The Influence of the Internal Forces of the Buckling Modes on the Load-Carrying Capacity of Composite Medium-Length Beams under Bending

The distribution of the internal forces corresponding to the individual buckling modes of lip-channel (LC) beams is investigated using the Semi Analytical Method (SAM) and the Finite Element Method (FEM). Channel section beams made of 8-layered GFRP (Glass Fiber Reinforced Polymer) laminate with three different layer arrangements were considered. The effect of the internal forces on the non-linear first-order coefficients corresponding to the interactive buckling was also studied. Moreover, distributions of the internal forces corresponded to the loading, leading to structure failure for which the load-carrying capacity was determined. The results indicated a high influence of the Nx internal force component on the buckling loads and load-carrying capacity of the LC-beams.


Introduction
Thin-walled beam structures are widely used as primary structural components in engineering practice. I-section and C-section beams are often utilised as basic elements in many applications. In most thin-walled beams, the load-carrying capacity of these structures is limited not only by their strength but also by their stability. Thus, the problem of loss of stability is a subject of interest among researchers.
The phenomenon of unexpected loss of the load-carrying capacity of bent C-beams was observed during experimental research [1,2]. Kolakowski et al. [3] interpreted this behaviour as the effect of the secondary global distortional-lateral mode. Studies were performed on steel beams with channel and lip-channel cross-sections, and various beam lengths. To identify this phenomenon, the semi analytical method (SAM) within the second order of non-linear approximation was used. This method was found to be very useful for the interpretation of the interaction of different buckling modes in the full load range. The proposed method facilitates the understanding of the phenomena occurring during coupled buckling. In [4], the majority of conclusions from [3] were confirmed for steel LC-beams. The beam lengths for which the interaction of the secondary global buckling mode affects load-carrying capacity the most were also determined.
In the authors' view, there is a lack of works devoted to the influence of the secondary global distortional-lateral buckling mode on the interactive buckling and load-carrying capacity of LC-beams. A few studies can be found in [5] in which the authors used FEM to verify the results of SAM calculations. High agreement of results was obtained. In [6,7], numerical and experimental studies were carried out on steel channel section beams to analyse the distortional-global interaction buckling. and discussed. The most important advantage of this method is the possibility of describitng the thin-walled structures in the full range of behaviours from global to local stability [26,27].
In this paper, the effect of the interactive buckling of the distortional-global and local-distortional buckling modes on the post-critical equilibrium paths and load-carrying capacity was taken into account. The dominant character of the internal force components was indicated using SAM. In this analysis, the change of the internal force distribution in the post-critical regime is studied. The main purpose of the work is to indicate which components of internal forces and which buckling mode have a dominant influence on interactive buckling. According to the authors, this is a novel approach in the literature on the subject. The comparison was performed for the load-carrying capacity and failure (within the Hashin failure criterion [28] and the elastic limit load-carrying capacity of the structure) of composite LC beams. The studies were performed for three LC medium-length beams, which differed in laminate stacking sequence. The length of the beams was chosen based on the results of previous works [3,4], to obtain the strongest buckling interaction.

Formulation of the Problem
Thin-walled prismatic composite beams were taken into consideration. It was assumed that the beams were composed of plates, connected along longitudinal edges, subjected to uniform major-axis bending moment and simply supported at both ends [2][3][4][5]25,27].
To take into account all buckling modes (global, local and coupled buckling), a plate model (i.e., 2D) of thin-walled structures was used. The exact geometrical relationships (i.e., full Green's strain tensor) were adopted for each plate to analyse both out-of-plane and in-plane bending of the i-th plate [26,27]: and: κ xi = −w i,xx κ yi = −w i,yy κ xyi = −2w i,xy (2) where u i , v i , w i -the components of the displacement vector of the i-th plate in the x i , y i , z i axis direction, respectively, and the x i -y i plane overlaps with the mid-plane before it buckles.
To solve the non-linear problem of structure stability, Koiter's theory [26,27] was implemented. The displacement fields U and the sectional force fields N were extended to the power series for the dimensionless amplitude of the r-th mode deflection ζ r : U ≡ (u, v, w) = λU 0 + ζ r U r + ζ 2 r U rr + . . . N ≡ (N x , N y , N xy ) = λN 0 + ζ r N r + ζ 2 r N rr + . . .
where: λ-load factor, U 0 , N 0 -the pre-buckling (i.e., unbending) fields, U r , N r -the first-order non-linear fields, U rr , and N rr -the second-order non-linear fields of the displacement and the sectional force, respectively. The range of indices is [1,J], where J is the number of interacting modes. It is assumed that the summation is over the repeated indices.
In thin-walled structures with initial geometric imperfections,Ū (only the linear initial geometric imperfections corresponding to the shape of the r-th buckling mode, i.e.,Ū = ζ r * U r are taken into account), the total potential energy, can be described by the following equation: Next, the equilibrium equations corresponding to Equation (4) have the form: 1 − M M r ζ r + a pqr ζ p ζ q + b rrrr ζ 3 r = M M r ζ * r r = 1, . . . , J where: M-magnitude of the applied bending moment; M r -buckling moment of the r-th buckling mode; ζ r -the dimensionless amplitude of the r-th buckling mode; ζ r *-the dimensionless amplitude of the initial deflection of the r-th buckling mode.
The buckling modes U i are mutually orthogonal, and can be described in the following form: σ 0 l 11 (U I ,U K ) = 0, (I,K) = [1,J], I K where J is the number of all relevant buckling modes that are considered to be crucial in the structural response. The coefficientsā 0 ,ā r ,ā pqr and b rrrr can be calculated from the equations described in the literature [26,27].
The following notations are introduced in Equation (5): a pqr = a pqr /a r b rrrr = b rrrr /a r In the semi-analytical method (SAM), approximate values of the b rrrr coefficients (5) are determined based on the linear buckling problem. This approach makes it possible to estimate precisely the values of the a pqr coefficients (5), according to the applied non-linear Byskov and Hutchinson theory [26,27]. It should be mentioned that a pqr coefficients (5) are the key factors affecting interaction. For a more detailed analysis see Appendix A.
A relative angle of rotation of the girder ends was determined as a function of the M/M min , through differentiation of the expression of potential energy (4) for M/M min [26,27]: M min M r a r ζ r (0.5ζ r + ζ * r ) where: α min -the minimum buckling angle of rotation of the beam subjected to pure bending, corresponding to the minimum value of the buckling moment M min . For the initial geometrical imperfect structure with amplitude ζ r *, at the point M s , where the load parameter M obtains its maximum magnitude (the so-called theoretical load-carrying capacity), the Jacobian of the non-linear system of Equations (5) is equal to zero. The solution is based on controlling the angle of rotation increase (7) for the structures subjected to bending. In this work, a three-mode approach (i.e., J = 3 in Equation (5)) was applied. This means that a model with three degrees of freedom is assumed. Instead of the finite strip method FSM, the exact transition matrix method is used in SAM.
It was decided to continue the analyses using FEM to validate the proposed SAM model and verify the obtained post-buckling equilibrium paths and load-carrying capacities. This comparison allows the determination of the cases for which global secondary buckling should be considered. It is worth noticing that the cost (i.e., time) of the SAM calculation is more than 30 times lower than that of FEM. In addition, the semi-analytical method SAM allows for a much simpler analysis of the observed phenomena and their interpretation compared to FEM.
The calculations were also performed using the finite element method (FEM) with the 18.2 ANSYS ® software version [29]. The composite beam was modelled using 10 754 SHELL181 elements and 65 434 degrees of freedom. SHELL181 is a 4-node element with six degrees of freedom at each node [29]. The size of the element was determined based on the convergence analyses and set to 3 mm. The boundary conditions were applied to ensure pure bending and simple support. To fulfil the condition of simple support, the displacements in the transverse directions at the ends of the beam (u x = 0, u y = 0 in Figures 1 and 2) were removed. Moreover, the displacement in the longitudinal direction was set to zero at the mid-length point of the beam and the mid-width point of the web. The beam loading providing the pure bending was applied in two ways, as follows [5]: • BC I-the beam was loaded by the bending moment generated from normal forces located at the nodes in both beam ends with different force magnitude ( Figure 1). The force distribution corresponds to the stress distribution in the case of pure bending. • BC II-the bending moment was applied by the displacement of the beam ends. The angle of rotation was applied in the "maternode" located at the centre of gravity of the cross-section and transferred to all nodes lying at both ends of the beam cross-sections ( Figure 2). This method of load application corresponds to that used in SAM.
Materials 2020, 13, x FOR PEER REVIEW 5 of 20 direction was set to zero at the mid-length point of the beam and the mid-width point of the web. The beam loading providing the pure bending was applied in two ways, as follows [5]: • BC I-the beam was loaded by the bending moment generated from normal forces located at the nodes in both beam ends with different force magnitude ( Figure 1). The force distribution corresponds to the stress distribution in the case of pure bending. • BC II-the bending moment was applied by the displacement of the beam ends. The angle of rotation was applied in the "maternode" located at the centre of gravity of the cross-section and transferred to all nodes lying at both ends of the beam cross-sections ( Figure 2). This method of load application corresponds to that used in SAM.  The bending moment and angle of rotation were determined separately for each type of load application (BC I and BC II). For BC I, the bending moment M was determined as the sum of forces acting on the beam and the distance from the neutral axis, while the angle of rotation α was determined based on the displacement of the flange in the beam's support, according to the following equation: where: Uzn-displacement in the z-direction of the point located on the flange, on the y-axis; direction was set to zero at the mid-length point of the beam and the mid-width point of the web. The beam loading providing the pure bending was applied in two ways, as follows [5]: • BC I-the beam was loaded by the bending moment generated from normal forces located at the nodes in both beam ends with different force magnitude ( Figure 1). The force distribution corresponds to the stress distribution in the case of pure bending. • BC II-the bending moment was applied by the displacement of the beam ends. The angle of rotation was applied in the "maternode" located at the centre of gravity of the cross-section and transferred to all nodes lying at both ends of the beam cross-sections ( Figure 2). This method of load application corresponds to that used in SAM.  The bending moment and angle of rotation were determined separately for each type of load application (BC I and BC II). For BC I, the bending moment M was determined as the sum of forces acting on the beam and the distance from the neutral axis, while the angle of rotation α was determined based on the displacement of the flange in the beam's support, according to the following equation: where: Uzn-displacement in the z-direction of the point located on the flange, on the y-axis; The bending moment and angle of rotation were determined separately for each type of load application (BC I and BC II). For BC I, the bending moment M was determined as the sum of forces acting on the beam and the distance from the neutral axis, while the angle of rotation α was determined based on the displacement of the flange in the beam's support, according to the following equation: where: U zn -displacement in the z-direction of the point located on the flange, on the y-axis; y max -maximum distance from neutral axis to the outer layer; For BC II, the bending moment M was determined as the reaction (moment) of the applied load in the "masternode", while the angle of rotation α was determined directly from the applied load. The solution differs for the two FEM models: for FEM I BC it is based on controlling the bending moment increase, while for FEM II BC on controlling the angle of rotation increase.
The studies were carried out for thin-walled beams of medium length equal to L = 500 mm with 1 mm wall thickness. Eight-layer symmetrical GFRP laminates with a thickness of individual ply equal to t 1 = 0.125 mm (i.e., t = 8t 1 ) were considered.
The cross-section of the lip channel beam is presented in Figure 3 with dimensions listed in Table 1. Mechanical properties, such as elastic properties and strength limits, were determined in the main orthotropic directions, as listed in Table 2. Depending on fibre orientation, three different stacking sequences (i.e., three instances) were analysed (Table 3).
Materials 2020, 13, x FOR PEER REVIEW 6 of 20 ymax-maximum distance from neutral axis to the outer layer; For BC II, the bending moment M was determined as the reaction (moment) of the applied load in the "masternode", while the angle of rotation α was determined directly from the applied load. The solution differs for the two FEM models: for FEM I BC it is based on controlling the bending moment increase, while for FEM II BC on controlling the angle of rotation increase.
The studies were carried out for thin-walled beams of medium length equal to L = 500 mm with 1 mm wall thickness. Eight-layer symmetrical GFRP laminates with a thickness of individual ply equal to t1 = 0.125 mm (i.e., t = 8t1) were considered.
The cross-section of the lip channel beam is presented in Figure 3 with dimensions listed in Table  1. Mechanical properties, such as elastic properties and strength limits, were determined in the main orthotropic directions, as listed in Table 2. Depending on fibre orientation, three different stacking sequences (i.e., three instances) were analysed (Table 3).

Linear Buckling Analysis
In the first stage, eigenbuckling analysis was performed for LC beams to determine the buckling stresses of the considered buckling modes.
The following indexes were introduced: 1-the primary global distortional-lateral buckling mode for m = 1 (m corresponds to the number of half-waves in the longitudinal direction);

Linear Buckling Analysis
In the first stage, eigenbuckling analysis was performed for LC beams to determine the buckling stresses of the considered buckling modes. The following indexes were introduced: 1-the primary global distortional-lateral buckling mode for m = 1 (m corresponds to the number of half-waves in the longitudinal direction); 2-the secondary global distortional-lateral buckling mode for m = 1; 3-local distortional buckling mode for m = 2; 4-local buckling mode for m > 2; The interactive buckling of three modes was examined in this work. Table 4 presents the critical stress values σ i (i.e., eigenvalues) for the three layer arrangements considered, obtained from FEM and SAM. The lowest values of the bifurcation loads σ i were obtained for: LC-1 for m = 2, LC-2 for m = 2, and LC-3 for m > 2 (for SAM) and for m = 1 (for FEM). Comparing the three analysed instances (LC-1, LC-2 and LC-3), the lowest values of the bifurcation loads for local and global modes were obtained for LC-1, while for LC-2 and LC-3 they are similar and larger by about 30% than for LC-1. For LC-1, the lowest values of the critical load σ i were obtained for local and global modes, with σ 1 /σ 3 = 1.1. Respectively, for LC-2 σ 3 /σ 1 and LC-3 σ 4 /σ 1 (SAM), σ 1 /σ 4 (FEM) were below 1.05. Therefore, the magnitudes of critical load for each instance are close to one another. Based on these results, it can be concluded that the layer arrangement of the LC-1 beam is the least favourable for the considered instances due to the value of the critical loads. Table 4. Comparison of buckling stresses for considered buckling modes.

LC-1
LC-2 LC-3 In Table 5, the values of the critical moment M min corresponding to the lowest critical load for each of the three instances, and also those corresponding to the M min angle of rotation at the support α min are listed. Comparing the three applied methods (FEM BC I, FEM BC II and SAM), a very high result accuracy was obtained for the two FE models. Some deviations can be noticed when SAM outcomes are taken into account. The magnitude of buckling stresses, bending moments and angle of rotation determined from SAM are up to 10% greater than from FEM. This difference between FEM and SAM results is due to the differences in the models in both methods. The number of degrees of freedom, three for SAM and over 65,000 for FEM, reveals that the results obtained from SAM are higher than FEM. The buckling modes (the so-called eigenvectors) corresponding to the analysed instances obtained from FEM are presented in Table 6. For each instance and buckling mode, the cross-sections of the buckling modes for the maximum deflection amplitude are also listed in Table 6. It should be mentioned that eigenvectors are determined in the increments of a unit. Comparing the individual buckling modes for each of the three considered instances LC-1, LC-2 and LC-3, it can be said that the buckling modes are practically identical for the given i = 1, 2, 3 in increments to a unit (i.e., deflections outside or inside for i = 3). For i = 4 (i.e., for the local mode for m > 2), the obtained buckling modes differ in the number of half-waves in the longitudinal direction. Therefore, for LC-1, there are ten half-waves, and 12 and 11 for LC-2 and LC-3, respectively.   As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This  As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This  As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This  As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This  As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This  As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This  As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This  As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This

LC-2
Materials 2020, 13, x FOR PEER REVIEW 8 of 20 As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This  As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This  As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This  As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This  As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This  As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This  As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This  As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This

LC-3
Materials 2020, 13, x FOR PEER REVIEW 8 of 20 As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This  As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This  As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This  As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This  As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This  As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This  As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This  As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the Nx, Ny and Nxy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces Nx and Ny were determined in half of the half-wave length, while Nxy at the end of the beam. The distribution of internal forces Nx, Ny and Nxy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (Nx, Ny or Nxy) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force Nx, the "Relative internal force" is defined as Nx/Nxmax, Ny/Nxmax or Nxy/Nxmax. This As part of the analysis of the eigenvalue problems corresponding to the bifurcation loads, the N x , N y and N xy components of internal forces were determined. For a given critical load, the force components were determined in increments of a unit. Therefore, magnitudes of the internal forces for different bifurcation loads cannot be compared. Only force components for a given eigenvalue can be compared.
Internal forces were determined numerically for the two analysed boundary conditions (FEM BC I, FEM BC II) and compared with the SAM results. The internal forces N x and N y were determined in half of the half-wave length, while N xy at the end of the beam. The distribution of internal forces N x , N y and N xy in the beam cross-section is presented in the charts dimensionless, by introducing the parameter "Relative internal force". This parameter is defined as the ratio of a given internal force (N x , N y or N xy ) to the maximum value of the internal force in the given instance and buckling mode. Taking into account the fact that in each analysed instance, the highest value is always reached by an internal force N x , the "Relative internal force" is defined as N x /N xmax , N y /N xmax or N xy /N xmax . This procedure allows the observation of the share of each internal force. Figure 4 presents indicative results with the distribution of the N x internal force for the primary global distortional-lateral buckling mode in the LC-1 beam cross-section. The force components are presented in the graphs as a force in the function of the cross-sectional length, as shown in Figure 4. The parameter is equal to zero at the beginning of the tensioned edge of the stiffener (bottom stiffener in Figure 1 for example) and gains the maximum value at the end of the compressed edge of the stiffener (top stiffener). For all force components and all modes for the three considered instances, good compatibility between the methods (i.e., FEM and SAM) can be observed.
buckling mode in the LC-1 beam cross-section. The force components are presented in the graphs as a force in the function of the cross-sectional length, as shown in Figure 4. The parameter is equal to zero at the beginning of the tensioned edge of the stiffener (bottom stiffener in Figure 1 for example) and gains the maximum value at the end of the compressed edge of the stiffener (top stiffener). For all force components and all modes for the three considered instances, good compatibility between the methods (i.e., FEM and SAM) can be observed.  Tables 7-9. For each considered buckling mode and stacking sequence, the highest share of Nx internal force is observed. Internal forces Nxy constitute a maximum of 40% of the Nxmax force, while the contribution of the Ny internal force reaches a maximum value of 25% of the Nxmax internal force. It should be emphasised that the highest shares of Ny internal forces occur in the case of the lowest local buckling mode (i = 4). For other considered buckling modes, their contribution is negligible (up to 1.5% of internal force for the secondary global distortional-lateral buckling mode in the LC-2 beam). Similarities can be observed in the distribution of internal forces among the considered samples.  Tables 7-9. For each considered buckling mode and stacking sequence, the highest share of N x internal force is observed. Internal forces N xy constitute a maximum of 40% of the N xmax force, while the contribution of the N y internal force reaches a maximum value of 25% of the N xmax internal force. It should be emphasised that the highest shares of N y internal forces occur in the case of the lowest local buckling mode (i = 4). For other considered buckling modes, their contribution is negligible (up to 1.5% of internal force for the secondary global distortional-lateral buckling mode in the LC-2 beam). Similarities can be observed in the distribution of internal forces among the considered samples.

Non-Linear Analysis
Next, a non-linear stability analysis was performed to determine the influence of the N x , N y and N xy internal force components on the critical equilibrium paths and on load-carrying capacity. In SAM, the buckling interaction occurs only within the first non-linear order of approximation. First-order non-linear coefficients (4) and the corresponding (through dependence (6a)) coefficients (5) are determined based on their eigenvalue modes (i.e., eigenvectors). For a more detailed analysis, see Appendix A. Determining the three index coefficients makes it possible to analyse the influence of each of these modes on the values of the coefficients, which is not possible as part of the eigenvalue problem.
To determine the effect of the most significant buckling modes, the values of a pqr coefficients were determined for each of the equations of equilibrium (5). In this analysis, only the three-mode approach was considered taking into account the interaction of the following buckling modes: • Case 1-the interaction of buckling modes 1, 2 and 3; • Case 2-the interaction of buckling modes 1, 2 and 4;

Nx Nxy Ny
For three instances LC-1, LC-2 and LC-3. For each instance, the largest three index coefficient max(a pqr ) was selected from the three equations of equilibrium (5). Table 10 lists the index designations of these a pqr coefficients and (in parentheses) the elements dependent on dimensionless deflection parameters. Table 10 presents the coefficients which play a major role in the interaction, i.e., coefficients meeting the relationships 0.8max(a pqr ) < a pqr < max(a pqr ) and second-order coefficients satisfying the relationship 0.2max(a pqr ) < a pqr < 0.8max(a pqr ). For a pqr coefficients, the following rule stands: e.g., coefficient a 211 , according to Equation (4), is equal to a 211 = a 112 + a 121 + a 211 . Comparing the results presented in Table 10, it can be observed that for LC-1, the coefficients at ζ 1 ζ 1 2 and ζ 2 ζ 3 2 in energy always play the primary role (7), while coefficients ζ 1 ζ 3 2 are secondary. In the non-linear analysis of a pqr coefficients, it was observed that they depend on the sum of three elements of integration, which rely among others on the internal force components N x , N y , N xy (for a more detailed analysis see Appendix A). Thus, the influence of these internal force components N x , N y , N xy on a pqr coefficients was analysed. For this purpose, the numerical reset for the individual forces' component was calculated in a computer program. The following dimensionless ratios were introduced: I1=a pqr (IX = 0)/a pqr , I2 = a pqr (IY = 0)/a pqr , I3 = a pqr (IXY = 0)/a pqr , where the coefficients occurring in the counter are determined by adopting IX = 0, IY = 0 and IXY = 0 expressed by Equations (A4)-(A7), respectively.
The maximum dimensionless ratios I1, I2, I3 for three instances and two cases are shown in Table 11. It is evident that the N x component has the most significant influence on the values of coefficients a pqr , because its omission reduces the coefficients by at least 94%. If the N y component is omitted, this influence is at most 4%, and for N xy it is at most 7% of the value. This demonstrates the dominant influence of the N x component on the magnitude of a pqr . Table 11. Dimensionless coefficients I1, I2 and I3 for LC-1, LC-2 and LC-3.
In Table 12, the dimensionless load-carrying capacity related to the lowest critical load M s /M min for two cases (i.e., case 1 and case 2) with the three-mode approach and three calculation methods is presented. In SAM, the interaction of the buckling modes for case 1 always provides lower load-carrying capacity M s /M min than for case 2. This indicates a stronger distortional-local influence of mode i = 3 (for m = 2) on interactive buckling than the "pure" local mode i = 4 (for m > 2). In SAM, it is possible to separate all buckling modes as opposed to FEM, where it is virtually impossible. For three instances, lower magnitudes of the M s /M min were obtained for FEM BC I and FEM BC II. The value of M s /M min load-carrying capacity is always lower than the minimum critical moment. For SAM, the load-carrying capacity drops to about 20% for LC-1, LC-2 and LC-3, while for FEM and LC-1 it drops to 2%, for LC-2 to 9%, and for LC-3 to 11%.  Table 13 presents the post-buckling equilibrium paths for three instances where the lowest load carrying-capacity was obtained by the SAM method. According to Equation (7) on post-critical equilibrium paths, the lowest value of the bifurcation moment M min and the corresponding angle of rotation at the support α min should be taken as a reference value. According to Table 4, for LC-1 the lowest value of bifurcation moment corresponds to i = 3, for LC-2 it is also i = 3, and for LC-3, and it is i = 4 (for SAM). Table 13 presents a comparison of the M/M min post-critical equilibrium paths as a function of the rotation angle on the support α/α min for the three used methods. Moreover, in the equilibrium paths obtained from FEM, the load-carrying capacity was marked with dots, while the load leading to composite failure determined with the Hashin failure criteria was marked with "x".
The post-critical equilibrium paths for FEM BC I and FEM BC II are close to each other in the considered ranges of α/α min variability. The curves for SAM are well matched to the curves for FEM. Nevertheless, they reach their maximum value, i.e., the load-carrying capacity M s /M min , below the load-carrying capacity obtained from FEM, and then fall much more rapidly than in the case of FEM. It should be mentioned, however, that SAM is a lower-bound estimation of the structure's load-carrying capacity. The average difference reaches a value close to 10% for LC-1 and LC-2, and the value of 5% for the LC-1 beam.
Using the Hashin failure criteria for the material properties given in Table 2, the maximum bending moments corresponding to the destruction of the first composite layer M H were determined using FEM for the considered instances. In Table 14, as in Table 12 for LC-1, LC-2 and LC-3, the dimensionless maximum loads referred to the lowest critical load value M H /M min for two cases (i.e., case 1 and case 2) for three modal approaches. The post-critical equilibrium paths are very flat near the load-carrying capacity. This results in the M H /M min values for the corresponding cases being similar to M s /M min (Tables 12 and 13). For LC-2, layer failure occurs immediately after reaching load-carrying capacity, while for the other two instances composite failure occurs before reaching the ultimate load. For each of the considered beams, failure occurred as a result of the matrix cracking under compression. Tables 14-17 show the distribution of internal forces N x , N y and N xy corresponding to the maximum moments M H /M min (Table 14) and the load-carrying capacity M s /M min (Table 12) obtained for FEM BC I and FEM BC II. Solid lines denote the distributions for the maximum moment M H /M min , and dashed lines those for M s /M min . For each considered instance, the largest share of N x internal force is observed. Internal force N xy constitutes a maximum of 20% of the N xmax force, while the contribution of the N y internal force reaches a maximum value of 3% of the N xmax internal force. In this case, the N x , N y and N xy internal force distributions are, as expected, more non-linear than in the case of distributions corresponding to a linear eigenvalue problem. The N x and N xy internal force distributions corresponding to the first failure (detected with the Hashin failure criteria) and the load-carrying capacity are similar. Larger differences can be observed for the N y internal force distributions. The graphs of the N x internal force distribution, considering the three-mode approach for modes i = 1, 2, 3 and i = 1, 2, 4 are very similar, as expected for FEM. LC-2 may raise some doubts. It should be noted; however, that when the absolute values of forces are taken into consideration, these graphs are coincident, which corresponds to the critical symmetrical equilibrium path.
SAM, as used in this study, explained the phenomenon of rapid loss of load-carrying capacity for medium-length LC-composite beams. FEM enables the verification of the observed phenomena. This phenomenon is always associated with the influence of the secondary global distortional-lateral buckling mode with the primary global distortional-lateral and local modes. This interaction was analysed taking into account the components of internal forces N x , N y and N xy . FEM was used to verify the observed phenomenon and analyse it more precisely. In addition, FEM made it possible to compare the determined load-carrying capacity with the maximum moment determined based on the Hashin failure criteria and to illustrate the distribution of the internal forces for these load values. The linear and non-linear analysis shows a very large influence of the internal force component N x on the bifurcation loads and load-carrying capacity of the LC-beams.      (Table 14) and the load-carrying capacity Ms/Mmin (Table 12) obtained for FEM BC I and FEM BC II. Solid lines denote the distributions for the maximum moment MH/Mmin, and dashed lines those for Ms/Mmin. For each considered instance, the largest share of Nx internal force is observed. Internal force Nxy constitutes a maximum of 20% of the Nxmax force, while the contribution of the Ny internal force reaches a maximum value of 3% of the Nxmax internal force. In this case, the Nx, Ny and Nxy internal force distributions are, as expected, more non-linear than in the case of distributions corresponding to a linear eigenvalue problem. The Nx and Nxy internal force distributions corresponding to the first failure (detected with the Hashin failure criteria) and the loadcarrying capacity are similar. Larger differences can be observed for the Ny internal force distributions. The graphs of the Nx internal force distribution, considering the three-mode approach for modes i = 1, 2, 3 and i = 1, 2, 4 are very similar, as expected for FEM. LC-2 may raise some doubts. It should be noted; however, that when the absolute values of forces are taken into consideration, these graphs are coincident, which corresponds to the critical symmetrical equilibrium path.  (Table 14) and the load-carrying capacity Ms/Mmin (Table 12) obtained for FEM BC I and FEM BC II. Solid lines denote the distributions for the maximum moment MH/Mmin, and dashed lines those for Ms/Mmin. For each considered instance, the largest share of Nx internal force is observed. Internal force Nxy constitutes a maximum of 20% of the Nxmax force, while the contribution of the Ny internal force reaches a maximum value of 3% of the Nxmax internal force. In this case, the Nx, Ny and Nxy internal force distributions are, as expected, more non-linear than in the case of distributions corresponding to a linear eigenvalue problem. The Nx and Nxy internal force distributions corresponding to the first failure (detected with the Hashin failure criteria) and the loadcarrying capacity are similar. Larger differences can be observed for the Ny internal force distributions. The graphs of the Nx internal force distribution, considering the three-mode approach for modes i = 1, 2, 3 and i = 1, 2, 4 are very similar, as expected for FEM. LC-2 may raise some doubts. It should be noted; however, that when the absolute values of forces are taken into consideration, these graphs are coincident, which corresponds to the critical symmetrical equilibrium path.  (Table 14) and the load-carrying capacity Ms/Mmin (Table 12) obtained for FEM BC I and FEM BC II. Solid lines denote the distributions for the maximum moment MH/Mmin, and dashed lines those for Ms/Mmin. For each considered instance, the largest share of Nx internal force is observed. Internal force Nxy constitutes a maximum of 20% of the Nxmax force, while the contribution of the Ny internal force reaches a maximum value of 3% of the Nxmax internal force. In this case, the Nx, Ny and Nxy internal force distributions are, as expected, more non-linear than in the case of distributions corresponding to a linear eigenvalue problem. The Nx and Nxy internal force distributions corresponding to the first failure (detected with the Hashin failure criteria) and the loadcarrying capacity are similar. Larger differences can be observed for the Ny internal force distributions. The graphs of the Nx internal force distribution, considering the three-mode approach for modes i = 1, 2, 3 and i = 1, 2, 4 are very similar, as expected for FEM. LC-2 may raise some doubts. It should be noted; however, that when the absolute values of forces are taken into consideration, these graphs are coincident, which corresponds to the critical symmetrical equilibrium path.  (Table 14) and the load-carrying capacity Ms/Mmin (Table 12) obtained for FEM BC I and FEM BC II. Solid lines denote the distributions for the maximum moment MH/Mmin, and dashed lines those for Ms/Mmin. For each considered instance, the largest share of Nx internal force is observed. Internal force Nxy constitutes a maximum of 20% of the Nxmax force, while the contribution of the Ny internal force reaches a maximum value of 3% of the Nxmax internal force. In this case, the Nx, Ny and Nxy internal force distributions are, as expected, more non-linear than in the case of distributions corresponding to a linear eigenvalue problem. The Nx and Nxy internal force distributions corresponding to the first failure (detected with the Hashin failure criteria) and the loadcarrying capacity are similar. Larger differences can be observed for the Ny internal force distributions. The graphs of the Nx internal force distribution, considering the three-mode approach for modes i = 1, 2, 3 and i = 1, 2, 4 are very similar, as expected for FEM. LC-2 may raise some doubts. It should be noted; however, that when the absolute values of forces are taken into consideration, these graphs are coincident, which corresponds to the critical symmetrical equilibrium path.  (Table 14) and the load-carrying capacity Ms/Mmin (Table 12) obtained for FEM BC I and FEM BC II. Solid lines denote the distributions for the maximum moment MH/Mmin, and dashed lines those for Ms/Mmin. For each considered instance, the largest share of Nx internal force is observed. Internal force Nxy constitutes a maximum of 20% of the Nxmax force, while the contribution of the Ny internal force reaches a maximum value of 3% of the Nxmax internal force. In this case, the Nx, Ny and Nxy internal force distributions are, as expected, more non-linear than in the case of distributions corresponding to a linear eigenvalue problem. The Nx and Nxy internal force distributions corresponding to the first failure (detected with the Hashin failure criteria) and the loadcarrying capacity are similar. Larger differences can be observed for the Ny internal force distributions. The graphs of the Nx internal force distribution, considering the three-mode approach for modes i = 1, 2, 3 and i = 1, 2, 4 are very similar, as expected for FEM. LC-2 may raise some doubts. It should be noted; however, that when the absolute values of forces are taken into consideration, these graphs are coincident, which corresponds to the critical symmetrical equilibrium path.  (Table 14) and the load-carrying capacity Ms/Mmin (Table 12) obtained for FEM BC I and FEM BC II. Solid lines denote the distributions for the maximum moment MH/Mmin, and dashed lines those for Ms/Mmin. For each considered instance, the largest share of Nx internal force is observed. Internal force Nxy constitutes a maximum of 20% of the Nxmax force, while the contribution of the Ny internal force reaches a maximum value of 3% of the Nxmax internal force. In this case, the Nx, Ny and Nxy internal force distributions are, as expected, more non-linear than in the case of distributions corresponding to a linear eigenvalue problem. The Nx and Nxy internal force distributions corresponding to the first failure (detected with the Hashin failure criteria) and the loadcarrying capacity are similar. Larger differences can be observed for the Ny internal force distributions. The graphs of the Nx internal force distribution, considering the three-mode approach for modes i = 1, 2, 3 and i = 1, 2, 4 are very similar, as expected for FEM. LC-2 may raise some doubts. It should be noted; however, that when the absolute values of forces are taken into consideration, these graphs are coincident, which corresponds to the critical symmetrical equilibrium path.

LC-3 Nx Nxy Ny
Interactio n i = 1,2,3 Interactio n i = 1,2,4 SAM, as used in this study, explained the phenomenon of rapid loss of load-carrying capacity for medium-length LC-composite beams. FEM enables the verification of the observed phenomena. This phenomenon is always associated with the influence of the secondary global distortional-lateral buckling mode with the primary global distortional-lateral and local modes. This interaction was analysed taking into account the components of internal forces Nx, Ny and Nxy. FEM was used to verify the observed phenomenon and analyse it more precisely. In addition, FEM made it possible to compare the determined load-carrying capacity with the maximum moment determined based on

LC-3 Nx Nxy Ny
Interactio n i = 1,2,3 Interactio n i = 1,2,4 SAM, as used in this study, explained the phenomenon of rapid loss of load-carrying capacity for medium-length LC-composite beams. FEM enables the verification of the observed phenomena. This phenomenon is always associated with the influence of the secondary global distortional-lateral buckling mode with the primary global distortional-lateral and local modes. This interaction was analysed taking into account the components of internal forces Nx, Ny and Nxy. FEM was used to verify the observed phenomenon and analyse it more precisely. In addition, FEM made it possible to compare the determined load-carrying capacity with the maximum moment determined based on

LC-3 Nx Nxy Ny
Interactio n i = 1,2,3 Interactio n i = 1,2,4 SAM, as used in this study, explained the phenomenon of rapid loss of load-carrying capacity for medium-length LC-composite beams. FEM enables the verification of the observed phenomena. This phenomenon is always associated with the influence of the secondary global distortional-lateral buckling mode with the primary global distortional-lateral and local modes. This interaction was analysed taking into account the components of internal forces Nx, Ny and Nxy. FEM was used to verify the observed phenomenon and analyse it more precisely. In addition, FEM made it possible to compare the determined load-carrying capacity with the maximum moment determined based on

LC-3 Nx Nxy Ny
Interactio n i = 1,2,3 Interactio n i = 1,2,4 SAM, as used in this study, explained the phenomenon of rapid loss of load-carrying capacity for medium-length LC-composite beams. FEM enables the verification of the observed phenomena. This phenomenon is always associated with the influence of the secondary global distortional-lateral buckling mode with the primary global distortional-lateral and local modes. This interaction was analysed taking into account the components of internal forces Nx, Ny and Nxy. FEM was used to verify the observed phenomenon and analyse it more precisely. In addition, FEM made it possible to compare the determined load-carrying capacity with the maximum moment determined based on

LC-3 Nx Nxy Ny
Interactio n i = 1,2,3 Interactio n i = 1,2,4 SAM, as used in this study, explained the phenomenon of rapid loss of load-carrying capacity for medium-length LC-composite beams. FEM enables the verification of the observed phenomena. This phenomenon is always associated with the influence of the secondary global distortional-lateral buckling mode with the primary global distortional-lateral and local modes. This interaction was analysed taking into account the components of internal forces Nx, Ny and Nxy. FEM was used to verify the observed phenomenon and analyse it more precisely. In addition, FEM made it possible to compare the determined load-carrying capacity with the maximum moment determined based on

LC-3 Nx Nxy Ny
Interactio n i = 1,2,3 Interactio n i = 1,2,4 SAM, as used in this study, explained the phenomenon of rapid loss of load-carrying capacity for medium-length LC-composite beams. FEM enables the verification of the observed phenomena. This phenomenon is always associated with the influence of the secondary global distortional-lateral buckling mode with the primary global distortional-lateral and local modes. This interaction was analysed taking into account the components of internal forces Nx, Ny and Nxy. FEM was used to verify the observed phenomenon and analyse it more precisely. In addition, FEM made it possible to compare the determined load-carrying capacity with the maximum moment determined based on

LC-3 Nx Nxy Ny
Interactio n i = 1,2,3 Interactio n i = 1,2,4 SAM, as used in this study, explained the phenomenon of rapid loss of load-carrying capacity for medium-length LC-composite beams. FEM enables the verification of the observed phenomena. This phenomenon is always associated with the influence of the secondary global distortional-lateral buckling mode with the primary global distortional-lateral and local modes. This interaction was analysed taking into account the components of internal forces Nx, Ny and Nxy. FEM was used to verify the observed phenomenon and analyse it more precisely. In addition, FEM made it possible to Materials 2020, 13, x FOR PEER REVIEW 17 of 20 Table 16. Internal forces in LC-2 beam under the load corresponding to the first failure and the loadcarrying capacity.

LC-3 Nx Nxy Ny
Interactio n i = 1,2,3 Interactio n i = 1,2,4 SAM, as used in this study, explained the phenomenon of rapid loss of load-carrying capacity for medium-length LC-composite beams. FEM enables the verification of the observed phenomena. This phenomenon is always associated with the influence of the secondary global distortional-lateral buckling mode with the primary global distortional-lateral and local modes. This interaction was analysed taking into account the components of internal forces Nx, Ny and Nxy. FEM was used to verify the observed phenomenon and analyse it more precisely. In addition, FEM made it possible to Materials 2020, 13, x FOR PEER REVIEW 17 of 20 Table 16. Internal forces in LC-2 beam under the load corresponding to the first failure and the loadcarrying capacity.

LC-3 Nx Nxy Ny
Interactio n i = 1,2,3 Interactio n i = 1,2,4 SAM, as used in this study, explained the phenomenon of rapid loss of load-carrying capacity for medium-length LC-composite beams. FEM enables the verification of the observed phenomena. This phenomenon is always associated with the influence of the secondary global distortional-lateral buckling mode with the primary global distortional-lateral and local modes. This interaction was analysed taking into account the components of internal forces Nx, Ny and Nxy. FEM was used to verify the observed phenomenon and analyse it more precisely. In addition, FEM made it possible to  Table 16. Internal forces in LC-2 beam under the load corresponding to the first failure and the loadcarrying capacity.

LC-3 Nx Nxy Ny
Interactio n i = 1,2,3 Interactio n i = 1,2,4 SAM, as used in this study, explained the phenomenon of rapid loss of load-carrying capacity for medium-length LC-composite beams. FEM enables the verification of the observed phenomena. This phenomenon is always associated with the influence of the secondary global distortional-lateral buckling mode with the primary global distortional-lateral and local modes. This interaction was analysed taking into account the components of internal forces Nx, Ny and Nxy. FEM was used to verify the observed phenomenon and analyse it more precisely. In addition, FEM made it possible to Materials 2020, 13, x FOR PEER REVIEW 17 of 20 Table 16. Internal forces in LC-2 beam under the load corresponding to the first failure and the loadcarrying capacity.

LC-3 Nx Nxy Ny
Interactio n i = 1,2,3 Interactio n i = 1,2,4 SAM, as used in this study, explained the phenomenon of rapid loss of load-carrying capacity for medium-length LC-composite beams. FEM enables the verification of the observed phenomena. This phenomenon is always associated with the influence of the secondary global distortional-lateral buckling mode with the primary global distortional-lateral and local modes. This interaction was analysed taking into account the components of internal forces Nx, Ny and Nxy. FEM was used to verify the observed phenomenon and analyse it more precisely. In addition, FEM made it possible to Materials 2020, 13, x FOR PEER REVIEW 17 of 20 Table 16. Internal forces in LC-2 beam under the load corresponding to the first failure and the loadcarrying capacity.

LC-3 Nx Nxy Ny
Interactio n i = 1,2,3 Interactio n i = 1,2,4 SAM, as used in this study, explained the phenomenon of rapid loss of load-carrying capacity for medium-length LC-composite beams. FEM enables the verification of the observed phenomena. This phenomenon is always associated with the influence of the secondary global distortional-lateral buckling mode with the primary global distortional-lateral and local modes. This interaction was analysed taking into account the components of internal forces Nx, Ny and Nxy. FEM was used to verify the observed phenomenon and analyse it more precisely. In addition, FEM made it possible to These coefficients are determined based only on linear eigenvalues. This is an important observation simplifying solutions to problems limited to the first approximation. Coefficientsā pqr (A3) are determined from the following integral relationships: σ p · l 11 (U q , U r ) = n i = 1 [N px (u q,x u r,x + v q,x v r,x + w q,x w r,x ) + N py (u q,y u r,y + v q,y v r,y + w q,y w r,y ) +N pxy (u q,x u r,y + u q,y u r,x + v q,x v r,y + v q,y v r,x + w q,x w r,y + w q,y w r,x )]dxdy = IX + IY + IXY (A4) where the abbreviated designation of integer expressions corresponding to the internal forces is introduced: N px (u q,x u r,x + v q,x v r,x + w q,x w r,x ) dxdy (A5) N py (u q,y u r,y + v q,y v r,y + w q,y w r,y ) dxdy (A6) N pxy (u q,x u r,y + u q,y u r,x + v q,x v r,y + v q,y v r,x + w q,x w r,y + w q,y w r,x ) dxdy (A7) and n-the number of component plates. As is evident, the internal forces correspond only to the buckling modes marked with the index p, while the displacement derivatives correspond to the buckling modes marked with indices q, r.