NONLINEAR FINITE ELEMENT ANALYSIS OF DAMAGED AND STRENGTHENED REINFORCED CONCRETE BEAMS

. Bending test of seven reinforced concrete beams are modeled in finite element program to validate the modeling strategies by comparing the structural response of the beams. Three beams in the set are pre-damaged and strengthened with fiber reinforced composites before the bending tests. Cracks are implemented into the model by inserting geometrical discontinuities to represent the pre-damaged beams. Parametric variables such as crack width, length and interval are chosen to simulate different pre-damage levels. Once the proposed modeling strategies are validated by real experimental tests then 196 finite element models are created to study the effects of pre-damage levels on the moment capacity of reinforced concrete beams repaired with CFRP. Results indicate that inclusion of pre-damage levels by means of cracks into the cross sections have significant effect on beams moment capacity. Reference to this paper should be made as follows: Aktas, M.; Sumer, Y. 2014. Nonlinear finite damaged and strengthened concrete beams, Journal of Civil Engineering and Management 20(2): 201–210.


Introduction
The use of externally bonded fiber reinforced composites (FRP) for rehabilitation and strengthening of reinforced concrete (RC) beams has been studied widely in civil engineering discipline. The material properties of FRP composites such as their high strength to weight ratio, good corrosion resistance and ease in application make FRP composites attractive for various industries. Many studies, either numerically or experimentally, have been done to investigate the behavior of externally bonded reinforced concrete elements with FRP in the last 30 years.
Most studies of FRP-rehabilitation have been conducted on virgin, uncracked specimens to estimate the ultimate strengthening capacity. A few experimental studies exist on pre-damaged specimens. However, neglecting the effects of cracks due to the existing pre-damage will lead the designers to incorrect retrofit or strengthening. Sharif et al. (1994) conducted a test program with specimens preloaded up to 85% of their ultimate capacity prior to unloading. They reloaded the specimens after bonding FRP composites. Arduni, Nanni (1997) used a similar procedure with a preloading of 30% of ultimate capacity. Benjeddou et al. (2007) created various damaged beams with different preloading values and with different concrete strength class. They repaired them by varying the amount of carbon fiber reinforced (CFRP) laminates and the dimension of laminate widths. Obaidat et al. (2010) used the finite element method to model the behavior of preloaded CFRP-strengthened beams without simulating cracks. Numerical studies lacks for predamaged beams with existing cracks. Experimental data are crucial to understand the behavior of strengthened RC beams and they are needed to verify the numerical modeling strategies. Nonlinear finite element analysis can provide a powerful tool to study the behavior of concrete structures (Pham, Al-Mahaidi 2005). However, in nonlinear finite element modeling of this kind of problem there are challenges such as: modeling of nonlinear response of concrete both in compression and tension, introducing existing pre-damages like cracks, defining interfacial bond behavior between materials and defining material models for CFRP and adhesive. In literature studies for modeling such behavior are available but there is very limited study for numerical modeling of pre-damaged RC beams, particularly with cracks. Malek, Saadatmanesh (1998) proposed an analytical model to simulate the composite behavior of concrete and steel reinforcement. While they modeled FRP composites, they neglected bond slip behavior between FRP and concrete. Niu et al. (2006) modeled interfacial bond behavior by using zero-thickness interface elements. Baky et al. (2007) also studied interfacial response of FRP strengthened RC beams. Monteleone (2008) modeled cracks by introducing cohesive elements in possible crack regions. Khalfallah (2008) proposed a numerical approach to determine the stress distribution between steel concrete interface near flexural cracks and emphasized the inclusion of tension stiffening in numerical modeling. Qin, Zhao (2009) developed a three dimensional finite element model through Abaqus software to simulate externally prestressed RC simple beams. Numerical simulations were verified by comparing the experimental results. Yang et al. (2009) investigated the flexural performance of FRP strengthened reinforced concrete beams by employing different FRP bonding lengths and by applying different prestressing methods into their experimental and numerical studies. All these studies mentioned above have been concentrated into the modeling techniques and prediction of the load-deflection behavior and ultimate load carrying capacities of FRP-strengthened undamaged beams. Modeling techniques should be developed to capture the effects of existing cracks on the load-deflection behavior of FRP-strengthened damaged beams. Thus, more realistic results can be captured through nonlinear finite element modeling.

Test cases for verification study
The experimental literature has been searched in order that an appropriate research program could be identified to use as a test case for the current verification study. Three different studies found in the literature centered on both FRPstrengthened undamaged and pre-damaged beams. Geometrical layouts of these three different test groups appear in Figure 1. All these experiments were tested under fourpoint loading to achieve a pure bending region. Material properties for concrete, FRP and epoxy used in these experimental program are presented in Table 1.  There is a total of seven experiments in these three groups. Three of these seven experiments are for unstrengthened RC beam and used for verifying the behavior of RC beam only. Only one experiment is for FRPstrengthened undamaged RC beam while the rest three experiments are for FRP-strengthened pre-damaged RC beams (Table 2).

Finite element modeling approach
The commercial multipurpose finite element software package ABAQUS is employed in this research. All the analyses were carried out under displacement controlled loading scheme and all the models are constructed in two dimensions. Modeling existing cracks due to damages of pre-loaded beams was at the heart of the research work reported on in the current paper. The following subsections explains the modeling techniques to simulate load-deflection behavior of FRP-strengthened damaged beams. The fundamental step was to model the nonlinear constitutive relations for concrete, rebar, FRP and epoxy separately.

Uniaxial tension and compression behavior of concrete
Since the compression and tension stress-strain relation of the test specimens were not reported in the test reports these relations were created by using mathematical models from literature. Stress-strain curve of concrete under uniaxial compression was obtained by employing Hognestad parabola (Hognestad 1951) along with linear descending branch. Some modifications were made to this parabola according to CEB-FIP MC90 due to the effects of closed stirrups. Thus, stirrups were not modeled individually but their effects were included in the properties of concrete. Figure 2a represents the uniaxial compressive response of concrete that is consistent with Eqn (1) in which σ is the compressive stress, f cu is the ultimate compressive stress, ε c * is the peak compressive strain, E is the elastic modulus and f c * is the modified compressive strength. More details about these modified values can be found in Arduini, Nanni (1997).
(1) a) Hognestad concrete compressive behavior b) Bilinear tensile behavior Bilinear model was adopted as displayed in Figure 2b for tensile behavior of concrete (Coronado, Lopez 2006). Crack opening was used instead of tensile strain and calculated as a ratio of the total external energy supply (G F ) per unit area required to create crack in concrete. In Eqn (2) tensile fracture energy of concrete, (G F ), is defined as a function of concrete compressive strength, f c *, and a coefficient, G fo , which is related to the maximum aggregate size (CEB-FIP MC90 1993). Several values for G fo are given in Table 3. (2)

Nonlinearity of concrete
In general, the nonlinearity of concrete under compression can be modeled by approaches based on the concept of either damage or plasticity, or both (Yu et al. 2010). Plasticity is generally defined as the unrecoverable deformation after all loads removed. Damage is generally characterized by the reduction of elastic stiffness. For Concrete Damage Plasticity model (CDP) two main failure mechanisms; tensile cracking and compressive crush-ing of the concrete are assumed. The CDP model used herein uses concepts of isotropic damage in combination with isotropic tensile and compressive plasticity to represent the inelastic behavior of concrete (Lubliner et al. 1989;Lee, Fenves 1998). The evolution of the yield surface is controlled by tensile and compressive equivalent plastic strains. As discussed in Lubliner et al. (1989) the unloaded response of concrete specimen seems to be weakened because the elastic stiffness of the material appears to be damaged or degraded. The degradation of the elastic stiffness on the stress-strain curve is characterized by two damage variables for tension and compression, d t and d c , which can take any values from zero to one. Zero represents the undamaged material while one represents total loss of strength (Lee, Fenves 1998). This response of concrete can be depicted as in Figure 3.
E 0 is the initial (undamaged) elastic stiffness of the material and pl ε are compressive plastic strain, tensile plastic strain, compressive inelastic strain and tensile cracking strain respectively. The evolution of the yield surface was defined by these hardening variables (equivalent tensile and compressive plastic strains). The equivalent tensile and compressive plastic strain can be automatically calculated by ABAQUS after the definitions of elastic material behavior. The stressstrain relations under uniaxial tension and compression were taken into account by considering damage variables in Eqns (3) and (4) (Lubliner et al. 1989).
Interface behavior between rebar and concrete was modeled by implementing tension stiffening in the concrete modeling to simulate load transfer across the cracks through the rebar. Tension stiffening also allows to model strain-softening behavior for cracked concrete (ASCE 1982). Thus, it is necessary to define tension stiffening in CDP model. Tension stiffening can be modeled by defining post failure stress-strain relation or by applying a fracture energy cracking criterion (Hu et al. 2004). When there is no reinforcement in significant regions of the model and cracking failure in not distributed evenly, mesh sensitivity problem exists. To overcome mesh sensitivity problem Hillerborg's (1976) fracture energy approach can be used instead of post failure stress-strain relation. In this approach; the amount of energy (G F ) required to open a unit area of crack was assumed as a material property. This postfailure stress displacement relation, taken from ABAQUS, is depicted in Figure 4. Thus, brittle behavior of concrete was defined by stresscracking displacement response rather than a stress-strain response (Fig. 4a).
As an alternative, G F can be implemented directly as a material property. However in this case, a linear loss of strength after cracking was assumed (Fig. 4b). From CDP perspective, both plastic displacement values can be calculated by using the Eqns (5) and (6). In these equations, the specimen length, l 0 , is assumed to be one unit length (l 0 = 1) and ck t u is the cracking displacement.
Finally, by employing these equations "effective" tensile and compressive cohesion stress , t c σ σ can also be defined as: The effective cohesion stresses determines the size of the yield (or failure) surface (Fig. 5). As mentioned before the yield function proposed by Lubliner et al. (1989) and modified by Lee, Fenves (1998) was employed in this study. Yield surface can be defined by introducing four parameters. The Poisson's ratio controls the volume changes of concrete for stresses below the critical value which is the onset of inelastic behavior. When the concrete reaches its critical stress value then it exhibits an increase in its plastic volume under pressure (Chen 1982). This behavior is taken into account by defining a parameter called the angle of dilation. In CDP model dilation angle (ψ) is measured in the p-q plane (p: hydrostatic pressure stress, which is a function of the first stress invariant, q: second deviatoric stress invariant) at high confining pressure and in this study it is accepted 30° as recommended in the literature. ∈, is an eccentricity that defines the rate at which the flow potential function approaches the asymptote. For instance the flow potential tends to a straight line as the eccentricity tends to zero. The ratio of initial biaxial compressive yield stress to initial uniaxial compressive yield stress is defined by σ bo /σ co . Finally, K c is the he ratio of the second stress invariant on the tensile meridian to compressive meridian at initial yield The parameter K c should be defined based on the full triaxial tests of concrete, moreover biaxial laboratory test is necessary to define the value of σ bo /σ co . This paper does not discuss the identification procedure for parameters; ∈, σ bo /σ co , and K c because tests to be verified in this study did not have such information. Since the experimental data were missing for these parameters, 0.1, 1.16 and 2/3 can be adopted for ∈, σ bo /σ co and K c respectively as mentioned in Yu et al. (2010). As mentioned above, damage parameters are required to complete the CDP model. But these damage parameters are not reported for the experiments to be verified in this study. Most of the RC flexural test reports in the literature failed to provide this information. In order to define these damage parameters, some laboratory tests in material level should be done in advance. In this study damage parameter for concrete compressive behavior was obtained from a verification problem given in ABAQUS verification manual. By applying curve fitting method to this example a 3 rd degree polynomial curve was obtained for compressive damage parameter (Fig. 6a). Then the very same equation was applied to get the damage parameters for the test cases to be verified in this study. Moreover tension damage parameters were obtained by checking the gradual loss in intensity of tensile strength of concrete. Tension damage models created for each test group are plotted in Figure 6b.

Steel and CFRP material model
The stress-strain curve of the reinforcing bar is assumed to be elasto-perfectly-plastic as shown in Figure 7a. The parameters needed to specify this behavior are the modulus of elasticity (E s ), Poisson ratio (υ) and yield stress (f y ). However, in the literature, most studies of RC structures strengthened by FRP have assumed the behavior of FRP to be linear.
The behavior of FRP composites was assumed to be linear elastic until reaching failure strain (ε u ) (Fig. 7b).

Adhesive: FRP-concrete interface model
Interface elements were used to model the bond mechanism between CFRP and concrete to simulate the adhesive. The constitutive behavior of interface elements was defined in terms of traction-separation behavior. This model includes initial loading, initiation of damage, and propagation of damage leading to eventual failure at the bonded interface. The material model of the FRP-concrete interface was simulated as a bond-slip relationship between local shear stress, τ, and the relative displacement, δ, at the interface. Lu et al. (2005) studied three bond-slip models using 253 pull tests from literature. They suggest that bilinear model give better results. Figure 8

displays a schematic representation of such behavior with what is used in the present work.
On the other hand three different bond-slip models; the precise, the simplified, and the bilinear models; have been recommended in their study. They suggest that bilinear model gives better results than the others do. Thus, in the current study, the bilinear model, as shown in Figure 8 was adopted.  Table 4. In these equations α 1 is a coefficient and it is recommended as 1.5 (Lu et al. 2005).

Finite element model
After defining material models geometry of the beam was modeled as plotted in Figure 9. For simplicity, stirrups were not considered as a geometrical entity but its effect was considered by introducing a confined concrete model for the RC. The element types used for constructing the finite-element models are listed in Table 5. Interaction between concrete, epoxy and FRP were achieved by the surface tie definition. Steel bars were embedded in concrete with the same degrees of freedom which also means that there is a perfect bond between concrete and steel. The concrete beam was modeled with four-node plain strain element. The reduced integration scheme was adopted to overcome shear locking. All the beams were loaded by displacement control in the vertical direction.

Crack modeling
Cracks were modeled as a geometric entity. Thus, the amount of crack width, length and interval as depicted in Figure 10 were calibrated to catch the experimental results.

Fig. 10. Variables in crack definition
The following matrix was created to explore the crack width, crack length and crack interval (Table 6). Values for these entities were selected after searching experimental literature. However, different type of mesh must be tried to define the crack tips where stress concentration is very high and Type 3 mesh (Fig. 11) was selected among the others. Another main concern in modeling the cracks is to neglect the bond slip between the primary rebar and surrounding concrete. Closely spaced flexure cracks impose a higher force demand on the reinforcement so the bond strength is expected to diminish with loading. When the stress plots are examined it is clear that the rebar at the crack intervals are under more stress than the rest. This can be seen from Figure 12. However, the intensity of the stress is less than the bond stress. Further, the presence of crack widths in the model would be the primary source of delamination of FRP upon loading. As Figure 13 clearly depicted that shear stresses are high at these regions. Delamination will start here as expected. Cover thickness is very important factor for preventing concrete cover splitting especially at these interfaces where stresses are high. Thus, proper cover thickness was provided in order not to allow cover splitting or delamination in this study. Fig. 13. Shear stresses at the interface of CFRP and concrete in the vicinity of the cracks After running 27 models for each test group and comparing the results of numerical models with those from experiments 1 mm, 80 mm and 85%xh were selected for crack width, interval and length, respectively.
These quantities are capable of catching the structural behavior of beams pre-loaded up to 90% of their ultimate capacity. Therefore, no direct crack width, length and interval are defined for any other given damage level.

Verification results
For the credence of the study, it is important to establish the accuracy of the modeling strategies used to conduct this work.
Employing the above modeling strategies the seven constant moment flexural cases were reproduced in the computer in order that suitable structural response results can be generated and compared with the results from experimental program mentioned in Figure 1 and Table 2.
Load deflection curve for un-strengthened RC beams are plotted in Figure 14 (a, b, c). Results for undamaged but FRP strengthened RC beam is given in Figure 14d.
Fig. 14. Verification results for un-strengthened and strengthened beams both without damage All these plots show that finite element modeling techniques applied herein, are valid for uncracked RC beams. Load deflection relation for initially cracked and FRP strengthened damaged RC beams appears in Figure 15 (a, b, c). Based on these results, it appears that the present modeling techniques are sufficiently robust to undertake the outlined parametric study to investigate the effects of existing cracks in RC beams.

Parametric study
In order to present the effects of existing cracks on strength of RC beams, RC beams were created by considering different parameters such as beam length (L/d), aspect ratio of cross section (d/b), stirrup spacing (s), tension reinforcement ratio (ρ tension ), concrete class and aggregate size (d max ). Numerical test matrix for this parametric study is as given in Table 7.  All the numerical specimens were created with the same cross-section width of 120 mm, compression reinforcement of 2 bars with 8 mm diameter and the ratio of shear span to effective depth of 4. Also all the specimens were strengthened with the same amount of CFRP and epoxy as used in Group 2 verification set. 96 numerical specimens were constructed by using the above mentioned parameters.
In order to investigate the effects of cracks on moment capacity, all numerical specimens were also created with cracks (Fig. 16). Dimensions of cracks were observed from studies mentioned in sub-section named crack modeling. Thus, total of 192 analyses were carried out and moment rotation relation for each beam was obtained.
Ratios of maximum moment capacity for both cracked (Mc max ) and uncracked (M max ) beams were calculated to understand the effects of modeling the existing cracks on strengthened RC beams (Table 8

Parametric results and discussion
The Strength Reduction Factor (SRF) is defined as the ratio of maximum moment capacity of uncracked beams to the cracked beams. Maximum value of strength reduction factor is found as 35.78%, which means the existence of cracks in pre-damaged RC beam decreased the moment value from 27.28 kNm to 17.52 kNm. Strengthening calculations for such beams should be done with correct initial capacity values in order to provide adequate retrofit to RC beam. Aggregate dimension, which is an important parameter for defining tensile fracture energy of concrete (G F ), does not have a significant effect on moment capacity for cracked concrete beams. Stirrup spacing, reinforcement ratio and concrete class values do not affect SRF values significantly. The most important parameter that has an effect on SRF value is found to be cross section aspect ratio (d/b). While the value of SRF vary from 0.2% to 23.86% with an average value of 6.86%, for square beams (d/b = 1), it does, however, vary from 3.48% to 35.78% with an average value of 18.4%, for rectangular beams (d/b = 2).

Conclusions
This study demonstrates that the load carrying behavior of a FRP reinforced damaged concrete beam should include the preexisting cracks in the model. Among the important findings are: 1. Experimentally observed loading trends and magnitudes for entire loading range of both RC beams and RC beams strengthened with FRP can be captured by employing CDP modeling approach. However, each numerical model has to be verified with real experiment for further parametric studies. 2. Mesh density, dilation angle and concrete fracture energy are needed to be investigated carefully for proper modeling. 3. Existing cracks can be modeled by introducing crack width, crack length and crack interval into the model. Results from parametric study to investigate these values indicated that 1 mm crack width, 85% crack length and 80 mm crack interval can capture the pre-damaged properties of RC beams. 4. Existing cracks in pre-damaged RC beams have to be considered in numerical models since moment capacity for pre-damaged beam is different from those undamaged beams. In this study, up to 35.78% difference between moment capacities is found for some cases. 5. RC beams with rectangular cross section are more sensitive than RC beams with square cross section in terms of having initial cracks modeled in finite element model. Along with pre-damages in reinforced concrete beams poor material quality, poor design and inadequate stirrup spacing are other fundamental problems for determining the moment capacity of RC beams. However, this study proves that the above mentioned problems are not as important as the geometrical discontinuities like cracks.
When the designers calculate the moment capacity of the beam they need to consider the existing cracks or damages. Failing to do this will estimate moment capacity more than the actual real case. Thus, the amount of CFRP will be less than the required one. Strengthening the damaged beams with this unsufficient CFRP amount will be unsafe. Geometrical damages should be implemented into the finite element analysis of damaged beams to avoid such case.