NUMERICAL INVESTIGATION OF REINFORCED CONCRETE BEAM DUE TO SHEAR FAILURE

This paper investigates the possibility of using a multi-surface plasticity model to predict shear failure in reinforced concrete beams. The analysis is carried out using the in-house software called 3D-NLFEA. The constitutive model for the concrete material is based on the plasticity-fracture model, which had previously developed to simulate the behavior of concrete cover spalling in reinforced concrete columns. To obtain the asymmetric shear failure pattern, random material properties imperfection for each meshed element is used. Two beams available in the literature are investigated and compared with the analysis results using 3D-NLFEA. From the comparisons, excellent agreement between the analysis and the test result was obtained. 3D-NLFEA can predict the peak load accurately. The peak load prediction only varies 2.19% for beam OA1 and 3.28 % for beam OA2, and it was lower than the test results. The failure crack patterns also show a typical diagonal crack extension from the support to the loading steel plate.


INTRODUCTION
Shear failure in the beam without shear reinforcement is still one challenging task in numerical simulation. A proper concrete constitutive model holds the key to capturing the load-deflection curve of the tested specimen and the cracking pattern. Typically, reinforced concrete (RC) beams without reinforcement that fails in shear have their strength and stiffness firstly degraded when the localized flexural crack occurs. The energy released from the flexural crack is then transferred to the longitudinal reinforcement, creating localized stress distribution on the longitudinal bars. As the load increases, the diagonal shear crack will occur, followed by a sudden drop in the load-carrying capacity of the beam. This diagonal shear crack can be avoided by providing shear reinforcement, and thus the shear strength and ductility of the beam can be improved significantly. The ductile behavior of the RC beams can also be obtained if the flexural reinforcement was designed to be yielded first.
For design purposes, the use of empirical equations was found to be adequate. It is a fast and robust method to estimate the load-carrying capacity of the beams. However, it should be noted that this empirical equation is derived from hundreds or more tested beam specimens with and without reinforcement failing in shear. Hence, the accuracy of the empirical formulation may vary significantly. On the other hand, using the non-linear finite element method to predict the behavior of RC beam failing in shear was more preferred as it can give a more accurate prediction of the peak load and show the cracking pattern in real-time. However, a concrete constitutive model that can capture the localized flexural and shear cracks is necessary. Although some researchers use the Modified Compression Field Theory (MCFT) [1][2][3] , which is specialized for predicting shear behavior of RC beams, in this paper, a constitutive model that is based on the multi-surface plasticity-fracture model [4] is used. The MCFT model in Vecchio and Collins [1] accurately predicted the load-deflection curve and the cracking patterns for both the flexural and shear cracks. However, the reinforcing bars are modeled as smeared reinforcement rather than embedded formulation, which requires further justification to determine the reinforcement ratio. For a complex reinforcement arrangement, it would not be practical to precompute the reinforcement ratio for the input in the modeling.
In applying non-linear finite element analysis of RC element, typically for the one that utilizes plasticity-fracture model, the smeared crack model is one of the most used methods to capture the crack propagation in concrete. This crack band theory fits well to be used together with the embedded truss formulation of the steel bar or even with the smeared reinforcement model. Inside the smeared crack model, there is a two cracks theory often discussed. One of them is called the fixed-crack model, and the other one is a rotating-crack model. Both fixed-and rotating-crack models are an extension of the classical isotropic fracture model [5] . Although the application of the smeared crack model with fixed-and rotating-crack model have been discussed widely elsewhere, the use of isotropic fracture model can also be used.
In this paper, an isotropic fracture model is used to investigate the RC beam failing in shear. The prediction of peak load, first cracking load, and crack propagation at some points in the load-deflection curve are presented. Typically, due to the symmetry of the modeled specimen, some researchers modeled only a half portion of the beam. With a half portion of the beam [6] , the diagonal cracks occurred in the center-line left and right, which is not true based on the experimental investigation. The localized shear crack was asymmetric from the experimental investigation and only occurred at one-half side of the beam. Another possible reason for using half-symmetry of the model is to cut the runtime of the program. It is also possible that if a full-beam model were used, the symmetric localized shear crack might occur, although it would be unlikely due to numerical precision error during the computation.
Hence, to induce an asymmetric crack pattern, a full model of the beam specimen with imperfections is required. These imperfections can be in the form of material or geometric imperfections. For RC structures, geometric imperfection may not significantly affect the thickness of the member is quite thick. On the other hand, material imperfection can be more successful for RC structures as it will create weak elements to initiate localized cracks in the model prematurely. Hence, it is possible to get an asymmetric cracking pattern. For that purpose, in this paper, a simple random material imperfection that generates a random number with some margin around the mean value is used. The random material properties are applied to generate a random value for the concrete compressive strength. Inside the 3D-NLFEA, the input material properties of concrete compressive strength will also affect the concrete tensile strength. As the concrete tensile strength varies for each of the meshed elements, the starting point of the localized cracks could start elsewhere. For each different random input, different starting crack locations are expected.

MATERIAL AND METHOD
The non-linear finite element package used in this research is called 3D-NLFEA. 3D-NLFEA is an in-house finite element package developed by the Piscesa [4,7] . 3D-NLFEA uses an Initial Stiffness Method (ISM) to maintain the equilibrium between the internal and external forces during the global iteration. To accelerate the convergence, the so-called Tangent Stiffness Like Projection Method (TSLPM) is used along with the Process Modification (PM) [7] or with the accelerated ISM [8] . To capture sudden drops in the load-deflection curve when the peak load reaches, the arc-length method is used [9] . 3D-NLFEA uses the SALOME platform as the preprocessor and PARAVIEW as the post-processor [10] . The concrete constitutive model used in this paper is the Piscesa et al. [41112] plasticity-fracture model. This plasticity-fracture model consisted of a multi-surface plasticity model, the Menetrey-Willam failure surface for concrete under compression, modified previously [12] to predict better the peak and residual loads for confined concrete and Rankine tension cut-off failure surface. Non-associative flow rule with a non-constant plastic dilation rate is used for concrete under compression and an associative flow rule for the tension cut-off failure surface. The returned mapping used is based on the implicit integration method [13,14] . When both failure surface actives, a nested return mapping outlined in Piscesa et al. [4] is used. The concrete constitutive model in 3D-NLFEA is similar to that of the plasticity-fracture model in the ATENA software package [5,15] . One can also find some examples of using the 3D-NLFEA finite element package in modeling RC structures in Christianto et al. [16] , Prasetyo et al. [17] , Alius et al. [18] .
The constitutive model for the steel reinforcing bar is modeled using an elastic-perfectly plastic model. The perfect bond model assumption is used between the bars and the solid parent element. Hence, there are no slips occur between the bars and the concrete. The steel reinforcing bar is modeled using an embedded element formulation [19] . The concrete element is modeled using the hexahedral element with BBar element technology [20] which uses different integration for the hydrostatic and deviatoric parts of the element B-Matrix. A reduced integration is used for the hydrostatic part, while for the deviatoric part, full integration is used.
There are two beam specimens available in the literature that are investigated. Beam OA1 and OA2 were tested by Bresler and Scordelis [21] . These beams are often used as the benchmark in numerical simulation to study the shear failure in RC beams. In 2004, the almost similar beam configurations of OA1 and OA2 were retested by Vecchio and Shim [6] . Figure 1 shows the cross-section for Beam OA1 and OA2. Beam OA1 has a width of 310 mm and a height of 556 mm.
On the other hand, beam OA2 has a width of 305 mm and a height of 561 mm. Beam OA1, prepared by Vecchio and Shim [6] , was constructed using two M25 and two M30 rebars, but it has almost the same reinforcement ratio as the Beam OA1 in Bresler and Scordelis [21] . Beam OA2 was constructed using two M25 and three M30 rebars. The arrangement of the longitudinal bars for both beams was shown in Figure 1 . The total lengths for beam OA1 and OA2 are 4100 mm and 5010 mm, respectively. The beam span (measured from the center-to-center supports) for beam OA1 is 3660 mm, and for beam, OA2 is 4570 mm.
The concrete compressive strength (fc) for Beam OA1 is 22.6 MPa, while for beam OA2 is 25.9 MPa. To consider the different strengths between the concrete cylinder and the concrete beam due to different size and curing conditions, a factor of 0.85 is used to reduce the actual concrete strength of the beam. The elastic modulus for the concrete material was computed using the equation in ACI 318-19. This yields the concrete modulus elasticity for beam OA1 and OA2 are 18137 and 18046 MPa, respectively. The concrete tensile strength (ft) used as an input is taken about 80% from the estimated concrete tensile strength without silica fume as outlined in Attard and Setunge [22] , Samani and Attard [23] . Hence, for beam OA1 and OA2, the concrete tensile strengths are 1.669 MPa and 1.829 MPa, respectively. In Vecchio and Shim [6] , the concrete tensile strengths are estimated from the concrete compressive strength, which gives 1.568 and 1.679 MPa for beam OA1 and OA2. The yield strength for the M25 bar is 440 MPa, while for the M30 bar is 436 MPa.
In the plasticity-fracture model, which utilizes the crack band approach, both the fracture energy and the internal length scale of the meshed element play an important role. The fracture energy used for both RC beams is computed using the CEB-FIP  2010 [24] , with the maximum aggregate size was set to 16 mm. The input for the base fracture energy (GF0) is 0.03 N/mm. The fracture energy input (GF) in the plasticity fracture model can be computed using Eq. 1.
where is the uniaxial concrete compressive strength in MPa. During the preprocessor stages, the maximum mesh size of the hexahedral element is set to 25 mm. Hence, theoretically, we should set the internal length scale to 25 mm for both beams. However, after several trials and run to get a closer prediction of the peak load with the test result, we found the optimum value for the internal length scale to be 22 mm. This can be attributed to higher fracture energy stored in the actual beam compared to the code-based value. It should also be noted that code prediction for the fracture energy should be conservative for design purposes, proving that it was conservative. Figure 2 shows the modeled beam specimens in 3D-NLFEA. The beam OA1 was modeled using an eight-node hexahedral solid element with a total of 50,258 elements (includes the steel loading plate and steel supports). The rebar element was modeled using a two-node truss element with a total of 656 elements. Beam OA2 has a larger size than beam OA. Hence, with the same maximum size of the concrete elements. The solid element used to model beam OA2 was 60,710 elements (including steel loading plate and steel supports). The rebar element was modeled using a two-node truss element with a total of 1005 elements. The loading is given with a displacement control of -0.1 mm for every load step. A monitoring displacement point at the bottom of the beam in the midspan was added for comparison purposes with the available test results.

RESULTS AND DISCUSSION
The following section explains the results and aalysis of our observations. The results are presented in four different aspects, i.e. load-deflection behavior, parameter identification, crack propagation, and 3d-NLFEA simulation.    Figure 3 shows the load-deflection curve for beam OA1 and OA2. In Figure 3 , both test results from Bresler and Scordelis [21] and Vecchio and Shim [6] . As for the numerical simulation, in addition to the 3D-NLFEA results, numerical simulation results from VecTor 2, which utilize the MCFT approach, were also shown in Figure 3 . An important observation from Vecchio and Shim [6] is that the stiffness of the retested RC beams with almost similar configurations is not the same. The beam tested by [21] shows higher stiffness compared to the beam tested by Vecchio and Shim [6] once the flexural cracks are formed.

Load-Deflection Behavior
As shown in Figure 3 , The initial stiffness prediction using 3D-NLFEA for both beam OA1 and OA2 agrees well with the initial stiffness of the test result. This indicates that the inputted elastic modulus for both beams is correct. By back calculating the concrete density as the required input when computing the elastic modulus, it returned the values for the concrete density equal to 2100 and 2000 kg/m3 for beams OA1 and OA2, respectively. The normal concrete density without the reinforcement, on average, is about 2200 kg/m3. In Vecchio and Shim [6] , we found the input for the elastic modulus to be considerably high. It was reported that for beam OA1, the concrete elastic modulus was 36,500 MPa, while for beam OA2, it was 32,900 MPa. By back calculates the concrete density for beam OA1 and OA2 concerning each beam concrete compressive strength, we get the value for the concrete density equal to 3,170 kg/m3 and 2,827 kg/m3 for beam OA1 and OA2, respectively.
Once the onset of localized flexural cracks occurred, the beam stiffness lowered. The onset location of the localized flexural cracks for both beams was close to the one investigated from the test result. This indicates that the inputted tensile strength for concrete material also agrees well with the test result. However, the concrete tensile strength input in 3D-NLFEA was adjusted to only equal 80% than the prediction as shown in Attard and Setunge [22] . In Vecchio and Shim [6] , the concrete tensile strength was even lower than the user input in 3D-NFLEA. Both 3D-NLFEA and VecTor 2 showed almost similar stiffness once the localized flexural cracks occurred.
After the localized flexural cracks are formed, the longitudinal bars in the cracked element have their stresses increase significantly. The shear stress along the perimeter of bars that are not cracked will increase until it forms another flexural crack when the tensile stress in concrete is reached. This explains the flexural crack propagation in concrete beams. The stiffness of the beam, as the flexural crack progressed, degraded at a lower rate. Small load capacity drops occurred in the load-deflection curve. This condition continues further up to the peak load.
For beam OA1, when the axial load reaches about 250 kN, VecTor 2 load-deflection curve softens and deviates from its previous loading path while the 3D-NLFEA curve is steadily increasing until the peak load is reached. The test result also showed similar behavior with 3D-NLFEA. For beam OA2, 3D-NLFEA and VecTor 2 performed about the same up to the axial load around 200 kN. After that point, VecTor 2 response is degraded, and the load-deflection curve deviates from its loading path. On the other hand, the 3D-NLFEA loading path keeps on increasing up to failure. At peak, the initiation of shear band failure, which is typical in the concrete beam without shear reinforcement, occurred. The crack bandwidth is larger than the flexural cracks. The energy stored in the beam was released, and a sudden drop in the peak load carrying capacity occurred. Table 1 shows the comparisons of the peak loads from the test results and the numerical models. Table 2 shows the comparisons of the mid-span displacement at peak load between the test results and the numerical models. As shown in Table 1, for beam OA1, the test result from Bresler and Scordelis [21] and Vecchio and Shim [6] showed that the peak load does not vary significantly, but the mid-span displacement was found to be significantly different. From the comparison for beam OA1, 3D-NLFEA shows an excellent axial load prediction, while VecTor 2 was more conservative. However, the mid-span deflection for 3D-NLFEA was higher than the beam tested by Bresler and Scordelis [21] but lower for the beam tested by Vecchio and Shim [6] . On the other hand, VecTor 2 mid-span displacement was found. to be the largest and was in close agreement with the beam tested by Vecchio and Shim [6] .
For Beam OA2, the test result from Bresler and Scordelis [21] showed a higher axial load than Vecchio and Shim [6] . 3D-NLFEA prediction was found to be in close agreement with the test result from Vecchio and Shim [6] and VecTor 2, again, showed conservative peak load prediction. As for the mid-span deflection for 3D-NLFEA was in good agreement with the test result from Vecchio and Shim [6] but lower if compared with the test result from Bresler and Scordelis [21] . For the mid-span deflection, VecTor 2 prediction was in close agreement with Bresler and Scordelis [21] . Overall, both 3D-NLFEA and VecTor 2 can be used to predict the response of RC beams failing in shear.

Parameter Identification
To gain an insight into the difference between the test results and the numerical models, the authors tried to correlate some parameters that may or may not affect the results. For the longitudinal bar properties, there was barely a difference. Hence, the cause of the significant difference caused by the reinforcement may not be significant. On the other hand, the concrete material properties can vary significantly. The aggregate size and composition may affect the fracture energy, initial micro-cracks, and the strength of the concrete material. Different curing conditions and loading parameters may also affect the behavior of the tested beam. As for the difference in numerical results, 3D-NLFEA utilizes different concrete constitutive models, incorporating a plasticity-fracture model with isotropic fracture formulation and embedded rebar formulation.
On the other hand, VecTor 2 uses the MCFT model, smeared reinforcement, and smeared cracking model. The fracture energy provided by the concrete, which is the function of the aggregate size, shape, and strength that was inputted in the concrete constitutive model, can result in different observed behavior. Hence, this author thought that for the test results, despite different measurement equipment accuracy, the concrete compressive strength, tensile strength, and fracture energy were one of the causes of the significant difference in the load-deflection curve. However, for the numerical results, the parameters were also the possible causes but more importantly, the used concrete constitutive model and modeling technique which these authors believe also contributed to the differences in the load-deflection curve.

Crack Propagation
To better understand the crack's propagation in the beam specimen, the maximum strains at three displacement points in the load-deflection curve are investigated. Figure 4 shows the maximum strains for beam OA1. As shown in Figure 4 , the cracks are mainly flexure at displacement monitored 4.0 and 6.0 mm. When the displacement is, However, at displacement 6.0 mm, flexural cracks start to incline and propagates to the loading point. At the peak point (at displacement 7.9 mm), shear cracks occurred with larger bandwidth than the flexural cracks.
It is also worth noting that the shear failure pattern was asymmetric. Hence, modeling half a beam should not be used, and we can use random material properties imperfection to trigger the asymmetric failure pattern. For beam OA2, the cracking pattern was shown in Figure

CONCLUSION
This paper has presented a non-linear simulation of an RC beam failing in shear. A plasticity-fracture model was used for the concrete constitutive model. An in-house 3D-NLFEA software package that utilizes a full three-dimensional model of the beam and random material imperfection to break the model's symmetry was used in the numerical simulation. The prediction of the peak load for both beams using 3D-NLFEA was better than the prediction using VecTor 2. This showed that the isotropic fracture model could be used to predict an RC beam that fails in shear. However, the mid-span displacement at peak predicted using 3D-NLFEA was lower than one series of the test result. From the cracking pattern investigation, the localized flexural cracks were well observed in 3D-NLFEA. The localized shear bands were also captured well, and it was asymmetric as expected due to the use of random material imperfection. This showed that the 3D-NLFEA software package was excellent in modeling concrete with various confinement devices and can be used to predict the shear behavior of RC elements. Further research should be carried out to investigate more RC beams that fail due to shear and failure in flexure and torsion. RC columns with high shear forces should also be investigated to ensure the confinement reinforcing bar is adequate to sustain the axial load under large lateral deformation.