Distortional-Global Buckling Interaction Relevance in Cold-Formed Steel Lipped Channel Columns

BUCKLING MODE INTERACTION IN COLD-FORMED STEEL (CFS) MEMBERS MUST BE CONSIDERED FOR THE STRUCTURAL DESIGN, SINCE THE MODAL INTERACTION MAY LEAD TO A SIGNIFICANT REDUCTION OF THE STRUCTURAL STRENGTH, USUALLY RECOGNIZED AS EROSION OF THE LIMIT LOAD. SO FAR, THE DISTORTIONAL-GLOBAL (DG) BUCKLING INTERACTION IS NOT COVERED BY CODES, E.G. BRAZILIAN CODE NBR 14762 (2010). THE PRESENT INVESTIGATION IS AIMED AT THE DG INTERACTION OF CFS LIPPED CHANNEL (LC) COLUMNS, WHICH IS THE MOST USUAL SECTION IN MANY APPLICATIONS. THE METHODOLOGY EVOLVES A PARAMETRIC ANALYSIS, OVER A SINGLE LC COLUMN UNDER DG BUCKLING INTERACTION. FIRST, A STUDY OF INITIAL GEOMETRIC IMPERFECTION (IGI) SENSIBILITY IS PERFORMED, WITH THE PURPOSE OF UNDERSTANDING THE INFLUENCE OF THE IGI ON THE COLUMN’S BEHAVIOR. MOREOVER, THE PARAMETRIC ANALYSIS IS EXTENDED TO A SET OF YIELDING STRESS AND COLUMN LENGTHS, IN ORDER TO UNDERSTAND THE ULTIMATE LOAD UNDER DIFFERENT TYPES OF DG BUCKLING INTERACTION NATURES. CONCLUSIONS ON THIS RESEARCH HAVE BEEN SHOWN THAT THE ACTUAL GLOBAL BUCKLING EQUATION FROM THE DIRECT STRENGTH METHOD (DSM) IS ALREADY SUITABLE TO COVER THE DG BUCKLING INTERACTION FOR THE CASE OF LIPPED CHANNEL MEMBERS UNDER AXIAL COMPRESSION.


Steel cold-formed thin-walled structures and buckling modes
Steel cold-formed members are widely applied in structural solutions, e.g. light steel frame construction, storage rack systems, mezzanine floors and building roofing systems. A large variety of cold-formed steel sections, CFS, can be manufactured with the help of light gage cold forming machines, supported by CAD/CAM facilities allowing performing simple and complex geometries, including holes and slotted perforations. Besides the competitive advantages of CFS structural systems, designers must take into account their particular structural behavior on the basis of the Vlasov (1961) theory of thin-walled members, with special attention to the warping torsional behavior and the buckling modes. In this context, the Brazilian standard ABNT NBR 14762 (ABNT, 2010) includes procedures for the structural design of CFS members, and the buckling effect on the strength of the members is of utmost importance.
CFS members with open sections display local (L), distortional (D) and global (G) buckling (Figure 1), with the respective axial compression critical loads crL P , crD P and crG P . The identification of the critical buckling loads and the associated modes can be obtained with the help of a numerical solution based on the finite strip method (FSM), with practical advantages when compared with the finite element method (FEM), due to a significant reduction in the degrees of freedom of the structural model. In addition, the FSM allows automatic buckling analysis along a defined member length range, allowing defining the variation of the buckling loads and the associated buckling modes, which is clearly presented to the user as the "signature curve" of the section. This arrangement is not easily available for FEM and would be much more time consuming if compared with the FSM. The computational program CUFSM, is a well-established reference by Li and Schafer (2013) and Ádany and Schafer (2006), and offers access to the signature curve (namely cr P vs. length L of the member) of the CFS section, allowing the identification of the buckling modes of members under axial compression, bending or any kind of initial stress distribution. Consequently, the identification of the single buckling modes, L, D and G, is an important step for the structural design of thin-walled CFS members. In addition, the buckling modes interaction must be considered for the structural design, since the modal interaction may lead to an important reduction of the structural strength, usually recognized as erosion of the limit load.

Distortional-global interaction: relevance & concepts
In the literature so far, the coupled instability phenomenon can be found as compound buckling, simultaneous buckling, interaction buckling and many other denominations. Historically, Koiter (1945) and Budiansky (1974) were the first two researchers that included the coupled instability phenomenon.
The nature of the phenomenon is complex and difficult to predict. Due to that, many authors try to classify the buckling interaction in various ways. Gioncu (1994) shows an overall of the couple instability state-of-art from over 230 papers until 1994, explaining concepts and classifications for the nature of the coupled phenomenon.
The distortional-global (DG) interaction can be classified by its type and nature. As reported by Martins, Camotim and Dinis (2018a), the type of distortional-global interaction can change depending on the type of the global buckling mode: (i) distortional/major-axis flexural-torsional (DMFT) and (ii) distortional/minor-axis flexural (DmF). However, for lipped channel columns, the DMFT interaction is more common to occur, while the DmF interaction happens more often on zed-sections columns. In addition, related to the nature of the DG buckling interaction, Martins, Camotim and Dinis (2018a) states a classification based on the ratio between global and distortional critical loads, , which gives the basis for the identification of the presence of the DG buckling interaction. Following the same idea, in the present study the slenderness ratio is considered, i.e.
, since the slenderness factor λ is commonly applied to design procedures and, in addition, gives a more clear identification of the effect of buckling behavior of CFS members. The categories proposed by Martins, Camotim and Dinis (2018a), with the slenderness parameter GD R and GD R λ , are classified as: (i) True DG interaction (TI): when the distortional and global critical buckling are close; this scenario may always happen for .
. GD 0 95 R 1 05 λ < < . According to Martins, Camotim and Dinis (2018a), this state can behave differently in three different groups, depending on the critical slenderness: (i.1) abrupt collapse for stocky columns, (iii) Secondary-global bifurcation DG interaction (SGI): also occurs for yield stress sufficiently high, enabling the buckling interaction to arise. Nonetheless, this one is more likely to develop due to its high post-critical strength reserve caused by the distortional buckling. This condition is more commonly observed in the range , and it is possible to observe other buckling modes interaction of the lipped channel CFS column, as LDG or just LD.

Design approaches for cold-formed steel columns under DG interaction
The direct strength method (DSM) is the most accepted design approach because of its simple application for civil engineers' design of thin-walled structural steel elements, based on the limit states associated with local, distortional, global and the local-global interactive buckling modes. The DSM is based on the Winter-type equation (Winter, 1968) & (Winter, 1947) for the local and distortional strength curves, and the "classical" column design curve for the global mode, taken from specifications of hot-rolled steel columns (e.g. ABNT NBR 8800 (ABNT, 2008) and ANSI/AISC 360-16 (2016) where y P is the squash load, defined as the gross sectional area times the steel yielding stress; G λ is the global slenderness; crG P is the elastic critical global buckling load. The DSM nominal axial strength for the distortional buckling, as described in the DSM, is expressed in Eq. (2).
where D λ is the distortional slenderness and crD P is the elastic critical distortional buckling load. One advantage of the DSM procedure is the capacity to express straightforward interactive equations, as shown by Schafer (2002). In the past few years, additional equations have been proposed considering interactive modes and also recalibration of the Winter-type equation (1968) coefficients for global and distortional equations, e.g. Dinis et al. (2020), Liu et al. (2020), and Santos et al. (2020).
Unlike the LG coupled phenomenon, the DSM-based nominal axial strength for the distortional-global (DG) interactive buckling has not been included in the standards. In the past few years, authors have been studying the DG Latin American Journal of Solids and Structures, 2020, 17(9), e332 4/30 coupled phenomenon behavior, as shown by  and Martins et al. (2018a). However, there is still a lack of studies of this phenomenon, mainly laboratory experiments. Schafer (2002) proposed an approach to the DSM that considers the DG interaction, shown in Eq. (3). A recent study carried out by Martins, Camotim and Dinis (2018a) tested the procedure proposed by Schafer (2002) on the basis of a parametric study, showing that this procedure is quite conservative, when compared with FEM results. Martins, Camotim and Dinis (2018a) also studied an additional approach that includes the DG interaction, as shown in Eq. (4).

GD
where nG P is the DSM nominal axial strength for the global buckling presented by Eq. (1); DG λ is the distortional slenderness based on the global strength nG P ; nD P is the DSM nominal axial strength for the distortional buckling described by Eq. (2); GD λ is the global slenderness based on the distortional strength nD P .
The present investigation is aimed at the DG interaction of CFS columns and one must observe that the buckling mode interactions develop the closer are critical loads: crD P and crG P must be close, at the same time crL P must be high enough to avoid the development of a combination of modes that includes the local buckling one. In the present case of DG, the investigated CFS column is the lipped channel, the most usual section in many applications, which is illustrated in Figure 1 with the identification of the fundamental buckling modes L, D and G.

NUMERICAL MODEL DESCRIPTION
For a numerical investigation of DG buckling interaction, a Finite Element Method (FEM) analysis is performed in order to understand the structural behavior (with assistance of ANSYS Mechanical APDL (ANSYS, 2016)). The FEM is addressed to capture and observe the nonlinear equilibrium path and detect the strength of the structural element.
Latin American Journal of Solids and Structures, 2020, 17(9), e332 5/30 The adopted model is described in six different groups: (i) discretization, (ii) end boundary condition, (iii) loading, (iv) material model, (v) initial geometric imperfections (IGI) and (vi) analysis method. For each group, the model is explained in detail, according to the ANSYS Mechanical APDL reference (ANSYS, 2016) and considerations retrieved from the literature. Figure 2 shows an illustration of a CFS column with a lipped channel section and both ends fixed. This figure illustrates some of the main components of the numerical model, such as discretization, end boundary condition and loading. Figure 2-a illustrates the finite element meshes, as well as the actual thickness of the end plate and the structural CFS lipped channel member.

Discretization
The adopted finite element type is the ANSYS SHELL281 (ANSYS, 2016). According to the ANSYS Theory Reference (ANSYS, 2016), the SHELL281 has 8 nodes, with 6 degrees of freedom per node, and is appropriate for linear, large rotation and large strain nonlinear situations. In addition, the SHELL281 formulation is placed on logarithmic strain and true stress measures. For this research, the finite element option contemplates the shell structural stiffness with bending and membrane considerations.
The mesh generation included quadrilateral and triangle-shaped elements with element size of 5 mm. This mixed mesh occurs resulting from the predefined geometric initial imperfections defined with the help of FStr Computer Application Lazzari (2020) (described in section 2.5) and is also due to the automatic mesh generation of the end plates of the column. Although the element shape chosen for the whole column is quadrilateral, some triangle-shaped elements may emerge, when the mesh is generated. This happens because the mesh generation of the structural member is formed by non-planar surfaces. According to ANSYS Theory Reference (ANSYS, 2016), SHELL281 provides more reliable results for triangular elements, while SHELL181 is not recommended for triangular-shaped elements. Figure 3 shows the mesh at the end plate of the column, detailing both types of possible mesh occurrences, quadrilateral and triangular shell element.

End boundary conditions
With respect to the boundary conditions, the column has a fixed-fixed end condition. The reason for this assumption is the actual condition of experimental tests, with fixed-fixed end sections condition, since fully pinned-pinned supports are more difficult to obtain in the laboratory (free bending in both minor and major axis). In addition, a fixed-fixed warping condition would always occur at the ends of the column under compression test. The constraints at the ends of the column are designed with a stiff plate rigidly fixed to the column cross-section, as illustrated in Figure 2-a.
In order to follow these support conditions, the FEM model includes constraints of the displacements in x and y directions, as well as the rotations around x, y and z direction at both ends of the column (see Figure 2-b and Figure 2-e). Another displacement constraint is placed at mid length of the column, at the center of the web, which prevents displacements in the longitudinal direction z (Figure 2-d). This additional constraint in the middle of the columns is needed to avoid free body translation, considering that symmetric behavior develops in the longitudinal direction of the column, imposed by the symmetric end boundary condition.

Loading
Equal compressive loads are applied at the ends of the column (see Figure 2-b and Figure 2-e). The external loads are applied at both ends as concentrated loads over the external face of the thick end plates, aligned with the centroidal axis of the column. The non-linear analysis, up to the column strength capacity, is based on loading steps. The maximum load is set up to be 1 N/N, and the first load step is 0.05 N/N for the arc-length method (described in section 2.6 Analysis methods).

Material model
The material model is defined as a balance of simplicity and strict accordance with the actual mechanical performance of the steel, in order to obtain an accurate parametric study. In this case, the material model with bilinear isotropic hardening (elastic-perfect plastic model) is considered, which applies the Von Mises yield criteria with an isotropic work hardening rule. This material model considers the initial slope of the stress-strain curve with the expected elastic Young modulus ( E 200 GPa = ) in the elastic strain region and a tangent slope modulus with zero magnitude, in the plastic strain region ( Et 0 MPa = ). The yield strength of the material varies, depending on the objective of the numerical test of the column, allowing variations of the column slenderness including all the buckling modes L, D and G. Also, the plates rigidly fixed to the ends of the column have the same mechanical properties as the CFS column. In addition, the major Poisson's ratio is n = 0.3 and an engineering stress-strain model.
With respect to the effects of the residual stresses and the cold-forming of the corners of the steel lipped channel member, the major part of investigations in the literature have neglected these effects. Dinis and Camotim (2015) has shown an insignificant impact in the numerical ultimate column load. Also, according to Ellobody and Young (2005), the small membrane residual stresses have demonstrated an irrelevant effect on the numerical ultimate load, stiffness of the column, load-shortening behavior and in the failure mode. In , a recent finite element modeling of local-distortional coupled phenomena in lipped channel CFS columns, the residual stresses and rounded corners effects were also neglected and the numerical FEM results indicate an accurate agreement with experimental results.

Generation of the initial geometric imperfections with FStr computational program
For the present research the original perfect geometry and imperfections of the column are created with the help of a previous finite strip method (FSM) analysis. The FStr computer application Lazzari (2020) performs previous elastic buckling analysis and, with the obtained modal critical shape, generates a node-mesh that is inserted into an APDL routine code as KEYPOINTS. The modal critical shape is inserted as the initial geometric imperfection (IGI), with a maximum amplitude depending on the buckling mode. In this case, the initially imperfect geometry of the column is composed of non-planar surfaces, connecting 4 nearby KEYPOINTS. Figure 4 shows the structural geometry generation with predefined imperfections given by KEYPOINTS produced by the FStr Computer Application.
The maximum amplitude parameters for the initial geometric imperfections for global and distortional buckling shapes are distinct. As proposed by Martins, Camotim and Dinis (2018a), the present research adopts / L 1000 for global imperfection and . 0 94t for distortional imperfection. The latter was originally proposed by Schafer and Peköz (1998a), as corresponding to a 50% probability that a random imperfection amplitude is below or above this value. Moreover, Santos (2017) has shown, in a study of lipped channel CFS columns under distortional buckling, that the ultimate load has not changed excessively (failure load variation below 5%) with the maximum amplitude of the distortional buckling geometric shape ranging from . 0 1t up to . 1 0t .

Analysis methods
For the present case the material and geometry nonlinearities are considered in a post-buckling behavior. The ANSYS built-in analysis method is based on the arc-length method, known as "The modified Riks method" (introduced by Riks (1979), Riks (1972) and Wempner (1971) with modifications being carried out later by many authors, e.g. Crisfield (1981), Crisfield (1983)). The method has a displacement control, during the loading, in order to find the fundamental paths before and after the limit point. This method is convenient for solutions of unstable problems that has a nonlinear static equilibrium.

Global buckling mode validation
The global (flexural-torsional) buckling validation is performed using three columns tested by Heva (2009), described in Table 1. This author has performed laboratory tests of multiple columns under different temperatures developing the flexural-torsional buckling. Updated information about the experiments can be found in Gunalan et al. (2014). For the present research, only the results for the specimen with room temperature (20 o C) are considered. Table 1 shows all the parameters and results of the numerical and laboratory tests performed by Heva (2009) and also the results obtained in the present research. Also, Table 1 includes all the measured geometries of the specimen, the material properties given by standard tensile tests and the measured specimen imperfections. Detailed information about the stress-strain curves and mechanical properties of CFS described in Table 1, can be found in Kankanamge and Mahendran (2011), for 250 and 450 steel grades, and in Ranawaka and Mahendran (2009), for 550 steel grade.
The actual stress-strain curves are obtained with a strain gauge at ambient temperature, as illustrated in Figure 5. Based on the stress-strain tests results, the material models adopted in the present investigation are bilinear and multilinear isotropic hardening, also graphically illustrated in Figure 5. The bilinear model is adopted with 0 y f σ = and t E 0 = . According to ANSYS APDL Theory Reference (ANSYS, 2016) the multilinear model is described by parts of the stress-strain curve, starting at the origin and defined by sets of positive stress and strain values with always positive slopes of the stress-strain curve.  It should be noticed that the author's present research numerical ultimate load is similar to the numerical model and laboratory tests presented by Heva (2009). The proposed numerical bilinear model has a relative difference of 6.4%, 6.0% and 4.0% with regard to tests of Heva (2009) and 3.2%, -1.1% and 1.2% differences with regard to with the numerical solution of Heva (2009), respectively for the models G250-1.95-1800, G450-1.90-1800 and G550-0.95-1800. Regarding the multilinear model, the relative difference is 6.7%, 1.1% and 0.0% when compared with tests of Heva (2009) and 3.4%, Latin American Journal of Solids and Structures, 2020, 17(9), e332 9/30 -5.6% and -2.7% differences with the numerical solution of Heva (2009), respectively for the models G250-1.95-1800, G450-1.90-1800 and G550-0.95-1800. For more convincing results, the graphs load versus displacement out of plane in the middle top flange at mid span (displacement position exemplified in Figure 6-b) for the 3 column specimens in Table 1 are shown in Figure 6 and Figure  7. These graphs show the displacement behavior of a specific point during the loading procedure.  It can be noted in Figure 6-b that the author's FEM has a similar equilibrium path as Heva's (2009) test, with a better fit with the multilinear model. On the other hand, in Figure 7-b, the laboratory tests are translated into the displacement direction, with a smooth equilibrium path. However, the bilinear numerical model examined in this research has a similar behavior as the numerical model examined by Heva (2009). Another important observation is concerned with the model G250-1.95-1800, shown in Figure 6-a. Note that comparing with the experiment, the present numerical model does not fit well with the experimental records. However, comparing the author's present study, both bilinear and multilinear, with the numerical model from Heva (2009), an almost exact fit of the equilibrium paths can be noticed.
Although, the multilinear material model has indicated better results for the column strength, it is more convenient to adopt the bilinear isotropic hardening model for the parametric study, because it makes the finite element analysis simpler, with satisfactory results. The main advantage of using a multilinear material model is to acquire more precise numerical results. However, the bilinear material model seemingly not influence significantly in the numerical results (equilibrium paths and column strength are somewhat close, with a maximum difference of 5% in the column strength), and it is a straight forward approach to perform in a parametric study. To sum up, the obtained results allow to conclude that the developed finite element model is validated for flexural-torsional buckling analysis.

Distortional buckling mode validation
The distortional buckling validation is performed using a single lipped channel column tested by Salles (2017), with the measured geometric dimensions referred to the out-of-section measurements illustrated in Figure 8-a. The test was performed at the COPPE Laboratory of Structures and Materials Professor Lobo Carneiro (LabEST). Updated information about the tested column can be found in . The lipped channel CFS column specimen was 2529 mm long, with material properties of 342 MPa of yield strength, 179.468 GPa of elastic modulus (quite low average Young modulus extracted from a set of standard tensile tests) and 9.845 GPa of tangent slope modulus ( t E ). As a result, the specimen reached its experimental column strength of 33.4 kN, while the finite element model in this research, lead to a numerical ultimate load of 35.3 kN. A maximum amplitude of 0.1 for the distortional buckling mode as initial geometric imperfection and a bilinear model for the strain-stress curve were considered. Figure 8-b shows a comparison of the equilibrium path of the web extremities at mid span, for the experimental test and the FEM solution. It can be noted in Figure 8-b that the equilibrium path of the FEM is close to that of the experimental results for small displacements. For larger displacements, the FEM model tends to follow a stiffer linear path until it collapses, while the laboratory test shows a long and smooth plateau before the collapse. Figure 9 shows the shape of the flange displacement out of plane measured by displacement transducer D4 along the column length, and Figure 10 for the displacement transducer D5, both for 5 loading steps. It can be observed that the same behavior of distortional buckling is clearly shown in the FEM and in the experimental test, with 3 half-waves of distortional modal shape. In addition, the buckling mode deformation increases with the load increment, until the column collapse. A fine agreement between numerical and experimental results is observed in Figure 9 and Figure 10.  The presented results indicate that the FEM model developed in the present investigation is able to accurately follow the equilibrium paths of both, global flexural torsional and distortional buckling modes of lipped channel coldformed columns.

PARAMETRIC ANALYSIS ON DG BUCKLING INTERACTION
First, a numerical parametric study is addressed to a lipped channel cold-formed steel column under true DG interaction, with different combinations of initial geometric imperfections (IGIs). Additionally, the study included a variation of the yield stress and column length, in order to understand the nature of behavior and strength of the DG interaction in a slenderness range, including Secondary-Global bifurcation DG interaction (SGI), True DG Interaction (TI) and Secondary-distortional bifurcation DG interaction (SDI). Figure 11 displays the critical load versus GD R λ of the LC 100x70x15x2.70 mm (out-to-out dimensions), performed by FStr and CUFSM (pure modes only), in a length range from 1000 to 3000 mm, i.e. . . GD 0 60 R 1 50 This ratio measures how far the critical loads of global and distortional modes are from each other. According to Martins, Camotim and Dinis (2018a), when (i) .
, the latter is the variable adopted by Martins, Camotim and Dinis (2018a).

Sensibility analysis of the IGIs' combination
This study consists of analyzing the structural behavior of a cold-formed lipped channel out-to-out section mm 2 , with different initial geometric imperfections (IGI) and yield strength. The goal of this study is to identify the column's strength and behavior IGIs sensibility using "impure" buckling mode combinations. The geometry is determined for a column experiencing strong DG interaction, for a length L characterizing True DG interaction (TI). The main reason of combining the buckling modes as IGI, emerged due to the difficulty in finding the DG buckling interaction in a simple elastic buckling analysis.
The combination of elastic buckling modes, as an alternative including the distortional and global buckling modes as IGI in the nonlinear analysis, arose because of the difficulty in obtaining a true distortional-global buckling shape through an elastic buckling analysis. Therefore, combining the global and distortional shapes as IGI, enabled the possible occurrence of DG coupled phenomenon. Figure 11 shows the critical load versus GD R λ for the proposed geometry and, as can be observed, the derivative of the curve for the first mode before L=1850 mm) and after this slenderness factor, changes drastically, as well as the critical modal shape. In this case, it is difficult to get a shape that presents a clear DG interaction. Because of this obstacle of getting a DG interaction shape mode, the first (Figure 12, for 0 θ =°) and second mode (Figure 12, for 90 θ =°) are combined. The buckling modal combination is performed using the first and the second mode shape, for a length of 1850mm. The buckling modes are illustrated in Figure 12  ), and details about mode composition and critical loads are given in Table 2. The critical load and the modal shapes are provided by the FStr Computer Application Lazzari (2020). The modal participation percentages are given by the Constrained Finite Strip Method, cFSM, with the help of CUFSM computer program Schafer (2020), confirming, in addition, the same accurate agreement between FStr and CUFSM results (see in Lazzari (2020)).

Figure 12
Analogy of modal shape initial imperfection combination in function of θ initial imperfection parameter (modal shapes amplified 10 times). It can be noted that the first mode (Table 2) is clearly a "global" one, which is more precisely defined as a flexuraltorsional mode. In this case, the DG buckling interaction is mainly focused in the distortional-flexural-torsional (D-FT) interaction of lipped channel cross sections. In order to study additional DG buckling interactions, such as the flexuraldistortional and torsional-distortional coupled phenomena, another geometric cross section is required that exhibits the flexural or torsional buckling mode, in transition with the distortional mode. Also, the first and second buckling modes (Table 2) are not pure modes, since they are obtained with the signature curve. Table 3 shows slenderness factors, ratios between slenderness factors and column strength computed with the direct strength method for both distortional and global buckling, for each value of the yield stress, i.e. corresponding to   Table 3: L λ , D λ and G λ are the local, distortional and global slenderness factor; is the ratio between G λ and D λ ; is the ratio between D λ and L λ ; / y cr P P is the ratio between the squash load over the minimum critical load (in this case of True DG buckling interaction, crD crG P P ≈ ); nG P and nD P are the DSM nominal axial strength given by DSM for distortional and global buckling, given by Eq. (1) and Eq. (2), respectively.
According to , for values of DL R λ higher than 1.20, the LD buckling interaction will not develop. Also, the proposed LD interaction equation by    Table 3,  , .
Since the displacement field is normalized, with a maximum displacement of 1.0, the coefficients 1 β and 2 β are the parameters of amplification of these shapes. These parameters are also defined in Eq. (5), where: ( ) C2 sin θ = ; θ is an angle for changing the buckling modal combination; 1 A and 2 A are the maximum amplitude for IGI, for the first and second mode respectively.
It should be noted that the parameters 1 β and 2 β (Eq. (5) A β = . Basically, θ is a single parameter that allows changing the IGIs and consequently the buckling modal combination shape. In order to illustrate the modes combination, Figure 12 shows the product 1 1 C A and 2 2 C A varying in a "trigonometric ellipse". This idea of combining the buckling modes in a "trigonometric ellipse" is inspired in a similar procedure proposed by Martins, Camotim and Dinis (2018a).
The maximum amplitude for the IGI, which is 1 A L 1000 = ⁄ and  Schafer and Peköz (1998a). With respect to these amplitudes, Table A.1 (Appendix) provides the percentage of each initial imperfection, in function of the θ initial imperfection parameter. This table shows the mode contribution changing from 100% to -100% for each mode, but with a total modulus summation of 100%. Basically, the negative percentage of the mode represents a buckling modal combination in opposite modal shape, i.e. multiplied by -1. Thereby, all the possible cases are combined. Figure 13-a shows clearly a cyclic behavior of the ultimate load and its behavior. Additionally, the symmetry of the behavior can be seen in the first quadrant. Basically, the behavior of the first and second quadrant are the same as the third and fourth. Another way to depict this cyclic behavior can be with the deformed columns at the ultimate load step displayed in Figure 13-b.  Note that in Figure 13-b the maximum displacement at ultimate load step changes based on the IGI shape. This shows the sensitivity of the column structural behavior depending on the initial geometric imperfection shape. Also, the cyclic behavior of the maximum displacement, from a global failure to a distortional one, is even clearer to notice. The maximum displacement cycle starts at mid span at the top flange and moves to the nearest distortional maximum halfwave deformation (around 40% of the maximum length). From that, the maximum displacement "jumps" to the bottom flange, with a distortional shape collapse. From that point, the maximum displacement moves to the mid span again, but at the bottom flange and, after that, the cycle repeats. Those modal shapes on the column strength step show that an analysis from is already sufficient to perform a complete DG interaction investigation. However, another cyclic and symmetric behavior for θ in first and second quadrant can be noticed. To conclude, a complete description of the relation between buckling modes combination (definition of the initial geometric imperfection) and the collapse mechanism, a more thorough analysis based on accurate observation of the equilibrium paths of the columns is required to avoid erroneous conclusions. In this case, the post-buckling equilibrium paths for  . Basically, this symmetric behavior evidences a "mirroring" behavior of the initial imperfection with θ from the first and second quadrant. This conclusion supports the parametric analysis of a large set of columns, to be conducted only with the initial imperfection for

IGI Combination with different yield stresses
The next step is to investigate columns with a higher yield strength. The same analysis conducted for columns with 345 MPa ( /  . This type of analysis is useful when examining the behavior of columns in a much larger elastic buckling range. The results for these analyses are shown in Lazzari (2020), on the other hand, the table with slenderness, ratios and DSM strength, is shown in Table 3 of the present investigation.
In order to visualize and compare the influence of each yield strength in the column structural behavior, the postbuckling equilibrium paths are presented. These results are shown in Figure 17 , Figure 17-c.2), the behavior is quite unique, with a smaller maximum displacement. This singular behavior is probably affected by a considerable absence of the global mode participation. Furthermore, for θ initial imperfection combination parameter close to o 90 , the ultimate load occurs on larger displacements, and for o 90 θ = (0G+1D), the ultimate load occurs in small displacements. Unfortunately, the finite element model with initial imperfection of 0G+1D did not converge for a higher squash load.
Evidently, all the elastic columns behave in global mode for larger displacements (except for columns with 0G+1D IGI, Figure 17-c.1). For columns with lower yield strength, on the other hand, the distortional mode develops at the beginning of the loading steps, when the distortional initial imperfection has a significant influence (see Figure 17- Moreover, it is noticed that for columns with a predominance of global initial geometric imperfection, the yielding of the material is symmetric and more concentrated at mid span. Whereas for columns with more distortional initial imperfection, the yielding process is more widely distributed. Comparing the strength of these columns with different IGI combination and different yield strength, it is possible to notice which initial imperfection influences more the column's load capacity. Figure 18 shows the column strength vs. squash load, normalized with a critical load in both axes.
As it was expected, and seen in Figure 18, for more elastic columns, the ultimate load increases. For

Analysis of DG buckling interaction nature
This study consists of investigating columns with different slenderness factors, to better understand the influence of the interaction nature on the column's strength. For this study the section LC 100x70x15x2.70 mm is considered with same geometry employed in the IGI combination analysis. However, the length of the column is modified, ranging from 1500 to 2200 mm, with an increment step of 50 mm (total of 15 lengths). This changing of column's length permits a modification of the global and distortional slenderness ratios and its contribution to the column behavior.
Basically, the investigation is performed for columns with three different natures: secondary-global bifurcation DG interaction (SGI), true DG interaction (TI) and secondary-distortional bifurcation DG interaction (SDI). In other to illustrate the DG buckling interaction investigation, Figure 11 shows the signature curve in contrast with the pure local, distortional, and global buckling mode curves.
Notice in Figure 11 that is associated with a length of 2200 mm. As can be observed, the range of SGI and SDI is not that extensive, because the analysis is limited to only one cross-section geometry. For a deeper understanding of the DG buckling interaction nature, it would be necessary to examine a wider variety of cross-section geometries (larger GD R λ range). The analysis performed in this research is focused on the TI nature, with fewer columns in the SGI and SDI regions.
Before starting the parametric analysis, it is important to define the IGI. This investigation is concerned with two different types of IGI: 100% Global (1G+0D) buckling mode IGI and 50% Distortional + 50% Global (0.5G+0.5D) buckling mode IGI. The reason of adopting those two initial imperfections is based on the results from the last study in the subsection 4.1 (Sensibility analysis of the IGIs' combination). Based on that study, it is concluded that for columns with a higher yield strength, 1G+0D IGI presents the lower ultimate column capacity. However, for columns with a lower yield strength, the 0.5G+0.5D IGI combination can give the lower column strength. These results revealed that a parametric study considering both IGI could indicate a difference in the column's strength, and consequently, some unconformity with the Direct Strength Method (DSM).
As a result of the IGI combination analysis, the following parametric study consists of using two different IGI (1G+0D and 0.5G+0.5D) for three different yielding stresses, 345 MPa, 508 MPa and 1016 MPa, i.e. applying 1G+0D with three yield stresses, and 0.5G+0.5D with also three yield stresses. For each of those yield stresses, the column strength capacity is reached, for 15 lengths, and compared with the actual distortional (Eq. (2)) and global direct strength equation (Eq. (1)). In addition, the strength is also compared to the proposed DG buckling interaction equations, proposed by Schafer (2002) (Eq. (3)) and by Martins, Camotim and Dinis (2018a) (Eq. (4)).
To sum up, 15 columns with different lengths, 3 types of yield stresses and 2 different IGI, totalizing 90 columns, were analyzed. The critical loads, slenderness factors, ultimate loads, and nominal axial strength for the 90 columns are reported in Table A.2 (appendix). With the purpose of a clear data visualization, Figure 19 illustrates graphically all the columns strength compared with DSM-based equations, for each yielding and IGI. Each graph shows a comparison with the global (Eq. (1)), distortional (Eq. (2)), distortional-global (Eq. (3)) and global-distortional (Eq. (4)) equations. More Latin American Journal of Solids and Structures, 2020, 17(9), e332 20/30 specifically, Figure 19-a shows the results for 345 MPa, Figure 19-b displays the results for 508 MPa and Figure 19-c reveal the results for 1016 MPa. Several observations can be made from Figure 19-a. First of all, the strength equation global-distortional, nGD P , (Eq. (4)) provides the most conservative approach, compared to the distortional-global (Eq. (3)) and global (Eq. (1)) equations, respectively nDG P and nG P . Secondly, the DSM nominal axial strength for global buckling (Eq. (1)), gives explicitly the lowest nominal strength, when compared to the DSM nominal strength for distortional buckling, nD P , (Eq. (2)), because all the / u nD P P data is below the / u nG P P data, and since the u P value is constant for both cases, it means that nD P is higher than nG P . Thirdly, in Figure 19-a.1 a small discontinuous gap can be observed for all the / u n P P ratios ( ). This gap refers to the half-wave switch from 4 to 5 waves, in the distortional buckling mode as IGI. It was realized that this half-wave changes interfered in the ultimate load, resulting in a slight reduction of the column strength. With respect to the type of initial imperfection, the 1G+0D initial imperfection has shown a smoother changing of the ultimate load, without discontinuous gaps. Finally, another interesting observation is the variation of the ratio / u nD P P , which displays a continuous decreasing related to the ratio GD R λ . This illustrates the fact that the column is approaching a global critical buckling region, since it is "leaving" the distortional critical buckling region.
There are a few remarks concerning the results for yield stress of 508 MPa, displayed in Figure 19-b. As reported by yielding of 345 MPa (Figure 19-a), the DSM-based equation global-distortional nGD P (Eq. (4)) provides the most conservative approach. Moreover, the DSM nominal axial strength for global buckling, nG P , is lower than the nD P equation, which according to the DSM, the global buckling mode is the one which governs the instability limit state of the structural element in all GD R λ range. Related to the IGI, the 1G+0D combination also has shown a softer changing of the ultimate load, without discontinuous gaps. As a final remark, the ratio / u nD P P has also demonstrated a similar behavior, as well as for yielding of 345 MPa. This result undoubtedly emphasizes the assumption that the columns reach a global critical buckling region. ). These gaps are related to the change of the half-wave number of the IGI for the distortional buckling mode. It was noticed that the first gap corresponds to the half-wave switch from 3 to 4 half-waves, and the second gap from 4 to 5 half-waves. With even numbers of half-waves, there is a slight increase in the ultimate load. The possible reason for this strength increase is due to the initial imperfection maximum displacement. When the half-wave number is odd, the maximum displacement from global and distortional modes coincides in the same cross-section, which is at mid span. Consequently, the maximum initial imperfection amplitude from each mode amplifies the maximum total displacement in the column, which end up resulting in a lower column strength capacity.
Even though the data results for yielding of 345 and 508 MPa are similar in some aspects, the results for yield stress of 1016 MPa, included in Figure 19-c, has shown a particular response. Since the columns with yield strength of 1016 MPa became further elastic, the behavior of the numerical ultimate load changed. It was noticed that the ratio / u nD P P decreases faster when GD R λ increases. Furthermore, the / u nG P P begins with 0.84 and speedily rises to 1.06 in both cases of IGI, showing that the ultimate load from the FEM approaches the DSM global equation quickly. On the other hand, the distortional-global nDG P (Eq. (3)) and global-distortional nGD P (Eq. (4)) interaction equations seem to stabilize on a plateau, with no abrupt changes. So far, the FEM data results are analyzed with the four possible equations separately. In order to compare the numerical ultimate load with the available standard DSM approach, Figure 20   It is noticeable that the larger the elastic behavior range, corresponding to higher yielding stress, the lower is the with 1G+0D IGI combination. Additionally, these columns with lower values of GD R λ and higher yielding, seem to diverge from the original DSM equation, due to a possible D-G buckling interaction in a secondary-global bifurcation D-G interaction (SGI) region. Lastly, comparing the two cases of initial imperfection combination, it is noticed that the mean, maximum and minimum value of / u nDSM P P ratio, are similar to each type of imperfection, while the 0.5G+0.5D initial imperfection combination presents the lower standard deviation and coefficient of variation.
Another approach to visualize the FEM results u P in Figure 20, is to display the corresponding data in the traditional column strength design curve ( / (1)), restricted to the global buckling mode. This graphical data is illustrated in Figure 21-b. This figure basically shows all the 90 columns strength over squash load ratio versus the global buckling slenderness factor, in addition with the global DSM equation and the Euler column curve for reference / 2 G 1 λ . Additionally, the same data is shown in Figure 21-a, but in comparison with the design approaches , , , nG nD nGD P P P and nDG P . In Figure 21-a it is shown how dispersed the DSM equations are from the ultimate load from FEM. Notice that the nD P procedure is clearly not corresponding to the column's strength. Furthermore, the , nG nGD P P , and nDG P appears to have more accurate approaches, even though the nGD P and nDG P procedures have shown to be quite conservative (points along the line -5% of n u P P = ).
Latin American Journal of Solids and Structures, 2020, 17(9), e332 22/30 In addition, notice in Figure 21-b that the difference of the 0.5G+0.5D and 1G+0D IGI in the global DSM equation is negligible. One important observation is the influence of the yielding stress. For yielding of 345 and 508 MPa, the / u y P P ratios are above the nG P equation, while for yielding of 1016 MPa, the / u y P P ratios are mostly below the nG P equation (see dashed circle in Figure 21-b). These results indicate that the global DSM equation handle well the columns under different types of DG buckling interaction nature (i.e. TI, SDI and SGI).  In general, Figure 22 has shown that the proposed equations by Schafer (2002) (Eq. (3)) and by Martins, Camotim and Dinis (2018a) (Eq. (4)) show a good agreement, i.e. lower standard deviation and coefficient of variation. Even though the standard deviation and coefficient of variation for both equations and both IGI combinations are relatively low, the mean, maximum, and minimum values have demonstrated that the proposed formulations are conservative. This means that the DSM nominal axial strength for the global buckling presented in the codes (Eq. (1)) is accurate enough for DG buckling interaction. These results are obviously limited by a group of the available column results, while the proposed DG buckling interaction equations handle quite conservatively the DSM-based nominal strength comparing to the ultimate load from the FEM. More specifically, note that the nDG P equation is less conservative than the nGD P , which can be justified by comparing the mean, standard deviation, coefficient of variation, minimum and maximum values ( Figure  22 and Table 4). Note that nDG P provides a lower mean, minimum and maximum value while the nDG P and nGD P equations give values of standard deviation and coefficient of variation relatively low (less than 0.044 and 3.6%, respectively).
In addition, it is noticed that the influence of the IGI type is not a relevant factor for the FEM performance of the columns. Even though the 1G+0D initial imperfection combination (Figure 22-b) has shown a sparser data for the nGD P equation, compared to the results from the nGD P equation with the 0.5G+0.5D initial imperfection (Figure 22-a), the influence on the numerical ultimate load is minimal. On the other hand, the nDG P equation revealed to have a minor impact on the data dispersion, when comparing both initial imperfections standard deviation and coefficient of variation.
Finally, comparing the nDG P , nGD P and nG P approaches, the conclusion can be drawn that the nG P procedure, already presented in the codes, handles the DG coupled phenomenon properly, for lipped channel columns tested in this study. The nominal axial strength equation for global buckling from the DSM was calibrated with many experiments of columns under global buckling. The coefficients calibrated for this equation, were revised by many authors, and it is a quite consolidated formulation (from low to intermediate slenderness columns, as reported by Dinis et al. (2020)).

FINAL REMARKS
The present research mainly provided an improvement in the comprehension of the distortional-global interaction buckling phenomenon, including both (i) its post-buckling behavior and (ii) the column ultimate load (strength). The investigation was aimed at CFS lipped channel columns, since these are CFS widely applied in steel construction and are prone to develop a DG buckling interaction.
Based on previous results and the numerical/computational support, a parametric study on distortional-global buckling interaction was carried out. The parametric analysis involved the FStr computational program for the definition of the initial geometric imperfections to be considered by the FEM non-linear analysis. Basically, only one cross-section of a lipped-channel was considered, with many combinations of yield stress y f and column's length L . The study was divided into two different investigations: (i) initial geometric imperfection combination and (ii) distortional-global buckling interaction nature, namely Secondary Distortional, Secondary Global and True Interaction, respectively SDI, SGI, and TI.

The importance of the imperfection combination with different yield stress
With respect to the initial geometric imperfection combination, the following points must be made: (i) Combining the initial geometric imperfection as first and second mode (which were classified as global and distortional mode), was an approach to include the distortional-global coupled phenomenon behavior into the post-buckling analysis, which was not clearly predictable by the elastic buckling analysis; (ii) Using an initial imperfection combination parameter θ , the initial geometric imperfection was modified, combining the global (flexural-torsional) with the distortional mode in different participation levels. Through this combination, the column strength analysis was performed and it was concluded that when θ belongs to the first and second quadrant (see Figure 12), a symmetric behavior was observed with the third and fourth quadrant. Indeed, the columns with an initial imperfection combination θ ≤ < , the predominance of the global mode with larger displacements was evident. On the other hand, columns with lower yield strength (less slender columns) developed a pronounced distortional buckling mode at the very beginning of the loading process, showing that the distortional initial imperfection significantly influences at the early steps of the loading; (vi) It has been shown that in columns with a lower yield strength, the distortional initial geometric imperfection (0G+1D) provides a lower ultimate load. For columns with intermediate slenderness, controlled in the present study by the variation of the yield stress, the combined initial geometric imperfection with 50% of global and 50% of distortional mode (0.5G+0.5D) gives the lower ultimate load. While for columns with very high yielding, the most detrimental ultimate load takes place with only global initial geometric imperfection (1G+0D). In conclusion, this analysis has shown that the initial geometric imperfection may affect the ultimate load and the column strength in different manners, according to the nature of DG buckling interaction.

The role of the nature of the DG buckling interaction in the behavior and strength of the columns
Regarding the investigation of the distortional-global buckling interaction nature, the following observations were reported: (i) It has been noticed that for high yield strength columns in the secondary-global bifurcation DG interaction region (which means lower values of GD R λ ) appears to diverge from the original direct strength method, which indicates a strong evidence of weak-to-moderate DG buckling interaction; (ii) Comparing the results with the two types of initial geometric imperfection considered (0.5G+0.5D and 1G+0D), it was concluded that the influence of the initial imperfection form was negligible in the ultimate load. Thereby, using an initial geometric imperfection with only a global (flexural-torsional) mode shape is sufficient to perform a parametric study on the distortional-flexural-torsional coupled phenomenon; (iii) The proposed equations by Schafer (2002) (Eq. (3)) and Martins, Camotim and Dinis (2018a) (Eq. (4)) addressed rather conservative results for the strength of the lipped channel columns described in the present research. More specifically, between the two approaches investigated here, the nGD P equation (Eq. (4)) has demonstrated to be more conservative than the nDG P (Eq. (3)) approach; (iv) The direct strength method equation addressed to global buckling can handle quite well the columns under the different DG nature (of the buckling interaction). However, the limits of the present investigation must be pointed, addressed only to lipped channel columns. Additional results of columns with different cross-section shapes and a wider range of GD R λ are needed to strengthen this assumption as a rule; (v) The authors believe that the new calibrated equations for global buckling, i.e. Dinis et al. (2020), might be a good approach for the case of DG buckling interaction for high slenderness columns. Additionally, for different types of end boundary condition, the DG buckling interaction approach, might be considered combining the new calibrated equations Latin American Journal of Solids and Structures, 2020, 17(9), e332 25/30 (i.e. Dinis et al. (2020) for global buckling and Liu et al. (2020) for distortional buckling) based on an idea proposed by Schafer (2002), i.e. replace the y P from the distortional equation (Liu et al. (2020)) by the nG P (Dinis et al. (2020)).
The conclusions obtained from the DG buckling interaction behavior are basically aligned with the ones found by and Martins, Camotim and Dinis (2018a),  and Martins et al. (2018d). However, more investigations on this topic are still needed, to better understands the phenomenon, since there is a lack of experimental tests of columns experiencing the distortional-global buckling interaction, in different types of nature.
To sum up, this work made possible a more thorough understanding of cold-formed lipped-channel columns under DG buckling interaction, with assistance of an elastic buckling analysis by the finite strip method, complimented by nonlinear FEM analysis. As initially proposed, the goal of this study has been achieved, with possible open topics that may be investigated in future research activities.