Aspects of Finite Element Modeling of Punching Shear Behavior of Reinforced Concrete Flat Slabs

A series of nonlinear FE analyses was conducted in this paper to model the punching shear failure of RC flat slabs. The analyses were carried out in the software DIANA using three-dimensional continuum elements. These analyses involve the test of different modeling choices in order to evaluate their influences on punching shear failure modeling. The parameters examined are the size of the FE mesh, the convergence methods and the concrete material input parameters. These parameters were investigated by comparing the results of load carrying capacity, load-deflection response and crack patterns from the FE analyses with a reliable experimental test. From the results obtained, it could be concluded that the numerical model was able to reproduce accurately the behavior of the tested slab and capture the punching shear failure.


INTRODUCTION
Punching shear failure is a complex phenomenon that occurs in reinforced concrete flat slab as a consequence of the stress concentration on slab-column region.A large number of experimental tests and analytical studies have been conducted on the last decades in order to determine the parameters that influence the punching shear behavior.Also, several models have been proposed, e.g.Kinnunen and Nylander (1960), Moe (1961), Nielsen and Braestrup (1976), Shehata (1990), Menétrey (1995) and more recently the Critical Shear Crack Theory (CSCT) proposed by Muttoni (2008).However, there is a still discussion about which one is the most accurately approach for the design of flat slabs.Besides, most of the design codes recommendation related to punching shear are based on empirical methods, e.g.Eurocode 2 ( 2004), ACI 318-14 (American Concrete Institute, 2014), and provided different methods that may differ significantly.
On the other hand, with the development of computational facilities, nonlinear finite element analysis (NLFEA) have become an interesting alternative for the assessment of concrete structures.These analyses may account nonlinear behaviors as concrete cracking and reinforcement yielding.Also, numerical parameter studies can be conducted to study the sensibility of each variable.In the case of modeling RC-flat slab, different investigations have been conducted.Early numerical studies implemented two-dimensional (2D) models using axisymmetric elements, e.g.Hallgreen (1996) and Menétrey (1995) for modeling circular slabs with ring reinforcement.However, this same approach was not suitable for slabs with orthogonal reinforcement.For that reason, tridimensional elements were adopted by other authors in order to represent in a more realistic way the punching shear behavior, e.g.Mamede et al. (2012), Belletti et al. (2015) and Shu et al. (2016).Those mentioned works have in common the use of the Smeared Crack approach for take into account the concrete cracking.The mean feature of this model is that the cracked concrete is assumed to remain as continuum material and the mechanical properties are modified to account for the effect of cracking, according to the evolving states of strain and stress (Cervera and Chiumenti, 2006).Despite of all the numerical investigations, modeling shear failure modes in slabs elements is still a challenging task using NLFEA models, mainly because of the lack of a general methodology and the quantity of input parameters.Few guidelines about nonlinear finite element modeling can be found in the literature (fib-Bulletin 45 (CEB-FIP, 2007), Hendriks et al. (2016)) but no one deals with the particular case of punching shear in flat slabs.Recently, Shu (2017) presented a "Multi-level Assessment Strategy" aiming to predict the punching shear behavior.This study included the evaluation of different modeling alternatives, involving meshing and input material parameters.
In this paper, a series of NLFEA were carried out in order to development a methodology for assess punching shear capacity in slab-column connection.It was selected a slab without shear reinforcement tested by Ramos (2003).Different model alternatives were analyzed, including meshing, input material parameters and also convergence techniques.The obtained results and the comparison with the experimental ones allowed to evaluate the accuracy of the development methodology.

EXPERIMENTAL MODEL
The studied flat slab was tested by Ramos (2003) and it had dimensions of 2300×2300×100 mm, supported by 8 steel plates with 100×100×50 mm and a square column of 200×200 mm located in its center.According to the nomenclature used by Ramos (2003), the slab was named AR2.The flexural reinforcement on the tension face of the slab was provided by 10 mm diameter steel bars in each direction with a spacing of 60 mm.The average yield strength was 639 MPa and the effective depth was 80 mm ( = 0.0169).On the compression face, the reinforcement was provided by 80 mm with a spacing of 200 mm approximately.The average value of concrete strength was measured using cube specimens.For this slab, the concrete strength was 48.9 MPa and the cylinder compression strength was considered of 80% of that value ( = 39.1MPa).Figure 1 illustrates the studied slab in this paper.Full details of the reference model could be found in Ramos (2003).
During the test the displacement on the top surface of the slab were measured using LVDT's.Also, the deformation of the bars close to the column was measured.The cracking pattern was monitored at the top of the slab.The slab had a punching shear failure with an ultimate load of 258 kN.

FINITE ELEMENT MODELING
In this paper, the finite element software TNO DIANA 9.4.4 was used to model the slab (DIANA, 2012).The main objective was to analyze the model behavior varying some properties of the mesh discretization, the convergence parameters and the nonlinear material constitutive models.For all the analyses only one quarter of the slab was modeled due to the symmetry.The loading and support plates were also included in the model.To simulate properly the behavior between the steel plates and concrete, interface elements were also included.The vertical displacement at the supports was fixed.In addition, on the symmetry axis, the perpendicular displacements were also restrained.A typical model with the boundary conditions is shown in the Figure 2. The bars were modeled using standard reinforcement embedded elements, which does not take into account the dowel effect and perfect bond between concrete and bars is assumed.The constitutive behavior of reinforcement is modeled using an elasto-plastic stress strain relationship without hardening.Finally, the loading was applied using an incremental displacement control through the central steel plate.

Baseline analysis
The concrete constitutive model and the convergence parameters used for the baseline analysis are listed in Table 1 and Table 2.These properties were used to carry out a mesh density study as will be shown in the next section.Also, a series of analyses were performed for examine the influence of the concrete properties and convergence parameters.In those analyses, one parameter of the mentioned tables was varied at the time with all the other parameters held constant.All the NLFEA were conducted in the same personal computer with the following configuration: Intel(R) Core(TM) i7-7700 CPU @ 3.6GHz Windows 10 Operational System 64bits.
Latin American Journal of Solids and Structures, 2018, 15(10), e120 4/22  The following sections will discuss the influence of each parameters analyzed in this work.

MESH DISCRETIZATION
The first part of this work was executed to study the influence of mesh discretization in the numerical solution.The main objective was analyze the effect of the mesh dimension in the accuracy of the model and the computation process time.Also, it was compared the behavior of the eight-node and twenty-node isoparametric solid brick elements, named in the software as HX24L and CHX60, respectively.The properties of eight different meshes tested are shown in the Table Latin American Journal of Solids and Structures, 2018, 15(10), e120  5/22 3. The Figure 3 illustrates the HX24L and CHX60 used in the mesh of the numerical models.The concrete material and convergence parameters were the same presented in the Tables 1 and 2. Figure 3. Finite element types and nomenclature used in Table 3.
The Figure 4 shows the comparison of load-displacement curves from the meshes A to F. The mesh type A analysis was performed in order to investigate the influence of the element selection.In this case, the HX24L (eight nodes brick element) was selected.The predicted load obtained was 209.374 kN (0.81P , where  is the load-carrying capacity of the experimental test by Ramos, 2003).The load-displacement curves for the Mesh A and D are shown in the Figure 4a.The behavior of both curves is similar, however the Mesh A response did not reach convergence after 14.5 mm.In addition, the crack pattern was less clear in comparison with other analyses in which were selected the CHX60 element (see Figure 5).Due to the poor results obtained, the HX24L element was dismissed for the others analyses.Latin American Journal of Solids and Structures, 2018, 15(10), e120  6/22 The analysis of meshes B to F were carry out to study the relation between the mesh density and accuracy of the numerical solution.The idea was to prove that upon a certain mesh refinement, the structural response remains the same.A course mesh size was proposed with only two layers of elements on the slab thickness (Mesh B).The number of elements was gradually increased until obtain the finest mesh (Type F, see Table 3).The Figure 6 shows the comparison between mentioned meshes.The load-displacement curves for the mesh B to F are presented in the Figure 4b.It was verified that the numerical results converged as the mesh was refined, therefore is demonstrated that a mesh independency solution is obtained through the use of the Smeared Crack approach with Total Strain Fixed formulation.The Table 4 shows the elapsed time in each analysis.As expected, the Mesh F required more time; this is explained by the number of elements and nodes, which means a larger system of equations to be solved.Nevertheless, it is important to notice that the increased of the ultimate load capacity (P) is small in comparison with the Mesh E.
Latin American Journal of Solids and Structures, 2018, 15(10), e120 7/22 The crack pattern obtained from the meshes B and F are show in the Figure 7.As expected, the coarse mesh has a less clear pattern, however it was possible to see the punching cone failure at the last convergence step at the mid cut section.The crack pattern obtained from the Mesh F allows differentiate the tangential and radial cracks on the top surface of the slab.
The last study related to mesh discretization was conducted to analyze a variable mesh size distribution.The aim was to reduce the number of elements and nodes of the mesh without affecting the accuracy of the numerical solution.An optimized discretization was proposed based on the layout used in the meshes D and E (see Figure 8a).Around the column, at a distant of three times the slab thickness, the mesh size was kept the same and outside that region was adopted a bigger element size as it is illustrated in Figure 8b.Two configurations with 4 and 5 layers of elements were tested (Mesh G and Mesh H respectively). Figure 9 shows the load-displacement curves obtained from the analysis.It can be noted that in both cases the numerical response between the original mesh and the optimized one is approximately the same.
Latin American Journal of Solids and Structures, 2018, 15(10), e120 8/22 In addition, the Table 5 shows that the elapsed time of each model was reduced.For that reason, was verified that the optimized meshes achieved the objective of reducing the computation time process without affecting the accuracy of the numerical solution, both for the Mesh G and for the Mesh H.For next analyses described in this paper the mesh H was adopted.

MATERIAL PARAMETERS STUDY
The material input parameters are important choices in order to provide a numerical analysis behavior as experimental tests.The overall aim of this parametric study was to investigate their sensibility of some parameters and its influences on the punching shear failure.This study includes the influence of the Young's modulus, fracture energy, tension softening behavior, tensile strength and shear retention factor.
This parametric study was conducted by varying the main parameters of the concrete and investigate its sensibility on load-deflection relationship and compressive strain behavior.All parameters mentioned were analyzed one at the time, and the other parameters remain as described in the baseline section of this paper (see Tables 1 and 2).Regarding the results, the load-deflection relationship describes the behavior of the model in the loading procedure.Also, the compressive strain behavior provides information regarding punching shear failure of the tests.
Figure 10 shows the general behavior of the mentioned load-radial compressive strain curves.According to Shu (2017), when the load applied in the slabs reaches 0.7 from the load-carrying capacity, the compressive radial strain by the side of the column reaches a peak and begins to decrease, indicating a change of internal forces.On the other hand, Guandalini (2006) assumes this radial strain peak at 0.8-0.9P .In general, from the results of this paper, it was observed this internal force change above 0.9P .Figure 12 shows the results from the materials input parameters that are associated to shear behavior.These parameters provide a relevant change in the response of the model in the loading procedure.It was observed a flexural failure behavior of the test by increasing the constant shear retention factor.This can be observed in the  = 0.2 curve (see Figure 12a and Figure 12b).The increase of the shear retention factor provides the model an increase of the punching shear strength leading the structure to a flexural failure.
Figure 12 (continued).Results from the materials input parameters associated to shear behavior.
The reduction of the shear retention factor () does not achieve the load-carrying capacity observed in the experimental results (see Figure 12a,  = 0.01).The model with a smaller shear retention factor does not consider the aggregate interlock on crack opening phenomena leading to a smaller load-carrying capacity.
This also was observed in the comparison between the Fixed and Rotate crack model (Figure 12c and 12d).The orientation of the cracks in the rotate smeared crack model changes according to the principal strains in each loading step, leading a crack surface only with normal components.As a consequence of the absence of shear stresses in the failure surface, an inaccurate load-carrying capacity that does not reach the experimental results could be observed (see Figure 12c).In the fixed approach, the crack orientation is also according to the principal strain, but once developed, the crack orientation remains the same, not depending of the principal strains in the following steps.
Figure 13 shows the results from the material parameters related to compressive behavior of the parametric analysis.By comparing the results from confinement and compression reduction due to lateral cracking (Figure 13a and 13b), it observed that the model in which these effects are not included reaches a smaller load-carrying capacity and strain peak in comparison with the reference model.According to Selby and Vecchio (1997), the confinement phenomena on concrete provides an increase of the compressive strength and strain peak.This leads to an increase of the load-carrying capacity observed in Figure 13a.Besides the confinement, the compressive fracture energy was another analysis related to compressive stresses.From the results observed in Figure 13c and 13d, the fracture energy does not provide a significant change in the behavior of the radial strain.This also was observed in the load-deflections curves.
According to the formulations related to the confinement and compressive fracture energy, both parameters provide a change in the ultimate compressive strain.In addition, the confinement model provides an increase of the compressive strength in order to remain the same compressive fracture energy value along the loading procedure.This leads to conclude that the slabs tested did not reach the final compressive strain of the constitutive law along the loading procedure, reaching the same strain peak.
The same compressive strain peak mentioned was associated with the tensile strain from the equilibrium procedure.The modification of the tensile behavior of the concrete provides the slab to achieve a different equilibrium state and, consequently, different compressive strain peak.This can be observed in the softening behavior analysis (see Figure 14a  and 14b).The linear softening modeled with the same tensile strength and tensile fracture energy has a different ultimate tensile strain in comparison with the reference Hordijk softening model.This leads to an equilibrium state that provides an increase of the compressive strain peak.Regarding load-carrying capacity, it was observed that the linear softening model reaches almost the same results in comparison with the reference Hordijk softening model and the curve of loaddeflection relationship present a good agreement each other (See Figure 14a).
Figure 14c and 14d illustrates the results from the tensile strength parametric analysis, it was observed that the tensile strength influences in the load-deflection relationship of the tests.An increase of the tensile strength provides an increase of the load-carrying capacity of the tests.However, the compressive radial strain peak remained almost the same indicating that the tensile strength does not influences in the results from the compressive strain.
In relation to the tensile fracture energy, the formulation of CEB-FIP MC 1990 leads to a tensile fracture energy almost 50% smaller than the fib MC 2010 (0.77 N/mm against 0.141 N/mm, respectively).The tensile fracture energy was directly proportional to the load-carrying capacity obtained from the model.This also conducts the model to a ductile behavior reaching larger strains at the end of the loading procedure.
Figure 14.Results from the materials input parameters associated to the tensile behavior.
On the other hand, a different result was observed in the Young's modulus parametric study (see Figure 15a and 15b).In general, an increase of the Young's modulus in comparison with the reference model provides a decrease of the loadingcarrying capacity.In addition, it was observed that a higher Young's modulus provides a smaller strain by the side of the column and, consequently, a brittler behavior by analyzing the compressive peak strain in Figure 15b.

CONVERGENCE ANALYSIS
In order to investigate the influence of the convergence parameters in the behavior of the studied slab, a parametric study was performed.The results from the convergence study was conducted to evaluate the elapsed time until the convergence of all the loads steps, number of iterations converged and load-displacement response.The material parameters from the baseline section remained the same for all analysis presented in this section.First of all, the influence in the behavior of the numerical response of the iteration methods Regular Newton-Raphson, Modified Newton-Raphson and the Secant Method were investigated.In this first analysis, the convergence criteria was based on force, energy and displacement simultaneously and the displacement load was applied in 200 steps with the size increment step of 0.2 mm, see CM01, CM02 and CM06 in Table 6.
In the second analysis, the Secant method was chosen to carry out the numerical model varying the convergence criteria by energy, force and displacement one at the time.The incremental size of the load step model remained 0.2 mm applied in 200 steps, see CM03, CM04, CM05 in Table 6.
CM03, CM07 and CM08 were studied in order to investigate the influence of the load step incremental size on the numerical analysis.The Secant convergence model was chosen in CM03, CM07 and CM08 based on energy convergence criteria with a tolerance of 10 .CM03 was carried out with a load step incremental size of 0.2 mm, in order to CM07 and CM08 have its load step incremental size of 0.5 mm and 1.0 mm, respectively.
Then, CM03, CM09 and CM10 were carried out to evaluate the response of the model by varying the tolerance of the convergence criteria.All these FE models use the Secant iterative method based on energy criteria.CM03, CM09 and CM10 were carried out with a tolerance of 10 , 10 and 10 , respectively.See Table 6 for more details about the convergence parametric study.All the other numerical, physical and geometrical parameters remained the same as described in section 3.1 (Baseline analysis).The maximum number of iterations per steps for all models were 200.From the analysis of the influence of the iteration convergence method, Figures 16 shows the results obtained from the tests.Figure 16a illustrates the load-deflection curves obtained from the analysis of the CM01, CM02 and CM06 models and the Figure 16b shows the step number versus the number of iteration required to reach the convergence of each step.
Figure 16.Results from convergence method analysis.
From the comparison of the load-deflections curves show in Figure 16a, the models stiffness has a good agreement but the load-carrying capacity and the final displacement reached are not the same.The load-carrying capacity reached by the Secant Method was 261.64 kN (1.01 , where  is the load-carrying capacity of the experimental test by Ramos, 2003), against 244.82 kN (0.95 ) from the Regular Newton-Raphson Method and 234.39 kN (0.91 ) from the Modified Newton-Raphson Method.
However, Figures 16b shows that the Secant Method requires more iterations per step than the other methods studied.According to the mathematical formulations, the Regular and the Modified Newton-Raphson methods reach the solution by changing the stiffness matrix of the nonlinear problem based on mathematics expression that do not take into account the previous solutions steps, while the Secant Method changes the matrix based on its current step.This requires more iteration per steps and, consequently, more time to reach the solution of each step.
The Modified Newton-Raphson changes the stiffness matrix only in the first iteration of each step, in some cases this could lead to a less computational work and less time to reach the solution on each step.However, it is important to remind that more iterations are needed to solve the problem when compared to Regular Newton-Raphson.The Modified Newton-Raphson Method reaches the last convergence step in 38 minutes, while The Regular Newton-Raphson Method and The Secant Method reaches in more than 1 hour.Nevertheless, by using the same matrix in all the iterations may cause problems in the convergence procedure, not solving the problem and reaching the maximum number of iterations.The number of the steps that converged using the Modified Newton-Raphson Method was 85, while the Regular Newton-Raphson Method was 95 and Secant Method was 111.This leads to a maximum displacement reached by the model based on the Modified Newton-Raphson Method of 16.8 mm, while the experimental results by Ramos (2003) was about 20 mm.In the tests presented in this paper, the Secant Method was able to best represent better the experimental results in comparison with the other iteration convergence method based on its load and displacement results, reaching 261.64 kN (1.01 ) and 22 mm (1.1 ).
The second analysis of the parametric convergence study was conducted to evaluate the Secant Method based on different convergence criteria.The criteria studied in the test was based on Energy (CM03), Force (CM04) and Displacement (CM05).Figures 17 shows the results from this analysis.The Figure 17a shows the load-displacement relationship and Figure 17b show the number of the steps that converged versus the number of iterations needed to reach the convergence.Regarding load-deflection relationship, Figure 17a shows that the test based on energy and force criteria had similar results.The results from the Secant Model based on displacement criteria did not have a good agreement with the other models studied in this paper, including the experimental results from Ramos (2003).From the displacement of 15 mm to the last step presented in Figure 17a, the model based on displacement criteria presented an instability in their results and an untrue equilibrium (see Figure 17a).
Figure 17b indicates that test based on displacement needed less iterations to converge each step.This can be observed in the first steps of the analysis and give a good agreement with the tests based on Force and Energy (see Figure 17a).After the 80th step (16 mm, Figure 17a) the number of iterations needed to converge each step and the solution presented were not according to the expected results from a punching shear failure, for that reason, the FE analysis of the model based on displacement criteria was interrupted by the authors in the 115 step.From these results, it could be concluded that the displacement criteria should be avoided.This recommendation is also suggested by Hendriks et al. (2016).This fact is explained by the use of incremental displacement vectors in the displacement criteria that could leads to an untrue equilibrium when small values of those vectors are reached between iterations (Claus, 2009).
The next study was focused on the convergence tolerance.Too loose norm criteria may lead to inaccurate and unreliable results, on the other hand, a too strict norm will lead to an excessive computational work, without improving on the results (Claus, 2009).The tolerance value must be determined carefully, independent of the criteria chosen to carry out the nonlinear analysis and its strictness.
In order to investigate the differences of a loose and a strict analysis, the results from CM03, CM09 and CM10 were compared.The models differs each other only by the tolerance of the convergence criteria, using the Secant Method to solve the nonlinear problem (See Table 6).Figure 18 illustrates the results obtained.Figure 18a illustrates the load-deflection relationship and Figure 18b the number of steps converged versus the number of iteration needed to reach the equilibrium.
Figure 18.Results from the convergence tolerance analysis.
Latin American Journal of Solids and Structures, 2018, 15(10), e120 16/22 From the results displayed in Figure 18, all the models were similar each other in terms of stiffness on the loaddeflection relationship.Nevertheless, the model with the tolerance of 10 (CM10) had some divergences along the load incremental steps.Figures 18b show that the model with the tolerance of 10 (CM10) needed less iterations than the other models studied.The loose tolerance of CM10 leads the analysis to have an oscillating behavior in the load-displaceresponse along the steps.
Comparison between the CM03 (tolerance of 10 ) and CM09 (tolerance of 10 ) presents another inaccurate solution along the incremental load steps in the analysis of the model.From the CM09 results (Figures 18b) it was observed that the number of iterations to reach the equilibrium in each step increases with the number of incremental steps converged, reaching the maximum number of iterations in the 83th incremental step (16.6 mm) of the analysis and achieved 0.9 .
In the other hand, CM03 need more iterations to reach the equilibrium on each step, but the equilibrium obtained it is a more reliable solution and the convergence of the subsequent steps leads to a final behavior more accurate than the other models.Regarding the time elapsed to finish the analysis, CM09 and CM10 elapsed about 25 minutes against 54 minutes by CM03.
In order to analyze the influence of the size of the load incremental steps of the models, CM03, CM07 and CM08 were carried out (see Table 6 for more details).Figure 19 illustrates the results from the analysis of the models.Figure 19a shows the results from load-carrying capacity versus max deflection of the model and Figure 19b shows the steps converged normalized by the total number of steps reached versus the number of iterations needed to reach the equilibrium of each steps.
CM03 was the reference model in this analysis.The size of the load incremental step carried out by CM03 was 0.2mm, reaching 40mm after convergence of all the 200 load steps.CM07 and CM08 were carried out with a larger incremental step size: 0.5mm and 1.0mm, respectively.In order to compare the models each other, the x-axis that corresponds of the incremental steps on Figure 19b was normalized by the maximum number of steps converged reached in the end of the FE analysis.
Figure 19.Results from the load incremental step size analysis.
Regarding load-deflection relationship, all the models studied to investigate the size of the incremental load steps are in a good agreement each other.The load-carrying capacity reached by the models are between 258kN and 265kN and the maximum deflection are between 21.2 mm and 23 mm (see Figure 19a).
The major difference on the results are the number of steps converged and the time elapsed on the analysis.CM03 (reference with 0.2 mm/step) needed 54 minutes to reach the final solution of the test against 29 minutes from the CM07 (0.5 mm/step) and 22 minutes from CM08 (1.0 mm/step).It also was observed that a larger incremental load steps size will lead to a higher number of steps needed to reach the equilibrium of each step, see Figure 19b.
In general, a small incremental size will provide a more stable and accurate results along the analysis procedure (Claus, 2009).Furthermore, a smaller incremental step size on the nonlinear procedure converges more easily.This can be observed in the number of iterations need to converge each steps in the Figure 19b.However, this leads to a large computational expense, observed in the large amount of steps converged and in the comparison of the time elapsed to finish the analysis of CM03 (109 steps and 54 minutes) and CM10 (25 steps and 22 minutes).
In reinforced concrete problems, a hybrid of large and small incremental load size may be a good alternative choice to conduct an extensive nonlinear analysis.A larger incremental size may be used in the beginning of the analysis, covering Latin American Journal of Solids and Structures, 2018, 15(10), e120 17/22 the elastic linear behavior and the small incremental size may be used after.This will provide more accurate results in the post-cracking behavior stage and a less computational expense in the linear elastic behavior stage of the model.The last parametric analysis studied in this paper was the influence of the line search technique on the steps convergence procedures.This analysis was compared by CM03 and CM11, which are the same numerical model but only with the difference of the CM03 was carried out with line search technique and CM11 was carried out without it.Figure 20 illustrates the comparison from CM03 and CM11.From the results shown in Figure 20, there not significant differences observed in the models.The load-carrying capacity reached is almost the same value (258 kN to CM03 against 257 kN to CM11) and the elapsed time differs only in 2 minutes (54 minutes to CM03 against 56 minutes to CM11.The line search is a mathematical technique to provide a solution in less time, furthermore the secant method has a search technique included that provides almost the same results, justifying the observed results.To summarize, the results from the convergence parametric study are presented in Table 7.

NUMERICAL AND EXPERIMENTAL RESULTS COMPARISON
The obtained results from the NLFEA were compared with the experimental results in order to validate the parametric study conducted in this paper.The mesh H, proposed in the section 4, modeled with the parameters implemented in the baseline model were used for this comparison.
Latin American Journal of Solids and Structures, 2018, 15(10), e120 18/22 7.1 Load-displacement curve The Figure 21 shows the load-displacement comparison for the Positions 1 and 2. It is important to remind that in the experimental scheme the load was applied in the plates around the column and the column itself was the only support of the slab.Thus, in order to compare the numerical results, the model reference system was changed to coincide with the test scheme.From the Figure 21 it can be observed that the load-displacement behavior of the experiment was similar in comparison to the numerical analysis results, especially at the position 1.In the case of the position 2 a stiffener behavior was obtained, however the general behavior of the experimental load-displacement curve is represented properly.

Load-strain curve for reinforcement bars
A good adjustment between the numerical and experimental results was found for all the bars measured (See Figure 22).At the beginning of the loading stage, the bars had a linear behavior that was correctly captured with the numerical model.In the last stages a yielding behavior was observed.In the reinforcement bar 1 (Figure 22a), was verified a huge increment of the strain, that was not possible to be represented in the numerical curve.For all the other bars, a better adjustment was obtained, even in the last stages of loading.
Figure 22 (continued).Comparison from load-strain on reinforcement bars between numerical models and experimental test.

Cracking Pattern
The crack pattern evolution is shown in the Figure 23 at three different stages.This evolution allows to examine the development of the radial and tangential crack around the column.Also, the post-failure crack pattern at the cross section is compared with the experimental results (See Figure 24).The numerical crack pattern shows that the model was able to capture the punching shear failure.It can be observed that the inclined shear crack have a similar angle and position in comparison with the experimental punching cone.

CONCLUSIONS
In this paper, several choices on NLFEA of punching shear slabs problems were carry out in order to investigate theirs influences in the overall behavior of the model.The proposed study aimed to investigate a variety of mesh discretization of the FE model in the accuracy of the results, the influence of the main RC materials input parameters and convergence parameters to reach the equilibrium.From the result obtained the following conclusion can be drawn: Based on the analysis of the mesh discretization, it was verified that the numerical solution converged as the mesh was refined.By analyzing the result of the different configurations tested, it is recommended the use of twenty-node brick element with quadratic interpolation and a minimum of 5 layers of elements along the slab thickness.In order to reduce the number of nodes and elements, a less fine mesh could be used outside the punching failure zone, without affecting the accuracy of the numerical solution.
Regarding the materials input, it was observed that the parameters that most influence the behavior of the slabs was the parameters associated to shear.The shear retention factor showed a relevant importance in the failure mode and the fixed smeared crack model presented a good agreement with the experimental results.This can be justified by the shear stresses that developed on the crack surfaces.The parameters associated to compressive and tensile behavior provides changes in the radial strain and in its final equilibrium state.
According to the results from the material input parametric analysis, the Young's modulus was inversely proportional to the load-carrying capacity and compressive radial strain peak.On other hands, the tensile strength, the shear retention factor and the tensile fracture energy were directly proportional to the load-carrying capacity.The compressive fracture energy does not provide significant changes on the overall behavior.Nevertheless, the confinement provides an increase of the compressive strength along the loading procedure, consequently increasing the load-carrying capacity.
From the convergence analysis, the Secant method provides more computational expense and a larger elapsed time to reach the final solution in comparison to the Regular and Modified Newton-Raphson.However, The Secant Method provides less convergence problems in comparison to the other methods.The displacement criteria demands less iterations to reach the convergence on each step, but it could leads to an inaccurate solution.
Based on the analysis of this paper, to reach a strict solution without problems of convergence and untrue equilibrium, it is recommended a Secant method using simultaneously force, energy and displacement criteria with a small load incremental step size.A Secant Method based by force criteria and a hybrid of a small and larger load incremental size may be used for a fast preliminary solution.

Figure 1 :
Figure 1: Dimensions and rebar layout of the experimental slabs.Adapted from Ramos (2003).

Figure 2 .
Figure 2. Typical mesh of the model with the boundaries conditions.

Figure 4
Figure 4. a) Comparison of load-displacement curves of A and D meshes.b) Load-displacement curves obtained from the mesh refinement analysis (Meshes B to F).

Figure 5 .
Figure 5.Comparison of crack pattern (principal tensile strain  ) obtained from the mesh A (left) and mesh D (right) models, both meshes had the same size but different type of element.

Figure 6 .
Figure 6.Discretization of a) coarse mesh B and b) fine mesh F.

Figure 7 .
Figure 7.Comparison of crack pattern (principal tensile strain  ) obtained from the mesh B (left) and mesh F (right) models.

Figure 8
Figure 8. a) Typical element size used for meshes D and E b) Optimized size element distribution used for meshes G and H.The meshes D and G had 4 layers of elements on the slab thickness while the E and H had 5 layers.

Figure 9 .
Figure 9. Load-displacement curves obtained from the mesh optimization analysis a) Comparison between the meshes D and G (4 element layers on slab thickness) b) Comparison between the meshes E and H (5 element layers on slab thickness).

Figure 10 .
Figure 10.Position where was measured the radial strain on the numerical model (left) and typical load-radial strain (right).

Figure 11 .
Figure 11.Varied Material input parameters.The bold alternatives are the baseline values.
Figure 17.Results from the convergence criteria analysis.

Figure 20 .
Figure 20.Results from the line search analysis.

Figure 21 .
Figure 21.Comparison from load-deflection relationship between numerical model and experimental test from Ramos (2003).

Figure 23 .
Figure 23.Evolution of the crack pattern in the numerical model.

Table 1 :
Baseline materials input parameters.

Table 2 :
Baseline convergence parameters

Table 3 :
Mesh configurations and their properties.The nomenclature used for the element dimension (abc) is illustrated in the

Table 4 :
Results obtained from the Meshes B to F varying the element size.

Table 5 :
Result obtained from the mesh optimization analysis.