# Dimensionless Analytical Solution and New Design Formula for Lateral-Torsional Buckling of I-Beams under Linear Distributed Moment via Linear Stability Theory

Authors: Wen-Fu Zhang; Ying-Chun Liu; Ke-Shan Chen; Yun Deng
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

Even for the doubly symmetric I-beams under linear distributed moment, the design formulas given by codes of different countries are quite different. This paper will derive a dimensionless analytical solution via linear stability theory and propose a new design formula of the critical moment of the lateral-torsional buckling (LTB) of the simply supported I-beams under linear distributed moment. Firstly, the assumptions of linear stability theory are reviewed, the dispute concerning the LTB energy equation is introduced, and then the thinking of Plate-Beam Theory, which can be used to fully resolve the challenge presented by Ojalvo, is presented briefly; secondly, by introducing the new dimensionless coefficient of lateral deflection, the new dimensionless critical moment and Wagner’s coefficient are derived naturally from the total potential energy. With these independent parameters, the new dimensionless analytical buckling equation is obtained; thirdly, the convergence performance of the dimensionless analytical solution is discussed by numerical solutions and its correctness is verified by the numerical results given by ANSYS; finally, a new trilinear mathematical model is proposed as the benchmark of formulating the design formula and, with the help of 1stOpt software, the four coefficients used in the proposed dimensionless design formula are determined.


---

## Body

## 1. Introduction

Structural stability has always been a key point in the design of steel structures [1–6]. It is well known that an I-beam has various buckling phenomena such as local buckling, distortion buckling, and lateral-torsional buckling, in which the lateral-torsional buckling (LTB) is also known as the global buckling. This may occur if the applied service loads exceed its LTB critical moment of the I-beam and hence, in practice design, it is important for the designers to obtain the accurate solution of the critical moment as the upper limit of buckling strength of the steel I-beams. However, until now, even for the simple-supported I-beams with doubly symmetric sections subjected to unequal end moments, the design formulas given by the codes/specifications of different countries are quite different (Figure 1) that makes the designers feel very confused. In fact, such problems also confuse those involved in the formulation of the codes/specifications. In order to avoid greater controversy, there is a trend in the specification that is the initiative to delete the relevant provisions, such as the new version of the EC3 specification [7]. However, the deletion of the relevant provisions does not mean that such problems do not exist. Therefore, it is one of the objectives of this paper to provide a more accurate design formula based on the dimensionless analytic solution, which is also the main motivation of this paper.Figure 1
Discrepancies of EUMF for a simply supported beam under linear distributed moment.Due to the complexity of the LTB phenomenon of I-beams under nonuniform distributed moment, so far only approximate analytical solutions or numerical solutions have been published.Some approximate analytical solutions can be found in the classical text books such as Chajes [2], Chen and Atsuta [3], Trahair [4], Timoshenko and Gere [5], and Bleich [6].The numerical solutions can be obtained from finite difference method [2, 3, 6], finite integral method [4], and finite element method (FEM) [8]. More recently, Suryoatmono and Ho [9] and Serna et al. [10] have reported the results of elastic LTB load analysis using finite difference technique to solve the governing differential equation; Lim et al. [11] have conducted an extensive investigation into the elastic LTB of I-beams using both the Bubnov-Galerkin approach [1] and the finite element method [8]. In their papers, Vlasov’s differential equation governing [1] for lateral-torsional buckling of beams was used; Park et al. [12] use the finite element program developed by Lim et al. [11] to perform parametric studies to investigate the lateral-torsional buckling behavior of singly symmetric I-beams, in which I-beams are modeled by the beam element with two nodes per element and seven nodal degrees of freedom; Greiner and Lindner [13], Greiner et al. [14], and Mohri et al. [15] also use finite element method to present a quite complete numerical study on lateral-torsional buckling of beams; Kim et al. [16] adopt the variation of the total potential energy to derive the stability equations for the lateral buckling analysis of an arbitrarily laminated thin-walled composite beam, and the exact analytical stiffness matrix is presented based on power series expansions.It is noteworthy that, with the popularization and widespread use of finite element software, since the 2000s, some researchers have tried to use the general-purpose finite element software to simulate the buckling behavior or to assess the related design formula of the steel I-beams. For example, Mohri et al. [17] adopt the beam elements (B31OS) and shell elements (S8R5) of ABAQUS [18] to model the lateral buckling phenomena. The authors have found that shell elements could consider local buckling, section distortion, and the local effects of concentrated loads and boundary conditions, which were all ignored in the beam element theory; Dolamune Kankanamge and Mahendran [19] use ABAQUS [18] to investigate the LTB behavior of simply supported cold-formed steel lipped channel beams subjected to uniform bending. The authors have found that European design rules are conservative, while Australian/North American design rules are unconservative; Tong [20] uses ANSYS [21] to regress the design formula for the singly symmetric I-beams under linear distributed moment; Sweedan [22] adopts ANSYS to numerically investigate the lateral stability of cellular steel beams subjected to equal end moments, mid-span concentrated loads, and uniformly distributed loads; Serna et al. [10] use Cosmos/M27 [23] to simulate the lateral-torsional buckling of the I-beams.In theory, FEM can be used not only to do research works but also to formulate design formulas. However, when you really try to formulate the design formulas through a lot of FEM analysis, you will have to deal with some real trouble. Firstly, the finite element analysis must be a dimensional analysis [8, 18, 21, 23]; that is, the section size and span (or slenderness) of the steel beams must be specified in advance. Hence there is an enormous amount of workload in the early stage, such as the design of hundreds of model beams and the relevant data preparation in the finite element modeling; secondly, one has to spend a lot of time and effort in later stages to select from a large number of data outputs and determine which results are useful LTB results; finally, and most importantly, our experience shows that, because the meshing of FEM model will greatly affect the accuracy of finite element analysis, hence the corresponding FEM analysis results lack of precision consistency and regularity. This may be the reason why the design formulas of different countries will be quite different.On the contrary, not only does the dimensionless analytic solution have the advantages of small input data amount and fast calculation speed, but also the formula formulated from these dimensionless results is more general than those obtained from FEM analysis. To the author’s knowledge, research on the dimensionless analytical solution is scare and began in 1970s with the Anderson and Trahair’s works [24], followed by Kitipornchai et al. [25] and Kitipornchai and Wang [26, 27] to simulate the buckling behavior of singly symmetric I-beams and Tee beams under moment gradient. However, the dimensionless Wagner’s coefficient and beam parameter defined in their work lack definite physical and geometric meaning, and more importantly, none of the dimensionless parameters is independent; that is, they are not independent variables and cannot be used to get a reasonable design formula. Therefore, it is a very important task to develop new dimensionless analytic solution and use it for design purpose.According to the theory and discussion in [28], it can be proved that an exact dimensionless solution of the LTB problem of the simply supported I-beam can be obtained by the standard trigonometric series. Recently, based upon the author’s work [29, 30], we have carried out a series of research work on the dimensionless analytical solutions and design formulas of the critical moment of the LTB of the I-beams [31–34]. Recently, Zhang et al. [35] derived a dimensionless analytical solution and design formula of the critical moment of the LTB of the cantilever beam with tip lateral elastic brace under uniform and concentrated load.As a continuation of the work of the author, this study shall present a dimensionless analytical solution via linear stability theory and design formula of the critical moment of the LTB of the simply supported I-beams under linear distributed moment. In Section2, the basic assumptions of the traditional linear stability theory and the dispute concerning the LTB energy equation are introduced, and then the attempts to resolve the dispute and the Plate-Beam Theory are described briefly; in Section 3, with the modal trial function expressed by the trigonometric series, the energy variation method is adopted in the derivation of the dimensionless analytical solution of the critical moment. Moreover, the definition of dimensionless critical moment and its applicability is discussed; next, the 1st-order and 6th-order approximate analytical solution are presented in Section 4. In Section 5, the numerical solution and convergence performance of the dimensionless analytical solution are discussed and its correctness is verified by ANSYS software in Section 6. This part also presents a new technique to satisfy the rigid section hypothesis. Finally, a new mathematical model of the design formula is presented and, with the help of 1stOpt software, the four coefficients used in the proposed dimensionless design formula are determined.


## 2. Theory Background of LTB Energy Equation for I-Beams via Linear Stability Theory

As is known, for the slender (L≥10d), thin-walled (d≥10t) I-beams, the LTB (i.e., overall buckling) will generally control the final design. In this case, the simplified mechanical models, that is, differential equation model and energy equation model, of the LTB problem of the thin-walled beams can be established by introducing appropriate assumptions.For the sake of completeness, a brief review and comments of the basis of linear stability theory are presented hereafter.


### 2.1. Basic Assumptions of Traditional Linear Stability Theory

(1)
Rigid section hypothesis: that is, the cross section is very rigid that its original shape is retained during buckling. This is the well-known Vlasov’s rigid section hypothesis, which means that local and distortional buckling are excluded in the LTB analysis.
(2)
Euler-Bernoulli’s hypothesis: that is, there is no shear deformation in the middle surface of the cross section. This assumption is used to describe the in-plane deformation of a shell/plate.
(3)
Kirchhoff’s hypothesis: that is, the shear deformation in the planes normal to the middle surface of the cross section is small and can be neglected. This assumption is used to describe the out-of-plane deformation of a shell/plate.
(4)
The material is an ideal isotropic material and follows Hooke’s law.
(5)
Displacements and twist angle are assumed to be small enough.The first two hypotheses are first explicitly proposed by Vlasov [1], and the last two are implicitly found in the theoretical derivation of LTB problems [1, 20, 36–38].Invoking assumption(1), the following transverse displacements were proposed by Vlasov:(1)u-y,z=uz-y-y0θz;v-x,z=vz-x-x0θz,where uz, vz are the transverse displacements of the shear center, θz is the twist angle of the cross section, and x0, y0 are the coordinates of the shear center.Invoking assumption(2), that is, by setting the shear strain of the middle surface to zero, the following longitudinal displacement was deduced: (2)w-x,z=wz-xu′z-yv′z-ωθ′z,where ω is the “warping function” [36, pg. 11], but Vlasov called it the sectorial area [1, pg. 16].In Vlasov’s monograph, (2) is called the law of sectorial area. Obviously, since the condition about the shear strain of the middle surface is used in his derivation, (2) is only applicable to describing the longitudinal displacement of an arbitrary point in the middle surface rather than that outside the middle surface. As a result, the effect of St. Venant’s torsion could not be deduced correctly by the law of sectorial area. Therefore, in essence, (2) cannot be directly used to derive the energy equation for LTB problem of thin-walled structures.In short, even though Vlasov’s theory is such invaluable that almost all the literature on the theory of thin-walled bars will refer to his work, thus making him an outstanding figure among the pioneers in the field [36–38], his theory is incomplete, since in his work only the LTB differential equations were established by the fictitious load method, and no LTB energy equation was presented, since his proposed displacement field could not be used to deduce correct energy equation. This is precisely because of such basic defects and the abuse of the law of sectorial area, leading to a lot of controversy in future generations.


### 2.2. Dispute concerning the LTB Energy Equation

The first correct LTB energy equation is proposed by Bleich [6] in 1952 (called traditional LTB energy equation). However, his theoretical derivation was not rigorous. That is, it is theoretically incomplete. The defects include the following: (1) the longitudinal displacement of the cross section is not clearly defined; (2) the strain energy corresponding to the St. Venant torsion cannot be derived naturally at the same time and thus should be additionally added by invoking the St. Venant formula; (3) the Wagner effect caused by the transverse load cannot be obtained naturally, because only the longitudinal potential energy is considered in his derivation.In addition, both Vlasov and Bleich’s theories do not apply to the combined cross sections (e.g., open-closed/solid-open) and composite cross sections (e.g., steel-concrete).Because the LTB energy equation is the foundation of the modern mechanics analysis, such as finite element method and approximate analytic solution, extensive research has been carried out since the 1960s. The related theoretical derivation on the LTB energy equation for I-beams has been detailed in [2–5, 20, 36–38]. It is found that since the correctness of the LTB energy equation is often subject to a variety of academic challenges, most of these efforts were aimed at trying to explain the rationality of the following LTB energy equation proposed by Bleich, that is, (3)Π=12∫LEIyu′′2+EIωθ′′2+GJk+2Mxβxθ′2+2Mxu′′θ-qaqθ2dz,where u is the lateral deflection of the shear center, θ is the twist angle of the cross section, EIy is the minor axis flexural rigidity, GJk is the St. Venant torsion rigidity, EIω is the Vlasov warping torsion rigidity, Mx is the moment function about the major axis, L is the length of the beam, q is the distributed load acting in the vertical plane of the I-beams, aq is the load position parameter, which is obtained by subtracting the y-coordinate of the load acting point from that of the shear center, and βx is the Wagner coefficient, and it is defined as(4)βx=12Ix∫Ayx2+y2dA-y0,where x, y is the coordinates of an arbitrary point with respect to the origin (the centroid of the cross section) and y0 is the coordinate of the shear center of the cross section.The most fatal challenge comes from the question of the validity of the Wagner hypothesis, which is a part of the theoretical basis for the traditional LTB theory. Since 1981, this problem has been continuously questioned and challenged by Ojalvo [38–41] and considerable ongoing efforts [42–47] have been made by some famous scholars, such as Trahair, Chen, and Knag et al., to confirm the validity of the Wagner hypothesis in the framework of Vlasov theory. Although these in-depth series of discussions have led to a further understanding of the Wagner hypothesis and the related LTB theory, the basis of the argument is still relatively inadequate and did not get everyone’s recognition. Concerning this, Ziemian in his book [48] concludes that “it should be noted that a challenge to part of the theory of monosymmetric beams, called the Wagner hypothesis, was presented by Ojalvo (1981). At the present time (2009), this challenge still has not been fully resolved.”


### 2.3. Attempts to Resolve the Dispute and the Plate-Beam Theory

From the appearance point of view, the problem of the Wagner coefficient is only reflected in process of deriving the load potential; hence some attempts have been made to improve the Vlasov’s displacement filed in order to obtain the correct Wagner coefficient. For example, based upon a more complex geometric analysis, a longitudinal displacement with two additional nonlinear terms was proposed by Trahair (see [4, pg. 384]):(5)w-x,z=w-xu′-yv′-ωθ′+-xv′θ+yu′θ.In addition, even the concept of the rotation matrix, which was suitable for the nonlinear analysis with large rotation [49], was introduced into the linear LTB analysis, and the transverse displacements with more additional nonlinear terms were proposed by Pi et al. [50] as follows:(6)u-y,z=uz-yθz-xu′2+θ22-yu′v′2-ωu′θ′,v-x,z=v-xθz-xu′v′2-yv′2+θ22-ωv′θ′,w-x,z=w-xu′-yv′-ωθ′+-xv′θ+yu′θ.Obviously, incorrect linear strain energy will be obtained if the above displacements of Trahair and Pi et al. were used. Therefore, the above improvements are cumbersome and unnecessary.In fact, The LTB problem is similar to the vibration problem of the prestressed concrete beam, and hence this problem belongs to the category of linear eigenvalue problem in essence. In modern mechanics terms, this problem is a small displacement finite strain problem, namely, in such analysis, assumptions(3) and (4) hold, thus increasing the nonlinear term contradiction with the above basic assumptions. In addition, the second variation criterion was used in their derivation, which has been discarded in modern variational principle [51–54]. These problems were also criticized by Ojalvo [41] as follows: “...the underlying theorem on which the derivation must be based (the theorem of stationary potential energy) is misused: Components of a finite strain tensor are used in the strain energy expression when infinitesimal strain expressions are clearly called for by the theorem.”In the author’s opinion, so far the problem of how to derive the linear strain energy correctly based only upon the commonly used engineering mechanics theory has not been resolved in both the traditional theory and Ojalvo’s new theory [38], since, in the whole process of the theoretical debates, this problem has been deliberately avoided by both sides of the debaters.In fact, there is a common theoretical basis in both Bleich’s theory and Ojalvo’s new theory; that is, they all tried to use engineering theory (beam or rod) to derive their LTB energy equation. For example, Navier’s plane section hypothesis is used in Bleich’s derivation of the linear strain energy. But these thoughts were rarely mentioned in the published LTB literature.This study considers a simply supported steel I-beam under linear distributed moment as shown in Figure2. Obviously, this I-beam is an open thin-walled beam composed of three flat plates (i.e., two flanges and one web plate). Therefore, how to more rationally describe in-plane and the out-of-plane deformation of a flat plate by the classical engineering mechanics is the essence of solving the debate of the LTB problem reasonably. Based on such a basic understanding, the author created the Plate-Beam Theory [30, 55–58] recently.Figure 2
Calculation diagram of simply supported beam under linear distributed moment.In this new theory, Vlasov’s hypotheses are preserved, but the warping function and the Vlasov’s longitudinal displacement are discarded. The longitudinal displacement, geometrical equation, strain energy, and initial stress potential energy of the in-plane bending and the out-plane bending and torsion of each flat plate are described by the Euler beam theory and Kirchhoff-Plate theory, respectively (Figure3). It is shown that this new theory combines the advantages of Vlasov and Bleich’s theory and applies not only to open [56]/closed cross sections [57] but also to open-closed/solid-open cross sections [58].Figure 3
Illustration of section deformation assumption and plate deformation analogy.In short, the Plate-Beam Theory provides a new way to resolve the dispute concerning the Wagner hypothesis completely, in which the rational elements of Bleich’s, Vlasov’s, and Ojalvo’s thoughts have been inherited and developed, while the warping function of Vlasov was abandoned. Furthermore, the correctness of the traditional LTB energy equation, such as (3), is reaffirmed via the commonly used plate and beam theory.


## 2.1. Basic Assumptions of Traditional Linear Stability Theory

(1)
Rigid section hypothesis: that is, the cross section is very rigid that its original shape is retained during buckling. This is the well-known Vlasov’s rigid section hypothesis, which means that local and distortional buckling are excluded in the LTB analysis.
(2)
Euler-Bernoulli’s hypothesis: that is, there is no shear deformation in the middle surface of the cross section. This assumption is used to describe the in-plane deformation of a shell/plate.
(3)
Kirchhoff’s hypothesis: that is, the shear deformation in the planes normal to the middle surface of the cross section is small and can be neglected. This assumption is used to describe the out-of-plane deformation of a shell/plate.
(4)
The material is an ideal isotropic material and follows Hooke’s law.
(5)
Displacements and twist angle are assumed to be small enough.The first two hypotheses are first explicitly proposed by Vlasov [1], and the last two are implicitly found in the theoretical derivation of LTB problems [1, 20, 36–38].Invoking assumption(1), the following transverse displacements were proposed by Vlasov:(1)u-y,z=uz-y-y0θz;v-x,z=vz-x-x0θz,where uz, vz are the transverse displacements of the shear center, θz is the twist angle of the cross section, and x0, y0 are the coordinates of the shear center.Invoking assumption(2), that is, by setting the shear strain of the middle surface to zero, the following longitudinal displacement was deduced: (2)w-x,z=wz-xu′z-yv′z-ωθ′z,where ω is the “warping function” [36, pg. 11], but Vlasov called it the sectorial area [1, pg. 16].In Vlasov’s monograph, (2) is called the law of sectorial area. Obviously, since the condition about the shear strain of the middle surface is used in his derivation, (2) is only applicable to describing the longitudinal displacement of an arbitrary point in the middle surface rather than that outside the middle surface. As a result, the effect of St. Venant’s torsion could not be deduced correctly by the law of sectorial area. Therefore, in essence, (2) cannot be directly used to derive the energy equation for LTB problem of thin-walled structures.In short, even though Vlasov’s theory is such invaluable that almost all the literature on the theory of thin-walled bars will refer to his work, thus making him an outstanding figure among the pioneers in the field [36–38], his theory is incomplete, since in his work only the LTB differential equations were established by the fictitious load method, and no LTB energy equation was presented, since his proposed displacement field could not be used to deduce correct energy equation. This is precisely because of such basic defects and the abuse of the law of sectorial area, leading to a lot of controversy in future generations.


## 2.2. Dispute concerning the LTB Energy Equation

The first correct LTB energy equation is proposed by Bleich [6] in 1952 (called traditional LTB energy equation). However, his theoretical derivation was not rigorous. That is, it is theoretically incomplete. The defects include the following: (1) the longitudinal displacement of the cross section is not clearly defined; (2) the strain energy corresponding to the St. Venant torsion cannot be derived naturally at the same time and thus should be additionally added by invoking the St. Venant formula; (3) the Wagner effect caused by the transverse load cannot be obtained naturally, because only the longitudinal potential energy is considered in his derivation.In addition, both Vlasov and Bleich’s theories do not apply to the combined cross sections (e.g., open-closed/solid-open) and composite cross sections (e.g., steel-concrete).Because the LTB energy equation is the foundation of the modern mechanics analysis, such as finite element method and approximate analytic solution, extensive research has been carried out since the 1960s. The related theoretical derivation on the LTB energy equation for I-beams has been detailed in [2–5, 20, 36–38]. It is found that since the correctness of the LTB energy equation is often subject to a variety of academic challenges, most of these efforts were aimed at trying to explain the rationality of the following LTB energy equation proposed by Bleich, that is, (3)Π=12∫LEIyu′′2+EIωθ′′2+GJk+2Mxβxθ′2+2Mxu′′θ-qaqθ2dz,where u is the lateral deflection of the shear center, θ is the twist angle of the cross section, EIy is the minor axis flexural rigidity, GJk is the St. Venant torsion rigidity, EIω is the Vlasov warping torsion rigidity, Mx is the moment function about the major axis, L is the length of the beam, q is the distributed load acting in the vertical plane of the I-beams, aq is the load position parameter, which is obtained by subtracting the y-coordinate of the load acting point from that of the shear center, and βx is the Wagner coefficient, and it is defined as(4)βx=12Ix∫Ayx2+y2dA-y0,where x, y is the coordinates of an arbitrary point with respect to the origin (the centroid of the cross section) and y0 is the coordinate of the shear center of the cross section.The most fatal challenge comes from the question of the validity of the Wagner hypothesis, which is a part of the theoretical basis for the traditional LTB theory. Since 1981, this problem has been continuously questioned and challenged by Ojalvo [38–41] and considerable ongoing efforts [42–47] have been made by some famous scholars, such as Trahair, Chen, and Knag et al., to confirm the validity of the Wagner hypothesis in the framework of Vlasov theory. Although these in-depth series of discussions have led to a further understanding of the Wagner hypothesis and the related LTB theory, the basis of the argument is still relatively inadequate and did not get everyone’s recognition. Concerning this, Ziemian in his book [48] concludes that “it should be noted that a challenge to part of the theory of monosymmetric beams, called the Wagner hypothesis, was presented by Ojalvo (1981). At the present time (2009), this challenge still has not been fully resolved.”


## 2.3. Attempts to Resolve the Dispute and the Plate-Beam Theory

From the appearance point of view, the problem of the Wagner coefficient is only reflected in process of deriving the load potential; hence some attempts have been made to improve the Vlasov’s displacement filed in order to obtain the correct Wagner coefficient. For example, based upon a more complex geometric analysis, a longitudinal displacement with two additional nonlinear terms was proposed by Trahair (see [4, pg. 384]):(5)w-x,z=w-xu′-yv′-ωθ′+-xv′θ+yu′θ.In addition, even the concept of the rotation matrix, which was suitable for the nonlinear analysis with large rotation [49], was introduced into the linear LTB analysis, and the transverse displacements with more additional nonlinear terms were proposed by Pi et al. [50] as follows:(6)u-y,z=uz-yθz-xu′2+θ22-yu′v′2-ωu′θ′,v-x,z=v-xθz-xu′v′2-yv′2+θ22-ωv′θ′,w-x,z=w-xu′-yv′-ωθ′+-xv′θ+yu′θ.Obviously, incorrect linear strain energy will be obtained if the above displacements of Trahair and Pi et al. were used. Therefore, the above improvements are cumbersome and unnecessary.In fact, The LTB problem is similar to the vibration problem of the prestressed concrete beam, and hence this problem belongs to the category of linear eigenvalue problem in essence. In modern mechanics terms, this problem is a small displacement finite strain problem, namely, in such analysis, assumptions(3) and (4) hold, thus increasing the nonlinear term contradiction with the above basic assumptions. In addition, the second variation criterion was used in their derivation, which has been discarded in modern variational principle [51–54]. These problems were also criticized by Ojalvo [41] as follows: “...the underlying theorem on which the derivation must be based (the theorem of stationary potential energy) is misused: Components of a finite strain tensor are used in the strain energy expression when infinitesimal strain expressions are clearly called for by the theorem.”In the author’s opinion, so far the problem of how to derive the linear strain energy correctly based only upon the commonly used engineering mechanics theory has not been resolved in both the traditional theory and Ojalvo’s new theory [38], since, in the whole process of the theoretical debates, this problem has been deliberately avoided by both sides of the debaters.In fact, there is a common theoretical basis in both Bleich’s theory and Ojalvo’s new theory; that is, they all tried to use engineering theory (beam or rod) to derive their LTB energy equation. For example, Navier’s plane section hypothesis is used in Bleich’s derivation of the linear strain energy. But these thoughts were rarely mentioned in the published LTB literature.This study considers a simply supported steel I-beam under linear distributed moment as shown in Figure2. Obviously, this I-beam is an open thin-walled beam composed of three flat plates (i.e., two flanges and one web plate). Therefore, how to more rationally describe in-plane and the out-of-plane deformation of a flat plate by the classical engineering mechanics is the essence of solving the debate of the LTB problem reasonably. Based on such a basic understanding, the author created the Plate-Beam Theory [30, 55–58] recently.Figure 2
Calculation diagram of simply supported beam under linear distributed moment.In this new theory, Vlasov’s hypotheses are preserved, but the warping function and the Vlasov’s longitudinal displacement are discarded. The longitudinal displacement, geometrical equation, strain energy, and initial stress potential energy of the in-plane bending and the out-plane bending and torsion of each flat plate are described by the Euler beam theory and Kirchhoff-Plate theory, respectively (Figure3). It is shown that this new theory combines the advantages of Vlasov and Bleich’s theory and applies not only to open [56]/closed cross sections [57] but also to open-closed/solid-open cross sections [58].Figure 3
Illustration of section deformation assumption and plate deformation analogy.In short, the Plate-Beam Theory provides a new way to resolve the dispute concerning the Wagner hypothesis completely, in which the rational elements of Bleich’s, Vlasov’s, and Ojalvo’s thoughts have been inherited and developed, while the warping function of Vlasov was abandoned. Furthermore, the correctness of the traditional LTB energy equation, such as (3), is reaffirmed via the commonly used plate and beam theory.


## 3. Dimensionless Analytical Solution

### 3.1. Moment Function and Modal Trial Functions

For the simply supported I-beam shown in Figure2, it is easy to obtain the moment function as follows:(7)Mxz=kM1+1-kM1zL,k∈-1,+1,where k=M2/M1 is the end moment ratio and M2, M1 are the left and right end moment, respectively, and M2≤M1.For the case of the simply supported I-beams, the following modal trial functions are proposed for the lateral deflectionuz and twist angle θz [29–35], which can be expressed as the sum of trigonometric series as follows:(8)uz=h∑n=1∞Amsin⁡mπzL;m=1,2,3,…,∞,(9)θz=∑n=1∞Bnsin⁡nπzL;n=1,2,3,…,∞,where Am and Bn are the undetermined dimensionless coefficients for the lateral deflection u(z) and twist angle θ(z), respectively; h is the distance between the centroid of top and bottom flange.From the perspective of structural dynamics,Am and Bn can be regarded as the generalized coordinates or generalized degrees of freedom. However, their dimensions are different, so in order to remove the dimension of Am so that Am and Bn are both undetermined dimensionless coefficients, we introduce h into the expression of u(z) in this work. This idea is proposed by the author, which is completely different from other previous research [25–27] as discussed in Section 3.4, and it has been used in LTB analysis of various steel beams [29–35].Obviously, the above trial functions are orthogonal and satisfy the following geometric boundary conditions of the simply supported beams; that is,(10)u0=u′′0=0;uL=u′′L=0,θ0=θ′′0=0;θL=θ′′L=0.Moreover, according to the theory and discussions given by literature [28], it can be proved that trigonometric series is the exact solution of the modal trial function for the problem shown in Figure 2, and hence the exact analytical expressions of the total potential energy and the one for buckling equation can be derived directly.


### 3.2. Dimensionless Total Potential Energy

Substituting the moment function (7) and the modal trial functions (8)-(9) into the potential energy functional (3) and integrating over the beam length L, one can get the integration results of the total potential energy.Without loss of generality, the following new dimensionless parameters [29–35] are introduced:(11)M~0=M1π2EIy/L2h;β~x=βxh;η=I1I2;K=π2EIωGJkL2,where M~0 is the dimensionless end moment; h is the distance between the centroid of the top and bottom flange; K is the dimensionless torsional stiffness parameters; η is the ratio of the moment inertia about minor axes of the top flange to that of the bottom flange; I1, I2 are the 2nd moment of inertia about minor axes of the top and bottom flange, respectively; β~x is the dimensionless Wagner’s coefficient.In the following derivation, the energy equation of (3) is first multiplied by a common factor L3/h2EIy, and then the substitution relationships of (11) are introduced. For the sake of clarity, the final integration results are listed as follows in three parts:(12)Π1=12∫0LEIyu′′2+EIωθ′2+GJkθ′2dz=∑m=1∞Am2m4π44+∑n=1∞Bn21+n2K2n2π4η4K21+η2,Π2=12∫L2Mxβxθ′2dz=M~0∑n=1∞Bn21+km2π4β~x4+M~0∑s=1∞∑r=1r≠s∞BsBr4-1+kπ2rsr2+s2β~xr2-s22;s±r=even,Π3=12∫L2Mxu′′θ′dz=-M~0∑m=1∞AmBm1+km2π44-M~0∑s=1∞∑r=1r≠s∞AsBr4-1+kπ2rs3r2-s22;s±r=even.The total potential energy of this problem is the sum of the above equations; that is,(13)Π=∑i=13Πi.


### 3.3. Dimensionless Analytical Solution of Buckling Equation

According to the principle of energy variation method, the following stationary conditions that minimize the total potential energy are required:(14)∂Π∂Am=0,m=1,2,3,…,∞;∂Π∂Bn=0,n=1,2,3,…,∞.Note that, in the above equations, m and n are arbitrary and independent.From the condition∂Π/∂Am=0, we can get the equilibrium equation with respect to the lateral deflection as follows:(15)∑m=1∞Amm4π42-M~0∑m=1∞Bm1+km2π44-M~0∑m=1∞∑r=1r≠mm±r=even∞Br4-1+kπ2rm3r2-m22=0.For each fixedm (corresponding to Line m), herein, one homogeneous equation corresponding to Am, Bm, and Br will be obtained.With the knowledge of the internal data structure of the above infinite series equation and using the following the notations listed as follows:(16)A=A1A2A3⋯T;B=B1B2B3⋯T,we can rewrite (15) in the following matrix form:(17)R0A+S0B=M~0R1A+S1B,(18)Rm,m0=m4π42,m=1,2,…,∞,Rs,r0=0,r≠s,s=1,2,…,∞,r=1,2,…,∞,(19)Ss,r0=0,s=1,2,…,∞,r=1,2,…,∞,(20)Rs,r1=0;s=1,2,…,∞,r=1,2,…,∞,(21)Sm,m1=1+km2π44,m=1,2,…,∞,Ss,r1=4-1+kπ2rs3r2-s22,s≠r,s±r=even,s=1,2,…,∞,r=1,2,…,∞,where A and B are the column vectors constituted by all undetermined coefficients Am and Bm, respectively; R0 and R1 are the coefficient matrixes for A in the left and right hand of the above equilibrium equation, respectively; S0 and S1 are the coefficient matrixes for B in the left and right hand of the above equilibrium equation, respectively.Similarly, from the condition∂Π/∂Bn=0, we can get the equilibrium equation with respect to the twist rotation as follows:(22)-M~0∑n=1∞An1+kn2π44-M~0∑n=1∞∑r=1r≠nn±r=even∞Ar4-1+kπ2nr3n2-r22+M~0∑n=1∞Bn1+kn2π4β~x2+∑n=1∞Bn1+n2K2n2π4η2K21+η2+M~0∑n=1∞∑r=1r≠nn±r=even∞Br4-1+kπ2rnr2+n2β~xr2-n22=0.With the notations in (16), the above equation can also be arranged in the following matrix form:(23)T0A+Q0B=M~0T1A+Q1B,(24)Ts,r0=0,s=1,2,…,∞,r=1,2,…,∞,(25)Ts,r1=Sr,s1,s=1,2,…,∞,r=1,2,…,∞,(26)Qm,m0=1+m2K2m2π4η2K21+η2,m=1,2,…,∞,Qs,r0=0,r≠s,s=1,2,…,∞,r=1,2,…,∞,(27)Qm,m1=-1+km2π4β~x2,m=1,2,…,∞,Qs,r1=-4-1+kπ2rsr2+s2β~xr2-s22,r≠s,s±r=even,s=1,2,…,∞,r=1,2,…,∞,where T0 and T1 are the coefficient matrixes for A in the left and right hand of the above equilibrium equation, respectively; Q0 and Q1 are the coefficient matrixes for B in the left and right hand of the above equilibrium equation, respectively.If (17) and (23) are combined together, the buckling equation can be expressed as follows in the form of block matrix:(28)R0S0T0Q0AB=M~0R1S1T1Q1AB,where M~0 is the dimensionless critical moment, that is, the dimensionless buckling strength; R, S, T, and Q are the subblock coefficient matrixes constructed from analytical expressions included in (18)–(21) and (24)–(27), respectively.Equation (28) is the dimensionless analytical buckling equation of the LTB problem of a simply supported steel beam under linear distributed moment obtained in the present work.Obviously, from the mechanical point of view, the matrix at the left-hand side of the equation represents the linear stiffness matrix of the steel girder and is independent of the load, while the matrix on the left side of the equation depends on the load pattern, which in the present case can be considered as a geometric stiffness matrix due to the linear distribution of moments.It is noted that the dimensionless analytical solution of the buckling equation presented in the paper not only has a clear data structure and but also is easy to program and solve. In addition, unlike the FEM solutions published in literature, this solution can easily be verified by anyone.


### 3.4. Remarks on Definition of Dimensionless Critical Moment

There are different ways to define the dimensionless critical moment, in which the well-known definition is the one proposed by Kitipornchai et al. [25] and Kitipornchai and Wang [26, 27] and can be expressed as follows:(29)M~cr=McrLEIyGJk.Corresponding to this definition, the following dimensionless Wagner’s coefficient, the load position parameter, and the beam parameter are defined by Kitipornchai:(30)δ=βxLEIyGJk;ε=aLEIyGJk;K~=π2EIyh2GJkL2=πhLEIyGJk.Obviously, all these parameters lack clear physical and geometric meanings. More importantly, since there is a common factorEIy/GJk in (30), none of the parameters defined in (30) is independent; that is, they are not independent variables; thus it is difficult to use them to regress a rational design formula.Instead, the following dimensionless critical moment is used in this paper:(31)M~cr=Mcrπ2EIy/L2h=McrhPEy.This definition is first proposed by Professor Zhang [29] since 2008 and later used in the LTB analysis of various steel beams [30–35].It is the first finding in this paper that this definition can be derived naturally from the total potential energy based upon the new definition of displacement function, that is, (8); secondly, it has a clear physical meaning as shown in the last term in (31); that is, the product of Euler’s bending buckling load around the minor axis and the distance between the centroid of the top and bottom flange is used to normalize the actual critical moment.In addition, in the case of the LTB of a doubly symmetric beam under pure bending, the definition of (31) also has a definite geometric meaning; that is,(32)M~cr=McrhPEy=uθ=121+K-2.That is, the dimensionless critical moment is the dimensionless distance between the fixed axis of rotation and the shear center for a cross section. It is clear that, for the doubly symmetric beams, the fixed axis of rotation always lies beyond the tension flange due toK>0.Finally, and most importantly, the other dimensionless parameters defined in (11) also have definite physical meanings and are independent of each other and are therefore suitable for the formulation of a new design formula.


## 3.1. Moment Function and Modal Trial Functions

For the simply supported I-beam shown in Figure2, it is easy to obtain the moment function as follows:(7)Mxz=kM1+1-kM1zL,k∈-1,+1,where k=M2/M1 is the end moment ratio and M2, M1 are the left and right end moment, respectively, and M2≤M1.For the case of the simply supported I-beams, the following modal trial functions are proposed for the lateral deflectionuz and twist angle θz [29–35], which can be expressed as the sum of trigonometric series as follows:(8)uz=h∑n=1∞Amsin⁡mπzL;m=1,2,3,…,∞,(9)θz=∑n=1∞Bnsin⁡nπzL;n=1,2,3,…,∞,where Am and Bn are the undetermined dimensionless coefficients for the lateral deflection u(z) and twist angle θ(z), respectively; h is the distance between the centroid of top and bottom flange.From the perspective of structural dynamics,Am and Bn can be regarded as the generalized coordinates or generalized degrees of freedom. However, their dimensions are different, so in order to remove the dimension of Am so that Am and Bn are both undetermined dimensionless coefficients, we introduce h into the expression of u(z) in this work. This idea is proposed by the author, which is completely different from other previous research [25–27] as discussed in Section 3.4, and it has been used in LTB analysis of various steel beams [29–35].Obviously, the above trial functions are orthogonal and satisfy the following geometric boundary conditions of the simply supported beams; that is,(10)u0=u′′0=0;uL=u′′L=0,θ0=θ′′0=0;θL=θ′′L=0.Moreover, according to the theory and discussions given by literature [28], it can be proved that trigonometric series is the exact solution of the modal trial function for the problem shown in Figure 2, and hence the exact analytical expressions of the total potential energy and the one for buckling equation can be derived directly.


## 3.2. Dimensionless Total Potential Energy

Substituting the moment function (7) and the modal trial functions (8)-(9) into the potential energy functional (3) and integrating over the beam length L, one can get the integration results of the total potential energy.Without loss of generality, the following new dimensionless parameters [29–35] are introduced:(11)M~0=M1π2EIy/L2h;β~x=βxh;η=I1I2;K=π2EIωGJkL2,where M~0 is the dimensionless end moment; h is the distance between the centroid of the top and bottom flange; K is the dimensionless torsional stiffness parameters; η is the ratio of the moment inertia about minor axes of the top flange to that of the bottom flange; I1, I2 are the 2nd moment of inertia about minor axes of the top and bottom flange, respectively; β~x is the dimensionless Wagner’s coefficient.In the following derivation, the energy equation of (3) is first multiplied by a common factor L3/h2EIy, and then the substitution relationships of (11) are introduced. For the sake of clarity, the final integration results are listed as follows in three parts:(12)Π1=12∫0LEIyu′′2+EIωθ′2+GJkθ′2dz=∑m=1∞Am2m4π44+∑n=1∞Bn21+n2K2n2π4η4K21+η2,Π2=12∫L2Mxβxθ′2dz=M~0∑n=1∞Bn21+km2π4β~x4+M~0∑s=1∞∑r=1r≠s∞BsBr4-1+kπ2rsr2+s2β~xr2-s22;s±r=even,Π3=12∫L2Mxu′′θ′dz=-M~0∑m=1∞AmBm1+km2π44-M~0∑s=1∞∑r=1r≠s∞AsBr4-1+kπ2rs3r2-s22;s±r=even.The total potential energy of this problem is the sum of the above equations; that is,(13)Π=∑i=13Πi.


## 3.3. Dimensionless Analytical Solution of Buckling Equation

According to the principle of energy variation method, the following stationary conditions that minimize the total potential energy are required:(14)∂Π∂Am=0,m=1,2,3,…,∞;∂Π∂Bn=0,n=1,2,3,…,∞.Note that, in the above equations, m and n are arbitrary and independent.From the condition∂Π/∂Am=0, we can get the equilibrium equation with respect to the lateral deflection as follows:(15)∑m=1∞Amm4π42-M~0∑m=1∞Bm1+km2π44-M~0∑m=1∞∑r=1r≠mm±r=even∞Br4-1+kπ2rm3r2-m22=0.For each fixedm (corresponding to Line m), herein, one homogeneous equation corresponding to Am, Bm, and Br will be obtained.With the knowledge of the internal data structure of the above infinite series equation and using the following the notations listed as follows:(16)A=A1A2A3⋯T;B=B1B2B3⋯T,we can rewrite (15) in the following matrix form:(17)R0A+S0B=M~0R1A+S1B,(18)Rm,m0=m4π42,m=1,2,…,∞,Rs,r0=0,r≠s,s=1,2,…,∞,r=1,2,…,∞,(19)Ss,r0=0,s=1,2,…,∞,r=1,2,…,∞,(20)Rs,r1=0;s=1,2,…,∞,r=1,2,…,∞,(21)Sm,m1=1+km2π44,m=1,2,…,∞,Ss,r1=4-1+kπ2rs3r2-s22,s≠r,s±r=even,s=1,2,…,∞,r=1,2,…,∞,where A and B are the column vectors constituted by all undetermined coefficients Am and Bm, respectively; R0 and R1 are the coefficient matrixes for A in the left and right hand of the above equilibrium equation, respectively; S0 and S1 are the coefficient matrixes for B in the left and right hand of the above equilibrium equation, respectively.Similarly, from the condition∂Π/∂Bn=0, we can get the equilibrium equation with respect to the twist rotation as follows:(22)-M~0∑n=1∞An1+kn2π44-M~0∑n=1∞∑r=1r≠nn±r=even∞Ar4-1+kπ2nr3n2-r22+M~0∑n=1∞Bn1+kn2π4β~x2+∑n=1∞Bn1+n2K2n2π4η2K21+η2+M~0∑n=1∞∑r=1r≠nn±r=even∞Br4-1+kπ2rnr2+n2β~xr2-n22=0.With the notations in (16), the above equation can also be arranged in the following matrix form:(23)T0A+Q0B=M~0T1A+Q1B,(24)Ts,r0=0,s=1,2,…,∞,r=1,2,…,∞,(25)Ts,r1=Sr,s1,s=1,2,…,∞,r=1,2,…,∞,(26)Qm,m0=1+m2K2m2π4η2K21+η2,m=1,2,…,∞,Qs,r0=0,r≠s,s=1,2,…,∞,r=1,2,…,∞,(27)Qm,m1=-1+km2π4β~x2,m=1,2,…,∞,Qs,r1=-4-1+kπ2rsr2+s2β~xr2-s22,r≠s,s±r=even,s=1,2,…,∞,r=1,2,…,∞,where T0 and T1 are the coefficient matrixes for A in the left and right hand of the above equilibrium equation, respectively; Q0 and Q1 are the coefficient matrixes for B in the left and right hand of the above equilibrium equation, respectively.If (17) and (23) are combined together, the buckling equation can be expressed as follows in the form of block matrix:(28)R0S0T0Q0AB=M~0R1S1T1Q1AB,where M~0 is the dimensionless critical moment, that is, the dimensionless buckling strength; R, S, T, and Q are the subblock coefficient matrixes constructed from analytical expressions included in (18)–(21) and (24)–(27), respectively.Equation (28) is the dimensionless analytical buckling equation of the LTB problem of a simply supported steel beam under linear distributed moment obtained in the present work.Obviously, from the mechanical point of view, the matrix at the left-hand side of the equation represents the linear stiffness matrix of the steel girder and is independent of the load, while the matrix on the left side of the equation depends on the load pattern, which in the present case can be considered as a geometric stiffness matrix due to the linear distribution of moments.It is noted that the dimensionless analytical solution of the buckling equation presented in the paper not only has a clear data structure and but also is easy to program and solve. In addition, unlike the FEM solutions published in literature, this solution can easily be verified by anyone.


## 3.4. Remarks on Definition of Dimensionless Critical Moment

There are different ways to define the dimensionless critical moment, in which the well-known definition is the one proposed by Kitipornchai et al. [25] and Kitipornchai and Wang [26, 27] and can be expressed as follows:(29)M~cr=McrLEIyGJk.Corresponding to this definition, the following dimensionless Wagner’s coefficient, the load position parameter, and the beam parameter are defined by Kitipornchai:(30)δ=βxLEIyGJk;ε=aLEIyGJk;K~=π2EIyh2GJkL2=πhLEIyGJk.Obviously, all these parameters lack clear physical and geometric meanings. More importantly, since there is a common factorEIy/GJk in (30), none of the parameters defined in (30) is independent; that is, they are not independent variables; thus it is difficult to use them to regress a rational design formula.Instead, the following dimensionless critical moment is used in this paper:(31)M~cr=Mcrπ2EIy/L2h=McrhPEy.This definition is first proposed by Professor Zhang [29] since 2008 and later used in the LTB analysis of various steel beams [30–35].It is the first finding in this paper that this definition can be derived naturally from the total potential energy based upon the new definition of displacement function, that is, (8); secondly, it has a clear physical meaning as shown in the last term in (31); that is, the product of Euler’s bending buckling load around the minor axis and the distance between the centroid of the top and bottom flange is used to normalize the actual critical moment.In addition, in the case of the LTB of a doubly symmetric beam under pure bending, the definition of (31) also has a definite geometric meaning; that is,(32)M~cr=McrhPEy=uθ=121+K-2.That is, the dimensionless critical moment is the dimensionless distance between the fixed axis of rotation and the shear center for a cross section. It is clear that, for the doubly symmetric beams, the fixed axis of rotation always lies beyond the tension flange due toK>0.Finally, and most importantly, the other dimensionless parameters defined in (11) also have definite physical meanings and are independent of each other and are therefore suitable for the formulation of a new design formula.


## 4. Approximate Analytical Solution

### 4.1. First-Order Approximation

Ifm=n=1, the following simple expression can be obtained from (28):(33)π42001+K2π4η2K21+η2A1B1=M~001+kπ441+kπ44-1+kπ4β~x2A1B1.Under the conditions that the coefficient determinant is zero, then the first-order approximation of the dimensionless critical moment of the I-beams under linear distributed moment can be written in a compact form:(34)M~cr=C1-C2a~+C3β~x+-C2a~+C3β~x2+η1+η21+K-2,where(35)C1=21+k;C2=0;C3=1or(36)M~cr=21+kβ~x+β~x2+η1+η21+K-2.Obviously, this 1st-order approximation only applies to cases wherek is greater than zero, and, at k=-1, the critical moment is infinite, which is unreasonable.Ifk=1 which means that the two end moments are equal, that is, the beam is under pure bending, then its dimensionless critical moment can be expressed in terms of (36) as follows:(37)M~cr=β~x+β~x2+η1+η21+K-2.It is easy to find that this expression is identical to the exact analytical solution in the textbook [1–6]. This conclusion indicates that the solution of the buckling equation presented in the paper is indeed an accurate solution in this special case.


### 4.2. Sixth-Order Approximation

Ifm=n=6, the following analytical expression can be obtained from (28):(38)R0S0T0Q06×6AB6×1=M~0R1S1T1Q16×6AB6×1,where(39)S0=T0=R1=0,R0=π420000008π400000081π42000000128π4000000625π42000000648π4,Q0=1+K2π4η2K21+η20000004+16K2π4η2K21+η20000009+81K2π4η2K21+η200000016+256K2π4η2K21+η200000025+625K2π4η2K21+η200000036+1296K2π4η2K21+η2,Ts,r1=Sr,s1=141+kπ489-1+kπ2016225-1+kπ2024-1+kπ21225329-1+kπ21+kπ49625-1+kπ20160441-1+kπ20021625-1+kπ2941+kπ443249-1+kπ2089-1+kπ2256225-1+kπ2076849-1+kπ241+kπ4128081-1+kπ2001000441-1+kπ20200081-1+kπ22541+kπ43000121-1+kπ2864-1+kπ212250329-1+kπ204320121-1+kπ291+kπ4,Q1=-121+kπ4β~x-409-1+kπ2β~x0-272225-1+kπ2β~x0-888-1+kπ2β~x1225-409-1+kπ2β~x-21+kπ4β~x-31225-1+kπ2β~x0-1160441-1+kπ2β~x00-31225-1+kπ2β~x-921+kπ4β~x-120049-1+kπ2β~x0-409-1+kπ2β~x-272225-1+kπ2β~x0-120049-1+kπ2β~x-81+kπ4β~x-328081-1+kπ2β~x00-1160441-1+kπ2β~x0-328081-1+kπ2β~x-2521+kπ4β~x-7320121-1+kπ2β~x-888-1+kπ2β~x12250-409-1+kπ2β~x0-7320121-1+kπ2β~x-181+kπ4β~x.Equation (38) gives the sixth-order approximation of the dimensionless critical moment of a steel beam under the linear distributed moment.It can be proved that the above result is exactly the same as that obtained by using trigonometric series with six terms as the modal function. This shows that the infinite series solution, that is, the exact analytical solution derived in previous section, is correct.Our research results [34, 35] have shown that, in most cases, the sixth-order approximation can achieve satisfactory results. However, if one wants to obtain more accurate results of the dimensionless critical moment, then more terms should be included in the numerical computations. Hence we shall present the numerical algorithm to solve the dimensionless analytical buckling equation, and the convergence performance is discussed in the next section.


## 4.1. First-Order Approximation

Ifm=n=1, the following simple expression can be obtained from (28):(33)π42001+K2π4η2K21+η2A1B1=M~001+kπ441+kπ44-1+kπ4β~x2A1B1.Under the conditions that the coefficient determinant is zero, then the first-order approximation of the dimensionless critical moment of the I-beams under linear distributed moment can be written in a compact form:(34)M~cr=C1-C2a~+C3β~x+-C2a~+C3β~x2+η1+η21+K-2,where(35)C1=21+k;C2=0;C3=1or(36)M~cr=21+kβ~x+β~x2+η1+η21+K-2.Obviously, this 1st-order approximation only applies to cases wherek is greater than zero, and, at k=-1, the critical moment is infinite, which is unreasonable.Ifk=1 which means that the two end moments are equal, that is, the beam is under pure bending, then its dimensionless critical moment can be expressed in terms of (36) as follows:(37)M~cr=β~x+β~x2+η1+η21+K-2.It is easy to find that this expression is identical to the exact analytical solution in the textbook [1–6]. This conclusion indicates that the solution of the buckling equation presented in the paper is indeed an accurate solution in this special case.


## 4.2. Sixth-Order Approximation

Ifm=n=6, the following analytical expression can be obtained from (28):(38)R0S0T0Q06×6AB6×1=M~0R1S1T1Q16×6AB6×1,where(39)S0=T0=R1=0,R0=π420000008π400000081π42000000128π4000000625π42000000648π4,Q0=1+K2π4η2K21+η20000004+16K2π4η2K21+η20000009+81K2π4η2K21+η200000016+256K2π4η2K21+η200000025+625K2π4η2K21+η200000036+1296K2π4η2K21+η2,Ts,r1=Sr,s1=141+kπ489-1+kπ2016225-1+kπ2024-1+kπ21225329-1+kπ21+kπ49625-1+kπ20160441-1+kπ20021625-1+kπ2941+kπ443249-1+kπ2089-1+kπ2256225-1+kπ2076849-1+kπ241+kπ4128081-1+kπ2001000441-1+kπ20200081-1+kπ22541+kπ43000121-1+kπ2864-1+kπ212250329-1+kπ204320121-1+kπ291+kπ4,Q1=-121+kπ4β~x-409-1+kπ2β~x0-272225-1+kπ2β~x0-888-1+kπ2β~x1225-409-1+kπ2β~x-21+kπ4β~x-31225-1+kπ2β~x0-1160441-1+kπ2β~x00-31225-1+kπ2β~x-921+kπ4β~x-120049-1+kπ2β~x0-409-1+kπ2β~x-272225-1+kπ2β~x0-120049-1+kπ2β~x-81+kπ4β~x-328081-1+kπ2β~x00-1160441-1+kπ2β~x0-328081-1+kπ2β~x-2521+kπ4β~x-7320121-1+kπ2β~x-888-1+kπ2β~x12250-409-1+kπ2β~x0-7320121-1+kπ2β~x-181+kπ4β~x.Equation (38) gives the sixth-order approximation of the dimensionless critical moment of a steel beam under the linear distributed moment.It can be proved that the above result is exactly the same as that obtained by using trigonometric series with six terms as the modal function. This shows that the infinite series solution, that is, the exact analytical solution derived in previous section, is correct.Our research results [34, 35] have shown that, in most cases, the sixth-order approximation can achieve satisfactory results. However, if one wants to obtain more accurate results of the dimensionless critical moment, then more terms should be included in the numerical computations. Hence we shall present the numerical algorithm to solve the dimensionless analytical buckling equation, and the convergence performance is discussed in the next section.


## 5. Numerical Solution and Convergence Performance

### 5.1. Matrix Expression of Eigenvalue Problem

In order to obtain a numerical solution to the dimensionless analytical solution, we shall take the number of trigonometric series in (8) and (9) as a finite number, such as m=n=N.If the common-used symbols used in traditional finite element method [8] are adopted, then (28) can be abbreviated as follows:(40)K02N×2NU2N=M~0KG2N×2NU2N,where(41)U2N=ABT,K02N×2N=R0T0S0Q0;KG2N×2N=R1T1S1Q1,where U2N is the buckling mode which is consisted of dimensionless undetermined coefficients (generalized coordinates); K02N×2N is the linear stiffness matrix of the steel beam; KG2N×2N is the geometric stiffness matrix of the steel beam under the linear distributed moment; M~0 is the dimensionless critical moment.Mathematically, the problem of solvingM~0 can ultimately be attributed to solving the generalized eigenvalue problem of (38), where M~0 is the smallest eigenvalue and U2N is the corresponding eigenvector (i.e., buckling mode).Obviously, the final form of the dimensionless buckling equation (see (38)) has not only a clear data structure that is easy to program but also a well-known form of the traditional finite element method and therefore has a definite physical meaning, while that given by Kitipornchai et al. [25–27] and McCann et al. [59] is a Hessian matrix, in which the dimensionless critical moment is implicit in that matrix. Therefore, their expressions lack clear physical meaning and are not easily solved by the existing mathematical software.


### 5.2. Convergence Performance and Numerical Algorithm

First, the convergence performance of the trigonometric series is presented. Typical graphs of the relationship between the dimensionless buckling moment and the number of terms of trigonometric series used in the modal trial functions are shown in Figure4. It can be seen that, in general, the trigonometric series is basically monotone convergence; the convergence rate for the negative end moment ratio is obviously lower than that for the positive end moment ratio.Figure 4
Convergence performance.
(a)
k
=
-
1
(b)
k
=
-
0.5
(c)
k
=
0.5
(d)
k
=
1Whenk changes from −1 to 1, the convergence rate is gradually accelerated, and the number of terms needed to obtain accurate numerical solutions is gradually reduced. Moreover, the results show that when k=-1, the convergence rate is the slowest, and more than six terms are needed to obtain accurate numerical solutions. However, when k=1, the convergence rate is the fastest, and only one term is needed to get the exact solution.Next, the numerical algorithm is discussed. In fact, any method used to solve the generalized eigenvalue problem is applicable to LTB problem of this paper, because it is different from the large-scale finite element analysis program, and the number of degrees of freedom involved in this paper is not large.It is noted that theeig(∗,∗) function provided in the MATLAB [60] can be used to obtain the eigenvalue and eigenvector of a given matrix. Even for matrixes with large dimensions, very high accuracy results can be obtained effectively. Therefore, this paper uses the eig(∗,∗) function to develop the developed MATLAB program.In order to automatically obtain a numerical solution of any desired accuracy, it is necessary to perform an iterative process in the developed MATLAB program in which the following convergence criteria are used:(42)M~0i+1-M~0iM~0i≤tol,where tol is the iteration tolerance. When the eigenvalue M~0 is required to 2s-digit accuracy [8], then tol=10-2s is used in the subsequent analysis.In order to facilitate the dimensionless parameter analysis and to perform the comparative study between the numerical solution from the exact dimensionless analytical solution and the ANSYS finite element analysis, a MATLAB program is developed according to the aforementioned numerical algorithm.


## 5.1. Matrix Expression of Eigenvalue Problem

In order to obtain a numerical solution to the dimensionless analytical solution, we shall take the number of trigonometric series in (8) and (9) as a finite number, such as m=n=N.If the common-used symbols used in traditional finite element method [8] are adopted, then (28) can be abbreviated as follows:(40)K02N×2NU2N=M~0KG2N×2NU2N,where(41)U2N=ABT,K02N×2N=R0T0S0Q0;KG2N×2N=R1T1S1Q1,where U2N is the buckling mode which is consisted of dimensionless undetermined coefficients (generalized coordinates); K02N×2N is the linear stiffness matrix of the steel beam; KG2N×2N is the geometric stiffness matrix of the steel beam under the linear distributed moment; M~0 is the dimensionless critical moment.Mathematically, the problem of solvingM~0 can ultimately be attributed to solving the generalized eigenvalue problem of (38), where M~0 is the smallest eigenvalue and U2N is the corresponding eigenvector (i.e., buckling mode).Obviously, the final form of the dimensionless buckling equation (see (38)) has not only a clear data structure that is easy to program but also a well-known form of the traditional finite element method and therefore has a definite physical meaning, while that given by Kitipornchai et al. [25–27] and McCann et al. [59] is a Hessian matrix, in which the dimensionless critical moment is implicit in that matrix. Therefore, their expressions lack clear physical meaning and are not easily solved by the existing mathematical software.


## 5.2. Convergence Performance and Numerical Algorithm

First, the convergence performance of the trigonometric series is presented. Typical graphs of the relationship between the dimensionless buckling moment and the number of terms of trigonometric series used in the modal trial functions are shown in Figure4. It can be seen that, in general, the trigonometric series is basically monotone convergence; the convergence rate for the negative end moment ratio is obviously lower than that for the positive end moment ratio.Figure 4
Convergence performance.
(a)
k
=
-
1
(b)
k
=
-
0.5
(c)
k
=
0.5
(d)
k
=
1Whenk changes from −1 to 1, the convergence rate is gradually accelerated, and the number of terms needed to obtain accurate numerical solutions is gradually reduced. Moreover, the results show that when k=-1, the convergence rate is the slowest, and more than six terms are needed to obtain accurate numerical solutions. However, when k=1, the convergence rate is the fastest, and only one term is needed to get the exact solution.Next, the numerical algorithm is discussed. In fact, any method used to solve the generalized eigenvalue problem is applicable to LTB problem of this paper, because it is different from the large-scale finite element analysis program, and the number of degrees of freedom involved in this paper is not large.It is noted that theeig(∗,∗) function provided in the MATLAB [60] can be used to obtain the eigenvalue and eigenvector of a given matrix. Even for matrixes with large dimensions, very high accuracy results can be obtained effectively. Therefore, this paper uses the eig(∗,∗) function to develop the developed MATLAB program.In order to automatically obtain a numerical solution of any desired accuracy, it is necessary to perform an iterative process in the developed MATLAB program in which the following convergence criteria are used:(42)M~0i+1-M~0iM~0i≤tol,where tol is the iteration tolerance. When the eigenvalue M~0 is required to 2s-digit accuracy [8], then tol=10-2s is used in the subsequent analysis.In order to facilitate the dimensionless parameter analysis and to perform the comparative study between the numerical solution from the exact dimensionless analytical solution and the ANSYS finite element analysis, a MATLAB program is developed according to the aforementioned numerical algorithm.


## 6. Verification by FEM Software

In this section, the general-purpose finite element software ANSYS is used to verify the correctness of the exact solution and 6th-order approximations given previously.


### 6.1. Description of FEM Model

It is well known that the BEAM189 in ANSYS software cannot correctly simulate the critical moment of LTB for the singly symmetric beams. Therefore, the SHELL63 element is used in this paper to simulate the buckling problem of the steel beams.The second reason is that since a typical steel I-beam is composed of three thin plates, it is theoretically correct using the SHELL elements to simulate the buckling problem of the I-beams. This is in accordance with Vlasov’s work [1] because his theory was established on the simplification of the shell theory and the concept of “reduced modulus” was used in his monograph (see equation (5.5) in [1]).It is shown that only in this way all kinds of buckling phenomena such as local buckling, distortion buckling, and lateral-torsional buckling (LTB) can be simulated and captured. However, only the LTB phenomena is concerned in this investigation. Consequently, when the SHELL63 element is used to simulate the desired lateral-torsional buckling of a steel beam, there is a problem; that is, for beams simulated by SHELL63 element, the Vlasov’s rigid section hypothesis cannot be automatically satisfied as the BEAM elements. Therefore, special measures have to be taken in the FEM simulation to ensure that the overall buckling; that is, LTB rather than the local buckling or the distortional buckling appears first.Some approaches have been proposed, such as the method of adding stiffeners or the method of treating the stiffeners as membrane elements [20] to approximate the rigid section hypothesis. However, it is found that, for k=-1, no matter how many stiffeners are added, no satisfactory results can be obtained. In order to overcome this difficulty, this paper presents an innovative simulation technique for the ANSYS software to satisfy the rigid section hypothesis.After exploration and modeling practice, we found that, in the ANSYS finite element model, the “CRIG” command is such a very effective command that you can be more confident to simulate the rigid section features of the buckled beams, without worrying about the failure of FEM numerical simulation. Figure5 shows the ANSYS finite element model after the CERIG command is used. Moreover, each end of the steel beam is simulated as the fork support as usual (Figure 6), and the end moment is simulated by applying a concentrated torque at the centroid of the cross section at each end of the steel beam (Figure 7).Figure 5
Simulation of rigid section hypothesis by “CERIG” command.Figure 6
Simulation of boundary conditions (fork supports).Figure 7
Simulation of end moments.The material properties are taken as follows: Young’s modulusE=206 GPa, Poison’s ratio υ=0.3, and the shear modulus of elasticity is given by G=E/21+υ.


### 6.2. Verification of FEM Model

Although the analytical and numerical solutions given in the previous section are dimensionless, the finite element analysis must be a dimensional analysis in which the section size and span (or slenderness) of the steel I-beams must be specified in advance.Considering that the finite element analysis here is of a verification nature, only three typical cross sections are selected for the study purpose. Section A is a doubly symmetric section, and both Section B and Section C are singly symmetric sections, in which Section B has a larger tension flange, whereas Section C has a larger compression flange. All pictures of the selected sections along with their geometrical properties are shown in Table1.Table 1
Section dimensions and geometric properties for the selected beams.
Section A
Section B
Section C
H
=
400 mm; h=372 mm;
H
=
400 mm; h=372 mm;
H
=
400 mm; h=372 mm;
h
w
=
344 mm; bf=400 mm;
h
w
=
344 mm; bf1=200 mm;
h
w
=
344 mm; bf1=400 mm;
t
f
=
28 mm; tw=18 mm;
b
f
2
=
400 mm; tf=28 mm;
b
f
2
=
200 mm; tf=28 mm;
I
y
=
2.987
×
10
-
4 m4;
t
w
=
18 mm; Iy=1.680×10-4 m4;
t
w
=
18 mm; Iy=1.680×10-4 m4;
J
k
=
6.523
×
10
-
6 m4;
J
k
=
5.059
×
10
-
6 m4;
J
k
=
5.059
×
10
-
6 m4;
I
ω
=
10.333
×
10
-
6 m6.
I
ω
=
2.296
×
10
-
6 m6.
I
ω
=
2.296
×
10
-
6 m6.The beam spans are chosen such that the span-to-depth ratios lie between 15 and 30. Table2 lists the dimensionless parameters associated with the calculation of the LTB critical moments.Table 2
Critical moments for the selected beams under pure bending.
Section
K
L/m
η
k
β
x
N = 6 McrTheory
N = 30  McrTheory
M
c
r
S
H
E
L
L
Diff.1/%
M
c
r
B
E
A
M
Diff.2/%
A
1.063
6
1
1
0
0.687
0.687
0.680
0.91
0.689
−0.35
0.797
8
1
1
0
0.802
0.802
0.801
0.22
0.808
−0.66
0.531
12
1
1
0
1.066
1.066
1.068
−0.23
1.076
−0.96
B
0.569
6
0.125
1
−0.322
0.391
0.391
0.383
1.85
0.641
61.35
0.427
8
0.125
1
−0.322
0.541
0.541
0.535
1.15
0.809
46.51
0.284
12
0.125
1
−0.322
0.871
0.871
0.867
0.52
1.162
30.43
C
0.569
6
8
1
0.322
1.035
1.035
1.048
−1.28
0.641
−39.11
0.427
8
8
1
0.322
1.185
1.185
1.204
−1.55
0.809
−33.11
0.284
12
8
1
0.322
1.515
1.515
1.539
−1.53
1.162
−25.02
Note. Diff.1 = (McrTheory-McrSHELL)/McrSHELL  ×  100%; Diff.2 = (McrTheory-McrBEAM)/McrBEAM  ×  100%.The above-mentioned FEM model is used to calculate the critical moment of the I-beams. Since, for a steel I-beam under pure bending, the critical moment has exact solution, that is, (37), then the FEM model of the beam with the sections listed in Table 1 is verified by this exact solution. Furthermore, the sixth-order approximate analytical solution and numerical solution of the critical moments are also calculated. All the results are listed in Table 2.It is observed that (1) for all cases where doubly symmetric I-beams are studied, the deviations from the FEM results are within 0.91%; (2) for all cases where singly symmetric I-beams are studied, the deviations from the FEM results are higher than those of the doubly symmetric beams, but the largest deviation is not more than 1.85%; (3) if three decimal places are retained, the results of the 6th-order approximate analytical solution and the numerical solution with 30 terms are exactly the same. This proves the appropriateness of the ANSYS FEM model developed in this paper.In addition, since no exact solution exits for the beams under linear distributed moment, the ANSYS FEM model is used to compare with the published data which is obtained from newly developed FEM model. Table3 listed the results of ABSYS and those of B3DW [15].Table 3
Critical moments for the beams of HEA-200 under linear distributed moment.
Number
k
SHELL/103 N⋅m
B3Dw/103 N⋅m
Diff.1/%
Theory/103 N⋅m
Diff.2/%
BEAM/103 N⋅m
Diff.3/%
L-1
1
82.064
82.69
0.76
81.872
−0.23
82.687
−0.99
L-2
0.75
93.574
94.79
1.30
93.358
−0.23
94.285
−0.98
L-3
0.5
108.090
108.94
0.79
107.853
−0.22
108.916
−0.98
L-4
0.25
126.419
127.45
0.82
126.175
−0.19
127.393
−0.96
L-5
0
149.142
150.45
0.88
148.935
−0.14
150.318
−0.92
L-6
−0.25
175.909
177.61
0.97
175.823
−0.05
177.353
−0.86
L-7
−0.5
204.151
206.32
1.06
204.317
0.08
205.944
−0.79
L-8
−0.75
225.915
228.38
1.09
226.436
0.23
228.096
−0.73
L-9
−1
219.686
221.97
1.04
220.378
0.31
221.905
−0.69
Section and its properties
H
=
190 mm; h=180 mm;
h
w
=
170 mm; bf=200 mm;
t
f
=
10 mm; tw=6.5 mm;
L
=
8 m; A=51.05 cm2;
I
y
=
3559.23 cm4; Iz=1333.33 cm4
J
k
=
14.8895 cm4; Iω=108×103 cm6
E
=
210 GPa
Note. Diff.1 = (B3Dw − SHELL)/SHELL × 100%; Diff.2 = (Theory − SHELL)/SHELL × 100%; Diff.3 = (Theory − BEAM)/BEAM × 100%.It can be observed that (1) all the results of B3DW are higher than those of the ANSYS model developed in this paper, and the largest deviation is less than 1.3%; (2) all the results given by the numerical solutions (referenced “theory”) are the lowest and closer to the those of ANSYS with the maximum difference within 0.31%. This provides an additional evidence of the appropriateness of the ANSYS FEM model, and the analytical and numerical solutions developed in the current work are verified preliminarily.It should be pointed out that the BEAM189 in ANSYS software can only be used to predict the critical moment for a doubly symmetric beam but will yield erroneous results for a singly symmetric beam.In summary, all the results show that the “CRIG” command in ANSYS can be used to simulate the Vlasov’s rigid section hypothesis, and its simulation effect is more effective and hence better than the method of adding stiffeners.


### 6.3. Verification of the Analytical and Numerical Solutions

This paragraph will further to verify the accuracy of the analytical and numerical solutions developed in this paper.For the steel I-beams with the section listed in Table1 under linear distributed moment, their critical moments and buckling modes are calculated by ANSY FEM model, sixth-order approximation, and numerical solution, respectively. The results of the calculated critical moments are listed in Table 4. The predicted buckling modes for the beams with doubly and singly symmetric sections are pictured in Figures 8 and 9, respectively.Table 4
Critical moments for the selected beams under linear distributed moment.
Section
K
L/m
η
k
β
x
N = 6 McrTheory
N = 30 McrTheory
M
c
r
F
E
M
Diff.1/%
A
1.063
6
1
0.5
0
0.906
0.906
0.897
0.99
1.063
6
1
0.1
0
1.179
1.179
1.164
1.30
1.063
6
1
0
0
1.265
1.265
1.247
1.45
1.063
6
1
−0.1
0
1.357
1.357
1.336
1.63
1.063
6
1
−0.5
0
1.766
1.766
1.720
2.69
1.063
6
1
−1
0
1.872
1.872
1.809
3.48
B
0.569
6
0.125
0.5
−0.322
0.512
0.512
0.502
1.92
0.569
6
0.125
0.1
−0.322
0.649
0.649
0.635
2.18
0.569
6
0.125
0
−0.322
0.689
0.689
0.673
2.29
0.569
6
0.125
−0.1
−0.322
0.730
0.730
0.713
2.41
0.569
6
0.125
−0.5
−0.322
0.908
0.908
0.883
2.88
0.569
6
0.125
−1
−0.322
1.075
1.075
1.088
1.19
C
0.569
6
8
0.5
0.322
1.365
1.365
1.385
−1.50
0.569
6
8
0.1
0.322
1.767
1.767
1.784
−0.93
0.569
6
8
0
0.322
1.889
1.889
1.900
−0.56
0.569
6
8
−0.1
0.322
2.017
2.017
2.029
−0.60
0.569
6
8
−0.5
0.322
2.263
2.262
2.228
1.55
0.569
6
8
−1
0.322
1.126
1.126
1.090
3.29
Note. Diff.1 = (McrTheory-McrFEM)/McrFEM  × 100%.Figure 8
Buckling mode graph of the beams with doubly symmetric section (Section A).
(a)
k
=
1
(b)
k
=
0.5
(c)
k
=
-
0.5
(d)
k
=
-
1Figure 9
Buckling mode graph of the beams with singly symmetric section (Section C).
(a)
k
=
1
(b)
k
=
0.5
(c)
k
=
-
0.5
(d)
k
=
-
1It can be found that (1) if three decimal places are retained, the results of the 6th-order approximate analytical solution and the numerical solution with 30 terms are exactly the same as shown in Table2; (2) when k changes from 1 to -1, the deviations from the FEM results gradually increase, and the largest deviation occurs in the case of k=-1 and is less than 3.54%; (3) there is close agreement between the buckling modes of the present theory and those of ANSYS. This validates both the analytical and numerical solutions developed in the current work, and it can be concluded that the dimensionless analytical solution presented provides accurate solutions for the LTB problem of the simply supported steel I-beams.It is noted that, for all I-beams studied in this paper, the results given by the 6th-order approximation are of high accuracy. Our other research results also support this conclusion [34–36]. Therefore, if you are able to accept the calculation results with three-decimal precision, the sixth-order approximation proposed in this paper is a simple and quick calculation method. However, in order to obtain a high-precision design formula, we propose to use the analytical solution with more terms; that is, use the numerical solution proposed in previous sections.Finally, it is noted that, for the case of Section A and Section B, SHELL solutions of ANSYS will provide lower critical moment predictions than the exact solution as expected, but, for Section C, the result is just the opposite. Therefore, for the case of the beam under linear distributed moment, FEM cannot be applied to formulate the design formula, because it lacks precision consistency in the analysis.


## 6.1. Description of FEM Model

It is well known that the BEAM189 in ANSYS software cannot correctly simulate the critical moment of LTB for the singly symmetric beams. Therefore, the SHELL63 element is used in this paper to simulate the buckling problem of the steel beams.The second reason is that since a typical steel I-beam is composed of three thin plates, it is theoretically correct using the SHELL elements to simulate the buckling problem of the I-beams. This is in accordance with Vlasov’s work [1] because his theory was established on the simplification of the shell theory and the concept of “reduced modulus” was used in his monograph (see equation (5.5) in [1]).It is shown that only in this way all kinds of buckling phenomena such as local buckling, distortion buckling, and lateral-torsional buckling (LTB) can be simulated and captured. However, only the LTB phenomena is concerned in this investigation. Consequently, when the SHELL63 element is used to simulate the desired lateral-torsional buckling of a steel beam, there is a problem; that is, for beams simulated by SHELL63 element, the Vlasov’s rigid section hypothesis cannot be automatically satisfied as the BEAM elements. Therefore, special measures have to be taken in the FEM simulation to ensure that the overall buckling; that is, LTB rather than the local buckling or the distortional buckling appears first.Some approaches have been proposed, such as the method of adding stiffeners or the method of treating the stiffeners as membrane elements [20] to approximate the rigid section hypothesis. However, it is found that, for k=-1, no matter how many stiffeners are added, no satisfactory results can be obtained. In order to overcome this difficulty, this paper presents an innovative simulation technique for the ANSYS software to satisfy the rigid section hypothesis.After exploration and modeling practice, we found that, in the ANSYS finite element model, the “CRIG” command is such a very effective command that you can be more confident to simulate the rigid section features of the buckled beams, without worrying about the failure of FEM numerical simulation. Figure5 shows the ANSYS finite element model after the CERIG command is used. Moreover, each end of the steel beam is simulated as the fork support as usual (Figure 6), and the end moment is simulated by applying a concentrated torque at the centroid of the cross section at each end of the steel beam (Figure 7).Figure 5
Simulation of rigid section hypothesis by “CERIG” command.Figure 6
Simulation of boundary conditions (fork supports).Figure 7
Simulation of end moments.The material properties are taken as follows: Young’s modulusE=206 GPa, Poison’s ratio υ=0.3, and the shear modulus of elasticity is given by G=E/21+υ.


## 6.2. Verification of FEM Model

Although the analytical and numerical solutions given in the previous section are dimensionless, the finite element analysis must be a dimensional analysis in which the section size and span (or slenderness) of the steel I-beams must be specified in advance.Considering that the finite element analysis here is of a verification nature, only three typical cross sections are selected for the study purpose. Section A is a doubly symmetric section, and both Section B and Section C are singly symmetric sections, in which Section B has a larger tension flange, whereas Section C has a larger compression flange. All pictures of the selected sections along with their geometrical properties are shown in Table1.Table 1
Section dimensions and geometric properties for the selected beams.
Section A
Section B
Section C
H
=
400 mm; h=372 mm;
H
=
400 mm; h=372 mm;
H
=
400 mm; h=372 mm;
h
w
=
344 mm; bf=400 mm;
h
w
=
344 mm; bf1=200 mm;
h
w
=
344 mm; bf1=400 mm;
t
f
=
28 mm; tw=18 mm;
b
f
2
=
400 mm; tf=28 mm;
b
f
2
=
200 mm; tf=28 mm;
I
y
=
2.987
×
10
-
4 m4;
t
w
=
18 mm; Iy=1.680×10-4 m4;
t
w
=
18 mm; Iy=1.680×10-4 m4;
J
k
=
6.523
×
10
-
6 m4;
J
k
=
5.059
×
10
-
6 m4;
J
k
=
5.059
×
10
-
6 m4;
I
ω
=
10.333
×
10
-
6 m6.
I
ω
=
2.296
×
10
-
6 m6.
I
ω
=
2.296
×
10
-
6 m6.The beam spans are chosen such that the span-to-depth ratios lie between 15 and 30. Table2 lists the dimensionless parameters associated with the calculation of the LTB critical moments.Table 2
Critical moments for the selected beams under pure bending.
Section
K
L/m
η
k
β
x
N = 6 McrTheory
N = 30  McrTheory
M
c
r
S
H
E
L
L
Diff.1/%
M
c
r
B
E
A
M
Diff.2/%
A
1.063
6
1
1
0
0.687
0.687
0.680
0.91
0.689
−0.35
0.797
8
1
1
0
0.802
0.802
0.801
0.22
0.808
−0.66
0.531
12
1
1
0
1.066
1.066
1.068
−0.23
1.076
−0.96
B
0.569
6
0.125
1
−0.322
0.391
0.391
0.383
1.85
0.641
61.35
0.427
8
0.125
1
−0.322
0.541
0.541
0.535
1.15
0.809
46.51
0.284
12
0.125
1
−0.322
0.871
0.871
0.867
0.52
1.162
30.43
C
0.569
6
8
1
0.322
1.035
1.035
1.048
−1.28
0.641
−39.11
0.427
8
8
1
0.322
1.185
1.185
1.204
−1.55
0.809
−33.11
0.284
12
8
1
0.322
1.515
1.515
1.539
−1.53
1.162
−25.02
Note. Diff.1 = (McrTheory-McrSHELL)/McrSHELL  ×  100%; Diff.2 = (McrTheory-McrBEAM)/McrBEAM  ×  100%.The above-mentioned FEM model is used to calculate the critical moment of the I-beams. Since, for a steel I-beam under pure bending, the critical moment has exact solution, that is, (37), then the FEM model of the beam with the sections listed in Table 1 is verified by this exact solution. Furthermore, the sixth-order approximate analytical solution and numerical solution of the critical moments are also calculated. All the results are listed in Table 2.It is observed that (1) for all cases where doubly symmetric I-beams are studied, the deviations from the FEM results are within 0.91%; (2) for all cases where singly symmetric I-beams are studied, the deviations from the FEM results are higher than those of the doubly symmetric beams, but the largest deviation is not more than 1.85%; (3) if three decimal places are retained, the results of the 6th-order approximate analytical solution and the numerical solution with 30 terms are exactly the same. This proves the appropriateness of the ANSYS FEM model developed in this paper.In addition, since no exact solution exits for the beams under linear distributed moment, the ANSYS FEM model is used to compare with the published data which is obtained from newly developed FEM model. Table3 listed the results of ABSYS and those of B3DW [15].Table 3
Critical moments for the beams of HEA-200 under linear distributed moment.
Number
k
SHELL/103 N⋅m
B3Dw/103 N⋅m
Diff.1/%
Theory/103 N⋅m
Diff.2/%
BEAM/103 N⋅m
Diff.3/%
L-1
1
82.064
82.69
0.76
81.872
−0.23
82.687
−0.99
L-2
0.75
93.574
94.79
1.30
93.358
−0.23
94.285
−0.98
L-3
0.5
108.090
108.94
0.79
107.853
−0.22
108.916
−0.98
L-4
0.25
126.419
127.45
0.82
126.175
−0.19
127.393
−0.96
L-5
0
149.142
150.45
0.88
148.935
−0.14
150.318
−0.92
L-6
−0.25
175.909
177.61
0.97
175.823
−0.05
177.353
−0.86
L-7
−0.5
204.151
206.32
1.06
204.317
0.08
205.944
−0.79
L-8
−0.75
225.915
228.38
1.09
226.436
0.23
228.096
−0.73
L-9
−1
219.686
221.97
1.04
220.378
0.31
221.905
−0.69
Section and its properties
H
=
190 mm; h=180 mm;
h
w
=
170 mm; bf=200 mm;
t
f
=
10 mm; tw=6.5 mm;
L
=
8 m; A=51.05 cm2;
I
y
=
3559.23 cm4; Iz=1333.33 cm4
J
k
=
14.8895 cm4; Iω=108×103 cm6
E
=
210 GPa
Note. Diff.1 = (B3Dw − SHELL)/SHELL × 100%; Diff.2 = (Theory − SHELL)/SHELL × 100%; Diff.3 = (Theory − BEAM)/BEAM × 100%.It can be observed that (1) all the results of B3DW are higher than those of the ANSYS model developed in this paper, and the largest deviation is less than 1.3%; (2) all the results given by the numerical solutions (referenced “theory”) are the lowest and closer to the those of ANSYS with the maximum difference within 0.31%. This provides an additional evidence of the appropriateness of the ANSYS FEM model, and the analytical and numerical solutions developed in the current work are verified preliminarily.It should be pointed out that the BEAM189 in ANSYS software can only be used to predict the critical moment for a doubly symmetric beam but will yield erroneous results for a singly symmetric beam.In summary, all the results show that the “CRIG” command in ANSYS can be used to simulate the Vlasov’s rigid section hypothesis, and its simulation effect is more effective and hence better than the method of adding stiffeners.


## 6.3. Verification of the Analytical and Numerical Solutions

This paragraph will further to verify the accuracy of the analytical and numerical solutions developed in this paper.For the steel I-beams with the section listed in Table1 under linear distributed moment, their critical moments and buckling modes are calculated by ANSY FEM model, sixth-order approximation, and numerical solution, respectively. The results of the calculated critical moments are listed in Table 4. The predicted buckling modes for the beams with doubly and singly symmetric sections are pictured in Figures 8 and 9, respectively.Table 4
Critical moments for the selected beams under linear distributed moment.
Section
K
L/m
η
k
β
x
N = 6 McrTheory
N = 30 McrTheory
M
c
r
F
E
M
Diff.1/%
A
1.063
6
1
0.5
0
0.906
0.906
0.897
0.99
1.063
6
1
0.1
0
1.179
1.179
1.164
1.30
1.063
6
1
0
0
1.265
1.265
1.247
1.45
1.063
6
1
−0.1
0
1.357
1.357
1.336
1.63
1.063
6
1
−0.5
0
1.766
1.766
1.720
2.69
1.063
6
1
−1
0
1.872
1.872
1.809
3.48
B
0.569
6
0.125
0.5
−0.322
0.512
0.512
0.502
1.92
0.569
6
0.125
0.1
−0.322
0.649
0.649
0.635
2.18
0.569
6
0.125
0
−0.322
0.689
0.689
0.673
2.29
0.569
6
0.125
−0.1
−0.322
0.730
0.730
0.713
2.41
0.569
6
0.125
−0.5
−0.322
0.908
0.908
0.883
2.88
0.569
6
0.125
−1
−0.322
1.075
1.075
1.088
1.19
C
0.569
6
8
0.5
0.322
1.365
1.365
1.385
−1.50
0.569
6
8
0.1
0.322
1.767
1.767
1.784
−0.93
0.569
6
8
0
0.322
1.889
1.889
1.900
−0.56
0.569
6
8
−0.1
0.322
2.017
2.017
2.029
−0.60
0.569
6
8
−0.5
0.322
2.263
2.262
2.228
1.55
0.569
6
8
−1
0.322
1.126
1.126
1.090
3.29
Note. Diff.1 = (McrTheory-McrFEM)/McrFEM  × 100%.Figure 8
Buckling mode graph of the beams with doubly symmetric section (Section A).
(a)
k
=
1
(b)
k
=
0.5
(c)
k
=
-
0.5
(d)
k
=
-
1Figure 9
Buckling mode graph of the beams with singly symmetric section (Section C).
(a)
k
=
1
(b)
k
=
0.5
(c)
k
=
-
0.5
(d)
k
=
-
1It can be found that (1) if three decimal places are retained, the results of the 6th-order approximate analytical solution and the numerical solution with 30 terms are exactly the same as shown in Table2; (2) when k changes from 1 to -1, the deviations from the FEM results gradually increase, and the largest deviation occurs in the case of k=-1 and is less than 3.54%; (3) there is close agreement between the buckling modes of the present theory and those of ANSYS. This validates both the analytical and numerical solutions developed in the current work, and it can be concluded that the dimensionless analytical solution presented provides accurate solutions for the LTB problem of the simply supported steel I-beams.It is noted that, for all I-beams studied in this paper, the results given by the 6th-order approximation are of high accuracy. Our other research results also support this conclusion [34–36]. Therefore, if you are able to accept the calculation results with three-decimal precision, the sixth-order approximation proposed in this paper is a simple and quick calculation method. However, in order to obtain a high-precision design formula, we propose to use the analytical solution with more terms; that is, use the numerical solution proposed in previous sections.Finally, it is noted that, for the case of Section A and Section B, SHELL solutions of ANSYS will provide lower critical moment predictions than the exact solution as expected, but, for Section C, the result is just the opposite. Therefore, for the case of the beam under linear distributed moment, FEM cannot be applied to formulate the design formula, because it lacks precision consistency in the analysis.


## 7. Proposed Dimensionless Design Formula

Nowadays, the concept of equivalent uniform moment coefficient (EUMF) is widely accepted in design codes in many countries. However, it is found that, even for the case of the doubly symmetric I-beams under linear distributed moment, the design formulas given by the codes/specifications of different countries are quite different. The corresponding comparison between the design formulas [61–65] with the exact solutions developed in current work is depicted in Figure 1. Due to space constraints, this paper does not intend to make too many evaluations of these results.This paper attempts to use a new mathematical model as the benchmark to solve this long-standing engineering problem with the above-mentioned dimensionless analytical solutions.


### 7.1. Trilinear Mathematical Model

Since the concept of equivalent uniform moment factor (EUMF) proposed by Salvadori [66] is widely recognized in the engineering community, similar to others, at first we have tried to use EUMF to describe the relationship between the critical moments and the gradient moment factork, but we ultimately failed.In order to analyze the cause of failure, through the parametric analysis, we first study its change law of the data. Figures10, 11, and 12 show the trend of the curves of the critical moments versus the gradient moment factor k with varied η and K, respectively. It can be seen that the trend of all the curves are very different that no rules can be followed. This means that it is not possible to describe such a complex change with a single EUMF. In other words, a single EUMF does not exist in the case of I-beams with unequal flanges subjected to unequal end moments; hence, we must break through the shackles of traditional concepts, so that we can get a breakthrough.Figure 10
Distribution pattern of critical moments with variedK and k.Figure 11
Critical moments versus gradient moment factor with variedη (K=2).Figure 12
Critical moments versus gradient moment factor with variedK (η=0.5).With the aid of the idea of using a piecewise linear approximation to simulate an actual curve, this paper proposes a new trilinear mathematical model used as the benchmark of formulating the design formula. This new idea is represented by the graphics shown in Figure13. There are four control points in this new model; that is, Mcr-1, Mcr-0.5, Mcr0, Mcr+1 are the four control moments. Since Mcr+1 can be expressed by the existing exact solution, that is, (37), only the other three moments are undetermined parameters, which should be obtained by the nonlinear regression techniques.Figure 13
Proposed trilinear model as the benchmark of design formula.


### 7.2. Parameter Regression Techniques and Comparison

As we all know, there are three key issues in regression technology, namely, the choice of target formula, how to generate data, and how to regress the unknown coefficient in the target formula.Firstly, according to the results of the previous parametric analysis, after a number of analyses, observations, and tries, we finally selected the following target formula:(43)M~cr=C1C3β~x+C3β~x2+C4ηC5+C6η21+K-2,where C1, C3, C4, C5, and C6 are undetermined regression coefficients and the later three are newly defined in this paper.In addition, the study also found that these undetermined regression coefficients for the cross sections with lager top flange are completely different from those of other cross sections in the range (-1.0≤k<-0.5), so we cannot expect to describe all the cases with a set of unified regression coefficients by using the target formula (43). Therefore, two sets of regression coefficients have to be used for C1, C3, C4, C5, and C6 in the proposed design formula for the case where k=-1.0 and k=-0.5.Next, the MATLAB program formulated in the paper can be employed to generate the above-mentioned dimensionless data. In the following data analysis, the end moment ratiok is fixed for each control moment, while β~x, η, β~xK are varied. K is changed from 0.2 to 5.0 with the varied step size of 0.02~0.05; η is changed from 0.1 to 10.0 with the varied step size of 0.05~1.0; β~x is changed from −0.4 to +0.4 with step size of 0.01.Thirdly, due to the highly nonlinear nature of the target formula, we chose the 1stOpt [67] software to do this regression [30–35]. 1stOpt is a software package for the analysis of mathematical optimization problem. The nonlinear curve fit of 1stOpt is powerful compared to any other similar software package available today. The greatest feature is that it no longer needs the end-user to provide or guess the initial start-values for each parameter, but, randomly generated by machine, the probability for finding the correct solution is higher than the others.Taking a cross section having a larger top flange and a negative moment (k=-1) as an example, based on the 475 sets of dimensionless data (β~x, η, K, M~cr) obtained from the dimensionless analytical solution, using Levenberg-Marquardt algorithm in the 1stOpt software, then after 39 iterations, the resulting regression coefficients are obtained as shown in Table 5 and the corresponding correlation coefficient is 0.99545.Table 5
Parameter
Best estimate
c
1
2.220165826
c
3
−0.458694891
c
4
7.259270297
c
5
1.482366765
c
6
2.735838541The resulting relationship between the target value and the calculated value is depicted in Figure14. In this figure, the red line represents the calculated curve, and the blue line represents the target curve. It is found that quite satisfactory regression results have been obtained.Figure 14
Relationship of the target value and calculated value.Using a similar approach, we can get all the undetermined coefficients and the relevant results are listed in Table6.Table 6
Formula and parameters for I-beams under linear distributed moment.
k
=
-
1
k
=
-
0.5
k
=
-
0
k
=
+
1
C
1
0.234 (2.220)
1.387 (1.095)
5.983
1
C
3
4.674 (−0.459)
0.822 (0.922)
0.287
1
C
4
72.762 (7.259)
0.615 (10.860)
11.468
1
C
5
0.861 (1.482)
0.587 (1.629)
11.713
1
C
6
0.566 (2.736)
0.273 (1.154)
10.446
1
M
~
0
=
M
1
π
2
E
I
y
/
L
2
h;
β
~
x
=
β
x
h;
η
=
I
1
I
2;
K
=
π
2
E
I
ω
G
J
k
L
2
M
~
c
r
=
C
1
C
3
β
~
x
+
C
3
β
~
x
2
+
C
4
η
C
5
+
C
6
η
2
1
+
K
-
2
Note. (1) Parameters in parentheses for I-beams of singly symmetric sections with larger top flange only; (2) applicable scope of the formula: K = 0.2~5.0; η = 0.1~10.0; β~x = −0.4~+0.4.It must be pointed out that the whole regression process takes only a few minutes, of which the calculation process of nearly 500 sets of data only takes a few seconds. If you intend to use the FEM to obtain such a large number of data sets, it may take at least a month. This fully demonstrates that the dimensionless solution has an unparalleled advantage over the finite element method in terms of efficiency and quality in obtaining a large number of regression data.Finally, the proposed dimensionless design formula is used to calculate the critical moments of the I-beams with the sections listed in Table1 and the results are compared with the theoretical solution. All results are plotted in Figure 15 and listed in Table 7.Table 7
Comparison of the results of design formula with those of theory for I-beams under linear distributed moment.
Section
K
L/m
η
k
β
x
N
=
6 McrTheory
N
=
30 McrTheory
M
c
r
R
e
g
r
e
s
s
i
o
n
Diff./%
A
1.063
6
1
1
0
0.687
0.687
0.687
0.00
1.063
6
1
0
0
1.265
1.265
1.255
−0.72
1.063
6
1
−0.5
0
1.766
1.766
1.737
−1.68
1.063
6
1
−1
0
1.872
1.872
1.919
2.54
B
0.569
6
0.125
1
−0.322
0.391
0.390
0.391
0.00
0.569
6
0.125
0
−0.322
0.689
0.689
0.664
0.51
0.569
6
0.125
−0.5
−0.322
0.908
0.908
0.915
5.43
0.569
6
0.125
−1
−0.322
1.075
1.075
1.194
11.08
C
0.569
6
8
1
0.322
1.035
1.035
1.035
0.00
0.569
6
8
0
0.322
1.961
1.961
1.958
−0.19
0.569
6
8
−0.5
0.322
2.263
2.262
2.286
0.88
0.569
6
8
−1
0.322
1.126
1.126
1.150
7.00Figure 15
Comparison between theory and the proposed formula.The results show that the proposed formula is of high accuracy and is applicable to both the doubly symmetric section and the singly symmetric section under positive and negative gradient moment.It must be noted that the trilinear mathematical model can also be converted to a cubic curve model or by adding more points to obtain a quadratic or higher-order curve model. Due to the limited space, these will be discussed elsewhere.


## 7.1. Trilinear Mathematical Model

Since the concept of equivalent uniform moment factor (EUMF) proposed by Salvadori [66] is widely recognized in the engineering community, similar to others, at first we have tried to use EUMF to describe the relationship between the critical moments and the gradient moment factork, but we ultimately failed.In order to analyze the cause of failure, through the parametric analysis, we first study its change law of the data. Figures10, 11, and 12 show the trend of the curves of the critical moments versus the gradient moment factor k with varied η and K, respectively. It can be seen that the trend of all the curves are very different that no rules can be followed. This means that it is not possible to describe such a complex change with a single EUMF. In other words, a single EUMF does not exist in the case of I-beams with unequal flanges subjected to unequal end moments; hence, we must break through the shackles of traditional concepts, so that we can get a breakthrough.Figure 10
Distribution pattern of critical moments with variedK and k.Figure 11
Critical moments versus gradient moment factor with variedη (K=2).Figure 12
Critical moments versus gradient moment factor with variedK (η=0.5).With the aid of the idea of using a piecewise linear approximation to simulate an actual curve, this paper proposes a new trilinear mathematical model used as the benchmark of formulating the design formula. This new idea is represented by the graphics shown in Figure13. There are four control points in this new model; that is, Mcr-1, Mcr-0.5, Mcr0, Mcr+1 are the four control moments. Since Mcr+1 can be expressed by the existing exact solution, that is, (37), only the other three moments are undetermined parameters, which should be obtained by the nonlinear regression techniques.Figure 13
Proposed trilinear model as the benchmark of design formula.


## 7.2. Parameter Regression Techniques and Comparison

As we all know, there are three key issues in regression technology, namely, the choice of target formula, how to generate data, and how to regress the unknown coefficient in the target formula.Firstly, according to the results of the previous parametric analysis, after a number of analyses, observations, and tries, we finally selected the following target formula:(43)M~cr=C1C3β~x+C3β~x2+C4ηC5+C6η21+K-2,where C1, C3, C4, C5, and C6 are undetermined regression coefficients and the later three are newly defined in this paper.In addition, the study also found that these undetermined regression coefficients for the cross sections with lager top flange are completely different from those of other cross sections in the range (-1.0≤k<-0.5), so we cannot expect to describe all the cases with a set of unified regression coefficients by using the target formula (43). Therefore, two sets of regression coefficients have to be used for C1, C3, C4, C5, and C6 in the proposed design formula for the case where k=-1.0 and k=-0.5.Next, the MATLAB program formulated in the paper can be employed to generate the above-mentioned dimensionless data. In the following data analysis, the end moment ratiok is fixed for each control moment, while β~x, η, β~xK are varied. K is changed from 0.2 to 5.0 with the varied step size of 0.02~0.05; η is changed from 0.1 to 10.0 with the varied step size of 0.05~1.0; β~x is changed from −0.4 to +0.4 with step size of 0.01.Thirdly, due to the highly nonlinear nature of the target formula, we chose the 1stOpt [67] software to do this regression [30–35]. 1stOpt is a software package for the analysis of mathematical optimization problem. The nonlinear curve fit of 1stOpt is powerful compared to any other similar software package available today. The greatest feature is that it no longer needs the end-user to provide or guess the initial start-values for each parameter, but, randomly generated by machine, the probability for finding the correct solution is higher than the others.Taking a cross section having a larger top flange and a negative moment (k=-1) as an example, based on the 475 sets of dimensionless data (β~x, η, K, M~cr) obtained from the dimensionless analytical solution, using Levenberg-Marquardt algorithm in the 1stOpt software, then after 39 iterations, the resulting regression coefficients are obtained as shown in Table 5 and the corresponding correlation coefficient is 0.99545.Table 5
Parameter
Best estimate
c
1
2.220165826
c
3
−0.458694891
c
4
7.259270297
c
5
1.482366765
c
6
2.735838541The resulting relationship between the target value and the calculated value is depicted in Figure14. In this figure, the red line represents the calculated curve, and the blue line represents the target curve. It is found that quite satisfactory regression results have been obtained.Figure 14
Relationship of the target value and calculated value.Using a similar approach, we can get all the undetermined coefficients and the relevant results are listed in Table6.Table 6
Formula and parameters for I-beams under linear distributed moment.
k
=
-
1
k
=
-
0.5
k
=
-
0
k
=
+
1
C
1
0.234 (2.220)
1.387 (1.095)
5.983
1
C
3
4.674 (−0.459)
0.822 (0.922)
0.287
1
C
4
72.762 (7.259)
0.615 (10.860)
11.468
1
C
5
0.861 (1.482)
0.587 (1.629)
11.713
1
C
6
0.566 (2.736)
0.273 (1.154)
10.446
1
M
~
0
=
M
1
π
2
E
I
y
/
L
2
h;
β
~
x
=
β
x
h;
η
=
I
1
I
2;
K
=
π
2
E
I
ω
G
J
k
L
2
M
~
c
r
=
C
1
C
3
β
~
x
+
C
3
β
~
x
2
+
C
4
η
C
5
+
C
6
η
2
1
+
K
-
2
Note. (1) Parameters in parentheses for I-beams of singly symmetric sections with larger top flange only; (2) applicable scope of the formula: K = 0.2~5.0; η = 0.1~10.0; β~x = −0.4~+0.4.It must be pointed out that the whole regression process takes only a few minutes, of which the calculation process of nearly 500 sets of data only takes a few seconds. If you intend to use the FEM to obtain such a large number of data sets, it may take at least a month. This fully demonstrates that the dimensionless solution has an unparalleled advantage over the finite element method in terms of efficiency and quality in obtaining a large number of regression data.Finally, the proposed dimensionless design formula is used to calculate the critical moments of the I-beams with the sections listed in Table1 and the results are compared with the theoretical solution. All results are plotted in Figure 15 and listed in Table 7.Table 7
Comparison of the results of design formula with those of theory for I-beams under linear distributed moment.
Section
K
L/m
η
k
β
x
N
=
6 McrTheory
N
=
30 McrTheory
M
c
r
R
e
g
r
e
s
s
i
o
n
Diff./%
A
1.063
6
1
1
0
0.687
0.687
0.687
0.00
1.063
6
1
0
0
1.265
1.265
1.255
−0.72
1.063
6
1
−0.5
0
1.766
1.766
1.737
−1.68
1.063
6
1
−1
0
1.872
1.872
1.919
2.54
B
0.569
6
0.125
1
−0.322
0.391
0.390
0.391
0.00
0.569
6
0.125
0
−0.322
0.689
0.689
0.664
0.51
0.569
6
0.125
−0.5
−0.322
0.908
0.908
0.915
5.43
0.569
6
0.125
−1
−0.322
1.075
1.075
1.194
11.08
C
0.569
6
8
1
0.322
1.035
1.035
1.035
0.00
0.569
6
8
0
0.322
1.961
1.961
1.958
−0.19
0.569
6
8
−0.5
0.322
2.263
2.262
2.286
0.88
0.569
6
8
−1
0.322
1.126
1.126
1.150
7.00Figure 15
Comparison between theory and the proposed formula.The results show that the proposed formula is of high accuracy and is applicable to both the doubly symmetric section and the singly symmetric section under positive and negative gradient moment.It must be noted that the trilinear mathematical model can also be converted to a cubic curve model or by adding more points to obtain a quadratic or higher-order curve model. Due to the limited space, these will be discussed elsewhere.


## 8. Conclusions

Because there is no reliable analytical solution available, this study is initially intended to develop an exact analytical solution to study the applicability of the design formulas published in literature. However, we found that results given by the existing design formulas vary widely, especially in the range-1.0≤k<0, so, according to the idea that scientific research should serve for engineering, a new dimensionless design formula is proposed for practical engineering and the following conclusions can be drawn.(1) There is a challenge to part of the theory of monosymmetric beams, called the Wagner hypothesis, was presented by Ojalvo (1981). This challenge can be fully resolved by the Plate-Beam Theory put forward by the author, in which only the commonly used plate and beam theory was used and the Vlasov’s warping function was discarded.(2) For the case of Section A and Section B, SHELL solutions of ANSYS will provide lower critical moment predictions than the exact solution as expected, but, for Section C, the result is just the opposite. Therefore, for the case of the beam under linear distributed moment, FEM cannot be applied to formulate the design formula, because it lacks precision consistency in the analysis.(3) By introducing a new dimensionless coefficient of lateral deflection, new dimensionless critical moment and dimensionless Wagner’s coefficient are derived naturally from the total potential energy. These new dimensionless parameters are independent of each other and thus can be used as independent variables in the process of regression formula.(4) The results show that, in the case of the simple-supported I-beams under linearly distributed moments, the dimensionless analytical solution presented in current work is exact in nature and has an unparalleled advantage over the finite element method in terms of efficiency and quality in obtaining a large number of regression data. Therefore, it provides a more convenient and effective method for the formulation of a more accurate design formula than finite element method.(5) It is found that, for all I-beams studied in this paper, the results given by the 6th-order approximation are of high accuracy. Our other research results also support this conclusion [30–35]. Therefore, if you can accept the calculated results of three-decimal precision, then the sixth-order approximation presented in this paper is a simple and quick calculation method.(6) The dimensionless parameters defined by Kitipornchai lack clear physical and geometric meanings. More importantly, none of these parameters is independent; thus it is difficult to use them to regress a rational design formula.(7) The dimensionless analytical buckling equation is expressed as the generalized eigenvalue problem as that used in the traditional finite element method. This not only facilitates the solution in MATLAB programming but also makes the analytical solution have a clear physical meaning, more easy to understand and apply.(8) To satisfy Vlasov’s rigid section hypothesis, the “CRIG” command is used in the FEM model of ANSYS. The results show that its simulation effect is more effective and hence better than the method of adding stiffeners.(9) BEAM189 in ANSYS software can only be used to predict the critical moment for a doubly symmetric beam but will yield erroneous results for a singly symmetric beam.(10) It is found that the concept of equivalent uniform moment factor (EUMF) is not applicable for the simple-supported I-beams with unequal flanges subjected to unequal end moments; that is, a single EUMF does not exist in this case.(11) This paper proposes a new trilinear mathematical model as the benchmark of formulating the design formula. The results show that the proposed design formula has a high accuracy and is applicable to both the doubly symmetric section and the singly symmetric section under positive and negative gradient moment. This model provides a new idea for the formulation of more reasonable design formulas, so it is suggested that the designers of the design specifications can formulate the corresponding design formulas with reference to this model.



---

*Source: 4838613-2017-12-31.xml*