Validation of Stainless-Steel CHS Columns Finite Element Models

Stainless-steel elements are increasingly used in a wide range of load-bearing structures due to their strength, minimal maintenance requirements, and aesthetic appearance. Their response differs from standard steels; therefore, it is necessary to choose a different procedure when creating a correct computational model. Seven groups of numerical models differing in the used formulation of elements integration, mesh density localization, nonlinear material model, and initial geometric imperfection were calibrated. The results of these advanced simulations were validated with published results obtained by an extensive experimental approach on circular hollow sections columns. With regard to the different slenderness of the cross-sections, the influence of the initial imperfection in the form of global and local loss of stability on the response was studied. Responses of all models were validated by comparing the averaged normalized ultimate loads and the averaged normalized deflections with experimentally obtained results.


Introduction
In comparison with the most commonly used carbon steel material, stainless steel is more recent in the field of structural engineering. For example, the UK's first stainlesssteel road bridge over the River Eamont in the village named Pooley Bridge was officially opened in October 2020. The first composite stainless-steel vehicular bridge in Europe (most probably also worldwide) was built in Cala Galdana on the island of Menorca, Spain, and was opened in June 2005 [1]. This duplex stainless-steel arch bridge over Algendar Creek replaced the existing reinforced concrete bridge which suffered severe corrosion issues due to the marine environment. Despite the higher cost of the stainless-steel material, the usage of this material is most beneficial in cases of structures with heavy traffic volumes or those exposed to such aggressive environment as was discussed for the case of the steel ASTM A1010 in a recent cost-efficiency study by Daghas et al. [2]. Usage possibilities of the stainless-steel material in combination with different materials have been recently studied by Pauletta et al. [3], who investigated the bond-slip behavior between stainlesssteel reinforcement bars and concrete, or by Corradi et al. [4], who reviewed the use of stainless-steel profiles to reinforce or repair historical wooden structures. Research is also conducted also on the enhancement of the stainless-steel corrosion resistance by Dinu et al. [5], who proposed a certain improvement of the grade AISI 304 stainless steel.
In addition to higher cost and strong corrosion resistance, another significant difference, compared to ordinary carbon steel, from the aspect of mechanical properties is the nonlinear stress-strain behavior. The stress-strain curve of a stainless-steel material possesses no pronounced yield plateau and exhibits an early declination from the linear elastic behavior with significant strain hardening. Along with higher strength and higher ductility, these are the most significant differences compared to carbon steel [6]. Various analytical equations describing this material behavior are found in the literature, e.g., in the study by Hradil et al. [7]. A constitutive model of stainless steel in a high-strain range was investigated by Peng et al. [8]. The stress-strain curves were also compared by Real et al. [9]. These descriptions are based on the expressions originally proposed by Ramberg and Osgood [10], subsequently modified by Hill [11]. The use of 0.2% proof stress, as the equivalent of the yield stress, has become a standard industry practice [12]. The variability of material parameter values is significantly high among different stainless-steel grades and different groups (types according to metallurgical structure) [13]. The most commonly used types of stainless steel are austenitic, ferritic, and duplex ones.
Material parameters can be feasibly determined through the process of parameter identification, where either the analytical stress-strain curve or the stress-strain curve as the result yielded from the finite element numerical analysis is fitted to be in the best possible match with the experimentally measured data. Such material model calibration has been conducted by Jindra et al. [14]. Values of the material parameters could also be obtained, based on the informative Annex C of the European standard EN 1993-1-4 [15], from the tables in AS/NZS 4673 [16] or SEI/ASCE-8 [17]. Recently, a statistical study of not only stainless-steel material parameter values, based on the results consolidated over the last decades, has been conducted by Arrayago et al. [18].
The design guidance to determine the flexural buckling capacity of stainless-steel CHS (circular hollow section) members is provided in EN 1993-1-4 [15]. This approach is consistent with the same task for carbon steel structural elements described in EN 1993-1-1 [19]. Many existing data points (based on either physical experiments or properly validated geometrically and materially nonlinear finite element analyses with imperfections -GMNIA) lie below the current EN 1993-1-4 [15] flexural buckling curve, as described by Young et al. [20], Rasmussen et al. [21], Ashraf et al. [22], Theofanous et al. [23], or Shu et al. [24]. Originally, due to the insufficient amount of stainless-steel CHS experimental data at the EN standard creation time, the buckling curve was calibrated mainly by considering cold-formed rectangular hollow section (RHS) and square hollow section (SHS) results of the column buckling tests [21]. The buckling curve based on cold-formed RHS and SHS experimental results may not be appropriate for CHS due to increased material strength in hardened corner areas of RHS and SHS cross-sections [25].
Physical experiments on stainless-steel CHS elements exposed to compressive loading have been conducted by Rasmussen et al. [26], who investigated austenitic stainlesssteel stub columns. Authors of similar research include Burgan et al. [27], Talja [28], Rasmussen [29], Young et al. [20], Uy et al. [30], Zhao et al. [31], Lam et al. [32], Gardner et al. [33], Bardi et al. [34], and Paquette et al. [35]. A wide range of loading eccentricities in the compression of CHS columns was investigated by Buchanan et al. [36]. It can be noted that the present article studies the static resistance of the column as a supporting structural member, and the heat and flow load cases [37] inside the tube are not considered.
The benchmark data presented in this study are obtained from a comprehensive experimental program conducted at Imperial College London (ICL) and Universitat Politècnica de Catalunya (UPC) and a previous research conducted by Buchanan et al. [38]. Tests of the two austenitic (A) cross-sections, 104 × 2 CHS and 106 × 3 CHS, and the duplex (D) cross-section 88.9 × 2.6 CHS were undertaken at ICL, and those of the two ferritic (F) cross-sections, 101.6 × 1.5 CHS and 80 × 1.5 CHS, were carried out at UPC [38].
To conduct further geometrically and materially nonlinear analyses of CHS stainlesssteel columns in compression, the parametrical numerical finite element models created in ANSYS Classic(v 19.0.) utilizing APDL (ANSYS parametric design language) [39] need to be properly validated.
This contribution aims at the validation of the numerical FE models of CHS columns buckling tests and comparison of the validation results with the research of Mr. Buchanan, Mrs. Esther Real, and Mr. Leroy Gardner [38], where another FE software, Abaqus/CAE 2016 [40], was used, and a slightly different mesh size [38] was used. Moreover, this study compares the results of various modeling approaches, where either different mesh sizes, different elements, or element formulations are considered. This approach comparison is based on several selected model cases.

Experimental Testing Program
The flexural stability tests consisted of 47 concentrically loaded pin-ended specimens with a wide range of local and global slenderness. Sample lengths were chosen to provide a range of global slenderness values λ.
A detailed description of the comprehensive experimental program is well documented and described in the study of Mr. Buchanan et al. [38] pp. 298-303.

Numerical Finite Element Model
A parametrical numerical finite element model has been created in the ANSYS Classic technology (v 19.0.) [39] using APDL macros. The variable input parameters were all the material parameters (Tables 1 and 2), geometrical parameters (Tables 3-7), as well as the initial global imperfection amplitude (ω 0 + e 0 ).
Two sets of material property values, referred to as SCP and TP, have been considered in this study. Both sets are based on the averaged values from all of the relevant available data: for the stub-column material properties (SCP, Table 1, based on the data from Table 4 in [38]) and for the tensile coupon material properties (TP, Table 2, based on the data from the Table 2 in [38]), respectively. E 0 is the material elastic Young's modulus, σ 0.2 is the material 0.2% proof stress, and n is a strain-hardening exponent, σ 1.0 is the 1% proof stress of the material, n' 0.2,1.0 is a strain-hardening exponent, and σ u is the ultimate tensile stress. Stress-strain relations defined by these parameters are elaborately described in Section 3.3.
Geometrical parameters, D (cross-section outer diameter), L (effective structural length, with the inclusion of the additional knife edge lengths), t (wall thickness), and the imperfection amplitudes (ω 0 + e 0 ) have been considered in accordance with the measured values of the respective specimens provided (Tables 3-7 below, based on the data from the Tables 5-9 in [38], respectively).  Table 3. Measured and calculated geometric properties of the cross-section set 106 × 3 CHS (circular hollow section) pin-ended columns (austenitic steel).

Specimen [38]
D (mm) t (mm) L (mm) L/(ω 0 + e 0 ) (-) ω 0 + e 0 (mm) D/t (-) A (mm 2 ) I (mm 4 ) All the 74 finite element analyses (which differ in material and geometrical parameters) to validate the numerical model results have been conducted only considering only the modeling approach #A. The other model variants (#B-#F) have been conducted only for chosen geometries, and the results are documented and compared as well.

Modeling Approach #A
Four-node structural shell elements (SHELL 181) with 3 translational degrees of freedom (DOF) and 3 rotational DOF per node have been used to model the CHS columns. The elements possess bending and membrane stiffness (Mindlin-Reissner theory). A reduced integration with hourglass control has been considered with one integration point (three through the thickness). The elements include the linear effect of transverse shear deformation. An assumed shear strain formulation of Bathe-Dvorkin is used to alleviate shear locking [39]. The geometrical shapes of all the shell elements used to model the CHS column are rectangles, with a maximal edge size of 8 mm in the longitudinal direction and 5 mm in the tangential direction (along the circumference) ( Figure 1a). This differs from the mesh size adopted by Buchanan [38], where the longitudinal and circumferential dimensions of the mesh are considered the same as the wall thickness of the CHS tube, t.
deformation. An assumed shear strain formulation of Bathe-Dvorkin is used to alleviate shear locking [39]. The geometrical shapes of all the shell elements used to model the CHS column are rectangles, with a maximal edge size of 8 mm in the longitudinal direction and 5 mm in the tangential direction (along the circumference) ( Figure 1a). This differs from the mesh size adopted by Buchanan [38], where the longitudinal and circumferential dimensions of the mesh are considered the same as the wall thickness of the CHS tube, t.

Modeling Approach #B
In this study, only a few selected cases have been reanalyzed considering a finer mesh (3 mm × 3 mm), and the results are discussed. The difference in mesh sizes is depicted in Figure 1a,b (depicted on the geometry of the specimen 88.9 × 2.6-400-P [38]).

Modeling Approach #C
For these cases, initial local geometrical imperfections have also been considered (Section 3.2). The global imperfections and the mesh size are the same as in approach #A.

Modeling Approach #D
Several model cases have been recalculated by adopting the solid structural finite elements (SOLID 185) for the CHS tube instead of the shell elements ( Figure 1c). The mesh sizes in the longitudinal and tangential directions are the same as in approach #A. Three elements have been considered in the radial direction, along the cross-section thickness t. The default key option (0) of the element technology has been considered; therefore, there is a full integration with the B-bar method [39] (selective reduced integration) which helps to prevent volumetric locking in nearly incompressible cases.

Modeling Approach #E
For few cases, the uniform reduced integration of the solid finite elements with hourglass control has been adopted. One integration point per solid element is considered in this formulation. Note: key option No. 2 is set up as 1 [39].

Modeling Approach #F
To compare the results, the enhanced strain formulation of solid elements has been adopted for selected cases. Shear locking, as well as volumetric locking in nearly incompressible cases, are both prevented if this formulation is considered [39]. Certain numbers of internal degrees of freedom are introduced, and therefore, this option is less efficient. Note: key option No. 2 is equal to 2 [39].

Geometrical Imperfections
A global imperfection was incorporated into the FE model using the form of the lowest global buckling modal shape obtained from prior modal analysis (an example is given in Figure 2a). A global imperfection amplitude was used to simulate an initial global imperfection, as well as eccentricity. The amplitudes of ω 0 + e 0 (Tables 3-7) were considered. Local imperfections were generally neglected. Only in a few selected cases (#C), a local imperfection amplitude of t/10 was considered (where t is the section thickness). The local imperfection took a shape of the lowest local buckling modal shape (an example of 88.9 × 2.6-400-P [38] is in Figure 2b,c. Note: shape according to "b" was considered).

Material Model
To describe the behavior of the stainless-steel material, a stress-strain relation proposed by Ramberg and Osgood [10], modified by Hill [11] has been adopted: where σ and ε are engineering stress and strain, respectively, E0 is the material elastic Young's modulus, σ0.2 is the material 0.2% proof stress, and n is a strain-hardening exponent. At strains higher than σ0.2 value, this model overestimates the stress values [41]. A two-stage compound stress-strain curve devised by Mirambell and Real [42] provides a better agreement with experimental data for stress values above the 0.2% proof stress [41]. For the case of the compressive loading, a certain modification of the second stage is proposed by Gardner [33]: where σ1.0 is the 1% proof stress of the material, n'0.2,1.0 is a strain-hardening exponent, and E0.2 is the stiffness (tangent modulus) at the 0.2% proof stress given as: For the finite element (FE) numerical analyses, a multilinear material model with isotropic hardening has been considered (Mises plasticity). A more detailed description is presented in the author's previous study [14]. In this study, however, the stress-strain behavior only up to the value of 0.1 σ0.2 is considered as ideal elastic (to neglect plasticity at low strains), instead of the previously considered value of 2/3 σ0.2 [14] (which was too high for some cases). The engineering (nominal) stress-strain material curves have been transferred into true stress and logarithmic strain dependences to match the results of geometrically nonlinear FE analyses:

Material Model
To describe the behavior of the stainless-steel material, a stress-strain relation proposed by Ramberg and Osgood [10], modified by Hill [11] has been adopted: where σ and ε are engineering stress and strain, respectively, E 0 is the material elastic Young's modulus, σ 0.2 is the material 0.2% proof stress, and n is a strain-hardening exponent. At strains higher than σ 0.2 value, this model overestimates the stress values [41]. A twostage compound stress-strain curve devised by Mirambell and Real [42] provides a better agreement with experimental data for stress values above the 0.2% proof stress [41]. For the case of the compressive loading, a certain modification of the second stage is proposed by Gardner [33]: where σ 1.0 is the 1% proof stress of the material, n' 0.2,1.0 is a strain-hardening exponent, and E 0.2 is the stiffness (tangent modulus) at the 0.2% proof stress given as: For the finite element (FE) numerical analyses, a multilinear material model with isotropic hardening has been considered (Mises plasticity). A more detailed description is presented in the author's previous study [14]. In this study, however, the stress-strain behavior only up to the value of 0.1 σ 0.2 is considered as ideal elastic (to neglect plasticity at low strains), instead of the previously considered value of 2/3 σ 0.2 [14] (which was too high for some cases). The engineering (nominal) stress-strain material curves have been transferred into true stress and logarithmic strain dependences to match the results of geometrically nonlinear FE analyses: where σ nom is the nominal engineering stress, ε nom is the nominal engineering strain, and ε true is the true total (mechanical) strain. For the compressive material properties, ε nom have been introduced with negative values. As it is impossible to define a negative tangent of the stress-strain relation while adopting isotropic hardening [39], the stress-strain relation has been defined as ideal plastic (with its tangent slope close to 0 but positive) instead of any kind of softening after the peak stress. An example of material model verification (parameter values of 106 × 3-400-FR in Table 4 in [38]) is depicted in Figure 3.
) ε ln( ε nom true where σnom is the nominal engineering stress, εnom is the nominal engineering strain, and εtrue is the true total (mechanical) strain. For the compressive material properties, εnom have been introduced with negative values. As it is impossible to define a negative tangent of the stress-strain relation while adopting isotropic hardening [39], the stress-strain relation has been defined as ideal plastic (with its tangent slope close to 0 but positive) instead of any kind of softening after the peak stress. An example of material model verification (parameter values of 106 × 3-400-FR in Table 4 in [38]) is depicted in Figure 3. In cold-formed CHS, small values of the membrane residual stresses have been observed, and as such can be neglected [43]. The residual stresses through thickness are implicitly incorporated by considering the measured values of the material properties [44].

Boundary Conditions and Loading
To simulate a pin-ended column, all the circumferential nodes at the end of the CHS tube are connected (in the radial direction) with a single node located on the section axis. The offset of the node is either 50 (40 + 10) or 77 mm in dependence on the experimental setups described in [38] p. 301. The nodal connections are modeled using rather stiff beam elements. Such a boundary condition is considered at both ends of the CHS tube (see Figure 1). All translational and 2 rotational degrees of freedom of the bottom node are constrained. The loading during the FEA has been conducted by a prescribed displacement of the upper node (in the direction of the CHS tube axis). Two translational and 2 rotational degrees of freedom of the upper node were constrained.

Results
The results of all the 74 FEA (finite element analyses) of the approach #A (shell elements), the ultimate axial loads Nu,FE, and corresponding ultimate midheight lateral deflections ωu,FE, together with the experimental results (Nu,exp and ωu,exp) [38], are documented in Tables 8-12.
The results of variants in accordance with the approach #B or #C (shell elements, finer mesh variant, or local imperfection variant) are not documented in tables but rather only graphically (Figures 4-6).
In the case of the modeling approach #D (solid elements, selective reduced integration), two whole sets of lengths (two chosen cross-sections) have been reanalyzed, and the results are summarized in Tables 13 and 14. The difference between approaches #D and In cold-formed CHS, small values of the membrane residual stresses have been observed, and as such can be neglected [43]. The residual stresses through thickness are implicitly incorporated by considering the measured values of the material properties [44].

Boundary Conditions and Loading
To simulate a pin-ended column, all the circumferential nodes at the end of the CHS tube are connected (in the radial direction) with a single node located on the section axis. The offset of the node is either 50 (40 + 10) or 77 mm in dependence on the experimental setups described in [38] p. 301. The nodal connections are modeled using rather stiff beam elements. Such a boundary condition is considered at both ends of the CHS tube (see Figure 1). All translational and 2 rotational degrees of freedom of the bottom node are constrained. The loading during the FEA has been conducted by a prescribed displacement of the upper node (in the direction of the CHS tube axis). Two translational and 2 rotational degrees of freedom of the upper node were constrained.

Results
The results of all the 74 FEA (finite element analyses) of the approach #A (shell elements), the ultimate axial loads N u,FE , and corresponding ultimate midheight lateral deflections ω u,FE , together with the experimental results (N u,exp and ω u,exp ) [38], are documented in Tables 8-12.
The results of variants in accordance with the approach #B or #C (shell elements, finer mesh variant, or local imperfection variant) are not documented in tables but rather only graphically (Figures 4-6).
In the case of the modeling approach #D (solid elements, selective reduced integration), two whole sets of lengths (two chosen cross-sections) have been reanalyzed, and the results are summarized in Tables 13 and 14. The difference between approaches #D and #A was very negligible in the matter of monitored data, and therefore, the remaining geometries have not been reanalyzed considering solid model (the presented results in Table 8 are almost the same as those in Table 13, and the results of Tables 9 and 14 likewise). Statistical comparison of the monitored data (N u,FE /N u,exp and ω u,FE /ω u,exp values) between approaches #A and #D are summarized in Table 15. For approach #A, only the results which are presented in Tables 8 and 9 are applicable for this statistical  evaluation (as for approach #D, only the results of corresponding cross-section geometries  are available-Tables 13 and 14).
FE models of approach #A were validated by comparing the averaged normalized ultimate loads N u,FE /N u,exp and the averaged normalized deflections ω u,FE /ω u,exp (Table 16). The validation results are also compared with the results of Mr. Buchanan et al. [38].
Furthermore, the selected model cases of the different solid element formulations considered in approaches #E and #F have not seemed to yield significantly different results in the matter of the observed outputs and are therefore documented only graphically in this study. Global slenderness λ is calculated in dependence on the compressive cross-section class (according to EN 1993-1-4 [15]) considering Equation (6) for section class (cl.) 1-3, and Equation (7) for cross-section class 4: where A is the cross-section area, σ 0.2 is the 0.2% proof stress, E is Young's modulus, I is the second moment of area, L is the effective length, and A eff is the effective cross-sectional area determined in accordance with the formula from BS 5950-1 [45]: Compressive class slenderness limits have been considered in accordance with EN 1993-1-4 [15]. The global slenderness λ and the compressive classes (cl.) have been calculated for all the model geometries for both sets of material properties (SCP and TP) and are documented in Tables 8-14 (values of λ in Tables 13 and 14 are, of course, the same as in Tables 8 and 9, respectively).
The N-ω (axial load-midheight deflection) curves are provided in Figures 4-11. For each figure, the "a" part is used to depict the whole force-deflection curve, and the "b" (zoomed) part of each Figure focuses on the area near the peak forces N u,exp, and N u,FE , so the negligible differences of the N-ω curves in these areas between considered modeling approaches are also visible. In some cases, the "a" part of the figure is described as "whole*"-in these cases, certain N-ω curves of the graph are not depicted until the very last converged or analyzed step of that case.
The differences in the matter of the force-deflection curves between modeling approaches #A, #B, and #C (shell numerical model with default mesh size, finer mesh size, and local imperfection, respectively) are depicted for selected specimens of the crosssectional set 88.9 × 2.6 in Figures 4-6.
The differences in the matter of the force-deflection curves between modeling approaches #A and #D (shell vs. solid numerical model) are depicted for the selected specimens of the cross-sectional set 106 × 3 and 104 × 3 in Figures 7-11. Moreover, in Figures 10 and 11, the results of the approaches #E and #F (different solid element formulations) are also documented and available for comparison. each figure, the "a" part is used to depict the whole force-deflection curve, and the "b" (zoomed) part of each Figure focuses on the area near the peak forces Nu,exp, and Nu,FE, so the negligible differences of the N-ω curves in these areas between considered modeling approaches are also visible. In some cases, the "a" part of the figure is described as "whole*"-in these cases, certain N-ω curves of the graph are not depicted until the very last converged or analyzed step of that case.    each figure, the "a" part is used to depict the whole force-deflection curve, and the "b" (zoomed) part of each Figure focuses on the area near the peak forces Nu,exp, and Nu,FE, so the negligible differences of the N-ω curves in these areas between considered modeling approaches are also visible. In some cases, the "a" part of the figure is described as "whole*"-in these cases, certain N-ω curves of the graph are not depicted until the very last converged or analyzed step of that case.      The differences in the matter of the force-deflection curves between modeling approaches #A, #B, and #C (shell numerical model with default mesh size, finer mesh size, and local imperfection, respectively) are depicted for selected specimens of the cross-sectional set 88.9 × 2.6 in Figures 4-6.    The differences in the matter of the force-deflection curves between modeling approaches #A and #D (shell vs. solid numerical model) are depicted for the selected specimens of the cross-sectional set 106 × 3 and 104 × 3 in Figures 7-11. Moreover, in Figures 10  and 11, the results of the approaches #E and #F (different solid element formulations) are also documented and available for comparison.  As described in [38], three specimens from those of shorter lengths (104 × 2-550-P, 104 × 2-750-P, and 106 × 3-550-P) changed the initial direction of the midheight lateral deflection (induced by the applied eccentricity) before reaching the peak force Nu,exp (hence the minus values of ωu,exp in Tables 8 and 9). One specimen (80 × 1.5-300-P) changed the deflection direction after the peak force Nu,exp had been reached (hence positive value of the ωu,exp for this case in Table 11). Note: the deflection reference data of the specimen 104 × 2-550-P was multiplied by −1 only for the purpose of depiction in Figure 8b. Original data were considered for the statistical evaluation. In the case of modeling approach #F (solid elements, enhanced strain formulation), the convergence of the numerical solution was more sensible to the step size of the analysis. When the step size in approach #F was adopted, the same as in all the other approaches, the last converged substeps were at the values of midheight lateral deflection ω approx. Six mm and 11.3 mm for specimens 104 × 2-950-P and 104 × 2-1650-P, respectively, (orange dash-dotted curves in Figures 10 and 11; note: the last point of the orange curve in Figure 11 is also a converged state). To obtain further analysis data, a finer step size has been adopted for the case 104 × 2-1650-P, and the results are described by the pink curve in Figure 11. However, in the postpeak region of this N-ω curve, certain oscillations have As described in [38], three specimens from those of shorter lengths (104 × 2-550-P, 104 × 2-750-P, and 106 × 3-550-P) changed the initial direction of the midheight lateral deflection (induced by the applied eccentricity) before reaching the peak force N u,exp (hence the minus values of ω u,exp in Tables 8 and 9). One specimen (80 × 1.5-300-P) changed the deflection direction after the peak force N u,exp had been reached (hence positive value of the ω u,exp for this case in Table 11). Note: the deflection reference data of the specimen 104 × 2-550-P was multiplied by −1 only for the purpose of depiction in Figure 8b. Original data were considered for the statistical evaluation.
In the case of modeling approach #F (solid elements, enhanced strain formulation), the convergence of the numerical solution was more sensible to the step size of the analysis. When the step size in approach #F was adopted, the same as in all the other approaches, the last converged substeps were at the values of midheight lateral deflection ω approx. Six mm and 11.3 mm for specimens 104 × 2-950-P and 104 × 2-1650-P, respectively, (orange dash-dotted curves in Figures 10 and 11; note: the last point of the orange curve in Figure 11 is also a converged state). To obtain further analysis data, a finer step size has been adopted for the case 104 × 2-1650-P, and the results are described by the pink curve in Figure 11. However, in the postpeak region of this N-ω curve, certain oscillations have developed. The amplitudes of these oscillations are not significantly higher than the oscillation amplitudes of the reference data (red curve). Further analyses of the approach #F have not been conducted due to much lower computational efficiency, as well as not a very significant difference in the matter of the monitored results (ultimate axial load N u,FE, and the corresponding midheight lateral deflection ω u,FE ).
The differences in the results for approach #E are rather negligible near the areas of the ultimate axial loads N u,FE and noticeable only in the areas of large deformations. Therefore, further analyses considering this approach have not been conducted, either.
Contour plots of the equivalent plastic strains and, more importantly, deformed geometries for selected cases are depicted in Figures 12-14. For each contour plot, the values of deformation load u z and the corresponding axial load N z are provided. Furthermore, the modeling approach is stated, and the link of the case with the corresponding F-u curve and figure number of the graph is provided.
Global buckling was the most common failure mode in the case of longer specimens (e.g., Figure 12c). Shorter specimens have rather developed a local buckle near the midheight of the compressed side of the tubular cross-section after the peak load has been reached (e.g., Figure 12a,b). The buckling modes are very similar to those described in the study by Buchanan et al. [38].  (a) 88.9 × 2.6-400-P (SCP) #C (b) 88.9 × 2.6-400-P (TP), finer step #B (c) 88.9 × 2.6-3080-P (SCP) #A Figures 13 and 14 depict the differences in the deformed shape at the compressed face near the midheight of the selected specimens when adopting the modeling approach of #A (shell elements), #D (solid elements), or #E (solid elements with uniform reduced integration). The area near the local buckling is deformed differently in these cases. Figure 15 depicts energy ratios of the case 104 × 2-1650-P (SCP) modeled by approaches #A and #E. The energy ratios are monitored at states of the last analyzed steps (Figure 15a,b), and also at the states of the maximal peak axial loads Nu (Figure 15c,d). In the course of the computation, the artificial energy is introduced to sustain hourglass control of the reduced integration elements, either shells from case #A or solids from case #E. In the case of modeling approaches #D (selective reduced integration) or #F (full integration), no artificial energy is being introduced. The solution is generally acceptable if the ratio of artificial energy to the total energy is less than 5% [39]. Mesh refinement is recommended otherwise. Figure 15 depicts the ratio of the artificial to total energy for each element. In the case of approach #A, the recommended value of 5% has been exceeded in the last step of the analysis only locally, with the maximal value of 23% in two shell elements (Figure 15a). However, in the modeling approach #E, the limit value has been exceeded in a significant number of the finite elements (Figure 15b). In the step where the peak loads Nu have been achieved (Figure 15c for approach #A and Figure 15d for the approach #E), the energy ratios are far below the recommended value of 5%, with generally lower values for the modeling approach case #A (shell elements). Note: pay attention to the different contour scales.  Figures 13 and 14 depict the differences in the deformed shape at the compressed face near the midheight of the selected specimens when adopting the modeling approach of #A (shell elements), #D (solid elements), or #E (solid elements with uniform reduced integration). The area near the local buckling is deformed differently in these cases. Figure 15 depicts energy ratios of the case 104 × 2-1650-P (SCP) modeled by approaches #A and #E. The energy ratios are monitored at states of the last analyzed steps (Figure 15a,b), and also at the states of the maximal peak axial loads N u (Figure 15c,d).
In the course of the computation, the artificial energy is introduced to sustain hourglass control of the reduced integration elements, either shells from case #A or solids from case #E. In the case of modeling approaches #D (selective reduced integration) or #F (full integration), no artificial energy is being introduced. The solution is generally acceptable if the ratio of artificial energy to the total energy is less than 5% [39]. Mesh refinement is recommended otherwise.   Figure 15 depicts the ratio of the artificial to total energy for each element. In the case of approach #A, the recommended value of 5% has been exceeded in the last step of the analysis only locally, with the maximal value of 23% in two shell elements (Figure 15a). However, in the modeling approach #E, the limit value has been exceeded in a significant number of the finite elements (Figure 15b). In the step where the peak loads N u have been achieved (Figure 15c for approach #A and Figure 15d for the approach #E), the energy ratios are far below the recommended value of 5%, with generally lower values for the modeling approach case #A (shell elements). Note: pay attention to the different contour scales.

Validation Data (Approach #A)
The average and COV (coefficient of variation) values of normalized ultimate load N u,FE /N u,exp and deflection ω u,FE /ω u,exp are depicted in Table 16 (based on the values from  Tables 8-12). These results are in good agreement with the research of Mr. Buchanan et al. [38] (also depicted in Table 16). The mean value of the normalized ultimate load N u,FE /N u,exp is close to one for both sets of the material properties, TP and SCP (tensile properties, stub-column properties), and the COV values are close to zero; therefore, it is possible to consider the accuracy as sufficient. COV is smaller when the material set of SCPs has been adopted in the FEA, the same as observed in the study by Buchanan et al. [38]. The value of the normalized ultimate deflection ω u,FE /ω u,exp is closer to one when the SCP set of material properties has been considered, again in accordance with the conclusions in [38]. The values of ω u,FE /ω u,exp are 0.627 and 0.660 for TP and SCP, respectively, and are in accordance with the result values of 0.612 and 0.708 from [38]. These values seem to be rather far from the accurate match (the value of 1), with quite a high coefficient of variation: 1.006-0.825 in this study and 0.920-0.834 in the study [38]. This is firstly due to the flat shape of the load-deflection curves near their maximum. Secondly and more importantly, as a consequence of the lateral deflection in the direction opposite to the initial geometrical imperfection before reaching the peak value of the axial force, which was observed in three (out of 37) reference cases, hence considered with the negative value (also see [38]). Note: in the fourth reference case, the deflection evolved into the opposite direction after reaching the peak load, and therefore, is not considered with the negative value of the lateral deflection [38].
Statistical values, the average (mean) value of the ultimate load N u,FE /N u,exp , and deflection ω u,FE /ω u,exp along with the standard deviations bars are graphically depicted in Figure 16 (based on Table 16). These values are based on the results of this study and of the study by Mr. Buchanan [38]. of 37) reference cases, hence considered with the negative value (also see [38]). Note: in the fourth reference case, the deflection evolved into the opposite direction after reaching the peak load, and therefore, is not considered with the negative value of the lateral deflection [38].
Statistical values, the average (mean) value of the ultimate load Nu,FE/Nu,exp, and deflection ωu,FE/ωu,exp along with the standard deviations bars are graphically depicted in Figure 16 (based on Table 16). These values are based on the results of this study and of the study by Mr. Buchanan [38].
In general, the validation data are very similar to the values discussed in [38], and therefore, the considered numerical FE model is possible to be considered as validated with no dispute. Consistency between FEA results and the experimental data has been achieved (especially for the case of SCP material properties depicted in Figure 8).
(a) (b) Figure 16. Average values with standard deviation bars for: (a) normalized ultimate loads; (b) normalized ultimate deflections; both depicted for the results considering SCP and TP material properties, based on this study and also the study by Mr. Buchanan [38]. In general, the validation data are very similar to the values discussed in [38], and therefore, the considered numerical FE model is possible to be considered as validated with no dispute. Consistency between FEA results and the experimental data has been achieved (especially for the case of SCP material properties depicted in Figure 8).

Mesh Size Influence (Approach #A vs. #B)
The difference between these two mesh sizes, coarser (8 mm × 5 mm, approach #A) or finer (3 mm × 3 mm approach #B), has been observed for the geometries of three different slenderness values from the set of the same cross-section (Figures 4-6, crosssection 88.9 × 2.6), considering the lowest structural length from the set (400 mm, Figure 4), the largest one (3080 mm, Figure 6), and the third value from (950 mm, Figure 5). This choice of various lengths for one selected cross-section is being considered as a sufficiently representative sample (all the cross-sections of this study are geometrically very similar).
Only in the case of little slenderness, the initial model stiffness is considerably larger when the finer mesh has been adopted ( Figure 4). This might be the reason for the adopted boundary conditions as well (the connection of the constrained nodes with the circumferential nodes at the CHS tube ends), with possibly larger influence in cases of the shorter specimens. The influence is greater for finer meshes due to the utilized automatized APDL macro (finer mesh = more circumferential nodes = more beams to connect the circumferences with the appropriate loading nodes). This issue is not considered a concern, as the real stiffness of the boundary conditions was not determined with certainty. In addition, the postpeak behavior is slightly different, when the finer mesh has been considered only in a few cases ( Figure 5). In some model cases, the analyses considering the finer mesh had more convergence difficulties while the same time step was considered (e.g., Figure 4, #B, SCP).
Based on the results of these selected geometries, the difference between various mesh sizes is very negligible in the matter of the monitored results (N u , ω u ), which are the principal interest of this study.

Local Imperfection Influence (Approach #A vs. #C)
Influences of the local initial imperfections are investigated for the same selected cases as described previously (Section 5.2). The difference between cases with adopted and neglected local imperfection is very small (Figures 4-6) for all the cases. Therefore, the local imperfections have not been implemented for the other cases to conduct the general validation of the numerical FE model. The effort invested in the implementation of the local initial geometrical imperfections is not meaningful for the objectives of this study. For different lengths of the column specimens, different modal shapes might be applicable to incorporate the local imperfection. Therefore, a manual check of the modal shapes themselves is required. On the other hand, the global imperfection shape has always been obtained from the first modal shape, which is much more feasible if the analysis process is to be automatized (e.g., via APDL macros).

Shell or Solid FE Model (Approach #A vs. #D)
The comparison between the shell model and the solid model in the matter of the globally monitored results (normalized ultimate load N u,FE /N u,exp and deflection ω u,FE /ω u,exp ) are depicted in Table 15. In the case of approach #D (solid model), two whole sets of lengths (considering two different cross-sections) have been reanalyzed. These results are summarized in Tables 13 and 14 (and are very similar to the corresponding results of case  #A, Tables 8 and 9, respectively). The mean values and the values of the coefficient of variations are practically the same for both approaches (Table 15). Load-deflection curves are compared in Figures 7-11. The behavior in the postpeak softening region (where the large strains are involved) is different for some modeled cases. However, this is not the objective of this study.

Full Integration of the Solid Elements (approach #F)
This approach has been investigated for only two selected cases, the cross-section of 104 × 2, and structural lengths of 950 mm ( Figure 10) and 1650 mm ( Figure 11). When the same step size (and the same maximal number of Newton-Rhapson iterations within one step) has been considered, the convergence problems occurred shortly after reaching the peak force N u,FE . A finer time step needed to be implemented to achieve the convergence after the peak load. This has been introduced for one case ( Figure 11). Although the convergence has been secured, certain oscillations have developed. The amplitude of these oscillations is not significantly higher than the oscillation amplitudes of the reference data, and therefore, the results are considered sufficiently good. Further analyses considering this approach, however, have not been conducted due to much higher computational demands. Moreover, the investigated results (N u , ω u ) are very similar to those obtained by much more efficient approaches.

Uniform Reduced Integration (Approach #E), Selective Reduced Integration (Approach #D), Shell Elements with Reduced Integration (Approach #A)
The differences in the monitored results are very negligible, again (Figures 10 and 11). A noticeable difference is in the shape of the local buckling area in the midspan of the column height (Figures 13 and 14). The distortion in this area seems to be rather smoother for cases where the #E approach has been considered (Figure 13c or Figure 14c), in comparison with the sharper area of local buckling when either the #D approach (Figure 13b or Figure 14b) or the default #A approach (Figure 13a or Figure 14a) have been adopted.
The energy ratio (artificial to total energy) [39] is depicted in Figure 15. In order to consider the results, where the elements with the reduced integration have been considered as good, this energy ratio needs to be below 0.05 (5%), as recommended in [39]. The ratio is below the recommended value at the time of the peak load for both cases, #A ( Figure 15c) and #E (Figure 15d), with slightly better (lower) values for approach #A. For the last observed step, however, the energy ratio exceeds 0.05 in a large number of the elements when approach #E has been adopted (Figure 15b), and only in few elements in the case of approach #A (Figure 15a). For the modeling approach #D, there is no additional artificial energy (this approach is considered to be more precise). Therefore, the postpeak results of approach #A are considered more convenient.

Other Observations
In all of the modeled cases, the region with the largest plasticity has developed in the midspan, on the compressed face of the columns. Examples are depicted in Figure 12. The shape of this area was slightly different when the finer mesh has been adopted (Figure 12b), compared to the coarser mesh of the same geometry (Figure 12a). In addition, the difference in the shape of the high plasticity area is between cases of various structural lengths (Figures 12 and 13). For the large lengths (Figure 12c), the global buckling occurred much sooner, and the analyses have not been conducted long enough to observe any local distortion of the cross-sectional shape.
Although proper statistical evaluation among all the discussed modeling approaches is not provided, further cases have not been reanalyzed, as the results of approach #A are very satisfactory, compared to the results documented in the study by Mr. Buchanan [38], which was the main objective of this study.

Further Application
In future work, the initial imperfections will be considered as random variables. This will make it possible to study their influence on the ultimate limit state using stochastic sensitivity analysis [46][47][48][49], where the main influences and interaction effects between imperfections can be the subject of research.

Conclusions
The performed simulations aimed at the comparison of advanced numerical models and the determination of their limitations in the simulation of the pressure test of thin-walled stainless-steel columns, considering the loss of stability. The study shows a significant effect of the correct description of boundary conditions on the overall response. In addition to the comparative (normalized ultimate load N u,FE /N u,exp and deflection ω u,FE /ω u,exp ), several generalizable conclusions could be obtained by cross-comparison, which is discussed in detail in Section 5. Experimentally determined data from extensive research were considered as reference values. Very good agreement was reached for all used models.
The mean (average) value of the normalized ultimate axial load N u,FE /N u,exp was 1.010 (standard deviation 0.074) in this study and 1.024 (standard deviation 0.083) in the study by Mr. Buchanan, for the numerical analyses based on the "tensile test material properties" (TP) values. In the case of numerical simulations based on the values of the "stub-column material properties" (SCP), the average normalized ultimate axial load was 0.998 (standard deviation 0.050) in this study and 1.013 (standard deviation 0.047) in the study by Mr. Buchanan et al. These statistical values are graphically depicted in Figure 16a (and also in Table 16). The values of normalized ultimate forces are very close to 1, and the standard deviation is also low. The validation results for the ultimate forces presented here are in good agreement with the experimental study.
In  Table 16). The values seem to be rather scattered (normalized average not close to 1.00, large deviation). However, these values are very similar to those from the study by Mr. Buchanan et al., where the finite element models are successfully validated. The large deviation of values is the reason for the flat nature of the load-deflection curve near its maximal value of the ultimate load; therefore, larger values of standard deviations are expected.
Another objective was to describe a robust design procedure with the appropriate use of numerical simulation. This was achieved by a comparison of several modeling approaches. The modeling approaches differ in the adopted type of finite elements (either shell or solid), in varying element formulations (e.g., full integration or reduced integration), in varying mesh sizes (finer, coarser) or local geometrical imperfection incorporation or absence. The effects of these modeling approaches on the monitored main objective model outputs (ultimate axial load N u and the corresponding midheight lateral deflection ω u ) are relatively very small and therefore negligible. Considering the shape of the whole force-deflection curve and the local buckling area, the type of adopted finite elements along with their formulation set up have the highest importance on the numerical model response. The most suitable, and also the most effective, are shell elements. The influence of the mesh size (considering two types of mesh, either the finer 3 mm × 3 mm or coarser 5 mm × 8 mm) is slightly more important for shorter specimens than for the longer ones. The most negligible influence on the results, in general, is caused by the effect of the local geometrical imperfection within the considered amplitude value of t/10 (where t is the section thickness-various value for various cross-sectional sets).
The validated numerical FE models presented here will be used for future analyses, where a study of initial imperfections influence on the ultimate limit state using stochastic sensitivity analysis will be conducted. These are simulations of large sets of samples, and therefore, it is necessary to choose the most concise computational model with the minimum time requirement.
Author Contributions: Conceptualization, methodology, analysis, and writing-original draft preparation, D.J.; writing-review and editing, supervision, J.K. and Z.K.; project administration, funding acquisition, Z.K. All authors have read and agreed to the published version of the manuscript.
Funding: This paper has been created with the financial support of Czech Science Foundation by project No.: 20−00761S "Influence of material properties of stainless steels on reliability of bridge structures". Another financial support was provided by the BUT fund for specific university research, the project number FAST-J-21-7169. Data Availability Statement: Data sharing not applicable.