# Two-Dimensional Mesoscale Finite Element Modeling of Concrete Damage and Failure
## 3 Automatic Generation of Cohesive Elements
 \tag{14}\]

where \(\sigma_{n}\) and \(\sigma_{t}\) represent normal and shear traction, respectively, and \(\delta_{n}\) and \(\delta_{s}\) are crack opening displacement in the normal and shear directions, respectively. In Equation (14), \(f\) shows the nonlinear function between the traction and crack displacement. When different fracture modes occur simultaneously on a crack tip, either by multiaxial loading of the geometry or induced by any geometry that favors multiple fracture modes under any uniaxial loading condition, then the geometry is under mixed-mode fracture states. Although the studied concrete samples in this study are loaded under the uniaxial loading condition, due to the high material heterogeneities, the cracks propagate in concrete under mixed-mode conditions rather than single-mode. Hence, its corresponding traction-separation laws are used in the finite element simulations [48]. Since bilinear TSL has a simple functional relationship and is utilized widely by many researchers in modeling heterogeneous materials such as composites, in this work, for each fracture mode, a bilinear TSL is utilized. Three model parameters must be utilized to fully characterize a bilinear TSL. These model parameters are the cohesive element penalty stiffness \(K\), the threshold stress \(\sigma_{max}\), and the critical fracture energy \(G_{C}\). Fracture energy can be interpreted as the required energy to form two new fractured surfaces [49]. The fracture energy can be calculated as follows:

\[G_{C}=\int_{0}^{\delta_{f}}\sigma\,d\delta\quad\text{for}\quad 0\leq\delta \leq\delta_{f}, \tag{15}\] \begin{table}
\begin{tabular}{l l l l l l} \hline \hline
**Parameter** & **Aggregate** & **Mortar** & **TSL-ITZ** & **TSL-Mortar** & **TSL-Aggregate** \\ \hline Maximum normal stress, \(\sigma_{n}^{0}\) (MPa) & - & - & 2.4 & 4 & - \\ Maximum shear stress, \(\sigma_{S}^{0}\) (MPa) & - & - & 10 & 30 & - \\ Normal mode fracture energy, \(G_{I}\) (N/mm) & - & - & 0.02 & 0.06 & - \\ Shear mode fracture energy, \(G_{II}\) (N/mm) & - & - & 0.4 & 1.2 & - \\ B-K criterion material parameter, \(\eta\) & - & - & 1.2 & 1.2 & - \\ \hline \hline \end{tabular}
\end{table} 
Table 2: _Cont._The cohesive element penalty stiffness is the stiffness of any cohesive element in the elastic region prior to the initiation of the damage. The penalty stiffness is a function of the elastic modulus of the interface and its thickness and is calculated as follows:

\[K=\frac{E}{t_{coh}}, \tag{16}\]

where \(t_{coh}\) is the cohesive element thickness, and \(E\) represent the elastic modulus. Accurate definition of the penalty stiffness \(K\) is crucial to appropriately predict the structural response of concrete parts since this parameter controls how easily the cohesive elements deform elastically in the interface region. A very large \(K\) would result in immediate initiation of the damage without preceding deformation when the threshold stress \(\sigma^{0}\) is reached. This would result in a large softening region where the amount of deformation (distance between \(\delta_{0}\) and \(\delta_{f}\)) leading to complete fracture would increase. On the other hand, for very small \(K\), with a given applied displacement, the model may never reach \(\sigma_{max}\) and, thus, no damage would initiate. This state in the finite element model would result in a relatively high sliding and penetration between the nodes in the interface region. Before initiation of the damage, the normal and shear stress can be calculated using the following equation:

\[ \begin{split}\sigma_{n}&=k_{n0}\langle\delta_{n} \rangle\quad\text{(normal direction)}\\ \sigma_{t}&=k_{s0}\,\delta_{s}\quad\quad\quad\text{ (shear direction)}\end{split}  \tag{17}\]

where \(k_{n0}\) and \(k_{s0}\) are the penalty stiffness in the normal and shear directions, respectively. The symbol \(<>\) used in Equation (17) represents the Macaulay bracket. The Macaulay bracket signifies that compressive deformation or stress state would not initiate damage in the model.