Finite element model and test results for punching shear failure of RC slabs

The nonlinear finite element analysis of pin-supported reinforced concrete slabs of moderate thickness is the main subject of this paper. This is an important issue of the mechanics of concrete structures. Thus, a nonlinear FE analysis of RC slab models subjected to punching shear was initiated and the comparison of the numerical and test results was made. The prediction engine of the crack pattern due to bending and extension is incorporated in the Mindlin-type (moderately thick) slab model together with other nonlinear features. This is in order to formulate an alternative model of RC slab in relation to the layered models or fully three-dimensional models. New formulas were applied for 3D constitutive relationships for concrete and for tension stiffening effect. Prediction of punching shear was facilitated by Podgórski’s failure criterion for concrete. On the one hand, a considerable advantage of the proposed approach is a relatively low numerical effort in comparison with the existing models, while on the other hand the applied model clearly describes the physical behaviour of a real slab. A supporting test programme for validation was run. Three RC slabs with a system of double-headed studs as the reinforcement against punching were tested by the authors in ITB Strength Tests Laboratory. The results for the units constructed as square slabs with a central short column subjected to full-scale tests were initially compared with assessments based on standard provisions and technical approvals. As a result of this approach, the overall prediction of the nonlinear behaviour of the test model, including the model of shear failure, is in accordance with the experimental data.


Introduction
Column-supported floor slabs are encountered in many types of structures, including RC tanks. Finite element modelling of punching shear failure of RC slabs is an important issue of the mechanics of concrete structures, and it was the subject of hitherto published papers of many authors. Abbasi et al. [1] utilised the layered model and analysed failure modes taking into consideration punching shear as a function of the reinforcement ratio and giving appropriate criteria of limit states. A series of important analyses with iterative FE algorithms were presented at the conference International Workshop on "Punching Shear Capacity of RC Slabs" held in Stockholm in June 2000. Polak and Guan [2,3] also used the layered element formulation with transverse shear. Guan and Loo [4] used this approach to analyse slab-corner column connection as well as slab-column systems. In solving similar problems, Wosatko et al. [5,6], as well as Genikomsou and Polak [7], implemented 3D elements utilising the damage plasticity model of concrete. Menetrey and Willam [8,9] performed an interesting analysis of punching shear in reinforced concrete, setting new criteria for process localisation. Vocke [10] thoroughly analysed the punching behaviour of RC flat slabs in edge and corner columns by using the nonlinear FE program and with full-scale testing.
The analytical method proposed in this paper has no direct reference to the methods described above, and it consists of two separate algorithms: the iterated crack patterns method (ICPM) previously successfully applied for the analysis of thin RC slabs [11], and the fracturing surface localisation method (FSLM) based on the application of a three-dimensional limit state condition for concrete. The first is an intermediate approach between the modified stiffness 1 3 36 Page 2 of 16 method and layered models of RC slab. This paper expands the previous investigations the authors carried out in scope of the FE analysis of thin RC planar structures [11] and the slabs of moderate thickness (allowing for the impact of shear deformation on the slab deflection) [12]. The finite element model for the nonlinear analysis of reinforced concrete slabs of moderate thickness is also an important subject of this paper. The prediction engine of the crack pattern due to bending is incorporated in the Mindlin-type slab model, together with other nonlinear properties, similar to the way presented by the authors of this paper at AMCM conference [12]. However, the novel formulae have been used herein for nonlinear constitutive relationships for concrete and for tension stiffening effect. The authors modified the constitutive model of concrete proposed by Kotsovos and Pavlović [13,14] to fit the material model for experimental data, which was important for the problem under consideration. The model was subjected to calibration procedures in the uniaxial and biaxial compression stress regime and was compared with the experimental data available in the literature, mainly with the test data of Kupfer et al. [15]. Additionally, in order to avoid overestimating slab capacity under brittle punching conditions, a special descending path was incorporated into the constitutive model for the range beyond the peak stress. The tension stiffening effect is taken into account, with the assumption that additional stress is carried by the reinforcement steel. In this paper, the corresponding model was derived according to the provisions of EC 2 [16] by using appropriate transformations. The constitutive matrix components are computed by means of numerical integration (incorporating Gauss quadrature) through the slab thickness. Ten different cracking patterns (all possible for slabs) are assumed in ICPM formulation. The use of each pattern implies the division of slab thickness into several layers (maximum three layers of concrete) and, consequently, a maximum of three sets of integration ranges. As a result, calculation of stiffness matrix elements takes a relatively small amount of computer time and, unlike in the case of commonly used layered model, free from the influence of an arbitrary chosen number of layers. The depth of the compression zone can be determined more accurately than in the case of the layered model. The effect of coupling between membrane deformations and flexure, which is caused by a non-symmetrical double reinforcement of the slab and cracking of the concrete is taken into account. As the approach analogous to deformation theory of plasticity is used for the material description and the stiffness matrix of Mindlin-type RC slab, a nonlinear analysis of RC slabs has been performed using the finite element direct iteration method. The failure mechanism of RC slabs (whether brittle punching or ductile bending) depends on the ratio of punching to bending strength. In the proposed model, both phenomena are described by separate regions of the element's constitutive matrix, where in turns shear or bending mode failure can be predicted by the material model of Kotsovos [13,14] and limit surface for concrete of Podgórski [17].
Extensive experimental research has been conducted on the phenomenon of punching shear in reinforced concrete floor slabs with reinforcement against punching in the form of double-headed studs. Studies on the implementation of this kind of reinforcement have been conducted by Eligehausen, Hegger, Beutel and Vocke [18][19][20], as well as Polak et al. [21]. The researchers from Silesian University of Technology [22][23][24] and Łódź Technical University [25] have also carried out widespread investigation on the subject. Starosolski et al. [26,27] performed some valuable tests on slabs reinforced with studs. In addition, Urban [28] and Hulimka [29] recently published monographic papers on punching RC slabs.
In order to run the test programme for validation, the authors of this paper examined three RC slabs (test specimens) with a system of double-headed studs as a transversal reinforcement in the ITB Strength Tests Laboratory in Warsaw [12]. The results of the tests on units constructed as square slabs with central short columns subjected to forcecontrolled loadings acting from the bottom were initially compared with assessments based on the requirements of different standards and technical approvals. However, the calculation results obtained this way are not the subject of this paper. The test results were then used to confirm the accuracy of the analytical method proposed in this paper.

Main steps of the nonlinear FE analysis
The main steps of the nonlinear FE analysis are as follows. The procedure begins from the solution of elastic, doubly reinforced concrete slab. The direct iteration is adopted for the nonlinear analysis in accordance with the total strain material modelling. Following the computation of nodal displacements, strains and stresses in concrete and steel, the limit state condition for concrete is checked. Cracks in concrete are introduced in the direction perpendicular to the axis of principal tensile stress at points where the cracking criterion has been reached. Cracking is taken into account by defining a certain cracking pattern. Parameters of integration are determined. For each loading step, the stress state is checked at every Gaussian point, and the material properties at these points are updated accordingly. A new stiffness matrix is formed, and iterations are carried out to ensure that the convergence checked by the displacement test is achieved. Further cracking of concrete is taken into account by checking a stress state following a certain cracking pattern, and such a pattern is updated if necessary. The previous cracking is stored for the next loading step. When the convergence is reached, a subsequent load case is used. Iterations to obtain a convergence are performed until each loading set is exhausted. The limit state condition for concrete (together with the imaging technique) sets the threedimensional layout of punching shear failure.

Failure criterion for concrete
As the ultimate surface for concrete, the failure criterion proposed by Podgórski [17] allowing for the useful simplifications proposed by the first author of this paper [30] was implemented in the concrete material model. The so-called JP limit state surface (Podgórski [17]) almost does not differ from the limit state surface of Willam and Warnke (W&W) (see [26]), under the assumption for this model that the parameter e = 0.514 (in the case of the data adopted in the numeric example). The parameter e is responsible for the appearance of sharp edges in the deviatoric cross section of the W&W failure surface. In the criterion assumed herein, the function r(φ), which describes the radius of the failure surface cross section by the octahedral plane, should depend on two parameters to fit the experimental data. Such a function can be written in the following manner: where α and β are constants fulfilling the following inequalities: 0 ≤ α ≤ 1 and 0 ≤ β ≤ π/6. Parameters describing the proportion of radii r for different φ angles can be written as: The radius (r), the distance from 0 on the hydrostatic axis (h) and the angle (φ) are related to the stress tensor and deviator invariants in the following manner: where where in turn I 1 = the first invariant of stress tensor, s ij = the deviator of the stress tensor (s ij = σ ij − σ kk δ ij /3), J 2 , J 3 = the second and the third invariants of stress tensor deviator, respectively, and σ o , τ o = the octahedral normal and shear stresses, respectively (both described as above). (1) The octahedral strains denoted as ε ο and γ ο are defined in a similar way: The plots of the r (φ) function described by Eq. (1) are presented in Fig. 1 (see also Figs. 2 and 3). It is difficult to distinguish the plots of JP (Podgórski) and Willam and Warnke (W&W) limit surfaces shown in the presented figures, because they overlap. According to Podgórski [17], α and β coefficients from Eq. (1) can be derived by means of iteration.
As this method of deriving α and β coefficients is not very effective, equations for the direct calculation of these parameters were proposed by Lewiński [25]. By substituting Eqs. (1) into (2) for φ = 0°, 30°, 60°, we get the set of three trigonometric equations. They can be solved with regard to parameters α and β: J. Podgórski confirmed the usefulness of the above expressions [27]. Assuming that in the octahedral set of coordinates the meridian of ultimate surface is described by a parabolic equation, the following equation for the failure surface for concrete was used (see Fig. 3): The values of C 0 , C 1 , C 2 , α and β coefficients can be derived on the basis of experimental data, utilising the following strengths of concrete: uniaxial tensile strength f t , uniaxial compressive strength f c , biaxial compressive strengths f cc (σ 1 /σ 2 = 1) and f 0c (σ 1 /σ 2 = 2), triaxial tensile strength f t . According to the tests performed by Kupfler et al. [15] f cc = 1.1 f c and f oc = 1.25 f c were assumed, while the coefficient κ = f t /f c was assumed to be equal to 0.1. Both λ and ϑ parameters can now be expressed in the form of equations dependent on the tests results: The coefficients C 0 , C 1 and C 2 from Eq. (8) now can be written in the following way: = sin 3 0 , where (10)

Material model for concrete
Concrete behaviour is described by the mathematical model proposed by Kotsovos and Pavlović [13]. The incorporated model is based on the three-modulus approach modified on the basis of the internal-stress concept. Two secant moduli, bulk K cS and shear G cS , are functions of the deviatoric and volumetric components of stress state: where K 0 and G 0 are bulk and shear initial moduli, respectively, while A, b, C and d are model parameters depending on f c value (see [13]). Instead of the coupling modulus, an equivalent superimposed stress state is introduced: Again, k, l, m and n parameters can be expressed in terms of f c on the basis of experimental tests. Finally, the stress-strain relationship can be expressed as: where σ id accounts for the coupling between stress deviator and volume change. The above relationships do not describe the effect of dilatancy (Kelvin effect) which involves volume increase (without the presence of negative hydrostatic pressure) after reaching a specific level of octahedral shear stress τ oOUFP -at the onset of unstable fracture propagation (OUFP), described below for the present approach by Expression (18). This volume increase can be described by the increase of octahedral normal strain (δε o ) equal to [14]: where the increase of octahedral shear stress (δτ o ) is equal to Instead of the application of the ultimate surface for concrete given by Kotsovos and Pavlović [13,14], the failure criterion proposed by Podgórski [17] was implemented in concrete material model, because, according to the authors, the latter model fits better with the experimental data. After simplifying Eq. (8), with the assumption that f t = 0.1f c , the following expression for the ultimate strength was obtained: On the basis of the above failure criterion and the experimental results, the surface separating stress state that does not cause dilatation phenomenon (OUFP) has been proposed in the form: The octahedral shear stress (δτ ο ) and strain (δε o ) increments were obtained by the iterative numerical analysis. In order to verify the model, the comparison of the numerical analysis results and experimental data obtained by Kupfer et al. [15] was made. For the value of mean compressive strength adopted for analysis: f c = 33.6 MPa, the values of the model parameters are as follows: The assumed constitutive model was tested for f c = 32.4 MPa and the subsequent relations of stresses in concrete: σ 2 = 0, σ 2 = 0.226σ 1 , σ 2 = 0.226σ 1 , σ 2 = 0.525σ 1 , σ 2 = σ 1 and also for f c = 19.1 MPa and the same stress relations. An example of stress-strain relationships for the proposed model for the relation σ 2 = 0.226σ 1 is presented in Fig. 4. The grey lines represent experimental data, the black ones-calculation results. In these comparisons, the value of σ 3 principal stress is equal to 0 and the compressive strength of concrete f c = 32.4 MPa. The constitutive model of concrete proposed by Kotsovos and Pavlović [13] has been modified to fit the model to the experimental data valid for the problem of brittle failure. In order to avoid the overestimation of slab capacity in the conditions of brittle punching, a special descending path was incorporated in the material model for the range beyond the peak stress.
A more recent work of Zisopoulos et al. [28] on the behaviour of concrete under different boundary conditions confirms that concrete can be described as a brittle material. Since, as the frictional forces at the sample base reduce, concrete suffers a faster loss of its load carrying capacity than it could result from test data obtained by Kupfer et al. [15].

Models for reinforcing steel and tension stiffening
A hypo-elastic model of reinforcing steel has been assumed for the purpose of numerical analysis because of the adoption of such constitutive models for concrete. The stress-strain relationship for this model is shown in Fig. 5. The coefficient ψ spi has the following form for The concept of modelling the tension stiffening by the modification of reinforcing bars' stiffness is not quite new [29]. The influence of the reinforcement bond in cracked concrete on slab deformation has been accounted for by the tension stiffening approach according to the provisions of Eurocode 2 [16]. Expression (7.18) given in this standard can also relate to the strains in reinforcing steel in the arbitrary i direction and take the form: where ε smi = the mean strain in the reinforcement (according to EC2 [16]) between the cracks in i direction (for intermediate bond of rebar), ε ucri , ε si = the values of the strains in the reinforcement calculated for the uncracked and fully cracked conditions (in case of unbonded reinforcement), respectively, ζ = a distribution coefficient allowing for tension stiffening (see Eq. 7.19 given in EC2 [16]). Therefore, the tension stiffening coefficient can be introduced to control the effective elasticity modulus of the reinforcing bars ( E si = E s ∕ si ) as follows: Substituting the coefficient ζ by Expression (7.19) (see EC2 [16]), after some transformations, we get the formula: where in turn σ si and σ sri are the steel stresses in cracked section, wherein σ sri concerns the first cracking. In order to take into account the mean stress in the tension reinforcement between the cracked sections, the steel stress in cracked section σ si can be replaced by the ratio of its mean value and the tension stiffening coefficient itself: σ smi /ψ si , as the stress σ smi is directly available from the iterative numerical analysis. The ratio of the strains ε ucri /ε si can be also determined in a numerical way and coefficient β can be assumed according to EC2 [16]. After some rearrangements, the coefficient ψ si takes the form:  where σ sri can be assume as equal to f ct (1 + α e ρ si )/ρ si , where α e = E s /E cm and ρ si = A si /A c . Since the tension stiffening effect has been taken into account in accordance with EC2 [16], so, in accordance with this Eurocode, the participation of concrete in tension in directions perpendicular to cracks was omitted. On the basis of Eqs. (19) and (23), the generalised coefficient of steel stiffness variation can be formulated as:

Crack propagation in concrete
The first step of the nonlinear procedure is to find a solution for the elasticity problem for a doubly reinforced concrete . slab, where coupling between membrane deformations and flexure is caused only by the lack of symmetry of the reinforcement layers in the slab. Ten iteratively defined cracking patterns, which can be identified step by step in case of bending and membrane deformations of the plate in elastoplastic phase of RC planar structures, shown in Fig. 6 were used [11] in this formulation.
The cracking patterns represent all types of cracking sets that may arise as a result of loads acting in the plane of the RC plate and as a result of its bending with twisting. The adoption of each pattern implies the division of slab thickness into several layers (a maximum of three layers of concrete). This enables the replacement of the summation of constituent stiffnesses over the layers (as in the case of the layered method) by the appropriate analytical or simple numerical integration through the total thickness. This approach, called ICPM above, allows for nonlinear analysis of relatively low numerical effort in comparison with layered analysis, and furthermore, the tips of the considered cracks need not lie on the interfaces of the previously imposed layers. The depth of a layer of cracked concrete is calculated at each step of iteration from the condition of zero principal tensile stress perpendicular to the crack. A maximum of three layers of concrete may Fig. 6 Cracking patterns of arbitrary RC slab. Hatched areas mean compression zones, while blank areas-cracked regions exist simultaneously: uncracked, with one-way cracks and with two-way cracks. In case of assumed constitutive model, there is no need to include a fourth layer for the crushed concrete. However, such an approach is possible and may be adopted when necessary.

Finite element solution procedure
The program for numerical calculations has been written in FreeFem++ language developed at the Université Pierre et Marie Curie ( [30]) and distributed by GPL licence. For the finite element solution, a plane triangular element with six nodes is used. Problem description in FreeFem++ requires formal variational form, which was formulated for moderate thickness Mindlin slab. The generalised stress-strain relationship takes the form: where B αβδγ , C αβδγ , D αβδγ and H αβ , are secant constitutive tensor components that depend on the displacement field. For the purposes of the FEM, the constitutive tensor is converted into the matrix form. All components of the secant constitutive matrix consist of two parts connected with the participation of concrete and steel reinforcement. The contribution of concrete in the secant constitutive matrix defined in the rotated nt set of coordinates has the form: in which the respective terms can be expressed as where β j is a shear-retention factor (taking into account the dowel actions and aggregate interlock in rough cracks), where in turn j = 1 for i = 2 and 3 (one-way cracks) and j = 2 for i = 4 (two-way cracks). Parameters of integration: Δz 1 ÷ Δz 8 are given in Table 1, where Δz l (l = n, t) are ranges of concrete layers with respect to the cracking pattern (Fig. 6) and can be determined numerically based on the condition of zero stress in a cracked cross section in the nt coordinate system depending on the strain and curvature from the previous iteration. The nt coordinate system is rotated relative to the xy system by an α angle in counterclockwise direction, and the direction n is normal to the first appearing crack. Similarly, the participation of steel reinforcement in the secant constitutive matrix defined in xy set of coordinates (for a given orthogonal reinforcement system) can be presented as follows in which the respective terms can be expressed as where where in turn where sq is the generalised coefficient of steel stiffness variation according to Eq. (24) and A sq is the cross-sectional area of reinforcement per unit width in the direction q, at a distance z sq from the midsurface of the slab, where in turn q = x, y. The index ( ′ ) in case of A sq, z sq and sq denotes the top reinforcement.
Let us investigate the ′ c component of constitutive matrix that is related to concrete behaviour (26) in a rotated (nt) set of coordinates and can be expressed as: Coefficients k i t and k i n can be expressed as: k n 3 = k n 4 = k t 2 = k t 4 = 0 and k n 1 = k n 2 = k t 1 = k t 3 = 5∕6. The contribution of steel reinforcement in the constitutive matrix component of a Mindlin slab is proportional to the ratios of reinforcement in particular directions. Taking into account the coefficients due to tension stiffening and plasticity, it is possible to write that where For the purpose of brittle punching phenomenon description, the internal bonds were used by means of the imposed locking of shear strains. These bonds were adopted in vertical plane stripes along the sets of the aforementioned reinforcement.
The constitutive matrix of RC structure consists of contributions of steel reinforcement as well as cracked or uncracked concrete where the matrix [D s ] is given by (30)

Supporting test programme for validation
The authors tested three transversely reinforced RC slabs (test specimens) in ITB Strength Tests Laboratory. The external dimensions of the plate of each test specimens were 2650 × 2650 mm, and the plate thickness 200 mm. The dimensions of a short column were: height 600 mm, 42.5R per 1 m 3 of dry mixture, 55 kg of additives (lime milk, flying ash), 639 kg of gravel 8/16, 522 kg of gravel 2/8, 649 kg of sand 0/2, 164 kg of water and 1.6 kg of FM 10 J admixture. The mechanical properties of the concrete, such as elasticity modulus (E cm ) and compressive strength (f c ) as well as tensile strength, were obtained from the tests results of cube and cylinder specimens. Six samples were tested for each strength parameter given in Table 2. The elasticity modulus (E cm ) and concrete tensile strength were evaluated only for the slab model № PI.
A comparatively low value of tensile strength of concrete (f ctm ) was achieved, which amounted of 2.13 MPa according to Brazilian split tests on concrete cylinders. Based on the results obtained on six samples, it was found that the minimum value of concrete tensile strength obtained during the tests was 1.78 MPa, while the maximum value of this strength was 2.48 MPa. The values of the reinforcement yield point (f y ) and tensile strength (f t ) were obtained from tensile tests results. Despite a rigid construction of the test stand, the displacements at specimen perimeter were measured during the tests. This allowed for the subsequent compensation of unpredictable deformations of the supporting structure. Square specimens simply supported around the perimeter were loaded at the centre by three 1000-kN hydraulic actuators (see Fig. 8). Horizontal and vertical views of test set-up configurations are presented in Fig. 8. The test in progress is visible in Fig. 9. For the purpose of the measurement of displacements, a number of dial gauges (numbered from 1 to 9) and the cathetometer (target point № 10) were used. The dial gauges layout is shown in Fig. 7.

Tests results
The obtained results, showing both the load bearing capacity of the tested specimens and values obtained from the material testing are listed in Table 2. Basically, two different modes of punching shear are encountered in case of RC slabs supported on columns: a typical mode for the slab without shear reinforcement (Fig. 10a) and for the slab specimen with shear reinforcement (Fig. 10b), considered herein. The real mode of the punching shear with shear reinforcement could be observed after the test and the cutting of the tested slab. The photograph of slab specimen № PI after the cut along the axis of symmetry is presented in Fig. 10b), and the location of inclined cracks caused by shear force is clearly visible. The cracking pattern on the upper surface of tested slab model № PIII is shown in Fig. 11. These photographs have been taken after the tests finish.
Comparison of test and analytical charts showing the relationship between force acting on the short column from the bottom and the slabs midpoint deflections is given in Fig. 12.

Analysis results and discussion
The primary aim of this paper is to present the computational model of RC slab-column connection and to compare the results with data obtained from the tests described herein.
The computational examples are presented, in which the stiff slab-column connection is used at the base of the floor slabs. Three numerical models have been analysed; the initial FE Fig. 9 Test in progress. The measurement of displacements by the cathetometer  Fig. 13.
Finite element meshes adopted for calculation of ultimate capacity (after a sensitivity analysis) for the models without and with transverse reinforcement are presented in Fig. 14. Midpoint deflections obtained from both experimental tests and nonlinear analysis of RC slab with and without transverse reinforcement are presented in Fig. 12. Finally, the fracture zones (dark colour), localised numerically according to the Podgórski failure criterion, are shown in Fig. 15. The crack propagation simulation on the top surface of the slabs shows that the real cracking mode presented in Fig. 11 is similar as the analysis result for the initial cracking around the column (Fig. 13a), while the final cracking modes (Fig. 13b, c) form the strips similar as the results of the yield line method for RC slabs. Comparison of the load bearing capacity of such modelled transversally reinforced and unreinforced test specimens indicates a correct response of the slabs in both cases. The post-peak parts of the load-deflection curves are not demonstrated in Fig. 12 as all the curves (both experimental and analytical) were obtained under the force (not displacement) control, in which authors' opinion is more realistic. It is visible, that the force-displacement relationship for numerical models of the slabs with and without transverse reinforcement is almost identical up to the 260 kN load (Fig. 12). At this point, the model without transverse reinforcement suffers an  abrupt loss of load carrying capacity (which corresponds to the brittle nature of the failure mechanism due to punching). A part of the load-deflection curve for the model with transversal reinforcement continues to ascend up to about 900 kN load, where the failure occurs, due to the ultimate flexural strength of the RC slab. The relation between the experimental data (ultimate force 980 kN-mean from the tests) and the numerical results shows a little underestimation of the slab load carrying capacity. In the first range of displacement values (from 3 to 16 mm), the discrepancies between experimental and numerical results are a consequence of the differences between the values of the existing tensile strength of concrete and the tensile strength assumed in computations, which leads to different stress redistributions after cracking. However, the slope character of the plots does not differ much as it results from the value of the elasticity modulus of concrete which is less susceptible to random variation in the properties of this material. The different shapes of the load-deflection curves in the final sections result from the computational application of the hypo-elastic steel model without hardening, while the hardening is observed in tests. However, due to simplicity of this model and being on the safe side, the authors of this paper have assumed that it is adequate for the computational purposes.
In the model proposed herein, additional internal bonds were used by means of the imposed locking of shear strains. These bonds are located in vertical plane stripes along the radial sets of transversal reinforcement (see Figs. 7 and 14). This modelling strategy of the double-headed studs stems from the very concept of their use consisting in such stiffening of the support zone that the failure occurs due to the bending, not punching shear (or both phenomena can occur simultaneously). Such assumption was confirmed by many tests performed at the German universities [10,18,19] as well as the calculations performed in accordance with the Eurocode 2 [16]. A series of slab specimens reinforced by the sets of double-headed studs was also tested at ITB Strength Testing Laboratory, and it turned out that the yielding of any stud due to stretching in the tests carried out as described above is practically impossible. It is visible that the results obtained for the model with transverse reinforcement are comparable to the real mode of punching shear shown in Fig. 10b), while without transverse reinforcement-to the real mode shown in Fig. 10a). Some disturbances can occur only in a small region around the column, because in the model of moderate thickness slab the vertical stresses and strains are neglected. Notwithstanding, in the assumed FE model with transverse reinforcement the computational destructive crack initiates close to the slab-column connection and runs through the studs up to the unreinforced region of the slab. The shape of the contour lines of the deflection surface at a failure for both slabs is presented on the same figure (Fig. 15). A significant increase in displacement values around the column perimeter, typical for punching failure mechanism, can be observed in the case of the model without transverse reinforcement. The failure mechanism of the RC slab depends on the punching to bending strength ratio. In case of a sufficiently high transversal  [17] for the models without and with transverse reinforcement reinforcement ratio, it is possible to observe ductile failure mechanism due to bending. However, the slabs with low transversal reinforcement ratio or unreinforced transversally suffer abrupt rupture due to brittle punching prior to reaching the ultimate load value corresponding to bending. Furthermore, the force-displacement relationship does not depend on the level of transversal reinforcement up to the point of bifurcation towards the brittle or ductile failure.

Summary and conclusions
Implementation of the described above-combined ICPM/ FSLM method in the nonlinear FE analysis of Mindlin-type slabs enabled us to devise the computer program for the analysis of moderate thickness RC slabs subjected to punching and the recognition of failure mechanism, with respect to the layout of transversal reinforcement. The computational results were confronted with the test results. Based on the results of such comparative analysis, the following conclusions are drawn: 1. The prediction of shear capacity was achieved by the above-described FSLM method. As a result, both types of damage were modelled separately but also, the combined flexure-shear mode failure can be considered by using the combined ICPM/FSLM method. This approach aided by the imaging technique allowed for the distinction between both phenomena of the failure. The numerical analysis results were compared with experimental results of the three slab specimens (reinforced by the sets of studs with connecting flat at the top) tested at ITB Strength Testing Laboratory obtaining fairly good accordance with the experimental data. 2. As the analytical method proposed in this paper has no direct reference to the methods described in the cited literature, the comparative analysis is limited in nature. However, it should be pointed out that a considerable advantage of the proposed approach is a relatively low numerical effort in comparison with layered or fully three-dimensional models because of the plane finite element discretisation of the floor slab and due to the division of the slab into a small number of concrete layers (usually two layers in cracked regions), which is additionally confirmed by a short computation time. Moreover, finite element plane discretisation is more suitable for use in the practical design of floor slabs. 3. In scope of surface crack modelling and load-deflection curves, a comparable accuracy was obtained as in the paper [3] and better than in the other cited publications. On the other hand, in the case of cited papers in which moderate thickness slab elements were used, none of them presented the computationally determined punching surface in a vertical cross section, taking into account the influence of double-head studs. From the point of view of the modelling accuracy of the loaddeflection curves, the presented method also gives more accurate results than in the cited papers in which threedimensional elements were used. However, the computational effort associated with the numerical implementation of the proposed method was much smaller compared to other available models using the threedimensional or layered FE analysis. 4. The novel formulas were used for nonlinear constitutive relationships for concrete and for tension stiffening effect. The constitutive model of concrete proposed by Kotsovos and Pavlović [13] was modified according to the limit state condition of Podgórski [17] to fit the material model to the experimental data valid for the considered problem. The model was subjected to calibration procedures for uniaxial and biaxial compression and was compared with the test data of Kupfer et al. [15]. Additionally, in order to avoid the overestimation of slab capacity in case of brittle punching, a special descending path was incorporated in the constitutive model for the range beyond the peak stress. 5. The tension stiffening effect is taken into account, with the assumption that additional stress is carried by the reinforcement. The adequate model was derived according to the provisions of EC 2 [16]; however, it was mathematically converted to take into account the mean stress in reinforcement between the cracks, without loss of physical meaning. 6. Therefore, the overall prediction of the nonlinear behaviour and the bearing capacity of the test slabs were improved in comparison with the former results [12] and show a good accordance with experimental data. The relation between experimental data (average ultimate force 980 kN) and numerical results (ultimate force 900 kN) shows little underestimation of the slab load carrying capacity. The reason is that plastic hardening was not taken into account in the implemented steel model (see Fig. 5). However, it was acceptable due to its simplicity and in order to stay on the safe side.