A discrete element model of concrete for cyclic loading

Computers and Structures 196 (2018) 173–185

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

A discrete element model of concrete for cyclic loading Sina Sinaie a, Tuan Duc Ngo a,⇑, Vinh Phu Nguyen b a b

Department of Infrastructure Engineering, The University of Melbourne, VIC 3010, Australia Department of Civil Engineering, Monash University, Melbourne, VIC 3800, Australia

a r t i c l e

i n f o

Article history: Received 23 June 2017 Accepted 27 November 2017

Keywords: Concrete Discrete element method Cyclic loading Damage model Residual strength

a b s t r a c t This paper takes advantage of the discrete element method to develop a model of concrete for cyclic simulations. For this purpose, a micro-mechanical damage model that also allows stress-reversals is formulated for inter-particle bonds. Moreover, a multi-phase implementation of the discrete element method is developed and used for two distinct reasons. First, to characterize aggregate and mortar particles separately. Second, to allow the effect of the interfacial transition zone to be taken into account. A strict validation approach is taken in this work, whereby the developed model is only calibrated against monotonic stress-strain curves and then evaluated for its performance under cyclic loading. Simulation results are constantly compared against experimental values. These comparisons illustrate the capability of the model to predict cyclic properties of concrete. Progression of damage is discussed in terms of numerical variables and also through the visualization of force chains and crack propagation. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction With the widespread use of concrete as a construction material in earthquake-prone regions, characterization of the cyclic properties of concrete has been the subject of many studies. These involve a range of experimental investigations carried out on plain concrete under cyclic loading [1–7]. As two of the pioneering contributors, Bahn and Hsu [1] characterized some important features of concrete under cyclic loading. What sets aside Bahn and Hsu’s [1] work is the fact that they also examined these features in relation to random cycles and evaluated the existing analytical models against such loadings. To make cyclic models of concrete applicable to the design of spirally reinforced concrete columns, Sakai and Kawashima [2] investigated the behavior of confined concrete under cyclic loading. As computational methods for simulating concrete structures were evolving, Sima et al. [3] aimed to improve the constitutive model of concrete that was being used in such computational methods, especially in relation to the smeared crack approach. Their constitutive model covered the cyclic behavior of concrete in the compressive and tensile regions and employed separate damage parameters for each region. The constitutive model of Sima et al. [3] has been improved by Breccolotti et al. [6] for a more accurate account of damage accumulation. ⇑ Corresponding author. E-mail addresses: [email protected] (S. Sinaie), [email protected] (T.D. Ngo), [email protected] (V.P. Nguyen). https://doi.org/10.1016/j.compstruc.2017.11.014 0045-7949/Ó 2017 Elsevier Ltd. All rights reserved.

More recently, Sinaie et al. [7,8] investigated the cyclic response of cylindrical concrete samples in the context of size effects and in combination with temperature effects. Aside from its standard use, concrete is also a main load-bearing component of many composite systems. Studies involving fiber-reinforced concrete [9] and concrete-filled tubes [10] have also dealt with the characterization of concrete under cyclic loading. When it comes to the simulation of concrete, meshless methods have shown to have certain advantages. This is especially true in relation to discontinuities (cracks), where such methods allow simple yet reliable formulations for the initiation and propagation of cracks. For example, Rabczuk and Belytschko [11] and Rabczuk et al. [12] incorporated their idea of representing cracks as a set of contiguous particles to simulate a range of 3D concrete problems involving static and dynamic loading with rate-dependent and independent material models. In their model, the crack surface is not restrained by the topology of the model, since the formation of a crack and its direction is determined by the failure criterion at the particle scale. Rabczuk et al. [13] later extended modified this approach by splitting the particles where cracking is detected. Doing so eliminated the need for additional variables at the location of the crack, therefore simplifying the formulation. Being a mesh-free particle-based method, the discrete element method (DEM) has also shown to be a reliable approach for the simulation of concrete. The DEM has also been coupled with other computational methods such as smooth-particle hydrodynamics (SPH) [14] for problems involving fluid-structure interaction and

174

S. Sinaie et al. / Computers and Structures 196 (2018) 173–185

computational fluid dynamics (CFD) for problems involving flow in porous material [15]. As one of the most earliest works, Cundall and Strack [16] evaluated a 2D DE model of packed disks. Although primarily aimed at granular assemblies such as sand, their work set the stage for the development of a computation model allowing interactions to be set at a particle scale, and therefore, applicable to various granular material. Due to its generality, aside from being extended into 3D and the addition of inter-particle bonds, the underlying concept of the discrete element method has not changed over the years. The main advancements made to the discrete element method have been in the direction of its applications, and the computational implementation, which in turn allows larger scale (or higher resolution) applications. At the time, Cundall and Strack [16] highlighted that their approach would allow a 1500-particle assembly to be simulated using no more that 64 kB of memory. Nowadays, simulations involving tens of thousands of particles is the norm and can easily be carried out on regular desktop computers. With the inclusion of inter-particle bonds into the DEM, such models have shown to be an effective solution to problems involving fracture [17–20]. This includes the simulation of concrete, as fracture and crack propagation play an important role in the behavior of this material. Hentz et al. [21,22] defined an interaction range as a criterion to establish inter-particle bonds (referred to interaction links in their work) between neighboring particles. With this method, they created particle assemblies to investigated the effectiveness of DE models in predicting the static [21] and dynamic [22] properties of concrete samples. Using a similar approach, Shiu et al. [23] used the discrete element method to simulate the penetration of missiles into reinforced concrete slabs. Aside from the particles that represented the concrete, they also included the reinforcement bars as a series of aligned particles and the impacting missile as a rigid assembly of particles. Note that the properties of these different materials are implemented into the model through the mechanical properties of the inter-particle bonds, which they referred to as parallel bonds. Up to this day, most researchers take a single-phase approach to their DE simulations, meaning that they have one single type of particles and one single type of bond between the particles. In contrast to this, Azevedo et al. [24] implemented a multi-phase approach to make a distinction between the aggregates and the cement matrix. Doing so gives the possibility to have a different sets of mechanical parameters (or even constitutive relations) depending on the type of the particles that exists at the two ends of the bond. Specifically for concrete, this means that the model can take into account three different phases corresponding to aggregate-aggregate, aggregate-mortar, and mortar-mortar bonds. Using this multi-phase approach, Sinaie et al. [25] investigated the thermal degradation of concrete from the standpoint of incompatibility between aggregates and the mortar matrix in terms of their coefficients of thermal expansion. Studies involving DE models of concrete have mostly dealt with monotonic loading, whereby micro-mechanical properties (defined at particle-scale) have been shown to properly translate into macro-mechanical characteristics (observed at sample-scale). However, this translation is yet to be shown for loading histories that involve stress reversals. To fill this gap, the present study explores the advantages/disadvantages of the DEM when it comes to simulating concrete under cyclic loading. To this end, a damagebased inter-particle bond model is developed to accommodate stress reversals and to take into account the degradation of micro-mechanical properties in compression, tension and shear. A multi-phase approach is also taken here, whereby the interfacial transition zone (ITZ) and the material weaknesses associated with it are introduced into the model. The numerical samples created and analyzed here are based on a series of laboratory experiments

carried out by the author in a previous study [7]. To avoid the imposition of an strenuous calibration requirement, material properties of the multi-phase DE model are calibrated against a fraction of the experimental results, and the remaining experimental data are used to validate the calibrated model. These validations clearly demonstrate the capability and reliability of the DEM for the simulation of concrete specimens under cyclic loading. It is important to note that the present paper aims to extend our previous work on DEM simulations of concrete samples [25,26]. For completeness, this paper is, however, treated independently and therefore, a complete description of the DEM model is given. Note that the micro-mechanical model that is formulated here is a refined version of the one used in Refs. [25,26]. At the macroscale (sample scale), the results of this study lay down the foundation for using the DEM to predict the response of concrete-like materials where stress reversals are involved. At the micro-scale, the findings from this and our previous works [25,26] provide a better understanding of the inter-particle interactions of such materials. 2. Multi-phase discrete element model This section describes the multi-phase DE model used in this study. The term multi-phase indicates that distinction is made between different types of particles (aggregates and mortar) and the interactions that can exist between them (aggregate-mortar and mortar-mortar bonds). Characterizing aggregate-mortar and mortar-mortar bonds separately allows the interfacial transition zone (ITZ) to be taken into account, which is a major factor influencing the strength and behavior of concrete. The overall mechanical properties of a particle assembly is the net result of the interactions between the assembly’s individual particles. These interactions are categorized as either contact or cohesive bond. Contact interactions appear when two particles touch each other, while cohesive bond interactions are only created if two particles are in close proximity at the beginning of the analysis. From a simulation standpoint, the main difference between the two categories is that contact interactions can dynamically appear and disappear during an analysis, however, if an existing cohesive bond interaction breaks, it is permanently lost. As mentioned before, this work uses a multi-phase implementation, which makes a distinction between aggregate particles and particles representing the mortar matrix [24,27]. Among the two, only mortar particles can form a cohesive bond with their neighbors. This is in contrast to single-phase models [28,29] where all particles are taken to be of the same material and all particles have the potential to form cohesive bonds. In the following text, properties associated with aggregate particles and mortar particles are distinguished by using the subscripts ‘a’ and ‘m’. 2.1. Kinematics of cohesive bonds As mentioned before, cohesive bonds are created only once, that is, at the start of the simulation. For a cohesive bond to exist between two particles, the two particles have to be in close proximity, and at least one of them has to be of mortar type. The spatial condition for a bond to exist between two particles i and j is:

jdij j 61þc ri þ rj

ð1Þ

where ri and r j are the radii of particles i and j, respectively, as illustrated in Fig. 1. The length of the distance vector jdij j ¼ jxj  xi j is calculated from the center of particle i (xi ) to the center of particle j (xj ). In Eq. (1), c > 0 is a parameter that allows a bond to be established when the two particles are in close proximity, but not

175

S. Sinaie et al. / Computers and Structures 196 (2018) 173–185

to a degradable-stiffness formulation. This is demonstrated later in Section 4 by comparing the results of the 2 approaches. As a result, in this study, stiffness values K n and K s are set to be independent of damage, and are calculated using:

dij

c

ri

Kn ¼

rj

i

dcj

necessarily in perfect contact. A value of c ¼ 0:25 is used in this study to allow each particle to be connected to its neighboring particles through an average of 10–11 bonds. The distance between the centroid of particle i and the contact point c is denoted by the vector dci and is equal to:




 ri dij ri þ rj

ð2Þ

The smaller of the two particles is used to calculate the cross section of the bond area: 2

Aij ¼ p½minðr i ; rj Þ

ð3Þ

Increments of displacements and rotations are calculated at every time-step, whereby Dx ¼ xt  xtDt and Dh ¼ ht  htDt . Therefore, the displacement of the contact point on particle i is:

Dxci ¼ Dxi þ ðdci  Dhi Þ

ð4Þ

and the relative displacement between the contact points on the two particles becomes Dxc;rel ¼ Dxcj  Dxci . The unit normal vector of a cohesive bond is defined as:

^ ij ¼ dij =jdij j n

ð5Þ

and consequently, the unit tangent vector is:

^sij ¼

^ ij Þn ^ ij Dxc;rel  ðDxc;rel  n ^ ij Þn ^ ij j jDxc;rel  ðDxc;rel  n

ð6Þ

2.2. Force-deformation relation for cohesive bonds Using the unit vectors defined in Eqs. (5) and (6), the force acting on particle i as a result of interacting with particle j (Fig. 1) is decomposed into normal and shear components:

^ ij þ F s ^sij Fij ¼ F n n

ð7Þ

where F n and F s are the magnitudes of the normal and shear components. Increments to these forces are calculated for each timestep using:

DF n ¼ K n Dun DF s ¼ K s Dus

Ei Aij jdci j

ð9Þ

where Ei is the elastic modulus of the material. The parameter ns ¼ 1=½2ð1 þ mÞ determines the ratio between the normal and shear stiffnesses, and by assuming m ¼ 0:25 for concrete, takes the value of ns ¼ 0:4 in this study [25,26]. Normal and shear deformations (Dun and Dus ) appearing in Eq. (8) are found using:

Fig. 1. Cohesive bond between a mortar matrix particle (left) and an aggregate particle (right).

dci ¼

Ki ¼

K s ¼ ns K n

j

dci

KiKj ; Ki þ Kj

ð8Þ

The graphical representation of these relationships is illustrated in Fig. 2. Note that the stiffness values K n and K s (in the normal and shear directions, respectively) are constant, regardless of the level of damage. This is in contrast to the formulation previously used by the authors [25,26], whereby both stiffnesses degrade as damage progresses. The present paper focuses on a constant-stiffness model as it gives far better results under cyclic loading compared

^ Dun ¼ Dxc;rel  n Dus ¼ Dxc;rel  ^s

ð10Þ

Moments are also developed as interacting particles move and rotate:

^ ij  Fij þ K r hij Tij ¼ r i n

ð11Þ

where the first term on the right hand side is the moment caused by the shear force defined by Eq. (7). The second term is developed by 2

the rotational stiffness K r ¼ nr K n ½minðr i ; rj Þ acting on the relative rotation hij ¼ hj  hi between particles i and j. The parameter nr is a constant ratio in this relation and is set to 0.4 for the purpose of this study. 2.3. Damage model for cohesive bonds Forces calculated by Eq. (8) must always fall into the admissible region illustrated in Fig. 3. This region is defined by:

8 > < F n 6 dt F t F n 6 dc F c > : F s 6 ds F s  li F n

ðun P 0:0Þ ðun < 0:0Þ

ð12Þ

where F t ¼ f t Aij ; F c ¼ f c Aij and F s ¼ f s Aij are limiting forces in the tensile, compressive and shear directions, and li is the coefficient of internal friction (Fig. 3). Note that f t ; f c and f s represent maximum tensile, compressive and shear stresses, respectively. The level of damage is incorporated into Eq. (12) through the parameters dt ; dc and ds which respectively correspond to tensile, compressive and shear strengths. When force components reach their limits, cohesive bonds degrade and consequently, strength values are reduced through the damage parameter d in Eq. (12). The extend of damage is governed by the maximum deformation experienced in each direction. For normal deformation in the tensile direction:

8 < dt ¼ 1:0 : dt ¼

utf udt f

ut ult

ðudt 6 ult Þ < 1:0 ðudt > ult Þ

ð13Þ

and in compression:

8 < dc ¼ 1:0

: dc ¼ 1:0 þ b



udc ulc

ðudt 6 ult Þ   1 > 1:0 ðudc > ulc Þ

ð14Þ

where ult and ulc denote linear limits for tensile and compressive deformations, respectively. Moreover, udt and udc are the maximum tensile and compressive deformations that have been experienced during the analysis, and utf is the tensile deformation limit at which the bond breaks. As shown in Fig. 2b, the value of b represents the slope of the compressive response after the initial linear region.

176

S. Sinaie et al. / Computers and Structures 196 (2018) 173–185

−Fn

Fn

Fs

Ks

Kn F¯t

δc F¯c F¯c

Kn

δt F¯t

Kn

Kn

Ks

ψ F¯s

βKn

ψδs F¯s F¯s −un

un ult

udt

uft

ulc

us

udc

uls

ufs uds

ψufs

ψuls (a)

(b)

(c)

Fig. 2. Micro-mechanical damage models under (a) tension, (b) compression and (c) shear using the constant-stiffness formulation.

2.5. Equations of motion, time integration and damping

|Fs |

Each particle in the assembly is independently governed by the equations of motion:

8 P P €i ¼ Fi  at j Fi j jxx__ i j < mi x i : Ii €hi ¼ P Ti  ar j P Ti j h_ i _ jh j

μi F¯s

ð16Þ

i

μf |Fc |

δs F¯s |Ft | δt F¯t F¯t

Fig. 3. Mohr-coulomb law for admissible stress states. Blue and red areas represent un-damaged and damaged states while the green area represents a pure friction state between two particles. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The relationship for shear damage is defined as:

8 < ds ¼ 1:0 : ds ¼

f wus uds f wus wuls

ðuds 6 wuls Þ < 1:0 ðuds > wuls Þ

ð15Þ

where w ¼ 1  li F n =F s . Note that F n can be tensile (> 0) or compressive (< 0).

P €i is where Fi is the resultant of all forces acting on particle i and x P the acceleration of the particle. Similarly, Ti is the net moment acting on particle i and € hi represents the angular acceleration of the particle. In Eq. (16), mi and Ii respectively denote the mass and the moment of inertia of particle i. The second term on each line of this equation represents the damping component, whereby at and ar are the translational and rotational damping coefficients, respectively. These are in the form of global damping, and serve the purpose of stabilizing the dynamic response of the simulations [29,30]. By monitoring the kinetic energy of the entire system and the maximum kinetic energy among all particles, it was found that setting at ¼ ar ¼ 0:1 properly dissipates kinetic energy while not affecting the results. It should be noted that increasing the damping coefficients causes fewer bonds to break under loading, resulting in higher strength values. This can be attributed to the suppression of the shock that follows the initial formation of a crack, therefore, inhibiting a realistic propagation of that crack. In this study, Eq. (16) is integrated over time using the Verlet method [31] with a constant time-step size. 3. Numerical setup

2.4. Contact interaction In contrast to cohesive bond interactions, contact interactions can appear or disappear at any time during the analysis between any two particles. Contact interactions appear whenever the distance between two particles, regardless of being of aggregate or mortar type, is less than or equal to the sum of their radii, or jdij j 6 r i þ rj . These interactions are dynamically detected at each time step. In this work, the formulation for contact interactions and cohesive bond interactions are quite similar. Both interactions have the same formulation, except that contact interactions are simplified by setting f t ¼ 0; f s ¼ 0 and f c ¼ 1, and omitting any type of damage. This can be seen as the green (darker) region depicted in Fig. 3. For contact interactions, the elastic moduli of aggregates and mortar particles are respectively set to typical values of EA ¼ 70 GPa and 24 GPa, while Poisson’s ratio is assumed to be m ¼ 0:25 for both materials.

This section describes the process of creating the DE models. This is based on the experiments carried out by Sinaie et al. [7] involving quasi-static displacement-controlled cyclic loading on normal strength concrete samples. Readers are referred to the original paper for details regarding the experimental setup, however, information required for the numerical setup will be pointed out here. The numerical setup is carried out in two steps: the first step involves the creation of particle assemblies while the second step takes care of calibrating the micro-mechanical properties. The outcome of this process is 4 individual numerical samples made up of the same material, but different in size. An important feature of this research is that material parameters are calibrated against the monotonic response alone, however, in the next section, the same parameters are used to validate the numerical samples against experimental cyclic results. This makes up a strict calibration-validation approach aiming to separate the model’s calibration requirements from its application.

S. Sinaie et al. / Computers and Structures 196 (2018) 173–185

criterion 1: High compacity (least empty space) criterion 2: Minimum particle overlaps

3.1. Creating particle assemblies The geometry of the numerical samples created here is based on the sample sizes tested by Sinaie et al. [7]. These involve cylindrical concrete samples with 38- and 63-mm diameters and aspect ratios of 1 and 2. These are denoted by DR, where the first blank field () is filled in by the diameter of the sample and the second field by its aspect ratio. Images of the cracked samples are shown in Fig. 4. The aggregate size distribution of the concrete samples was found in the laboratory and is given in Fig. 5a. Accordingly, the particle size distribution of the numerical setup is calculated and shown in Fig. 5b. All four numerical samples constitute the same particle size distribution. The aim is to create assemblies of randomly distributed aggregate particles (according to Fig. 5b) that are embedded in mortar particles. It is important to note that mortar particles are all set to be 2-mm in diameter, while aggregate particles take diameters larger than 2 mm. To ensure that the mortar particle size is small enough, the smallest sample (D38R1) was created and tested with mortar particles that were 1-mm in diameter. The resulting stress-strain curve was not very different from when 2-mm particles were used, however, having a larger number of particles resulted in a longer simulation run time. The computer code used here to carry out the sample creation process is the same code that is used later to simulate the numerical samples under cyclic loading. However, at this stage, for the purpose of sample creation, interaction between the particles is of contact type only and friction and gravity forces are disabled. Disabling these features allows a more uniform random distribution of the particles. This paper uses a ‘growing particle’ approach for sample creation. For this purpose, a boundary is first established using elastic half-space walls which span outwards. This boundary determines the permissible region for particles, and takes the form of the internal volume of a cylinder corresponding to the size of the sample under consideration. Once this permissible region is established, the following steps are carried out: 1. Aggregate and mortar particles are created and randomly dispersed within the permissible region. 2. These particles are programmed to ‘grow’ from a single point (zero-radius) into their final size over a finite duration. As illustrated in Fig. 6, during their growth, particles move around as they come into contact and push away from each other, all within the permissible region. 3. Once particles have grown to their full size, the simulation is allowed to continue running until the entire assembly reaches a state of zero kinetic energy (Fig. 6). For each sample, the number of aggregate particles is calculated using Fig. 5b. However, the exact number of mortar particles is found through trials by repeating steps 1–3. A trial is considered successful, if the following conditions are met:

D63R2.0 D63R1.0

177

The first criterion is satisfied by increasing the number of mortar particles, while the second one is satisfied by decreasing it. These criteria are indirectly monitored for the assembly by measuring the net force exerted by the particles on the top and bottom walls of the prescribed boundary. If the net force exerted by the particles on the top (or bottom wall) is less than 0.1 N at the end of a trial, a higher number of mortar particles is selected for the next run. Conversely, if the net force value is higher than 1.0 N, less mortar particles are used for the next trial. Note that a good initial estimate for these trials is to calculate the number of mortar particles according to the space that is not occupied by aggregate particles. Once these conditions are satisfied and the entire assembly reaches a state of zero kinetic energy, internal bonds are created according to Eq. 1. A value of c ¼ 0:25 is used for Eq. (1) in order for each particle to be connected to others through an average of 10–11 bonds. This concludes the sample creation process, whereby particle coordinates and their bonds are saved to file for the next stage. Using the approach described above, the 4 particle assemblies shown in Fig. 7 are created. In these images, aggregate and mortar particles are differentiated by red and blue colors, respectively. It should be mentioned that the approach taken here for sample creation is, in concept, similar to the ‘material-genesis’ method of Potyondy and Cundall [30]. However, while they satisfied criteria 1 and 2 by downscaling the size of all particles, the present paper takes a trial-error approach by changing the number of mortar particles. As a result, the method used here will allow the predetermined particle sizes of Fig. 5 to be retained.

3.2. Model calibration An important feature of this study is the aim to calibrate the model by only using a small portion of the experimental data and then using the remaining data to assess the calibrated model. As a result, out of 4 samples, only samples D38R1 and D38R2 are considered in the calibration process, and only by using their monotonic peak stress (f c ) and corresponding strain (c ) values. The remaining data, which includes the cyclic response of all 4 samples plus the monotonic response of samples D63R1 and D63R2 are used to validate the calibrated model. Making the distinction between aggregate (A) and mortar (M) particles allows separate material parameters to be set for aggregate-mortar (AM) and mortar-mortar (MM) bonds (listed in Table 1). Doing so enables the model to take into account the effect of the ITZ on the mechanical properties of aggregate-mortar bonds. Mechanical properties of the ITZ have been investigated in the literature and indicate that as a result of higher porosity, transition zones exhibit a relatively lower modulus and tensile strength in comparison to the rest of the cement paste [32,33]. In this study, AM bonds are the ones encompassing the ITZ between aggregates and the mortar matrix. Therefore, considering the values reported in the literature, the lower stiffness and strength of AM bonds (relative to MM bonds) are taken into account by:

D38R2.0 D38R1.0

Fig. 4. Images of concrete samples after being subjected to cyclic loading [7].

reduction 1: 20% on the elastic modulus (EM ) reduction 2: 30% on tensile strength (f t ) while the remaining parameters of AM and MM bonds are taken to share the same features. In order to find the material parameters of AM and MM bonds, an inverse approach is used here. This involves back to back iterations of running a simulation, comparing the results with

S. Sinaie et al. / Computers and Structures 196 (2018) 173–185

Passing (%)

100

Particles Mortar Aggregate

50

0 0.1

1

2

10

12k

Particle count per liter (n)

178

Mortar particles: d = 2 mm

10,339

n = 52 ∼ 57k 8k

4,360 4k

0

519 300 188 127 3

4

5

6

7

8

9

Particle diameter (d, mm)

Sieve size (mm)

(b)

(a)

Fig. 5. Aggregate size distribution for experimental and numerical samples. (a) Sieve analysis of experimental concrete samples. (b) Particle count per 103 m3 volume of numerical samples.

t = 0.001 s Fw = 0.0 kN Ek = 0.0 Nm

0.005 0.0 2.7(10-4 )

0.010 469.3 8.4(10-1 )

0.050 122.9 1.8(10-2 )

1000.0 0.388 3.9(10-8 )

Fig. 6. Creating a randomly dispersed 3-dimensional particle assembly from ‘growing’ particles. The variation of the force exerted on the imaginary top wall (F w ) and the total kinetic energy (Ek ) of the system over time (t) is also given.

experimental values of f c and c and then tweaking the material parameters listed in Table 1. This iterative process is continued until the simulated values of f c and c (for D38R1 and D38R2) are within 10% vicinity of the median of the experimental values. Final parametric values resulting from the calibration process is given in Table 1. Fig. 8 compares the calibrated numerical model against experimental stress-strain curves. It should be noted that the time spent to calibrate the model, makes up a large portion of the total time spent on this work. This is considered a drawback of the DEM, especially in comparison to finite element methods. As a result, the development of an efficient calibration approach has become the focus of researchers in the field [34]. GPU-based parallelized implementations of the DEM have demonstrated impressive performance gains, due to the high thread-count of graphics cards and the nature of DEM algorithms [35,36]. Here, simulations are carried out using an in-house DE computer program developed by the authors. This program implements parallel procedures based on the OpenCL standard. Using an AMD RadeonTM R7-260x GPU (an affordable, consumer-grade graphics card), simulations take 2–96 hours of run-time, based on the number of particles and the number of loading cycles. The overall computation process of the computer code is presented as a flowchart in Fig. 9. 4. Simulation results The numerical samples created using the procedure described in the last section are analyzed here. Simulations involve displacement-controlled loading cycles applied through the movement of the top and bottom loading platens. Semi-infinite Hertzian elastic bodies are use to represent these platens, whereby Young’s

modulus and Poisson’s ratio are set equal to those of steel (200 GPa and 0.3, respectively). In addition to normal forces, these halfspaces can also exert tangential shear forces to the particles that they are in contact with. These shear forces follow a static/kinetic friction model using a single coefficient of friction lf . As long as a particle stays in contact with the platen, increments of its relative tangential displacement are accumulated, whereby the tangential force is calculated and compared to the maximum friction force. However, once the contact breaks, the displacement is reset to zero. Both top and bottom platens are displaced at a constant rate of 0:5  103 ms1 . This is the case for all samples except D63R2 where a value of 1:0  103 ms1 is used to avoid a very long runtime. These rates were found to fall in a range where simulations are not sensitive to the resulting strain rate. Although these rates might seem high considering the strain-rate dependency of concrete observed in the laboratory [22], they are not high enough to have a significant effect on the numerical models of this study. This is due to the fact that the only rate-dependent parameters used here are those of global damping (at and ar ) described in Section 2.5. Loading-unloading cycles are controlled by changing the direction of motion of the top and bottom platens. Compressive loading transitions into unloading once preset strain values are reached in the sample. These preset values are multiples of the peak strain value (c ) as shown in Fig. 11. Unloading reverts to compressive loading when the force exerted by the top platen becomes less than 1.0 kN. Time increments are taken to be a constant value of Dt ¼ 107 s throughout the simulations. This time increment is small enough to avoid numerical instability. Moreover, preliminary simulations

179

S. Sinaie et al. / Computers and Structures 196 (2018) 173–185

Aggregate: 685 Mortar: 2319 AM: 7345 MM: 5526

R = 1.0

Aggregate: 3115 Mortar: 10922 AM: 34340 MM: 26438

R = H/D

DR

Aggregate: 1367 Mortar: 4644 AM: 14584 MM: 11194

R = 2.0

H

Aggregate: 6224 Mortar: 22204 AM: 68816 MM: 54327

38

63

D

Fig. 7. Dimensions of the four cylindrical samples considered in this study with their particle and bond count (AM: aggregate-mortar, MM: mortar-mortar). Red particles represent aggregates and blue particles represent the mortar matrix. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1 Micro-mechanical properties for the final calibrated model (AM: aggregate-mortar and MM: mortar-mortar).

AM MM

f t (MPa)

EA (GPa)

EM (GPa)

ns

nr

f s =f t

f c =f t

utf =ult

usf =uls

b

li

lf

5.15 7.35

70.0 NA

18.0 22.5

0.40

0.40

3.50

6.00

7.00

6.00

0.05

0.45

0.50

60

40

stress (MPa)

20 D38R1

0

D63R1 Exp. Sim.

indicated that using smaller time increments has negligible effect on the results. All damage parameters described through Eqs. (8)–(15) are set to an undamaged state at the beginning of the analysis (that is, d ¼ 1:0; k ¼ 1:0). In addition, global damping parameters at and ar of Eq. (16) are both set to 0.1. In the following text, the macro-mechanical response of the numerical samples are presented and compared to experimental results of Sinaie et al. [7]. The properties used for comparison are depicted in Fig. 10 and involve the unloading strain u , plastic strain p , reloading strain r , reloading stress f r , and reloading tangent modulus Er . These are commonly used parameters when it comes to the cyclic behavior of concrete [1,3,4]. 4.1. Qualitative properties

40

Before presenting quantitative simulation results and comparing them with experimental data, a brief examination of some qualitative properties of the model is given.

20 D38R2

0.0

0.5

1.0

D63R2

0.0

0.5

1.0

1.5

strain (%) Fig. 8. Stress-strain curves of the final calibrated model against experimental data ranges [7].

4.1.1. Calibration and validation As mentioned before, only samples D38R1 and D38R2 are used to calibrate the material parameters given in Table 1, whereby samples D63R1 and D63R2 are added for further validation of those parameters. As a result, it is expected that samples D38R1 and D38R2 show better agreement with their experimental

180

S. Sinaie et al. / Computers and Structures 196 (2018) 173–185

start read data (particles, materials)

initialize bonds

parallel

initialize timestep detect contacts: particle-wall (PW)

detect contacts: particle-particle (PP)

if bond active

calculate PW contact forces

calculate PP contact forces

calculate bond forces /check failure criteria

loop over all contacts/bonds

parallel

accumulate forces update positions/velocities loop over all particles update time loop over time steps end

fc

1.5

fr

1.0

c p

u r

strain Fig. 10. Definition of mechanical properties for cyclic stress-strain curve.

counterpart (Fig. 8). It should be mentioned that the ‘size-effect’ characteristics of DE models of concrete specimens have been studied in detail by the first author in another work [26]. The results of that work indicate that such models of concrete are capable of capturing size effects resulting form sample height, and boundary condition effects resulting from platen friction. Moreover, although size effects resulting from sample diameter are qualitatively predicted, results do not quantitatively match experimental data. These features are immediately seen in Fig. 8, and can be attributed to the fact that irregularly-shaped aggregates are simplified to spherical particles in the DEM.

normalized stress (σ/fc )

stress

Fig. 9. Flowchart for the computation process of the model.

D38R1

D63R1

D38R2

D63R2

Mon. Cyclic

0.5

0

1.0

0.5

0.0

1.0

2.0

0.0

1.0

2.0

3.0

normalized strain (/c ) Fig. 11. Simulation results for numeric samples under cyclic loading, and comparison with simulations under monotonic loading. Stress and strain values have been normalized for each individual sample using their respective f c and c values.

181

S. Sinaie et al. / Computers and Structures 196 (2018) 173–185

1.5

D38R1

1.5

D63R1

Exp. Sim.

Exp. Sim.

1.0

normalized reload stress (fr /fc )

normalized stress (σ/fc )

1.0

0.5

0 D38R2

D63R2

1.0

0.5

0.5 D38R1

0

D63R1

1.0

0.5 D38R2

0.0

1.0

2.0

0.0

1.0

2.0

3.0

0.0

1.0

2.0

D63R2

0.0

1.0

2.0

Fig. 12. Comparison of cyclic stress-strain curves obtained from simulation and those observed in experiments [7]. Stress and strain values have been normalized for each individual sample using their respective f c and c values.

Fig. 14. Comparison of experimental [7] and numerical results for the variation of reloading stress (f r ) against unloading strain (u ) for concrete samples subjected to cyclic loading.

3.0

3.0 Exp. Sim.

2.0

2.0

1.0 D38R1

0

D63R1

2.0

1.0 D38R2

1.0

2.0

normalized reload strain (r /c )

normalized plastic strain (p /c )

Exp. Sim.

0.0

3.0

normalized unload strain (u /c )

normalized strain (/c )

1.0 D38R1

0

2.0

1.0

D63R2

0.0

1.0

2.0

3.0

normalized unload strain (u /c ) Fig. 13. Comparison of experimental [7] and numerical results for the variation of plastic strain (p ) against unloading strain (u ) for concrete samples subjected to cyclic loading.

4.1.2. Monotonic and cyclic response Fig. 11 shows the simulated cyclic stress-strain curves of all 4 samples in relation to their simulated monotonic curves. As seen in this figure, the cyclic stress-strain response never exceeds the monotonic curve. In other words, the monotonic curve envelopes the cyclic curve. This property is a primary observation made from laboratory tests [1,4] and indicates an innate restriction for the reloading stress (f r ) at different levels of damage. While the monotonic curve represents an upper-bound to the cyclic response, the question of an existing lower bound still exists. Moreover, throughout the literature, different independent variables have been used as the cyclic-damage index [3,6,7,37]. However, finding the best indicator for the level of damage needs further investigation. Such studies, although possible in the laboratory, can potentially be carried out more effectively using

D63R1

D38R2

0.0

1.0

2.0

D63R2

0.0

1.0

2.0

3.0

normalized unload strain (u /c ) Fig. 15. Comparison of experimental [7] and numerical results for the variation of reloading strain (r ) against unloading strain (u ) for concrete samples subjected to cyclic loading.

numerical models. The present study creates the framework for these future investigations.

4.2. Quantitative properties Fig. 12 shows cyclic stress-strain curves obtained from simulation against experimental curves. An advantage of numerical simulations over laboratory experiments is that for the former, the same sample can be subjected to different loading schemes. As a result, for simulations, the peak strain (c ) of numerical samples is found beforehand through monotonic loading and is then used to prescribe load reversals under cyclic loading. This is not possible in the laboratory, and therefore, the unloading-reloading of numerical and experimental curves of Fig. 12 are not concurrent.

182

S. Sinaie et al. / Computers and Structures 196 (2018) 173–185

normalized reload modulus (Er × c /fc )

6.0 Exp. Sim.

4.0

2.0 D38R1

0

D63R1

4.0

2.0 D38R2

0.0

1.0

2.0

D63R2

0.0

1.0

2.0

3.0

normalized unload strain (u /c ) Fig. 16. Comparison of experimental [7] and numerical results for the variation of reloading tangent modulus (Er ) against unloading strain (u ) for concrete samples subjected to cyclic loading.

Comparing stress-strain curves can be considered circumstantial since experimental data obtained from different laboratory samples (under the same loading) always has some degree of scatter. Nevertheless, having numerical and experimental curves next to each other can give insight to their most prominent similarities and differences. In addition to that, more information can be gained when numerical and experimental results are compared in terms of damage level. Here, the variable representing damage level is chosen to be the unloading strain (u ). Figs. 13–16 show the variation of the cyclic properties defined in Fig. 10 against unloading strain. Simulation results presented in this section all belong to the single set of numerical samples created in the last section (Fig. 7), and are compared to the scatter plots resulting from experimentation. 4.2.1. Plastic strain (p ) Fig. 13 illustrates the variation of plastic strain (p ) in relation to unloading strain (u ). According to this figure, simulation results

show reasonable correspondence with experimental measurements for samples D38R1 and D38R2. However, for samples D63R1 and D63R2, plastic strain (p ) appears to be underestimated, with simulation results diverging from experiments as the unloading strain (u ) increases. Under-estimation of p for samples D63R1 and D63R2 implies that the model is over-estimating the amount of deformation recovered upon unloading. By comparing the numerical and experimental curves shown in Fig. 12, one can notice the difference between the slopes of their unloading paths. This unloading slope is higher for the experimental curves, to the extent that the unloading slope of the first cycle is even higher than the initial modulus. This variation of the unloading slope can be attributed to the compaction that a concrete sample experiences under compression. Furthermore, the unloading paths of experimental curves experience much less non-linearity in comparison to the unloading paths of simulated curves. This also contributes to the discrepancy of the plastic strain (p ) between simulation and experiments. It is quite difficult to give a detailed explanation for the degree of nonlinearity exhibited in the stressstrain response during unloading (Fig. 11), and not much information is given in the literature. Including a stiffness hardening behavior into the micro-mechanical model or modifying the contact model to allow for compaction can potentially reduce the gap between experimental and numerical values of p . 4.2.2. Reloading stress (f r ) Fig. 14 shows the variation of reloading stress (f r ) against unloading strain (u ). For f r , the correspondence between simulation and experimentation is quite well for all samples. This indicates that at different levels of u , the damage induced into concrete samples is properly represented by the numerical model. At the same time, by viewing reloading stress (f r ) as an indicator of residual strength, the correlation between damage-level and residual-strength is shown to be reasonably predicted by the model. 4.2.3. Reloading strain (r ) The variation of reloading strain (r ) against unloading strain (u ) is shown in Fig. 15. Similar to reloading stress (f r ), the agreement between simulation and experimentation exists for all of the samples. With f r and r reflecting the main residual properties of a damaged concrete specimen, the reliability of DE models to

/c

3.0

D38R1

2.0 1.0 c = 0.008947

0.0

σ/fc

1

2

fc = 53.5 MPa

3

1.0

4 5

0.5

6

0.0

n/n0

n0 = 12, 871

1.0 0.5 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

pseudo-time Fig. 17. For sample D38R1, variation of strain, stress and active (unbroken) bond against pseudo-time. The labels (1–6) shown over the stress-time curve indicate moments of peak stress and correspond to the force-chain snapshots illustrated in Fig. 19.

183

S. Sinaie et al. / Computers and Structures 196 (2018) 173–185

/c

3.0

D38R2

2.0 1.0 c = 0.004211

0.0

σ/fc

1

2

fc = 45.3 MPa

3

1.0 4 5

0.5 0.0

n/n0

n0 = 25, 778

1.0 0.5 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

pseudo-time Fig. 18. For sample D38R2, variation of strain, stress and active (unbroken) bond against pseudo-time. The labels (1–5) shown over the stress-time curve indicate moments of peak stress and correspond to the force-chain snapshots illustrated in Fig. 19.

10 mm

1

2

3

4

5

6

Fig. 19. Development of internal force-chains over the cyclic loading history. Snapshots were taken at moments of peak-stress. Compression and tensile forces are shown in blue and red, respectively, while gray lines represent bonds that have been broken. Labels shown over the samples correspond to the same ones given in Figs. 17 and 18. In order to avoid an otherwise cluttered image, these force-chains only include a longitudinal 10-mm strip of each sample. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

reflect the correlation between damage and strength is demonstrated through these results. Moreover, it should be reminded that the DE model was only calibrated against monotonic properties, and the fact that it is predicting cyclic properties so well, is evidence to the soundness of its underlying principles.

experimental data. This implies potential application of DE simulations for studies involving non-destructive testing using the dynamic modulus of elasticity [38] and health monitoring techniques where the elastic modulus plays an important role in the propagation of acoustic emissions [39].

4.2.4. Reloading modulus (Er ) Fig. 16 illustrates the variation of reloading modulus (Er ) against unloading strain (u ). It is observed that for all samples, there is good correspondence between simulation results and

4.3. Failure-progression of cohesive-bonds Figs. 17 and 18 plot stress and strain variations of samples D38R1 and D38R2 (from Fig. 11) against a pseudo-time variable

184

S. Sinaie et al. / Computers and Structures 196 (2018) 173–185

−Fn

Fn

Fs

Ks

Kn F¯t λt Kn δt F¯t

δc F¯c F¯c

βKn λc Kn

Kn

ψ F¯s

λs Ks

ψδs F¯s F¯s −un

un ult

udt

uft

ulc

(a)

udc

(b)

us uls

ufs ψuls

uds

ψufs

(c)

Fig. 20. Degrading-stiffness micro-mechanical damage models under (a) tension, (b) compression and (c) shear.

through the sample and the load-bearing region has taken an hourglass shape. This hour-glass shape can be seen in the force-chain snapshots shown in Fig. 19.

normalized stress (σ/fc )

D38R1

4.4. Degrading-stiffness micro-mechanical model

D38R2

1.0

0.5

0.0

1.0

2.0

normalized strain (/c ) Fig. 21. Simulation results for numeric samples under cyclic loading using degradable-stiffness micro-mechanical models.

The results presented so far have all been attained using the constant-stiffness micro-mechanical model described in Section 2. This section aims to show the result of allowing the stiffness of the bond to degrade as a function of damage. Such a model is depicted in Fig. 20 where the damage parameters kt ; kc and ks account for the degradation of stiffness in tension, compression and shear, respectively. Similar to other damage parameters, these are typically defined in terms of the maximum deformation in their respective direction. In order to illustrate the shortcoming of using the degradingstiffness models of Fig. 20, samples D38R1 and D38R2 are simulated under cyclic loading using the same material properties given in Table 1. The result of these simulations are shown in Fig. 21. These results clearly demonstrate the poor performance of degrading-stiffness models in relation to load reversals (compare to Fig. 11).

5. Conclusions (t). Note that the strain rate is constant in terms of this time variable, whereby positive and negative slopes respectively correspond to loading and unloading phases of the cyclic history. These figures also show the number of unbroken (active) cohesive bonds over the same time span. It is observed that the macro-mechanical compressive strength of a numerical sample is directly related to the number of unbroken bonds of the sample (n). In other words, the number of unbroken bonds is an indication of residual strength. However, the relation is not linear. Indicated by the negative slope of the n  t curve, these figures show that cohesive-bonds fail during both loading and unloading phases. However, as expected, the rate of failure is higher during loading. These plots also show that the rate of bond-failure is not the same for all cycles, and the rate tends to decline as the cycles progress and the residual strength decreases. Fig. 19 shows snapshots of the force-chains developed inside the sample. These snapshots are taken at the same moment when the stress reaches its maximum value during each loading cycle. The labels 1–6 in Fig. 19 correspond to the same labels depicted in Figs. 17 and 18. Illustrated through Figs. 11, 17 and 18, by the time the loading sequence enters the 4th cycle, residual strength f r is almost at 0:5f c . At this point, major cracks have propagated

This study evaluated the performance of a DE model of concrete in the context of cyclic loading. The proposed numerical model took advantage of a multi-phase implementation, allowing aggregates and the mortar matrix to be differentiated using different types of particles. A set of four different numerical samples with different sizes were created based on the concrete specimens tested by Sinaie et al. [7]. A single set of material parameters was determined for the model using the monotonic results of only two (out of four) of the laboratory specimens, allowing all the cyclic data (for all four specimens) to be used for further validation of the numerical model. The reason behind using a small subset of the experimental data for calibration was to avoid making a strenuous calibration process a requirement for the outcomes. With the micro-mechanical material model being rateindependent, to ensure that simulations represented a quasistatic condition similar to the experiments, loading rates were kept small while a global damping model was employed to dissipate kinetic energy. Monotonic and cyclic characteristics of the numerical samples were obtained under displacement-controlled loading and compared to their experimental counterparts. The following conclusions can be drawn from these comparisons:

S. Sinaie et al. / Computers and Structures 196 (2018) 173–185

1. A DE model is capable of reproducing the cyclic characteristics of concrete samples observed in the laboratory, both qualitatively and quantitatively (Figs. 13–16). 2. The model performs well under cyclic loading, even though the material parameters were only calibrated against a monotonic stress-strain response. 3. The variation of the reloading strain, stress and modulus (r ; f r and M r ) with unloading strain is reasonably predicted by the DE model. This indicates that such models can be reliably used to predict residual strength in terms of damage. 4. Combined with the visualization of internal forces, the number of unbroken inter-particle bonds gives a good indication of the residual load bearing capacity of a sample (Figs. 17–19). Successfully validating the capability of the DEM for the cyclic response of concrete specimens lays the foundation for assessing the cyclic-damage response of concrete under more complex loading schemes.

Acknowledgment The financial support from the Australian Research Council (ARC) via Discovery project DP170100851 (Tuan Duc Ngo) and DECRA project DE160100577 (Vinh Phu Nguyen) is gratefully acknowledged. The first author also wishes to express his gratitude to the Civil Engineering Department at Monash University where he started working on the discrete element method.

References [1] Bahn BY, Hsu CTT. Stress-strain behavior of concrete under cyclic loading. ACI Mater J 1998;95(2). [2] Sakai J, Kawashima K. Unloading and reloading stress–strain model for confined concrete. J Struct Eng 2006;132(1):112–22. [3] Sima JF, Roca P, Molins C. Cyclic constitutive model for concrete. Eng Struct 2008;30(3):695–706. [4] Osorio E, Bairán JM, Marí AR. Lateral behavior of concrete under uniaxial compressive cyclic loading. Mater Struct 2013;46(5):709–24. [5] Fathi H, Farhang K. Effect of cyclic loadings on heated self-compacting concrete. Constr Build Mater 2014;69:26–31. [6] Breccolotti M, Bonfigli MF, D’Alessandro A, Materazzi AL. Constitutive modeling of plain concrete subjected to cyclic uniaxial compressive loading. Constr Build Mater 2015;94:172–80. [7] Sinaie S, Heidarpour A, Zhao XL, Sanjayan JG. Effect of size on the response of cylindrical concrete samples under cyclic loading. Constr Build Mater 2015;84:399–408. [8] Sinaie S, Heidarpour A, Zhao X. Effect of pre-induced cyclic damage on the mechanical properties of concrete exposed to elevated temperatures. Constr Build Mater 2016;112:867–76. [9] Caratelli A, Meda A, Rinaldi Z. Monotonic and cyclic behaviour of lightweight concrete beams with and without steel fiber reinforcement. Constr Build Mater 2016;122:23–35. [10] Zhou F, Xu W. Cyclic loading tests on concrete-filled double-skin (SHS outer and CHS inner) stainless steel tubular beam-columns. Eng Struct 2016;127:304–18. [11] Rabczuk T, Belytschko T. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Comput Methods Appl Mech Eng 2007;196(29):2777–99. [12] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A geometrically non-linear threedimensional cohesive crack method for reinforced concrete structures. Eng Fract Mech 2008;75(16):4740–58.

185

[13] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A simple and robust threedimensional cracking-particle method without enrichment. Comput Methods Appl Mech Eng 2010;199(37):2437–55. [14] Wu K, Yang D, Wright N. A coupled SPH-DEM model for fluid-structure interaction problems with free-surface flow and structural failure. Comput Struct 2016;177:141–61. [15] Sobieski W, Zhang Q. Multi-scale modeling of flow resistance in granular porous media. Math Comput Simul 2017;132:159–71. [16] Cundall PA, Strack OD. A discrete numerical model for granular assemblies. geotechnique 1979;29(1):47–65. [17] Ergenzinger C, Seifried R, Eberhard P. A discrete element model predicting the strength of ballast stones. Comput Struct 2012;108:3–13. [18] Guo L, Latham JP, Xiang J. Numerical simulation of breakages of concrete armour units using a three-dimensional fracture model in the context of the combined finite-discrete element method. Comput Struct 2015;146:117–42. [19] Chen X, Chan AH, Yang J. Simulating the breakage of glass under hard body impact using the combined finite-discrete element method. Comput Struct 2016;177:56–68. [20] Nguyen TT, Bui HH, Ngo TD, Nguyen GD. Experimental and numerical investigation of influence of air-voids on the compressive behaviour of foamed concrete. Mater Des 2017;130:103–19. [21] Hentz S, Daudeville L, Donzé FV. Identification and validation of a discrete element model for concrete. J Eng Mech 2004;130(6):709–19. [22] Hentz S, Donzé FV, Daudeville L. Discrete element modelling of concrete submitted to dynamic loading at high strain rates. Comput Struct 2004;82 (29):2509–24. [23] Shiu W, Donzé FV, Daudeville L. Penetration prediction of missiles with different nose shapes by the discrete element numerical approach. Comput Struct 2008;86(21):2079–86. [24] Azevedo NM, Lemos J, de Almeida JR. Influence of aggregate deformation and contact behaviour on discrete particle modelling of fracture of concrete. Eng Fract Mech 2008;75(6):1569–86. [25] Sinaie S, Heidarpour A, Zhao X. A micro-mechanical parametric study on the strength degradation of concrete due to temperature exposure using the discrete element method. Int J Solids Struct 2016;88:165–77. [26] Sinaie S. Application of the discrete element method for the simulation of size effects in concrete samples. Int J Solids Struct 2017;108:244–53. [27] Nitka M, Tejchman J. Modelling of concrete behaviour in uniaxial compression and tension with DEM. Granular Matter 2015;17(1):145–64. [28] Tran V, Donzé FV, Marin P. A discrete element model of concrete under high triaxial loading. Cem Concr Comp 2011;33(9):936–48. [29] Oñate E, Zárate F, Miquel J, Santasusana M, Celigueta MA, Arrufat F, et al. A local constitutive model for the discrete element method. Appl Geomater Concr Comput Part Mech 2015;2(2):139–60. [30] Potyondy D, Cundall P. A bonded-particle model for rock. Int. J. Rock Mech. Min. Sci. 2004;41(8):1329–64. [31] Radhi A, Behdinan K. Contemporary time integration model of atomic systems using a dynamic framework of finite element lagrangian mechanics. Comput Struct 2017;193:128–38. [32] Scrivener KL, Crumbie AK, Laugesen P. The interfacial transition zone (ITZ) between cement paste and aggregate in concrete. Interface Sci 2004;12 (4):411–21. [33] Mondal P, Shah SP, Marks LD. Nanoscale characterization of cementitious materials. ACI Mater J 2008;105(2):174–9. [34] Behraftar S, Torres SG, Scheuermann A, Williams D, Marques E, Avarzaman HJ. A calibration methodology to obtain material parameters for the representation of fracture mechanics based on discrete element simulations. Comput Geotech 2017;81:274–83. [35] Peng L, Xu J, Zhu Q, Li H, Ge W, Chen F, et al. GPU-based discrete element simulation on flow regions of flat bottomed cylindrical hopper. Powder Technol 2016;304:218–28. [36] Zheng J, An X, Huang M. GPU-based parallel algorithm for particle contact detection and its application in self-compacting concrete flow simulations. Comput Struct 2012;112:193–204. [37] Li P, Wu YF, Gravina R. Cyclic response of FRP-confined concrete with postpeak strain softening behavior. Constr Build Mater 2016;123:814–28. [38] Amini K, Jalalpour M, Delatte N. Advancing concrete strength prediction using non-destructive testing: Development and verification of a generalizable model. Constr Build Mater 2016;102:762–8. [39] Behnia A, Chai HK, Shiotani T. Advanced structural health monitoring of concrete structures with the aid of acoustic emission. Constr Build Mater 2014;65:282–302.