A Novel 3D Reinforced Particle Model for Reinforced Concrete Fracture Assessment: Formulation and Validation

Abstract

Rigid particle models (PMs) that explicitly consider the influence of the aggregate structure and its physical interaction mechanisms have been used to predict cracking phenomena in concrete. PMs have also been applied to reinforced concrete fracture, but the known studies have adopted simplified reinforcement and reinforcement/particle interaction models. In this work, a novel 3D explicit discrete element formulation of reinforcement bars discretized through several rigid cylindrical segments is proposed, allowing the 3D reinforced particle model (3D-RPM) to be applied to reinforced concrete fracture studies, namely for shear failure. The 3D-RPM is evaluated using known three-point and four-point bending tests on reinforced concrete beams without stirrups and on known shear transfer tests due to dowel action. The 3D-RPM model is shown to reproduce the crack propagation, and the load displacement response observed experimentally for different steel contents under three-point bending, for different beam sizes, under four-point bending, and for different bar diameters, under shear. The validation examples highlight the importance of including a nonlinear reinforcement/particle interaction model. As shown, an elastic model contact leads to higher vertical loads in three-point and four-point bending tests for the same set of contact properties and, in the shear tests, leads to an overestimation of the maximum shear strength and to an increase in the model initial stiffness.

1. Introduction

The assessment of reinforced concrete structures requires numerical tools that can represent the formation, propagation and localization of cracks, the concrete/reinforcement interactions and their effect on the overall behavior. It is difficult to characterize this complex behavior using a continuum displacement formulation, as a suitable stress–strain law may not exist, or the constitutive law may be excessively complex, especially under shear loading [1].

Detailed rigid particle models (PMs) based on the discrete element method (DEM) have been shown to be an interesting alternative for fracture assessment given that the development of failure surfaces and their interaction occurs naturally as the structure is initially idealized as a discontinuous media. By considering the material grain structure and the physical mechanisms of grain interaction, PMs include in a direct way the material randomness and internal length. PMs based on the DEM were initially applied to rock fracture [2,3,4,5,6] and later extended to other quasi-brittle materials such as concrete [7,8,9,10,11,12,13], reinforced concrete [14,15], old masonry walls [16,17,18,19] and asphalt mixtures [20,21,22].

Several enhancements have been proposed in both 2D and 3D to address some of the limitations of the initial PM formulation for fracture studies [2], namely more complex particle shapes closer to the real grains [4,5,10,11,23], heterogeneity [6,8,17,19], adopting contact models that allow moment transmission [11], increasing particle contact interaction range [3,7,9,11], considering the real material internal structure through microcomputed tomography [10,24,25,26], adopting complex contact constitutive models including softening branches or strain rate effects [9,11,12,13,14,15], adopting viscoelastic contact models [20,22], including particle deformability [5,23,27] and adopting hybrid discrete element/finite modeling approaches [18,21].

When compared with meshless-based particle models based on a continuum approach, such as smoothed particle hydrodynamics (SPH) [28,29], the material point method [30,31], the extended finite element method [32] and the peridynamics approach [1,33], which allow discontinuities to occur and propagate, DEM-based PMs have the advantage of considering in a direct way the heterogeneous internal structure of concrete and the particle geometry, which is important to model the aggregate interlock mechanism present in reinforced concrete.

To model reinforced concrete structures, it is necessary to consider the reinforcement model and to account for the reinforcement/concrete interaction mechanisms. Three different approaches to represent reinforcement bars, depending on the type of problem being assessed, are usually adopted in continuum-based finite element (FEM) models [34]: (i) a discrete representation, usually adopting 1D bars that are connected to the concrete through linkage elements (springs) that may have bond/slip laws in the axial direction; (ii) a distributed representation, where the reinforcement bars are smeared over the concrete elements, a perfect bond is assumed and a composite stiffness matrix is usually built; and (iii) an embedded representation, where the reinforcement bars are directly modeled with their axial stiffness and a perfect bond is assumed.

In DEM, reinforcement bars are first modeled in 2D as a special contact between two discrete entities, using the same procedure as adopted for the traditional block-to-block contact [35]. A similar procedure has been adopted in the 3D numerical modeling of reinforced concrete tension stiffening tests [36]. For 2D DEM-based PMs, reinforcement has been considered by adopting in the particles that are intersected by the reinforcement bar an additional contact that adopts similar principles to the single contact point usually adopted in PMs [37]. In [15], a 2D reinforcement model that discretizes the reinforcements bars with 1D rigid lines and adopts a 1D reinforcement line/particle interaction through axial and shear springs was proposed and later applied to old masonry strengthening analysis [19]. In 3D DEM-based PMs, the most straightforward way to include the reinforcement is to discretize the reinforcement bars using lines of spherical particles of the same diameter as the reinforcement bar diameter [14,38], with mechanical contact properties that match the reinforcement properties. This approach may lead to a higher number of particles to be adopted given the particle size (reinforcement diameter), and the reinforcement/particle interaction inherits an artificial boundary roughness.

In similar mesoscale models based either on the Lattice model [39], on the Rigid Body Spring Model (RBSM) [40] or on the Lattice Discrete Particle Model [41], steel reinforcement is usually represented by two different approaches: (i) one-dimensional bar elements, Euler or Timoshenko beams, that interact with the concrete particles via a spring interface that may include nonlinear bond/slip laws [42,43,44,45,46]; (ii) the reinforcement bar is represented using a finer particle discretization of the steel reinforcement geometry [47,48], adopting different material properties for the particles representing the reinforcement bars and for the particles representing concrete, and the reinforcement/concrete interface is handled as a usual particle/particle interaction.

Three-dimensional PMs’ applications to reinforced concrete fracture studies have adopted simplified representations of the reinforcement and of the reinforcement/particle interaction models [14,37,38]. At mesoscale studies, specially under shear loading conditions, the representation of the reinforcement with an array of spherical particles that have an associated artificial roughness [14,38], or with a 1D reinforcement line that interacts with particles through interfaces [37], needs to be improved. To tackle this issue, a novel 3D reinforcement/particle model (RPM) based on the 2D model [15] is proposed for reinforced concrete fracture studies. In the proposed 3D-RPM, the reinforcement bars are discretized with cylindrical rigid elements that interact at the end nodes. Each cylindrical rigid bar element is sub-discretized with spherical particles along its length, which are rigidly associated with the rigid cylindrical bar. The spherical particles are only used for contact interaction purposes. With the proposed modeling approach, the reinforcement bar model and the reinforcement bar discretization required for contact purposes with the surrounding concrete are made independent. The artificial roughness associated with sub-discretizing a cylindrical rigid element with spherical particles for contact interaction purposes is smoothed by removing from the reinforcement/particle contact unit normal the vector component parallel to the rigid element axial direction.

The 3D-RPM is validated using known three-point and four-point bending tests on reinforced concrete beams and on known shear transfer tests due to dowel action. As shown, the 3D-RPM can reproduce the crack formation, propagation and localization and the load displacement response observed experimentally in the different tests, for different steel contents (three-point and dowel tests), and for different beam sizes and ratios (four-point tests). The numerical studies presented also show the relevance of adopting a modeling approach that allows for bond/slip and where the mechanical reinforcement behavior is made independent of the contact interaction.

2. Three-Dimensional Reinforced Particle Model (3D-RPM) Formulation

2.1. Particle Model (PM) Based on the Discrete Element Method

In a 3D-PM based on the discrete element method (DEM), the structure is represented as an assemblage of particles that interact with the neighboring particles at the contact points or interfaces. Given computational restrictions, the particles are usually assumed to be rigid and to have a spherical shape [2].

In a 3D-PM DEM-based model, it is assumed that during a single timestep, the disturbances can only propagate from a particle to its immediate neighbors, following a soft contact approach where the overlap of the particles is accepted. This overlap is not a real overlap, but it intends to model, in an indirect way, the deformation of the interacting particles at the contact point. The motion of a single particle is governed by the resultant moments and forces acting upon it. Newton’s second law of motion is integrated twice to define the particle displacements and rotations. In a 3D-PM DEM, a centred-difference scheme is usually adopted leading to the following expressions for the particle velocities at time $t+\mathsf{\Delta }t/2$:

$$x_i^{t+\mathsf{\Delta }t/2}=x_i^{t-\mathsf{\Delta }t/2}+{\displaystyle \frac{F_i\left(t\right)+F_i^d\left(t\right)}{m}}\mathsf{\Delta }t$$

(1)

$${\omega }_i^{t+\mathsf{\Delta }t/2}={\omega }_i^{t-\mathsf{\Delta }t/2}+{\displaystyle \frac{M_i\left(t\right)+M_i^d\left(t\right)}{I}}\mathsf{\Delta }t$$

(2)

where $F_i\left(t\right)$ and $M_i\left(t\right)$ are the total applied force and momentum applied at instant $t$, $F_i^d\left(t\right)$ and $M_i^d\left(t\right)$ are the damping force and momentum applied at instant $t$, $m$ and $I$ correspond to the mass and inertia of the particle, $x_i^{t+\mathsf{\Delta }t/2}$ and ${\omega }_i^{t+\mathsf{\Delta }t/2}$ are the velocities at time $t+\mathsf{\Delta }t/2$, and $x_i^{t-\mathsf{\Delta }t/2}$ and ${\omega }_i^{t-\mathsf{\Delta }t/2}$ are the velocities at time $t-\mathsf{\Delta }t/2$.

The uniaxial tests here presented were carried out under displacement control, adopting a constant plate velocity that assures quasi-static loading conditions. A local non-viscous damping approach with a damping value of 0.70 was adopted [49]. The three-point and four-point bending tests follow an adaptive dynamic relaxation viscous damping (ADR) approach under, respectively, displacement and load control, allowing larger loading steps to be adopted and faster convergence rates [15]. In the shear transfer due to dowel action tests, a similar procedure to that adopted in the uniaxial tests was adopted.

2.1.1. Voronoi-Generalized PM (VGCM-3D)

The adopted particle interaction follows a 3D Voronoi-generalized contact model (VGCM-3D) [11]. The Laguerre–Voronoi tessellation of the gravity centres of the interacting spherical particles is adopted to set the contact surface and the contact location of the two interacting particles, given the common Voronoi facet, as shown in Figure 1. Compared with complex polyhedral-based PMs [23,36], a PM that adopts a VGCM-3D keeps the simplicity of spherical particle interaction [2], reducing the need for a substantial increase in computational effort.

For a given local contact point, with a local normal ($k_n^{\left[J\right]}$) and shear stiffness ($k_s^{\left[J\right]}$), the local normal and shear forces increments are obtained following an incremental linear law:

$${∆F}_n^{[J,N]}=-k_n^{\left[J\right]}{∆x}^{[J,N]}=-k_n^{\left[J\right]}\left({\dot{x}}_i^{\left[J\right]}∆t\right)n_i$$

(3)

$${∆F}_i^{[J,S]}=-k_s^{\left[J\right]}{∆x}_i^{\left[J,S\right]}=-k_s^{\left[J\right]}\left({\dot{x}}_i^{\left[J\right]}∆t{-{∆x}^{[J,N]}n}_i\right)$$

(4)

where ${∆x}^{[J,N]}$ is the normal contact displacement increment, ${∆x}_i^{\left[J,S\right]}$ is the contact shear displacement increment, ${\dot{x}}_i^{\left[J\right]}$ is the contact displacement increment, and $n_i$ is the contact unit normal. The approach adopted in [2] is followed, so the normal components (forces and displacements) are stored as scalars, and the shear components are stored as vectors.

The VGCM-3D contact stiffness is set given the adopted Young’s modulus of the equivalent continuum material ($\overline{E}$) and given a constant that relates the normal and the shear stiffness spring value ($\alpha$):

$$k_n^{\left[J\right]}={\displaystyle \frac{\overline{E}}{d}}{A}_c^{\left[J\right]}$$

(5)

$$k_s^{\left[J\right]}=\alpha k_n^{\left[J\right]}$$

(6)

where $A_c^{\left[J\right]}$ is the local point j-associated contact area, and $d$ is the distance between the particles centre of gravity. The local-contact-point-associated area is defined by the total sum of one-third of the triangular areas associated with the given local contact point, see Figure 1. The local contact strength values are defined based on the maximum contact tensile stress (${\sigma }_{n.t}$), the maximum contact cohesion stress ($\tau$), and the local contact point area ($A_c^{\left[J\right]}$):

$$F_{n.max}^{\left[J\right]}={\sigma }_{n.t}{A}_c^{\left[J\right]}$$

(7)

$$C_{max}^{\left[J\right]}=\tau A_c^{\left[J\right]}$$

(8)

Given the contact maximum cohesion strength, the maximum shear contact force at the local point is given as follows:

$$F_{s.max}^{\left[J\right]}{=C}_{max}^{\left[J\right]}+F_n^{\left[J\right]}{\mu }_c$$

(9)

where ${\mu }_c$ is the contact friction coefficient.

2.1.2. Model Generation

In the adopted model generation, the aggregate content comprising one or more particle diameter dimension ranges is inserted from the highest to the smallest particle diameter range. A spherical shape is initially adopted. The spherical particles that represent the mortar are subsequently introduced adopting a porosity value of 0.1 and a uniform distribution, featuring diameters usually lower than the minimum diameter adopted to represent the aggregate content [11]. After a compact spherical particle assembly is formed, the centres of gravity of the spherical particles are triangularized based on a weighted Delaunay algorithm. The dual Laguerre–Voronoi diagrams are defined considering the constructed weighted Delaunay tetrahedra. A VGCM-3D contact between two neighboring particles is established given the associated Laguerre–Voronoi facet.

Figure 2 shows the adopted PM generation procedure regarding the uniaxial tests that are used to calibrate the contact properties for the three-point bending case study. For the aggregate, two sieve sizes were adopted, the largest sieve size with an aggregate content of 24% of the volume for particle diameters between 8.0 mm and 12.7 mm, and the smallest sieve size for particle diameters between 6.0 and 8.0 mm for an aggregate content of 20% of total volume. The particles representing the mortar were subsequently introduced adopting a porosity value of 0.1 and a uniform distribution, featuring diameters ranging between 5.6 and 6.4 mm [11]. Due to the computational costs associated with PMs, a particle size distribution larger than the actual concrete aggregate and mortar sizes was usually adopted. Additional DEM-based concrete model generation approaches can be found in [50].

2.1.3. Vectorial Bilinear Weakening Model (BL)

A vectorial bilinear softening contact damage model (BL) is adopted in the normal and shear directions [11,15] (Figure 3). The damage value in each contact direction, tensile–normal and shear, is updated given the maximum contact displacement registered in that direction, and the contact damage is given by the sum of the tensile damage and the shear damage.

The BL contact model has been shown to be a reliable choice for PM fracture studies, offering several advantages in detailed 3D PM DEM-based models [11]. Compared with more sophisticated contact models [41], the BL contact model simplifies the calibration process by reducing the number of contact strength parameters and has lower associated computational costs, without comprising the predicted concrete macroscopic response for different loading scenarios [11].

2.2. Proposed Reinforcement Model

In this work, a new procedure is proposed that discretizes each reinforcement bar using several rigid cylindrical segments that interact with each other at the reinforcement bar ends, where the reinforcement elastic and strength properties are lumped. The novel explicit formulation for the reinforcement model here presented follows similar principles to that adopted in a 2D RPM [15]. A similar reinforcement 2D model can also be found in [51] within the RBSM framework, which is the first formulation that models 2D bars using rigid bodies connected by springs, for in- and out-of-plane frame vibration analysis and beam plastic collapse.

The contact interaction between two rigid cylindrical elements at a given interaction location is presented in Figure 4. As in the 3D-PM, each cylindrical rigid bar element has six degrees of freedom of displacement/rotation and force/moment at its centre of gravity.

The location of the interaction node ($x_c^{\left[C\right]}$) where the elastic and strength properties are lumped, is given by the coordinate averages of the end nodes of each rigid cylinder adopted in the discretization of the reinforcement element:

$$x_c^{\left[C\right]}=0.5\left(x_i^{\left[AJ\right]}+x_i^{\left[BI\right]}\right)$$

(10)

where $x_i^{\left[AJ\right]}$ is the location of the end node $J$ of rigid cylindrical element A, and $x_i^{\left[BI\right]}$ is the location of the beginning node $I$ of rigid cylindrical element B. For a given rigid cylindrical element, the locations of the beginning node and of the end node are defined given the rigid cylindrical element length ($L^{\left[A\right]}$) and the axial direction ($a_i^{\left[A\right]}$):

$$x_i^{\left[AI\right]}=x_i^{\left[AI\right]}-{\displaystyle \frac{L^{\left[A\right]}}{2}}{a}_i^{\left[A\right]}$$

(11)

$$x_i^{\left[AI\right]}=x_i^{\left[AI\right]}+{\displaystyle \frac{L^{\left[A\right]}}{2}}{a}_i^{\left[A\right]}$$

(12)

The axial direction at the contact location ($a_i^{\left[C\right]}$) is given by averaging the axial directions of the interacting cylindrical rigid discrete elements under analysis:

$$a_i^{\left[C\right]}=0.5\left(a_i^{\left[A\right]}+a_i^{\left[B\right]}\right)$$

(13)

At the interaction location, the elastic properties are concentrated into an axial stiffness $k_a$, a shear stiffness $k_s$, a bending stiffness $k_b$ and a torsional stiffness $k_t$ that are given as follows:

$$k_a={\displaystyle \frac{EA}{L}};k_s={\displaystyle \frac{GI}{L}};k_b={\displaystyle \frac{EI}{L}};k_t={\displaystyle \frac{EJ}{L}}$$

(14)

where L is the length associated with the interaction location given by $L=0.5\left(L^{\left[A\right]}+L^{\left[B\right]}\right)$, A is the cross-sectional area of the rigid cylindrical element given by $A={\displaystyle \frac{\pi D^2}{4}}$, D is the rigid cylindrical element diameter, I is the moment of inertia of the cross-section given by $I={\displaystyle \frac{\pi D^4}{64}}$, J is the polar moment of inertia of the cross-section of a circular section of diameter D given by $J={\displaystyle \frac{\pi D^4}{32}}$, $E$ is the Young’s modulus of the reinforcement bar, and $G$ is the shear modulus of the reinforcement bar.

At the contact location, the axial and shear contact force increments are obtained following an incremental linear law:

$${∆F}^{\left[a\right]}=-k_a\left({\dot{x}}_i∆t\right)a_i=-k_a{∆x}^{\left[a\right]}$$

(15)

$${∆F}_i^{\left[s\right]}=-k_t\left({\dot{x}}_i∆t{-{∆x}^{\left[a\right]}a}_i\right)=-k_t{∆x}_i^{\left[s\right]}$$

(16)

where ${∆x}^{\left[a\right]}$ is the contact displacement axial increment, ${∆x}_i^{\left[s\right]}$ is the contact shear displacement increment, and ${\dot{x}}_i$ is the contact displacement increment. Similarly, the contact location torsional increment and bending increment are also obtained following an incremental law:

$${∆M}^{\left[t\right]}=-k_{\theta }\left({\dot{\theta }}_i∆t\right)a_i=-k_t{∆\theta }^{\left[t\right]}$$

(17)

$${∆M}_i^{\left[b\right]}=-k_b\left({\dot{\theta }}_i∆t{-{∆\theta }^{\left[t\right]}a}_i\right)-k_b{∆\theta }_i^{\left[b\right]}$$

(18)

In the numerical examples presented in Section 3, a trilinear stress–strain diagram was adopted at the interaction node, see Figure 5.

2.3. Proposed Reinforcement/Particle Interaction Model

One possibility for handling the reinforcement bar/concrete interaction would be to adopt a cylindrical/spherical particle detection/interaction scheme. Instead, in this work, each cylindrical rigid bar element is sub-discretized with spherical particles along its length, that are rigidly associated with the rigid cylindrical bar, Figure 6.

With the adopted contact interaction scheme, it is possible to adopt for the interaction between the rigid cylindrical particles (reinforcement bars) and the rigid spherical particles (concrete) the simplified spherical particle/spherical particle interaction model as defined in Section 2.1. For the reinforcement/concrete interaction, only one local contact point is adopted. With the adopted approach, the reinforcement bar discretization with rigid cylindrical elements that control the reinforcement bar’s mechanical behavior is made independent of the particle discretization adopted for contact purposes.

The unit normal of the concrete spherical particle/rigid element spherical particle contact is corrected considering the axial direction of the cylindrical rigid element to which the rigid element spherical particle belongs. In this way, the roughness associated with sub-discretizing a cylindrical rigid element with spherical particles is smoothed, avoiding the numerical influence of artificial particle interlocks.

The elastic and strength properties of the reinforcement cylindrical rigid element/concrete spherical particle contacts are defined using the same methodology adopted for the particle/particle contacts, see Section 2.1. Similarly, a BL contact model is adopted for reinforcement/concrete interactions. The adopted algorithm also ensures that a given concrete particle can only have one single contact point with a rigid cylindrical element sub-discretized with spherical particles.

3. Validation Examples

3.1. Methodology

The 3D-RPM with reinforcement and reinforcement/particle interaction capabilities is validated using a three-point bending reinforced concrete beam model [52] and a four-point bending reinforced concrete beam model [53]. A known shear transfer test due to dowel action only was also assessed [54].

The three-point and shear-point experimental tests that were taken as a reference [52,53] were chosen because they are well-known studies of reinforced concrete beams with only longitudinal reinforcement, that make it simpler to validate the 3D-RPM, but also introduce a required complexity, given that it handles different steel contents and beam geometries. The shear test taken as reference [54] was chosen because it is a well-known study that addressed shear transfer in concrete panels by first assessing aggregate interlock and dowel action on its own, and later assessed both phenomena contributions in reinforced concrete cracked panels [55].

In all numerical examples presented, the VGCM-3D contact elastic and strength parameters were calibrated using similar geometries to those adopted in the experimental studies that were taken as a reference, following a trial-and-error procedure, which relies on the user experience. It is important to mention that there are contact calibration approaches that can ease the calibration process and reduce the user expertise requirements, such as those based on design of experiment method (DOE), on machine learning methods (ML), or on evolutionary approaches (GA) [56,57,58].

Given the adopted particle size, it was decided not to represent the concrete heterogeneity; for this reason, the same contact properties were adopted independently of the particle type (aggregate/mortar).

3.2. Three-Point Bending Test

The 3D-RPM with reinforcement and reinforcement/particle interaction capabilities was applied in the assessment of a three-point bending test geometry that was adopted in the analysis of minimum reinforcement area in high-strength concrete [52]. Only the smallest beam size from the experimental study carried out in [52] was numerically assessed, Figure 7.

The concrete, made of crushed aggregate with a maximum diameter size of 12.7 mm, had an average maximum compressive stress of 91.2 MPa, obtained using cubic specimens measuring 160 mm, and a modulus of elasticity of 34.2 GPa. A Poisson’s coefficient of 0.18 and a maximum tensile strength of 4.6 MPa were assumed for contact calibration purposes.

As pointed out in Section 2.1.2, two sieve sizes were adopted, the largest sieve size representing particles with diameters ranging between 8.0 mm to 12.7 mm had an aggregate content of 24% of the total volume, and the smallest sieve size representing particles with diameters ranging between 6.0 and 8.0 mm had an aggregate content of 20% of the total volume. A uniform distribution with particle diameters ranging between 5.6 and 6.4 mm was adopted to represent the mortar. The best-fit values of the contact elastic and strength properties were obtained through a trial-and-error calibration process on 160 mm cubic specimens, using uniaxial compression and tensile tests, see

Table 1. Three-point bending test: 3D-RPM best-fit contact properties ¹.

Elastic		Strength
$\overline{\mathit{E}}$ $\left(\mathbf{M}\mathbf{P}\mathbf{a}\right)$	$\mathit{\alpha }$	${\mathit{\sigma }}_{\mathit{n}.\mathit{t}}$ $\left(\mathbf{M}\mathbf{P}\mathbf{a}\right)$	$\mathit{\tau }$ $\left(\mathbf{M}\mathbf{P}\mathbf{a}\right)$	${\mathit{\mu }}_{\mathit{c}}$	${\mathit{G}}_{\mathit{f}.\mathit{n}}$ $\left(\mathbf{N}/\mathbf{m}\mathbf{m}\right)$	${\mathit{G}}_{\mathit{f}.\mathit{s}}$ $\left(\mathbf{N}/\mathbf{m}\mathbf{m}\right)$
52.88	0.25	6.50	32.40	0.3	0.0082	6.084

¹ Macroscopic numerical predictions: Maximum tensile strength (4.6 MPa), Maximum Compression strength (90.6 MPa), Elastic modulus (34.2 GPa) and Poisson’s coefficient (0.18).

.

In all numerical examples, the reinforcement bars with a cross-sectional area equal to that adopted in the experimental studies were discretized with 32 mm long cylindrical rigid elements. For interaction reasons, a 0.50 radius overlap was adopted for the sub-discretization of the cylindrical rigid elements with spherical particles. A Young’s modulus of 200 GPa was assumed for the steel bars. As in the experimental study [52], the distance from the reinforcement bars to the lower beam face was set equal to one-tenth of the beam depth. In the numerical examples, the steel content and reinforcement bar sizes considered (1Φ4, 2Φ5, 2Φ8, 2Φ10) were similar to those adopted in the experimental studies [52].

As a reference, the three-point bending tests with a 5 mm reinforcement diameter had around 11,495 particles representing the aggregates and 46,495 particles representing the mortar, corresponding to a total of approximately 385,000 VGCM-3D contacts (≈2,320,000 local contacts). The reinforcement/particle interaction was represented by approximately 5300 single-contact-point VGCM-3D contacts.

Regarding the reinforcement/particle interaction, three different numerical approaches were adopted: (i) an elastic model (EL); (ii) a BL model with the same strength properties as those adopted for the concrete contacts (NL); and (iii) a BL model with a 50% reduction in the contact strength properties adopted for concrete (NLR).

Figure 8 shows the NLR approach for each steel content and adopted reinforcement bar size, the predicted numerical final deformation magnified 20 times, and the predicted final crack patterns. In all numerical tests, cracking first occurred at the bottom zone close to the zone of maximum bending moment (beam midspan). Later, the cracked surface developed towards the upper plate where the vertical load was applied. As shown in Figure 8, for higher reinforcement steel ratios, secondary diagonal cracking also occurs for higher loading values. Higher reinforcement steel ratios allow the load to be transmitted from the upper plate at midspan to the lower supporting plates through an arch effect; the higher the steel ratio, the higher is the number of diagonal cracks that can be formed.

Figure 9 shows the vertical load–midspan displacement numerical predictions for the different adopted reinforcement steel content: without steel (WS); 1Φ4; 2Φ5; 2Φ8; and 2Φ10. For each reinforcement steel content, the numerical results obtained with each adopted reinforcement/particle interaction approach are compared: (i) elastic model (EL); (ii) BL model with the same strength properties as those adopted for the concrete contacts (NL); and (iii) BL model with a 50% reduction in the contact strength properties adopted for concrete contacts (NLR). Also presented are the experimental curves [52] that were obtained for each tested reinforcement content (Exp (1990)).

For the case without reinforcement, Figure 9a, the numerical prediction (WS) is shown to have a good agreement with the experimental response (WS-Exp), the numerical response being slightly more brittle than that observed experimentally. If a more refined PM was adopted and the concrete heterogeneity represented, it would be possible to numerically predict a more ductile post-peak behavior.

When compared with the elastic reinforcement/particle contact model (EL), the predictions obtained with a BL contact model (NL and NLR) are in better agreement with the experimental tests, for all the assessed reinforcement scenarios (Figure 9b–e). As shown, a BL contact model with a 50% strength reduction (NLR) predicted the best agreement with the known experimental results for all the steel reinforcement contents (Figure 9b–e).

Figure 9f shows the 3D-RPM predictions with a 50% strength reduction for the reinforcement steel contents under assessment (NLR). The numerical tests clearly indicate that the transition from a brittle to a more ductile response occurs for a level of reinforcement corresponding to 2Φ5, similarly to that obtained in the experimental studies [52].

3.3. Four-Point Bending Test

A four-point bending configuration, as shown in Figure 10, was also assessed using the proposed 3D-RPM. The test geometry was adopted to study the size effect on diagonal shear failure of reinforced concrete beams without stirrups [53]. The micro-concrete, composed of aggregates with a maximum size of 4.8 mm, exhibited an average compressive strength of 46.8 MPa, as measured on cylindrical specimens with a diameter of 76 mm and a height of 152 mm. A Poisson’s coefficient of 0.18, a Young’s modulus of 36.2 GPa, and a maximum tensile strength of 3.2 MPa were assumed for contact calibration purposes. Like in the experimental setup, the numerically tested geometries had a constant thickness of 38.1 mm. In the numerical models, a ${\sigma }_1$ value of 790.0 MPa was adopted, which is within the range of yield strengths reported, 690.0 MPa to 890.0 MPa [53].

In the PM, the aggregate was represented with a single sieve size that had an aggregate content of 35% of the volume representing particles with diameters ranging between 5.0 mm and 7.2 mm. A uniform distribution with particle diameters ranging between 4.6 and 5.2 mm was adopted to represent the mortar. The best-fit values of the contact properties were obtained through a trial-and-error calibration process on 76 mm × 76 mm × 152 mm specimens, using uniaxial compression and tensile tests,

Table 2. Four-point bending test: 3D-RPM best-fit contact properties ¹.

Approach	Elastic		Strength
	$\overline{\mathit{E}}$ $\left(\mathbf{M}\mathbf{P}\mathbf{a}\right)$	$\mathit{\alpha }$	${\mathit{\sigma }}_{\mathit{n}.\mathit{t}}$ $\left(\mathbf{M}\mathbf{P}\mathbf{a}\right)$	$\mathit{\tau }$ $\left(\mathbf{M}\mathbf{P}\mathbf{a}\right)$	${\mathit{\mu }}_{\mathit{c}}$	${\mathit{G}}_{\mathit{f}.\mathit{n}}$ $\left(\mathbf{N}/\mathbf{m}\mathbf{m}\right)$	${\mathit{G}}_{\mathit{f}.\mathit{s}}$ $\left(\mathbf{N}/\mathbf{m}\mathbf{m}\right)$
NL	55.10	0.25	4.55	13.3	0.5	0.0035	0.8835
EL			3.90	13.8	0.5	0.0026	0.9511

¹ Macroscopic numerical predictions: Maximum tensile strength (NL: 3.2 MPa; EL: 2.75 MPa), Maximum Compression strength (46.8 MPa), Elastic modulus (36.2 GPa), and Poisson’s coefficient (0.18).

. Also presented are the best-fit contact properties that give the best maximum load predictions when adopting an elastic model for the reinforcement/particle interaction (EL).

Two different beam sizes were numerically tested, the smallest size adopted in the experimental study (d = 40.64 mm) and an intermediate size (d = 60.96 mm) that corresponds to an intermediate size between the smallest size (d = 40.64 mm) and the following size adopted in the experimental study taken as reference (d = 81.28 mm) [53].

The adopted reinforcement bars were discretized with 10.16 mm long cylindrical rigid elements. For the intermediate beam size that was only tested numerically, the beam geometry and the steel reinforcement were defined given the average values of the experimental tests (d = 40.64 mm and d = 81.28 mm). The cross-sectional areas are equal to those adopted in the experimental studies (d = 40.64 mm: 2Φ6.19 and d = 60.96 mm: 2Φ8.34). For interaction reasons, a 0.50 radius overlap was adopted for the sub-discretization of the cylindrical rigid elements with spherical particles. A Young’s modulus of 200 GPa was assumed for the reinforcement steel bars.

As a reference, the shear-point bending tests with the highest numerically adopted beam depth (d = 60.96 mm) had around 3198 particles representing the aggregates and 13,894 particles representing the mortar, corresponding to a total of approximately 105,500 VGCM-3D contacts (≈633,000 local contacts). The reinforcement/particle interaction was represented by approximately 2550 single-contact-point VGCM-3D contacts. Regarding the reinforcement/particle interaction, two different numerical approaches were adopted: a BL model with the same strength properties as adopted for the concrete contacts (NL), and an elastic contact model (EL) with slightly different contact strength properties that were required to be adjusted to match the maximum load predictions.

Figure 11 shows the predicted crack pattern evolution and the final displacement (magnified 10 times), for the smallest beam size numerically tested (d = 40.64 mm), when an NL reinforcement/particle model is adopted. As shown, several tensile cracks occur at the lower face in the zone of maximum bending moment (central part midspan in between the upper plates). Following this, some of the cracks grow from the lower face towards the upper plates, initially in a vertical direction and later with a diagonal curved shape. At higher loads, additional diagonal cracks due to tensile loading occur. Up to the maximum load, there is some symmetry in the final cracking patterns, but at a later stage, a diagonal failure tends to localize in one of the sides, see also Figure 11f.

Figure 12 shows the predicted crack pattern evolution and the final displacement (magnified 10 times), for the largest beam size numerically tested (d = 60.96 mm), when an EL reinforcement/particle model is adopted; a similar crack propagation and localization occurs.

Table 3. Maximum vertical load: 3D-RPM best-fit contact properties and experimental results [44].

$\mathbf{Beam}\mathbf{Size}\left(\mathbf{m}\mathbf{m}\right)$	Pexp (kN)	PNum.NL (kN)	PNum.EL (kN)
40.64	5.92–6.49	6.20	6.20
60.96	8.45–8.83	8.50	8.50

presents, for each beam numerically assessed (40.64 mm and 60.96 mm), the maximum load numerical predictions for the reinforcement/particle interaction when adopting a BL model with the same strength properties as adopted for the concrete contacts (P_Num._NL) and when adopting a reinforcement/particle interaction elastic model (P_Num._EL). Also presented are the range of maximum loads obtained experimentally (P_Exp) [44]. As shown in Table 3, the 3D-RPM model’s maximum load predictions for each beam size are in good agreement with the maximum values obtained experimentally for both EL and NL reinforcement/particle contact models.

Figure 13a,b show for each numerically assessed beam size (40.64 mm and 60.96 mm) the maximum load versus beam size predicted diagrams and the total damage evolution for, respectively, the reinforcement/particle NL contact model (NL) and for the reinforcement/particle EL contact model (EL). The numerical predictions of load displacement and total contact damage evolution are very similar for both reinforcement/particle models (NL and EL) for both beam sizes (40.64 mm and 60.96 mm). For both geometries and contact model approaches, the total contact damage at peak value is of the order of 30%, indicating that extensive cracking occurs in the reinforced concrete beams at peak load as also shown in Figure 11e and Figure 12e.

3.4. Shear Transfer Due to Dowel Action Test

The proposed 3D-RPM has also been used for the analysis of shear transfer tests in reinforced concrete panels due to dowel action [54]. In the dowel tests with the plate geometry presented in Figure 14, the aggregate interlock effects were mitigated by constructing a smooth low-friction crack passing through the center of the specimen and adopting a two-stage casting procedure [54]. The mean concrete cube strength measured in 100 mm cubes varied between 27.6 MPa and 37.6 MPa. A Poisson’s coefficient of 0.18, a Young modulus of 30.0 GPa, a maximum tensile strength of 3.5 MPa and a maximum compressive strength of 32.50 MPa were assumed for contact calibration purposes. In the numerical models, a ${\sigma }_1$ of 400 MPa was adopted, which is within the range of the yield strengths reported in the experimental studies that were carried out in reinforced concrete panels by the same authors in a follow-up study that assessed shear transfer in reinforced concrete panels [55].

In the PM, the aggregate was represented with two sieve sizes, a larger one that had an aggregate content of 25% of the volume representing particles with diameters ranging between 10.0 mm and 12.0 mm, and a smaller one that had an aggregate content of 30% of the volume representing particles with diameters ranging between 8.0 mm and 10.0 mm. A uniform distribution with particle diameters ranging between 7.0 and 9.0 mm was adopted to represent the mortar. The best-fit values of the contact properties were obtained through a trial-and-error calibration process on 100 mm cubes, using uniaxial compression and tensile tests,

Table 4. Shear transfer due to dowel action test: 3D-RPM best-fit contact properties ¹.

Elastic		Strength
$\overline{\mathit{E}}$ $\left(\mathbf{M}\mathbf{P}\mathbf{a}\right)$	$\mathit{\alpha }$	${\mathit{\sigma }}_{\mathit{n}.\mathit{t}}$ $\left(\mathbf{M}\mathbf{P}\mathbf{a}\right)$	$\mathit{\tau }$ $\left(\mathbf{M}\mathbf{P}\mathbf{a}\right)$	${\mathit{\mu }}_{\mathit{c}}$	${\mathit{G}}_{\mathit{f}.\mathit{n}}$ $\left(\mathbf{N}/\mathbf{m}\mathbf{m}\right)$	${\mathit{G}}_{\mathit{f}.\mathit{s}}$ $\left(\mathbf{N}/\mathbf{m}\mathbf{m}\right)$
48.15	0.20	5.40	10.95	0.2	0.0092	1.411

¹ Macroscopic numerical predictions: Maximum tensile strength (3.5 MPa), Maximum Compression strength (32.5 MPa), Elastic modulus (30.0 GPa) and Poisson’s coefficient (0.18).

. In the 3D-RPM, the smooth joint effect was modeled by imposing to the VGCM3D contacts in the smooth joint vicinity a constant horizontal unit normal and a zero shear stiffness spring value, only allowing normal contact forces to be transmitted through the smooth joint.

The reinforcement steel bars were discretized with cylindrical rigid elements with a length of 2.25 mm. Like in the previous examples, for interaction reasons, a 0.50 radius overlap was adopted for the sub-discretization of the cylindrical rigid elements with spherical particles. The shear tests due to dowel action with a 16 mm reinforcement diameter had approximately 4100 particles representing the aggregates and 46,495 particles representing the mortar, corresponding to a total of 320,700 VGCM-3D contacts, which, in turn, corresponded to a total of 1,924,200 local contacts. The reinforcement/particle interaction was represented by around 2400 single-contact-point VGCM-3D contacts.

Regarding the reinforcement/particle interaction, two different reinforcement/particle models were adopted: an elastic contact model (EL) and BL model with the same strength properties as adopted for the concrete inter-particle contacts (NL). Additionally. for the NL reinforcement/particle interaction model, the influence of the adopted ${\sigma }_1$ for the steel strength was also assessed by adopting a lower ${\sigma }_1$ value of 380 MPa (NL-S1) and a higher ${\sigma }_1$ value of 420 MPa (NL-S2). The reinforcement steel bar level of discretization was also assessed for the NL reinforcement/particle model by also evaluating a 3D-RPM that adopts cylindrical rigid element discretization with a length of 15.0 mm (NLd).

The predicted numerical shear displacement/shear stress diagrams for a reinforcement steel bar diameter of 12 mm (Φ12) and for a reinforcement steel bar diameter of 16 mm (Φ16) are presented in Figure 15, for the EL contact model (EL), for the NL contact model (NL), for the NL contact model with a lower ${\sigma }_1$ value (NL-S1), for the NL contact model with a higher ${\sigma }_1$ value (NL-S2), and for the NL contact model that adopts cylindrical rigid elements discretization with a length of 15.00 mm (NLd). Also shown are the experimental curves [45]. As shown in Figure 15, a reasonable agreement is obtained between the 3D-RPM predictions and the experimental results for the two reinforcement bar diameters evaluated. For both reinforcements, it is shown that an elastic reinforcement/particle model, Φ12 (EL) and Φ16 (EL), leads to an increase in the model initial stiffness and in the predicted maximum shear strength.

Figure 15a,d also show that the reinforcement discretization, Φ12 (NLd) and Φ16 (NLd), influences the numerical predictions, namely for stiffer reinforcements with larger diameter, where the predicted numerical response switches from a ductile response to a more brittle response, Φ16 (NLd) in Figure 15d. This shows the relevance of adopting a reinforcement/particle interaction model that independently handles the contact interaction and the bar mechanical behavior. In zones where a very high reinforcement deformation is expected, such as a plane joint, it is important to adopt a reinforcement bar discretization length much lower than the adopted minimum particle radius. For both adopted reinforcement diameters, the predicted maximum shear stress is almost proportional to the adopted steel bar ${\sigma }_1$ value (Figure 15b,e).

Figure 15c,f show, for both reinforcement diameters and for the NL reinforcement/particle model, Φ12 (NL) and Φ16 (NL), the total damage evolution (dashed line) and the reinforcement/particle contact damage evolution (dashed-point line). For both reinforcement diameters, the damage at the reinforcement/particle contacts occurs almost from the onset of loading, and closer to the shear yield plateau, the contact damage at the reinforcement/particle contacts tends to stabilize (10% for Φ12 and 12% for Φ16). Figure 15d,f also show that a higher contact damage is associated with a larger bar diameter, similar to what occurred in the experimental tests [54], and that at the shear yield maximum strength plateau, the damage evolution also tends to stabilize. The total damage results also indicate that the damage is localized with only 2.6% total damage predicted with the Φ16 (NL) model.

Table 5. Maximum dowel strength: Theoretical and numerical values ¹.

$\mathsf{\varphi }$ $\left(\mathbf{m}\mathbf{m}\right)$	F Dowel.theory (MPa)		τ Theory (MPa)		τ NL (MPa)
	${\mathit{\sigma }}_1$	${\mathit{\sigma }}_2$	${\mathit{\sigma }}_1$	${\mathit{\sigma }}_2$
12.0	21.34	23.86	1.42	1.59	1.41
16-0	37.95	42.42	2.53	2.83	2.59

¹ F _Dowel.theory = $1.3{\mathsf{\varphi }}^2{\sigma }_y^{0.5}{\sigma }_c^{0.5}$ [54] where ${\sigma }_c$ is the concrete compressive strength measured in cubes, and ${\sigma }_y$ is the yield stress of the reinforcing steel.

presents the maximum theoretical shear dowel strength (F _Dowel.theory) and the corresponding maximum theoretical shear stress (τ _Theory) [54], given the adopted numerical concrete macroscopic properties, Table 3, and the adopted values of the trilinear steel reinforcement model, ${\sigma }_1$ and ${\sigma }_2$ and the predicted numerical values with an NL reinforcement/particle interaction model (τ _NL), Φ12 (NL) and Φ16 (NL). As shown, the numerical predicted maximum shear strength, for each reinforcement size, is within the range of the maximum theoretical shear strength values.

Figure 16 shows the predicted numerical final deformation, magnified 10 times, and the predicted final crack patterns for the Φ16 model and for the NL reinforcement/particle interaction model. As shown, cracking occurs in the vicinity of the reinforcement bars due to compression leading to a splitting failure that grows towards the supported faces and towards the specimen interior; a similar phenomenon was recorded in the experimental tests [54].

Figure 16b also shows that the proposed reinforcement model discretized with rigid cylinders connected with springs at the end nodes can consider the reinforcement bar flexure, direct shear and kinking, which are known to be the most relevant factors in the contribution of the reinforcement steel to the shear strength.

4. Discussion

A novel 3D reinforced PM is proposed to model reinforced concrete fracture where (i) the concrete structure (aggregate/mortar) is approximately given by polyhedral shaped particles; (ii) the reinforcement bars are discretized with several rigid cylindrical segments that interact with each other at the reinforcement bar ends where the elastic and strength properties of the reinforcements are lumped; (iii) the cylindrical rigid elements are sub-discretized with spherical particles along their length to ease the contact interaction process, adopting a simplified spherical particle/spherical particle interaction model.

Compared with the usual 3D PM approach that discretizes the reinforcement bars using lines of spherical particles of diameter equal to the bar diameter [14,38] or with a 1D reinforcement line that interacts with particles through interfaces [37], the proposed novel 3D reinforced PM approach allows the following:

• The elimination of the outer artificial roughness associated with the spherical particle discretization by removing the cylindrical rigid element axial direction from the reinforcement/particle contact unit normal. This outer roughness is more complex to eliminate if a line of spheres is adopted to model concrete [14,38].
• Additional flexibility in the reinforcement bar discretization; finer discretizations lower than the reinforcement bar diameter can be adopted if large deformations are expected or larger discretization than the bar diameter can be adopted if lower deformation gradients are expected, with computational gains, giving an additional flexibility that is not possible if an array of particles is adopted to represent the reinforcement [14,38].
• Bond/slip behavior to be considered by taking into account, at the reinforcement/particle interface, the proper state of stress, including the confinement stress, which cannot be accomplished in a straightforward manner if a 1D line reinforcement bar that interacts with particles through interfaces is adopted to represent the reinforcement and the reinforcement/particle interaction [37].

The 3D-RPM crack pattern predictions in the three-point bending test are similar to those observed experimentally [32,59,60]. With a 2D-RPM [15], the predicted crack pattern is much more localized even for higher reinforcement contents; the crack pattern predicted with the novel 3D-RPM for higher reinforcement contents is in much better agreement with the crack pattern observed experimentally under three-point bending, where additional cracking and diagonal cracking are observed [32,59,60,61]. In the three-point bending test, it is shown that the reinforcement/particle interaction BL model with a 50% reduction in the contact strength properties adopted for concrete contacts (NLR) leads to a better agreement with the known experimental results for different steel contents [52].

In the four-point bending test, the predicted crack patterns and failure modes are in good agreement with those observed in experimental tests for the same beam sizes [53], in 2D numerical studies [15] and in similar experimental tests of beams with only longitudinal reinforcement [33,62]. The results show the relevance of adopting a 3D-RPM that includes bond/slip at the reinforcement/particle interfaces. In the four-point bending test, to have similar predictions, it is necessary to reduce the macroscopic tensile strength (15% reduction) when adopting an elastic reinforcement/particle model. The presented results also show that for this type of geometry, loading and support conditions, the discretizaton adopted for the reinforcement bars can use a length higher than the reinforcement bar diameter, reducing the modeling requirements, when compared with PMs that adopt an array of spherical particles [14,38].

Regarding the shear transfer due to dowel action tests, the results presented show that the proposed reinforcement model discretized with rigid cylinders connected with springs at the end nodes can consider the reinforcement bar flexure, direct shear and kinking, which are known to be the most relevant factors in the contribution of the reinforcement steel to the shear strength. The results presented show that the proposed 3D-RPM can be used as a local model to study and evaluate the effect of passive anchors as a possible reinforcement solution for small gravity dam stability [63] and to study the effect of the dowel bar load transfer efficiency in jointed plain concrete pavements [64]. Note that with a 2D-RPM [15] or with a 3D PM that adopts a 1D line reinforcement bar that interacts with particles through interfaces [37], it is not straightforward to model the dowel effect.

5. Conclusions and Further Developments

The results presented show that the 3D-RPM can be applied to fracture analysis of reinforced concrete geometries at the mesoscale to further understand the fracture process and to validate new reinforcement materials. The validation examples presented in Section 3, based on three-point and four-point bending tests on reinforced concrete beams and shear transfer tests due to dowel action, show the following:

• That with the proposed 3D-RPM, it is possible to calibrate the contact properties in simple tests (uniaxial compression and uniaxial tension) to numerically predict responses close to those observed experimentally in reinforced concrete specimens, for different steel contents, structural geometries and sizes and boundary and loading conditions.
• The relevance of adopting a reinforcement/particle contact model that considers bond/slip behavior. An elastic model contact leads to higher vertical loads in three-point and four-point bending tests for the same set of contact properties. In the shear transfer due to dowel action tests, an elastic reinforcement/particle model leads to an overestimation of the maximum shear strength and to an increase in the model initial stiffness.

Further work is underway to better calibrate the reinforcement/particle interaction contact properties with known tension stiffening experimental results and to adopt the 3D-RPM as a local model within a hybrid discrete element/finite element model to assess passive anchors as a reinforcement solution for small gravity dam stability analysis.