Thin-Walled Stress recovery from one dimensional models for tapered bi-symmetric thin-walled I beams: De ﬁ ciencies in modern engineering tools and procedures

This paper highlights several issues of the procedures nowadays adopted for the recovery of cross-sections stress distribution within tapered thin-walled I beams. In particular, de ﬁ ciencies are evident even considering bi-symmetric structural elements behaving under the assumption of plane stress. In fact, analytical results available in the literature since the ﬁ rst half of the past century highlight that the continuous variation of the height of a in ﬁ nite long wedge induces shear stress distributions substantially di ﬀ erent from the ones occurring in prismatic beams. Unfortunately, this peculiarity of non-prismatic beams is neglected or treated with coarse approaches by most of the modern engineering tools and procedures, leading to inaccurate descriptions (and also severe un- derestimations) of the real stress magnitude. After a comprehensive literature review on this speci ﬁ c topic, the paper compares most common stress-recovery procedures with a new, simple, and e ﬀ ective tool derived from a recently proposed non-prismatic planar beam model. The numerical examples reported in the paper highlight that the approaches available in the literature and widely used in practice estimate the parameters of interest for practitioners with errors bigger than 50% leading therefore to unreliable results. Conversely, the herein proposed tool leads to errors smaller than 5% in all the considered cases, paving the way to a new generation of e ﬀ ective tools that practitioners can use for the design of such structural elements.


Introduction
Structural elements with variable height have been widespread in several engineering fields since the nineteenth century, because they allow to optimize strength and stiffness of structures with significant material savings. This practice became extremely popular in steel constructions since this workable material allows to manufacture structural elements with complex geometry without a significant increase of the production costs. Nowadays, on the one hand, highly challenging structures like skyscrapers, bridges, and large-span pavilions lead the structural optimization to be increasingly important whereas, on the other hand, modern production technologies like high precision cutting and automatic welding allow a further reduction of the production costs. As a consequence, the usage of non-prismatic structural elements is becoming an increasingly diffused practice with significant economical benefits. Unfortunately, the modeling of nonprismatic structural elements is a non trivial task of the design process that could impede effective usage of non-prismatic structural elements.
First of all, also considering a simple example like a staticallyloaded non-prismatic beam behaving under the assumption of plane stress, the most stressed cross-section generally does not coincide with the cross-section subjected to the maximal internal force [19]. Therefore, the design of non-prismatic beams requires a stress analysis more extensive and accurate than the one usually adopted for prismatic beams. Even more important, recent papers have highlighted that the continuous variation of the cross-section geometry of a homogeneous beam deeply influences several aspects of the beam behavior like (i) the cross-section stress distribution [8], (ii) the beam's constitutive relation [7], (iii) the stiffness [24], (iv) the stability, (v) the post-buckling behavior, (vi) the dynamic response, etc. As a consequence, the analysis of non prismatic beams can not be done using tools developed for prismatic beams but it requires specific-purpose models that carefully consider all the phenomenon occurring within the body of this special class of structural elements.
Within the so far introduced list, the most elementary but maybe also the most critical problem is the distribution of stresses within the cross-section. In fact, an accurate stress analysis is the first mandatory step for the evaluation of the cross-section resistance and, furthermore, it represents the core of the beam models together with the kinematics [6], and, therefore, it deeply influences all the other following aspects. Aiming at a better understanding of the problematic so far introduced, we start by considering several analytical results that allow us to highlight main critical points.

The infinite long wedge
Let us consider the infinite long wedge -depicted in Fig. 1-behaving under the hypothesis of plane stress. Neglecting the self-weight of the material and considering the concentrated loads applied in the wedge's vertex, it is possible to calculate the analytical solution for the 2D equilibrium equations i.e., the stress distribution within the body [ (1) As discussed by Bresler et al. [12], the analytical solution (1) highlights that only radial stresses σ r resist the concentrated forces applied at the vertex, without any contribution of the tangential shear τ rθ . Conversely, both the radial stress σ r and the tangential shear τ rθ resist the bending moment applied at the vertex.
Representing the stress with respect to the Cartesian coordinate system Oxy depicted in Fig. 1, a linear distribution of horizontal stress σ x and a quadratic distribution of shear stress τ with respect to y provide a reasonable approximation of the analytical solution (1) for small values of α [35,13], as illustrated qualitatively in Table 1. In greater detail • the radial stress distribution resulting from a load parallel to the wedge bisector leads to a constant distribution of horizontal stress and a linear (odd) distribution of shear, • the radial stress distribution resulting from a load perpendicular to the wedge bisector leads to a linear (odd) distribution of horizontal stress and a parabolic (even) distribution of shear, with vanishing magnitude on the wedge bisector, • the stresses resulting from a bending moment lead to a linear (odd) distribution of horizontal stress and a parabolic (even) distribution of shear with vanishing mean value.
The results so far introduced can be generalized to a tapered beam as illustrated by Krahula [20], Russo and Garic [33], Hodges et al. [17], Auricchio et al. [5], Beltempo et al. [8], and Balduzzi et al. [7]. Table 1 clearly highlights that the stress distribution in tapered beams is strikingly different compared to the one in a prismatic beam, where • a constant distribution of horizontal stress and a vanishing shear resist the axial load, • a linear (odd) distribution of horizontal stress and a parabolic (even) distribution of shear, with maximum on the beam axis, resist the shear-bending load, • a linear (odd) distribution of horizontal stress and a vanishing shear resist the pure bending load.
In greater detail, the horizontal stress distributions on prismatic and tapered beams coincide whereas the shear distributions are substantially different. In fact, tapered beams lead to non-vanishing shear for every load condition and, also for the shear-bending load, the maximal shear magnitude and its location within the cross-section are substantially different. Specifically, as clearly stated by Trahair and Ansourian [37], according to prismatic beams design formulas for shear-bending load, "the distribution of the shear stress is parabolic, with zero stresses at the top and bottom edges, and 1.5 times the average at the axis. For the wedge, the shear stress variation is also parabolic, but with zero stress at the axis and 3 times the average at the top and bottom edges." This means that, even in the simplest case, the design formulas for prismatic beams underestimate the maximal shear by 50%, locate the maximal shear in the wrong position, and this error occurs for every taper angle different from zero.
Even if limited to an extremely simple problem, the results so far discussed clearly demonstrate that the procedures usually adopted for the stress analysis of prismatic beams can not be effective for non prismatic beams, but even lead to dangerous underestimations.  Extending the discussion to more realistic situations we consider thin-walled non-prismatic beams, limiting our attention to bi-symmetric I shaped cross-sections with variable web height.
Despite the stress distribution within non-prismatic beams has a crucial role in determining their stiffness and strength, its importance is often underestimated both in engineering practice and research. In fact, most of literature on thin-walled tapered and non-prismatic beams focuses on the second order effects, buckling, and post-buckling behaviors that, for thin-walled structures, represents the main failure mechanism [9,2,31,3,23,38,36,26,14]. Unfortunately, in most of the so far cited papers, it is assumed that, in tapered beams, stress distributes within the cross-section as in prismatic beams, neglecting therefore also the effects of non-trivial shear stress distribution on constitutive relations, stiffness, stability, and post-buckling behavior. Consistently, also national and international design standards provide specific indications for the buckling limit loads of tapered structures, but do not provide specific indications about the presence of non-trivial shear stresses [16]. As an example, introducing the problem of the non-prismatic beam buckling and referring to the differences of stress distribution within tapered and prismatic beams, Marques [23] states that "for small tapering angles <°( 15 ) this difference is negligible and, as a result, regarding member design, the design formulas for prismatic have been extended for the case of tapered members". Clearly, this statement disagrees with the analytical results reported in Section 1.1, leads to extremely coarse and unsafe cross-section strength estimations, as will be demonstrated in the following, and maybe undercuts the effectiveness of any attempt for buckling and post-buckling analysis.
An even more critical modeling procedure widely adopted in practice is the usage of the so called "stepped" Finite Element (FE) models [30,21,22]. It consists in modeling the tapered beam with a certain number of prismatic elements (usually 10 or 20) in which the variable parameters (the area and the inertia) are approximated with their mean value within the element length. Clearly, this procedure ignores the non trivial constitutive relations and the resulting stiffness [24]. As a consequence, stepped FE leads to heavy errors in the evaluation of both the structural element stiffness and the internal forces acting in static-indeterminate structures [26]. Last but not least, the prismatic beam stress recovery is used in post-processing, leading to severe underestimations of the real stresses acting within the beam already discussed in Section 1.1 and deeply examined in the following sections.
Looking at the literature, the former approach for the recovery of stress distribution within non-prismatic beams was developed by Bleich [10] and consists in a generalization of the Jourawasky theory [18]. Aiming at evaluating the cross-section resistance, the author evaluates the shear stress in the cross-section's geometric center assuming that the shear magnitude reaches its maximum there, as in prismatic beams. Unfortunately, the assumption on the location of maximal shear is wrong, as can be seen looking at the results reported in Table 1. As a consequence, the procedure turned out to be very imprecise and has strongly been criticized by Paglietti and Carta [28,29].
Searching for an effective approach capable to determine the crosssection resistance, several books propose simplified procedures [11,12]. The main idea behind all the so far proposed approaches is that the boundary equilibrium forces principal stresses within flanges to be oriented as the flanges' plane. As a consequence, the forces resulting from stresses acting on the flanges have a vertical component that contributes to the vertical equilibrium of the cross-section. In particular, Blodgett [11] proposed a simple analysis tool which assumes that flanges resist to the bending moment, subsequently estimates the flanges contribution to the shear equilibrium according to the idea so far introduced, and finally evaluates the mean value of the shear stress within the web. Similar principles were used more recently by Rubin [32], Redmond [31] whereas Trahair and Ansourian [37] exploit the analytical solution for the infinite long wedge, but recover results similar to the ones coming from the procedures so far introduced.
More recently, Vu-Quoc and Léger [39] proposed a more refined approach for the recovery of shear stress distribution within the web of tapered I beams subjected to shear bending loads within the web plane. The proposed method accounts for the boundary equilibrium and the variation of the cross-section area, but neglects derivatives of the moment of inertia and the flanges contribution to the shear equilibrium. Nevertheless, the obtained shear stress recovery formulas was recently used by Chiorean and Marchis [14], within an higher order model for tapered beams.
Finally, Balduzzi et al. [6] have recently proposed a simple planar non-prismatic beam model capable to tackle all the aspects so far introduced. In particular, the paper proposes a planar stress representation based on a rigorous generalization of the Jourawsky theory, following a derivation path not different from the one proposed by Bleich [10]. Furthermore, it also provides instruments capable to manage accurately the effects of non trivial stress distribution on the beam stiffness and displacements. The numerical results reported in the paper indicate that the proposed model evaluates accurately both stresses and displacements.
The only remaining modeling approach available in literature is the usage of 3D models e.g., realized with highly refined FE based on 3D continuum mechanics [4,40], shell, or 3D FE. This approach is extremely accurate but also expensive from the computational point of view [27,1]. As a consequence, it is rarely used for the global modeling of real structures where engineers prefer to use 1D FEs. The only situations in which 3D models may be adopted are the enhanced checks like local stability and the enhanced design of structural details like stiffeners. To the authors knowledge, the only field in which shell or 3D FE are used for the modeling of a whole structure is the research, mainly for the validation of less expensive analysis approaches [15,39,25,3,23,40].

Paper aims and outline
According to the remarks discussed in Section 1.2, this paper aims at comparing performances and verifying effectiveness of the most important procedures for the stress recovery so far introduced. Specifically, the paper compares the following methods: • METHOD 1: prismatic beam stress recovery formulas i.e., the method adopted by the most of the commercial software for structural analysis, • METHOD 2: the simplified stress recovery proposed by Blodgett [11] that assumes constant stresses within flanges and web, but tackle the shear contribution coming from inclined flanges, • METHOD 3: the stress recovery proposed by Vu-Quoc and Léger [39] that accounts for web boundary equilibrium but neglects derivative of the cross-section inertia function and the flanges' contribution to the vertical equilibrium, • METHOD 4: the enhanced stress recovery proposed by Balduzzi et al. [6], • METHOD 5: the stress analysis performed trough 3D FE that will be used as reference solution in the discussion of numerical results.
It is worth highlighting that METHODS 2, 3, and 4 assume that the beam behaves under the assumption of plane stress state. Despite this assumption could seem highly restrictive, it is widely accepted in engineering practice since the first-design of most of the steel structures can be done analyzing planar frames. Furthermore, METHODS 2 and 3 provide the stress recovery only for bi-symmetric linearly tapered beams subjected to shear-bending loads. Therefore, due to the preliminary nature of the considered investigation and aiming at performing a fair comparison, the present paper considers only beam's geometry and loads that could be tackled by al the considered methods.
The paper is structured as follows.
• Section 2 provides formulas for the cross-sectional stress analysis approaches introduced so far. Specifically, (i) it provides a synopsis of formulas used for prismatic beams, (ii) it briefly resumes the analysis approach proposed by Blodgett [11], (iii) it reports formulas derived by Vu-Quoc and Léger [39], and (iv) it derives closedform formulas for I shaped bi-symmetric tapered beams starting from the stress representation proposed by Balduzzi et al. [6].
• Section 3 reports numerical results obtained by the application of the formulas so far introduced to a simple problem. Specifically, it considers the stress distribution within several cross-sections of a tapered beam subjected to shear bending load.
• Section 4 discusses the numerical results.
• Section 5 reports final remarks and delineates future research.
The simple numerical example considered in Section 3 will highlight that all the engineering tools and procedures for the recovery of stresses from one dimensional models nowadays used both in everyday engineering practice and research suffer from heavy deficiencies and are not reliable. Conversely, the procedure proposed by Balduzzi et al. [6] seems more accurate and promising for the development of new and effective instruments.

Analysis of cross-section stress distribution
Let us consider the homogeneous thin walled tapered beam depicted in Fig. 2(a). The cross-section has two symmetry axes as illustrated in Fig. 2(b) and only the web height h x ( ) w varies along the beam axis as illustrated in Fig. 2(c). Specifically, we assume that it varies linearly with respect to the x coordinate where ′ h represents the slope of the flanges and h 0 is the web height at x = 0.
In the following we are going to assume that (i) the beam behaves under the assumption of small displacements and strain, (ii) the material is linear-elastic and isotropic, (iii) cross-sections behave as a rigid body, and (iv) the beam behaves under the assumption of planar stress state, as usually accepted in engineering practice. As a consequence of assumption (iv), all parameters are constant with respect to the depth coordinate z, simplifying the 3D beam geometry to the planar problem depicted in Fig. 2(c) in which the mechanical properties of layers are proportional to the flange width and the web thickness, respectively. Finally, two bending moments M (0) and M l ( ) and two shear forces V (0) and V l ( ) (satisfying equilibrium equations) act on the beam extremities, as illustrated in Fig. 2(c). For convenience in the notation, we introduce the piece-wise function b x y ( , ) defining the variation of flange width and web thickness with respect to the y coordinate

METHOD 1 Prismatic beam's stress recovery
According to the formulas proposed in standard literature for prismatic beams the horizontal stress distribution within the cross-section can be evaluated as where, for the bi-symmetric cross-section the cross-section inertia J x ( ) reads Furthermore, the shear stress distribution within the cross-section can be evaluated as where the first order of area about the centroid axis Q x y ( , ) V results as Fig. 2. Thin-walled non-prismatic beam, Cartesian coordinate system, dimensions, and load definitions. Assuming that only the flanges resist the bending moment, the magnitude of the force F f resulting from flange's horizontal stress reads and the resulting horizontal stress in the flanges is defined as ( ) x f f f f (9) where, given a positive bending moment the stress is positive in the bottom flange and negative in the top one.
As illustrated in Fig. 3, the horizontal equilibrium on the flanges's boundary induces a non vanishing shear which magnitude is Therefore, the resulting flange's vertical force is defined as leading to the following distribution of the shear stress

METHOD 3 (Vu-Quoc and Léger [39])
The beam model proposed by Vu-Quoc and Léger [39] assumes the same distribution of the horizontal stress of the prismatic beam. Therefore, Eq. (4) is valid also for the herein proposed approach.
Conversely, according to [39,Formula (26)] and consistently with the notation and sign convention introduced at the beginning of Section 2 the shear stress distribution within the web is defined as

METHOD 4 (Balduzzi et al. [6])
The beam model proposed by Balduzzi et al. [6] assumes the same distribution of the horizontal stress of the prismatic beam. Therefore, Eq. (4) is valid also for the herein proposed approach.
Conversely, the recovery of the shear stress distribution within the cross-section is done according to the Jourawsky theory derivation [18]. Considering the horizontal equilibrium of an infinitesimal portion of the beam of size × dx dy and after few calculations the shear stress distribution results to be defined as where C is a constant of integration. Substituting the horizontal stress distribution (4), Eq. (14) becomes Thereafter, the boundary equilibrium (see Fig. 3) and the horizontal equilibrium between flange and web allows to determine the value of the constant of integration. Readers may found further details about the so far resumed procedure in [6, Section 3.3].
Introducing the herein considered beam geometry (see Fig. 2) in Eq. (15), the shear stress distribution within the cross-section is defined as It is worth recalling that both the functions Q x y ( , ) is minimal for y = 0, whereas Q x y ( , ) V is maximal. As a consequence, Formula (16) suggests that, when the shear load is prevailing, the maximal shear will be located at the cross-section geometrical center, as usual for prismatic beams. On the contrary, when the bending moment is prevailing, the maximal shear will be on the web boundary or in the flanges.

METHOD 5 3D FE (reference solution)
In order to provide a reference solution describing the mechanical behavior of the considered body as accurately as possible we modeled the whole beam body with 3D FE, through the commercial software  [34]. Specifically, we exploited the problem symmetry with respect to the plane z = 0 and we modeled only the left hand half of the domain, imposing that displacements with respect to the z direction vanish on the symmetry plane. Furthermore, we assumed that, at the beam ends, shear forces result from a uniform surface vertical shear stress over the web and the bending moment results from a uniform surface traction/compression over the flanges. In order to obtain a numerical solution, we set the following mechanical properties of the material: = E 210 000 MPa and = ν 0.3 and we impose a vertical and two horizontal constraints on two nodes in the initial cross-section. Finally, we adopted a structured mesh of × × 2000 3 23 trilinear elements for the flanges and a structured mesh of × × 2000 160 2 trilinear elements for the web. Globally, 916,000 elements constitute the 3D FE model. It is worth highlighting that a so refined mesh is usually not required in engineering practice, but is herein adopted in order to ensure that errors affecting the reference solutions are negligible. Finally, since the proposed analysis of cross-section stress distributions are done considering a plane stress state and aiming at performing a fair comparison, we further elaborated the 3D FE results calculating the mean value of all the quantities of interest with respect to the z direction.

Numerical results
Aiming at discussing the proposed methods capabilities, this section provides the stress recoveries done with methods introduced in Section 2 for a simple example.
According to the notations introduced in Fig. 2, we assign the following numerical values It is worth noticing that the inclination angle of the flanges is equal to 2.3 deg, significantly smaller than the limit of 15 deg indicated by Marques [23] as the limit for the validity of prismatic stress recovery. Furthermore, the web height is the prevailing dimension in the left hand part of the domain whereas the flanges depth is the prevailing dimension on the right hand part of the domain. Finally, we set the forces acting on the beam ends as follow It would be noted that the applied forces are in equilibrium. Furthermore, the bending moment is positive in the right hand part of the domain, vanishes at = x 7000 mm and becomes negative in the left hand part of the domain. As a consequence, in the left hand part of the domain the beam height decreases according to the bending moment magnitude whereas in the right hand part the beam height and the bending moment have opposite derivatives.
Aiming at discussing the stress analysis tools introduced in Section 2, we consider 5 cross-sections A x ( ) i with = x 1000, 3000, 5000, 7000 i , and 9000 mm, respectively.
It is worth recalling that the maximal Mises's stress is a crucial parameter for practitioners and often used to check the resistance of cross-sections. For the cross-section stress recovery methods introduced in Section 2, the Mises's stress is calculated as whereas the reference solution considers the results coming from the 3D FE software ABAQUS. Obviously, the reference solution takes into account all the components of the 3D stress tensor resulting therefore significantly more refined than the information provided by Eq. (21). Anyway, the comparison of the simplified methods with accurate 3D FE results will highlight their effectiveness in tackling real structures. For each quantity ζ, we consider the relative error = − e ζ ζ ζ ζ ref ref (22) expressed in percentage. Differently from usual error definitions, the absolute-value operators are omitted in Eq. (22). This choice depends on the fact that we would like to highlight when simplified formulas introduced in Section 2 overestimate (i.e., lead to a positive error) or underestimate (i.e., lead to a negative error) the reference values. In the authors' opinion, this information is crucial since it allows to determine if the predictions are on the safe side or not. Table 2 reports the values of shear stress evaluated in the middle and on the boundary of the web. In fact, as discussed in Section 1 and due to the cross-section symmetries, the location of the maximal shear in tapered beams can be both at the middle or on the boundary of the web. Table 3 reports the values of the Mises's stress evaluated in the middle and on the boundary of the web and at the highest/lowest point of the cross-section. It is worth recalling that, in analogy with prismatic beams, the chosen points are the ones that could be critical.
Both Tables 2, 3 indicate if the cross-section stress-analysis methods predict the location of the maximal stress correctly (✓ right, wrong). Finally, bold text in Table 3 highlights the location of the maximal Mises's stress in the reference solution, the estimations of this value provided by the different methods, and their relative errors.

Discussion of the results
We start by noticing that, among the five considered cross-sections, . This simple observation confirms what is stated by Kim et al. [19] i.e., the location of maximal horizontal stress does not coincide with the location of maximal bending moment.
In greater detail, for prismatic beams the cross-section inertia I x ( ) is a constant parameter. Therefore the search for maximal + σ x h h ( , /2 ) x w f and for maximal M x ( ) are equivalent (see Eq. (4)). Conversely, for tapered beams, finding the maximal requires to solve the following equation Obviously, the solution of Eq. (23) is non trivial and could require specific numerical tools, usually not needed for prismatic beams.
A few simple analytical calculations would allow to generalize the so far introduced considerations to the maximal shear and the maximal Mises's stress but this does not seem to be necessary to make the point clear.

Horizontal stress
Formula (4) (METHOD 1 prismatic beam's stress recovery, METHOD 3 (Vu-Quoc and Léger [39]), and METHOD 4 Balduzzi et al. [6]) results to be accurate in most of the cross-sections. Specifically, a significant difference to the reference solution can only be noticed in A x ( ) 4 ( Fig. 7(a)). Nevertheless, in this cross-section the magnitude of horizontal stress is negligible if compared to the shear stress magnitude in the same cross-section ( Fig. 7(b)) and also when compared to the horizontal stresses in other cross-sections ( Fig. 4(a), Fig. 5(a), Fig. 6(a), and Fig. 8(a)). In particular, it is worth noticing that the difference so far highlighted depends on the fact that the proposed methods can   tackle only first order effects, whereas the 3D FE model has the capability to catch higher order effects that result to be significant in this cross-section. Formula (9) (METHOD 2 Blodgett [11]) results to be more coarse than Formula (4). In particular, it usually overestimates the maximal horizontal stress in flanges with errors up to 20% (Fig. 4(a)), but, unfortunately, errors are not always on the safe side ( Fig. 8(a)). Furthermore, it is not able to evaluate horizontal stress in the web, leading therefore to an extremely coarse description of the stress distribution within the web.

Shear stress
Formulas (6) (METHOD 1 prismatic beam's stress recovery) and (13) (METHOD 3 Vu-Quoc and Léger [39]) are not reliable in evaluating the real shear stress distribution, as illustrated by Fig. 4(b), Fig. 5(b), Fig. 6(b), and Fig. 8(b). In detail, the only case in which they lead to accurate results is the cross-section subjected to a pure shear load ( Fig. 7(b)). Furthermore, as reported in Table 2, METHOD 1 (prismatic beam's stress recovery) leads to errors ranging from − 45% to + 115% whereas METHOD 3 (Vu-Quoc and Léger [39]) leads to errors ranging from − 70% to + 190%. Last but not least, both Formulas (6) and (13) are  also not able to individuate the right position of the maximal shear (see Table 2 A x ( ) 1 and A x ( ) 2 ). In fact, since the bending moment is the prevailing internal force in A x ( ) 1 and A x ( ) 2 , the maximal shear is not located at the middle of the web (see Fig. 4(b) and Fig. 5(b)), but Formula (6) completely ignores this important effect. Finally, it is worth noticing that the prismatic beam formula always predicts quasi-vanishing shear in the flanges, whereas the reference solution indicates that it has generally a non-negligible magnitude.
Formula (12) (METHOD 2 Blodgett [11]) allows to provide a reasonably accurate prediction of the shear stress mean value in flanges and is generally more accurate than both Formulas (6) (METHOD 1 prismatic beam's stress recovery)) and (13) (METHOD 3 Vu-Quoc and Léger [39]) since it leads to errors ranging from − 25% to 20%. Nevertheless it is not able to individuate the maximal shear location and its estimations do not guarantee a conservative estimation of the quantities of interest for practitioners. Formula (16) (METHOD 4 Balduzzi et al. [6]) results to be the most effective and accurate since it is able to individuate correctly both the maximal shear and its location, leading to errors smaller than 1% in all the cross-sections.

Mises's stress
Looking at Table 3, it is worth noticing that (i) METHOD 1 (prismatic beam's stress recovery) leads to errors ranging from − 45% to 115%, (ii) METHOD 2 (Blodgett [11]) leads to errors ranging from − 90% to 20%, (iii) METHOD 3 (Vu-Quoc and Léger [39]) leads to errors ranging from − 70% to 190%, and (iv) METHOD 4 (Balduzzi et al. [6]) leads to errors always smaller than 5%. Furthermore, it is also possible to notice that using both METHOD 1 (prismatic beam's stress recovery) and METHOD 3 (Vu-Quoc and Léger [39]) the more consistent errors occur at y = 0, whereas using METHOD 2 (Blodgett [11]) the more consistent errors occurs at the intersection between web and flanges. Finally, METHODS 1, 2, and 3 sometime are not able to individuate the right position of the maximal Mises's stress.
Since the maximal Mises's stress defines the cross-section resistance and the usage of METHODS 1, 2, and 3 could lead to underestimate this fundamental parameter more than 50%, we can conclude that they are unreliable. Furthermore, since they are not able to individuate the right location of the maximal Mises stress -i.e., the location of the failure point-they lead designers to uselessly increase the size of certain beam components with any immediate benefit for both the cross-section and the structure strength. As a consequence, literature methods turn out to be not only ineffective but also extremely dangerous.
Conversely, METHOD 4 (Balduzzi et al. [6]) results to be the most effective and accurate since it is able to individuate correctly both the maximal Mises stress and its location, leading to errors smaller than 5% in all the cross-sections.

Final remarks
This paper tackles the problem of the recovery of stress distribution within the cross-section of bi-symmetric tapered steel beams. Analytical results (available in the literature since the first half of the past century) highlight that the continuous variation of the beam height induces shear stress distribution substantially different from the one occurring in prismatic beams.
In order to predict the stress distribution within tapered beams, a certain number of methods were proposed in literature and several of them were compared within this paper trough the evaluation of their performances in a simple practical example. Specifically the paper considers • METHOD 1 stress recovery usually adopted for prismatic beams, • METHOD 2 stress recovery proposed by Blodgett [11], • METHOD 3 stress recovery proposed by Vu-Quoc and Léger [39], • METHOD 4 stress recovery proposed by Balduzzi et al. [6]. METHOD 1 can underestimate the maximal shear up to 45% and lead to underestimation up to 25% for the Mises's stress. Similarly, METHOD 2 can underestimate the maximal shear up to 25% and the Mises's stress up to 55%. METHOD 3 can underestimate the maximal shear up to 70% and the Mises's stress up to 10%. As a consequence, METHODS 1, 2, and 3 are not reliable and should be avoided in design practice. Conversely, METHOD 4 leads to errors always smaller than 5% in all the considered cross-sections, representing therefore an effective and accurate tool for the cross-section stress analysis.
Further developments of this work will include the generalization of the METHOD 4 to 3D thin-walled beams and the investigation of the effects that non-trivial stresses have on beam stiffness, buckling, and dynamic behavior. Table 3 Evaluation of Mises's stress in critical points, correctness of the maximal Mises's stress position (✓ right, wrong), and relative errors. Bold text highlight the maximal Mises's stress within the cross-section in reference solution. MET.1 (prismatic beam's stress recovery), MET.2 (Blodgett [11]), MET.3 (Vu-Quoc and Léger [39]), MET.4 (Balduzzi et al. [6]), and MET.5 3D FE (reference solution).