A discrete element implementation for concrete: particle generation, contact calculations, packing under gravity and modeling material response

Abstract

A three-dimensional discrete element modelling capability for concrete based on rigid particles of arbitrary shape and size has been developed. The novel particle generation algorithm allows control of particle size, angularity and flakiness. General rigid body kinematics including finite rotations is accounted for, and an explicit time integration algorithm that conserves energy and momentum is implemented. An efficient contact algorithm with several features to increase the efficiency of the contact computations has been developed. This enables the gravity packing problem for arbitrary shaped particles to be solved in reasonable run time. The proposed procedure is used to generate assemblies of concrete specimens of various sizes that are homogeneous and isotropic in the bulk, and can capture the wall effect due to the formwork. The calibrated specimens are seen to be capable of accurately capturing experimentally observed macro stress strain response and failure patterns. The influence of aggregate shape on texture formation in the packed specimen, and on macro strength and failure patterns in hardened concrete, is demonstrated and is seen to be consistent with experimental results.

Appendices

Appendix A: Contact algorithm: details of modifications

1.1 A.1 Finding the intersection between the valid regions for any two vertices in 3D

The valid region for vertex i is formed by the ‘bounding planes’ of the vertex, each bounding plane being formed by the normals of the adjacent faces of particle A that meet at vertex i. These normals thus comprise the ‘bounding vectors’ of the bounding planes. The valid region for vertex j can be similarly determined from the normals of the faces of particle B that meet at vertex j (Fig. 29a, b). The intersection of the valid regions of vertex i belonging to particle A and vertex j belonging to particle B may just be the valid region of vertex i if the valid region of vertex i is contained entirely within the valid region of vertex j, or vice versa. This is the case illustrated in the top figure in Fig. 29c.

Another possibility is that the two valid regions actually intersect (bottom figure in Fig. 29c). In such a case, the lines of intersection between the bounding planes defining the valid region of vertex i and the bounding planes defining the valid region of vertex j have to be found. For an intersection to be valid, the line of intersection has to lie within the valid region of both the planes (Fig. 30). Once the vector defining the line of intersection has been determined, this vector replaces one of the two bounding vectors for both the bounding planes being considered. Thus when a bounding plane of j finds valid intersections with a bounding plane of i, it has to be determined which of the two bounding vectors defining the bounding plane of j will remain within the intersection region and which of the two bounding vectors has to be discarded. Simple geometric tests were developed to do this.

Fig. 29

Intersection region between the valid regions of vertex i, particle A and vertex j, particle B. a Valid region of A at i, b valid region of B at j, c intersection between valid regions at i at j

Fig. 30

Lines in darker outline in b show boundaries of valid region obtained using line of intersection and retained bounding vectors of planes ‘a’ and ‘b’

1.2 A.2 Kinetic energy based local tracking algorithm

A kinetic energy based algorithm was developed to determine the optimum interval $t_{2}$–$t_{1}$ between global trackings. It involves performing detailed checks, as described in the flow chart in Fig. 31, to determine the optimum size of the interval $[t_{1}, t_{2}]$ independently for each particle involved in a contact interaction. Once this has been done for both particles, the minimum of the two independently determined optimum intervals is taken to be the actual time interval between successive global trackings for that contact.

Since the flow chart is self-explanatory only basic features of the algorithm are discussed here. If the total kinetic energy of the system at time t is substantially less than the peak value of kinetic energy attained over the entire loading history up to time t, then the system as a whole is close to a state of static equilibrium. In such a situation, in order to determine whether the configuration of a particular particle is changing sufficiently slowly, one only needs to examine the history of the kinetic energy of that particle (Conditions I and II in the flow chart). On the other hand when the total kinetic energy of the system is significant, one has to compare the kinetic energy of a particular particle with the kinetic energy of the fastest moving particle in the system in order to determine whether that particle is truly a slow moving particle (Conditions III and IV in the flow chart). Examining the time history of the kinetic energy of individual particles is then not sufficient, since the kinetic energy of a particle may have reduced substantially from its peak kinetic energy, but compared to the other particles in the system, it may still be a relatively fast moving particle. In addition, special provision is made for interactions between particle pairs that are ‘dormant’. Dormant interactions are interactions where two particles are sufficiently close to each other to require that their configurations be tracked to determine whether they are coming in contact. However, their relative velocities are so small that it is unlikely that this would happen at any proximate time. Many such interactions are found to exist during the packing phase of the computations, and monitoring such interactions at every increment to determine if they are near to establishing contact is a significant computational burden. This can be reduced if the dormant interactions are identified. If the kinetic energies of both particles in an interaction are very low compared to the kinetic energy of the fastest moving particle in the system, and the gap between the particles is significant, then such particle pairs are flagged as dormant and contact detection computations are only performed after a significantly large number $(n_{{\textit{dormant}}})$ of time increments. The kinetic energy of each dormant pair is however monitored every increment; if there is a substantial increase in the kinetic energy of either particle in the dormant pair then the pair is reactivated i.e. contact detection computations are performed every increment.

The default values of the parameters used in the algorithm are given in Fig. 31. These values were found to be robust for initial particle impact velocities in the range of 0.25-3 m/s and particle masses in the range of $1\times 10^{-5}$ to $1\times 10^{-1}$ kg, which were thought to be typical of the concrete casting (packing) process. However the parameter values are independent of the logical structure of the algorithm, and may be modified if required.

Fig. 31

Kinetic energy based algorithm to find intervals between global tracking

Fig. 32

Contact point calculation. a Multiple penetrations at nodes i, j, k and l, b increased penetration at nodes k and l due to penetration reduction at i and j, c single contact point at centroid of all penetrated nodes

1.3 A.3 Determining a single contact point between contacting particles

Once the common plane has been determined, all vertices of the slave particle that are within a distance $d_{{\textit{gap}}}$ (Sect. 3.3) of the common plane are found, and are treated as vertices that have positive penetration. One possibility is to enforce contact at the vertex of the slave particle A, say vertex i, that has the largest penetration. This works reasonably well when there are no other vertices of particle A that have positive penetration. However in case there are additional vertices of particle A that have also penetrated, but have penetrations that are smaller than that of vertex i, this leads to undesirable behavior. This can be seen from Fig. 32, where four vertices of slave particle A, say vertices i, j, k and l have positive penetrations. Vertices i and j have penetrated by the same amount and their penetrations are larger than that of vertices k and l (Fig. 32a). If the contact point is taken to be the centroid of vertices i and j, and the contact force on particle A applied at this location, particle A rotates, resulting in increased penetration at vertices k and l (Fig. 32b). Since the objective is to reduce the penetration of all the vertices i, j, k and l, this is not desirable. To prevent this from happening, the contact point is chosen to be the centroid of all vertices of slave particle A that have positive penetrations i.e. the centroid of the vertices i, j, k and l (Fig. 32b). This ensures that all these vertices experience a translational force pushing them apart from the common plane.

To understand the differences with [11] one may consider the situation depicted in Fig. 32a. Chang and Chen [11] would treat this case as an ‘edge’ contact between slave edge ij and the master. Equal contact forces would act at the vertices i and j pushing the two particles apart. In the current algorithm however, a single contact force would act at the centroid (Fig. 32c). If, as in Fig. 32c, the penetration of the centroid is smaller than that of the vertices i and j, the magnitude of the contact force computed using the current algorithm would be somewhat smaller. Hence it may take a larger number of increments to resolve the penetration of edge ij. However, unlike [11], the current algorithm would not result in a reduction in the gap between the common plane and the vertices k and l (Fig. 32b).

While the above comparison is for a particular situation, the current approach seems to work well overall. There is no increase in penetration in the long run as compared to [11] (cf. Fig.13), and the savings in computational time are enormous.

Appendix B: Rotary scaling factor for very flaky/elongated particles

The following procedure ensures that flaky particles in the assembly, with very small values of rotary inertia about their axes of maximum elongation, do not give rise to very low values of $\Delta t_{{\textit{stable}}}$.

The largest eigen value for the undamped interactions of each two particle system is found by solving the eigen problem:

$$\begin{aligned}&{[}\tilde{{{\varvec{M}}}} -\omega ^{2} {\tilde{{{\varvec{K}}}}}{]} {\varvec{\uppsi }}=\mathbf{0 }\nonumber \\&{\tilde{{{\varvec{M}}}}} =\left[ {{\begin{array}{cccc} \mathbf{M}_\mathrm{A} &{} \mathbf{0 }&{} \mathbf{0 }&{} \mathbf{0 } \\ \mathbf{0 }&{} \mathbf{M}_\mathrm{B}&{} \mathbf{0 }&{} \mathbf{0 } \\ \mathbf{0 }&{} \mathbf{0 }&{} \mathbf{I}_\mathrm{A}&{} \mathbf{0 } \\ \mathbf{0 }&{} \mathbf{0 }&{} \mathbf{0 }&{} \mathbf{I}_\mathrm{B} \\ \end{array}}}\right] {\tilde{{{\varvec{K}}}}} =\left[ {{\begin{array}{cccc} \mathbf{K}&{} {-\mathbf{K}}&{} {\mathbf{K}{\hat{\mathbf{r}}}_\mathrm{A}^\mathrm{T} }&{} -\mathbf{K}{\hat{\mathbf{r}}}_\mathrm{B}^\mathrm{T} \\ {-\mathbf{K}}&{} \mathbf{K}&{} {-\mathbf{K} {\hat{\mathbf{r}}}_\mathrm{A}^\mathrm{T} }&{} \mathbf{K}{\hat{\mathbf{r}}} _\mathrm{B}^\mathrm{T} \\ \mathbf{K}{\hat{\mathbf{r}}}_\mathrm{A} &{} -\mathbf{K}{\hat{\mathbf{r}}} _\mathrm{A} &{} {\hat{\mathbf{r}}}_\mathrm{A} \mathbf{K}{\hat{\mathbf{r}}}_\mathrm{A}^\mathrm{T} &{} -{\hat{\mathbf{r}}} _\mathrm{A} \mathbf{K}{\hat{\mathbf{r}}}_\mathrm{B}^\mathrm{T} \\ -\mathbf{K}{\hat{\mathbf{r}}} _\mathrm{B} &{} \mathbf{K}{\hat{\mathbf{r}}} _\mathrm{B} &{} \hbox {}-{\hat{\mathbf{r}}} _\mathrm{B} \mathbf{K}{\hat{\mathbf{r}}} _\mathrm{A}^\mathrm{T} &{} {\hat{\mathbf{r}}} _\mathrm{B} \mathbf{K}{\hat{\mathbf{r}}}_\mathrm{B}^\mathrm{T} \\ \end{array}}} \right] .\nonumber \\ \end{aligned}$$

(B.1)

$\tilde{{{\varvec{M}}}}$ and $\tilde{{{\varvec{K}}}}$ are the generalized mass and stiffness matrices of dimension $12\times 12$ and ${\varvec{\uppsi }}$ and $\omega$ are the eigenvectors and eigen values. The generalized mass matrix $\tilde{{{\varvec{M}}}}$ contains contributions from the mass as well as the rotary inertias of the two particles, A and B, and is evaluated at their respective centers of mass. $\mathbf{M}_{\mathrm{A}}$ and $\mathbf{M}_{\mathrm{B}}$ are $3\times 3$ diagonal matrices. $\mathbf{M}_{\mathrm{A}}$ is the mass matrix of particle A, with mass $M_{A}$ and $\mathbf{M}_{\mathrm{B}}$ is the mass matrix of particle B with mass $M_{B}$. $\mathbf{I}_{\mathrm{A}}$ and $\mathbf{I}_{\mathrm{B}}$ are the $3\times 3$ rotary inertia matrices for particles A and B. The generalized stiffness matrix $\tilde{{{\varvec{K}}}}$ involves $3\times 3$ skew symmetric matrices ${\hat{\mathbf{r}}}_{\mathrm{A}}$ and ${\hat{\mathbf{r}}}_\mathrm{B}$, whose axial vectors are the vectors $\mathbf{r}_{\mathrm{A}}$ and $\mathbf{r}_{\mathrm{B}}$ (Fig. 33) as well as the diagonal matrix K whose diagonal term is the normal spring stiffness $K_{n}$.

Fig. 33

Spring-mass-rotary inertia System for two interacting particles

The equation governing the free vibration of the spring-mass-rotary inertia system with only translational degrees of freedom active is given by:

$$\begin{aligned} \left[ {\begin{array}{cc} \mathbf{M}_\mathrm{A} &{} 0 \\ 0&{} \mathbf{M}_\mathrm{B} \\ \end{array}} \right] \left[ {{\begin{array}{cc} {\ddot{\mathbf{x}}}_\mathrm{A} \\ {\ddot{\mathbf{x}}}_\mathrm{B} \\ \end{array}}}\right] +\left[ {{\begin{array}{cc} \mathbf{K}&{}-\mathbf{K}\\ -\mathbf{K}&{}\mathbf{K}\\ \end{array}}}\right] \left[ {{\begin{array}{cc} \mathbf{X}_\mathrm{A}\\ \mathbf{X}_\mathrm{B}\\ \end{array}}}\right] =\mathbf{0 } \end{aligned}$$

(B.2)

where $\mathbf{x}_{\mathrm{A}}$ and $\mathbf{x}_{\mathrm{B}}$ are the spatial degrees of freedom at the center of mass of particles A and B respectively. The eigen value problem corresponding to Eq. (B.2) is solved for all interacting particles to obtain an estimate of the stable time increment of the translation-only system. If this time increment is denoted as $\Delta t_{stable}^{trans-only}$, the goal is to scale the rotary inertias $\mathbf{I}_\mathrm{A}$ and $\mathbf{I}_\mathrm{B}$ such that the scaled rotary inertias $\mathbf{I}_\mathrm{A}^{scaled}$ and $\mathbf{I}_\mathrm{B}^{scaled}$ are such that the largest eigen value of the scaled system, $\omega _{max}^{scaled}$, is less than or equal to the largest eigen value of the translation-only system, $\omega _{\max }^{{\textit{trans}}-{\textit{only}}}$, i.e.

$$\begin{aligned} \omega _{max}^{scaled} \le \omega _{max}^{trans-only} =2/{\Delta }t_{stable}^{trans-only} \end{aligned}$$

(B.3)

In order to obtain a usable bound on the rotary inertias using the above condition, the following conservative simplifications are made. ${\hat{\mathbf{r}}}_{\mathrm{A}}$ and ${\hat{\mathbf{r}}}_\mathrm{B}$ in $\tilde{{{\varvec{K}}}}\hbox { are}$ replaced by the following upper bound matrix r:

$$\begin{aligned} \mathbf{r}=\hbox {}\left[ {{\begin{array}{ccc} 0&{}\quad {-d}&{}\quad d \\ d&{}\quad 0\quad &{}\quad {-d} \\ -d&{}\quad d&{}\quad 0 \\ \end{array}}} \right] ;\Vert \mathbf{r}\Vert _\infty \ge \hbox {max}\left( \Vert {\hat{\mathbf{r}}} _\mathrm{A}\Vert _\infty ,\,\,\Vert {\hat{\mathbf{r}}}_\mathrm{B}\Vert _\infty \right) \nonumber \\ \end{aligned}$$

(B.4)

and $\mathbf{I}_\mathrm{A}$ and $\mathbf{I}_\mathrm{B}$ in $\tilde{{{\varvec{M}}}}$ are replaced by the following lower bound:

$$\begin{aligned} \mathbf{I}=\left[ {{\begin{array}{lll} I&{}\quad 0&{}\quad 0 \\ 0&{}\quad I&{}\quad 0 \\ 0&{}\quad 0&{}\quad I \\ \end{array} }} \right] ;\quad \Vert \mathbf{I}\Vert _\infty \le \mathbf{min} \left( \Vert {{\varvec{I}}}_{{\varvec{A}}}\Vert _\infty ,\Vert {{\varvec{I}}}_{{\varvec{B}}}\Vert _\infty \right) \end{aligned}$$

(B.5)

In Eq. (B.4) d is the effective radius of the larger of the two particles A and B, and in Eq. (B.5) the isotropic rotary inertia, I, is defined such that its scalar component $I\hbox { is equal to }{\textit{min}}\left( {{\varvec{\omega }}_{\mathbf{I}_\mathrm{A} } , {\varvec{\omega }}_{\mathbf{I}_\mathrm{B}}}\right)$ where ${\varvec{\omega }}_{\mathbf{I}_\mathrm{A}}$ and ${\varvec{\omega }}_{\mathbf{I}_\mathrm{B}}$ represent the spectrum of $\mathbf{I}_\mathrm{A}$ and $\mathbf{I}_\mathrm{B}$ respectively. $\omega _{max}^{scaled}$ can then be calculated by solving the eigen value problem given by Eq. (B.1) with $\tilde{{{\varvec{K}}}}$ and $\tilde{{{\varvec{M}}}}$ modified as above:

$$\begin{aligned} \omega _{max}^{scaled} =\sqrt{\frac{K_n (M_A +M_B )}{M_A M_B }+\frac{6K_n d^{2}}{I}} \end{aligned}$$

(B.6)

The corresponding bound on I is then obtained by enforcing Eq. (B.3):

$$\begin{aligned} I\ge 6K_n d^{2}\big /\left[ 4/\left( {\varDelta t_{stable}^{trans-only}} \right) ^{2}-\hbox {}K_n (M_A +M_B )/M_A M_B \right] \nonumber \\ \end{aligned}$$

(B.7)

The rotary inertias $\mathbf{I}_\mathrm{A}$ and $\mathbf{I}_\mathrm{B}$ are scaled to ensure that their smallest eigen values satisfy this lower bound, thereby ensuring that the stability of the system is governed by $\Delta t_{stable}^{trans-only}$.

Inclusion of viscous damping leads to a reduction in the value of $\Delta t_{stable}^{trans-only}$ since the largest frequency of the constrained spring mass system described by Eq. (B.2) is affected by the presence of the viscous damping term $\mathbf{C}\left[ {{\begin{array}{l} {\dot{\mathbf{x}}}_\mathrm{A} \\ {\dot{\mathbf{x}}}_\mathrm{B} \\ \end{array}}}\right]$. The reduced stable time increment for the damped system is then given by:

$$\begin{aligned} \left( \Delta t_{stable}^{trans-only}\right) _{damped} =\Delta t_{stable}^{trans-only} \left( {\sqrt{1+\xi ^{2}}-\xi }\right) \nonumber \\ \end{aligned}$$

(B.8)

where $\xi$ is the ratio of the damping coefficient to its critical value for the spring-mass-damper system with only translational degrees of freedom active. $(\Delta t_{stable}^{trans-only} )_{damped}$ is then used instead of $\Delta t_{stable}^{trans-only}$ in Eq. (B.7).

Appendix C: Calibration of meso elastic parameters for hardened concrete

The meso moduli $K_n^*\hbox { and }K_s^*$ must be calibrated to replicate the elastic response at the macro level, determined by the macro deformability parameters, E (Young’s modulus) and $\upsilon$ (Poisson’s ratio). The functions $E\left( {K_n^*,K_s^*}\right)$ and $\nu \left( {K_n^*,K_s^*}\right)$ are not known for arbitrary polyhedral particles. Trying to obtain them through regression analysis from experimental data requires performing, for every DEM assembly, a large number of experiments to obtain reasonable fits. An alternative calibration procedure was therefore developed. The calibration procedure starts with assuming values of the meso-level material parameters $K_n^{*} \hbox { and }K_s^*$, say $(K_n^*)^{0}\hbox { and }(K_s^*)^{0}$. Two simple load tests that are independent of each other, e.g. uniaxial compression and triaxial compression, are then performed on the DEM assembly. The elastic strain energy of the DEM assembly $W_{DEM}^I ,I=1,2$ being the number of independent tests, is then calculated using the assumed values $(K_n^*)^{0}\hbox { and }(K_s^*)^{0}$, and the normal and shear contact forces $F_k^n$ and $F_k^s$ at each contact interaction k:

$$\begin{aligned} 2W_{DEM}^I= & {} \left( {1/(K_n^{*} )^{0}} \right) \mathop \sum \nolimits _{k=1}^N \left[ {(F_k^n )^{2}/A_{\mathrm{int}}^k } \right] \nonumber \\&+\left( {1/(K_s^{*} )^{0}} \right) \mathop \sum \nolimits _{k=1}^N \left[ {(F_k^s )^{2}/A_{\mathrm{int}}^k } \right] \end{aligned}$$

(C.1)

It is assumed in Eq. (C.1) that there are N contact interactions, each with area of interaction equal to $A_{\mathrm{int}}^k$. After performing the summations Eq. (C.1) can be written more compactly as:

$$\begin{aligned} W_{DEM}^I =a_I /(K_n^{*})^{0}+b_I /(K_s^{*})^{0} \end{aligned}$$

(C.2)

where:

$$\begin{aligned} a_I =\mathop \sum \nolimits _{k=1}^N (F_k^n )^{2}/A_{\mathrm{int}}^k \hbox { and }b_I =\mathop \sum \nolimits _{k=1}^N (F_k^s )^{2}/A_{\mathrm{int}}^k \end{aligned}$$

(C.3)

The strain energies of the corresponding macro experiments, $W_{Macro}^1\hbox { and }W_{Macro}^2$, are also found in terms of the components of the macro stress tensor ${\varvec{\sigma }}$ and macro modulus E and Poisson’s ratio $\nu$. If $I = 1$ corresponds to a uniaxial test, with the uniaxial stress $\sigma _{zz}$ applied in the z direction,

$$\begin{aligned} W_{Macro}^1 =\sigma _{zz}^2 /2E \end{aligned}$$

(C.4)

If $I = 2$ corresponds to a triaxial test, with triaxial stress components $\sigma _{xx} ,\sigma _{yy} \hbox { and }\sigma _{zz}$ applied in the x, y and z directions,

$$\begin{aligned} W_{Macro}^2= & {} (1/2E)\left[ \sigma _{xx}^2 +\sigma _{yy}^2 +\sigma _{zz}^2 \right] \nonumber \\&-\left( {\nu /E} \right) \left[ {\sigma _{xx} \sigma _{yy} +\sigma _{yy} \sigma _{zz} +\sigma _{zz} \sigma _{xx} } \right] \end{aligned}$$

(C.5)

The goal of the calibration procedure is to correct the assumed values of $(K_n^*)^{0}\hbox { and }(K_s^*)^{0}$ to ensure that the strain energy of the DEM model equals that of the macro model, for both $I = 1, 2$. Denoting the calibration factor for $(K_n^*)^{0}$ as $c_{K_n^*}$ and the calibration factor for $(K_s^*)^{0}$ as $c_{K_s^*}$ this requires that the following two equations must be satisfied:

$$\begin{aligned}&a_1 \big /\left( {c_{K_n^{*}}(K_n^{*})^{0}} \right) +b_1 \big /\left( {c_{K_s^{*} } (K_s^{*} )^{0}} \right) =W_{macro}^1 \end{aligned}$$

(C.6)

$$\begin{aligned}&a_2 \big /\left( {c_{K_n^{*}} (K_n^{*} )^{0}} \right) +b_2 \big /\left( {c_{K_s^{*}}(K_s^{*} )^{0}} \right) =W_{macro}^2 \end{aligned}$$

(C.7)

The above two simultaneous equations are solved for $c_{K_n^{*}}$ and $c_{K_s^*}$. This yields the calibrated values $K_n^*=c_{K_n^*} (K_n^*)^{0}\hbox { and }K_s^*=c_{K_s^*} (K_s^*)^{0}$. The calibrated values of $K_n^*\hbox { and }K_{s}^*$ also depend on the average thickness of the mortar layer between the aggregates; the calibration factors $c_{K_n^*}$ and $c_{K_s^{*}}$ incorporate this dependence as well.

While ideally, only two tests are necessary to calibrate $K_n^*\hbox { and }K_{s}^*$, in practice it was observed that improved results could be obtained by performing three tests (e.g. uniaxial compression, triaxial compression and triaxial tension) instead of two. This resulted in three equations for the two unknowns $c_{K_n^{{*}}}$ and $c_{K_s^*}$. The values of $c_{K_n^\mathrm{*}}$ and $c_{K_s^*}$ were then found by least squares minimization. The stress states for the three tests are described in Table 7. It is necessary for the stress state corresponding to each test to differ significantly from the other tests, as otherwise the resulting equations may be ill-conditioned or near-singular. In such a case, it may not be possible to obtain accurate corrected values of $K_n^*\hbox { and }K_{s}^*$ from the equations.

The elastic calibration procedure was verified for assemblies which differed only in the arrangement of packed particles but had the same meso geometry and same aggregate composition and mortar mix. Twelve such assembles were considered. The calibrated values of $K_{n}^*\hbox { and }K_{s}^*$ obtained from a single (benchmark) assembly were used in all eleven remaining assemblies. All the assemblies were identically loaded in uniaxial compression. The Young’s modulus and Poisson’s ratio were calculated for each assembly by plotting the macro axial stress vs. macro axial strain, and evaluating the ratio of macro lateral strain to macro axial strain. For each assembly the E and $\upsilon$ values obtained from the DEM simulations were compared with the target values of E and $\upsilon$, denoted as $\bar{E}$ and $\bar{\upsilon }$. The $E/\bar{E}$ ratios were found to vary between 0.94 to 1.13 while the $\upsilon /\bar{\upsilon }$ ratios ranged between 0.97 and 1.04. The coefficient of variation for the $E/\bar{E}$ ratios was 4.936% while that for the $\upsilon /\bar{\upsilon }$ ratios was 2.458%. The dispersions in $E/\bar{E}$ and $\upsilon /\bar{\upsilon }$ ratios compare favorably with values reported in the literature (cf. [24]).

Table 7 Loading states used for elastic calibration

If the same cement and aggregate type as well as concrete mix proportions are used for assemblies with different meso geometries, the volume of mortar remains the same but the mortar gets distributed over very different total aggregate surface areas, since the total surface area of the aggregates varies with meso geometry. As a consequence, the average thickness of the mortar layer between the aggregates changes with meso geometry. Therefore if the calibrated values of $K_n^*\hbox { and }K_s^*$, obtained for an assembly with meso geometry ‘A’, are used in an assembly with a different meso geometry ‘B’, the macro Young’s modulus and Poisson ratio values may not match. However, according to experiments reported by previous researchers (e.g. [25, 56]), changes in meso geometry should not lead to any significant differences in the initial elastic properties of concrete. To account for this additional scalar correction factors $\gamma _{n}$ and $\gamma _{s}$, representing corrections to $K_{n}^*\hbox { and }K_s^*$ due to differences in average mortar thickness between the two assemblies, have to be determined. These factors are found on a trial and error basis. For assemblies with new meso geometries, a couple of trials were found sufficient to determine modified meso moduli $\gamma _n K_{n}^*\hbox { and }\gamma _s K_s^*$ that resulted in the $E/\bar{E}$ and $\upsilon /\bar{\upsilon }$ ratios becoming sufficiently close to 1.0 (Table 8).

Table 8 Calibrated meso elastic parameters for meso geometries with varying aggregate angularity but identical concrete mix proportions, aggregate size gradation, cement and aggregate type

Appendix D: Calibration of meso damage parameters for hardened concrete

(a) Choice of damage initiation stress: Experimental data on interfacial fracture properties of concrete is limited. Recently, Rao and Prasad [45] performed experiments to determine damage parameters for mortar-rock interfaces. The damage initiation stress values reported in this paper are used as the mean values of the respective Weibull distributions. These damage initiation stress values are not significantly different from those proposed in [61].

(b) Determination of the interface fracture energy and damage evolution parameters: No experimental data is available for the meso level interface fracture energy and hence these values are determined by performing several trials. However some guidance can be obtained from the results of macro experiments [31, 46] which suggest that Mode II fracture energy can be up to 25 to 30 times higher than the Mode I fracture energy. It is assumed that this holds true at the meso scale as well. Hence if a DEM sample, with a large Mode II to Mode I meso fracture energy ratio, is loaded under uniaxial tension, it is expected that interface failure due to tensile cracking in Mode I would be the dominant mode of failure. Several DEM simulations of cubic specimens under uniaxial tension were therefore performed, in which the Mode II fracture energy was kept fixed at a high value and the Mode I fracture energy and Mode I damage evolution parameters varied until the macro stress–strain curve obtained from the DEM simulation showed a close match with the experimental stress strain curve. The Mode I fracture energy and Mode I damage evolution parameters were then held fixed at these values. Subsequently multiple DEM simulations of uniaxial compression tests were performed, in which the Mode II fracture energy and Mode II damage evolution parameters were varied. The Mode II fracture energy and Mode II damage evolution parameters which resulted in best match between the DEM macro stress–strain curve and the experimental curve for uniaxial compression were then selected.

(c) Final check: Finally a further DEM uniaxial tension test was simulated with Mode I and Mode II fracture parameters as obtained in (b) above. If the DEM macro stress–strain curve matched the experimental uniaxial tension test curve, and the crack patterns were along expected lines, the calibration of the meso fracture parameters was regarded as complete.

In case of significant differences, the Mode II fracture parameters were held fixed at their values obtained at the end of (b), and the Mode I fracture parameters varied to obtain a good match with the expereimental uniaxial tension stress strain curve. The entire cycle described in (b) was then repeated. A consistent set of fracture parameters satisfying the final check was usually obtained after a few iterations.