8-MW wind turbine tower computational shell buckling benchmark. Part 1: An international âround-robinâ exercise

An assessment of the elastic-plastic buckling limit state for multi-strake wind turbine support towers poses a particular challenge for the modern finite element analyst, who must competently navigate numerous modelling choices related to the tug-of-war between meshing and computational cost, the use of solvers that are robust to highly nonlinear behaviour, the potential for multiple near-simultaneously critical failure locations, the complex issue of imperfection sensitivity and finally the interpretation of the data into a safe and economic design. This paper reports on an international ‘round-robin’ exercise conducted in 2022 aiming to take stock of the computational shell buckling expertise around the world which attracted 29 submissions. Participants were asked to perform analyses of increasing complexity on a standardised benchmark of an 8-MW multi-strake steel wind turbine support tower segment, from a linear elastic stress analysis to a linear bifurcation analysis to a geometrically and materially nonlinear buckling analysis with imperfections. The results are a showcase of the significant shell buckling expertise now available in both industry and academia. This paper is the first of a pair. The second paper presents a detailed reference solution to the benchmark, including an illustration of the Eurocode-compliant calibration of two important imperfection forms.


Introduction
Thin-walled metal shells exhibit among the most complex behaviour of any classical structural form, owing to the interaction of various sources of nonlinearity and a limit state behaviour often controlled by catastrophic buckling. These nonlinearities include plasticity, changes of geometry and the interactions with both geometric and material imperfections, and the reader will find a treasure trove of information relating to these in the 'blue book' of Teng and Rotter [1] and the 'red book' of ECCS EDR5#2 [2] Recommendations. Given the many phenomena which must be considered in the design of a shell structure, where even a simple treatment requires the use of the mathematically two-dimensional 'membrane theory', it is perhaps not surprising that the Eurocode on the strength and stability of metal shells EN 1993-1-6 was the first in world to regulate structural engineering design assisted by finite element analysis. The first draft was published in 1999 as ENV 1993-1-6 [3] and eventually in finalised form as EN 1993-1-6 [4] in 2007. In the years that followed, the Standard received an amendment EN 1993-1-6:A1 in 2017 [5] which introduced updated rules and resistance equations as well as an entirely new design concept named Reference Resistance Design [6,7]. Between 2017 and 2020, a Project Team was active which evolved the Standard further, culminating in a significantly revised prEN 1993-1-6 in 2021 [8] which is due to go out for the stage of 'CEN Enquiry' in March 2023 when it is made available for public use although it is not yet a final publication.
Amongst many innovations, a key legacy of the earliest incarnation of EN 1993-1-6 is a suite of acronyms LA, LBA, MNA, GNA, GMNA and GMNIA whose meaning is explained shortly. These had in fact been devised as early as 1994 during a meeting of the committee led by Prof. J. Michael Rotter (University of Edinburgh) working in close collaboration with Prof. Herbert Schmidt (University of Duisburg-Essen) and Prof. Richard Greiner (TU Graz) on what would later become the draft ENV 1993-1-6 (and which by 1995 had already adopted its future structure of being organised around four distinct ultimate limit states), though the idea for using acronyms as a shorthand for different computational analyses appears to originate from the doctoral thesis of Werner Guggenberger in 1992 [9] at TU Graz and the Reader may find an early incarnation of these in Greiner and Guggenberger [10]. These acronyms have gained significant acceptance in both research and industry as a standardised taxonomy to describe computational analyses of plastic collapse, elastic-plastic and elastic buckling limit states. Their proliferation in journal publications over time is illustrated in Fig. 1, where the first modest usage appears after the launch of the ENV in 1999, followed by a significant increase with the release of the EN in 2007 and again in 2017. The GMNIA acronym, representing the most 'advanced' computational analysis that includes both geometric and material nonlinearities and explicitly considers imperfections of any kind, enjoys particular popularity and appeared in over 40 papers annually since 2019. This reflects the fact that the taxonomy has now found widespread adoption in other applications of steel structural engineering, including a formal incorporation into the upcoming entirely new prEN 1993-1-14 in 2022 [11] on the design of general steel structures assisted by finite element analysis.
The Lead Authors Adam J. Sadowski and Marc Seidel have both taken an active role in the development of EN 1993-1-6 under the leadership of Prof. J. Michael Rotter, with Sadowski having been a UK delegate to CEN/TC250/SC3/WG6 since around 2016 (eventually taking over from Prof. Rotter as its Convenor in 2021) and Seidel having been a German delegate since 2019. Sadowski had additionally been a member of the Project Team CEN/TC250/SC3/PT5, also initially led by Prof. Rotter, tasked with developing the new prEN1993-1-6, to which Seidel had also contributed detailed technical commentary on the various intermediate drafts. After CEN/SC3 had accepted the final draft to go out for international 'CEN Enquiry' at the general meeting in March 2022 at which point it obtained the prEN designation, and given the substantial changes that were introduced into this Eurocode since the 2007 or 2017 versions, the Lead Authors determined that it would be valuable to take stock of the state of computational shell buckling expertise around the world and assess the extent of adoption of some of its key innovations.
To this end, over 40 academics and industry practitioners worldwide were approached with the invitation to participate in a 'roundrobin' exercise. Each participant in this exercise was requested to complete the same exercise 'blindly' without any guidance on the part of the organisers and entirely independently of the other participants (who remained unknown to one another until just before submission of this paper for peer review). Out of 49 colleagues approached by the Lead Authors in April 2022 (identified on the basis of their professional networks), 37 replied (~76% response rate) of which 30 agreed to participate, even if tentatively (~61% acceptance rate). The brief was then circulated and the exercise initialised, with submissions requested to be made in a standardised spreadsheet template with a deadline of early October 2022. In the end, 27 submissions were received (~55% conversion rate) in addition to the two performed independently by the Lead Authors of which one is treated here as a 'reference solution' presented in full in the companion paper. The standard of submissions was consistently very high. Analysts were invited to be co-authors on this study in recognition of the significant and entirely voluntary effort involved and to encourage participation. They are presented in alphabetical order by surname from the third co-author onwards in the list of authors, though a full list of participants and a world map may be found in the Acknowledgements. Readers of this paper should be aware that while its focus is chiefly on the accurate computational prediction of the buckling ultimate limit state and the possible pitfalls that an analyst may encounter along the way, it is also a social Table 1 Geometry of the top segment of an 8 MW wind turbine support tower. overview of the international shell buckling community as of 2022. Other successful 'round-robin' exercises where a technical community was invited to work anonymously on a standardised computational benchmark in structural engineering and which have contributed significantly to the understanding of the state of the relevant disciplines include that of Holst et al. [12,13] on granular media pressures in silos using finite element and discrete element methods, Dyke and Caicedo [14,15] on the seismic response of cablestayed bridges, Hinton et al. [16] on polymer composite failure theories and more recently by Lange and Boström [17] on loaded steel beams under standard fire conditions.

Geometry and other inputs
Participants were provided with the geometry of the top segment of a realistic design for a 36-m high 8-MW multi-strake metal wind turbine support tower as detailed in Table 1 and illustrated in Fig. 2a with geometrical parameters defined in Fig. 2b. Only the midsurface diameters of the top and bottom of the strakes (d 1 and d 2 respectively), thicknesses t and heights h were released. A 'strake' is defined in prEN 1993-1-6 as a zone of constant thickness within a shell constructed with a stepped wall, and for wind turbine support towers it corresponds to what is often called a 'can'. The additional derived data shown for information purposes for each strake in Table 1 includes the meridional angle of incline to the vertical β, the average slenderness ρ avg / t expressed in terms of the average circumferential radius of curvature ρ avg (≈ r = d/2 for near-cylindrical geometries) and approximate dimensionless length parameters ω and Ω about which more will be said later. The participants were not instructed in how to treat the bottom boundary condition other than being told that the bottom flange provided 'restraint'. The material was nominally a S355J0 grade steel, with an elastic modulus of 210 GPa, a Poisson ratio of 0.3, a yield strength of g345 MPa and a density of 7850 kg/m 3 .

Load cases
Participants were asked to perform a series of increasingly sophisticated computational analyses (presented shortly) under two load cases (LCs), with all non-gravity loads stipulated as acting at the centroid through the top of the top flange (Fig. 2c, which was also issued to the participants) although not how they should be specifically modelled: • LC1 without global torsion: A transverse shear force Q = 1.76 MN, a bending moment M = 33 MNm and a vertical downwardsacting force V = 4 MN (originating largely from the dead weight of the rotor and nacelle assembly). Participants were also advised that gravity-induced 'self-weight' of the tower section was present, although not what to do with this information. • LC2 with global torsion: Similar to LC1 with V = 4 MN and self-weight, but with slightly reduced Q = 1.6 MN and M = 30 MNm, and an additional torsional moment T = 22 MNm about the centroidal axis in the direction shown in Fig. 2b.
The values given are design loads according to IEC 61400-1 [18] and no further partial factors needed to be applied.

Requested analyses
The following quasi-static computational analyses were requested for both LCs: • LA (linear elastic stress analysis of the perfect geometry) to compute the mass of the tower, the 6 reactions at the base of the tower and the 3 displacements and 3 rotations at the load application centroid. This is to act as a basic sanity check for the integrity of the finite element modelling and to identify obviously erroneous submissions which would then be removed from further consideration. • LBA (linear bifurcation analysis of the perfect geometry) to compute the 'first ten' linear bifurcation buckling eigenvalues (scaling factors on the applied set of loads) including clear images of their associated eigenmodes and a best estimate of the locations of the peak amplitudes. • MNA (materially nonlinear plastic collapse analysis of the perfect geometry) to compute the evolution of the load proportionality factor (LPF) on the applied loads and the 6 displacement / rotation components of the load application centroid. Participants were asked to identify the critical LPF corresponding to the small-displacement plastic collapse mechanism, and to provide a clear image of this mechanism. • GMNA (geometrically and materially nonlinear analysis of the perfect geometry) to compute similar information as for the MNA but relating to the elastic-plastic buckling mode. • GMNIA, a GMNA with explicitly modelled imperfections, with the sole additional specification that it must represent an "'excellent' (but not perfect)" construction quality.
Participants were not advised to adopt any particular solver settings, material laws or imperfection models. However, they were invited to inform their modelling choices using the provisions of 'EN 1993-1-6 ′ (no specific version was specified) or the Standards in their local jurisdiction. For those familiar with EN 1993-1-6, they were additionally invited to perform a subsequent assessment of the buckling resistance of the tower segment using the LBA-MNA design method [19] and then again using the Reference Resistance Design (RRD) method [6], with the caveat that they report the exact version of the Standard that they used for this purpose.

Requested results
Participants were invited to report the following in a standardised spreadsheet template: • Their affiliation as they would like it to be publicly displayed.
• The extent of their experience in the FE modelling of shell structures and wind turbine support towers specifically.
• The FE software that was used.
• Specific modelling details, such as how they treated the flanges, base boundary condition and load application.
• Motivation for mesh design, including element type, size, uniformity and whether it was changed across analyses or load cases. • Solver settings including incrementation choices and their motivation.

The community
A total of 29 submissions were processed, of which ~ 69% were submitted by colleagues from academia and ~ 31% by colleagues from industry. The imbalance of participants in favour of academia is a reflection of the relatively greater freedom available there to spend time engaging in tasks of this nature. Participants declared between 1 and 25 years' experience in the FE modelling of shell structures (median of ten years), but between 0 and 25 years' experience in the FE modelling of wind turbine support towers specifically (median of zero years). However, the overall quality of the submissions received suggests that FE modelling skills acquired working on one type of shell structure are very transferable to the analysis of other types.

Software choices
Of all submissions, ~66% were completed the exercise using the ABAQUS finite element software (with versions used including 6.14-1, 6.17 and 2017 to 2022), while the remaining ~ 34% used the ANSYS package (either Workbench or APDL, with versions used including 19.4, 2021 R1/R2 and 2022). The varying different versions of each software are unlikely to have influenced the results in any meaningful way, given that the finite element formulations and buckling-relevant solvers have been stable for many years and new versions offer either only largely cosmetic changes or focus on specialised capabilities such as XFEM that were not requested here. Academic participants exhibited a distinct preference for ABAQUS with ~ 80% opting for it out of this group, whilst ANSYS was broadly favoured by participants from industry with ~ 67% using it out of this group. Both software packages offer a powerful scripting functionality to automate the generation of complex models, and have a similar capability in terms of performing linear, eigenvalue and nonlinear path-tracing analyses.

Modelling of the tower geometry
Most submissions exploited the property that the tower segment is a shell of revolution to efficiently generate the midsurface by the ABAQUS or ANSYS pre-processor, with at least two submissions choosing to generate the model input file directly using their own custom-written mesh generator. The majority of submissions (~93%) chose to create a 'reference point' (RF) on the axis of revolution to simplify the introduction of centroidal point loads and the specification of boundary conditions. A range of valid approaches was used to treat the top Flange 101, such that almost all submissions had a top boundary that was either rigidly circular or had a sufficiently high finite stiffness to act as such. Maintaining the circularity of end boundaries is particularly important to restrict out-ofround distortions in shells of revolution subject to unsymmetrical loads [20,21]. It should be added that the approximate dimensionless length Ω ≈ 1.09 of present tower segment classifies it as being either inside (if to consider uniform bending conditions [22]) or on the cusp of being inside (if to consider cantilever bending conditions [23]) a 'transitional' length domain where bending-induced ovalisation begins to have a deleterious effect on the buckling resistance of a cylinder with circularity restrained at end boundaries. If circularity is not maintained, the tower may behave as a 'long' cylinder undergoing fully-developed pre-buckling ovalisation [21] under bending loads which would lead to a conservative GMN(I)A resistance.
A substantial portion of submissions (~60%) modelled the top flange as a thicker shell strake. The point loads were then applied to a centroidal RF that was connected to the top-most edge of shell nodes via a rigid body kinematic coupling that also served to maintain a rigid circular cross-section at this location. At least two submissions did not create a RF and applied the loads in a distributed manner along the top-most circumferential shell edge, and one submission modelled a transverse circular plate at the top edge of the same thickness as the flange (50 mm) to maintain circularity. Four submissions (of which one is the reference solution) modelled the top flange as a beam element representing a circular hollow section whose node was connected to the top edge of Strake 102 via a rigid body kinematic coupling. The point loads were then applied directly to the top-most node of this beam element. Since classical beam elements cannot normally undergo cross-sectional distortion, circularity was maintained. Lastly, at least five submissions did not physically model the top flange at all, representing its effect instead with a direct rigid body coupling between the top edge of Strake 102 and the centroidal RF.
The bottom boundary represents a more straightforward modelling detail. Some submissions (~39%) modelled the bottom Flange 116 as a shell strake linked kinematically to a centroidal RF whose 6 dofs were constrained to zero. Approximately the same number of submissions (~39%) did not use a RF at the base and instead applied boundary conditions directly to the nodes of the bottom edge of Flange 116, with choices varying between a fully restrained BC1r boundary or a pinned BC1f boundary. Other submissions modelled Flange 116 as a beam while the reference solution ignores it entirely and instead creates a fully-constrained RF at the base of Strake 115. Each one of these modelling options is entirely valid provided that the flange is stiff enough to prevent ovalisation (if it is not, the entire tower would need to be modelled rather than a single segment), and unlikely to be a source of any significant differences in the results.

Mesh design
Although it was not always clear, most submissions, including the reference solution, likely opted for a uniform mesh with no local refinement. This was variously motivated as having been done either for the sake of modelling speed and convenience, or because there was an awareness of the possibility of local buckling being critical within any of the strakes. At least two submissions used a 'doublebiased' graded mesh within each strake that was finer near the junctions and coarser in the 'membrane' region, a fair approach if the aim is to accurately capture the localised regions of compatibility bending near strake junctions but there is then a risk that a critical local buckling mode which would have occurred in a now coarsely-meshed 'membrane' regions may be missed. Most submissions did not vary the mesh from one analysis to another, again partly out of convenience but also because the eigenmode shapes of an LBA were often used as the source imperfection for the GMNIA for which consistency of the mesh across both analyses is vital. The majority of submissions (~75%) opted for a 4-node bilinear thin/thick shell element (S4/S4R in ABAQUS or SHELL181 in ANSYS), while the remaining ~ 25% opted for an 8-node biquadratic thin/thick shell element (S8R in ABAQUS or SHELL281). At least one submission used different meshes for both LCs.
Submissions invariably adopted one of two distinct strategies to establish an 'optimal' mesh resolution. The first is the classical approach of progressive refinement based on feedback from some convergence metric, taken variously as the highest LA or GMNIA displacement or von Mises stress or the critical LBA eigenvalue. Usually, a 'fine' mesh was first adopted and then made progressively coarser until the corresponding prediction for the metric exceeds that of the 'finest' mesh by some margin, taken variously in the range 1-5%. While this is a good approach in principle, unless special care is taken the analyst cannot know how 'fine' their 'finest' mesh really is and the possibility that such a mesh may not adequately capture the smallest-scale local meridional compression buckles (which control the buckling limit state and which can be very small relative to the dimensions of the shell [20,24];) cannot be discounted. The average 4-node element sizes for meshes optimised in this manner were ~ 119 mm and ~ 129 mm in the circumferential and meridional directions respectively (calculated using averages of element size ranges reported for ~ 32% of submissions). For 8node elements the corresponding size was ~ 193 mm in both directions (just two submissions).
In the second approach, analysts did not do any progressive refinement but instead established element sizes relative to some characteristic measure governing the scale of the most localised feature that the mesh must capture. This measure was typically taken as the linear bending half-wavelength of a cylindrical shell given by λ ≈ 2.444√(rt), which offers a very close approximation to the actual expression for a conical shell when it is near-cylindrical. Meridional element sizes were variously reported as being in the range λ/30 to λ/5, with similar circumferential sizes adopted in to maintain approximately square elements. For meshes designed in this manner, the overall average mean 4-node element sizes were ~ 91 mm and ~ 83 mm in the circumferential and meridional directions respectively (~36% of submissions), or ~ 115 mm and ~ 86 mm for 8-node elements respectively (~14% of submissions), on average finer than those achieved using the classical progressive refinement approach. In the remaining submissions the mesh design methodology was not explicitly stated. A visual comparison of the two mesh protocols is given in Fig. 3. This second approach is recommended as it guarantees that the mesh is appropriately designed for the scale of the smallest possible buckling feature that it must accurately represent and dispenses with the need for a laborious progressive mesh refinement study. It is the approach adopted by the reference solution (black cross in Fig. 3), and further information may be found in the companion paper and in Sadowski [25].

LA calculations
The combined results of all LA stress analysis calculations (tower mass, base reactions and tip dofs) are presented in Figs. 4 and 5, the aim being to filter out any submissions that may be grossly in error from further consideration. A number of submissions exhibited very obvious order-of-magnitude errors for individual results which were altered directly, but this was the only direct intervention made during processing and all subsequent results were taken at face value exactly as reported. All values presented below are absolute values to account for varying coordinate system conventions used.
The figures show 'box-and-whisker' plots of the distributions of the predictions for each entity according to the load case (LC). The blue 'boxes' are in places not easy to discern because they represent the 25th and 75th percentiles of the distribution which are in fact extremely close, a positive outcome that signifies that the majority of submissions computed the correct reactions and tip dofs and there is very little variability around the median (the red lines). The black crosses represent the results of the reference solution and are consistently in agreement with the predictions of the majority of submissions. However, the lower and upper 'whiskers' represent the 0th and 100th percentiles respectively, and they are on occasion very far away from the 'boxes' because of outliers (black red-filled circles) which appear to represent either calculation or reporting errors, and in at least one case a misinterpretation of the LCs due to an ambiguity in the way bending moments were presented in the original brief.
Any submission for which the absolute value of a non-negligible reaction differed by more than 25% from the overall median was flagged as an erroneous outlier and that submission was not included in subsequent analyses (these same submissions also happened to correspond to most of the outliers in the tip deflections seen in Fig. 5, with the remaining ones being attributable to obviously misreported individual entries). This was done so that the variability in the remaining results could be plausibly attributed to modelling choices only rather than to gross errors. Unfortunately, 7 out of 29 (~24%) submissions were removed from further consideration in this way, illustrating the vital importance of correctly performing and documenting an LA as a stepping-stone to the correct execution of more sophisticated analysis where the scope for error is even greater. It is, of course, not guaranteed that the remaining 23 submissions are free from either calculation or reporting error. In what follows, there may be fewer than 23 reported calculations because not all submissions presented a complete set of results.

LBA calculations
The combined results for the 22 valid remaining LBA calculations, a linear perturbation procedure, are presented here. The majority of the remaining submissions (~68%), including the reference solution, chose to scale all loads in a single eigenvalue analysis step for either load case (LC), including self-weight. While this has the physically unrealistic consequence of scaling the effect of gravity beyond 100% intensity, the resulting critical LBA eigenvalues (R cr ) are then slightly on the conservative side and the practice is codecompliant in that self-weight is then treated as a regular design load. Of the remaining submissions, six applied gravity to 100% intensity in a pre-loading step prior to the perturbation procedure in which only the tip loads were scaled, while at least one applied both gravity and the V load in a pre-load step and subsequently scaled only Q, M and T on the basis that the V load represents the gravityinduced weight force of the nacelle and similarly should not scale beyond 100% intensity.
The different predictions for the 'first ten' LBA eigenvalues for both LCs are illustrated with box-and-whisker plots in Fig. 6. For LC1, the majority of predictions are in the vicinity of the overall median for each eigenvalue, suggesting a consistent set of predictions. Those that are excessively high are for a submission with a particularly coarse mesh which gives a very stiff response, highlighting the importance of an adequate discretisation. One submission reported negative eigenvalues for the first five eigenmodes, about which more is said shortly. The eigenvalues all increase monotonically with the eigenmode number, representing increasingly higher-order buckling modes with increasingly localised features in the region of Strakes 112-114 (Fig. 7, top).
For LC2, the uniform torsional loading appears to dominate the pre-buckling stress state (in that the structure is energetically closer to bifurcating into a torsional mode than any other) which causes some non-uniqueness regarding the positioning of the corresponding torsional eigenmode (Fig. 7, bottom). For this reason, the eigenvalues come in groups (here, pairs) and represent very similar eigenmodes whose positions are merely shifted around the circumference. For the same reason, the direction of the dominant torsional moment does not matter and some solvers reported eigenvalues that were plausible in magnitude but negative in sign which contributes to the apparently severe observed 'scatter'. The plausibility of the magnitudes of the negative eigenvalues is entirely coincidental and potentially misleading, and should not be interpreted as meaning that the load set in LC2 is somehow sign-reversible. It is strongly recommended to mitigate against the eigenvalue solver reporting negative eigenvalues as these unrealistically represent a load set (including gravity) acting in the opposite direction to that intended, for example by specifying a minimum eigenvalue of zero, and negative eigenvalues should never be retained. Of those submissions that reported solver settings, the majority (~87%) used their software's implementation of the powerful Lanczos algorithm.

MNA calculations
The combined results for the 21 valid MNA calculations, a path-tracing procedure to compute the small-displacement reference plastic collapse load R pl and determine its corresponding plastic collapse mechanism are presented here. The majority of submissions (~57%, fewer than for the LBAs) scaled all of the loads within a load case to failure (including self-weight), but at least seven and two submissions respectively chose to apply gravity or gravity with the V load to 100% intensity in a separate pre-load step and not scale these to failure hereafter.
The different predictions for the equilibrium relationships between the load proportionality factor (LPF) on the scaled loads and the Q-direction displacement for LC1 and LC2 are presented in Figs. 8 and 9 respectively. These are representative of the scatter observed across all submissions for other tip displacements and rotations, which are not shown, and illustrate the very wide range of approaches taken to perform an MNA. A selection of corresponding images submitted in response to a request for the 'critical failure mode' is shown in Fig. 10, where this was variously represented by means of contour plots of yield flags, von Mises stresses and strains, or equivalent plastic strains. The majority correctly identify the junction between Strakes 111 and 112 as the location of extensive plastic deformation for LC1 and additionally the base of Strake 103 for LC2.
At least a part of the scatter in the critical LPF may be explained by the varying material laws employed, all of which assumed a bilinear stress-strain curve (Fig. 11) albeit with varying choices for the post-yield strain hardening modulus. The majority of submissions assumed either zero strain hardening (with the portion of the curve being either exactly horizontal or with a likely negligible gradient of the order of just 10 1 MPa, representing the '*' in Fig. 11), with several assuming E/100 = 2100 MPa and just a handful assuming E/1000 = 210 MPa. It should be noted that in strict code-compliance with prEN 1993-1-6 [8] an MNA should be performed with an ideal elastic-rigid plastic law, although a strain hardening modulus up to E/1000 is permitted in case of numerical convergence issues and prEN 1993-1-14 [11] recommends an even more restrictive E/10000 = 21 MPa. Another source of scatter in the equilibrium curves (also in those of GMN(I)A, but not in any critical LPF) may have been an unclear specification in the brief, with some submissions potentially reporting the max. obtained displacement components anywhere in the structure as opposed to those at the tip load application point. Of those submissions that reported solver settings, approximately half employed their software's implementation of a load-controlled Newton-Raphson static solver and half employed a load-controlled arc-length solver. As MNAs cannot undergo a loss of geometric stiffness and thus cannot exhibit a reduction in the load proportionality factor, either of these solver choices is equally adequate.
The actual reference R pl was usually, but not always, reported as the peak LPF of the portion of the equilibrium curve that was computed (see first boxplot in Figs. 8 and 9). Where it was not (second boxplot in both figures), it was either estimated as a future projection based on the tendencies of the computed portions of the equilibrium curve using a Modified Southwell Plot (two submissions) or its enhancement the Convergence Indicator Plot [26] (also two submissions), or it was variously taken as the load factor at which the max. equivalent plastic strain exceeded 5-20% or at which the von Mises stress exceeded a certain fraction of the yield stress. The former approach results in an R pl that is slightly higher than the peak LPF while the latter leads to R pl that are lower, although overall the distribution of values reported as the actual R pl is slightly lower than for the peak LPF (compare the second and first boxplots in both figures). In the Lead Authors' view, the concern about reporting a value of R pl so as to avoid an excessive computed strain or stress somewhere in the structure is unnecessary because the reference plastic collapse load anyway represents an idealised limit state. It may also be unconservative because it leads to a lower estimate of the dimensionless slenderness λ = √(R pl /R cr ) which, depending on how it is used, may lead to a dimensionless resistance χ(λ) that is higher than it should be. Equally, it is undesirable and potentially unsafe to overestimate R pl by specifying an excessively high strain hardening modulus which for an MNA should ideally be exactly zero.

GMNA calculations
The combined results of the 22 valid GMNA calculations are presented in the form of computed equilibrium curves of the LPF Fig. 8. Equilibrium plot of the computed LPF against the Q-direction tip displacement, with box-and-whisker plots of the peak LPF from these curves and the reported critical R pl for LC1. The black line and crosses represent the reference solution. Fig. 9. Equilibrium plot of the computed LPF against the Q-direction tip displacement, with box-and-whisker plots of the peak LPF from these curves and the reported critical R pl for LC2. The black line and crosses represent the reference solution.
against Q-direction displacement (with other tip displacement or rotation components again displaying a similar degree of scatter) in Fig. 12 for both LC1 & 2. The choice of material laws was slightly different to the MNAs (compare Figs. 11 & 13), with more submissions now choosing to employ either measured or codified complete engineering stress-strain curves (perhaps unnecessarily, given that design calculations should arguably be limited to nominal inputs). Of those submissions that reported solver settings, the usage was approximately evenly split between a load-controlled arc-length solver and a load-controlled Newton-Raphson solver (with one submission using an implicit dynamic solver). The arc-length solver is known to be capable of following a softening equilibrium path by being able to reduce the controlling LPF (see discussion in Section 6.3 of Sadowski et al. [7]). However, only a handful of submissions using an arc-length solver included any significant portion of a post-buckling path (Fig. 12) and most showed distinct signs of convergence problems in the vicinity of the first buckling event. The difficulty in predicting the post-buckling response appears to be an ongoing and significant limitation of modern finite element software and deprives the analyst of important insights into the complete nonlinear behaviour of their structural system. Lastly, while there is some scatter in the reported critical characteristic resistances R GMNA and still some worrying outliers (some analysts reported major convergence issues when using kinematic coupling to a RF at the top boundary which they attempted to overcome through various means), most submissions are relatively consistent in their predictions for both load cases.

GMNIA calculations
The combined results of the 22 valid GMNIA calculations are explored here in detail. The somewhat cryptic brief specification that the tower model should represent an "'excellent' (but not perfect) construction quality" was a hint to the participant to consider designing the amplitude of the adopted imperfection to be compliant with a Fabrication Tolerance Quality Class (FTQC) A of EN 1993-1-6. Equally, for participants who do not use EN 1993-1-6 in their jurisdiction this was a hint to consider a correspondingly lenient amplitude according to their Standard of choice if it permits differentiation of the design resistance according to anticipated quality of construction. As no hint was given as to the expected form of the imperfection, the variability in the choices made across all submissions, as was as the resulting scatter in the predicted GMNIA resistances, was substantial (compare the 'box' for GMNIAs in Fig. 14 with those for the other analyses). If to express this scatter through a dispersion coefficient calculated as the interquartile range (distance between the 25th and 75th percentiles) normalised by the median, the computed GMNIA resistances exhibit between ~ 30% and ~ 21% dispersion for LC1 and LC2 respectively compared to ~ 2-4% for the LBA and GMNA resistances (but ~ 13-14% for MNA resistances due to their sensitivity to the adopted material model).
Across all submissions, there is little consensus as to whether the tower design is actually safe on the basis of the GMNIA: ~42% and ~ 26% of submissions have (unfactored) characteristic R GMNIA resistances that are less than unity for LC1 and LC2 respectively (both counts would increase if a partial factor γ M1 were included). This is apparent from the complete set of GMNIAs for both LCs shown as empirical cumulative distributions in Fig. 15 which illustrate the frequency of occurrence of GMNIA resistances up to and Fig. 11. Plot of adopted material laws in terms of engineering stress-strain curves for MNA. The black line represents the reference solution. Some submissions reported this as the software input starting at a zero plastic strain rather than a total strain. including a particular value, and which participants may use to identify their own submission and placement relative to the others.
Most submissions (~64%) employed only LBA eigenmodes (Fig. 7) as the equivalent geometric imperfection for both load cases, shown on the left of Figs. 16 and 17 in the form of equilibrium curves of the LPF against the Q-direction displacement as well as boxand-whisker plots of distributions of the critical resistances. The choice of actual eigenmode settings was variously argued as having been made on the basis of their perceived realism, modelling convenience (both ABAQUS and ANSYS offer a particularly convenient mechanism for importing LBA results as mesh perturbations into subsequent analyses) or conservatism based on preliminary test calculations. For example, the LC1 critical eigenmode (Fig. 7 top) may be argued to represent a local (relative to the total dimensions of the tower) accidental dimple and thus a plausible imperfection for a GMNIA of LC1. However, a similar argument for the LC2 critical eigenmode (Fig. 7 bottom) is harder to justify since it does not resemble any likely accidental or systematic defect, although this was the imperfection adopted by most GMNIA submissions for LC2 as it is quite damaging. At least two submissions did not use the LBA eigenmode corresponding to the smallest (positive) eigenvalue as an imperfection, but instead identified a 'most critical LBA eigenmode' out of a selection of different eigenmodes trialled in different GMNIAs (a reminder that the LBA eigenmode corresponding to the smallest eigenvalue may not necessarily lead to the most conservative equivalent imperfection). Other than modelling skill, an important contribution to the high scatter observed is likely due to the varying ways that the amplitude of the imperfection was controlled. Submissions usually controlled the scaling of the eigenmode so as to achieve a target imperfection amplitude as specified by the dimple tolerance requirements for FTQCA in EN 1993-1-6 but variously using the 2007, 2017 or 2021 versions of the Standard (the Fig. 13. Plot of adopted material laws in terms of engineering stress-strain curves for GMN(I)A. The black line represents the reference solution. Some submissions reported this as the software input starting at a zero plastic strain rather than a total strain. provisions have changed over the years), and in one case using the out-of-roundness tolerance requirements. Most critical load factors were correctly reported as being the tip of the computed equilibrium curve, but some submissions conservatively reported this as the LPF at which yielding was first detected.
A smaller portion of submissions (~23%) employed the algebraic weld depression formulation of Rotter and Teng [27] to create an axisymmetric weld imperfection at every strake junction, controlling its amplitude to achieve a FTQC A according to the EN 1993-1-6 dimple tolerance requirement (2017 and 2021 versions only). These are shown in the middle of Figs. 16 and 17. Submissions using the weld depression exhibited a smaller dispersion coefficient than those using eigenmode (~9% vs ~ 30% respectively for LC1 and ~ 2% vs ~ 21% respectively for LC2). While this partly reflects the smaller number of submissions using weld depression imperfections (their generation is more involved and cannot realistically be done without scripting), it is primarily a reflection of the significantly reduced ambiguity in the resulting geometry since the Rotter and Teng formulation follows a mathematical formula whereas an LBA eigenmode depends on the modelling choices under which that LBA was performed. While the median prediction for GMNIA with weld depression imperfections is less conservative than for GMNIA with LBA eigenmode imperfections, it is cautioned that this should not be generalised as a universal tendency. Ideally, a design calculation using GMNIA should trial several candidate imperfections at codecompliant amplitudes corresponding to the given FTQC to establish the likely imperfection sensitivity of the system. More is said on this in the companion paper. A selection of submitted screenshots illustrating the computed buckling mode is shown in Fig. 18 for both LCs and any assumed imperfections. These were variously reported as being of different stress states (typically von Mises) at  differently scaled deformed states, peak equivalent plastic stress states or incremental displacements between the pre-and postbuckling states. The variability in these is a direct reflection of the variability in the modelling choices made.
Lastly, the remaining submissions adopted unique imperfection methodologies that cannot readily be generalised, although as these were not assessed as outliers on the basis of the LA results their predictions have been included in the 'all imperfections' distributions shown in Figs. 14 to 17. These broadly included imperfections consisting of superpositions of multiple LBA eigenmodes at an amplitude not governed by a Standard, combinations of weld depressions and LBA eigenmodes within a single model, the use of the Energy Barrier Criterion [28], the use of a cantilever-like global imperfection created by analysing the tower under a single tip load and importing it as mesh perturbation, and finally the representation of a single purposefully-placed weld depression by a cosine shape although at an amplitude compliant with EN 1993-1-6 dimple tolerances.

Optional LBA-MNA and RRD assessments
As an optional component of this 'round-robin' exercise, participants were additionally invited to perform LBA-MNA and RRD assessments of the tower design under both load cases. LBA-MNA is a semi-computational design method present since the earliest ENV 1993-1-6 in 1999 [3] that permits an analyst to use the LBA and MNA predictions for R cr and R pl respectively to estimate the dimensionless slenderness as λ = √(R pl /R cr ) and then to use an algebraic formulation of a 'capacity curve' to obtain a knockdown factor χ(λ) on the reference plastic collapse load R pl so as to estimate a characteristic resistance R k = χ(λ)R pl [2,19] which is then used together with a partial factor γ M1 to obtain a design resistance R d = R k /γ M1 . The Reference Resistance Design (RRD) method [6] introduced in EN 1993-1-6:A1 in 2017 [5] took this one step further to create a resistance check that is entirely algebraic for the end user. Both the LBA-MNA and RRD provisions have undergone significant re-development in the new prEN 1993-1-6 [8], as have the source equations for the capacity curves which reside in Annexes D and E therein.
It has always been challenging to interpret these provisions when applied to more complex multi-strake or multi-segment shells, in particular regarding the treatment of situations where critical eigenmode from an LBA and the plastic collapse mechanism from an MNA occur in very different locations and / or under very different controlling stress conditions (compare Figs. 7 and 10), and also in terms of the appropriate choice of capacity curve with which to translate the slenderness λ into a knockdown factor χ. The collected results of 14 (~63%) and 9 (~41%) submissions which carried out LBA-MNA and RRD assessments respectively are summarised in Tables 2 and 3 without further commentary so that the variability in interpretation and outcomes can speak for itself. There is clearly little consensus as to most appropriate way to apply the provisions nor whether the tower has safely passed either design check. Indeed, most attempts concluded that the tower is not safe for LC1 or LC2. Any recommendations for the appropriate interpretation of LBA-MNA and RRD are outside of the scope of this report, but clearly much development remains to be done to complete the coverage of load cases and geometries in the form of capacity curves, investigate the possibility of reducing the methods' conservatism and to illustrate their usage on complex multi-strake civil engineering shells.

Conclusions
Some invaluable conclusions can be drawn from this 'round-robin' exercise on computational shell buckling of a standardised metal wind turbine support tower design presented in this paper. These are a fair reflection of the state of the relevant expertise across academia and industry internationally at the end of 2022. In particular: • There were more submissions from academics than industry practitioners by a ratio of approximately 2.2 to 1. The software used were exclusively ABAQUS and ANSYS by a ratio of 1.9 to 1 and were more likely to have been used by academics and industry practitioners respectively, reflecting the traditions of both communities. Both software packages exhibited comparable capabilities as far as the analyses requested by the exercise presented here are concerned. • A finite element mesh purposefully designed to accurately represent the scale of the smallest possible buckling-relevant feature will on average lead to finer meshes than one obtained by a classical iterative mesh convergence study and with much less modelling effort. • A consensus should be established regarding the appropriate way of treating gravity loads, since the code-compliant treatment of scaling them to arbitrary intensity like any other design load is physically questionable. • Not performing a carefully-documented LA as a stepping-stone to more complex analyses is a false economy, as it is clearly a valuable sanity check of model correctness. For the determinate tower structure investigated here, reactions can even be verified manually. For the same reason LBA and MNA should always be performed prior to a GMNA, and a GMNIA should always be performed last. • LBAs should be performed with a minimum requested eigenvalue set at zero to avoid eigenmodes with negative eigenvalues corresponding to direction-reversed load cases. All negative eigenvalues should be discarded. Reverse-sign non-gravity loads should be modelled as a separate load case if these can occur. • MNAs should not be terminated early on the grounds that an equivalent plastic strain has been greatly exceeded somewhere in the model since the analysis is anyway aiming to compute an idealised reference plastic collapse mechanism. A low MNA result may lead to an unconservatively low dimensionless slenderness. • MNAs should not be executed with a tangent modulus larger than E/10000 (ideally a zero tangent modulus should be used) as this may lead to the converse situation where the yield limit is significantly exceeded. • The scatter exhibited by the collective GMNIA results was found to greatly exceed the scatter observed in the collective LBA, MNA or GMNA results. This large scatter is due to the very many different assumptions made regarding imperfection forms, calibration of amplitudes and varying modelling skill.
It is hoped that the present study will give the international community working with metal shell structures food for thought and, together with the reference solution offered in the companion paper, will inform their professional practices going forward. The significant scatter observed across all calculations, but GMNIA in particular, suggests that some sort of standardisation of modelling protocols is very desirable to enable reproducibility of calculations and consistency of assumptions and practices across different teams of analysts. Similarly, there is currently no consensus on the application of the LBA-MNA or RRD design concepts which have  consistently led to the conclusion that the benchmark tower is not safe for either load case, although particularly for the load case involving torsion. This, too, is food for thought for the developers of the Eurocodes.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
The data that has been used is confidential.