Punching shear strength of eccentrically loaded RC flat slabs without transverse reinforcement

The effects of moment transfer at the flat-slab-inner-column connection, and the influence of concrete strength on the punching resistance of slabs, are analysed in the paper. Seven full scale specimens are investigated experimentally and numerically. The results are presented as load-rotation curves and compared with the critical shear crack theory (CSCT) failure criterion, and with other comparable expressions. A detailed nonlinear FEM analysis, involving calibration of model with experimental data, is conducted to ensure better correspondence of numerically obtained loadrotation relationship with experimental results.


Introduction
Reinforced concrete (RC) flat slabs resting on columns are among the most commonly used structural systems for multistorey buildings. Lacking the moment frames, these systems are expected to be adequately stiffened in horizontal direction by RC walls, providing the slab-column connection does not contribute significantly to lateral strength; i.e. moments transferring from column to the slab are rather small. However, their presence is always unfavourable, and may significantly affect the punching strength of the slab. They could be induced by horizontal (seismic) lateral drift of structural system, by uneven distribution of gravity load, or due to varying spans of continuous slabs. Over the past two decades, the punching shear design of flat slabs has been mostly based on the critical shear approach. Relying on the research of Kinnunen and Nylander conducted in the 1960s [1], in which the punching resistance of slabs was related to development of a critical crack, Muttoni et al. formulated a mechanical model for assessing the punching strength of flat slabs [2][3][4], which is known as the critical shear crack theory (CSCT). According to this proposal, in addition to the concrete compressive strength, f c , the shear resistance of a slab is also greatly affected by the width of a critical crack, and decreases with crack development. The rational explanation of the theory is based on the idea that a critical shear crack ( Figure 1) propagates along a compressed inclined strut, which transfers the shear to a column, reducing its shear strength. The width of the critical crack, w, is assumed to be proportional to the slab radial rotation, ψ (Figure 1, eqn. (1)).

Figure 1. Critical shear crack model
In addition to crack width, the maximum aggregate size also exerts -via the interlocking effect (crack roughness) -a considerable influence on the level of shear to be transferred via cracks. According to Walraven [5] and Vecchio and Collins [6], crack roughness and its capacity to carry shear forces may be accounted for by dividing crack width by the sum (d g0 +d g ), where d g is the maximum aggregate size, while d g0 is the reference size. In this way, the punching resistance of the slab becomes a function of the following factor (d is the effective slab depth): , w = ψ ·d, d g0 = 16 mm (1) On that basis, Muttoni [3] proposed (for the first time in 2003) a hyperbolic formulation for the failure criterion, which shows a rather good correspondence with experimental results (although the dispersion of results around the failure curve is significant). The criterion is written in two forms: the first one targets a mean value ("mean" criterion) of the analysed experimentally obtained punching strengths, while the second one [7], the "safe" criterion, targets a 5 %-fractile value, tending to include various sources of irregularities: where b 0 is the perimeter of the critical section located at a distance d/2 from the column edge ( Figure 1). Such an approach has been adopted in Swiss code [8], fib pre-normative Model Code 2010 [7] and in the "new" Eurocode 2 [9]. The above criteria have been derived for the centric case, in which the distribution of shear forces along the perimeter is close to uniform. However, the punching shear resistance is reduced because of non-uniform shear distributions, which may occur due to concentrations at the corners of large loaded areas, slab discontinuities, or, which is the case in this investigation, moment transfer between the slab and the column. These effects may be accounted for by reduction of control perimeter: generally, the shear-resisting control perimeter is related to the maximum shear force per unit length perpendicular to the perimeter. To take into account a variable shear distribution due to moment transfer, the control perimeter may be assumed reduced in the following way [7]: where e u is the eccentricity of the resultant of shear forces with respect to the centroid of the basic perimeter, and b u is the diameter of a circle with the same area as the region inside the basic parameter. For a square column section having an edge length b c , the diameter of an equivalent circle is shown in Figure 2. Punching shear strength of eccentrically loaded RC flat slabs without transverse reinforcement Knowing the failure criterion, the punching shear strength and the related deformation capacity are determined by the point of intersection between the failure criterion and the load-rotation curve, which mainly represents flexural behaviour of the slab, i.e. an increase in rotation with an increase in load ( Figure 3).

Figure 3. Punching shear strength and rotation capacity
The load-rotation relationship may be determined in a number of ways. Kinnunen and Nylander developed an analytical formulation for an axis-symmetric case [1] using the bilinear moment-curvature relationship. Muttoni [3,4] improved it by accounting for the tension stiffening effect, using the quadrilinear relationship. For practical purposes, a load-rotation curve may be approximated by the parabolic function having an exponent of 3/2 (like in Level II or Level III approximation in fib Model Code 2010, [7,10]): where r s denotes the position in which the radial bending moment is zero with respect to the support axis, f y is the yield strength of the reinforcement, and E s is the reinforcement modulus of elasticity. The moment m is an average moment per unit length in the support strip, while m R is an average flexural strength per unit length in the support strip. The width of the support strip, b s , is determined as follows ( Figure 4): In eqn. (5), e u is the resultant eccentricity of shear forces with respect to the centroid of the basic control perimeter, V/8 is the average moment acting in support strip without moment transfer, while the moment transferred to a column (V·e u ) acts across the width of the support strip, b s , half on each side. The coefficient of proportionality k m equals 1.5 for Level II of approximation, while it may be replaced with 1.2 if a more refined analysis is taken for designing flexural reinforcement (Level III approximation according to fib Model Code 2010). However, it should be noted that expression (5) is a parabolic approximation of an analytical solution where the tensile strength of the concrete and tension stiffening effects are neglected. As a consequence, a relatively large distance between the experimental and the approximate solution may be expected for lower rotations and for concrete of higher tensile strength ( Figure 5). Finally, the relationship may be obtained through a nonlinear analysis of the structure, taking into account the cracking, tension-stiffening effects, yielding of reinforcement or any other nonlinear effect relevant for good assessment of the system (Level IV approximation in fib Model Code 2010). The research presented in this paper is related to the punching resistance of the RC flat slabs of various concrete strengths subjected to the combined vertical load and bending moment of variable intensity. It is limited to slabs without transverse reinforcement. Such flat slabs constitute a large portion of existing and newly constructed ones, and they may be strengthened against punching. Furthermore, the exclusion of transverse reinforcement as a relevant factor provides conditions for separate treatment of analysed phenomena. The experimental program involving seven full-scale specimens differing in concrete strength and transferred moment intensity was conducted. The unbalanced moment was introduced by GRAĐEVINAR 70 (2018) 9, 757-770 Zoran Brujić, Danijel Kukaras, Radomir Folić, Sohela Ali, Arpad Čeh an eccentric vertical load at constant eccentricity (with load change), simulating the moments induced by unequal spans or an uneven load distribution within real structures. The acquired data were reviewed in the light of the CSCT failure criteria (eqns. (2) and (3)) and approximate load-rotation proposals, eqn. (5). A more refined nonlinear FEM analysis of the tested specimens is carried out using ANSYS Mechanical 14.5 software. The restriction involving the constant shear modulus, which is set by the software, has been overcome through the calibration process in which a numeric model is calibrated against test results taken as reference values.

Experimental analysis of eccentrically loaded flat slab specimens
The experimental research was performed by the authors in 2015 and 2016 at the Laboratory of the Faculty of Civil Engineering in Subotica. The aim of the experimental research was to determine the dependency of the punching shear strength of the slab with respect to, on the one hand, the bending moment that is transferred from column to slab and, on the other hand, the compressive strength of concrete.
In total, seven specimens of the same geometry (marked S1 to S7) were prepared and tested to specimen failure. A single interior column and the surrounding part of the slab were isolated as a full-scale specimen, and the load was applied in the form of an eccentric vertical force, thus enabling a combined transfer of force and moment. Such setup, where the eccentricity remains constant, is suitable for the analysis of transferring moments induced primarily by uneven spans or load distribution. For the moments induced by horizontal forces, the deformation of the slab outside the modelled part of the slab is significant and must also be analysed [15]. The shape and dimensions of the specimens ( Figure 6) were chosen so as to roughly correspond to the part of the flat slab structure of a span common for multi-storey structures (approximately 4 m) in the vicinity of the column.
The specimen represents the part of the slab (1) inside a circularly shaped section of contraflexure of radial moments. For elastic slabs of constant stiffness, such a section was approximately located at a distance 0.22L from the column axis, where L is the slab span. The slab depth was 180 mm and it was cast in an octagonal shape providing for a simpler formwork. The column of the specimen (2) was of square cross-section, with dimensions 250x250 mm. At the bottom, it ended with the horizontal extension in the form of a corbel (3), which enabled eccentric load application. To avoid unwanted local effects of reinforcement anchoring, the column reinforcement was anchored straight into the concrete block (4) cast over the top of the slab. The reinforcement of all specimens was identical. The slab was reinforced with Ø14/100 mm rebars placed in two orthogonal directions of the upper slab region, and with Ø10/200 mm rebars similarly placed in the lower slab region (see Figure  21). The column was strengthened with 14Ø16 bars that were evenly distributed across its circumference. The overall quantity of the main slab reinforcement was chosen so as to ensure a "moderate" reinforcement ratio of approximately 1 %. Specimens differed only with respect to concrete class: first three specimens were made of the concrete mixture for concrete class C30/37, while the mixtures for the remaining specimens were designed for concrete classes ranging from C60/75 to C80/95. Five different concrete mixtures were used in total, while the maximum aggregate size, d g , was 16 mm. The specimen dimensions and reinforcement were selected so as to ensure that the punching failure occurs prior to flexural failure. The testing was conducted in such a way ( Figure 8, Figure 9) that a specimen (1) with its corbel rested on hydraulic jacks, which were placed within the adjustable box frame (3). The frame was fixed on top of the star shaped steel profiles (2) where two steel profiles formed one spoke. The total of eight high-strength bolts (4) were fixed at the tip of each spoke and threaded through the corresponding openings on the slab, where they were fixed at the upper slab surface with a pair of 8mm thick steel plates (5). The load eccentricity was ensured by accurate positioning of the adjustable box frame at the specified location with respect to vertical axes of the column. The initial position of the specimen was achieved by manually tightening high-strength bolt nuts onto the upper surface of the slab until the perfect horizontal positioning.

Figure 8. Experiment setup: 1 -concrete specimen, 2 -steel support, 3movable hydraulic jack box, 4 -threaded rods, 5 -steel plates
Simultaneously with the casting of the concrete, testing specimens, cubes and cylinders, were performed in order to determine real mechanical properties of concrete at the moment of punching shear testing. The age of specimens used in the testing ranged between 137 and 188 days. Age information for individual samples is given in Table 1, together with the measured concrete compressive strength at the moment of testing, and the corresponding characteristic compressive cylinder-based strength, calculated according to the Eurocode strength increase function for the Class N cement. A reasonably good correspondence with the projected values was achieved. The behaviour of each specimen during testing was monitored using an appropriate measuring equipment. Displacement sensors were used to register vertical displacements of the upper slab surface (Figure 10.a), as well as its horizontal movement. A total of eleven strain gauges for measuring dilation of steel were installed at appropriate locations on top reinforcement bars ( Figure 10.b): five gauges (S1 to S5) along the middle x-direction bar, three gauges along the side x-direction bar, and three along y-direction middle bar. The same strain gauges were installed on each high-strength bolt (out of the eight bolts in total) which, together with additional two (used for control measurements), makes a total of ten gauges. Strain gauges for concrete were installed on the lower surface of the slab, near the column. Measurement data acquisition was conducted with the 72 channel data logger ( Figure 11) at the sampling rate of 5 Hz. The loading protocol was similar for all specimens. The loading protocol for specimen S4 is shown with the stepped line in Figure 12. Very similar procedures were used for other samples: the load was gradually increased up to approximately 40 % of the projected punching resistance. Then the specimen was almost completely unloaded and, finally, the load was increased monotonically, in small increments followed by pauses of approximately 3 minutes, until failure. The intensity of the force was determined in an approximate manner. However, the strain gauges set on threaded bolts provided, after calibration, forces values in bolts, acting as slab supports. On the basis of these forces, the exact values of the applied force and its eccentricity were calculated and plotted in Figure 12. The eccentricity of the applied force remains practically constant, while eccentricity in y-direction is negligible. As an illustration, deformations measured in top reinforcement of Specimen S2, at four selected positions, are plotted in Figure 13. Calculated eccentricities of applied force, and force values causing punching failure, are presented in Table 2. The eccentricities determined along y-axis were smaller than 5mm. All seven specimens failed due to brittle punching, after development of plastic deformations in top reinforcement of the slab.  Punching shear strength of eccentrically loaded RC flat slabs without transverse reinforcement which was obtained after removing the branches -that correspond to the unloading and reloading (to the previously reached value) -from the set of measured data. Although this resulted in lower plot accuracy in initial parts of the curves, it is the upper part that is of primary interest in the analyses. Having in mind that, due to eccentric load, the deformed shape of the specimen is asymmetric, and that the tangential change of the radial rotation can be approximated with the sine function (as indicated for instance in [11], and confirmed by experimental data), the reference slab rotations were determined as the arithmetic mean of rotations in the direction of the eccentricity (rotation ψ 0 ) and in the opposite direction (rotation ψ 180 ), according to the rules shown in Figure 14: in order to eliminate displacements due to support deformation, rotations were determined on the basis of relative displacements -difference in measured displacements (u i ) at points D2-D1, and D4-D5: ,  The change of experimentally determined slab rotations with an increase in load is shown in Figure 15 for all seven specimens. In order to put the results in the form needed for further analysis, the curves were plotted as monotonically increasing functions,   As expected, due to higher stress values, the curves that correspond to a larger eccentricities were less steep when compared to those that correspond to smaller ones.

Experimental results vs. design proposal
The experimentally obtained load-rotation relationships are shown in the diagram in Figure 16 as dashed curves, while their ends, the values which correspond to the specimen punching failure, are marked with circles and compared to two CSCT failure criteria, presented in eqns (2) and (3). The chart is plotted in a relative form in which failure curves are not dependent on eccentricity or concrete strength, i.e. the same failure curves are valid for all specimens. This is achieved by dividing the load-rotation relationship by the coefficient of eccentricity, k e (eqn. (4)). The resulting punching forces are in the domain of expected values defined by the "mean" failure criterion, eqn. (2), except for specimen S3 where failure occurred earlier. However, even in this case, the punching failure strength corresponds to the "safe" criterion, eqn. (3). Besides, it was observed that specimens with higher concrete strengths, specimen S7 in particular, exhibit a somewhat higher punching strength, i.e. up to 10 % higher compared to that estimated using the "mean" criterion.
In the following analysis, the experimentally obtained loadrotation relationships are compared to the approximate expression corresponding to the approximation level III of fib Model Code 2010 (eqn. (5); factor of proportionality, k m , equals 1.2). It has already been stated that such an approximation primarily targets the rotational capacity of the slab, and so significant deviations of approximate load-rotation curves from the test data are expected. A brief analysis of the expression, eqn. (5), shows the relationship is determined by concrete properties only to a small extent. Namely, concrete strength affects only the average flexural strength per unit length, m R , and this effect is almost negligible, as shown in Table 3 where the bending resistance is calculated according to the rigid-plastic theory, as in [10], but put in a non-design form: The effective height, d, is 146 mm (Figure 21), the yield strength of flexural reinforcement, f y , equals 500 MPa, and the flexural reinforcement ratio, ρ, is 1.055 %.
The results are shown in the following charts ( Figure 17 and Figure 18) where significant differences between the compared curves may be observed. They increase with load eccentricity and especially with concrete strength (tensile strength in particular). Good correspondence between compared curves is achieved only for sample S1, where the values of concrete strength and load eccentricity are relatively low. Due to the nature of approximation involved, where the contribution of tensile strength of concrete and tension-stiffening effects are neglected, all deviations are on the side of safety. Furthermore, the chosen specimen geometry is characterized with a relatively low column-size-to-slab-depth ratio which, as concluded in [11], leads towards lower rotation capacities of slabs, thus additionally emphasizing the conservative nature of the approximation. Consequently, the results of the applied approximation appear as overly conservative for practical purposes. Thus, in the analysed domain of flat slabs, considering moderate reinforcement ratios, higher concrete strengths and relatively large moments transferring at columnslab connection, it is advisable to use either more complex analytic predictions where the effects related to the concrete tensile strength are involved (for instance, proposals based on the quadrilinear moment-curvature relationship

Numerical modelling of specimens
The numerical modelling of experimental work is performed using the ANSYS software (v14.5). Concrete is modelled with the SOLID65 element, defined by eight nodes and having three degrees of freedom at each node. It is commonly used for the 3D modelling of concrete with or without reinforcement. Its most important aspect is treatment of nonlinear material properties, i.e. it is capable of modelling: cracking in tension, crushing in compression, plastic deformation, and creep. Although up to three different smeared reinforcement (rebar) specifications may be defined in the element, in order to achieve a more realistic numerical simulation of the experimental program no rebar specifications are used, because punching failure is largely affected by the shear resistance of the system and the rebar element is not capable of shear. Instead, a discrete line element is used for reinforcement. The possibilities of modelling the shear-stress-strain relationship, or shear modulus change, for cracked concrete in ANSYS are limited to the constant shear modulus approach: The shear transfer coefficients β t and β c are introduced as the shear strength reduction (shear retention) factors for those subsequent loads that include sliding across the crack face. The coefficient β t refers to an open crack, while the coefficient β c is used for a closed crack. Both factors range between 0 and 1. Generally, the use of a shear retention factor is not consistent with the available experimental results in that the shear-stressstrain relationship is nonlinear and it is dependent on the crack width or the principal tensile strain [13]. Furthermore, there is no rational basis for selecting the value of the retention factor. Minding this, and being limited by the software capabilities, a variable shear retention factor approach is introduced in this work. The nonlinear concrete model, assigned to all SOLID65 elements, comprises a linear-elastic model, a multilinear isotropic hardening model (MISO) defined by uniaxial stressstrain relationship, and a CONCRETE material model [12,13] with cracking and crushing capabilities. If a material fails in uniaxial, biaxial, or tri-axial compression at any integration point, the material is assumed to crush at that point (complete loss of structural integrity). The quasi-plastic behaviour of concrete is introduced through a rate-independent plasticity. The multilinear isotropic hardening model (MISO) is chosen despite the fact that only the monotonic load change is analysed. As this model can not simulate crack formation, it is paired with the CONCRETE model that is responsible for controlling the cracking and crushing of concrete only. The cracking is modelled by a plane of weakness in the direction perpendicular to the crack face, through an adjustment of material properties, which effectively treats the cracking as a smeared band of cracks. Although introduced (Figure 19b) , the corresponding cylinderbased strengths at 28 days, and at the age corresponding to the time of testing, are determined in accordance with the Eurocode 2 concrete aging function and using an appropriate conversion ratio. The Eurocode 2 stress-strain relationship for non-linear structural analysis is used to determine an increasing part of the stress-strain relationship, which is presented as multilinear relationship defined by six points in the stress-strain plane (Figure 19a). In order to avoid possible numerical problems, the corbel element of the specimen, and the top parts of the slab that are in contact with the supporting steel plates, are modelled with the SOLID185 elements, to act as linearly isotropic. The linearelastic properties assigned to these elements correspond to those used for nonlinear concrete model ( Table 4). The reinforcement is modelled in a discrete form, through the use of a 2-node 3D BEAM188 element based on the Timoshenko beam theory that includes shear-deformation effects. The element has six degrees of freedom at each node, while the provided stress-stiffness terms enable it to analyse flexural, lateral and torsional problems in nonlinear manner. All reinforcement sections are defined as circular solid ones, having the diameter matching that of the reinforcement applied. The nonlinear steel model comprises the linear-elastic and bilinear isotropic hardening model (BISO). Linear-elastic model properties involve the modulus of elasticity, E, equalling to 200.000 N/mm 2 , and the Poisson's coefficient, ν, amounting to 0.3. The BISO model is utilized to describe nonlinear behaviour of reinforcing steel. It is defined by a bilinear effective stress versus effective strain curve. The initial branch of the curve is elastic; its slope corresponds to the elastic modulus of the material. The plastic strain develops beyond the specified initial yield stress, f y , and the second branch of the and stress vs. totalstrain curve is inclined at an angle defined by the userspecified tangent modulus, E T . The von Mises yield surface is used, the flowing rule is associative, perpendicular to the surface, while the hardening rule is determined by the tangent modulus E T . The yield stress is assumed to be 500 N/ mm 2 , while the tangent modulus of the post-elastic branch is 10000 N/mm 2 . Supporting steel plates are modelled using  Punching shear strength of eccentrically loaded RC flat slabs without transverse reinforcement SOLID185 elements that are defined with eight nodes, and with three degrees of freedom in each node: translations in x, y and z directions. It allows for prism, tetrahedral and pyramid degenerations when used in irregular regions. Since no plastic behaviour of steel plates is expected, only linear elastic material properties are assigned to the elements. The symmetry of the specimen about X-Y plane is used in order to reduce the model size (DOF). Only a half of the specimens is treated numerically (Figure 20), while the remaining half of the specimens is replaced with the corresponding boundary conditions (restrained Z translations in nodes belonging to the plane of symmetry). The concrete volume is modelled as a 3D solid volume that follows the exact dimensions of the specimens. In order to provide the collinear nodes of the mesh along the discrete bars of reinforcement (Figure 21), the solid model is divided into smaller volumes. The concrete volume is sliced whenever required by reinforcement position, while meshing of volumes is performed automatically, setting the desired length of a finite element edge to 25 mm. As a result, a total of 21838 solid elements (25389 nodes) are created (Figure 22a). The clear majority of created elements are 8-node parallelepipeds, while only a small number of elements (26 in total) in the zone of the supporting plates are wedge shaped. The longitudinal reinforcement of the slab and columns is torsional and is restrained at its ends by perpendicular elements connecting the reinforcement end nodes to the adjacent solid mesh nodes (Figure 22.b).
Although in experiments the specimens were loaded by applying the force from hydraulic jacks, in numerical simulations the load was applied in the form od displacement.
There is no significant difference in the final response of the model between two load application approaches, although a better convergence is achieved by nonlinear simulation. Thus, displacements, safely beyond the displacement referring to punching failure, are applied to nodes along the line perpendicular to the bending plane, positioned at the eccentricity which corresponds to the one measured during the experiment, at the bottom of the specimen. However, the measured eccentricities are variable (Table 2) and do not exactly match the solid elements mesh, and so intermediate nodes had to be provided through use of rigid beam elements spreading from solid element mesh nodes to appropriate intermediate nodes. Besides the applied displacements, a fully restrained displacement is provided in Z-direction for all nodes in the symmetry plane, in Y-direction for four nodes at the top of the supporting steel plates, and in X-direction for the nodes in close proximity of the loaded nodes ( Figure 23).

Figure 23. Boundary conditions
In this study, the Newton-Raphson procedure is applied in nonlinear finite element analyses to solve higher order differential equations. The accuracy and convergence of GRAĐEVINAR 70 (2018) 9, 757-770 Zoran Brujić, Danijel Kukaras, Radomir Folić, Sohela Ali, Arpad Čeh the solution depends on the factors such as mesh density, constitutive properties of concrete, convergence criteria, tolerance values, etc. In the problems for which a convergent solution might be hard to achieve, such as in the nonlinear reinforced concrete modelling, the use of the full Newton-Raphson solution procedure is suggested, paired with the load applied in very small increments.

Analysis of numerically obtained results
As an illustration (just for specimen S2), some of the obtained results are presented in a graphic form: the distribution of total mechanical strains shown in Figure 24a mostly implies the crack distribution layout, while Figure 24b shows the yielding parts of the tensioned top reinforcement, which roughly corresponds to the experimentally measured data shown in Figure 13. The distribution of deflections ( Figure 25) of the top surface along the X-direction linear path provides the possibility for determining characteristic and, consequently, reference (average) rotations in the same manner as if they were obtained experimentally ( Figure 14, eqn. (7)). In such a case, the loadrotation curve may be plotted for each analysed model and compared with the experimental one, as shown in Figure 15. However, as mentioned earlier (see Table 4), each specimen is modelled using several different values for the open crack shear transfer coefficient (hereinafter referred to as the shear retention factor or SRF). As an illustration, Figure 26 shows numerically obtained load-rotation curves, and an experimentally determined curve, for four different values of this coefficient, related to specimen 4 (chosen as illustrative; there is no significant difference when compared with other specimens). Despite the shortfall registered for small rotation values (which is probably due to underestimation of tensile strength of concrete in numerical model), the comparison of experimental and numerically obtained curves clearly shows that the chosen value of the shear retention factor may significantly determine the load-rotation curve, and that lower values of this factor are, consequently, related to a lower bending stiffness behaviour (see slope of plotted curves) of the modelled sample. It can be seen in experimental curve that the gradual decrease in stiffness is faster compared to the decrease registered in any numerical curve. Or, with an increase of rotations, the slope of the experimental load-rotation curve corresponds to the slope of the numerical curve with progressively lower SRF values. The above information implies the possibility of introducing a "variable" shear retention factor to control the change in model stiffness. Namely, if the change of a shear retention factor is known, the slope of the resulting numerical curve may be determined through interpolation between known slopes of the constant SRF curves. Consequently, the resulting load-rotation curve will match the experimental one much more closely than any of the constant SRF curves. After a detailed observation of the obtained curves and comparison with experimental ones, it can be concluded that the rate of SRF change is significantly influenced by concrete strength. Thus, the change in shear retention factor with an increase in rotations is analysed as monotonically and nonlinearly decreasing and dependent on the concrete compressive strength. For this particular analysis, it is assumed that the SRF value reaches zero at the rotation of 0.02, while nonlinear change is defined with the exponent α, as a function of concrete strength. Furthermore, through an extensive regression analysis, which involved "straightened" curves for all seven specimens (using moving average approximations), and which favoured better matching in "fully cracked" zones of the diagrams (being of greater interest here), the following simple form of dependence is obtained: Knowing the SRF value for a rotation, the slope of the curve is obtained using parabolic interpolation, and then numerical integration is used to determine the resulting numerical loadrotation curve. In the following charts (in order to avoid an overlapping of curves, the data are divided into two charts), the resulting numerical simulation curves are plotted against the experimentally obtained ones, cf. Figure 27 and Figure 28. A reasonably good correspondence between numerical curves and experimental results has been achieved. Finally, based on numerically obtained load-rotation relationships and two failure criteria ("mean" criterion defined by expression (2) and "safe" criterion defined by expression (3)), two values of punching shear strengths per specimen are determined, as intersection points between the numerical load-rotation and two failure criteria: V FC -the "mean" failure criterion strength, and V FCS -the "safe" failure criterion strength. For the selected set of specimens, these values are plotted in Figure 29 and compared to the corresponding experimental values. In Table 5, these values are listed for all seven specimens and compared in relative form to experimental values. With the exception of specimen S3, a very good correspondence has been achieved, with deviations of less than 10 %, mostly on the safe side. Even the results for specimen S3 fit the "safe" criterion.

Conclusions
The parameters analysed in this paper are: load eccentricity and concrete strength. Experimental specimens are related to flat slabs of constant geometry (column size, slab span, and slab depth), vertical load at fixed eccentricity, and constant reinforcement ratio and layout, which limits the scope of this analysis. The comparison between experimentally obtained punching shear resistance values and the predicted values determined by "mean" failure criterion, eqn. (2), proposed by Muttoni et al., which includes the effects of transferring moments through reduction of the control perimeter, eqn. (4), revealed a relatively good correspondence of values, which confirms the applicability of this proposal. Although the criterion slightly underestimates the punching resistance of samples which utilize higher concrete strengths, the registered differences are within the acceptable margin (difference of up to 10 %). However, the punching phenomenon is characterized by a wide dispersion of resistance values, which conservatively leads towards the use of "safe" failure criterion for design purposes, eqn. (3). Such an approach proved valid even in actually performed investigations, where a punching shear significantly lower than expected was obtained for specimen S3, but it still satisfied the reduced criterion. The approximation of load-rotation relationships using the parabolic function, eqn. (5), showed significant safe-side differences, primarily due to neglect of tensile effects. In the analysed group/domain of structures, which are characterized by relatively low rotation capacities, such simplification leads to results that are overly conservative for practical purposes.
If an appropriate FEM analysis is used instead of the above procedure, satisfactory load-rotation relationships can be obtained, but only after software limitations present in the approach based on constant shear modulus are overcome. By adopting experimental results as reference, after regression analysis, this has been achieved by establishing the relationship between shear retention and rotations of the slab. Such an approach opens possibilities for further numeric experiments, which would extend the analysed domain by including structures characterized by intermediate values of the varied parameters. Due to high conservativeness of approximate proposals and a small number of experimentally considered specimens, the applicability of the suggested method of accounting for the effects of load eccentricity, eqn. (5), cannot be perceived in a quantitative manner. Still, the deviations shown in Figure 17 indicate that the effects of load eccentricity are overestimated, again on the safe side.