A computational strategy to establish algebraic parameters for the Reference Resistance Design of metal shell structures

The new Reference Resistance Design (RRD) method, recently developed by Rotter [1] , for the manual dimensioning of metal shell structures effectively permits an analyst working with only a calculator or spreadsheet to take full advantage of the realism and accuracy of an advanced nonlinear finite element (FE) calculation. The method achieves this by reformulating the outcomes of a vast programme of parametric FE calculations in terms of six algebraic parameters and two resistances, each representing a physical aspect of the shell’s behaviour. The formidable challenge now is to establish these parameters and resistances for the most important shell geometries and load cases. The systems that have received by far the most research attention for RRD are that of a cylindrical shell under uniform axial compression and uniform bending. Their partial algebraic characterisations required thousands of finite element calculations to be performed across a four-dimensional parameter hyperspace (i.e. length, radius to thickness ratio, imperfection amplitude, linear strain hardening modulus). Handling so many nonlinear finite element models is time-consuming and the quantities of data generated can be overwhelming. This paper illustrates a computational strategy to deal with both issues that may help researchers establish sets of RRD parameters for other important shell systems with greater confidence and accuracy. The methodology involves full automation of model generation, submission, termination and processing with object-oriented scripting, illustrated using code and pseudocode


Introduction
Shells are widely recognised to exhibit the greatest complexity of any structural form, a consequence of a behaviour governed by local buckling, global nonlinear collapse, plasticity and sensitivity to various imperfection forms, in addition to a myriad of possible geometries and load conditions [2] . The power of modern computing, coupled with the recent development of a powerful framework to characterise the results in a uniform manner, have led to the development of Reference Resistance Design (RRD) [1] . Without this framework, every problem must be treated in an ad-hoc manner. However, RRD requires a huge volume of nonlinear calculations required to support it. This paper presents a new and highly efficient technique to achieve that goal.
Although early advances in the analysis of shell buckling problems were made through algebraic studies [3][4][5] , it did not take long for the set of shell structural systems that could be processed design of metal shells under the generic title 'design by global numerical analysis' [14] . The first and simplest method is 'LBA-MNA' [15] , where the designer must use a computer calculation to determine the reference elastic critical buckling and plastic collapse resistances of the perfect shell from a linear bifurcation (eigenvalue) analysis (LBA) and an ideal elastic-plastic materially nonlinear but small displacement analysis (MNA) respectively. The second method is termed 'GMNIA' [16,17] and involves full geometrically and materially nonlinear analyses with explicitly modelled imperfections. Analyses of intermediate complexity are also defined in this Standard, termed GNA, GMNA and GNIA, where 'G' and 'M' indicate the inclusion of geometric and material nonlinearity respectively, while 'I' indicates the presence of imperfections. Neither of these computer-aided methodologies come easily to many practitioners, who may lack the training, the time or the software necessary to undertake nonlinear FE analyses on anything but the most important of designs. Indeed, the ECCS European Design Recommendations on the Buckling of Metal Shells [14] offers a commentary on EN 1993-1-6 and devotes much discussion to how such analyses should be conducted, illustrating many of the potential pitfalls that can make a calculation go wrong or be badly misinterpreted.
The new method termed Reference Resistance Design (RRD) was recently devised by Rotter [1,2,18] to permit analysts to take advantage of advanced computational predictions without having to perform these themselves. RRD is built around a capacity curve which describes the base relationship between a shell's nonlinear characteristic resistance normalised by its reference plastic collapse resistance ( χ = R k / R pl ) and its dimensionless slenderness ( λ = Ý( R pl / R cr ), where R cr is the reference elastic critical buckling resistance). The entire relationship is characterised by six independent algebraic parameters α G , α I , β, η, λ 0 and χ h , each one containing information about a real physical effect as described shortly. The most recent formulation of the capacity curve [17] is given in Eq. 1 and illustrated schematically first as a classic buckling curve in Fig. 1 a and again in a modified form [19] in Fig. 1 b. A detailed account of background to the capacity curve and proposed further enhancements may be found in Rotter [15] and Doerich and Rotter [20] .
In Eq. 1 , α G and α I are separate adjustments on R cr to respectively account for the influences of geometric nonlinearity and imperfection sensitivity in slender shells that fail elastically. As most shell systems under simple load conditions suffer from prebuckling geometric softening and detrimental imperfection sensitivity, α G and α I are usually less than unity. However, more complex systems such as cylindrical silos under unsymmetrical eccentric discharge pressures may exhibit beneficial changes of geometry leading to geometric stiffening and/or strengthening 'imperfections' [21,22] , such that α G and α I may become greater than unity. In theory, this very general formulation supports both possibilities. Although α G ≥ 1 cannot be avoided if the beneficial geometric nonlinearity relates to the fundamental mechanics of the system, a situation where α I > 1 usually signifies that an inappropriate 'imperfection' has been chosen and it is likely that a more carefully-chosen deviation from the perfect geometry will lead to α I < 1 [22] .

The
'plastic range factor' β identifies the point ( λ, χ) = ( λ p , 1 -β) on the capacity curve at which plasticity first begins to have a noticeable effect on the behaviour and marks the onset of the 'elastic-plastic' region. The 'interaction exponent' η controls the curvature of the relationship in this region of Fig.  1 a, with η < 1, > 1 and ≡ 1 allowing convex (positive curvature, indicating a strong elastic-plastic interaction), concave (negative curvature, indicating a weak interaction) and idealised linear (zero curvature) relationships respectively in the χ vs λ plane. The 'squash limit' λ 0 defines the slenderness at which the reference full plastic conditions R pl is reached, while χ h ≥ 1 defines the projected intercept of the vertical axis for a fictitious shell with zero slenderness, allowing for resistances higher than R pl for very stocky shells due to geometric or strain hardening.
In effect, RRD permits an analyst working with only a calculator or a spreadsheet to take full advantage of the accuracy and realism of a GMNIA, on the understanding that the six algebraic parameters that account for the complex nonlinear phenomena have already been established a priori for all practically relevant parameter ranges. While it is possible in theory to establish these on the basis of laboratory testing, the expense involved in doing so for all systems and parameter ranges would be absolutely prohibitive and difficult to justify in the modern age of limited budgets, powerful software and cheap computing power. A computational programme using a validated and robust FE model remains the only defensible option. The significant reduction in uncertainty in the structural resistance model obtained through RRD should reduce the need for excessive compensation through high partial safety factors, permitting more economic designs.
The RRD method has now been accepted for use in design through an Amendment to EN 1993-1-6 [13] coupled with the creation of a new normative Annex E that will hold RRD expressions for algebraic parameters and other information [18,23] . Two systems have currently been approved based on available research, cylindrical shells under uniform global bending and spherical shells under uniform external pressure, but the call is out to the research community to gradually extend this set to cover new geometries and load cases. In particular, Rotter [2] lists a number of important systems that are ripe for implementation in the RRD format, including tanks and silos with stepped walls, lap-jointed construction, wind loading, stiffened and corrugated shells, localised loading, conical shells under external pressure, cylinders with openings and cutouts, shells with ringbeams and others.

Scope of the study
Although RRD and the capacity curve together form a robust framework for the manual dimensioning of shell structures, the challenge of establishing a database of algebraic parameters for the most important geometries and load cases is formidable. A conceptual strategy for doing so is presented in Section 6 of [2] , which readers are warmly invited to consult. The aim of this paper is to complement the above with a practical computational strategy on how to design parametric studies for RRD specifically and perform the necessary calculations. In particular, the authors take the view that in order to establish such a database, the international shell buckling research community should embrace the full power of automation now widely available in modern computing to aid in the generation, submission, termination and processing of the many thousands of necessary parametric LBA, MNA, GNA, GNIA and GMNIA studies. The powerful commercial ABAQUS 6.14-2 [11] finite element package is employed here which links with the object-oriented Python and compiled FORTRAN programming languages for extended functionality, although the logic contained in the code and pseudocode examples given may readily be adapted to any other software with a scripting facility. The use of ABAQUS for illustration is not inappropriate, as it has arguably been the most widely used finite element package in the past two decades for shell buckling research (e.g. [7,[24][25][26][27][28][29][30] , among many others).
An outline Python code is presented here, based on a novel conceptual analysis hierarchy, for scripts that automate the mass generation of FE model input files in an ordered manner. It is then shown how to arrange these model files within ordered subdirectories based on the varied system variable, permitting automated submission to the FE solver in simultaneous batches of approximately equal numbers of models. Fully-functional FORTRAN code for ABAQUS user subroutines that permit the automated termination of nonlinear equilibrium path-tracing 'Riks' analyses once specific criteria are met, a vital capability, is presented in its entirety in Appendix A . Finally, an object-oriented post-processing methodology for the automated extraction of the RRD algebraic parameters is presented. The intention at this stage is to show how these can be obtained as accurately and as completely as possible, even if ultimately only simplified relationships are processed into a RRD design rule. For illustration purposes, reference is made throughout to the structural system that has, to date, possibly received the most research attention in this regard, cylindrical shells under uniform bending. This system finds ubiquitous application in generic load-bearing circular hollow sections, piles, pipelines, chimneys, wind turbine support towers and tall slender silos that may be subject to transverse loading, wind, seismic or other actions that induce a global bending moment.

The analysis hierarchy
It is proposed here that any metal shell system being characterised can be treated abstractly as a data tree according to an ordered hierarchy of analysis Levels. A physical variable, varied within a particular Level, is held constant when ascending to a higher Level along the same branch. A flowchart illustrating this is shown in Fig. 2 , customised for the paper's recurring example system of cylinders under uniform bending. The capacity curve upon which RRD is based ( Fig. 1 ) is formulated in terms of the resistances R cr and R pl which are linear in the elastic modulus E and the yield stress σ y respectively. It may be argued that these can be kept constant at reasonable values for the entirety of the characterisation and constitute the base state, or Level 0. Changes to the structural form or the loading distribution constitute changes to the system and require a start from a different Level 0.
Recent research has begun to acknowledge and incorporate the beneficial effect of strain hardening on steel design through simple piecewise-linear characterisations of the stress-strain relationship [31][32][33][34] . Without strain hardening, a system can never attain the full plastic condition at a finite slenderness, unrealistically forcing λ 0 to be zero and χ h to be unity in any RRD characterisation ( Fig. 1 ). The innately slender nature of shells means that their behaviour is governed mainly by stability considerations, and plastic strains are unlikely to go far beyond the 'yield plateau' portion of the stress-strain curve. The simplest and most conservative material law is thus a bilinear relationship with an initial elastic region of modulus E followed by linear strain hardening region of modulus E h (expressed as a percentage of E ) which should form the Level 1 variable. It may be set to zero for an ideal elastic-plastic characterisation, though a recent study of 225 stress-strain curves by Sadowski et al. [34] found that the 'yield plateau' actually exhibits a statistically significant positive modulus E h of the order of 0.1% of E , a more accurate but still conservative lower bound on E h that supports realistic values of λ 0 and χ h . The Level 2 variable should relate to the relative amplitude δ/ t of the imperfection form, carefully chosen to ensure that it is relevant and appropriately deleterious to the shell system within the constraints of practical constructions [13,14,22] .
The suggested analysis hierarchy permits Levels 1 and 2 to be universal for all metal shell systems, since every shell is affected by plasticity and imperfections. By contrast, Levels 3 and above govern a system's dependency on geometric nonlinearity which is controlled by variables unique to its geometry such as the length, width, radius and thickness (or their ratio r / t for compactness) etc., whose Levels may be arranged in any logical order. While the hierarchy may be followed methodically to perform a vast number of GMNIAs and obtain all RRD parameters by brute force, it is perhaps better used to design sub-studies to obtain individual parameters in isolation, simplifying the overall task. For example, characterising R cr and R pl requires only LBAs and MNAs to be performed at Level 3 and automatically 'assume' that the branch corresponding to E h = 0 and δ/ t = 0 has been chosen at Levels 1 and 2 respectively. Obtaining α G requires Level 3 GNAs under the same variable choices but using large displacement theory, while α I requires Level 2 GNIAs. Different imperfection forms may be investigated, each resulting in a different α I . Since LBAs, GNAs and GNIAs are by definition elastic with plasticity 'deactivated', and since MNAs should not include strain hardening [16,17] , Level 1 is redundant for these sub-studies. However, the remaining plasticity-related RRD parameters β, η, λ 0 and χ h must be established using GM-NIAs performed at Level 1, though potentially on the basis of fewer analyses across smaller parameter ranges if an understanding of other aspects of the system has been developed. The hierarchy formally establishes that any RRD parameter at any Level is a function only of the physical variables at the same Level or above. For ex- ample, in Fig. 2 , R cr , R pl and α G (Level 3) are independent of the imperfection amplitude δ/ t , a lower Level 2 variable, but dependent on r / t and the length L (same Level or higher); α I (Level 2) is independent of E h , a Level 1 variable, but dependent on δ/ t, r / t and L (same Level or higher); while β, η, λ 0 and χ h (Level 1) are dependent on E h , δ/ t, r / t and L . A shell system is only fully characterised once the complete dependencies of all RRD parameters on the system's physical variables have been established across all practical ranges.
Lastly, the two reference resistances R cr and R pl are the anchors of the capacity curve in RRD ( Fig. 1 ) and should be determined as accurately as possible. R cr is fully defined as the lowest positive eigenvalue in a linear bifurcation analysis of a system with compressive stresses, while R pl is defined as the lowest collapse load arising from a rigid (elastic-) ideal plastic calculation [2,17] . Both are based on small displacement theory and are always calculable using LBA and MNA computations respectively. For systems with relatively simple stress states, such as shells of revolution under the fundamental loads of uniform compression, external pressure and torsion, classical analytical expressions may already exist for either limit state. These are often based only on shell membrane theory and a simple criterion of failure i.e. a relevant stress at a single point reaches the classical elastic critical buckling or first yield stress. Such simple predictions are usually only approximations to R cr and R pl since they neglect more complex influences, but they act as a valuable scaffold upon which enhancements can be made on the basis of computational LBAs and MNAs respectively. For thin cylinders under uniform bending, for example, the Algorithm 1 Outline Python code for a Level 3 master generating script with input file logistics, specialised for cylinders under uniform bending (analysis hierarchy in Fig. 2 ).
analytical expression for the classical elastic critical buckling moment is calculated using classical beam theory as: where W el = π r 2 t is the thin-walled elastic section modulus, σ cl = 0.605 Et / r is the classical elastic critical buckling stress, E is the Young's modulus, and r and t are the radius and thickness respectively. However, this simple expression neglects the edge boundary conditions in short cylinders which restrain local buckling in a linear bifurcation analysis, an effect that can only be captured if the full pre-buckling stress state is included in the analysis and a three-dimensional buckle is permitted to develop. Rotter et al. [30] performed hundreds of Level 3 LBAs to show computationally that this boundary effect may be accommodated by a simple adjustment to M cl , giving the following more complete expression for R cr : where ω = L / Ý( rt ) is a carefully-devised dimensionless length that enables R cr to be characterised independently of the r / t ratio. The significance of the ω parameter is explored in greater detail in Section 7 .

Automated model generation and input file logistics
Advanced modern finite element packages often interface with a scripting language to enhance functionality and aid in the automation of parametric studies (e.g. Python in ABAQUS, APDL in ANSYS, Matlab in COMSOL), an invaluable tool for establishing databases of RRD parameters. The analysis hierarchy shown in Fig. 2 lends itself directly to a computational implementation as a system of nested loops, with Level 1 forming the outer loop and each higher Level forming an inner loop. However, this naïve approach may quickly lead to a truly unmanageable number of calculations, as each variation of a low-Level variable spawns an entire branch of potentially thousands of additional GMNIAs. Prudent low-Level choices may allow substantial economies in computational effort with minimal loss of information on RRD parameter dependencies.
It is advantageous to take Level 3 as the outer loop within a master model generating script, with all higher Levels of the system forming the inner loops. In this way, the same master script may be used with only minor adjustments to automate the generation of LBAs, MNAs and GNAs (to individually establish R cr , R pl and α G ) and, with specialisations for (a potentially minimised number of) low-Level variable choices, also GNIAs (for α I ) and GM-NIAs (for β, η, λ 0 and χ h ). An outline of a generating script writ- ten in Python illustrating how the generation of .inp input files (which contain the entirety of the model information that is submitted to the ABAQUS/Standard FE solver for analysis) may be automated is shown in Algorithm 1 . When 'executed' using the example LVL3 vector for cylinders under uniform bending, the code will create 199 directories named L3_10 to L3_1000 (i.e. for r / t = 10 to 10 0 0 in steps of 5), each of which will contain as many .inp files for that Level 3 variable (i.e. r / t ratio) as requested in the LVL4 vector (i.e. length L ). These will all be of the same analysis type (e.g. LBA, MNA, GNA, GNIA, GMNIA) as specified in Analysis .
The script additionally includes commands to aid in the automated submission of the large number of generated .inp files. For a processor with 8 logical cores (4 physical cores at 2 threads per core), for example, setting Batches = 4 requests the creation of four new Python scripts within each Level 3 folder called L4_runme_X.py where X is 0 up to 3, each one with content similar to that shown in Fig. 3 . These may be submitted for analysis in parallel on four separate command consoles (via python runme_X.py ), and each Level 4 .inp file inside its Level 3 folder will be analysed using 2 processor cores (due to CPUs = 2 ). The commands additionally specify that a FOR-TRAN user subroutine named MN_Riks_killer.f or GN_Riks_killer.f be invoked. The vital purpose of these is explained in detail in Section 6 .
For cylinders under uniform bending, the Level 3 r / t ratio directly controls the resolution of points along the capacity curve (see Chen et al. [27] ) and thus the base data set upon which RRD parameters can be accurately established. In contrast with Levels 1 and 2, the Level 3 variable is thus likely to cover a substantial range for any shell system leading to a long LVL3 vector and many Level 3 folders, so it may be convenient to gather all sets of individual Level 4 submissions scripts into a single set at Level 3. The complete Python code given in Algorithm 2 shows how a L3_runme_X.py file can conveniently be built out of any number of L4_runme_X.py files, an example of which is shown in Fig. 4 . This code additionally suggests useful housekeeping operations, such as deleting unnecessary files once an analysis has been completed to preserve disk space. In particular, the capacious ABAQUS output database .odb files (containing the entirety of the computed nodal and element variables) are no longer required once an analysis has terminated, as the buckling loads may be obtained from the much smaller .dat and .sta text files that are created simultaneously.

Introduction
Automation would be of limited use if it were not possible to stop an analysis once a 'criterion of failure' has been met and proceed to the next one. While modern FE solvers are now able to follow nonlinear equilibrium paths deep into the post-buckling domain [35,36] , it is invariably only the first critical point that is of interest for 'design by global numerical analysis' [13] and RRD. Unfortunately, FE software often fail to provide a satisfactory facility for automatically terminating a nonlinear analysis at the first critical point. For shell systems, this is exacerbated by their complexity and the many conditions under which 'failure' may be deemed to occur [17] . This section presents a strategy for the automated detection of all relevant 'criteria of failure' for shells and illustrates its implementation in ABAQUS. It is hoped that eventually these will become standard facilities in all structural FE software with a nonlinear solver.

Automated termination of linear bifurcation analyses (LBAs) and materially nonlinear analyses (MNAs)
The capacity curve ( Fig. 1 ) is defined relative to the reference resistances R cr and R pl obtained from LBA and MNA computations respectively. As all RRD parameters are defined relative to these two reference results, it is important that they be established as accurately as possible. An LBA is a matrix eigenvalue calculation ( * BUCKLE step in ABAQUS) which automatically terminates once the desired number of eigenvalues have been extracted. Assuming that a robust meshing scheme is applied, LBAs require no further monitoring. An MNA, by contrast, is a small displacement analysis to determine the plastic limit state corresponding to a plastic collapse mechanism, manifest as a 'plateau' on the equilibrium path ( Fig. 5 ). This requires an incremental path-tracing analysis, most robustly performed in ABAQUS as a * STATIC step using the 'Riks' [37] modified arc-length algorithm with geometric nonlinearity deactivated and an ideal elastic rigid-plastic material law with no strain hardening.
Doerich and Rotter [38] give a detailed account of the challenges involved in performing MNAs and accurately determining R pl . For shell systems dominated by membrane yielding, the equilibrium path attains an obvious horizontal plateau relatively 'quickly' with only a modest number of analysis increments. However, for systems with localised plastic collapse mechanisms involving extensive local bending and surface yielding, the plateau may not be reachable within a tolerable number of increments ( Fig. 5 ). Doerich and Rotter [38] propose a method to accurately predict the load proportionality factor (LPF) at which the 'plateau' may eventually form (and thus obtain R pl ) based only on a partially-computed portion of a still-rising equilibrium curve of LPF vs |DOF|, the tracked degree of freedom ( Fig. 6 a). The equilibrium plot is first transformed into a Modified Southwell (MS) plot of LPF vs |LPF/DOF| ( Fig. 6 b), which consists of an initial vertical portion corresponding to the linear elastic response and a plateauing portion corresponding to the plastic response. The projection of the tip of the non-vertical portion of the curve on the vertical axis (i.e. a linear fit through the three most recent adjacent data points) gives a first estimate of R pl , called P MS . This is then used in a second transformation to construct a Convergence Indicator Plot (CIP) of P MS vs ϖ = ( P MS -LPF)/LPF ( Fig. 6 c). The vertical axis intercept of a linear least-squares fit through the CIP data gives an enhanced estimate of R pl based on only a modest number of increments, as illustrated in Doerich and Rotter [38] on several shells systems. The equilibrium path, MS plot and CIP are not unique for any system and the choice of tracked DOF requires a careful assessment prior to automation.
The CIP projection method was arguably conceived for application during post-processing to make the most of an MNA that was allowed to run for as long as was deemed tolerable and then terminated manually. Although an experienced user may become a good judge of how many increments are optimal for a good R pl estimate together with the CIP, this is hardly a satisfactory strategy for automation. This paper offers a FORTRAN implementation of an 'on the fly' CIP that can be constructed during an ongoing MNA and automatically terminate it once a Kill Condition (KC) is satisfied. It exploits the extended functionality of user subroutines supported by ABAQUS, compiled during pre-processing, a powerful tool although it requires the proprietary Intel® FORTRAN compiler [39] . This code is contained in a user subroutine referred to in Algorithm 1 as MN_Riks_killer.f and presented in full as Subroutine 1 in Appendix A , based on example code given in Section 1.1.49 of the ABAQUS User Subroutines Reference Guide for version 6.14-2 [11] . The code reads the current LPF and tracked DOF values at each increment (for which a 'Riks' analysis specifically is required) and proceeds through each of the stages shown in Fig. 6 , though it only builds the CIP once a significant deviation from verticality is detected on the MS plot (Coefficient of Variation  [38] CIP projection procedure to accurately estimate R pl from partial MNA data and terminate once a Kill Condition is met (programmed in FORTRAN in Subroutine 1). of the horizontal axis data greater than, say, 10%). A stop command is sent to the ABAQUS/Standard solver once the KC is satisfied that the absolute relative difference between any two adjacent CIP projections of R pl is less than a specified tolerance ( Fig. 6 d), say, 0.01%.

Automated termination of geometrically nonlinear analyses (GNAs, GNIAs, GMNIAs etc.)
Shells are widely known to exhibit unstable buckling behaviour characterised by an abruptly descending and highly nonlinear postbuckling equilibrium path. These may be computed by powerful arc-length 'path tracing' methods, of which the most robust implementation in ABAQUS is again the modified 'Riks' algorithm [37] . The 'Riks' procedure has been shown to be reliable in following geometrically nonlinear equilibrium paths deep into the postbuckling range for various applications [21,24,35,36,40,41] .
In its detailed description of rules governing 'design by global numerical analysis using GMNIA' which attempt to cover all possible metal shell systems, the EN 1993-1-6 [13] Standard defines four 'criteria' that represent a state of failure ( Fig. 7 ). Criteria C1 and C2 represent the most common manifestations of unstable postbuckling, a reversal in the direction of the equilibrium path either due to a limit point (Criterion C1; global snap-through buckling) or a bifurcation (Criterion C2 a ; local buckling). Criteria C3 and C4 represent more ambiguous cases where there is no drop in the LPF and the equilibrium path continues to rise indefinitely, although there may be an intermediate 'kink' where the curve abruptly changes slope. For want of a reliable way of detecting the 'kink', the conservative rule was written to terminate the GMNIA either once a 'largest tolerable deformation' is reached (Criterion C3) or the highest equivalent stress in the system reaches first yield (Criterion C4). Criterion C3 is already possible in FE software which allow a DOF to be tracked and terminate the analysis once the DOF reaches a pre-defined value (which EN 1993-1-6 would supply), while the very conservative Criterion C4 is a poor basis upon which to establish accurate RRD parameters. The remainder of this section thus focuses on the automation of Criteria C1 and C2, as well as on the automated detection of 'kinks' on the equilibrium curve.
The occurrence of unstable post-buckling is manifest by a decrease in the load proportionality factor (LPF), which occurs for both limit points and bifurcations on the equilibrium path ( Fig. 7 ). Unfortunately, there is currently no way through either the ABAQUS CAE interface or the .inp file to automatically terminate a 'Riks' analysis when the LPF first begins to decrease, and to the authors' understanding this facility is also missing by default from other FE software. For shells that exhibit limit point 'snap-through' buckling, a maximum number of increments before termination may be specified based on prior experience. However, this approach is not reliable for bifurcations since the LPF incre- ment may be greatly cut back by the solver in the vicinity of a bifurcation causing the analysis to proceed very slowly, exhausting the permitted number of increments without capturing the bifurcation. While this too is potentially countered by setting a high maximum number of increments, it is a very inefficient solution with no guarantee of success.
The automatic termination of an on-going geometrically nonlinear 'Riks' analysis is a crucial capability that is likely to find application far beyond the present context of RRD. To this end, this paper offers a second FORTRAN user subroutine for ABAQUS that implements Criteria C1 and C2 above in addition to a 'kink detector' . It is referred to in Algorithm 1 as GN_Riks_killer.f and presented in full as Subroutine 2 in Appendix A . Usable with any geometrically nonlinear 'Riks' analysis (GNA, GNIA or GMNIA), the code sends a stop command to the solver when the ongoing analysis is deemed to satisfy any one of the following three KCs: • KC1: current increment's LPF is less than the previous increment's LPF. This captures both limit point and bifurcation instabilities (Criteria C1 and C2), though for the latter only if the solver correctly switches path ( Fig. 7 , ' a ' curve for Criterion C2). • KC2: the change in the LPF from one increment to the next is smaller than a specified tolerance. The authors have found this to be a useful mechanism to prevent the analysis from becoming stuck in an unending loop of cutbacks and increases in the arc-length increment size due to repeated attempts at convergence in the vicinity of a numerically 'problematic' bifurcation. • KC3: the absolute relative difference between the current and previous incremental radii of curvature of the equilibrium curve exceeds a specified tolerance.
KC3, the 'kink detector', is perhaps the most contentious. It is known that some imperfect shell systems do not necessarily exhibit a monotonically decreasing relationship between the buckling resistance and imperfection amplitude, and a bifurcation point with a descending equilibrium path may, for deeper imperfections, turn into a smooth inflection point followed by an ascending path [42] . For example, Fig. 8 a shows equilibrium paths in terms of the LPF against the end rotation and peak inward radial deformation DOFs for a very short cylinder under uniform bending with a 'Type A' Rotter and Teng [43] weld depression imperfection with a deep amplitude at midspan (from Fajuyitan et al. [44] ). Fig. 8 b shows equilibrium paths of the LPF against the end shortening and peak inward radial deformation DOFs for a stepwise varying wall thickness cylindrical silo under eccentric discharge with regularlyspaced weld depression imperfections (the '80%' curve from Figure  13 in Sadowski and Rotter [41] ). Both of these systems exhibit unstable bifurcation buckling at lower imperfection amplitudes where the criterion of failure is obvious and easily detectable by KC1, but at higher amplitudes the bifurcations vanish (with no reported negative eigenvalues in the global tangent stiffness matrix) and there is instead a continuous transition from a pre-to post-buckled state but with a clear 'kink' on the equilibrium path. The deformed shapes both before and after the 'kink' are illustrated in Fig. 9 , and clearly show the development of localised buckling deformations. Detecting the 'kink' has hitherto required careful manual and visual study of the equilibrium path, and it is easily missed if it is not anticipated.
Any implementation of KC3 is complicated by the fact that the equilibrium path is not unique and depends on the tracked DOF, whereas KC1 and KC2 are assessed on the basis of the LPF alone. Tracking a different DOF may cause the smooth transition from pre-to post-buckling appear as a stronger or weaker 'kink' on the equilibrium path (compare left and middle columns in Fig. 8 ). A very simple 'kink detector' would require only the calculation of the second derivative (curvature) of the equilibrium path using a finite difference scheme, coupled with a criterion that tests for a change in the sign of the curvature at any given analysis increment. The change in sign then corresponds to an inflection point, i.e. a change from a softening to a hardening equilibrium path or vice-versa, and is consistent with the criterion used to determine the buckling load of imperfect shells by Yamaki [42] . However, the authors believe that in the most general case, a 'kink' may in fact occur without an inflection point where, for example, a softening path abruptly transitions to another softening path with the tangent stiffness remaining positive and no drop in the LPF. In such a case, a simple 'change in sign' criterion for the curvature would not be sufficient to detect the 'kink'.
The proposed simple implementation of a 'kink detector' computes the local radius of curvature at any increment by fitting a circle to the three most-recent points on the LPF vs DOF curve. In a geometrically nonlinear analysis, even the initial 'linear' portion of the equilibrium curve is never perfectly linear thus the radius of curvature is always finite, although use of double precision arithmetic is advised nonetheless. The proposed detector sends a stop command to the solver if the absolute relative percentage change of two adjacent radii of curvature is found to exceed a tolerance, say 80%. The choice of tolerance should be tailored to the system and the tracked DOF, and should be low enough to be sensitive to the 'kink' but high enough to avoid false positives. This detector makes no assumption about the sign of the curvature, it simply tests for disproportionate changes in it from one increment to another. However, a more sensitive detector could track incremental changes in higher derivatives computed using backward finite difference operators, as in some cases the 'kink' may be too smooth for detection using the curvature alone. In particular, automated termination poses a particular challenge for systems with a naturally stable post-buckling response, such as a rectangular plate under uniform axial load as shown in Fig. 8 c (Benchmark #3 in Sadowski and Rotter [45] ) where the transition from pre-to post buckling is often too gradual to qualify as a 'kink'.
A final issue deserves mention. During a geometrically nonlinear analysis of a perfect shell system, the FE solver may at times fail to switch to the secondary path at a bifurcation and instead continue erroneously on the fundamental path with an increasing LPF ( Fig. 7 ; b curve for Criterion C2). However, the tangent stiff-  ness matrix will cease to be positive-definite and the software may report negative eigenvalues beyond this point, hinting at the presence of a lower-energy state on an ignored secondary path [46] . This problem is traditionally countered by introducing a very minor mesh perturbation in the form of the critical buckling eigenmode (i.e. 1st LBA mode [24,25,47,48] ) which turns the bifurcation into a limit point and ensures a continuous response. Buckling then becomes easily detectable by KC1. However, this procedure does not always reliably ensure a switch at the lowest bifurcation, and certain systems may require such a deep perturbation that the initial tangent stiffness is affected. If the first bifurcation is missed, the analysis risks dangerously overpredicting the buckling resistance. The authors suggest that a further 'KC4' should be im-plemented which terminates the analysis once the tangent stiffness matrix stops being positive-definite and negative eigenvalues are detected for the first time, even if no switch in the load path occurs. However, this does not appear to be currently possible within ABAQUS or other software.

Identification of dimensionless groups governing geometric nonlinearity
Rotter [2,16] emphasises the importance of isolating the effect of geometric nonlinearity ( α G ) for any RRD characterisation, as this permits complete capacity curves to be established for 'constant geometric nonlinearity' where α G is manifest as a straight vertical line in R k / R pl vs R k / R cr space ( Fig. 1 b). The analysis hierarchy ( Fig. 2 ) shows that the effect of geometric nonlinearity may be isolated through a sub-study of GNAs performed at Level 3 and above, leading to a characterisation of α G that is a function of geometric variables only (e.g. r / t and L ). In a number of shell systems investigated so far, a wise application of background theory has enabled the construction of a single dimensionless group, effectively a 'multi-Level' variable, that captures the entirety of α G 's dependence on the shell geometry. Examples include: • Cylinders under axial compression [42,49] : both ω = L / Ý( rt ) and the Batdorf [50] parameter Z = ( L 2 /( rt )) Ý(1-ν 2 ) have been used as dimensionless variables to characterise this system.
It appears that, with the exception of Z , the above are all special cases of a general 'Level 3/4' dimensionless variable that controls where a = 1/2 gives ω, a = −1/2 gives , a = 4/7 gives ξ , while a = 0 or 1 give simply L / r or L / t respectively. It is suggested that α G for many cylindrical shell systems may be expressed compactly in terms of ζ a with an appropriately specialised value of a , and general forms akin to Eq. 4 are likely to exist for other shell geometries. Establishing the governing dimensionless geometric group is an important investment, as it reduces by one the number of variable dependencies on all RRD parameters. For example, Level 1 RRD parameters β, η, λ 0 and χ h then become a function of only three physical variables (i.e. E h , δ/ t and ζ a -directly representing material nonlinearity, imperfections and geometric nonlinearity) , rather than four (i.e. E h , δ/ t, r / t and L ).

Extracting capacity curves at 'constant geometric nonlinearity'
The outline of a script in pseudocode to process the data generated by a programme of GMNIAs performed according to the analysis hierarchy in Fig. 2 is presented in Algorithm 3 . In keeping with the preceding discussion and Algorithm 1 , it is advantageous to perform the processing at Level 3 to avoid the complexity of many nested loops. The code to read buckling load factors from text data files generated by the FE software is easily written and benefits from an ordered arrangement of Level 4 model files within Level 3 folders as suggested by Algorithm 1 and 2 . In ABAQUS, R cr may be obtained as the lowest positive eigenvalue read from a .dat file, while R pl and R k may be obtained as the highest LPF read from a .sta file created by a 'Riks' analysis.
Algorithm 3 suggests making use of the powerful data abstractions available in an object-oriented language such as Python or Matlab [53] to organise the significant amount of information generated during a programme of GMNIAs. The contents of each Level 3 directory of completed GMNIAs (e.g. for constant r / t ) can be stored in a custom-written Level_3 object. These consist of equal-length lists of R k load factors, the corresponding Level 4 variable value (e.g. varying L ) and the constructed Level 3/4 variable value (e.g. varying ζ a ). Continuing the illustration in Algorithm 1 , for example, there would be 199 instances of the Level_3 class associated with the branch choices made at Levels 0 to 2.
The capacity curve ( Fig. 1 ) is a conceptual abstraction that lends itself particularly well for implementation as a data object. Since these should be established for 'constant geometric nonlinearity' [2,16] , an instance of the Capacity_Curve class suggested in Algorithm 3 should be created at each desired value of the constructed Level 3/4 ζ a variable which controls α G . To gather the data necessary to create a Capacity_Curve object at a particular Level 3/4 variable value, it is sufficient to search all Level_3 objects within the same branch of Levels 0 to 2 for GMNIA R k load factors at that Level 3/4 value. This gives sufficient information to construct both a 'traditional' ( Fig. 1 a) and 'modified' ( Fig. 1 b) capacity curve, which may then be used to perform a constrained least-squares fit using Eq. 1 as the functional form and thus establish a set of RRD parameters specific to those Level 0, 1, 2 and 3/4 variable values. In this light, it is advantageous to design the LVL4 vector in Algorithm 1 to always generate the same vector of desired Level 3/4 variable values for any constant Level 3 variable value, as illustrated in the final section.

Illustration on perfect cylinders under uniform bending
The logic of Algorithm 3 was implemented in Matlab [53] and is used here to illustrate how RRD parameters may be established for perfect cylindrical shells under uniform bending. Rotter et al. [30] found that separate phenomena govern the behaviour at dif-  ferent length ranges with α G characterised by two separate Level 3/4 variables, ω ≡ ζ 1/2 and ≡ ζ −1/2 . For brevity, only the behaviour dominated by pre-buckling cross-sectional ovalisation controlled by ≥ 0.5 is illustrated. Fig. 10 a shows an excerpt from approximately three thousand GMNIA R k load factors arranged in a three-dimensional R k / R pl vs R k / R cr vs space, where the resistance R is the applied end bending moment M . These were generated automatically and computed over only 48 hours on a single unsupervised workstation powered by a modest Intel® Core TM i7-6700 processor (3.4 GHz clock speed) with 16 GB of RAM assuming Batches = 4 , CPUs = 2 ( Algorithm 1 ) and automated analysis termination using GN_Riks_killer.f . For the range ≥ 0.5, accurate analytical results exist for R cr and R pl [27] thus they did not need to be obtained using Level 3 LBAs or MNAs. The Level 4 vector of lengths was designed to maintain steps of 0.5 in . The coloured curves on Fig. 10 a parallel to the axis are saved by Matlab in memory as Level_3 objects. The solid black curves in the planes perpendicular to the axis represent 'modified' capacity curves ( Fig. 1 b) saved in memory as Capacity_Curve objects (shown in detail on Fig. 10 b for = 0.5, 1.0, 1.5, 2.0, 2.5 and 3.0). A typical plot of the dependency of the RRD parameters on the Level 3/4 variable in the range 0.5 ≤ ≤ 7 (extracted from 2 × 22 capacity curves) is shown in Fig. 11 for a perfect cylinder and two different initial strain hardening moduli; E h = 0 (ideal elastic-plastic) and 0.1% of E = 200 GPa (see Sadowski et al. [34] ).
An enhanced formulation of the Level 1 RRD parameter η was used [20] which assumes a linear variation with λ between η 0 at λ 0 and η p at λ p . The extracted α G parameter follows the relationship with [30,51] showing that 'Brazier' ovalisation is fully-developed for > 3.5 and responsible for an approximate 50% reduction in the nonlinear elastic buckling resistance. It is entirely unaffected by E h .
The β parameter, representing the first visible onset of the effects of plasticity, is only very slightly affected by E h . As these illustrations have been performed for a perfect shell only, there is as of yet no information to extract α I which would require a full range of Level 2 branch choices. For the case of no post-yield strain hardening ( E h = 0), the χ h and λ 0 RRD parameters were constrained to be fixed at unity and zero respectively. This has a significant effect on the shape of the curves for low λ, which exhibit very different η 0 and λ 0 values depending on whether strain hardening is present or not. A complete RRD characterisation of cylinders under uniform and various non-uniform bending moment distributions is currently being carried out using the methodology presented in this paper, and will be published at a later stage.

Conclusions
Reference Resistance Design (RRD) is a method of analysis recently devised by Rotter [1] and approved for implementation in the Eurocode on the design of metal shells EN 1993-1-6. It permits a designer to take full advantage of the accuracy and complexity of an advanced nonlinear FE analysis of a metal shell without bearing the significant expense of having to perform it. The information about the behaviour of any system is encoded within six physical algebraic parameters and two physical resistances built around a base 'capacity curve' functional form. Establishing complete databases of RRD parameters for the most important shells systems is a significant challenge in modern shell buckling research.
This paper has presented a computational strategy to aid researchers in designing a programme of parametric FE analyses to establish complete RRD characterisations of metal shell structural configurations and load cases. The research community is particularly urged to embrace the possibilities for mass automation already present in modern FE software to handle the many thousands of fully nonlinear calculations necessary to establish RRD characterisations. A strategy for automated generation, submission, termination and processing of analyses is presented, built around a novel formal analysis hierarchy. Annotated Python code fragments and object-oriented pseudocode are presented to illustrate the concepts, whose logic may be specialised for any shell system and in any programming language.
A crucial capability is the automated termination of an ongoing FE analysis once a 'criterion of failure' has been met. This paper explores the many different failure criteria possible for metal shells depending on the analysis type, and proposes computational procedures for the automated detection of each one.
Fully-functional implementations written in FORTRAN are presented for immediate use as ABAQUS user subroutines, but the logic may be implemented in any FE software with a scripting functionality. Makers of FE software may wish to carefully consider incorporating some of the proposed functionality into their products, as the applications extend far beyond RRD and shell buckling.