A non-linear procedure for the numerical analysis of crack development in beams failing in shear

In this work, a consistent formulation for the representation of concrete behavior before and after cracking has been implemented into a non-linear model for the analysis of reinforced concrete structures, named 2D-PARC. Several researches have indeed pointed out that the adoption of an effective modeling for concrete, combined with an accurate failure criterion, is crucial for the correct prediction of the structural behavior, not only in terms of failure load, but also with reference to a realistic representation of crack initiation and development. This last aspect is particularly relevant at serviceability conditions in order to verify the fulfillment of structural requirements provided by Design Codes, which limit the maximum crack width due to appearance and durability issues. In more details, a constitutive model originally proposed by Ottosen and based on non-linear elasticity has been here incorporated into 2D-PARC in order to improve the numerical efficiency of the adopted algorithm, providing at the same time an accurate prediction of the structural response. The effectiveness of this procedure has been verified against significant experimental results available in the technical literature and relative to reinforced concrete beams without stirrups failing in shear, which represent a problem of great theoretical and practical importance in the field of structural engineering. Numerical results have been compared to experimental evidences not only in terms of global structural response (i.e. applied load vs. midspan deflection), but also in terms of crack pattern evolution and maximum crack widths.


INTRODUCTION
he analysis of the response up to failure of reinforced concrete (RC) structures through numerical techniques requires the adoption of effective material constitutive models, able to correctly represent the behavior of concrete and steel at the element level.This represents a quite complex task, since as loading increases, RC behavior is influenced by several non-linear mechanisms which often interact with each other, such as concrete cracking and crushing, aggregate interlock, bond-slip behavior, dowel action and yielding of reinforcement.Therefore, realistic simulations of RC structural behavior require a correct description of each of these phenomena into the adopted material model, which can be subsequently implemented within the Finite Element (FE) framework.In this context, among the several constitutive laws proposed in the past within smeared crack formulations ( [1][2][3][4][5][6], just to mention some), 2D-PARC model [7] represents an effective tool to perform non-linear analyses up to failure of RC structures, since it is able to account for the most influencing aforementioned contributions.This model is structured in a modular framework, so that each mechanical phenomenon is indeed individually analyzed by using a proper constitutive law and then the corresponding contribution is inserted into a material stiffness matrix.In more detail, the present work focuses on the evaluation of concrete contribution before and after the development of cracking.As regards concrete modeling, the basic requirements concern the choice of an accurate failure criterion to be employed in conjunction with a proper non-linear stress-strain relationship.The Seventies have seen the development of several T failure criteria derived from different failure data (e.g., [8][9][10][11]), which have pointed out that concrete under a biaxial (or triaxial) state of stress exhibits different stiffness, strength and ductility than under uniaxial loading.Consequently, the strength characteristics of concrete under a general multi-axial state of stress cannot be reproduced into a material model by directly using experimental uniaxial stress-strain curves.To represent the problem, several mathematical material models have been proposed, which can be substantially divided into orthotropic, non-linear elastic, plastic or endochronic ones [12].With reference to orthotropic formulations, a fairly simple solution can be obtained by referring to the concept of equivalent uniaxial strain, as originally proposed by Darwin and Pecknold in [13] and implemented into 2D-PARC model [7].An interesting model based on non-linear elasticity was instead proposed by Ottosen [14].This formulation provides non-linear stress-strain relations for concrete by only properly changing the secant values of Young modulus and Poisson ratio.In this way, even if the model is able to realistically represent concrete behavior under a general stress state, its calibration is quite simple, only requiring experimental data obtained by standard uniaxial tests.For its flexibility and numerical feasibility, it has immediately appeared suitable for the implementation into FE programs, and applicable to the analysis up to failure of different types of RC structures [15,16].Basing on the approach proposed by Ottosen [14] and on the implementation performed by Barzegar [16], concrete modeling into 2D-PARC has been here revised.This modified formulation has been verified with reference to the analysis of RC beams without shear reinforcement [17,18], where concrete modeling assumes particular relevance both before and after crack pattern development.Furthermore, this represents a structural problem of significant theoretical and practical importance, subjected in the past to great experimental efforts [19,20], since the comprehension of the mechanisms of shear transfer across cracks and failure in beams without shear reinforcement can improve the knowledge of concrete contribution on shear strength also in beams with web reinforcement.

NUMERICAL MODEL
he behavior of reinforced concrete (RC) beams without shear reinforcement has been herein studied through a non-linear finite element (NLFE) procedure.In order to account for reinforced concrete mechanical non-linearity, a suitable constitutive model, named 2D-PARC [7], has been adopted.This constitutive model is based on a smeared-fixed crack approach and its theoretical basis, which has been deduced for a RC membrane element subjected to general in-plane stresses, can be found in details in [7] and in [21,22] with reference to its extension to the case of steelfiber reinforced concrete (SFRC) elements or to the 3D case.In the uncracked stage, concrete and steel are schematized like two material working in parallel, by assuming perfect bond between them.When the principal maximum stress violates the failure envelope in the cracking region, crack pattern is assumed to develop with a constant spacing am1.Afterwards, a strain decomposition procedure is adopted, by subdividing the total strain into two components, respectively related to RC between cracks and to all the resistant mechanisms that develop after crack formation (i.e.aggregate bridging and interlock, tension stiffening and dowel action).These resistant contributions are expressed as a function of two main variables, namely crack width w1 and sliding v1, and included into the crack stiffness matrix [D cr1 ].The behavior of RC between cracks is instead described by adopting the same approach used in the uncracked stage, even if a slight modification is operated on both concrete and steel stiffness matrices, [D c ] and [Ds], so as to account for the degradation induced by cracking.The cracked RC stiffness matrix is then expressed in the following form: where [I] is the identity matrix.This formulation has been successfully applied to the analysis of different types of structures (such as panels, beams, slabs, etc…), providing a good prediction of experimental evidences both in terms of load-deformation response up to the ultimate capacity of the considered element, and in terms of crack pattern evolution and failure mode (e.g.[7,23,24]).However, the proposed algorithm is quite complex and requires long calculation times, which are mainly related to the approach followed in the evaluation of concrete stiffness matrix [Dc].Both in the uncracked and cracked stage, concrete is indeed modeled as an orthotropic, non-linear elastic material subjected to a biaxial state of stress, which is duplicated by means of two equivalent uniaxial curves, following Darwin and Pecknold approach [13].These uniaxial curves for concrete in compression and in tension report the actual stress as a function of an equivalent uniaxial strain, which is in turn determined according to [25].These curves are characterized by maximum strength values derived from an analytical T biaxial strength envelope based on the one suggested by Kupfer et al. [8] (see also [7] for further details), and by material secant moduli in the two orthotropic directions, to be inserted into the matrix [Dc].Moreover, in the cracked stage the terms of matrix [D c ] are adequately softened through an empirical damage coefficient  related to crack width w 1 , whose expression can be still found in [7].The need of adopting a so refined model for the description of concrete behavior, instead of a simple linear-elastic matrix, is mainly aimed to a correct simulation of those elements made of plain concrete or reinforced in a single direction, such as, i.e., beams without shear reinforcement.In these cases, the structural behavior is indeed mainly governed by concrete performances, and consequently its correct constitutive modelling is mandatory for a realistic prediction of both element stiffness and strength.This work illustrates an alternative procedure that still allows a quite sophisticated representation of concrete behavior but requires a lower computational effort.Its implementation into the 2D-PARC model is also briefly discussed.

Modeling of concrete behavior in the uncracked stage
The constitutive relation herein adopted for concrete modeling represents a specialized 2D form, according to Barzegar [16,26], of the 3D non-linear elastic model originally proposed by Ottosen [14,15].The main advantages of this model lay on its simple definition, since the stress-strain relation is expressed as a function of only two parameters, i.e. the secant values of the Young modulus E c and of the Poisson coefficient , which are properly modified to account for material non-linearity.As a consequence, concrete stiffness matrix [Dc] can be written as: in the global x-y co-ordinate system.Hence, with respect to the previous formulation implemented into 2D-PARC relation, this model depends on a lower number of parameters and allows to bypass the evaluation of the two equivalent uniaxial strains, so reducing the required computational effort.Moreover, this model is very flexible, since it can be used in conjunction with any failure criterion, by simply modifying a single parameter, the so-called "nonlinearity index", as discussed in the following.Similarly to the original approach of 2D-PARC, the failure envelope proposed by Kupfer et al. [8,9] has been still considered, even if its analytical expression has been slightly modified according to [16] in the region corresponding to tension-compression, so as to avoid possible discontinuities in those points where the maximum principal tensile stress is close to zero.Fig. 1 shows the adopted failure envelope and the analytical expressions describing each considered region (tension-tension, tension-compression and compression-compression); the bold line indicates the part of the curve effectively implemented into the 2D-PARC model, according to its conventions (that is  1c ≥  2c ).Concrete secant elastic modulus E c is computed using a parameter called "nonlinearity index", , which depends on how far the current stress point is from failure and is then related to the amount of non-linearity in the stress-strain curves [14,16].In case of biaxial compression this parameter can be evaluated through the expression: where 2c is the maximum principal compressive stress and 2fin represents its corresponding value on the failure envelope, determined by keeping fixed the other principal compressive stress  1c (being  1c ≥  2c , and assuming compressive stresses as negative).Thus,  < 1,  = 1 and  > 1 respectively correspond to stress states located inside, on, and outside the considered failure curve.When tensile stresses are present, the nonlinearity index is instead computed in terms of effective stresses.To this aim, the actual state of stress ( 1c ,  2c ), where at least  1c is a tensile stress, is properly turned into an "equivalent compressive case", by superposing an hydrostatic pressure -1c to the existing stress field.In this way, a new state of stress ('1c, '2c) = (0, 2c -1c) is obtained and the nonlinearity index is evaluated as: where f c is the uniaxial compressive strength of concrete.By following this procedure, the value of  is properly reduced when tensile stresses occur; thus,  < 1 always holds in this case. - Tension-tension After having determined the nonlinearity index, concrete secant elastic modulus Ec can be then calculated as: where E ci is the initial value of concrete Young modulus, D is a compressive post-peak nonlinearity parameter that determines the degree of strain softening when concrete crushing occurs (see [14,16] for details), and E' cf is the secant modulus corresponding to peak stress.When a tensile stress is present, E'cf is simply evaluated as in case of uniaxial compression, i.e.E' cf = E cf = f c / c0 , while for biaxial compression the following relation is adopted: A being the ratio between the initial value of concrete Young modulus and the secant one corresponding to peak stress (Eci/Ecf), while the term x takes into account the dependence on the actual loading and is evaluated through the relation: The first addend of Eq. 7 represents the failure value of the invariant 2 c J f .Based on the definition of the nonlinearity index  (Eq.3), the following expression can be found: Similarly, 1 3 is the value of   2 c f J f for uniaxial compressive loading.
As regards the evaluation of the secant value of Poisson coefficient to be inserted into Eq.2, it has been experimentally observed that in presence of compressive stresses concrete tends first to compact, while after the appearance of microcracks it tends to dilate.As a consequence, its value should be properly adjusted during the analysis so as to correctly represent this behavior.Following the procedure proposed in [14,16], the Poisson coefficient is kept fixed until  reaches a "limit value"  a equal to 0.8, and afterwards it is updated by applying the following relation: where  i is the initial value of the Poisson coefficient (assumed equal to 0.2) and  f represents its secant value at peak (approximately equal to 0.36).
The so obtained secant values of concrete elastic modulus and Poisson coefficient have been inserted into the concrete stiffness matrix [D c ] (Eq.2); in this way, only the constitutive relation adopted for concrete modeling has been modified, by keeping unchanged the global structure of 2D-PARC algorithm.It should be observed that all the above described procedure is applied until either cracking or crushing occurs.Cracking takes place when the current stress state reaches or violates the failure envelope in the cracking region, i.e. 1c ≥ 1max in tension-compression or simply  1c ≥ f ct in case of tension-tension, being f ct the tensile strength of concrete.In this case, the fixed-smeared crack approach previously described is applied (see [7] for further details), even if the modeling of concrete between cracks has been properly modified as described in the following Subsection.On the contrary, when the stress state reaches or violates the failure envelope in the crushing region, i.e.  ≥ 1 in case of biaxial compression or 2c ≤ 2max in case of compression-tension, a very simplified procedure is applied herein.For the considered case studies, only small portions of the modeled structures can enter indeed in the post-crushing regime during the analysis (e.g.near concentrated loads or at supports) without affecting the global behavior, which is instead much more influenced by the appearance of tensile cracks.However, concrete crushing should be included in the model formulation so as to avoid numerical difficulties or wrong predictions of the ultimate failure load.To this aim, in 2D-PARC model a reduced Young modulus Ec is simply considered, equal to 25% that of undeformed concrete, as suggested also in [27].

Modeling of concrete behavior in the cracked stage
The above described formulation has been applied also in the cracked stage for the evaluation of the stiffness matrix of concrete between cracks, [D c ].However, in order to include damaging due to the presence of cracks, a proper reduction in the concrete compressive strength and stiffness has been operated.To this aim, the biaxial concrete failure envelope shown in Fig. 1 has been properly reduced, following the relation suggested in [26]: where max is the current maximum principal tensile strain in cracked concrete and * c f is the modified value of the uniaxial compressive strength at peak.The corresponding modified strain *  co can be in turn calculated as a function of the initial value of the strain c0 corresponding to peak stress in uniaxial compression, through the expression:

COMPARISONS WITH EXPERIMENTAL OBSERVATIONS
he capability of the proposed model to describe RC global behavior and crack pattern evolution has been verified against the results of two well-documented experimental programs concerning beams without shear reinforcement [17,18].Among several research projects, the one carried out by Vecchio and Shim [17] has been selected owing T to its effectiveness, being a duplicate of the well-known Bresler and Scordelis beam tests [20], always regarded as benchmark data.The choice of the experimental program undertaken by Podgorniak-Stanik [18] has instead been related to the availability of several experimental data monitored during test execution, mainly concerning the crack pattern evolution with increasing loads.
Description of experimental tests [17,18] The attention has been initially focused on three beams without stirrups, named OA1, OA2, OA3, tested by Vecchio and Shim [17].These specimens had the same rectangular cross section -305 mm wide and 552 mm deep -and a net span respectively equal to 3660 mm, 4570 mm and 6400 mm, corresponding to an increasing amount of tension reinforcement, heavy enough to make the beams critical in shear.The main geometrical details of the considered specimens and their reinforcement arrangement are summarized in Fig. 2a.The three beams were characterized by a progressively increasing concrete compressive strength f c (which was equal to 22.6, 25.9 and 43.5 MPa for specimen OA1, OA2 and OA3, respectively).The main characteristics of the adopted reinforcement, both in terms of geometrical details and steel properties, can be found in [17], to which reference is made.BN50D 2700 3000 2 M20, 1 M25, 10M10 Figure 2: Geometric dimensions (in mm) and reinforcement arrangement of the analyzed beams: (a) OA [17] and (b) BN series [18].
All the tests were performed under loading control, with a central point load, until the approaching of the ultimate stage, when the procedure was switched to displacement control so to allow the evaluation of the post-peak behavior.As already mentioned, the purpose of this experimental program was to recreate, as much as possible, the Bresler and Scordelis tests [20], in terms of geometrical dimensions, reinforcement details, material strengths and loading.Compared to these latter, the beams tested by Vecchio and Shim [17] exhibited indeed a very similar behavior, with only few minor differences; as a consequence, only the specimens described in [17] have been considered in the FE analyses reported herein.
In addition to this series, other four beams (originally named BN25, BN25D, BN50, BN50D) belonging to the extensive experimental program carried out by Podgorniak-Stanik [18], have been numerically analyzed.Series 25 and 50 differed from each other in terms of transverse cross-section dimensions, respectively equal to 300 mm x 250 mm and 300 mm x 500 mm, and in terms of net span, which was nearly doubled for the second one (Fig. 2b).Moreover, the two series were characterized by almost the same tension reinforcement ratio, which was equal on average to 0.85% for beams BN25 and BN50, and to 1.21% for beams BN25D and BN50D.These two last specimens contained indeed additional small longitudinal bars distributed along the web.More details about the geometric dimensions and reinforcement arrangement of the examined beams can be found in Fig. 2b.All the beams were cast by using concrete with a cylindrical compressive strength equal to 37 MPa.Mechanical properties of reinforcing steel and rebar details can be found in [18].All the tests were performed under loading control, by applying a central point load in several steps.At the end of each step, the load was lowered to approximately 90% of its current peak value and then increased again.During all test execution, crack pattern evolution and crack width were monitored in detail.

Numerical results vs. experimental observations
Numerical analyses have been carried out by implementing the above described 2D-PARC model into a commercial FE code (ABAQUS, [28]).Taking advantage of the symmetry of the problem, only one half of each beam has been simulated, by adopting a FE mesh constituted by quadratic, isoparametric 8-node membrane elements with reduced integration (4 Gauss integration points).Numerical analyses have been performed under displacement control, by applying an increasingly displacement at the loading point, in order to achieve a better numerical convergence and evaluate also the post-peak behavior.Numerical and experimental results have been first compared by considering the global response, in terms of applied load vs midspan deflection, as can be seen from Figs. 3, 4. In more details, Figs. 3 and 4a-b refer to the three specimens tested by Vecchio and Shim [17].On the same graphs, also the results obtained by Bresler and Scordelis [20] on nominally identical beams are reported.Fig. 3 shows a comparison between the experimental data relative to specimen OA2 and the numerical results obtained by adopting the 2D-PARC model, respectively implementing the non-linear constitutive relation described above or a simple linear elastic matrix for describing concrete behavior.These results highlight that the assumption of a constant value for concrete Young modulus during the analysis -equal to its initial value (E c = E ci )obviously provides a stiffer response, which is more pronounced for higher values of the applied load, where the effect of mechanical non-linearity is more significant.Moreover, the adoption of a linear elastic matrix for concrete leads in this case to a slight underestimation of the failure load.A more refined description of concrete behavior allows not only an improved modeling of the experimental response until failure, but also a much more stable solution, since the numerical analysis is characterized by a better convergence.On the contrary, both numerical analyses correctly predict the experimental shear failure mode and the corresponding final crack pattern.It should be here pointed out that the numerical response obtained through the original formulation of 2D-PARC model described in [7] is almost superimposed to that provided by adopting the concrete modeling proposed in this work.For this reason, the corresponding curve has not been plotted, so as to allow a better comprehension of the graphs of Fig. 3.In both cases, concrete is indeed treated as a non-linear elastic material subjected to a biaxial state of stress, but the here proposed approach is preferable since it requires a reduced computational effort and leads to an improved convergence of the algorithm.Figs.4a-b show similar comparisons for the other two specimens tested by Vecchio and Shim [17], namely OA1 and OA3 beams.In this case, experimental data have been only compared with the numerical response provided by the improved formulation of 2D-PARC model proposed herein, so to better appreciate its effectiveness.The graphs highlight the good capability of the model to represent the global behavior both at serviceability (cracking load) and at ultimate limit state.The peak load is accurately predicted, as well as the brittle shear failure characterized by no ductility.An accurate simulation of the experimental behavior, both in terms of stiffness and failure load, has been also obtained for specimens tested by Podgorniak-Stanik [18], as shown in Figs.4c-d.The adopted FE model is also able to predict that beams BN25D and BN50D (containing an additional distributed reinforcement on element sides, in addition to the main flexural bottom rebars) fail for higher shear stresses than the corresponding specimens BN25 and BN50, as provided by test results.Therefore, the satisfactory agreement between experimental and numerical responses proves that the proposed model is able to correctly describe the experimental behavior for different specimen geometries, as well as different reinforcement arrangements.Further comparisons between numerical and experimental results are also provided in terms of cracking development and crack widths, as depicted in Figs. 5 and 6.Fig. 5 shows the crack pattern at failure for the three specimens tested by Vecchio and Shim [17].As can be observed, the model exhibits a fine capability of reproducing the experimental diagonaltension crack, catching the very brittle and sudden failure typical of beams containing no shear reinforcement.Furthermore, maximum crack widths, which represent one of the most difficult parameter to predict in numerical analyses, are substantially comparable to experimental ones.Similar results have been also obtained for RC beams tested by Podgorniak-Stanik [18].More attention has been here devoted to the analysis of crack pattern evolution for increasing loads, in terms of crack distribution and widths.Some of the obtained results have been reported in Fig. 6 (for brevity, only for specimen BN50 [18]).As can be seen, crack patterns are reasonably well described during the entire loading history, highlighting that the proposed model can represent a powerful tool in the analysis of RC structures also at serviceability conditions, where crack control represents one of the fundamental issues to be checked.For the analyzed beam, the first numerical flexural crack occurred at approximately the same loading level registered during the experimental test.As proved by both numerical and experimental results, this first crack has been followed by the appearance of subsequent flexural cracks, until the attainment of the last loading stage, when a significant shear diagonal crack developed as extension of existing cracks.

CONCLUSIONS
n this paper, a constitutive non-linear model for the analysis of reinforced concrete structures, named 2D-PARC, has been modified so as to improve its computational efficiency, while maintaining its capability of providing correct predictions of the structural behavior.To this aim, the constitutive relation originally adopted in 2D-PARC for the modeling of concrete, both before and after cracking, has been properly substituted by implementing the one proposed by Ottosen, which is based on non-linear elasticity and is characterized by a high numerical feasibility.The accuracy of the proposed formulation has been tested through comparisons with experimental results on RC beams without stirrups failing in shear, which represent a quite difficult case study, particularly when modelled through two-dimensional membrane analysis.Based on the obtained results, it has been generally proved that 2D-PARC model is able to provide accurate predictions of strength, load-deformation response and failure mode.Moreover, also crack pattern evolution, both in terms of crack distribution into the considered element and in terms of maximum crack width, can be satisfactorily evaluated.Finally, it can be pointed out that the modular structure of 2D-PARC model makes it an attractive and versatile alternative approach for the analysis of RC structures, since the representation of each single resistant contribution can be easily changed when more refined and/or efficient constitutive relations become available.I

Figure 3 :
Figure 3: Comparison between experimental[17] and numerical results in terms of applied load vs. midspan deflection for specimen OA2.Numerical analyses have been repeated twice, by differently modeling concrete behavior.

Figure 4 :
Figure 4: Comparison between experimental and numerical results in terms of applied load vs. midspan deflection for specimens (a), (b) of the OA series[17], and (c), (d) of the BN series[18].

Figure 5 :
Figure 5: Experimental (left side, [17]) vs numerical (right side) crack patterns and crack widths at failure for specimens OA.