Quantification of the resistance modeling uncertainty of 19 alternative 2D nonlinear finite element approaches benchmarked against 101 experiments on reinforced concrete beams

Nineteen 2D nonlinear finite element analysis (NLFEA) solution strategies were benchmarked against a wide variety of 101 experiments on reinforced concrete beams failing in bending, flexural shear, or shear compression. The relatively high number of solution strategies was motivated by the conviction that choices for the constitutive models, the finite element kinematics and equilibrium settings will interact, and must therefore be tested in conjunction. Modeling uncertainty distribution parameters are presented for the 19 solution strategies, using all beams, and using beams with and without stirrups separately. The resulting statistics are discussed against the correctness of the simulated failure modes and failure loads, revealing that rotating crack models perform well for the relatively ductile failures in beams with stirrups, while fixed crack models perform better for the more brittle failures in beams without stirrups. The tailored solution strategies that predict failure modes correctly, imply a log‐normal distribution of the modeling uncertainty with relatively low coefficients of variation. The outlook is that these estimates of the statistical properties of the modeling uncertainties could serve as a basis within safety formats.


| INTRODUCTION
Two decades ago, Frank Vecchio published a milestone paper in this journal with the inspiring title "Nonlinear finite element analysis of reinforced concrete: at the Discussion on this paper must be submitted within two months of the print publication. The discussion will then be published in print, along with the authors' closure, if any, approximately nine months after the print publication. crossroads?" 1 : Applications of nonlinear finite element analysis (NLFEA) should be done with a healthy degree of caution and skepticism. Wherever possible, analysis models should be validated or calibrated against benchmark tests involving specimens of similar construction and loading details, dependent on mechanisms anticipated to be significant of the analysis problem at hand. The paper further included criticism towards the research community, and associated technical committees, that may have failed the profession in some respects: many working in the area have directed their efforts to developing sophisticated numerical models and methods of analysis, not directly suited to reinforced concrete.
We now observe a clear trend break. Standards organizations have recognized the potentials of NLFEA-but are also well aware that NLFEA is much more difficult to codify than resistance models based on lower levels of approximation. 2 Being aware of these concerns, standards committees involved in codifying NLFEA focus on performance-based codes in which the accuracy of the NLFEAs for a specific application field should be evaluated before being used for structural design or assessment. This accuracy is typically expressed by a modeling uncertainty, which acknowledges that no model can be perfect. The modeling uncertainty should be determined by systematically confronting a welldefined NLFEA approach with a series of experimental benchmark studies, all relevant for a well-defined specific application field. The upcoming fib Model Code 2020 3 will include a completely revised section on NLFEA which will exactly follow this line of thought.
Schlune et al. 4 reviewed round robin NLFEA modeling competitions, concerning structural elements which failed in compression, in bending (with under-and overreinforced sections) or in shear. Typically, the results, blindly predicted by the round robin participants, were distributed around the response of the laboratory tests with a coefficient of variation ranging from 5% (for ductile flexural failure of under-reinforced cross-sections) to a staggering 40% (for brittle shear failure, either in diagonal tension or shear compression). Engen et al. 5 termed this notion of modeling uncertainty as "between-model uncertainty" as it provides a relevant measure of the modeling uncertainty if an analyst would randomly pick a reasonable NLFEA approach for the problem at hand, without prior verification of this approach with similar benchmarks. Results from renowned competitions have been documented by Collins et al., 6 van Mier & Ulfkjoer, 7 Jaeger & Marti, 8 and Collins et al. 9 and more recently by Yang et al. 10 Indeed, usually these round robin competitions were based on challenging tests, for which little comparable tests were available. Consequently, if a design code would assume a default value for a modeling uncertainty, thus assuming that no proper benchmarking is available, applying the notion of "between-model uncertainty" makes sense.
Other researchers followed the above line of thought in which a well-defined NLFEA approach, often termed as a "solution strategy", is tested against series of experimental benchmark studies, all within a certain application field. Engen et al. 5 termed this notion of modeling uncertainty as "within-model uncertainty" and it is this notion that makes sense to be used as a measure of the accuracy for wellcalibrated solution strategies. Allaix et al., 11 Engen et al., 12 and Cervenka et al. 13 all tested specific solution strategies, against, respectively, 14, 38, and 33 benchmarks, and thus quantified the within-model uncertainties for the respective NLFEA approaches. Belletti et al., 14,15 Castaldo et al., 16,17 and Gino et al. 18 tested in their studies up to 18 approaches, thus resulting in multiple within-model uncertainties per study.
In this article, we follow this line rigorously by quantifying the modeling uncertainty for 19 different solution strategies against 101 experimental benchmarks, leading, following a full factorial set-up, to 1919 nonlinear finite element simulations. The set of benchmarks comprises reinforced and prestressed beams, published in literature; the focus is on the load resistance and on the failure modes.
The set of solution strategies comprises only standard approaches, popular in engineering practice, systematically categorized into constitutive models, kinematic formulations and equilibrium settings. They are all based on displacement-based, incremental-iterative nonlinear finite element analyses. They all use a 2D plane stress approach, adopt smeared cracking (and crushing) for the concrete and either use embedded truss reinforcements with perfect bond, or adopt separate truss or beams elements plus interface models to simulate bond-slip between concrete and steel reinforcement. All of them use arc-length control, with the purpose of avoiding displacement control, which is seen as more stable but less generally applicable in practice.
The different solution strategies are selected such that they address current dilemmas in engineering practice. First, with respect to the material modeling, a main variable is the use of a fixed crack or a rotating crack model for the concrete. Fixed smeared crack models describe the constitutive behavior along axes that are fixed upon crack or crush initiation. Progressive hardening and softening are then formulated in the direction perpendicular to the fixed crack or crush plane, with a shear retention factor or function being introduced to describe the reduced shear capacity of the cracked material as compared to the uncracked elastic material. Rotating smeared crack models describe the constitutive behavior in the continuously rotating principal stress-strain system. These models can be conceived to include an implicit shear term that guarantees coaxiality between principal stress and strain. Guidance for selecting these crack models is limited. Second, with respect to the kinematics, the strategies are distinguished by using different element sizes. A typical dilemma in engineering applications, which involves a trade-off between economy and accuracy, is the size of the model and thus the feasibility of the numerical model. Also, the presence of a bond-slip relationship and the inclusion of beam action in addition to the truss action in the reinforcing bars is a common issue because they are linked to the possibility for specific mechanisms to arise, which could be critical for the target failure mechanisms, but can also potentially complicate the analyses. By selecting bond-slip, it naturally follows that we should also focus on selecting an adequate bond-slip material model. Thirdly, with respect to equilibrium, it is noted that the, sometimes, "subtle" aspect of obtaining "fully" converged solutions is often overlooked by the users while it may clearly influence the prediction results, and is for this reason explicitly addressed in the selected strategies in this study. The strategies are distinguished by using different convergence norms and a different maximum number of equilibrium iterations per load step. It is noted that the paper is limited to 2D representations of constitutive and kinematic behavior. 3D effects like confinement and detailed localizations have not been explicitly modeled.
The main objective is to provide estimates of the statistical properties of the modeling uncertainties for each modeling strategy, in view of the application fields of reinforced and prestressed beams. First, we assess how well the simulated failure modes correspond to the experimental counterparts; next we quantify the accuracy of the simulated load resistance. For the latter, the modeling uncertainty is defined traditionally as the ratio of the experimentally found load resistance and the numerically simulated resistance. 12 For a limited number of analyses, we present detailed analysis results, which are envisioned to be a strict requirement for any practical nonlinear finite element application, in order to interpret the found statistical properties. The outlook is that these estimates of the statistical properties of the modeling uncertainties could serve as a basis within safety formats, based on either the partial factors method, the global factor method or the simulationbased method, adopting the terminally used in the upcoming fib Model Code 2020. 3

| SELECTION OF BENCHMARK EXPERIMENTS ON REINFORCED CONCRETE BEAMS
To provide objective and widely applicable results, a diverse selection of reported experiments on reinforced and prestressed beams was used to investigate the model uncertainties and accuracy of failure mode modeling. The diversity related to beam dimensions, materials and the presence or nonpresence of stirrups, leading to a wide range of failure modes varying from bending, to diagonal tension and shear compression. Beams, that are often used for benchmarking, have been included. Additionally, beams from other research projects were added to cover a wide application field, including beams that are likely to occur in daily engineering practice of design and assessment. The final set of benchmark beams is listed in Table 1.
This amounted to 101 benchmark experiments, see de Putter 33 for details. Their reinforcement ratio varied from 0.35% to 5.31%, the shear span over depth ratio (a/d) varied from 1.1 to 7.0. The total depth varied between 90 and 1600 mm and the concrete compressive strength varied between 18 and 126 MPa. The latter is important as all other concrete and bond-slip parameters were derived from correlations to the reported compressive strength. Mainly three-point and four-point loading was used. The selected set of experiments included statically indeterminate and cantilevered beams as well. The distributions of the main quantities over the benchmarks are depicted in Figures 1 and 2. It is noted that one beam was left out of this plot due to its depth of 1712 mm. An overview of the experimental failure modes and reinforcement configurations is presented in Table 2.

| NONLINEAR FINITE ELEMENT APPROACHES FOR THE ANALYSIS OF REINFORCED CONCRETE BEAMS
A NLFEA solution strategy comprises all choices and assumptions regarding constitutive relations, kinematic compatibility, and equilibrium conditions. 34 The selected solution strategies for the present study are based on guidelines issued by the Dutch Ministry of Infrastructure and the Environment 14,35 and include variations thereof, resulting in 19 solution strategies. These guidelines are software independent. It is noted though that in practice, a particular strategy may not be transferred one-to-one when switching to another software code: the specific software implementations may contain small subtle, differences. However, the emphasis in the current article is on the great importance of end-user decisions when making a model with a selected software; Diana was used for all analyses. The correctness of the constitutive and FE implementations for single-element and structural component models is shown in a verification manual, which partly reduces the implementation uncertainty. 36 The relatively high number of strategies were motivated by the conviction that selected choices from the constitutive, kinematic and equilibrium classes will interact and must therefore be tested in conjunction.
The common basis for the solution strategies is summarized in Table 3 failure load is retrieved from the analyses. In many practical applications of NLFEA of concrete structures, nonconverged steps are unavoidable. Local effects and successive cracking potentially introduce temporary convergence issues-but should not be interpreted as failure.
To avoid biased interpretations of numerical results, a secured interpretation was used in which the failure load was defined as the highest retained load with at least one of the convergence criteria satisfied. Second, in order to interpret the simulated failure mode, the following assumptions were made. A flexural failure displays either yielding of the longitudinal reinforcement and crushing of the compression zone or just excessive yielding of the longitudinal reinforcement. Shear compression failure was marked by crushing of the diagonal strut and flexural shear was identified as failure due to opening of the critical inclined crack without failure of the compression zone. In some cases, excessive yielding of the stirrups leading to divergence of the analysis was observed numerically, termed "stirrup failure" in this article, in which the stirrups cannot provide sufficient plastic deformation till the crushing of the concrete. All material parameters were derived from the reported compressive strengths, following the formulas from the fib Model Code 2010 44 and following Nakamura & Higai 45 for the relation between the compressive and tensile fracture energy. Finally, all analyses were performed using the DIANA finite element program, release 10.4. 36 Table 4 provides an overview of the 19 solution strategies. The 19 solution strategies differ from each other on six aspects. Where all strategies used a total-strain based smeared crack approach, which is rather popular among practitioners, see Table 3, the first distinction between the different strategies was the adopted concept for the crack orientation within the constitutive model. Already in the 80s, for instance in Rots & Blaauwendraad, 46 there were lively discussions on whether to use fixed or rotating crack modeling strategies for such models and to which extend this depends on the application at hand; a discussion that to date was never fully concluded. In this study a damage-based shear retention model is used for the fixed crack model, in which the shear modulus along the cracked surface is reduced proportionally to the reduction of the secant stiffness of the crack in opening direction. 47 Such reduction of the shear modulus with increasing crack normal strain represents physical reality, as tortuous cracks provide less shear stiffness once they open up wider. The particular damage-based function adopted has proven to serve well in a variety of RC cases, whereas other shear retention models such as the assumption of a constant shear retention factor may lead to larger sensitivities. The element sizes depended on the height h of the beams by using a fixed division of the height. The reinforcements either only included axial action, in embedded truss reinforcements and truss elements, or also included bending action, in reinforcement beam elements. The bond between reinforcement steel and concrete was either considered as fully bonded, or followed shear stress-slip relations by Shima 48 or the fib Model Code 2010. 44 The same bond models have been used after prestressing of the strands. For the stirrups, embedded truss reinforcements have always been used. For some strategies the maximum number of equilibrium iterations per load step were increased. As temporary nonconverged steps were accepted, the maximum number of iterations were expected to influence the results. The convergence norms were based on energy (E) and on the out-of-balance forces (F). 36 The tolerances for the energy and force norm ratios have been indicated in the table. The logical OR indicates that convergence is obtained if one of the two norms is reached, whereas AND denotes that both norms should be satisfied before convergence is obtained.

| QUANTIFICATION OF THE ACCURACY OF PREDICTED FAILURE MODES
As will be shown in Section 5, the concrete model selection has the most significant impact on the accuracy of the load resistance prediction. Moreover, Section 5 will reveal that the solution strategy F9 is one of the better strategies, simulating the load resistance most accurately. For this reason, before presenting these modeling uncertainties for all strategies in Section 5, this section will first show the accuracy of the simulated failures modes of strategies F9 and R8. Note that the only difference between these strategies is the crack model, which is fixed versus rotating. Table 5 compares the failure modes that were observed in their respective experiments, to the failure modes that were modeled in this study. The wellpredicted failure modes, along the diagonals for R8 and F9, have been set in bold. For the rotating crack strategy R8, 24 experimentally observed flexural failure modes were simulated as compression shear failure. For the fixed crack strategy, 13 experimentally observed bending failure modes were simulated as stirrup failure. These instances are underlined in the table and will be discussed in the remainder of this section.
First, the 24 wrongly predicted failure modes for the rotating crack strategy can all be exemplified by the predicted results of test B702A1 from Yang & Koekoek. 25 Figures 3 to 6 show numerical results. The forcedisplacement characteristic in Figure 3 shows the experimental behavior as a solid line and the NLFEA prediction, using strategy R8, with a dashed line. At a load level close to the reported failure load, a drop in the retained load (marked as "Critical event") was observed in the model. Figure 4 mainly shows the crack normal strains, indicating the local crack orientations per element. In order to arrive at crack widths, these strains should be integrated over the crack band width which is related to the dimensions of the finite elements, or even of several elements in case the strains are not fully localized. However, in this article, we do not focus on crack width predictions and the plot of the crack strains is thought to be sufficient to give an impression of the crack pattern. Before the drop in the load bearing, the crack pattern as depicted in Figure 4 on the left was simulated. Afterwards, at 210 kN load level, the crack pattern on the right was observed. It shows how the position of the critical shear crack changed from its original position, which crosses the compressive strut, to its final position being under the strut. The over-rotation of the crack pattern, on the right, enables the possibility of having an alternative shear force transfer mechanism through the compression strut, thus, increases the load bearing capacity. Figure 5 compares the predicted crack pattern at 174 kN from Figure 4 to the experimentally observed crack pattern at failure. 25 For completeness, the crack pattern predicted with strategy F9 is also presented. It shows that the initial prediction of the critical crack was correct in both strategies. Figure 6 monitors the simulated crack opening in the element marked in Figure 4.
After an instantaneous increase of the crack opening at 174 kN, the sustained load decreased after which the crack closed, and load increased again.
In the case of the rotating crack strategy R8, this type of behavior was observed in most models of beams without shear reinforcement. After opening of the flexural shear crack, rotation of the cracks accommodated an alternative equilibrium state leading to a shear compression failure at a substantially higher load level. For 24 of the 32 cases, Table 5, the models all developed a flexural shear crack, after which the formation of a strut was observed until a compression shear failure was found. Convergence behavior indicated problems around the steps where the equilibrium state changed, after which the steps converged as before. In the other eight cases, the opening of the critical crack led to numerical divergence and forced the analysis to terminate.
Second, the 13 wrongly predicted failure modes for the fixed crack strategy F9, Table 4, can all be exemplified by the predicted results of test B321 from Rashid & Mansur. 24 The fixed crack models of beams with stirrups all showed premature stirrup-or shear failure. As soon as diagonal cracks were formed, fixed in direction and quickly losing capacity due to the conservative shear retention model, the stirrups could not be utilized fully, resulting in a stirrup failure before the full shear or flexural failure could be reached. This is illustrated in Figures 7 and 8. Figure 8 shows that the load resistance is substantially under-predicted by the fixed crack solution strategy F9. Figure 7 illustrates the premature stirrup failure, which is exemplary for all investigated fixed crack solution strategies, indicating that beams with stirrups could better be modeled with rotating crack models. In Table 5, this behavior has been classified as "stirrup failure", because the main cause of the observed failure was just that. This specific failure mode was not reported by the experiments as it is a spurious numerical behavior. It is noted that better results for the beams with stirrups could also have been obtained by applying a less conservative shear retention model to the constitutive model for beams with shear reinforcement, thus increasing the shear capacity of the cracked concrete and preventing the premature failure modes from occurring. This was not in the scope of the present study.

| QUANTIFICATION OF MODEL UNCERTAINTIES FOR THE LOAD RESISTANCE
This study set out to provide estimates of the statistical properties of the modeling uncertainties and identify the uncertainty for different solution strategies. To that extend a total of 1919 instances of the modeling uncertainty were recorded, 101 for each of the 19 different solution strategies. The distribution parameters of the  210 kN (false and "over-rotated", right). A single crack orientation is shown per finite element within-model uncertainty for each solution strategy were determined, as in Engen et al. 12 In addition, as the previous section revealed issues with predicting the right failure modes for beams with shear reinforcements for rotating crack approaches and for beams without shear reinforcement for fixed crack approaches, a division was made in the population of benchmarks, also considering beams with and without shear reinforcement separately. Table 6 shows the resulting distribution parameters for all the 101 beams as a single data set as well as the results for the 57 benchmarks with and 44 without shear reinforcement as separate data sets. As discussed in the previous section, rotating crack models could predict an alternative equilibrium state, resulting in wrongly predicted failure modes. For this reason, solution strategy R8* was added, which uses the NLFEA results of strategy R8, but for benchmarks expected to experience a flexural shear failure, the failure load was no longer taken to be the maximal sustained load, but the load at moment where the flexural shear crack opened as identified in Figure 3. This makes strategy R8* as the only solution strategy which includes an explicit manual interpretation of results.
The mean values, μ θ , indicate that all solution strategies are to some extend biased, ranging from the nonconservative value 0.893 when rotating crack solution strategy R6 is used for beams without stirrups to the conservative value 1.322 when fixed crack solution strategy F1 is used for beams with stirrups. The reason for the conservative results of fixed crack models was already discussed in the previous section. In all cases, the mean values are considered as acceptable, especially since they are accompanied with explanations based on the predicted failure modes. Since a known bias can be easily accommodated in any safety format, the focus of our discussion on the results in Table 6 is now moved to the coefficients of variation, V θ .
If a single solution strategy for both beams with and without stirrups has to be recommended, the added solution strategy R8* clearly outperforms all other solution strategies with a coefficient of variation as low as 0.118 for the total data set. If a solution strategy for beams with stirrups is sought, notably all the rotating crack strategies perform better than the fixed crack strategies. And finally, if a solution strategy for beams without stirrups is sought, the fixed crack models in general perform better than the rotating models. For this division of beams without stirrups, the relevance of carefully interpreting the failure mode, when using rotating crack models, is clearly demonstrated from the result of strategy R8*. Only in case of strategy R8* the coefficient of variation reaches the relatively low value of 0.141.
Upon closer inspection of Table 6, the influence of finer meshes and tighter convergence criteria can be analyzed. Comparing strategies F2, F5, and F7, with, respectively, 20, 30, and 40 elements over the height, shows that the fixed crack model was somewhat mesh-sensitive, where a finer mesh resulted in higher resistance predictions, especially for deep beams. Comparing strategies F3 and F6 shows the relevance of the convergence criterion. The price of using a loose convergence criterion, as in strategy F6, is a relatively high coefficient of variation of 0.235. F I G U R E 6 Crack opening (and closure) in the element marked in Figure 4 for rotating crack strategy R8. The force on the y-axis is the load on the beam, used to indicate the progression of the analysis F I G U R E 5 Experimental result B702A1, without stirrups, 25 with the simulated crack pattern (overlay) for rotating cracks strategy R8 (left) at 174 kN and fixed crack strategy F9 (right) at 171 kN The bond slip models influenced the results in the following manner. For the beams without stirrups, the best results were found for solution strategy F9, using the fib MC2010 bond slip model. For solution strategy F8, in which using the Shima bond slip model is the only difference, both the bias and the coefficient of variation were found to be larger. With 40 iterations per step, solution strategies F1 (fully bonded), F2 (Shima bond slip model), and F4 (fib MC2010 bond slip model), the difference between the Shima model and fully bonded was less pronounced and the fib MC2010 bond slip model resulted in substantially better results for beams without shear reinforcements. The V θ of 1.40 for strategy F1 is however the highest coefficient of variation found within the family of fixed F-strategies for beams without shear reinforcement. For the rotating R-strategies for beams with shear reinforcement, the variation of the bond slip model did not result in remarkable differences in the coefficients of variation. It can be concluded that the modeled failure mechanisms in beams without shear reinforcement are sensitive to the bond slip model, while those of beams containing stirrups are not. This can be attributed to the different physical mechanisms observed in beams with and without stirrups. Bond slip modeling improves the crack pattern, crack spacing and the predicted crack openings. The flexural shear mechanism in beams without stirrups heavily depends on the opening and location of the critical crack, and thus modeling of bond-slip pays off. In beams with stirrups the ultimate limit state is governed by crushing of the compression struts and yielding of the stirrups rather than cracking. Consequently, a detailed prediction of crack locations using bond-slip models is less relevant for beams with stirrups.  Families of solution strategies can be distinguished based on their most important features; the choice of the crack model is a clear example of this. Solution strategies may be termed robust in case the results are only slightly sensitive to variations of parameters that are secondary in the solution strategy. In this study we consider the kinematic and equilibrium settings as secondary parameters. Table 6 shows that the family of fixed crack strategies for beams without stirrups is robust, with little variation of the obtained coefficients of variation. The same holds true for the family of rotating crack strategies for beams with stirrups. Figure 9 presents this graphically via boxplots. The boxplot "All data" is generated as follows. For each beam the 19 solution strategies lead to 19 modeling uncertainties θ, of which the standard deviation of ln(θ) is calculated. The boxplot shows the distribution of these standard deviations considering all beams. The other boxplots are generated in a similar way, considering only fixed or rotating crack models, combined with considering only beams without or with stirrups. Figure 9 underlines the robustness of rotating crack models for beams with stirrups and shows that selecting solutions strategies for beams without stirrups should be done with care.
To study the statistical distribution of the modeling uncertainty, Figure 10 shows the frequency distribution of all 1919 analyses with its log-normal fit. Most safety formats make use of the assumption that the modeling uncertainty follows a log-normal distribution. Because of the lack of large sample sizes, this assumption has never been thoroughly checked. In this study, because of the large amount of samples, we may verify the assumption making use of the Shapiro-Wilkes test. On these 1919 values of ln(θ i ) the Shapiro-Wilkes test rejected the H 0 hypothesis of a log-normal distribution, using a p value <0.05, see Table 7. Also, when using the results of solutions strategies R8 and F9 separately, the Shapiro-Wilkes test on these subsets rejected the hypothesis. However, the hypothesis of a log-normal distribution could not be rejected in case of solution strategy R8* for all beams, or when using strategy F9 solely for the beams without stirrups, or when using strategy R8 solely for the beams with stirrups. This supports the findings on correctly modeling the failure modes for these solutions strategies, see Section 4.  Finally, Figure 11 relates the found instances of the ductility index and modeling uncertainty for the two strategies F9 and R8*. This ductility index was proposed by Engen et al. 12 and is computed as a by-product of every NLFEA. Defined as the ratio of the plastic dissipation of the reinforcement to the plastic dissipation of the system, at failure, the ductility index, takes values between 0 and 1, and indicates to which degree the failure mode is governed by the ductile reinforcements. The figure shows a higher uncertainty for a ductility index, X ductility , lower than 0.6. This supports the observation made in Engen et al. 12 For practical NLFEAs, where the identification of the failure modes can be more difficult, the output of the ductility index can thus be helpful when interpreting the results.

| CONCLUSIONS
Nineteen NLFEA solution strategies were benchmarked against a wide variety of 101 experiments on reinforced and prestressed concrete beams, with and without stirrups, failing in bending, flexural shear, or shear compression. The 1919 results were then used to investigate the accuracy of the modeled failure behavior and the resistance modeling uncertainty. The former was judged upon in a qualitative sense, discussing type and brittleness or ductility of the failure modes, while the latter was treated quantitatively, with means and coefficients of variations for ultimate loads computed.
The 19 solution strategies differ on three aspects, classified according to concrete constitutive modeling choices, kinematic modeling choices and equilibrium modeling choices. It was noted that although the solutions strategies were presented as software independent, a particular strategy may not be transferred one-to-one when switching to another software code; this was not further investigated and Diana 36 was used for all analyses. The emphasis in the current article was on the great importance of end-user decisions when making a model with a selected software.
The constitutive modeling choices for the concrete appeared as major, in particular the choice for a fixed or a rotating smeared crack/crush model, both being popular models among RC practitioners. One may argue that the fixed crack model is in general a superior solution strategy for both types of beams, with and without shear reinforcement, since the rotating crack model does not consider additional action, such as an explicit shear resistance of cracked concrete. On the other hand, it can be argued that the performance of a fixed crack model will depend heavily on the specific shear retention model adopted and the associated risk of over-estimating the load capacity. Selecting a suitable shear retention model was not part of this study. The results shown in this article reveal that for the benchmarks with shear reinforcement the rotating crack models were found to be superior, while fixed crack models performed better for the benchmarks without shear reinforcement. For beams with shear reinforcement, rotating crack models properly accommodate the rotation of compression struts and subsequent redistribution and reserve, while fixed crack models do not accommodate that and may lead to premature failure of stirrups. For beams without shear reinforcement it is the other way around; here, the failure is governed by a sudden brittle crack, properly captured by fixed crack models, while rotating crack models then falsely allow for over-rotation of cracks and struts, building up too much capacity and overestimating the failure load. This calls for modeling choices that are tailored to certain structural types, in this case structures with and without stirrups. These aspects were demonstrated to be in line with the quantitative uncertainty statistics, revealing means much closer to the experimental values and variations much smaller once tailoring was adopted as compared to the entire pool of solution strategies. It was furthermore demonstrated that the unwanted effect of the rotating crack model for beams without stirrups could be mitigated by carefully postprocessing the results in terms of manually identifying a "critical event" associated with brittle failure.
The kinematic and equilibrium modeling choices appeared to be secondary for the benchmarks considered. Kinematic assumptions like the inclusion of bond-slip F I G U R E 1 1 Modeling uncertainty versus ductility index for fixed crack solution strategy F9 (in blue) and rotating crack solution strategy R8* (in black). A square denotes that a brittle behavior was observed by interpreting the NLFEA results; a triangle denotes that a ductile behavior was simulated with the fib MC2010 bond slip model did improve especially the fixed crack model predictions for beams without stirrups, because of more accurate representation of the crack pattern. Mesh fineness did not have a major effect, probably because the softening curves for tension and compression were regularized against the element dimensions. Equilibrium assumptions like the choice of convergence norms and maximum number of iterations per load step appeared to affect the results, but again not to a similar extent as the rotating/fixed constitutive assumptions. In general, stricter norms and a larger number of maximum equilibrium iterations help to identify the proper failure modes and loads, as demonstrated by the summary table of uncertainty statistics (Table 6). However, the use of small steps and an arc-length procedure for tracing postpeak load paths, as adopted in all solution strategies herein, is of equal importance.
It was shown that specifically for the optimal rotating crack solution strategy (R8) for the subset of beams with stirrups, and for the optimal fixed crack solution strategy (F9) for the subset of beams without stirrups a Shapiro-Wilkes test did not reject the hypothesis that modeling uncertainty is a log-normally distributed random variable. The same held for the adjusted solution strategy of rotating cracks with manual identification of critical failure (R8*) for all beams. In other cases, this hypothesis was rejected.
It was confirmed that there is a relevant link between the model-uncertainty and the ductility index. This follows from the general observation that quasi-brittle shear behavior, especially for beams without stirrups, is more complex and therefore more sensitive to subtle choices in the solution strategy.
All conclusions point to the same direction. Beams with stirrups are well modeled with rotating crack models. These models are robust in the sense that the results are hardly influenced by changes that are secondary in the solution strategy, like the bond-slip model, the equilibrium conditions and the mesh size. Beams without stirrups are well modeled with fixed crack models, but require more attention.
The outlook is that the estimates of the statistical properties of the modeling uncertainties could serve as a basis within safety formats. Often, multiple modeling approaches are compared on a single benchmark (in Round-Robin tests), or a single modeling approach is compared for multiple benchmarks, but the systematic study of multiple modeling approaches (in this case 19) on multiple benchmarks (in this case 101) is thought to help further improve the reliability of NLFEA, optionally with tailored subfamilies of approaches for specific structural classes, covering both between-model and withinmodel uncertainties.