Eigenproblem Versus the Load-Carrying Capacity of Hybrid Thin-Walled Columns with Open Cross-Sections in the Elastic Range

The phenomena that occur during compression of hybrid thin-walled columns with open cross-sections in the elastic range are discussed. Nonlinear buckling problems were solved within Koiter’s approximation theory. A multimodal approach was assumed to investigate an effect of symmetrical and anti-symmetrical buckling modes on the ultimate load-carrying capacity. Detailed simulations were carried out for freely supported columns with a C-section and a top-hat type section of medium lengths. The columns under analysis were made of two layers of isotropic materials characterized by various mechanical properties. The results attained were verified with the finite element method (FEM). The boundary conditions applied in the FEM allowed us to confirm the eigensolutions obtained within Koiter’s theory with very high accuracy. Nonlinear solutions comply within these two approaches for low and medium overloads. To trace the correctness of the solutions, the Riks algorithm, which allows for investigating unsteady paths, was used in the FEM. The results for the ultimate load-carrying capacity obtained within the FEM are higher than those attained with Koiter’s approximation method, but the leap takes place on the identical equilibrium path as the one determined from Koiter’s theory.


Introduction
Experimental investigations [1] on the thin-walled columns made of modern materials with open cross-sections under compression have proved that the multi-mode buckling can be hazardous. When choosing the length of thin-walled elements, this effect has to be considered. These observations were confirmed in [2][3][4], where the buckling problem of laminated C-sections subject to pure bending in the flange plane [2,3] or in the web plane [4] was discussed. Attention was drawn to a rapid stability loss of beams for this type of loading. An unexpected phenomenon of stability loss of C-section beams of medium length subject to bending in the web's plane was observed, which was interpreted as a sudden influence of the secondary global distortional buckling mode characterized by a significantly higher value of the bifurcation stress than the primary mode [5]. The phenomenon was further investigated in [6,7] to define for which lengths this phenomenon is the most dangerous. The results of the experimental investigations and the calculations conducted with a semi-analytical method based on Koiter's asymptotic theory of thin-walled C-section beams in the web's plane of medium length were verified with the FEM calculations [8]. A significant effect of the higher global distortional mode on the load-carrying capacity of such beams was confirmed. With steel C-section and top-hat columns subject to compression, a catastrophic unexpected influence of the sec-ondary global distortional mode on the load-carrying capacity of the structure was indicated as well [3]. In [9], attention was paid to an effect of the distribution of inner forces, in particular, the global distortional longitudinal force, on an interaction of buckling modes and the load-carrying capacity of composite structures under bending. An influence of the secondary global mode on the load-carrying capacity of C-section beams of medium length under bending, made of a step-variable gradient material, was analyzed as well [10]. One of important conclusions from this work is such that a layer of glue connecting both layers of materials can be neglected in the modeling of such plates.
In stability of thin-walled structures, the determination of bifurcational values, which corresponds to a solution to the problem of eigenvalues, is one of fundamental issues. An eigenvector corresponds to eigenvalues. With structure stability, it is most often a buckling mode that is one of eigenvectors. It is also possible to determine inner sectional forces, including membrane forces. However, it should be remembered that eigenvalues are determined with an accuracy up to a certain constant. In [9], the key role of inner longitudinal forces corresponding to global distortional modes and of local distortional buckling modes was pointed out. These modes rarely correspond to the lowest bifurcational values of thin-walled structures [3,[6][7][8][9][10][11]. An effect of distortional modes on post-buckling equilibrium paths and the load-carrying capacity of structures was also investigated in [12][13][14][15][16][17][18].
Many bifurcational values, buckling modes, and components of membrane inner forces for various numbers of buckling half-waves along the longitudinal direction (or for the assumed length of the structure) can be determined for the linear eigenproblem, which is possible thanks to, among others, the semi-analytical method (SAM) based on Koiter's theory [6][7][8][9]11,19]. It should be mentioned that distributions of membrane inner forces are in balance. It is possible to indicate an eigenvalue and its corresponding eigenvectors, which can affect a rapid unexpected stability loss, because of the linear solution. This problem has become a starting point for the study presented here. Thin-walled hybrid columns in compression, with channel and top-hat cross-sections, were chosen as numerical examples to conduct the analysis. Column walls consist of two isotropic layers with step-variable material constants, i.e., an inner material A layer and an outer material B layer, as in [10]. Both layers are permanently integrated. In the numerical modeling, an effect of the layer connecting both the materials was not considered.
A multimodal approach to the phenomenon of stability loss of thin-walled columns in the elastic range was applied in the study. Numerical simulations were conducted with two methods, namely: the semi-analytical method (SAM) based on Koiter's theory [19][20][21] and the finite element method (FEM), which was to verify the results obtained with the first method. In the SAM, a plate model was used to analyze individual walls of the hybrid column employing the classical laminate plate theory (CLPT). Differential equations of the plate system under analysis were formulated with a variational calculus. A nonlinear solution to the stability loss problem was attained with a perturbation method, including additionally the full Green's strain tensor for thin-walled plates, the second Piola-Kirchhoff's stress tensor in Lagrange's description, and a numerical method of the transition matrix using Godunov's orthogonalization [8,19,21].

Formulation of the Problem
Multimodal buckling of hybrid columns with two types of cross-sections, namely: channel ( Figure 1a) and top-hat (Figure 1b), made of two isotropic materials of step-variable gradation of material properties, i.e., of material A and material B, in the elastic range and subject to compression, was considered. The cross-section does not alter along the whole length of the column. The thickness of each layer is the same and equal to 0.5 mm. The total thickness of the walls under analysis is 1 mm. The materials of the columns are made to satisfy Hooke's law. The material constants for the material A layer are: Young's modulus 71 GPa, Poisson's ratio 0.33, and for the material B layer 107 GPa and 0.34, respectively. The columns under study are simple supported at both ends.

Semi-Analytical Method (SAM)
The nonlinear problem of stability loss in thin-walled structures within Koiter's first-order approximation theory [19] was solved within the modified analytical-numerical method (ANM). In this approach, a precise transition matrix method is applied in contrast to the finite strip method. In the case under study, the analytical-numerical method should be additionally extended by the second-order approximation of Koiter's theory. In the proposed solution, second-order coefficients are estimated within the semi-analytical method (SAM) [8,21], in which it is postulated that values of these coefficients are defined based on the eigenproblem. It causes the SAM approach to yield a lower-bound estimation of post-buckling equilibrium paths. The SAM allows for solving the nonlinear problem of stability of thin-walled structures made of arbitrary materials, including hybrid ones.
The structures under consideration are prismatic thin-walled columns built of plates connected along longitudinal edges. In order to account for all buckling modes, from global through local to the coupled ones, a plate model (i.e., 2D) of thin-walled structures was employed. Columns were freely supported at their both ends. Exact geometrical relationships (i.e., full Green's strain tensor) were assumed for each plate [19,21]: where: u, v, w are displacements along the axes: x, y, z. According to Koiter's theory, the fields of displacement U and the fields of sectional forces N and M are distributed into a series regarding the dimensionless amplitude of deflection of the r-th buckling mode normalized to the thickness of the first component plate (i.e., ζr) [19][20][21]: where: λ is dimensionless load, components of the fields of displacement, and sectional forces N and M: U0, N0 in the pre-buckling state, Ur, Nr, Mr in the post-buckling state (first-order approximation), Urr, Nrr, Mrr second-order approximation. The index r = 1, …, J corresponds to the assumed number of the eigenvalue, whereas J is a number of eigenvalues considered in the analysis. It should be added that in the proposed formulation, a summation principle over the repeating index holds.
In thin-walled structures with the geometrical imperfections Ū (i.e., only linear initial geometrical imperfections corresponding to the r-th buckling mode are accounted for, i.e., ), the total potential energy can be expressed in the following form [19][20][21]: where the indices assume the values: p = 1, …, J, q = 1, …, J, r = 1, …, J. For the potential energy written in such a way (4), a system of equilibrium equations can be written shortly as: are the coefficients of the first-order (6) and second-order (7) nonlinear approximation, correspondingly [19][20][21]. The eigenvalues Ur are mutually orthogonal, which can be expressed as: σ0 l11(UI,UK) = 0, (I,K) = [1,J], I ≠ K, where J is a number of all buckling modes that are considered essential in the structure performance and were subject to the analysis. It allows the coefficients occurring in formula (6) to be written in the form: The above coefficients (Equations (6), (8) and (9)) depend only on eigenvalues. This significantly simplifies solving the nonlinear problem of stability loss in the first-order approximation. The first term in relationship (9) is described with the formula: r  y  q  y  r  x  q  x  r  y  q  y  r  x  q  x  r  y  q  y  r  x  q  xyp   y  r  y  q  y  r  y  q  y  r  y  q  yp  x  r  x  q  x  r  x  q  x  r  x  q (  [  )  ,  (   ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,  ,   ,  ,  ,  ,  ,  ,  , where n-number of component plates that form the column under analysis.
It follows from relationship (10) that the sectional forces corresponding to the eigenvalue denoted with the index p (i.e., ) exert a considerable influence on magnitudes of the three-index coefficients pqr a (Equation (9)). Because with these indices integrals of different signs are summed up, it should be remembered that the symmetrical modes regarding the axis of the cross-section can interact with an arbitrary number of these modes, whereas the anti-symmetrical modes interact only with anti-symmetrical pairs and with an arbitrary number of the symmetrical ones. With columns freely supported, solutions along the column length are developed into a sinus function, which causes subsequent eigenmodes to differ in the number of half-waves (further denoted as m) that appear along the column length. In the case of columns that are medium-long and long, global buckling occurs when one half-wave occurs, i.e., m = 1, and local buckling when the number of these half-waves is higher, i.e., m  1. To each buckling mode, for instance, denoted as r, a number of half-waves mr can be assigned. It follows from the conducted simulations that the three-index coefficients pqr a (Equation (6)) can be neglected when the number of half-waves differs considerably, i.e., 0   r q m m for the se-lected eigenvalues. Choosing the eigenvalues for analysis, one should follow the principle that these modes have an approximate number of half-waves that occurred along the column length. It has been found that these coefficients (Equation (6)) are equal to 0, when taking three eigenvalues to the analysis, a relationship holds that such that a sum of the numbers of their half-waves (i.e., mp + mq + mr) is an even number. Such an approach allows one to determine precisely the values of the coefficients apqr, according to the applied Byskov and Hutchinson's nonlinear theory [19,21]. The one-index coefficients r a (Equation (8)) depend on fractional values. Attention should be paid to the fact the coefficients apqr according to (6) introduced to equilibrium Equation (5) refer to subsequent values of bifurcational loads, whereas the coefficients pqr a are independent and determined from formula (9). The condition of zeroing of the Jacobian of system of Equation (5) corresponds to the boundary point (i.e., the ultimate load-carrying capacity or the maximal load) or to the point of secondary bifurcation. The dimensionless shortening of the column ends /min was determined as a function of the dimensional stress σ/σmin by differentiation of the expression for potential energy (4) with respect to σ/σmin [8,19,21] where: min-minimal shortening of the column corresponding to the minimal value of the bifurcational load σmin.

Finite Element Method (FEM)
The verifying simulations were conducted with the finite element method (FEM) employing a commercial package Abaqus [22]. A numerical model of the thin-walled structures under analysis was generated with four-node shell elements of eight degrees of freedom in each node (element type: S8R-eight-node doubly curved shell, reduced integration). The element size was determined through the convergence analysis to be equal to 2 mm. Boundary conditions ( Figure 2) were adapted to the requirements imposed by the SAM [23,24]: 1. In the center of the column, the displacement of the whole cross-section along the column length is equal to 0, i.e., uz = 0; 2. At both ends, it is assumed that the point of intersection of the web with the axis of symmetry of the cross-section (denoted as point 1 in Figure 2a,b) cannot displace along the y-axis, i.e., (uy)1 = 0. The web in segment 1-2 cannot displace along the x-axis, i.e., (ux)1-2 = 0. The arm in segment 2-3 displaces with respect to the y-axis in the same way as corner 2, i.e., (uy)2-3 = (uy)2. Reinforcement 3-4 displaces along the x-axis identical to corner 3 (Figure 2b), i.e., (ux)3-4 = (ux)3.
Another variant of boundary conditions, in compliance with the SAM requirements, was proposed by Szymczak and Kujawa [25][26][27][28]. Comparing the results of simulations of bifurcational loads and post-buckling behaviors for both variants of the boundary conditions, one could state that they agreed. During the verifying calculations, it was found that the determined eigenmodes and post-buckling equilibrium paths were identical in both the cases of the modified boundary conditions. It should be added that the nonlinear stability problem was solved both with the Riks algorithm and the Newton-Raphson algorithm. The results did not differ, thus the present work shows solely the results for the Riks algorithm, which allows for tracing not only stable equilibrium paths. To obtain the numerical convergence and a possibility to estimate the maximal load, the time-step was very fine (i.e., total time-step was 1, initial: 0.001, maximum: 0.01, and minimum: 1 × 10 −20 ). The SAM is much more efficient and allows for a much simpler analysis of the phenomena under investigation and their interpretation when compared to the FEM. Its main limitation lies in a necessity to select properly the eigenvalues, which affect con-siderably the post-buckling behavior of thin-walled structures. The second limitation of this method in comparison to the FEM is the fact that several eigenmodes are considered in the analysis (in general, fewer than 10).

Channel Column
In the first stage of the numerical simulations, it was assumed for the considered hybrid (i.e., materials A-B) channel of the geometrical dimensions of the cross-section shown in Figure 1a that the outer material B layer and the inner material A layer had the identical thickness of 0.5 mm each and, for the assumed material constants, the values of the bifurcational loads (i.e., eigenvalues) were determined in a broad range of variability of the column length from 25 up to 2000 mm, assuming that only one buckling half-wave (i.e., m = 1) formed along its length.
In Figure 3a, values of bifurcational stresses are given in MPa for the first three symmetrical buckling modes, whereas in Figure 3b the first three anti-symmetrical buckling modes as a function of the length of a single buckling half-wave are expressed in mm, correspondingly. The following numerical notations of the curves were assumed: 1. Primary modes corresponding to lowest bifurcational values when m = 1; 2. Secondary modes when m = 1; 3. Third modes, higher when m = 1, and the letter notations were used: S-denotes symmetrical modes with respect to the symmetry axis of the column cross-section, A-stands for anti-symmetrical modes.
The curve 1S in Figure 3a corresponds to the lowest symmetrical buckling modes when m = 1. The curve assumes the local minimum for the length of a single half-wave equal to 110 mm. For shorter lengths, it grows rapidly and reaches the local maximum when its length is 670 mm, then the curve descends gently for higher values of the length. The curve 2S has two local minima for the half-wave lengths: 70 and 670 mm, and it has one local maximum for the length of 380 mm. For the half-wave shorter than 1000 mm, the bifurcational stresses are lower than 1000 MPa. The local minimum for the curve 2S and the local maximum for the curve 1S occur for the same value of the half-wave equal to 380 mm. For the curve 3S, the local minimum occurs when the half-wave length is 390 mm, and the local maximum for the length of 210 mm. Bifurcational values for the curve 3S are one order of magnitude higher than for the curve 1S. In Figure 3b, the curve 1A (anti-symmetrical buckling mode) attains the local minimum for the half-wave length equal to 85 mm, the maximum for the length of 450 mm, then it descends for longer columns. The curve 2A has two local minima: for the half-wave lengths of 45 and 450 mm and two maxima for the lengths: 230 and 900 mm. In the half-wave's range length from 150 to 1100 mm, bifurcational values of the curve 2A do not exceed 1000 MPa. Bifurcational values for the curve 3A, within the whole range of length alternations, are in practice higher than 1000 MPa, and only for the length shorter than 100 mm are inconsiderably lower. Components of membrane forces of the first-order approximation were normalized by the normalization condition of the deflection amplitude of the i-th buckling mode for the first component plate thickness (i.e., ζi). In the further part of the study, according to (3b), only maximal absolute values of membrane components of the inner forces Nxi, Nyi, Nxyi, expressed in N/mm with accuracy up to the dimensionless amplitude of the i-th buckling mode ζi, will be considered. In Figure  As regards the plots of sectional forces (Figures 4 and 5), their maximal values and the corresponding lengths of buckling half-waves along the column lengths, in particular the longitudinal components Nxi, are of interest in practice as most significant in their interaction [9]. In Figure 4a, for the symmetrical modes, it is easy to see the dominant value of Nx3 ≈ 470 N/mm when the length is equal to 220 mm-the curve 3S (i.e., when the index i = 3). The remaining maximal values of Nxi for the curves 1S (i = 1) and 2S (i = 2) are over five times lower. For the curve 3S, the force Nx3 is higher than 50 N/mm in the length range from 190 to 500 mm. For the curve 1S, the values of Nx1 < 30 N/mm, whereas for 2S-Nx2 < 90 N/mm. In Figure 4b, maximal values of the forces Nyi are highest for very short channels for the length shorter than 30 mm, and they do not exceed the value of 10 N/mm for the first (i = 1), 15 N/mm for the second (i = 2), and 30 N/mm for the last buckling mode (i = 3), correspondingly. In Figure 4c  In Figure 5a, for anti-symmetrical modes, the maximal values of Nxi for the curves are: 1A-Nx1 < 25 N/mm, 2A-Nx2 < 220 N/mm, and 3A-Nx3 < 210 N/mm, correspondingly. The values of Nx2 and Nx3 are close to each other and twice as low as 3S. The plots of the curves 2A and 3A are by far more complex than for symmetrical modes. The curve 1A attains its maximum for the half-wave length of 600 mm. The curve 2A has its maximum for Nx2 ≈ 220 N/mm and the half-wave length of 270 mm and the local maximum for 60 N/mm and 1000 mm. It attains also two minima for the half-wave lengths: 50 and 670 mm. The curve 3A has two maxima for: Nx3 = 210 N/mm and the length 200 mm and Nx3 = 160 N/mm and the length of 470 mm, correspondingly. Its minima occur for Nx3 = 65 N/mm and the length 315 mm, and Nx3 = 15 N/mm and the length 50 mm. As can be easily noticed in Figure 5b, maximal values of the forces Ny2 are equal to for 1A-Ny1 < 10 N/mm, 2A-Ny2 < 25 N/mm, 3A-Ny3 < 35 N/mm. In Figure 5c, the maximal values of Nxy1 for the curve 1A are less than 3. The curve 2A has two maxima: Nxy2 = 30 N/mm and the length is 250 mm, and Nxy2 = 7 N/mm and the length is 1000 mm. Its minima occur when Nxy2 = 6 N/mm and the length is 50 mm, and Nxy2 = 2 N/mm and the length is 660 mm. By analogy, for the curve 3A, we have the maxima for: Nxy3 = 40N/mm and the length of 200 mm, and Nxy3 = 50 N/mm and the length 475 mm, and the minima for: Nxy3 = 21 N/mm and the length 285 mm, and Nxy3 = 7 N/mm and the length 50 mm. The maximal values of Nx2 and Nx3 (curves 2A and 3A in Figure 5a) occur for the same values as for Nx2 and Nx3 (curves 2A and 3A in Figure 5c).
The next stage consisted in a nonlinear analysis of interactive buckling, which allowed for determination of post-buckling equilibrium paths and, possibly, the load-carrying capacity (i.e., maximal loading) of the hybrid channel. Two total column lengths, for which absolute maximal values of the longitudinal forces Nxi were the highest (cf. Figures 4a and 5a), were selected for the analysis to observe strong interactions between the chosen eigenmodes. The following total lengths (further denoted as L) of the channel were assumed in the analysis: 210 and 250 mm. Table 1 lists values of bifurcational stresses for two assumed total lengths of the channel determined with two methods, namely: the SAM and the FEM. For the SAM, maximal absolute values of membrane sectional forces and several half-waves (in brackets) along the longitudinal direction that form along the length of the columns under analysis (denoted as the parameter m) and letter notations of the conditions on the axis of symmetry of the cross-section (i.e., S or A), are given. It should be added that for each variant denoted with the parameter m and the letter S or A, there are many eigenmodes, which are described: the lowest one is referred to as the primary eigenmode, whereas the subsequent ones, which are higher eigenmodes that differ from the primary one by deformations of the cross-section, can be found. The values of local bifurcational stresses are given for the buckling mode of two half-waves along the length (i.e., m = 2), as for this number of half-waves, the value of the bifurcational load is the lowest as well as for the mode of a single half-wave along the length (i.e., m = 1). It follows from a comparison of bifurcational loads that the SAM and FEM results are in high conformity. The maximal error for all the cases does not exceed 3%.
All the buckling modes from Table 1 for two lengths of the channel are shown in the next figures (Figures 6 and 7). For L = 210 mm, it is shown in Figure 6 that: (a) modes 1-3 are symmetrical modes, corresponding to subsequent modes with one half-wave along the column length (i.e., m = 1); (b) modes 4-6 are symmetrical modes, corresponding to subsequent modes with two half-waves along the column length (i.e., m = 2); (c) modes 7-9 are anti-symmetrical modes, corresponding to subsequent modes with one half-wave along the column length (i.e., m = 1); (d) modes 10-12 are anti-symmetrical modes, corresponding to subsequent modes with two half-waves along the column length (i.e., m = 2).
By analogy, all the modes from Table 1 when L = 250 mm are shown in Figure 7. Comparing the buckling modes for the total length 210 mm ( Figure 6), it can be easily noticed that for nearly all 12 buckling modes, the channel corner does not displace except for the symmetrical mode-mode 3-when one half-wave occurs along the column length (i.e., i = 3 and m = 1S) and for two anti-symmetrical modes-modes 8, 9-also, when one half-wave occurs along the column length (i.e., i = 8-9 and m = 1A). For these modes, the corners displace, and this is caused by an appearance of large longitudinal membrane forces Nxi when i = 3, 8, 9. These modes correspond to distortional buckling modes. Identical conclusions follow for the channel of the length L = 250 mm (Figure 7). Comparing the results for both the lengths and the maximal absolute values of the forces Nxi for buckling modes 3, 8, 9 (Table 1) to the corresponding buckling modes, one can see that there is a coupling between them, that is to say, on the basis of the displacement of the corner for the given mode, it can be stated that it is caused by a considerable inner longitudinal force. The reverse reasoning is correct as well. Table 1. Bifurcational stresses (i.e., σi in MPa) and maximal absolute values of membrane forces (i.e., Nximax/Nyimax/Nxyimax in N/mm) for selected buckling modes of channels with two lengths: L = 210 mm and L = 250 mm.

Mode Number
Denoted as i-Index  In the SAM, it is possible for the nonlinear problem to consider a strictly defined and finite number of buckling modes within an interaction of these modes. This allows for determining crucial modes that decide post-buckling equilibrium paths and the load-carrying capacity of the structure. When an interaction of J modes is considered, it means that the numerical model has such a number of degrees of freedom. In the FEM, it is not possible to determine which modes are to be considered in the numerical model. Both in the SAM and the FEM, a certain number of modes with inaccuracies and their amplitude should be assumed to solve a nonlinear problem. In this study, it was assumed for the SAM that the amplitude of the initial deflection was equal to: ζr *  = 0.2 (formula (5)), where r = 1, 2, …, J. The imperfection signs were chosen in the most dangerous way, i.e., to obtain the lowest value of the ultimate load-carrying capacity or the least slope equilibrium path. In the FEM, identical values of imperfections were assumed of course and their signs were chosen in the same way.
In Table 2, the dimensionless ultimate load-carrying capacities σS/σ4 obtained with the SAM and the FEM for both lengths of channels under consideration are presented. In the SAM, 3-, 4-and 6-modal approaches were considered, selecting various combinations of modes (Table 1-buckling modes [1][2][3][4][5][6][7][8][9][10][11][12]. For the first two, only symmetrical modes were considered, whereas when anti-symmetrical modes were also accounted for, a 6-modal approach was applied. This was caused by the necessity to consider an even number of anti-symmetrical modes in the interaction. Indices of buckling modes (i.e., i-index), which were considered in the SAM, and the dimensionless value of the ultimate load-carrying capacity referring to the minimal value of the buckling load, i.e., σs/σmin = σs/σ4 for the given total length of the column, were also given. In the case of the channel length L = 210 mm, for an interaction of three modes, the ultimate load-carrying capacity was not obtained. The post-buckling equilibrium path was an ascending curve. In a 4-modal approach, the dimensionless ultimate load-carrying capacity σs/σ4 is equal to 1.42 for an interaction of symmetrical modes. An increase in the number of modes up to six in the first case, i.e., when only anti-symmetrical modes with one half-wave along the length (m = 1) are accounted for, does not result in a change in the load-carrying capacity when compared to the 4-modal approach. When interactions of mode 10 (m = 2A) and modes 8 or 9 (m = 1A) are considered in the 6-modal approach, a significant lowering of the load-carrying capacity σs/σ4 ≈ 0.94 takes place. For modes 8 and 9 we have much higher values of longitudinal forces than for mode 7. In the FEM, 4 and 6 different imperfections corresponding to the modes of initial deflections considered in the SAM were assumed, and the lowest value of the dimensionless load-carrying capacity σs/σ4 equaled 1.49. In the FEM, it was not possible to confirm the SAM results obtained when anti-symmetrical modes were taken into consideration. The FEM load-carrying capacity corresponds to the SAM 4-modal approach. Similar differences in the evaluation of the load-carrying capacity with a coupled interaction of symmetrical modes in the SAM and an invisible effect of these modes on the load-carrying capacity in the FEM were found in [3].
For the channel of the total length L = 250 mm for each modal approach, the ultimate load-carrying capacity (i.e., σS/σ4) was attained. However, in a 3-modal approach, the load-carrying capacities differed 1.4 times from the other approaches. In a 4-modal approach for symmetrical modes, σs/σ4 = 1.39 was obtained similar to the 6-modal approach, in which two anti-symmetrical modes (modes 8, 9 when m = 1A, Table 1) are accounted for. For the 6-modal approach, a lower value of the load-carrying capacity was obtained at anti-symmetrical modes 10 (m = 2A) and 8 (m = 1A) than for modes 10 and 9. This difference was equal to 1.15 times. It results from the fact that Nx8>Nx9 and σ8 < σ9. From the FEM, σs/σ4 = 1.49 was obtained. Similarly, as for the length of L = 210 mm, an effect of the anti-symmetrical modes on the load-carrying capacity could not be confirmed in the FEM. The load-carrying capacities for two different column lengths from the FEM are identical. It should be remembered that the SAM yields a lower-bound evaluation of the load-carrying capacity and foresees the ultimate load-carrying capacity with accuracy below 10% regarding the FEM for symmetrical modes. Because of the time involved in numerical calculations, this outcome is very satisfactory. Considering the results from [3], the authors think that further investigations on an influence of anti-symmetrical modes on interactive buckling and load-carrying capacity should be carried out. In Figure 8, post-buckling equilibrium paths in the system: dimensionless compressive stress-dimensionless column shortening (i.e., (σs/σmin) − (Δs/Δmin), where σmin corresponds to the lowest buckling stress, whereas Δmin-denotes the shortening corresponding to this stress) obtained within the SAM and the FEM for the channel of the total length L = 210 mm are presented. An analogous plot for the length L = 250 mm is drawn in Figure 9.  In Figure 8, the curve representing the SAM 4-modal approach reaches its maximal value for the dimensionless shortening equal to 4.5, whereas for the FEM it is 2.0 at approximate values of the dimensionless ultimate load-carrying capacity σs/σmin. The FEM model is more rigid than the SAM model. It should be remembered that in the SAM only up to six degrees of freedom were considered, whereas the FEM model had four orders of magnitude more. A nonlinear analysis was conducted only for the elastic range. The dimensionless ultimate load-carrying capacity attained with the SAM equaled 0.9, and for the FEM, after a jump on the identical path as in the SAM simulations, it was approx. 1.2 at the point the calculations were disrupted by a loss of convergence. For the SAM 6-modal approach, according to the conclusions from [3], numerical calculations should be continued from the authors' viewpoint.
Similar conclusions can be drawn for the channel of the length equal to L = 250 mm ( Figure 9). In the case of the FEM simulations, the ultimate value occurs when the dimensionless shortening is 2.0. Thus, a slight alternation in the length (about 1.2 times) causes the influence of global modes to be much more considerable, as the lowest bifurcational values for local modes are practically the same for both lengths. The FEM post-buckling path falls down after it reaches the value of the ultimate load-carrying capacity and the curve clearly reaches the minimum for σ/σ4 ≈ 1.2, and then the FEM code disrupts the calculations. In the authors' opinion, new buckling modes appear at the points the FEM calculations are disrupted (i.e., at σ/σ4 ≈ 1.1) for both lengths of the channels under analysis. This causes the SAM and FEM results to come closer. They should be subject to further in-depth analysis.
The next figure shows a collapse mode that corresponds to the ultimate load-carrying capacity σs/σ4 of channels of the lengths: L = 210 mm and L = 250 mm (Figure 10). The FEM collapse mode for both the column lengths are identical. Two buckling half-waves appear along the channel length. In the half-length, a visible deviation from the straight line can be seen, which corresponds to an effect of the distortional buckling mode.

Top-Hat Column
The second thin-walled column subject to consideration is a hybrid column with a top-hat cross-section of the geometrical dimensions of the cross-section shown in Figure  1b. An addition of external boundary reinforcements causes stiffening of the structure compared to the channel and thus, enlarges bifurcational loads. Identically as in the channel previously discussed, it is assumed that the outer layer is made of material B, and the inner layer is made of material A. The thicknesses and the material constants are the same as in the previously discussed column.
In Figure 11, a variability in bifurcational loads as a function of the length of one buckling half-wave for symmetrical modes (Figure 11a) and for anti-symmetrical modes (Figure 11b) between lengths from 25 to 1000 mm, which is most interesting from the practical point of view, is shown. Like in Figure 3, the following notations of curves are assumed: 1-primary mode, 2-secondary mode, 3-third mode, S-symmetrical mode, A-anti-symmetrical mode. The curve 1S in Figure 11a  Next, in Figures 12 and 13, maximal absolute values of membrane sectional forces for symmetrical and anti-symmetrical buckling modes, respectively, are presented. As discussed in Section 3.1 for hybrid channels, the most important are maximal values of the longitudinal forces Nxi when i = 1,2,3. For the first symmetrical mode of the longitudinal force (curve 1S in Figure 12a), the maximal value is attained when the length of a single buckling half-wave is equal to 220 mm, and it exceeds the value of 100 N/mm when the half-wave length changes from 185 to 250 mm. Respectively, the curve 2S reaches its maximum when the half-wave length is 125 mm and exceeds 100 N/mm in the range from 60 to 225 mm, whereas the curve 3S attains the maximum when the half-wave lengths equal 80 and 330 mm and exceeds the value of 100 N/mm in the ranges: from 60 to 160 mm and from 260 to 440 mm. Membrane transverse forces (Figure 12b Maximal values of cross-sectional forces for anti-symmetrical modes are presented in Figure 13. Maximal longitudinal forces are (Figure 13a), correspondingly, for the curves: 1A has a maximal value of 300 N/mm and for lengths 60-250 mm are higher than 100 N/mm; the curve 2A has two maxima: 350 N/mm for the range 50-170 mm and 160 N/mm for 250-480 mm; the curve 3A has two maxima as well: 600 and 200 N/mm for the range 60-300 mm. As can be easily noticed in Figure 13b (Figures 12 and 13) to C-columns (Figures 4 and 5), maximal values of inner forces occur for shorter lengths of structures under compression, which is caused by boundary reinforcements of TH-columns in relation to C-columns.
Considering maximal absolute values of longitudinal forces in evaluation of their influence on the load-carrying capacity of hybrid TH-channels, the following lengths were selected: 80, 130, 225, and 330 mm. In Table 3, values of bifurcational stresses and maximal values of membrane sectional forces for columns with two total lengths of 80 and 130 mm, corresponding to three symmetrical and anti-symmetrical buckling modes regarding the axis of symmetry, are presented. The ratios of values of subsequent bifurcational loads to the lowest value of the bifurcational load σi/σ4 are included. For the TH-column of the total length 80 mm (Table 4), the lowest value of bifurcational load occurs for the symmetrical mode (S) with the number of half-waves along the longitudinal direction equal to m = 1S (mode 1 in Table 3). As one can easily see, values of longitudinal inner forces are significantly higher than 100 N/mm, except for the lowest value of bifurcational load σ1.  Table 4. Dimensionless overload coefficient σ/σmin as a function of the dimensionless shortening Δ/Δmin = 5.0 for TH-columns (ANM results).

Total Length L in mm J-Mode Approach
Mode Numbers (i-Index in Table 3 Buckling modes are depicted in Figure 14. For the two lowest symmetrical and anti-symmetrical buckling modes, the corners of the TH-column do not displace. For the remaining four modes, the corner joining the web with the flanges does not displace, whereas the second corner joining the web with the boundary reinforcements is subject to displacement, which corresponds to the distortional-local mode. Next, post-buckling equilibrium paths were determined, employing the SAM only and covering from a one-modal approach (i.e., J = 1-uncoupled buckling) up to a 4-modal approach (J = 4). In the calculations, identical to the channels, the amplitude of imperfections (5) was assumed: ζr *  = 0.2 where r = 1, 2, …, J and various combinations of their signs. Here, the post-buckling equilibrium path does not attain the maximal value. Therefore, in Table 5, for individual modal approaches for the assumed dimensionless value of shortening Δ/Δmin = Δ/Δ1 = 5.0, dimensionless values of overloads: σ/σmin = σ/σ1 are given. As seen, an interaction between the buckling modes under investigation is absent in practice. Attention should be drawn to the fact that for the 4-modal approach, a pair of anti-symmetrical modes was considered. While analyzing the TH-channel of the total length equal to 130 mm, it was found that the lowest value of the buckling load corresponded to the symmetrical mode of two half-waves along the column length (i.e., m = 2S-mode 7 in Table  3). In this table, values of bifurcational loads and maximal values of membrane forces for two types of modes and two values of half-wave numbers m = 1, 2 are listed. The ratios σi/σmin = σi/σ7 are given as well. Values of longitudinal sectional forces are higher than 100 N/mm in 10 cases, except for two lowest values of bifurcational loads for m = 1 and m = 2. Buckling modes are presented in Figure 15. As can be easily seen, the buckling modes: symmetrical modes 3 and 5 ( Figure 15a) and anti-symmetrical modes 2, 4, 6 for m = 1 (Table 3) Table 3 (a), and anti-symmetrical modes 2, 4, 6 in Table 3   The remaining modes are the local ones, for which the corners of the TH-column cross-section do not displace practically. Like the TH-column of the length 80 mm, post-buckling paths do not attain maximal values (i.e., ultimate load-carrying capacity). Thus, in Table 4, for the value of the dimensionless shortening Δ/Δmin = Δ/Δ7 = 5.0, the value of the load is σ/σmin = σ/σ7. For a one-modal approach (J = 1), the overload is equal to 3.67, whereas for 2-to 4-modal approaches, the overload decreases approx. 1.25 times. When only an interaction of symmetrical modes is accounted for, the overload decreases to 3.08. When two additional anti-symmetrical modes are included in the 4-modal approach, the overload diminishes up to 2.94.
In further considerations, TH-columns of total length equal to 225 mm were dealt with. In Table 5 Table 5 (i.e., σ19). For each number of half-waves, the values of sectional forces for three symmetrical and three anti-symmetrical buckling modes were analyzed.
The maximal bifurcational loads were attained for mode 6 when m = 1 (Table 5). In three cases (modes 11, 12, and 24), the FEM modes which would correspond to the ANM modes were not found. The most difficult problems consisted in assigning the FEM modes for m = 2, where the error in bifurcational loads reaches 28%. A very large number of degrees of freedom in the FEM causes these modes to differ significantly from the ANM "pure" modal modes. With an exception for the cases when the difference in bifurcational loads is over 10% (four cases) and in three cases lacking recording, the differences in bifurcational loads between the ANM and the FEM are lower than 7%. Among 24 modes under consideration, only for three modes (modes 7, 13, 19 in Table 5) are the values of maximal longitudinal forces lower than 50N/mm, and for 18 modes they are higher than 100 N/mm.
In Figure 16a Table 5) and anti-symmetrical modes (modes 2, 4, 6, Table 5) are distortional modes, because the corners displace along with the boundary reinforcements. Anti-symmetrical mode 2 is the only one where the corner displaces between the web and the flange. The distortional modes for m = 2 are modes 8-12, and for m = 3 only anti-symmetrical mode 18. In Table 6, values of the dimensionless ultimate load-carrying capacity referring to the lowest value of the bifurcational load σ19 (i.e., σS/σ19) for both the SAM and the FEM are listed. The value of the amplitude of the initial deflection was assumed as ζr *  = 0.2 (when r = 1, 2, …, J). An interaction of the primary mode when m = 1 (mode 1 in Table 5) and subsequent higher modes (i.e., modes 7-24 in Table 5, when m > 1) was considered. This is caused because modes of a various number of half-waves m higher than 1 forming along the length of the column under consideration (i.e., for m > 1) do not interact in practice. In the interaction of selected modes 1-6 when m = 1 (Table 5) and modes 19-24 when m = 4 (Table 5), 2-4-modal approaches were considered. The attained values of the dimensionless load-carrying capacity for these approaches differ slightly, and the lowest value of the dimensionless load-carrying capacity equal to 1.98 (Table 6) was obtained for the 4-modal approach. In Table 6 values of the dimensionless load-carrying capacity for an interaction of selected modes 1-6 when m = 1 (Table 5) and modes 13-18 when m = 3 (Table 5) are given. The lowest value of the bifurcational load σ13 when m = 3 is slightly higher than σ19. Here, 2-6-modal approaches were considered. With an increase in the number of degrees of freedom (i.e., for higher values of J), to which nonzero imperfections correspond as well, the dimensionless load-carrying capacity decreases less than 0.5%. For an interaction of selected modes 1-6 (m = 1, Table 5) and modes 7-12 (m = 2, Table 5), the lowest bifurcational value occurs for mode 7, for which the value of the bifurcational load referring to the lowest value of the bifurcational load σ7/σ19 equals 1.37. As already mentioned for modes 7-12 (m = 2, Table 5), in six modes under consideration, five are distortional modes, for which longitudinal inner forces are significant. In coupled interactions for modes 1-6 (m = 1, Table 5) and modes 7-12 (m = 2, Table 5), 2-and 4-modal approaches were considered. The obtained value of the dimensionless load-carrying capacity equals σS/σ19 = 2.01. It should be taken into account that the minimal value of the bifurcational load for mode 17 (m = 2, Table 5) is 1.37 higher than for mode 19 (m = 4, Table 5), and the value of the dimensionless load-carrying capacity determined from the equations of equilibrium was recalculated to the lowest bifurcational value σ19. Thus, distortional-local interactions when m = 2 can be significant, as already noticed in [9]. In Table 6 results of the dimensionless ultimate load-carrying capacity for the SAM are compared to the FEM. Also, in this case, the SAM outcomes are lower and differ by 1.34 times. The lowest post-buckling equilibrium paths are shown in Figure 17, whereas Figure 18 presents a collapse mode corresponding to the ultimate load-carrying capacity according to the FEM. For the column of the total length equal to 225 mm, post-buckling equilibrium paths up to the value σ/σ19 ≈ 2 overlap ( Figure 17). As seen, the buckling mode for σS/σ19 corresponds to the coupled buckling mode 1 and 19, when m = 1 and m = 4 ( Figure 18). Table 6. Values of the dimensionless ultimate load-carrying capacity for the TH-channel of total length 225 mm (the semi-analytical method (SAM) and the finite element method (FEM) results).

SAM FEM J-Mode Approach
Mode Numbers (i-Index in   Finally, TH-columns of total length 330 mm were subject to analysis. Taking into consideration the conclusions for the columns of total length 225 mm, analogous to the shorter lengths of the columns, three bifurcational values for the mode with a single half-wave occurring along the column length (i.e., m = 1) and the lowest local modes with five half-waves (i.e., m = 5) for symmetrical and anti-symmetrical conditions of the axis of symmetry along with their corresponding values of inner forces, are given in Table 7 for the ANM. The values of maximal longitudinal inner forces are higher than 30 N/mm for 11 of the 12 modes under consideration, and they are higher than 100 N/mm for seven modes as well. Also, the FEM results were quoted. For the two highest bifurcational loads corresponding to the anti-symmetrical modes, the buckling mode corresponding to the SAM was not found in the FEM. The maximal differences between the values obtained with the ANM and the FEM are lower than 3.5%.
The buckling modes for the ANM are shown in Figure 19. As can be easily noticed, all modes 1-6 for one half-wave along the column length (i.e., m = 1, Table 7) are distortional modes, and for modes 7-12 (m = 2, Table 7), they are symmetrical only. Attention should be drawn to modes 4, 5 (m = 1, Table 7) and mode 11 (m = 5, Table 7), for which displacements of the corner joining the web with the flange are high. (c) (d) Figure 19. Buckling modes for the TH-channel of total length 330 mm (Table 7): symmetrical modes 1, 3, 5 (a), modes 7, 9, 11 (b), and anti-symmetrical modes 2, 4, 6 (c), modes 8, 10, 12 (d).  Table 8 values of the dimensionless load-carrying capacity σS/σ7 for the SAM are given on the assumption of 2-6-modal approaches. The values of the load-carrying capacity are practically the same and equal to σS/σ7 = 1.57, whereas for the FEM they are equal to 2.0. Identical to shorter lengths of the TH-column, the dimensionless load-carrying capacity in the SAM is lower than in the FEM. In the FEM, modes that stabilize the post-bucking path contribute. The FEM model is more prone to this than the SAM model. Figure 20 shows post-buckling equilibrium paths for the SAM and the FEM. The next figure depicts a collapse mode for the load-carrying capacity according to the FEM ( Figure 21).  Figure 20). In Figure 20, equilibrium paths obtained from the SAM and the FEM are presented. The shape of the plots is close to the value of σ/σ7 ≈ 1.6. For higher values of the overload, the FEM post-buckling path grows up to the load-carrying capacity σS/σ7 = 2.04. The most interesting is the part of the curve obtained with the Riks method, for the curve descending to the load-carrying capacity value determined in the SAM (i.e., for σS/σ7 ≈ 1.6), then due to a loss of convergence, the FEM Abacus package stopped the computations. As opposed to the post-buckling curves in Figures 8 and 9 (FEM), the descending part of the curve (Figure 20) lies outside the path regarding the ascending part. Such a curve shape is very unusual for the descending part of the path, and it requires further detailed numerical analysis. However, this does not fall within the scope of the present study.
In Figure 21a, a buckling mode for the ultimate load-carrying capacity σS/σ7 = 2.04 is presented. Along the TH-channel length, one can see four buckling half-waves for the web, and two half-waves are very clearly visible. There is only one buckling half-wave on the flanges with the boundary reinforcements. Figure 21b shows a mode for the point when the computations were disrupted (i.e., for σ/σ7 ≈ 1.6). Here, one buckling half-wave occurs along the column length, both in the web and in the flanges. It indicates that a new buckling mode, which causes the column to be relieved and which leads to its destructions, appears.

Conclusions
A problem of rapid, unexpected exceeding of the load-carrying capacity was discussed using the example of hybrid thin-walled columns with open cross-sections under compression. The numerical simulations, conducted with the semi-analytical SAM based on Koiter's theory, allowed for solving the linear eigenproblem and determining bifurcational loads, their corresponding eigenvalues, and components of membrane sectional forces. Distributing membrane sectional forces that are in balance based on the linear solution, it was possible to indicate the eigenvalue and its corresponding eigenvectors, which can affect rapid, unexpected exceeding of the load-carrying capacity of the hybrid thin-walled structure under discussion. The results were verified with the FEM in the linear and nonlinear range.
The detailed simulations were conducted for simple, supported, hybrid C-channels and TH-channels of the assumed lengths. The walls of the columns under investigation were made of two isotropic materials with a step-variable gradation of mechanical properties in the elastic range. While selecting the column lengths, the authors followed the principle that the determined maximal absolute values of longitudinal forces should attain their maxima in the solutions to the eigenproblem, which causes strong interactions between the selected eigenmodes. Summing up the results obtained for the assumed lengths and the determined maximal absolute values of membrane longitudinal forces with their corresponding buckling modes, one can see that there was a coupling between them based on the displacement of the corner for a given mode caused by a strong longitudinal sectional force. The reverse reasoning is correct as well.
In the SAM, it was possible to consider a precisely defined and finite number of buckling modes, taking into account in the interaction of these modes. This allowed one to determine the key modes that decided post-buckling equilibrium paths and the load-carrying capacity of the structure. In the FEM, it was practically impossible to choose which modes were to be considered. Including two anti-symmetrical modes in the multimodal SAM approach caused a significant decrease in the load-carrying capacity, and the stronger it was, the higher were the values of longitudinal sectional forces observed for the selected modes. The most crucial interactions of anti-symmetrical modes took place for the distortional global and local mode. In the FEM, it was not possible to confirm the results obtained when the anti-symmetrical modes were analyzed. The ultimate load-carrying capacity determined with the FEM corresponded to the SAM only if symmetrical modes were considered. Similar differences in the evaluation of the load-carrying capacity when a coupled interaction of anti-symmetrical buckling modes in the SAM multimodal approach was included, and an invisible effect of these modes on the load-carrying capacity in the FEM was found in [3].
In the FEM numerical simulations, two algorithms were employed to solve the nonlinear problem of stability loss, namely, the Riks algorithm and the Newton-Raphson algorithm. The results for stable equilibrium paths in both algorithms were identical. However, the Riks algorithm allowed us to catch the effect of a jump between two stable equilibrium paths. The post-buckling equilibrium path determined with the Riks method under low overloads was the same as that obtained with the SAM. At higher overloads, it overestimated the ultimate load-carrying capacity, then it jumped onto another equilibrium path, coming closer to the SAM path, when only symmetrical buckling modes were accounted for in the interaction. In the authors' opinion, a jump between equilibrium paths results from the fact that new buckling modes appear. It causes the outcomes attained within both methods to become closer. These observations should be followed by further in-depth analysis.