Reinforced Concrete Finite Element Modeling based on the Discrete Crack Approach

The behavior of reinforced concrete elements is complex due to the nature of the concrete that is weak in tension. Among these complex issues are the initial cracking and crack propagation of concrete, and the bond-slip phenomenon between the concrete and reinforcing steel. Laboratory tested specimens are not only costly, but are limited in number. Therefore a finite element analysis is favored in combination to experimental data. The finite element technique involving the cracks inserting is one of the approaches to study the behavior of reinforced concrete structures through numerical simulation. In finite element modeling, the cracks can be represented by either smeared or discrete crack. The discrete crack method has its potential to include strain discontinuity within the structure. A finite element model (FEM) including the concrete cracking and the bond-slip was developed to simulate the nonlinear response of reinforced concrete structures.


Introduction
Reinforced concrete is a composite material consisting of concrete and reinforcing steel. Concrete compression strength is very high, with a low tensile strength. Furthermore, concrete is categorized to be brittle in nature so that the use of plain concrete as a structural element becomes not feasible. The use of reinforcing steel that has a very high tensile strength and a high ductility performance in the tensile areas of concrete, will lead to a versatile material that is both excellent in terms of strength and ductility. Since concrete is weak in tension, the cracking of concrete as well as the interaction between the reinforcing steel and the cracked concrete causes highly nonlinear behavior in reinforced concrete structures [1].
Numerical analysis of reinforced concrete structures can be carried out by various methods. Among the widely used is finite element analysis. Many numerical studies using finite element method have been conducted to describe the behavior of reinforced concrete structures.
To simulate the actual behavior of reinforced concrete structures precisely, cracks needs to be inserted in the finite element model (FEM). Previous research works [2][3][4][5] modeled the cracks by either the smeared or discrete crack model. Since it was introduced by Rashid [6], the smeared crack approach in the finite element method has gained popularity, and an extensive amount of research works had been conducted based on this model. A smeared crack assumes an evenly distributed crack spreading upon concrete failure. The crack direction is presumed perpendicular to the direction of the principal tensile stress. The model also assumes the behavior of the cracked material as a continuum, the material is considered isotropic before cracking, and orthotropic prior to the development of cracks. Furthermore, the interface between concrete and reinforcing steel is assumed to be fully bonded. The advantage of the smeared crack method is its ease in comprising the algorithms into a computer program since no changes in the topology of the structure are accounted for during the entire loading process.
Contradictory to the smeared crack model, the discrete crack is modeled as a node separation on the side of adjacent elements when principle tensile stresses of nodes have reached the tensile strength of concrete. This technique is recognized to be able to represent strain discontinuity of the structure. However, the node separation has its consequence of changing the topology of the structure. Additionally, the behavior of reinforced concrete structures is also influenced by the bond-slip mechanism that occurs between concrete and reinforcing steel. The slip that occurs between the concrete and reinforcing steel will introduce the strain incompatibility near the cracks [7].
In this research, an attempt has been made to accommodate the influence of the concrete cracking and bond-slip in the reinforced concrete structure by constructing a FEM with the discrete crack approach. Complimentary experimental test specimens were prepared and tested under the exact same loading condition as the FEM to analyze the accuracy and correctness of the developed program. The mechanical properties of the concrete and reinforcing steel were obtained from laboratory tested specimens, and functioned as input to the FEM. The Visual Basic language [8] was used to generate the program algorithms.

Finite Element Modeling The Concrete Element
Concrete was modeled as a plane-stress, isoparametric quadrilateral element having two-by-two Gauss points. The stress-strain constitutive relationship of the material was based on CEB-fib Model Code [9]. The stress-strain relationship in both compression and tension was modeled separately. The stress-strain relationship in uniaxial compression is linear up to 30% of the compressive strength, prior to this point the material is nonlinear up to failure ( Figure 1). A post peak is observed, and the stiffness of the material will decrease rapidly beyond the ultimate strength [1]. The stress-strain relationship in uniaxial tension tends to respond linearly to 90% of the tensile strength, followed by a slight decrease in material stiffness [9]. The tangent stiffness method was used for modelling the nonlinearity of the material in the finite element analysis.
The material matrix was constructed using a general formulation incorporating orthotropic behavior as proposed by Chen and Saleeb [1]. At early loading stages, the material is elastic, having a constant modulus of elasticity. The developed FEM takes into consideration the effect of the cracking as a discrete crack model. The failure criteria were evaluated based on the major principal stresses at the nodes based on the Kupfer-Hilsdorf-Rusch failure envelope [10]. Once the major principal tensile stress at a nodal point exceeded the boundary of the failure envelope, a discrete crack was inserted in the FEM followed by a node separation.

The Reinforcing Steel Element
The discrete element was used for modelling the reinforcing steel. The steel bar was modeled as a one-dimensional bar element and was positioned at the interface between the corresponding isoparametric concrete elements. This discrete element was chosen because of its ability to accommodate the influence of bond-slip in the interface between the concrete and reinforcing steel with high satisfaction [11][12][13].
The stress-strain relationship of the reinforcing steel was modeled as an idealized bi-linear curve representing the elastic-linear behavior, and yielding.

The Bond Element
The bond element was modeled as a linkage element as proposed by Ngo and Scordelis [11]. The linkage element represented the stress transfer from the concrete to the steel bar as a double spring having one direction parallel and one perpendicular to the axes of the reinforcing steel ( Figure 2). The spring in the direction parallel to the axes of the reinforcing steel (Kv) represents the shear behavior or bond-slip between the concrete and reinforcing steel, while the spring in the direction perpendicular to the axes of the reinforcing steel (Kh) counts for the dowel action.

Figure 2. Linkage Element
The beam in this research was loaded in bending, and therefore dominated by the bending mode. As a consequence, the effect of dowel action is considered small and could be neglected. The value for the dowel stiffness Kh was assigned a significantly large number to the power of six (10 6 ) and therefore considered sufficiently high to diminish the effect of the dowel action [14]. The nodes at the boundary between the reinforcing steel and the concrete were connected by two nodes with the linkage element as connecting media. Initially, these double nodes have identical coordinates but different numbering, created during the meshing process. The Eligehausen model [3] was chosen to model the bond-slip relationship.
Node I in Figure 1 represents the node in the reinforcing steel, while node J originally coinciding node I, was placed in the concrete. When loading was applied, a disjunction between the two nodes in the direction of Kv occurred. The displacements of the steel node I in the global coordinate system were d1 and d2 respectively, while the displacements in the concrete node J were d3 and d4. In the local coordinate system, the axes were determined based on the position of the reinforcing steel. The angle  indicated the direction between the local and the global coordinate system.

The Discrete Crack Mechanism
An increase in loading will result in an increase in nodal stresses. When a major principal tensile stress 1 at a particular nodal point exceeds the failure envelope mandated by the Kupfer-Hilsdorf-Rusch criteria, a discrete crack is inserted in the FEM, enabling node separation. Three cases of crack propagation were considered in this research. The first case was the initiation of a discrete crack occurring in the concrete's extreme fibers in tension (Figure 3), the second was the case of crack propagation intersecting the reinforcing steel ( Figure  4), and the last case represented the propagation of diagonal cracks in the structures ( Figure 5).

Case 1
Case 1 is the initiation of a discrete crack starting from within the concrete element in the extreme fibers in tension (Figure 3). Node A was set as a starting point of the crack propagation, for this specimen the node A was located at the bottom of the beam, where the fibers are in tension. When the major principal tensile stress at node A exceeded the concrete tensile strength (1 > ftm), a node separation was performed by duplicating node A with a new node (node B) originally coinciding with node A. The node separation resulted in an increase in number of nodes in the structure. The mechanism of Case 2 is illustrated in Figure 4. This case represents the cracking of concrete at the interface. If the major principal tensile stress in the concrete node of the link in the interface has reached the failure envelope, then the crack propagated to the next node (node G) is perpendicular to the interface. This cracking mechanism resulted in a crack length of lcr in the reinforcing steel element. In the first stage of this crack propagation, the node in the concrete (node C) and the node at the reinforcing steel (node D) were duplicated, resulting in the creation of nodes E and F, coinciding with nodes C and D. The stress transfer from the concrete to the steel bar resulted in an increasing force in the reinforcing steel element. The force transferred from the concrete to the reinforcing steel was calculated using Equation (1), while the bond stress which occured in the reinforcing steel element was obtained by Equation (2). Slip occurred opposite to the right and left side of the cracks due to the pull out effect. This sequence was followed by the formation of a debonded length (lel). The bond-slip phenomenon was determined through the bond-slip constitutive model from Eligehausen [3]. The dragging of the node is a representation of the change in the crack tip. In between the crack tip, a new reinforcing steel element with a length of lel was added to the model, adjacent to the original reinforcing steel element. This additional reinforcement steel element had the exact same material properties as the original steel bar.
where: fC = the major principal tensile stress at node C fG = the major principal tensile stress at node G lel = the length of the propagation crack b = the widht of the beam ρ = perimeter of the reinforcing steel l = the length of the reinforcing steel

Case 3
The Case 3 is a continuation of the crack propagation that passed the reinforcing steel bar. The crack at this stage propagated from node G to node I ( Figure  5). The direction of crack propagation was distinguished into two directions, vertical and horizontal to the axes of the reinforcing steel. Crack propagation in the horizontal direction is a representation of a diagonal crack in the structures. The criterion of the crack propagation is based on the magnitude of the crack angle (θcrk) in the global coordinate system. When in the first quadrant, a crack with a crack angle less than 45 degrees is considered propagating horizontally. But if the angle of the crack is between 45 and 90 degrees, then the crack is assumed to propagate vertically. The duplication of node H was to accommodate the crack propagation of node G to node I.  The experimental set-up is illustrated in Figure 7. The loading test was conducted using a loading frame for structural experimentation, and the load was measured by a load cell with a capacity of 500 kN. The load was transferred to the beam using a two-point loading system, 300 mm apart, transferred via a sufficiently rigid WF steel beam. As the load was applied, the LVDTs installed at the center measured the deflection of the beam specimen. Additionally, the strain at the extreme upper and lower part of the concrete in the maximum flexural region, and the strain in the tension steel were measured using strain gauges. While the experiment was in progress, the crack patterns at each loading stage was directly traced and marked on the beam.

Validation
The developed FEM was run for the identical data to the experimental specimen B1. To evaluate the correctness of the developed FEM, the response of the load-versus-displacement curve generated by the FEM program was compared to the experimental data, and presented in Figure 8. In general, the curves showed that the FEM program produced a close approximation to the actual load-displacement behavior in flexure. The response was linear with a constant slope up till the formation of first cracking in the beam's tension zone. Prior to cracking, the slope of the load-displacement curve decreased gradually due to crack propagation within the concrete, the bond-slip between concrete and reinforcing steel, and yielding of reinforcing steel. The program had tended to slightly underestimate the initial crack formation and the crack intersecting with the reinforcing steel when compared to the data obtained from the laboratory specimen. The reason for this deviation could originate from the fact that the program didn't accommodate for the self-weight of the beam.

Figure 8. Validation of Load-displacement Relationship
Visual comparison of the beam's cracking pattern was also accessed to achieve a better understanding to the FEM's flexure behavior of the beam. Figure 9 shows the cracking of the experimental beam, and the crack propagation as predicted by the FEM program. The crack propagation as predicted by the FEM program was in good agreement with the results as observed form the experimental beam. In the area between the points of loading where the shear stresses are relatively very low, it was shown that the cracking started at the most extreme fiber in tension and procreated along a line perpendicular to the horizontal axes. At the nodal points in the line of the applied load, the FEM predicted a vertical crack propagation pattern up till the beam's mid-height. Approaching the center of gravity of the beam's cross section, the shear stresses in the section became more pronounced, resulting in a horizontal cracking in this area suggesting shear failure. Figure 10 shows the stress response of the tensile reinforcing steel at mid-span, as a function of load.

Figure 10. Load versus Reinforcement Tensile Stress
For the same load levels, the experimental beam B1 exhibited wider crack-widths, in combination with less number of cracks, when combined to the FEM beam. The FEM resulted in smaller cracks, that are substantially larger in number. This disparity is seen in Figure 10. The experimental beam exhibited an almost horizontal path after the first crack, while the FEM beam's reinforcement response followed a theoretical path suggesting that cracks were formed uniformly over time. The crack formation for the experimental specimen responded in a different manner, the cracks were formed and widened over a relatively very short period, influencing the tensile stress transfer from the concrete to the steel. This outcome is due to the fact that the stress transfer process is very sensitive to load responses. A finer meshing in the concrete tensile area and a smaller load increment are required to enable accurate modeling of this behavior. Figure 11 shows the moment-curvature response of the beams. The results of the FEM are in good agreement with the experimental curve, despite the differentiation in crack formation and propagation pattern. The curves suggested that within the boundaries of the curvature analysis, the average tensile and compression strain difference between the FEM and experimental beam are relatively small.

Conclusion and Future Research Work
The reinforced concrete structure with the influence of bond-slip and crack propagation was modelled numerically. The modeling of the bond-slip as link element and the representation of crack propagation as discrete crack model was proven to be correct in simulating the nonlinear response of a reinforced concrete structure in flexure. Discrete crack model including the technique of node separation was therefore proven to be thorough in representing the strain discontinuity of the structure due to the crack formation.
Further research work needs to be conducted to expand the finite element model for accommodating the vertically placed reinforcing steel, to model the transverse reinforcement. The crack propagation was, to large extends, correctly predicted by the FEM. However, a finer meshing especially in the region with high stress concentrations will result in a more precise cracking pattern within the beam, a diagonal crack pattern could thus resulted. The FEM program could further be utilized with an import option to accommodate the output of a mesh generator that could create smaller elements in the area of first cracking.