Lateral–Torsional Buckling of Cantilever Steel Beams under 2 Types of Complex Loads

: Cantilever steel beams are an essential structural element in civil engineering ﬁelds such as bridges and buildings. However, there is very little research on the critical moment ( M cr ) of cantilever beams subjected to a concentrated load (CL) or a combination of concentrated load and uniformly distributed load (CUDL) when the concentrated load is not limited to the free end. Therefore, the focus of the current paper is the calculation of M cr for cantilever steel beams under CL and CUDL. This paper proposes a program and a simple closed-form solution for M cr that are applicable to the elastic buckling analysis of cantilever I-beams under CL and CUDL. Based on the Rayleigh–Ritz method, a matrix equation and the corresponding procedure about M cr under CL and CUDL are derived by using inﬁnite trigonometric series for the buckling deformation functions. The value of M cr and the corresponding mode of buckling can be obtained efﬁciently by considering the symmetry of the section, the ratio of two load values and the load action position. Experimental results and ﬁnite element calculations validate the numerical solutions of the procedure. A closed-form solution for M cr is derived according to the assumption of a small torsion angle and the speciﬁc values of each coefﬁcient in the closed-form solution of M cr are calculated by the proposed procedure. The results show that the procedure and closed-form solution for M cr presented in this paper have a high degree of accuracy in calculating the M cr of the cantilever beam under CL and CUDL. The deviations between the results calculated by the proposed procedure and data from existing literature are less than 8%. These conclusions are capable of solving the calculation problem of M cr for cantilever beams under CL or CUDL, which are both signiﬁcant load cases in engineering. The study provides a reference for the design of cantilever steel beams.


Introduction
Lateral-torsional buckling (LTB) is a common global instability phenomenon for slender structures, where structures that are bent in the plane of greatest flexural rigidity may buckle laterally and torsionally as the external load reaches the critical value.Therefore, it is necessary to consider LTB in the design of the beam when the flexural rigidity of the beam in the plane of bending is greater than the lateral flexural rigidity.Cantilevers are a common slender structure in civil engineering.Therefore, it is necessary to consider the limiting load of LTB in addition to deformation and stress analysis during the engineering design process.The American National Standard Institute/American Institute of Steel Construction (ANSI/AISC 360-16) proposed the lateral-torsional buckling modification factor, C b , when calculating the LTB of simply supported beams.C b had an explicit formula but it was only recommended directly taking the value of 1.0 when calculating the C b of cantilever beams [1]; EN 1993-1-1:2005 also had a specific formulas for lateral-torsional buckling of simply supported beams, but it did not mention the instability calculation of cantilever [2].
As early as 1960, Clark proposed a formula for calculating the critical moment (M cr ) using C 1 , C 2 and C 3 factors [3].Since then, the calculation method of C 1 , C 2 and C 3 has been increasingly refined.Many scholars have studied C 1 , C 2 and C 3 factors of lateral-torsional buckling of simply supported beams.Zhang et al. [4] studied the general formula of M cr under linear distributed moment.The study used dimensionless parameters and infinite series to calculate M cr and fitted the specific values of C 1 , C 2 and C 3 .It was revealed that the more terms the buckling deformation function takes, the more accurate the calculation results will be.However, due to the complexity of the external load, the solution to M cr also needs to depend on other methods.Bresser et al. [5] further studied the expression of C b of simply supported beams under support moments, point loads, uniformly distributed loads and a combination thereof.Kucukler et al. [6] presented a theory of stiffness reduction coefficient for the LTB design of welded web-tapered steel beams.Rossi et al. [7] conducted elastic LTB analysis on supported beams using ABAQUS and GBTul.Sahraei et al. [8] assessed the LTB mode of plane frames by finite element theory, which greatly reduced the amount of modeling and calculation.So far, numerous problems about the elastic LTB of simply supported beams are being gradually solved.
For cantilever steel beams, buckling deformation is more complex due to the characteristics of the boundary condition.Unlike simply supported beams, the maximum displacement and torsion angle of cantilever are located at the free end rather than near the mid-span.Therefore, the coefficients of C 1 , C 2 and C 3 are no longer fixed values.Timoshenko used infinite series to study the singly symmetric cross-section cantilever beam and obtained the critical load of the lateral-torsional buckling [9].Anderson and Attard [10,11] carried out LTB experiments on four types of I-shaped cantilever steel beams with different cross-sections.Andrade et al. [12] extended the domain of application of C 1 , C 2 and C 3 factors to cantilevers.In addition, Andrade evaluated the performance of a 1D model of elastic LTB behavior [13].Since then, many scholars have calculated M cr 's formula of cantilever beams with different sections under a large number of different loads.Zhang et al. [14] proposed a set of formulae to calculate M cr under concentrated load at the free end or uniformly distributed load.Ozbasaran et al. [15,16] proposed a closed-form expression for calculating elastic critical LTB loads and a design procedure for cantilevers.The design procedure considered elastic buckling and inelastic buckling and the calculation results of the procedure were highly consistent with Eurocode 3-2005 or AISC360-10.Because the boundary conditions of cantilever beams are more complex than those of simply supported beams, the analytical solutions of M cr of the cantilever are highly complex.Therefore, the solution to M cr requires more approaches, such as the finite difference [17], finite element [18][19][20], energy method [21,22], finite integral [10] or the lateral-torsional buckling modification factor [1].Among them, the energy method is one of the most common and basic methods.Therefore, the correct derivation or selection of total potential energy is the key to obtaining accurate calculation results.Additionally, Trahair et al. [23,24] studied the inelastic LTB of cantilevers and continuous steel beams in order to provide better design methods than current specifications.Similarly, Demirhan et al. [25] carried out inelastic buckling experiments on biaxial symmetric cantilever steel beams and proposed a numerical simulation method to accurately predict the displacement and inelastic critical buckling load of cantilever steel beams.
In addition to the study of ordinary simply supported beams and cantilever beams, many scholars also considered other factors, such as pre-stressed beams, initial imperfections, flange-web interaction and material properties.Lorkowski et al. [26] studied the M cr of LTB of pre-stressed two-chord columns by experiment and numerical analysis.Kim et al. [27,28] presented an LTB theory of pre-stressed steel H-beams.This theory considered the initial rotation angle and it was applicable to simply supported beams and cantilever beams.Zhang used the energy method to obtain the buckling load of pre-stressed steel I-beams [29].Lebastard et al. [30] proposed analytical formulations for supported beams whose support was restrained from warping, and some scholars [31][32][33] have studied the LTB of special sections such as built-up sections.Saoula et al. [34] proposed a new theory to predict the shape and amplitude of unique global and local initial imperfections for the LTB of doubly symmetric I-section cantilevers.Erkmen [35] and Kimura [36] studied the elastic LTB mode considering web-distortion or flange-web interaction.Jáger et al. [37][38][39] studied the elastic LTB mode of corrugated web beams.In addition, an increasing number of scholars have considered the influence of material properties on the LTB performance of beams [40][41][42].For example, Khalaj et al. [43,44] studied the effect of chemical composition of materials on their strength by using neural networks and artificial intelligence technology.Virgin et al. [45] studied the LTB performance of cantilever beams using 3D printing technology.
It could be seen from those studies on the lateral-torsional buckling of cantilever steel beams that most of their transverse loads were mainly a single load.The position of the load acting on the cross-section was restricted to the top flange (TF), bottom flange (BF) and the shear center (SC).The position of the load acting in the beam length direction was restricted to the free end.There was less research on lateral-torsional buckling of singly symmetric I-beams under CL or CUDL.In engineering practice, cantilever steel beams are usually affected by the concentrated load or the combination of concentrated load and uniformly distributed load, and the position of the concentrated load is not limited to the free end.For example, some cantilever beams will be subject to the concentrated load from the secondary beams and uniformly distributed load from the floor slabs.The position of the concentrated load may be in the middle of the span or other positions.In view of this, CL and CUDL are both very important load cases in the engineering field.Therefore, it is necessary to study this situation.
The innovation of this paper is the numerical method and the closed-form solution about M cr , which aims to efficiently obtain the M cr and corresponding mode of buckling of cantilever steel beams under CL and CUDL.A matrix equation with M cr is developed by using infinite trigonometric series for the buckling deformation functions and Rayleigh-Ritz method.Based on the matrix equation, a MATLAB program was developed for LTB analysis of cantilever I-beams.The MATLAB program considered the ratio of two loads, the symmetry of the section, the action position of the concentrated load in the direction of the beam length and the vertical position of two loads.It is worth noting that the matrix equation and MATLAB program can efficiently obtain the M cr and buckling mode of the cantilever beam under any single load except constant moment.Therefore, it has a wide range of applications.Moreover, it was able to solve M cr and the buckling mode of cantilever beams under CL or CUDL, which are both common load cases in engineering applications.Then M cr calculated by the MATLAB program was compared with results reported in existing references [10,11,16,18] to verify the correctness of the theory.Last, a simple closed-form solution for M cr under CL and CUDL is proposed by analyzing the factors affecting the LTB modes of cantilever beams, which enables computing M cr with relative ease.

Introduction and Comparison of Total Potential Energy Equation
Figures 1 and 2 illustrate the LTB of a singly symmetric cantilever I-beam.As shown in Figures 1 and 2, L is cantilever length.P y is concentrated load and q y is uniformly distributed load, their coordinates are respectively (0, y Py ) and (0, y qy ), the shear center's coordinates are (0, y s ), where a Py = y Py − y s , a qy = y qy − y s .λ is the ratio of distance between the action point of P y and the fixed end to the cantilever's length.
Zhang and Tong [14,46] thought that the LTB theory of steel beams should consider linear strain energy, nonlinear longitudinal strain energy, nonlinear shear strain energy and nonlinear transverse strain energy at the same time.Compared with the experimental results in reference [10] and reference [11], the total potential energy proposed by Zhang is accurate in calculating the M cr of the cantilever beam.At the same time, it is also confirmed that the total potential energy proposed by Zhang and Tong is applicable to the boundary conditions of cantilever beam through the discussion in reference [14].
where E is Young's modulus, G is the shear modulus, Iω is the warping constant, It is the torsional constant, Iy is the moment of inertia about the y-axis, βx is Wagner's coefficient, u and φ are the lateral deflection and twist angle of section, respectively, φ(λl) is the twist angle at the point where the transverse concentrated load acts during the lateral-torsional buckling of the beam, and Mx is the moment around the x-axis.The ratio of the uniformly distributed load to concentrated load in Figure 1 is assumed to be as follows: The magnitude of the moment at any position of cantilever (Mx) is expressed by Equation (3).where E is Young's modulus, G is the shear modulus, Iω is the warping constant, It is the torsional constant, Iy is the moment of inertia about the y-axis, βx is Wagner's coefficient, u and φ are the lateral deflection and twist angle of section, respectively, φ(λl) is the twist angle at the point where the transverse concentrated load acts during the lateral-torsional buckling of the beam, and Mx is the moment around the x-axis.The ratio of the uniformly distributed load to concentrated load in Figure 1 is assumed to be as follows: The magnitude of the moment at any position of cantilever (Mx) is expressed by Equation (3).For the load case shown in Figures 1 and 2, the total potential energy of Zhang and Tong can be expressed as Equation (1).
where E is Young's modulus, G is the shear modulus, I ω is the warping constant, I t is the torsional constant, I y is the moment of inertia about the y-axis, β x is Wagner's coefficient, u and ϕ are the lateral deflection and twist angle of section, respectively, ϕ(λl) is the twist angle at the point where the transverse concentrated load acts during the lateral-torsional buckling of the beam, and M x is the moment around the x-axis.
The ratio of the uniformly distributed load to concentrated load in Figure 1 is assumed to be as follows: The magnitude of the moment at any position of cantilever (M x ) is expressed by Equation (3).
Note that when λ is equal to 0, M x takes the second function in Equation (3).When λ is equal to 1, M x takes the first function in Equation (3).
According to Equation (3), boundary conditions and the properties of the integral, the following equations can be obtained. (5) If λ = 1, then M x (l) = P y and ϕ(0) = 0.The following equation can be obtained by Equations ( 5) and (6).
If λ < 1, the following equation can be obtained by Equations ( 5) and ( 6): Equation ( 9) can be obtained by substituting Equations ( 4), (7) and (8) into Equation (1).Equation ( 9) is the expression of the traditional total potential energy in reference [47].Therefore, two total potential energy equations are equivalent under the load case studied in this paper.For convenience of calculation, Equation (9) with the simpler form is selected as the basis for calculation.

The Buckling Deformation Functions
When using the Rayleigh-Ritz method for stability analysis, the assumption of the buckling deformation functions is very critical, which directly affects the calculation accuracy.First, the buckling deformation functions should match the actual deformation form as much as possible and can be expressed by infinite series.Secondly, the buckling deformation functions must meet the displacement boundary conditions.If it can also meet the mechanical boundary conditions, the accuracy will be higher.Based on the above situation, the following trigonometric series are used to describe the lateral deflection u(z) and twist angle ϕ(z) of cantilever steel beams [4,12].
where h is the distance between the centroids of the top and bottom flange, and A n and B n are mutually independent generalized coordinates.From the perspective of structural dynamics, it can be seen that the dimensions of A n and B n are different.Therefore, h is introduced into the expression of u(z) to eliminate the dimension of A n , so that A n and B n are both undetermined dimensionless coefficients.This idea is from reference [4] and it can make the following calculation easier.Equations (10) and (11) satisfy the displacement boundary condition, i.e., u(0) = u (0) = 0 and ϕ(0) = ϕ (0) = 0.For the free end, it satisfies u(l) = 0, u (l) = 0, ϕ(l) = 0 and ϕ (l) = 0.

Calculation of Total Potential Energy Equation
The maximum absolute value of moment (|M 0 |) is expressed as the following equation: On substituting the moment function (3) and the buckling deformation functions (10) and (11) into the total potential energy Equation ( 9) and integrating over the beam length, one can obtain the integration results of the total potential energy.In order to make the results more simple, new dimensionless parameters in reference [4] are introduced below.
where I 1 is the moment of inertia of the top flange around the y-axis, and I 2 is the moment of inertia of the bottom flange around the y-axis.Some equations are derived from Equation ( 13) and defined thus: The expression of the total potential energy equation is divided into four parts.Each part is integrated and multiplied by l 3 /h 2 EI y .Equations ( 13)- (15) are substituted to obtain the specific potential energy results.The specific results are as follows.16) The total potential energy equation can be written as the sum of Equations ( 16)-( 19): According to the principle of stationary total potential energy, the stability conditions of the system can be obtained as follows.
The following expressions are made available by Equation (21).
Appl.Sci.2023, 13, 5830 The condition for the linear equations to hold is that the coefficient determinant of the independent parameters A n and B n is zero.For convenience of writing, the following matrices can be defined.
Equation ( 34) is derived from Equation (1) through integration.In the integration process, the buckling deformation functions Equations ( 10) and ( 11) are substituted into Equation (1).Since Equations ( 1), ( 10) and ( 11) meet the boundary conditions of cantilever beams, Equation ( 34) is applicable to the calculation of the M cr of cantilever beams.
For convenience of writing, the following matrices can be defined.
It can be seen from Equation ( 34) that the solution of the dimensionless critical moment The flowchart given in Figure 3 elaborates on the idea of the procedure.

Verification of the Results and Discussion
The proposed matrix equation is verified according to the experimental data [10,11,1 and finite element analysis data [16] of cantilever beams in existing literature.The dimen sions of specimens are shown in Figure 4 and Table 1.

Verification of the Results and Discussion
The proposed matrix equation is verified according to the experimental data [10,11,18] and finite element analysis data [16] of cantilever beams in existing literature.The dimensions of specimens are shown in Figure 4 and Table 1.Because it is difficult to obtain the LTB load of a cantilever beam under the uniformly distributed load and other complex loads, most lateral-torsional buckling tests applied concentrated load at the free end of the cantilever beam.In order to verify the accuracy of the method in this paper and the applicability of the total potential energy equation based on the experimental data, it is necessary to make α = 0 and λ = 1.At this time, the load case of the cantilever beam became a concentrated load at the free end.The experimental (Mex) [10,11,18] and theoretically predicted (Mth) critical moment and the Mth/Mex ratio are given in Table 2 for all the selected specimens.The specimens in reference [18] are not numbered, so the vertical position of the concentrated load is used instead of the specimen numbers.At the same time, as the experimental data can only verify the load case of CL, it is necessary to extract the data of finite element analysis from existing literature to verify the load case of CUDL.The critical moment calculated by finite element analysis (MFEM) from reference [16] and theoretically predicted (Mth) critical moment and the Mth/MFEM ratio are given in Tables 3-5 for all the selected specimens.
In Table 2, compared with the experimental results, the maximum value of the Mth/Mex ratio is 1.04, and the minimum value is 0.92.In Tables 3-5  Because it is difficult to obtain the LTB load of a cantilever beam under the uniformly distributed load and other complex loads, most lateral-torsional buckling tests applied concentrated load at the free end of the cantilever beam.In order to verify the accuracy of the method in this paper and the applicability of the total potential energy equation based on the experimental data, it is necessary to make α = 0 and λ = 1.At this time, the load case of the cantilever beam became a concentrated load at the free end.The experimental (M ex ) [10,11,18] and theoretically predicted (M th ) critical moment and the M th /M ex ratio are given in Table 2 for all the selected specimens.The specimens in reference [18] are not numbered, so the vertical position of the concentrated load is used instead of the specimen numbers.At the same time, as the experimental data can only verify the load case of CL, it is necessary to extract the data of finite element analysis from existing literature to verify the load case of CUDL.The critical moment calculated by finite element analysis (M FEM ) from reference [16] and theoretically predicted (M th ) critical moment and the M th /M FEM ratio are given in Tables 3-5 for all the selected specimens.In Table 2, compared with the experimental results, the maximum value of the M th /M ex ratio is 1.04, and the minimum value is 0.92.In Tables 3-5, there are four sets of data with significant differences.These data came from SC under CL in Table 4.In reference [16], the distance between the SC and TF of the section II with a reinforced top flange is very small and the differences between their M cr should not be significant.However, the differences in Table 4 are particularly significant and it can be concluded from reference [16] that M th is more accurate than M FEM .Therefore, these four sets of data are not included in the analysis scope.In Tables 3-5, compared with the data of finite element analysis from reference [16], the maximum value of the M th /M FEM ratio is 1.07, and the minimum value is 0.97.The average of all the ratios is 0.98 with a standard deviation of 0.13.Therefore, critical moment can be accurately solved by the numerical method proposed in this paper.
The predictions of the matrix equation and the data from existing literature are compared in Figure 5.The comparison method in Figure 5 uses the modified version of demerit point classification (DPC) [40].DPC is used to evaluate the predictions of the matrix equation in terms of safety, accuracy and economic aspects.A specific explanation of the modified version of DPC can be found in Table 6.According to Table 6 and Figure 5, a penalty is assigned to each range of the M th /M ex or M th /M FEM ratio.The predictions based on Equation ( 34) are all safe in Figure 5.To sum up, the predictions of 100% of the specimens are within the appropriate safety range (total penalty = 0).Therefore, it can be concluded that Equation (34) resulted in a positive prediction in terms of safety, accuracy and economic aspects.In addition, the impact of the number of generalized coordinates (n) on the accuracy of calculation results is studied.The specimens of section II in Table 1 are selected and their beam lengths are 4 m.The specific results are shown in Table 7 and Figure 6.According to Table 7 and Figure 6, the maximum error between M th and M FEM is only 3.61%, and calculation accuracy is high when n is increased to the convergence term.The maximum error between M th and M FEM is only 3.96% when n is 10, which is within the allowable error range, and is only 0.35% higher than that corresponding to the convergence term.When n is equal to 1 or 2, the error between M th and M FEM is greater than 17%.
. Convergence term refers to the value of n when the M th unchang as n increases.
Through the comparative analysis of finite elements, it can be seen that the critical moment of cantilever beam under the load cases studied in this paper can be calculated within the allowable error range when n is greater than or equal to 10, and when n is equal to 1 or 2, the error is too large to calculate.Therefore, the program in Figure 3 needs to be modified.The value of n in this program can directly take any value greater than 10, such as 100, without the need for a loop statement.It can greatly improve the efficiency of the procedure in Figure 3.In addition, the procedure of this paper refences the principles of GBTUL software in reference [48].The specific content of the code can be found in Appendix A. At the same time, a user graphical interface is developed for users.

Comparative Analysis between Matrix Equation and Equations Proposed in Other References
In past studies, many scholars have also proposed the equation of the M cr of cantilever beam, such as references [12,17].
The method of reference [12] is the Rayleigh-Ritz method.It used a set of nondimensional parameters and deduced the governing equations of cantilever beams.Reference [12] solved the governing equation by discretizing the beam model.However, reference [12] only describes its method and did not derive the explicit expression of the discrete matrix equation.Because the discrete matrix had too many physical quantities and was different from the expression in this paper, it will not be introduced in this paper.See Section 3 in reference [12] for details.(a) (b) Through the comparative analysis of finite elements, it can be seen that the critical moment of cantilever beam under the load cases studied in this paper can be calculated within the allowable error range when n is greater than or equal to 10, and when n is equal to 1 or 2, the error is too large to calculate.Therefore, the program in Figure 3 needs to be modified.The value of n in this program can directly take any value greater than 10, such as 100, without the need for a loop statement.It can greatly improve the efficiency of the procedure in Figure 3.In addition, the procedure of this paper refences the principles of GBTUL software in reference [48].The specific content of the code can be found in Appendix A. At the same time, a user graphical interface is developed for users.

Comparative Analysis between Matrix Equation and Equations Proposed in Other References
In past studies, many scholars have also proposed the equation of the Mcr of cantilever beam, such as references [12,17].
The method of reference [12] is the Rayleigh-Ritz method.It used a set of non-dimensional parameters and deduced the governing equations of cantilever beams.Reference [12] solved the governing equation by discretizing the beam model.However, reference [12] only describes its method and did not derive the explicit expression of the discrete matrix equation.Because the discrete matrix had too many physical quantities and was different from the expression in this paper, it will not be introduced in this paper.See Section 3 in reference [12] for details.
The method used in reference [17] is the finite difference approach.It established the governing equation of cantilever beams using the Galerkin method.The finite difference method was used to discretize the beam model and solve Mcr.
In order to compare the advantages and disadvantages of the equations in reference [12,17] and this paper, section I in Table 1 is selected for calculating Mcr.The lengths of the cantilever is 1.5, 2, 3 and 4 m, respectively, and the other properties of the beam are shown in Table 1.The concentrated load or uniformly distributed load is applied at the TF, SC and BF of the beam.The specific calculation results are shown in Figure 7.In Figure 7, the horizontal coordinates represent different load cases.When the numbers are 1, 2, 3 and 4, The method used in reference [17] is the finite difference approach.It established the governing equation of cantilever beams using the Galerkin method.The finite difference method was used to discretize the beam model and solve M cr .
In order to compare the advantages and disadvantages of the equations in reference [12,17] and this paper, section I in Table 1 is selected for calculating M cr .The lengths of the cantilever is 1.5, 2, 3 and 4 m, respectively, and the other properties of the beam are shown in Table 1.The concentrated load or uniformly distributed load is applied at the TF, SC and BF of the beam.The specific calculation results are shown in Figure 7.In Figure 7, the horizontal coordinates represent different load cases.When the numbers are 1, 2, 3 and 4, the position of the load is TF and the lengths of the beams are 1.5 m, 2 m, 3 m, and 4 m, respectively.When the numbers are 5, 6, 7 and 8, the position of the load is SC and the lengths of the beams are 1.5 m, 2 m, 3 m, and 4 m, respectively.When the numbers are 9, 10, 11, and 12, the position of the load is BF and the lengths of the beams are 1.5 m, 2 m, 3 m, and 4 m, respectively.
It can be seen from Figure 7 that the calculation results of Equation ( 34) are similar to those of reference [12], but are slightly different from those of reference [17].Through comparison and analysis, it can be concluded that Equation (34) in this paper has the following advantages: (i) Equation ( 34) has higher accuracy than reference [17]; (ii) Equation ( 34) can calculate M cr under CUDL, CL or uniformly distributed load, while the equation proposed in reference [12] can only calculate M cr under a single load and the equation proposed in reference [17] can only calculate M cr under CL and CUDL when the concentrated load is limited to the free end; (iii) Equation ( 34) considers the ratio of uniformly distributed load to concentrated load, the symmetry of the section, the action position of the concentrated load in the direction of the beam length and the vertical position of two loads.It considers very comprehensive factors and can calculate M cr in all cases.The specific comparison results can be seen in Tables 8-11.In Tables 8-11, " √ " indicates that it can be calculated and "×" indicates that it cannot be calculated.
distributed load to concentrated load, the symmetry of the section, the action position of the concentrated load in the direction of the beam length and the vertical position of two loads.It considers very comprehensive factors and can calculate Mcr in all cases.The specific comparison results can be seen in Tables 8-11.In Tables 8-11, "√" indicates that it can be calculated and "×" indicates that it cannot be calculated.

(a) (b)
Figure 7.Comparison of calculation results of three equations(ʺAndrade+2007" is reference [12] and "Ozbasaran+2013" is reference [17]): (a) calculation results under concentrated load at the free end; (b) calculation results under uniformly distributed load.Comparison of calculation results of three equations("Andrade+2007" is reference [12] and "Ozbasaran+2013" is reference [17]): (a) calculation results under concentrated load at the free end; (b) calculation results under uniformly distributed load.When n is greater than 1, Equation ( 34) can be transformed into a 2n degree equation about ∼ M 0 , which cannot be easily solved by mathematical methods.The following derivation process reveals the fundamental reasons for determining the highest order of ∼ M 0 in Equation (34).
The following buckling deformation functions are used.
According to Equation ( 21), a matrix equation of n + r order similar to Equation (34) can still be obtained: where the physical meaning and expression of R 0 , S 0 , T 0 , Q 0 , R 1 , S 1 , T 1 and Q 1 are the same as before.
According to Equations ( 26) and ( 28), the coefficient determinant of independent parameter A n (n = 1, 2, 3 . . . ) and B n (n = 1, 2, 3 . . . ) can be expressed concisely as follows: If n is greater than or equal to r, according to Equation (40) and the definition of determinant, that is, the determinant of order n is the algebraic sum of the products of n elements taken from different rows and columns, the operation result of Equation ( 40) is a 2r degree polynomial approximating ∼ M 0 .If n is than r, the operation result of Equation ( 40) is a n + r degree polynomial approximating ∼ M 0 .Therefore, only when r is equal to 1 and n is greater than or equal to r can Equation (40) be a quadratic polynomial approximating ∼ M 0 , and then ∼ M 0 can be accurately expressed by mathematical analytic expression.
Figure 8 shows the twist angle of each point along the z-axis of the beam when λ = 1, a Py = a qy and the torsional stiffness coefficient K for each is 0.5 and 1, 2.5, respectively.
Appl.Sci.2023, 13, x FOR PEER REVIEW 18 of 28 Therefore, only when r is equal to 1 and n is greater than or equal to r can Equation (40) be a quadratic polynomial approximating  , and then  can be accurately expressed by mathematical analytic expression.
Figure 8 shows the twist angle of each point along the z-axis of the beam when λ = 1, aP y = aq y and the torsional stiffness coefficient K for each is 0.5 and 1, 2.5, respectively.
It can be seen from Figure 8 that the twist buckling mode of a cantilever beam is related to the torsional stiffness coefficient K.However, it is difficult to find a twist angle function containing only one generalized coordinate to uniformly describe the twist buckling mode of the cantilever beam according to the results of Figure 8 and the existing literature.Therefore, only when the number of terms of the twist angle function is greater than or equal to 10 can the twist angle deformation of the cantilever beam be accurately described according to the conclusions in Table 7 and Figure 6.However, it has been proven that the analytical expression of  can be obtained accurately only when r equals 1, that is, the number of terms of twist angle function takes as 1.The contradiction formed before and after means that the Mcr of the cantilever beam cannot be accurately expressed by the analytical expression.Therefore, this paper must use the numerical calculation method to obtain the Mcr's expression.

Derivation of Closed-Form Solution of Mcr Based on Twist Angle Function
With the assumption that torsional rotation is small, the second derivative of u with respect to z can be written as given below from [9].
The moment distribution function   , the twist angle function φ(z) and Mmax are defined as follows.It can be seen from Figure 8 that the twist buckling mode of a cantilever beam is related to the torsional stiffness coefficient K.However, it is difficult to find a twist angle function containing only one generalized coordinate to uniformly describe the twist buckling mode of the cantilever beam according to the results of Figure 8 and the existing literature.Therefore, only when the number of terms of the twist angle function is greater than or equal to 10 can the twist angle deformation of the cantilever beam be accurately described according to the conclusions in Table 7 and Figure 6.However, it has been proven that the analytical expression of ∼ M 0 can be obtained accurately only when r equals 1, that is, the number of terms of twist angle function takes as 1.The contradiction formed before and after means that the M cr of the cantilever beam cannot be accurately expressed by the analytical expression.Therefore, this paper must use the numerical calculation method to obtain the M cr 's expression.

Derivation of Closed-Form Solution of M cr Based on Twist Angle Function
With the assumption that torsional rotation is small, the second derivative of u with respect to z can be written as given below from [9].
The moment distribution function κ(z), the twist angle function ϕ(z) and M max are defined as follows.
where M 0 is the moment at the fixed end of the cantilever beam.The value of M 0 is less than 0. ϕ 0 (z) is the basic function of twist angle function ϕ(z) of the cantilever, and B is an undetermined coefficient.According to the conclusions obtained in Table 7 and Figure 6, ϕ 0 (z) could use the following formula.
Some physical quantities are defined as follows.
Then Equation( 1) can be written as follows.
a q y q y B 2 r 5 + 1 2 a P y P y B 2 r 6 (52) M cr 's expression can be obtained from ∂Π ∂B = 0.
When the load case is CUDL, α = 1 and a Py = a qy are considered.Some physical quantities are defined as follows.
R 1 = 2r 1 (54) The expression of M cr can be written as follows.

Practical Calculation Formula of M cr
For the load case of CL, this paper assumes that α = 0 and calculates the M cr 's value for some common load cases when λ is equal to 1/3, 1/2, 2/3 and 1, respectively.
For the load case of CUDL, since the vertical position of the concentrated load (a Py ) and the uniform load (a qy ) on the beam are usually the same under actual load cases, this paper assumes that α = 1 and a Py = a qy and calculates the M cr 's value for some common load cases when λ is equal to 1/3, 1/2, 2/3 and 1, respectively.
It can be seen from Figure 8 that the twist angle mode of the cantilever beam during lateral-torsional buckling is mainly related to the torsional stiffness coefficient K. Since most of the commonly used values of K are between 0.3 and 2.5, this paper takes 15 different values of K between 0.3 and 2.5 and calculates B n by substituting relevant parameters into the MATLAB program compiled according to Equation (34) so as to obtain the specific expression of ϕ 0 (z) in Equation (45).Then the specific values of R 1 ~R5 can be obtained as shown in Tables 12 and 13.According to reference [49], Wagner's coefficient has only a slight impact on the buckling mode.Therefore, the conclusions in Tables 12 and 13 are applicable to both uniaxial symmetric and biaxial symmetric I-beams.Corresponding R 1 ~R5 's values can be obtained by interpolation when K has other values.Therefore, the M cr under CL or CUDL can be solved through Tables 12 and 13.

Example Verification
Section I, Section II that reinforced bottom flange and their properties in Table 1 are selected for calculating M cr .The length of cantilever is 4 m.It is assumed that the material is elastic and the component is free from defects.M cr is obtained according to Equation (59) and the finite element analysis respectively when λ takes different values.
Abaqus software is used to generate finite element models of the cantilever beams.The Finite Element Method (FEM) is applied to this software.The models of cantilever beams were generated in the CAE module and calculations were carried out in the STANDARD module [50].Cantilever beams are divided to 13,000 finite elements.S8R5 shell elements are used in this model and results are obtained by eigenvalue.In order to keep the original shape of the section and better conform to the rigid perimeter assumption, one stiffener is set for each l/6 in the length direction and a total of 5 stiffeners are set.In order not to create unnecessary constraints on other deformations of the cantilever beams, the stiffener is only rigidly connected with the web, and only couples with the in-plane displacement of top and bottom flanges.Tables 14 and 15 show that the results obtained by Equation (59) and Abaqus are nearly identical at SC.For TF and BF, difference between Abaqus results and Equation (59) results increases as λ decreases.Andrade and Camotim [13] explained this phenomenon.They verified that 1D frame element assumption loses its validity as the beam length decreases.And the situation of λ decreases is similar to this situation.
Table 12.Specific values of R 1 , R 2 , R 3 , R 4 and R 5 under CL.

Conclusions
A novel numerical method and a new closed-form solution for the elastic LTB analysis of cantilever steel beams have been successfully developed and validated in this paper.The numerical method was developed based on postulated buckling deformation functions and the Rayleigh-Ritz method.The closed-form solution was developed based on the assumption that the torsion angle is small and based on the numerical method.The values of M cr predicted by the proposed numerical method are sufficiently comparable to the data from existing literature and the following findings can be drawn: (1) M cr and the generalized coordinates A n and B n in the buckling deformation functions of cantilever I-beams under CL and CUDL can be accurately calculated by the numerical method proposed in this paper.Compared with the experimental results and the finite element results of relevant reference, the maximum of the M th /M ex and M th /M FEM ratio is 1.07, and the minimum error is 0.92.The average of all the ratios is 0.98 with a standard deviation of 0.13.(2) The safety, accuracy and economic aspects of the prediction of the numerical method are assessed using DPC.The results show that the numerical method can effectively predict M cr of cantilever steel beams under CL and CUDL.(3) The number of generalized coordinates determines the accuracy of M cr .The number of generalized coordinates is defined as n.The obtained errors of M cr do not exceed 3.96% when n is greater than or equal to 10.The errors of M cr exceed 17% when n is equal to 1 or 2. (4) The degree of M cr in Equation (34) varies with the number of generalized coordinates of twist angle function.When the number of generalized coordinates of twist angle function is n, the degree of M cr in Equation ( 34) is 2n.(5) The torsional buckling mode of the cantilever steel beam is mainly determined by K, which is the torsional stiffness coefficient.
However, in order to provide reference for the calculation of the M cr of cantilever steel beams matching the specification of ANSI/AISC 360-16, a large number of numerical calculations and statistical evaluations should be performed to determine the expression of C b , which is the lateral-torsional buckling modification factor.In addition, further investigations are needed on the LTB performance of beams composed of other materials that are not the isotropic material, such as Glass Fiber Reinforced Polymer (GFRP), laminate composite material and pultruded fiber reinforced polymer (PFRP) material, which are increasingly used in beams.Wagner's coefficient u, ϕ the lateral deflection and twist angle of section u(z), ϕ(z), ϕ 0 (z) the buckling deformation functions ϕ(λl)

Nomenclature
the twist angle at the point where the transverse concentrated load acts during the lateral-torsional buckling of the beam M x the moment about the x-axis α the ratio of the product of q y multiplied by l to P y h the distance between the centroids of the top and bottom flanges A n , B n generalized coordinates n, r the numbers of terms of buckling deformation functions M 0 , |M 0 | the moment when z = 0, the maximum absolute value of moment I 1 , I 2 the moment of inertia of the top flange around the y-axis, the moment of inertia of the bottom flange around the y-axis η the ratio of I 1 to I 2 M 0 the ratio of |M 0 | to the product of Euler critical load multiplied by h β x , a Py , a qy the ratio of β x to h, the ratio of a Py to h, the ratio of a qy to h K torsional stiffness coefficient R 0 , S 0 , T 0 , Q 0 , R 1 , S 1 , T 1 , Q

(∼M
0 ) is the smallest positive real number of the generalized eigenvalues of the matrices D 0 and D 1 .The matrix A B formed by the generalized coordinates is the generalized eigenvector corresponding to the smallest positive real number of the generalized eigenvalues of the matrices D 0 and D 1 .Our MATLAB program can be compiled to solve M cr and A B .

Figure 3 .
Figure 3. Flowchart of the procedure.

Figure 3 .
Figure 3. Flowchart of the procedure.

Figure 6 .
Figure 6.Comparison of predicted critical moment corresponding to different n values and MFEM: (a) comparison of cantilevers with section II with reinforced top flange; (b) comparison of cantilevers with section II with reinforced bottom flange.

Figure 6 .
Figure 6.Comparison of predicted critical moment corresponding to different n values and M FEM : (a) comparison of cantilevers with section II with reinforced top flange; (b) comparison of cantilevers with section II with reinforced bottom flange.

Figure 7 .
Figure 7.Comparison of calculation results of three equations("Andrade+2007" is reference[12] and "Ozbasaran+2013" is reference[17]): (a) calculation results under concentrated load at the free end; (b) calculation results under uniformly distributed load.

Figure 8 .
Figure 8.The torsion angle φ varies with the coefficient of K.

Figure 8 .
Figure 8.The torsion angle ϕ varies with the coefficient of K.

Table 1 .
Dimensions of specimens.

Table 1 .
Dimensions of specimens.

Table 3 .
The critical moment predicted by finite element analysis(M FEM ) and theoretically predicted (M th ) critical moment for Section I.

Table 4 .
The critical moment predicted by finite element analysis (M FEM ) and theoretically predicted (M th ) critical moment for Section II-reinforced top flange.

Table 5 .
The critical moment predicted by finite element analysis (M FEM ) and theoretically predicted (M th ) critical moment for Section II-reinforced bottom flange.

Table 6 .
Modified version of the Demerit Points Classification (DPC) criteria.

Table 6 .
Modified version of the Demerit Points Classification (DPC) criteria.

Table 7 .
Comparison of predicted critical moment corresponding to different n values and M FEM .
Note: Error

Table 8 .
Comparison of the ratio of uniformly distributed load to concentrated load.

Table 9 .
Comparison of cross-sectional symmetry.

Table 8 .
Comparison of the ratio of uniformly distributed load to concentrated load.

Table 9 .
Comparison of cross-sectional symmetry.

Table 10 .
Comparison of applicable load cases.

Table 11 .
Comparison of applicable vertical position of all loads.

Deduction of the Closed-Form Solution of M cr
4.1.Analysis of the Reason Why M cr Cannot Be Expressed Accurately by Analytic Expression

Table 14 .
Comparison of critical moment under CL.

Table 15 .
Comparison of critical moment under CUDL.

Table 15 .
Cont. : M cr1 is the calculated value by Equation (59).M cr2 is the calculated value by ABAQUS. Note Cartesian coordinates P y , q y concentrated load and uniformly distributed load a Py , a qy the difference between the ordinate of load and the ordinate of shear center λ the ratio of distance between the action point of Py and the fixed end to the cantilever's length E Young's modulus G shear modulus I ω , I t , I y warping constant, torsional constant and moment of inertia about y-axis β x