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.
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.
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 Section 2, 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 (), thin-walled () 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 , the following transverse displacements were proposed by Vlasov:where , are the transverse displacements of the shear center, is the twist angle of the cross section, and , are the coordinates of the shear center.
Invoking assumption , that is, by setting the shear strain of the middle surface to zero, the following longitudinal displacement was deduced: 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, where is the lateral deflection of the shear center, is the twist angle of the cross section, is the minor axis flexural rigidity, is the St. Venant torsion rigidity, is the Vlasov warping torsion rigidity, is the moment function about the major axis, is the length of the beam, is the distributed load acting in the vertical plane of the I-beams, is the load position parameter, which is obtained by subtracting the -coordinate of the load acting point from that of the shear center, and is the Wagner coefficient, and it is defined aswhere , is the coordinates of an arbitrary point with respect to the origin (the centroid of the cross section) and 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]):
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:
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 and 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 Figure 2. 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.
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 (Figure 3). 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].
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 Figure 2, it is easy to obtain the moment function as follows:where is the end moment ratio and , are the left and right end moment, respectively, and .
For the case of the simply supported I-beams, the following modal trial functions are proposed for the lateral deflection and twist angle [29–35], which can be expressed as the sum of trigonometric series as follows:where and are the undetermined dimensionless coefficients for the lateral deflection and twist angle , respectively; is the distance between the centroid of top and bottom flange.
From the perspective of structural dynamics, and 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 so that and are both undetermined dimensionless coefficients, we introduce into the expression of 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,
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 , one can get the integration results of the total potential energy.
Without loss of generality, the following new dimensionless parameters [29–35] are introduced:where is the dimensionless end moment; is the distance between the centroid of the top and bottom flange; 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; , are the 2nd moment of inertia about minor axes of the top and bottom flange, respectively; is the dimensionless Wagner’s coefficient.
In the following derivation, the energy equation of (3) is first multiplied by a common factor , 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:
The total potential energy of this problem is the sum of the above equations; that is,
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:Note that, in the above equations, and are arbitrary and independent.
From the condition , we can get the equilibrium equation with respect to the lateral deflection as follows:
For each fixed (corresponding to Line ), herein, one homogeneous equation corresponding to , , and 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:we can rewrite (15) in the following matrix form:where and are the column vectors constituted by all undetermined coefficients and , respectively; and are the coefficient matrixes for in the left and right hand of the above equilibrium equation, respectively; and are the coefficient matrixes for in the left and right hand of the above equilibrium equation, respectively.
Similarly, from the condition , we can get the equilibrium equation with respect to the twist rotation as follows:
With the notations in (16), the above equation can also be arranged in the following matrix form:where and are the coefficient matrixes for in the left and right hand of the above equilibrium equation, respectively; and are the coefficient matrixes for 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:where is the dimensionless critical moment, that is, the dimensionless buckling strength; , , , and 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:
Corresponding to this definition, the following dimensionless Wagner’s coefficient, the load position parameter, and the beam parameter are defined by Kitipornchai:
Obviously, all these parameters lack clear physical and geometric meanings. More importantly, since there is a common factor 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:
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,
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 to .
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
If , the following simple expression can be obtained from (28):
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:whereor
Obviously, this 1st-order approximation only applies to cases where is greater than zero, and, at , the critical moment is infinite, which is unreasonable.
If 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:
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
If , the following analytical expression can be obtained from (28):where
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 .
If the common-used symbols used in traditional finite element method [8] are adopted, then (28) can be abbreviated as follows:wherewhere is the buckling mode which is consisted of dimensionless undetermined coefficients (generalized coordinates); is the linear stiffness matrix of the steel beam; is the geometric stiffness matrix of the steel beam under the linear distributed moment; is the dimensionless critical moment.
Mathematically, the problem of solving can ultimately be attributed to solving the generalized eigenvalue problem of (38), where is the smallest eigenvalue and 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 Figure 4. 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.
(a)
(b)
(c)
(d)
When 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 , the convergence rate is the slowest, and more than six terms are needed to obtain accurate numerical solutions. However, when , 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 the 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 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:where tol is the iteration tolerance. When the eigenvalue is required to 2-digit accuracy [8], then 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 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 , 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. Figure 5 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).
The material properties are taken as follows: Young’s modulus GPa, Poison’s ratio , and the shear modulus of elasticity is given by .
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 Table 1.
The beam spans are chosen such that the span-to-depth ratios lie between 15 and 30. Table 2 lists the dimensionless parameters associated with the calculation of the LTB critical moments.
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. Table 3 listed the results of ABSYS and those of B3DW [15].
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 Table 1 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.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
It 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 Table 2; (2) when changes from 1 to , the deviations from the FEM results gradually increase, and the largest deviation occurs in the case of 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 factor k, 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. Figures 10, 11, and 12 show the trend of the curves of the critical moments versus the gradient moment factor with varied and , 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.
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 Figure 13. There are four control points in this new model; that is, , , , are the four control moments. Since 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.
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:where , , , , and 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 (), 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 , , , , and in the proposed design formula for the case where and .
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 ratio is fixed for each control moment, while , , are varied. 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; 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 () as an example, based on the 475 sets of dimensionless data (, , , ) 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.
The resulting relationship between the target value and the calculated value is depicted in Figure 14. 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.
Using a similar approach, we can get all the undetermined coefficients and the relevant results are listed in Table 6.
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 Table 1 and the results are compared with the theoretical solution. All results are plotted in Figure 15 and listed in Table 7.
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 , 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.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work is supported by the National Natural Science Foundation of China (51578120 and 51178087) and the Scientific Research Foundation of Nanjing Institute of Technology (YKJ201617).