The effect of cross-section geometry on the lateral-torsional behavior of thin-walled beams: Analytical and numerical studies

Thin-Walled Structures 184 (2023) 110535 Contents lists available at ScienceDirect Thin-Walled Structures journal homepage: www.elsevier.com/locate/tws Full length article The effect of cross-section geometry on the lateral–torsional behavior of thin-walled beams: Analytical and numerical studies Muhammad Z. Haffara, Martin Horáčekb, Sándor Ádánya,∗ aDepartment of Structural Mechanics, Faculty of Civil Engineering, Budapest University of Technology and Economics, H-1111 Budapest, Műegyetem rkp. 3, Hungary bInstitute of Metal and Timber Structures, Faculty of Civil Engineering, Brno University of Technology, Veveří 331/95, 602 00 Brno, Czech Republic ARTICLE INFO Keywords: Lateral–torsional buckling Geometrically nonlinear analysis Geometric imperfections Mono-symmetric cross-sections ABSTRACT In this paper, the elastic lateral–torsional behavior of simple beams is discussed by presenting a novel analytical solution and performing numerical studies. The motivation of the presented research is the observation that classic analytical prediction and finite element prediction are, typically, significantly different when the secondorder nonlinear behavior of beams with initial imperfections is analyzed. To understand and explain the observed differences, a novel analytical model is worked out for the geometrically nonlinear analysis of beams with initial geometric imperfections. The advancement in the presented analytical solution is the explicit consideration of the changing geometry as the load increases. The most important steps of the derivations are summarized, and the resulting formulae are briefly discussed. The derivations are done for general crosssections, however, the bending is assumed to act in one of the principal planes. Numerical studies are also presented, focusing on mono-symmetric cross-sections. As part of the numerical studies, first, the results of the new analytical formulae are compared to those from shell finite element analysis. The results suggest that the new formulae can capture the most essential elements of the behavior observed in the shell finite element calculations, justifying that the cross-section shape might have a significant effect on the nonlinear lateral–torsional behavior of beams. Then the effect of the lateral–torsional buckling is predicted by calculating buckling reduction factors, using the results of the geometrically nonlinear finite element calculations. The capacity prediction results, again, justify that the cross-section shape, as well as the sign of the assumed geometric imperfection, might have a non-negligible effect on the buckling reduction factors, and on the capacity of the member. 1. Introduction If a structure or structural part is slender, a potential failure mode is buckling. Buckling may take place in various forms, but the phenomenon can usually be described as follows: when the load intensity is gradually increased on the structure, the displacements are slowly increasing (i.e. primary displacements), but as the load approximates a certain level, the structure starts to develop rapidly increasing displacements (i.e., secondary displacements) so that the nature and/or the direction of the secondary displacements is distinctly different from that of the primary ones. If the structure is a beam and it is subjected to uniaxial bending in one of the principal planes, the primary displacements are in the plane of the loading, while the secondary displacements involve lateral translation (i.e., translations perpendicular to the plane of loading) and twisting rotation. In the case of beams, this phenomenon is identified as lateral–torsional buckling (LTB). There are multiple ways to analyze the buckling phenomenon. If the structure is free from imperfections and its material is perfectly ∗Corresponding author. E-mail address: [email protected] (S. Ádány). elastic, the analysis is usually termed linear buckling analysis (LBA), which leads to the buckling shapes and critical load values (i.e., critical moments in the case of LTB). Closed-form analytical solutions for the critical moments are known and can be found in textbooks [1–3], at least for simpler cases. However, when the cross-section, loading, or boundary conditions of the beam are more complex, it is not easy to find analytical solutions; that is why the topic is under research till now [4–7]. Another way of buckling analysis is when some geometric imperfections are directly considered. If the analysis is elastic, it is popularly abbreviated as GNI analysis or GNIA (i.e., Geometrically Non-linear Analysis with Imperfections). The approach was first applied by Young for columns [8]; it was shown that the initial imperfection is amplified due to the compressive force. The formula for the amplification factor is known as the Young formula and is widely used even in design calculations. It can be shown that the formula is valid for any elastic structure with initial geometric imperfection, at least under certain conditions [9]. GNIA can also be used for capacity prediction, as first https://doi.org/10.1016/j.tws.2023.110535 Received 17 November 2022; Received in revised form 5 January 2023; Accepted 5 January 2023 Available online 17 January 2023 0263-8231/©2023 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). M.Z. Haffar, M. Horáček and S. Ádány Thin-Walled Structures 184 (2023) 110535 proposed by Ayrton and Perry for compressed columns [10]. This approach is the basis of the European buckling curves, proposed by [11], which are still used in the current Eurocode standards [12,13], both for column buckling and LTB. Further generalization of the approach is proposed in [14]. A more advanced version of buckling analysis considers not only initial imperfections but also material nonlinearity, known as GMNI analysis or simply GMNIA. While GNIA can be performed analytically, at least for simpler cases, GMNIA requires some numerical tool, which is usually the finite element method. Using GMNIA the load– displacement path can be directly established, and the maximum point of the load–displacement curve identifies the capacity of the structure. Both GNIA and GMNIA need some initial geometric imperfection. In the practice, it is either assumed in a predefined shape, or as a buckling shape from an LBA. The advantage of using a buckling shape as initial imperfection is that this approach is quite general and automatic, i.e., can be applied to virtually any structure and during the process, it requires no or little engineering decision. The proper determination of the imperfection magnitude is a crucial question that is not fully answered yet. (In many structural problems it is also a question which buckling shape is to be used, but in the case of simple beams the first buckling mode is a proper choice.) The classic analytical solution for GNIA is the Young formula, recently it was shown in [15,16], that there is a discrepancy between the results predicted by the classic analytical GNI solution (see e.g., [9,17]) and those calculated by shell FEM GNI analyses. The discrepancies can be important. Also, the discrepancies are not limited to the difference of certain numerical values, some basic features of the behavior are affected. To reveal the reasons of the experienced discrepancies, the authors developed an advanced analytical model to solve the GNIA problem for beams subjected to lateral–torsional behavior. A simpler version of the analytical model has been introduced and discussed in [18], where doubly symmetric I-section beams were considered. In this present paper a generalized version of the analytical model is presented. The aim is, still, to have an analytical solution for the nonlinear behavior of simple beams with initial geometric imperfections, however, here the cross-section is arbitrary. Since non-symmetric cross-sections are rarely used in the practice, the focus is on mono-symmetric crosssections. Moreover, it is discussed what practical consequences can be expected due to the differences between the classic and the more realistic solutions (e.g., FEM or the here-presented advanced analytical solutions). In the paper, first, the proposed analytical model is summarized (in Section 2), then the solution is discussed (in Section 3). In Section 4 numerical solutions from the introduced analytical model are compared to those from FEM GNI analyses. In Section 5, further numerical studies are included, where capacity predictions are presented to demonstrate the importance of the differences between the classic and more advanced solutions when GNIA or GMNIA is used for capacity prediction. Finally, conclusions are drawn. 2. The analytical model 2.1. General In Section 2the applied analytical model is described. The model is developed for simply supported beams subjected to two opposite endmoments, see Fig. 1, hence the moment is uniform along the length. The cross-section is arbitrary, but it is assumed that the loading is in a principal plane, i.e., the bending takes place around a principal axis of the cross-section. The coordinate system is adjusted to the principal axes, that is Yand Zare the two principal axes, and the Xaxis is aligned to the mass center of the cross-section. The mechanical features are as follows: (i) the thin-walled member is modeled as an assembly of plane plates (referred to also as Fig. 1. Simply supported beam in uniform bending. strips), (ii) global buckling is defined as the modes with no crosssection distortion, and with no in-plane transverse extension and no in-plane shear deformations, (iii) the stress–strain field is assumed according to Kirchhoff plate theory for the out-of-plane behavior and to a conventional 2D stress–strain state for the in-plane behavior of each plate element, (iv) the material model follows Hooke’s law, (v) the end moments are applied as distributed loading, linearly varying over the (end) cross-sections. It is to mention that mechanically identical analytical models have been successfully applied e.g., in [19,20]. 2.2. Kinematics Since global buckling is defined as displacements with no crosssection distortion, the global displacements of the reference line of the member (defined by the cross-section mass centers) are expressed as follows 𝑈=𝑈𝑚𝑓𝑈(𝑋) 𝑉=𝑉𝑚𝑓𝑉(𝑋) 𝑊=𝑊𝑚𝑓𝑊(𝑋) ∅=∅𝑚𝑓∅(𝑋) (1) where Um,Vmand Wmare global translational displacement amplitudes, ∅𝑚is the torsional displacement amplitude, and the ffunctions determine the displacements’ distributions in the longitudinal direction. Though the longitudinal displacements may have some effect, as discussed e.g., in [19], this effect is typically negligible for practical cross-sections and practical length ranges, therefore, the global longitudinal displacement of the cross-sections will be neglected here, i.e., Um= 0 is assumed. Moreover, since the loading is in the X-Z plane, the primary displacement is the W. Even though there might be an interaction between the various displacement components, it is reasonable to assume that Wis due, mostly, to primary loading. The primary loading is a uniform moment; therefore a quadratic function is considered for W, in accordance with the classic first-order solution. Finally, for Vand ∅half sinewaves are considered, in accordance with classic linear buckling solutions for LTB. Therefore, the assumed displacement functions are as follows: 𝑊(𝑋) = 𝑊𝑚 4 𝐿2(𝐿−𝑋)𝑋 𝑉(𝑋) = 𝑉𝑚sin (𝜋𝑋 𝐿) ∅(𝑋)=∅𝑚sin (𝜋𝑋 𝐿)(2) where Lis the member length. It is to underline that V,Wand ∅are the total displacements, measured from the perfect, straight position. The goal here is to derive an analytical solution for GNIA, hence initial geometric imperfections are considered, too. In accordance with classic solutions, half sinewaves are applied as follows: 𝑉𝑖𝑛𝑖(𝑋) = 𝑉𝑚,𝑖𝑛𝑖 sin (𝜋𝑋 𝐿),∅𝑖𝑛𝑖(𝑋)=∅𝑚,𝑖𝑛𝑖 sin (𝜋𝑋 𝐿)(3) 2 M.Z. Haffar, M. Horáček and S. Ádány Thin-Walled Structures 184 (2023) 110535 Fig. 2. Global/local coordinates and displacements. The global displacements of the member fully determine the displacements of the strips, too, because the cross-section is rigid. By considering that the global Xaxis is aligned to the mass center of the cross-section, and the local and global Xand xare parallel with each other, the displacements of the strips’ (longitudinal) mid-lines are expressed as follows: 𝑢𝑚,𝑖(𝑥)=−𝜕𝑉 𝜕𝑋 𝑌𝑚,𝑖 −𝜕𝑊 𝜕𝑋 𝑍𝑚,𝑖 +𝜕∅ 𝜕𝑋 𝜔𝑚,𝑖 𝑣𝑚,𝑖(𝑥) = 𝑉 𝑐𝑜𝑠𝛼𝑖+𝑊 𝑠𝑖𝑛𝛼𝑖+ ∅ ((𝑌𝑚,𝑖 −𝑌𝑆)𝑠𝑖𝑛𝛼𝑖−(𝑍𝑚,𝑖 −𝑍𝑆)𝑐𝑜𝑠𝛼𝑖) 𝑤𝑚,𝑖(𝑥)=−𝑉 𝑠𝑖𝑛𝛼𝑖+𝑊 𝑐𝑜𝑠𝛼𝑖+ ∅ ((𝑌𝑚,𝑖 −𝑌𝑆)𝑐𝑜𝑠𝛼𝑖+(𝑍𝑚,𝑖 −𝑍𝑆)𝑠𝑖𝑛𝛼𝑖) 𝜑𝑚,𝑖(𝑥)=∅ (4) where Y𝑚,𝑖 and Z𝑚,𝑖 are the global coordinates of the ith strip’s midpoint, YSand ZSare the global coordinates of the shear center of the cross-section, 𝜔𝑚,𝑖 is the sectoral coordinate (with respect to the shear center) at the location of the ith strip mid-point, and 𝛼𝑖of the local y-axis of the ith strip to the global Y-axis of the member, see Fig. 2. Note that formulae for the calculation of shear center and sectoral coordinates can be found in textbooks, as well as a good summary is given in the Eurocode for cold-formed steel, see Annex C of EN 1993-1-3:2006 [13]. By using the above displacements in the mid-point of the strips, the local displacement functions of the strips can be constructed as follows: 𝑢𝑖(𝑥, 𝑦, 𝑧) = 𝑢𝑚,𝑖 −𝑣𝑚,𝑖𝑦−𝑤𝑚,𝑖𝑧−𝜑𝑚,𝑖𝑦𝑧 𝑣𝑖(𝑥, 𝑦, 𝑧) = 𝑣𝑚,𝑖 −𝜑𝑚,𝑖𝑧 𝑤𝑖(𝑥, 𝑦) = 𝑤𝑚,𝑖 +𝜑𝑚,𝑖𝑦 𝜑𝑖(𝑥) = 𝜑𝑚,𝑖 (5) Thus, the u(x,y,z), v(x,y,z), w(x,y) and 𝜑(x) local displacement functions are expressed by the Vm,Wmand ∅𝑚global displacement amplitudes. 2.3. Total potential For the solution the energy method is followed. Therefore, the internal potential (i.e., accumulated elastic strain energy) and external potential (i.e., the negative of the work done by the loads) should be expressed. However, it is assumed that the load is applied in increments, hence the increment of the total potential is necessary, expressed as the sum of the internal potential increment and the external potential increment, as follows: 𝛥𝛱 =𝛥𝛱𝑖𝑛𝑡 +𝛥𝛱𝑒𝑥𝑡 (6) After ‘j-1’ increments, the load value is 𝑀𝑌 𝑎, and the corresponding displacement amplitudes are 𝑊𝑚𝑎,𝑉𝑚𝑎 and ∅𝑚𝑎. This state is an equilibrium state, and it is referred to as state ‘a’. The goal is to find the displacement increments 𝛥𝑊𝑚𝑗,𝛥𝑉𝑚𝑗 and 𝛥∅𝑚𝑗 in the 𝑗th incremental step as the load is further increased by 𝛥𝑀𝑌, that is when the member reaches the next equilibrium state, referred to as state ‘b’. The load and the displacements at the end of the incremental step are, therefore: 𝑀𝑌 𝑏 =𝑀𝑌 𝑎 +𝛥𝑀𝑌𝑊𝑚𝑏 =𝑊𝑚𝑎 +𝛥𝑊𝑚𝑗 𝑉𝑚𝑏 =𝑉𝑚𝑎 +𝛥𝑉𝑚𝑗 ∅𝑚𝑏 = ∅𝑚𝑎 +𝛥∅𝑚𝑗 (7) 2.3.1. Internal potential increment The strain energy can readily be expressed in both states ‘a’ and ‘b’, and then the energy increment is simply the difference. Mathematically, 𝛥𝛱𝑖𝑛𝑡 =𝛱𝑖𝑛𝑡,𝑏 −𝛱𝑖𝑛𝑡,𝑎 (8) with 𝛱𝑖𝑛𝑡,𝑎 =1 2 𝑛 ∑ 𝑖=1 ∫𝑉 𝐸𝑡𝑖(𝜕𝑢𝑖𝑎 𝜕𝑥 )2 +𝐸𝑡𝑖3 12 (𝜕2𝑤𝑖𝑎 𝜕𝑥2)2 +𝐺𝑡𝑖3 3(𝜕2𝑤𝑖𝑎 𝜕𝑥𝜕𝑦 )2 d𝐴 𝛱𝑖𝑛𝑡,𝑏 =1 2 𝑛 ∑ 𝑖=1 ∫𝑉 𝐸𝑡𝑖(𝜕(𝑢𝑖𝑎 +𝛥𝑢𝑖𝑗 ) 𝜕𝑥 )2 +𝐸𝑡𝑖3 12 (𝜕2(𝑤𝑖𝑎 +𝛥𝑤𝑖𝑗) 𝜕𝑥2)2 +𝐺𝑡𝑖3 3(𝜕2(𝑤𝑖𝑎 +𝛥𝑤𝑖𝑗) 𝜕𝑥𝜕𝑦 )2 d𝐴 (9) where 𝑢𝑖𝑎 and 𝑤𝑖𝑎 are the displacement functions of the ith strip at state ‘a’, and 𝛥𝑢𝑖𝑗 and 𝛥𝑤𝑖𝑗 are displacement increment functions of the ith strip in the jth load incremental step. These local displacement functions are linked to the global displacement functions through Eqs. (4)–(5). Moreover, in Eq. (9) the integral, in fact, means double integration with respect to xand y, for the whole surface area of the strip, i.e., xis taken from 0 to L, and yis taken from (-b𝑖/2) to (+b𝑖 /2). Moreover, Lis the member length, b𝑖and t𝑖are the width and thickness of the ith strip, respectively, nis the number of strips, and E and Gare the modulus of elasticity and the shear modulus, respectively. It is to observe that the strain energy increment is expressed with respect to the global displacement increments 𝛥𝑉𝑚𝑗,𝛥𝑊𝑚𝑗 and 𝛥∅𝑚𝑗 (and obviously, is also dependent on the current displacements). 2.3.2. External potential increment For the external potential, the generic expression is as follows 𝛱𝑒𝑥𝑡 = − 𝑛 ∑ 𝑖=1 ∫𝑉 𝜎𝑥,𝑖𝜀𝑥,𝑖d𝐴(10) In this expression 𝜎𝑥,𝑖 is the longitudinal normal stress function for the ith strip, 𝜀𝑥,𝑖 is the corresponding strain function (i.e., longitudinal normal strain). Hence, it is implicitly assumed here that the other stress and strain components are negligibly small, which is a reasonable and classic assumption for global buckling problems. Still, the exact meaning of 𝜎𝑥,𝑖 and 𝜀𝑥,𝑖 needs further consideration. Since we want a geometrically non-linear analysis, we need to consider linear and non-linear strain terms. At a certain state, e.g., state ‘a’ and ‘b’, the linear, first-order strain function for the ith strip is expressed as: 𝜀𝑥,𝑖𝑎𝐼=𝜕𝑢𝑖𝑎 𝜕𝑥 𝜀𝑥,𝑖𝑏𝐼=𝜕𝑢𝑖𝑏 𝜕𝑥 =𝜕𝑢𝑖𝑎 𝜕𝑥 +𝜕𝛥𝑢𝑖𝑗 𝜕𝑥 (11) Hence the strain increment is: 𝛥𝜀𝑥,𝑖𝑗𝐼=𝜀𝑥,𝑖𝑏𝐼−𝜀𝑥,𝑖𝑎𝐼=𝜕𝛥𝑢𝑖𝑗 𝜕𝑥 (12) For the non-linear part the Green–Lagrange strain tensor is applied. The strain in the ith strip can be expressed as: 𝜀𝑥,𝑖𝑎𝐼𝐼 =1 2[(𝜕𝑣𝑖𝑎 𝜕𝑥 )2 +(𝜕𝑤𝑖𝑎 𝜕𝑥 )2] 𝜀𝑥,𝑖𝑏𝐼𝐼 =1 2[(𝜕(𝑣𝑖𝑎 +𝛥𝑣𝑖𝑗) 𝜕𝑥 )2 +(𝜕(𝑤𝑖𝑎 +𝛥𝑤𝑖𝑗) 𝜕𝑥 )2](13) 3 M.Z. Haffar, M. Horáček and S. Ádány Thin-Walled Structures 184 (2023) 110535 from which the non-linear strain increment can be expressed as: 𝛥𝜀𝑥,𝑖𝐼𝐼 =1 2(𝜕𝛥𝑣𝑖𝑗 𝜕𝑥 )2 +𝜕𝛥𝑣𝑖𝑗 𝜕𝑥 𝜕𝑣𝑖𝑎 𝜕𝑥 +1 2(𝜕𝛥𝑤𝑖𝑗 𝜕𝑥 )2 +𝜕𝛥𝑤𝑖𝑗 𝜕𝑥 𝜕𝑤𝑖𝑎 𝜕𝑥 (14) where 𝑣𝑖𝑎 and 𝑤𝑖𝑎 are the displacement functions of the ith strip at state ‘a’, and 𝛥𝑣 and 𝛥𝑤𝑖𝑗 are displacement increment functions of the ith strip in the jth load incremental step; and these local displacement functions can be derived from the global displacement functions by using Eqs. (4)–(5). The stress in the member, therefore the stress in each strip is dependent on the loading, as well as on whether we consider the changed stress state as the member deforms. If the member is subjected to uniaxial bending, the primary, first-order stress is linearly varying with Z. For the ith strip, at state ‘b’, it can be expressed as: 𝜎𝑥,𝑖𝑏𝐼= − 𝑀𝑌 𝑎 +𝛥𝑀𝑌 𝐼𝑌(𝑍𝑚,𝑖 +𝑦⋅𝑠𝑖𝑛 (𝛼𝑖)+𝑧⋅𝑐𝑜𝑠 (𝛼𝑖)) (15) where 𝐼𝑌is the second moment of area calculated for the global Y-axis. Since a load incremental procedure is considered here, where in each incremental step the actual displacements are calculated and the stresses are updated, we need to consider the effect of changing geometry, i.e., the second-order stresses due to 𝑉and ∅. To determine these second-order stresses, we utilize classic differential equations, namely Eqs. (5–78) and (5–80) of [1]. Eqs. (5–78) and of [1] can be written (by using the notations of the actual paper) as follows: 𝐸𝐼𝑧 𝜕4(𝑉−𝑉𝑖𝑛𝑖) 𝜕𝑋4−𝑀𝑌 𝜕2∅ 𝜕𝑋2= 0 (16) where 𝐼𝑍is the second moment of area calculated for the global Z-axis. This equation expresses the moment equilibrium in the lateral direction (i.e., in the plane perpendicular to the plane of loading). The first term represents the stress resultant, which, if an initial geometric imperfection is present, is to be calculated from the (𝑉−𝑉𝑖𝑛𝑖)displacement. The second term represents the lateral component of the M𝑌as the cross-section is twisted. Considering the assumed displacement functions of Eqs. (2)–(3), and after simplification, we get: (𝑉𝑚−𝑉𝑚,𝑖𝑛𝑖)𝜋2 𝐿2𝐸𝐼𝑍sin (𝜋𝑋 𝐿)+𝑀𝑌∅𝑚sin (𝜋𝑋 𝐿)= 0 (17) The resultant of the stresses for the lateral bending can be expressed as: 𝑀𝑍=𝐸𝐼𝑍 𝜕2(𝑉−𝑉𝑖𝑛𝑖) 𝜕𝑋2(18) or, with considering the shape functions of Eqs. (2)–(3): 𝑀𝑍= − (𝑉𝑚−𝑉𝑚,𝑖𝑛𝑖)𝜋2 𝐿2𝐸𝐼𝑍sin (𝜋𝑋 𝐿)(19) Let us substitute the above equation into Eq. (17), and we get the lateral bending moment due to the twisting rotation of the cross-section as follows: 𝑀𝑍=𝑀𝑌∅(20) Thus, when the member is twisted, 𝑀𝑌generates 𝑀𝑍, i.e. the primary bending generates lateral bending. From the lateral bending the stress in the 𝑖th strip, at state ‘b’, is expressed as: 𝜎𝑥,𝑖𝑏𝐼𝐼,𝑙𝑎𝑡𝑒𝑟𝑎𝑙 = − (𝑀𝑌 𝑎 +𝛥𝑀𝑌)∅𝑚𝑎 𝐼𝑍(𝑌𝑚,𝑖 +𝑦⋅cos (𝛼𝑖) +𝑧⋅sin (𝛼𝑖))sin (𝜋𝑋 𝐿)(21) Similarly, if the member is laterally displaced, M𝑌generates bimoment. Eq. (5–80) of [1] can be written (by using the notations of the actual paper) as: 𝐸𝐼𝑤 𝜕4(∅ − ∅𝑖𝑛𝑖) 𝜕𝑋4−𝐺𝐼𝑡 𝜕2(∅ − ∅𝑖𝑛𝑖) 𝜕𝑋2− 2𝛽𝑍𝑀𝑌 𝜕2∅ 𝜕𝑋2−𝑀𝑌 𝜕2V 𝜕𝑋2= 0 (22) where 𝐼𝑤is the warping constant, 𝐼𝑡is the torsion constant, and 𝛽𝑍is a non-symmetry constant, defined by Eq. (38). Eq. (22) expresses the equilibrium of torsional moments. The first and second terms represent the stress resultants from warping torsion (i.e., resultant of warping normal stresses) and Saint-Venant torsion (i.e., resultants of shear stresses), respectively. If initial geometric imperfection is present, these terms should be calculated from the (∅ − ∅𝑖𝑛𝑖)displacement. The third and fourth terms represent the torsional moment component of the external M𝑌moment, due to the displacements. Both 𝑉and ∅are assumed to have a half sinewave longitudinally, and the same for the initial shapes. By considering these functions we get: 𝐸𝐼𝑤 𝜋2 𝐿2(∅𝑚− ∅𝑚,𝑖𝑛𝑖)+𝐺𝐼𝑡(∅𝑚− ∅𝑚,𝑖𝑛𝑖)+ 2𝛽𝑍𝑀𝑌∅𝑚+𝑀𝑌𝑉𝑚= 0 (23) from which: (∅𝑚− ∅𝑚,𝑖𝑛𝑖)=−𝑀𝑌𝑉𝑚− 2𝛽𝑍𝑀𝑌∅𝑚 𝐸𝐼𝑤 𝜋2 𝐿2+𝐺𝐼𝑡 (24) By introducing: 𝐹𝑡=𝐺𝐼𝑡and 𝐹𝑤=𝜋2𝐸𝐼𝑤 𝐿2(25) the Eq. (24) can be written as: (∅𝑚− ∅𝑚,𝑖𝑛𝑖)=−𝑀𝑌𝑉𝑚− 2𝛽𝑍𝑀𝑌∅𝑚 𝐹𝑤+𝐹𝑡 (26) It is known that the bimoment is the torsional resultant of the warping normal stresses, therefore: 𝐵=𝐸𝐼𝑤 𝜕2(∅ − ∅𝑖𝑛𝑖) 𝜕𝑋2(27) from which, by using Eqs. (26)–(27), the bimoment amplitude can be expressed as: 𝐵𝑚= −𝐸𝐼𝑤 𝜋2 𝐿2(∅𝑚− ∅𝑚,𝑖𝑛𝑖)= −𝐹𝑤(∅𝑚− ∅𝑚,𝑖𝑛𝑖)(28) Substituting Eq. (26) into Eq. (28) we get the bimoment amplitude which represents the influence of the second-order displacements, as follows: 𝐵𝑚=𝐹𝑤 𝑀𝑌𝑉𝑚+ 2𝛽𝑍𝑀𝑌∅𝑚 𝐹𝑤+𝐹𝑡 =𝑀𝑌(𝑉𝑚+ 2𝛽𝑍∅𝑚)1 1 + 𝐹𝑡∕𝐹𝑤 (29) From the bimoment the stress in the ith strip can readily be calculated, e.g., at state ‘b’ it is: 𝜎𝑥,𝑖𝑏𝐼𝐼,𝑡𝑤𝑖𝑠𝑡 = −(𝑀𝑌 𝑎 +𝛥𝑀𝑌)𝑉𝑚𝑎 + 2𝛽𝑍∅𝑚𝑎 𝐼𝑤 1 1 + 𝐹𝑡∕𝐹𝑤 𝜔𝑖(𝑦, 𝑧) sin (𝜋𝑋 𝐿) (30) where 𝜔𝑖is the sectoral coordinate function for the ith strip. Thus, with all the above considerations, the external potential increment can be expressed as follows: 𝛥𝛱𝑒𝑥𝑡 = − 𝑛 ∑ 𝑖=1 ∫𝑉(𝜎𝑥,𝑖𝑏𝐼)(𝛥𝜀𝑥,𝑖𝐼+𝛥𝜀𝑥,𝑖𝐼𝐼 )d𝐴− 𝑛 ∑ 𝑖=1 ∫𝑉(𝜎𝑥,𝑖𝑏𝐼𝐼,𝑙𝑎𝑡𝑒𝑟𝑎𝑙 +𝜎𝑥,𝑖𝑏𝐼𝐼,𝑡𝑤𝑖𝑠𝑡)(𝛥𝜀𝑥,𝑖𝐼𝐼 )d𝐴 (31) The first term shows the work of the primary loading, while the second term is due to the secondary stresses, i.e., the stresses induced by the secondary displacements. It is to observe that the work increment, again, is expressed with respect to the global displacement increments 𝛥𝑉𝑚𝑗,𝛥𝑊𝑚𝑗 and 𝛥∅𝑚𝑗, hence the increment of the total potential is expressed by the displacement increments. 4 M.Z. Haffar, M. Horáček and S. Ádány Thin-Walled Structures 184 (2023) 110535 2.4. Equilibrium equation In equilibrium the total potential is stationary, therefore: 𝜕𝛥𝛱 𝜕𝛥𝑊𝑚𝑗 = 0 𝜕𝛥𝛱 𝜕𝛥𝑉𝑚𝑗 = 0 𝜕𝛥𝛱 𝜕𝛥∅𝑚𝑗 = 0 (32) Eq. (32) is a system of three equations, which can be summarized into one single matrix equation as follows: 𝐊𝐞𝜟𝐝+𝑀𝑌 𝑏𝐊𝐠𝜟𝐝=𝜟𝐟(33) where the first matrix is the 𝐊𝐞elastic stiffness matrix, while the second one is the 𝐊𝐠geometric stiffness matrix of the problem: 𝐊𝐞=⎡⎢⎢⎢⎢⎢⎣ 𝐹𝑌 64 𝜋2𝐿0 0 0𝐹𝑍 𝜋2 2𝐿0 0 0 𝐹𝑋 𝜋2 2𝐿 ⎤⎥⎥⎥⎥⎥⎦ 𝐊𝐠=⎡⎢⎢⎢⎢⎢⎣ 0 0 −∅𝑚𝑎 2 𝐿 0 0 𝜋2 2𝐿 −∅𝑚𝑎 2 𝐿 𝜋2 2𝐿 𝜋2 𝐿𝛽𝑍+ ∅𝑚𝑎 4𝜋 3𝐿𝛽𝑌 ⎤⎥⎥⎥⎥⎥⎦ (34) Moreover, the displacement vector contains the displacement increments for the given jth load step 𝜟𝐝=⎡⎢⎢⎢⎢⎣ 𝛥𝑊𝑚𝑗 𝛥𝑉𝑚𝑗 𝛥∅𝑚𝑗 ⎤⎥⎥⎥⎥⎦ (35) and the load vector at the right-hand-side of the equation is defined as: 𝜟𝐟=⎡⎢⎢⎢⎢⎢⎢⎢⎣ −𝑀𝑦𝑏 8 𝐿+𝑀𝑦𝑏∅2 𝑚𝑎 2 𝐿−𝑊𝑚𝑎𝐹𝑌 64 𝜋2𝐿 −𝑀𝑦𝑏∅𝑚𝑎 𝜋2 2𝐿− (𝑉𝑚𝑎 −𝑉𝑚,𝑖𝑛𝑖)𝐹𝑍 𝜋2 2𝐿 −𝑀𝑦𝑏𝑉𝑚𝑎 𝜋2 2𝐿−𝑀𝑦𝑏∅𝑚𝑎 𝜋2 𝐿𝛽𝑍−𝑀𝑦𝑏∅𝑚𝑎 4𝜋 3𝐿𝛽𝑌 +𝑀𝑦𝑏𝑊𝑚𝑎∅𝑚𝑎 2 𝐿− (∅𝑚𝑎 − ∅𝑚,𝑖𝑛𝑖)𝐹𝑋 𝜋2 2𝐿 ⎤⎥⎥⎥⎥⎥⎥⎥⎦ (36) In the above matrices and vectors, the 𝐹𝑥,𝐹𝑦and 𝐹𝑧symbols are defined as follows: 𝐹𝑌=𝜋2𝐸𝐼𝑌 𝐿2, 𝐹𝑍=𝜋2𝐸𝐼𝑍 𝐿2, 𝐹𝑋=𝐹𝑡+𝐹𝑤=𝐺𝐼𝑡+𝜋2𝐸𝐼𝑤 𝐿2(37) where non-symmetry constants are defined as: 𝛽𝑌=𝑦𝑆− 0.5∫𝐴(𝑦2+𝑧2)𝑦d𝐴 𝛽𝑍=𝑧𝑆− 0.5∫𝐴(𝑦2+𝑧2)𝑧d𝐴 (38) where 𝑦𝑆, 𝑧𝑆are the coordinates of the shear center relative to the centroid. The above equation system can be solved analytically. By introducing  𝑀2 𝑐𝑟 =𝑀𝑐𝑟2− 2𝑀𝑐𝑟𝐹𝑍𝛽𝑍(39) the resulting displacement increments are given as in Box I 3. Discussion of the analytical solution 3.1. LBA Based on the above derivation, it is easy to get the LBA problem, by simplification. In LBA there is no initial imperfection, only the primary stresses and only the second-order strain terms are to be considered. Moreover, the first equation in Eq. (33) is not necessary, since in the solution the primary displacements have no role. Finally, no load increments are necessary, i.e., it is enough to focus on the first step of the incremental procedure; consequently: the 𝛥𝑉𝑚𝑗 and 𝛥∅𝑚𝑗 displacements are simply the total displacements 𝑉𝑚and ∅𝑚, and 𝛥𝑀𝑌is simply the total 𝑀𝑌. The equation system in Eq. (33) therefore simplifies to the following one: [𝐹𝑍0 0𝐹𝑋].[𝑉𝑚 ∅𝑚]+[0𝑀𝑌 𝑀𝑌2𝑀𝑌𝛽𝑍].[𝑉𝑚 ∅𝑚]=[0 0](43) The solution exists if 𝑀𝑌takes a specific value, the 𝑀𝑐𝑟 critical moment, which is expressed as: 𝑀𝑐𝑟 =𝐹𝑍𝛽𝑍±√(𝐹𝑍𝛽𝑍)2+𝐹𝑋𝐹𝑍(44) By back-substitution, we can find the buckled shape. Obviously, the 𝑉(𝑋)and ∅(X)displacement functions have a half-sinewave shape longitudinally. Since the 𝐊𝐞+𝐊𝐠matrix is rank-deficient in the critical state, i.e., when 𝑀𝑌=𝑀𝑐𝑟, the 𝑉𝑚and ∅𝑚amplitudes are dependent on each other, otherwise the amplitudes are arbitrary. The relationship between 𝑉𝑚and ∅𝑚is as follows 𝑉𝑚= −∅𝑚 𝑀𝑐𝑟 𝐹𝑍 (45) 3.2. Classic GNIA Classic GNIA solution can also be obtained by simplifying Eq. (33). Only the primary stresses and only the second-order strain terms are to be considered. Moreover, the first equation in Eq. (33) is not necessary. If the initial geometric imperfection is the buckling shape from the LBA, 𝑉𝑚,𝑖𝑛𝑖 and ∅𝑚,𝑖𝑛𝑖 are dependent on each other, see Eq. (45). Finally, no load increments are necessary, i.e., it is enough to focus on the first step of the incremental procedure. Here, to distinguish between the initial displacement and the displacement increment, we keep the notation for the 𝛥𝑉𝑚and 𝛥∅𝑚displacement increments, however, the moment increment 𝛥𝑀𝑌is simply the total 𝑀𝑌. Finally, the equation system of Eq. (33) simplifies to the following one: [𝐹𝑍0 0𝐹𝑋].[𝛥𝑉𝑚 𝛥∅𝑚]+[0𝑀𝑌 𝑀𝑌2𝑀𝑌𝛽𝑍].[𝛥𝑉𝑚 𝛥∅𝑚] =[−𝑀𝑌∅𝑚,𝑖𝑛𝑖 −𝑀𝑌𝑉𝑚,𝑖𝑛𝑖 − 2𝑀𝑌𝛽𝑍∅𝑚,𝑖𝑛𝑖](46) The solution for the displacement increments: 𝛥𝑉𝑚=𝑉𝑚,𝑖𝑛𝑖 1 𝑀𝑐𝑟∕𝑀𝑌− 1 𝛥∅𝑚= ∅𝑚,𝑖𝑛𝑖 1 𝑀𝑐𝑟∕𝑀𝑌− 1 (47) The total displacement is the sum of the initial one and the increment, which leads to the well-known formulae expressing the displacement amplification (also known as the Young formula): 𝑉𝑚=𝛥𝑉𝑚+𝑉𝑚,𝑖𝑛𝑖 =𝑉𝑚,𝑖𝑛𝑖 1 1 − 𝑀𝑌∕𝑀𝑐𝑟 ∅𝑚=𝛥∅𝑚+ ∅𝑚,𝑖𝑛𝑖 = ∅𝑚,𝑖𝑛𝑖 1 1 − 𝑀𝑌∕𝑀𝑐𝑟 (48) It is observed that the ratio of Vand ∅remains the same during the analysis. 3.3. Doubly symmetrical cross-sections The expressions for the displacement increments, see Eqs. (40)– (42) can be simplified if the cross-section is doubly symmetric, i.e., if 𝛽𝑌=𝛽𝑍= 0. This also means that:  𝑀𝑐𝑟 =𝑀𝑐𝑟 = ±√𝐹𝑋𝐹𝑍(49) It is to observe that the second-order lateral translation and twisting rotation are different from those predicted by the Young formula. 5 M.Z. Haffar, M. Horáček and S. Ádány Thin-Walled Structures 184 (2023) 110535 𝛥𝑊𝑚= −𝑀𝑌 𝑏 𝜋2 8𝐹𝑌[ 𝑀2 𝑐𝑟 𝑀𝑌 𝑏2− 1]−𝑉𝑚,𝑖𝑛𝑖∅𝑚𝑎 𝜋2𝐹𝑍 32𝐹𝑌+ ∅𝑚,𝑖𝑛𝑖∅𝑚𝑎  𝑀2 𝑐𝑟 𝑀𝑌 𝑏 𝜋2 32𝐹𝑌− ∅𝑚𝑎 𝜋𝐹𝑍𝛽𝑌 3𝐹𝑌−𝐹𝑍𝛽𝑍𝜋2 4𝐹𝑌  𝑀2 𝑐𝑟 𝑀𝑌 𝑏2− 1 + ∅𝑚𝑎 8𝐹𝑍𝛽𝑌 3𝜋𝑀𝑌 𝑏 +2𝐹𝑍𝛽𝑍 𝑀𝑌 𝑏 − ∅2 𝑚𝑎 𝐹𝑍 8𝐹𝑌 −𝑊𝑚𝑎 (40) 𝛥𝑉𝑚= 𝑉𝑚,𝑖𝑛𝑖  𝑀2 𝑐𝑟 𝑀𝑌 𝑏2+𝑉𝑚,𝑖𝑛𝑖 2𝐹𝑍𝛽𝑍 𝑀𝑌 𝑏 + ∅𝑚𝑎 𝑀𝑌 𝑏 2𝐹𝑌− ∅𝑚,𝑖𝑛𝑖 𝐹𝑋 𝑀𝑌 𝑏 +𝑉𝑚,𝑖𝑛𝑖∅𝑚𝑎 8𝐹𝑍𝛽𝑌 3𝜋𝑀𝑌 𝑏 −𝑉𝑚,𝑖𝑛𝑖∅2 𝑚𝑎 𝐹𝑍 8𝐹𝑌  𝑀2 𝑐𝑟 𝑀𝑌 𝑏2− 1 + ∅𝑚𝑎 8𝐹𝑍𝛽𝑌 3𝜋𝑀𝑌 𝑏 +2𝐹𝑍𝛽𝑍 𝑀𝑌 𝑏 − ∅2 𝑚𝑎 𝐹𝑍 8𝐹𝑌 −𝑉𝑚𝑎 (41) 𝛥∅𝑚= −𝑉𝑚,𝑖𝑛𝑖 𝐹𝑧 𝑀𝑌 𝑏 − ∅𝑚𝑎 𝐹𝑧 2𝐹𝑦+ ∅𝑚,𝑖𝑛𝑖 𝑀𝑐𝑟2 𝑀𝑌 𝑏2  𝑀2 𝑐𝑟 𝑀𝑌 𝑏2− 1 + ∅𝑚𝑎 8𝐹𝑍𝛽𝑌 3𝜋𝑀𝑌 𝑏 +2𝐹𝑍𝛽𝑍 𝑀𝑌 𝑏 − ∅2 𝑚𝑎 𝐹𝑍 8𝐹𝑌 − ∅𝑚𝑎 (42) Box I. Moreover, the formulae for 𝑉and ∅are different. Independently of whether the initial shape is buckling shape or not, the amplification for 𝑉and ∅are different. An important characteristic of these formulae is that if the sign of the initial geometry is reversed then (i) the vertical 𝑊displacement is unchanged, (ii) the sign of the lateral increment is reversed, and (iii) the sign of the twisting rotation increment is reversed. This also means that a symmetric bifurcation is predicted (as the initial displacement converges to zero). As the bending moment increases, the denominator of the formulae can decrease to zero, which identifies singularity. The corresponding moment can be expressed as 𝑀𝑌 ,𝑠𝑖𝑛𝑔 =𝑀𝑐𝑟 √1+∅2 𝑚𝑎 𝐹𝑍 8𝐹𝑌 (50) From the formulae, it can be seen that the singularity belongs to a bending moment smaller than 𝑀𝑐𝑟. The distance of the singularity to 𝑀𝑐𝑟 is largely dependent on the twisting rotation. However, since the whole theory assumes that the rotations are small, the bending moment where singularity happens is only marginally smaller than 𝑀𝑐𝑟. Still, an important difference is that while the classic analytical solution suggests that singularity happens only at infinitely large displacements, the new analytical formulae predict the singularity at finite displacements. 3.4. Mono-symmetric cross-sections, bending in the symmetry plane Now we consider a mono-symmetric cross-section. The axis of symmetry is the Z-axis, i.e., 𝛽𝑌= 0, and the bending is acting in the symmetry plane. The expressions in Eqs. (40)–(42) for the displacement increments can slightly be simplified. By analyzing the formulae for the lateral translation and for the twisting rotation, it can again be observed that they are different from the Young formula and from each other. It can also be understood that the sign of the initial geometry has no real effect, i.e., a symmetric bifurcation is predicted (as the initial displacement converges to zero) Singularity happens when the denominator is equal to zero. The bending moment which causes singularity is the solution of a quadratic equation as follows: 𝑀𝑌 ,𝑠𝑖𝑛𝑔2(1+∅2 𝑚𝑎 𝐹𝑍 8𝐹𝑌)−𝑀𝑌 ,𝑠𝑖𝑛𝑔2𝐹𝑍𝛽𝑍−(𝑀𝑐𝑟2− 2𝑀𝑐𝑟𝐹𝑍𝛽𝑍)= 0 (51) from which the bending moment can be expressed as: 𝑀𝑌 ,𝑠𝑖𝑛𝑔 = 2𝐹𝑍𝛽𝑍±√(2𝐹𝑍𝛽𝑍)2+ 4 (1+∅2 𝑚𝑎 𝐹𝑍 8𝐹𝑌)(𝑀𝑐𝑟2− 2𝑀𝑐𝑟𝐹𝑍𝛽𝑍) 2(1+∅2 𝑚𝑎 𝐹𝑍 8𝐹𝑌)(52) However, if the twist angle is small, the (1+∅2 𝑚𝑎 𝐹𝑍 8𝐹𝑌)term is close to one. With this approximation, the above formula can be simplified to: 𝑀𝑌 ,𝑠𝑖𝑛𝑔 ≅𝑀𝑐𝑟 (53) That is, the singularity belongs approximately to 𝑀𝑐𝑟. 3.5. Mono-symmetric cross-sections, bending perpendicular to the symmetry plane Now we consider mono-symmetric cross-sections where the axis of symmetry is the Y-axis, i.e., 𝛽𝑍= 0. The bending is still in the vertical plane, that is in a plane perpendicular to the symmetry plane. The expressions for the displacement increments in Eqs. (40)–(42) can be simplified. The critical moment formula is identical to that for doubly symmetric cross-sections, see Eq. (49). Partly similar observations can be made as above. However, now there is a significant difference: the expression in the denominators of the formulae is dependent on the sign of the initial displacement; and not only the sign of the expression, but also its absolute value. This means that the magnitudes of both the primary (vertical) displacement and the secondary displacement are influenced by whether the initial imperfection is considered with positive or negative sign. This identifies asymmetric bifurcation if the initial imperfection converges to zero. As the bending moment increases, the denominator of the formulae can decrease to zero, which identifies singularity. Mathematically, the corresponding bending moment is the solution of a quadratic equation: 𝑀𝑌 ,𝑠𝑖𝑛𝑔2(1+∅2 𝑚𝑎 𝐹𝑍 8𝐹𝑌)−𝑀𝑌 ,𝑠𝑖𝑛𝑔∅𝑚𝑎 8𝐹𝑍𝛽𝑌 3𝜋−𝑀𝑐𝑟2= 0 (54) from which the bending moment is 𝑀𝑌 ,𝑠𝑖𝑛𝑔 = ∅𝑚𝑎 8𝐹𝑍𝛽𝑌 3𝜋±√(∅𝑚𝑎 8𝐹𝑍𝛽𝑌 3𝜋)2 + 4 (1+∅2 𝑚𝑎 𝐹𝑍 8𝐹𝑌)𝑀𝑐𝑟2 2(1+∅2 𝑚𝑎 𝐹𝑍 8𝐹𝑌)(55) Even if the ∅2 𝑚𝑎 term is neglected, we have: 𝑀𝑌 ,𝑠𝑖𝑛𝑔 ≅ ∅𝑚𝑎 2𝐹𝑍𝛽𝑌 3𝜋±√(∅𝑚𝑎 2𝐹𝑍𝛽𝑌 3𝜋)2 +𝑀𝑐𝑟2(56) The positive solution can be obtained by using the positive sign in the above formula. Moreover, if ∅𝑚𝑎 = 0, then 𝑀𝑌 ,𝑠𝑖𝑛𝑔 =𝑀𝑐𝑟. If ∅𝑚𝑎 >0,𝑀𝑌 ,𝑠𝑖𝑛𝑔 is clearly larger than 𝑀𝑐𝑟; the larger the ∅𝑚𝑎, the larger the difference between 𝑀𝑌 ,𝑠𝑖𝑛𝑔 and 𝑀𝑐𝑟. If ∅𝑚𝑎 <0, then 𝑀𝑌 ,𝑠𝑖𝑛𝑔 can be larger or smaller than 𝑀𝑐𝑟, depending on the magnitude of ∅𝑚𝑎. The dependence of singularity on the sign of the twisting rotation again justifies that the bifurcation problem is asymmetric in the case of mono-symmetric cross-sections loaded perpendicularly to the plane of symmetry. 6 M.Z. Haffar, M. Horáček and S. Ádány Thin-Walled Structures 184 (2023) 110535 Fig. 3. Displacement increments for various initial imperfections shapes. 3.6. The effect of the imperfection shape The above presented solution is valid for any imperfection shape (with a half-sinewave longitudinal distribution). It is worth looking at the effect of imperfection shape. However, to make the comparison of the various imperfection shapes, we focus on the first incremental step only, i.e., 𝑉𝑚𝑎 =𝑉𝑚,𝑖𝑛𝑖 and ∅𝑚𝑎 = ∅𝑚,𝑖𝑛𝑖 and 𝑀𝑌 𝑏 =𝛥𝑀𝑌. Moreover, to further simplify the final formulae, we assume a doubly symmetric cross-section, and we neglect the higher-order displacement terms in the formulae. Four cases are considered, and the obtained displacement increment formulae are summarized in Fig. 3. The important observations are as follows: •The formulae, even in this simplified case, are (with one exception) different from the Young formula. •Even if the buckling shape is used as an imperfection, the displacement amplifications are different for the lateral translation and the twisting rotation. •If one of the components of the initial geometry (i.e., either the lateral translation or the twisting rotation) is zero, still, nonzero increments are resulted even in the first load incremental step. Therefore, a single imperfection component (i.e., either the lateral translation or the twisting rotation) generates both lateral translation and twisting rotation during the non-linear analysis. It is to note that the above observations are valid for more general cases, too, without the applied simplifications, (i.e., general crosssections, and/or considering higher-order displacement terms in the formulae,) however, the resulted formulae are much more complex, that is why not shown here. 4. Numerical studies: comparison of analytical and shell FEM results 4.1. Summary of the numerical studies Numerical studies are completed by using the above-described new analytical model, comparing the results to those from shell FEM. Four cross-sections are considered, CS1 to CS4, as shown in Fig. 4, where the plate dimensions are given in mm, and the widths are interpreted as middle-line dimensions. The member length is the same in each case: Fig. 4. Cross-sections. L= 5 m. The buckling shape from LBA is applied as initial geometric imperfection. The analytical model is used in a very simple way. The load is imposed in small increments. The procedure starts from 𝑊𝑚𝑎 = 0, 𝑉𝑚𝑎 =𝑉𝑚,𝑖𝑛𝑖 and ∅𝑚𝑎 = ∅𝑚,𝑖𝑛𝑖, then 𝛥𝑀𝑌is applied, and the displacement increments are calculated by using Eqs. (40)–(42), to have the displacement values 𝑊𝑚𝑏,𝑉𝑚𝑏 and ∅𝑚𝑏, using Eq. (6). In the next incremental step these values are the starting values, and new increments are calculated, and then these steps are repeated. Since the formulae are simple, a large number of incremental steps can be applied. The experience is that 𝛥𝑀𝑌=𝑀𝑐𝑟∕100 is enough, but in the presented examples 𝛥𝑀𝑌=𝑀𝑐𝑟∕200 load increment has been applied. The shell FEM calculations have been performed by using the commercial software Abaqus [21]. Linear 4-sided shell elements (type S4R) based on the Reissner–Mindlin plate theory are used. A regular mesh is generated with a mesh size of 10 mm. Standard isotropic steel material is used, 𝐸= 210 GPa,𝐺= 80.8 GPa. The boundary conditions applied at the end sections of the beam simulate simple supports, i.e., globally and locally hinged supports but without restricting the warping. Practically, the flanges were supported vertically only, the web was supported horizontally only. At one end of 7 M.Z. Haffar, M. Horáček and S. Ádány Thin-Walled Structures 184 (2023) 110535 Fig. 5. Sample deformed shapes from shell FEM GNIA. (Only half of the member is shown.). the beam, the member was supported in the longitudinal direction at the center of gravity. The cross-sections are defined in a specific way so that the 𝑀𝑐𝑟∕𝐹𝑍 values – for the given length – are the same. The importance of this is that the ratio of the two components of the buckling shape is constant: i.e., 𝑉𝑚,𝑖𝑛𝑖∕∅𝑚,𝑖𝑛𝑖 is constant, as an immediate consequence of Eq. (45). This means that if, for example, the buckling shape is scaled so that the initial twist would be the same in each case, this automatically guarantees the equivalence of the initial lateral out-of-straightness values, too. Hence, exactly the same initial imperfection can be applied for all the analyzed members; if differences are found in the nonlinear behavior, these differences cannot be caused by the differences in initial shapes. In fact, we have considered an imperfection value with a lateral out-of-straightness amplitude equal to 20 mm, i.e., 𝐿∕250, which can be considered as a reasonable value for an equivalent geometric imperfection in a design situation. 4.2. Results Stresses and displacements are calculated for the whole load– displacement path for each cross-section. Sample deformed shapes from FEM analyses are shown in Fig. 5, where the half of each beam is shown at a load level approximately equal to the critical moment. Though the displacements are large, the figures prove that cross-section distortion is negligible, hence the obtained results represent lateral–torsional behavior. Fig. 6. Load–stress plots, doubly symmetric I-section (CS1). The stresses are plotted in Figs. 6,8and 10, for locations of the cross-sections: BL is the bottom flange left side, BR is the bottom flange 8 M.Z. Haffar, M. Horáček and S. Ádány Thin-Walled Structures 184 (2023) 110535 Fig. 7. Load–displacement plots, doubly symmetric I-section (CS1). Fig. 8. Load–stress plots, mono-symmetric I-sections (CS2 and CS3). right side, TL is the top flange left side, TR is the top flange right side. The displacements are plotted Figs. 7,9and 11. It is to note that in the vertical axis of each plot the moment is normalized by the critical moment, and for this normalization the critical moment calculated by the shell FEM model is employed. 4.3. Observations from GNI results The vertical displacement 𝑊is plotted in Fig. 7 for the doubly symmetric I section. As can be observed, the new analytical model predicts a nearly linear load–displacement relationship for the entire range of loading. For lower loads the calculated vertical translation values are practically identical to those from classic first-order theory, but very similar even for larger loads. On the contrary, the shell FEM values are slightly different, and the load–displacement plot is clearly nonlinear. The small differences at lower loads are not due to nonlinearities; they show that the global flexural stiffness of the shell model is slightly different from a classic beam model (mostly due to shear deformations which are neglected in the beam models but considered in the shell models). These observations on the primary displacements are little affected by the cross-section shape, which is why 𝑊is shown only for one cross-section. By looking at the stress plots, one can observe that the stresses are quite similar for all the cases, up to approx. 80%–90% of the critical moment. This proves that the analytical model reasonably well considers the secondary stresses due to the secondary displacements, see Eqs. (20) and (29). Similarly, the load–displacement paths calculated by the new analytical model and by the shell FEM analysis are similar in all the cases. Though non-negligible differences exist, the analytical results and shell FEM results show similar and significant deviation from the classic analytical solution (which is shown by the ‘Young’ curves). As the formulae of the new analytical model suggest (see e.g. Fig. 3), the amplification of 𝑉and ∅are not identical. The difference in the amplifications clearly depends on the cross-section shape. Practically, this difference can be regarded as small, but still, the numerical results prove that the deformed shape during the geometrically non-linear analysis is not merely the amplification of the initial imperfect shape, even if the imperfection is taken as the buckling shape. It is to underline that perfect agreement between the analytical model and the shell FEM cannot be expected due to several reasons. 9