Software implementation of section class and resistance calculation for general loading case *

In this paper, the problems arising in implementation of the cross-section class and resistance according to Eurocode, especially for classes 1 and 2 of the steel sections, are presented for general loading case. As Eurocode assumes full plastification of the section, regardless of corresponding strain in the material, it is inevitable to find the position of the plastic neutral axis for ultimate limit state of the section. But, for this purpose, one cannot use the Eurocode’s expressions for section resistance. Moreover, solution and strategies used in the steel design module of the Matrix 3D are presented in the paper.


Introduction
Trends in modern technical regulations, to which Eurocode belongs, are certainly on the side of almost complete automation of the calculations in the appropriate software. Unfortunately, many of the commercial software don't follow the rules and logic of regulations to their full extent, either due to the complexity of the necessary algorithms or their cumbersome change, frequent changes in regulations, or often insufficient analysis and knowledge of the regulations themselves. In Serbia, in addition to one widespread commercial software, there have been a couple of attempts to develop smaller academic software to the level of comprehensiveness of complete structural calculation with design, but no such project has been completed. At the Faculty of Technical Sciences, the academic software Matrix 3D [5] [6] [7] is being developed, which implements a module for dimensioning steel elements according to Eurocode. The initial problem that arises in this kind of code is described here in more detail and a solution is proposed, mostly for the purpose of encouraging and assisting other similar endeavors. Another reason is that the literature related to this particular problem is scarce, mostly because it is considered as an elementary problem. Furthermore, Eurocode does not give its equations for the resistance of the "U" cross-section, and the assumptions made by [2] [3] papers, regarding the position of the neutral axis, have been proven wrong in this paper. The strategy used in Matrix 3D to efficiently calculate the cross section class but also the resistant surface of the cross section is presented in following chapters, along with the assumptions made during the solution process. It is well documented [2] [3] [4] that the Eurocode expressions for cross section resistance are not very precise, hence there is a justification to use as simple method as possible. Surely, a * Awarded as best young researcher's paper, presented at ASES 2020 Symposium, Arandjelovac, Serbia, May 2021. * Corresponding author: E-mail address: nozica.tanja@uns.ac.rs more precise solutions are possible, but for the price of more time-consuming computational effort.

Problem definition
During the design of each steel element, according to Eurocode, it is first necessary to determine the class of the cross-section, and then the resistance of the cross-section. The class of the cross section depends on the stress distribution in each individual part of the cross-section. For classes 1 and 2, it is necessary to find the position of the plastic neutral axis for each load combination. To do this, it is necessary to know the limit strain state in order to extrapolate strains due to the design values of internal forces. Considering that Eurocode allows for the total plastification of cross-section it is impossible to know the strain state. It is illogical because it is implicitly allowed for strain to reach infinity. This is why the stress state must be found, without firm theoretical background. In other words, it is necessary to know the actual cross-section plastic resistance to determine the position of the plastic neutral axis.
2.1 Resistance of cross-section due to bending and axial force according to en 1993 In order to determine the plastic limit state of the section due to biaxial bending and axial force, Eurocode 3 provides expressions in EN 1993-1-1 [1]. The parameters α and β are included in this approximate expression, which represent the influence of the type of cross section (U, I, L etc.). Based on the expression, an interaction surface that represents the cross-sections plastic resistance, can be created ( Figure 1).

Figure 1. Biaxial bending and axial force interaction surface (EC3) for IPE300 cross-section at axial load levels n=NEd/NRd
The surface defined by Eurocode does not match the real surface because it is a simplified one. Another problem is that there are only expressions for "I" sections and hollow sections, but not for "U" sections.
Anyhow, the internal forces in the section for a given load combination can assume one of the three different position in the force space: 1) inside the resistance surface, 2) on the resistance surface or 3) outside the surface. The case when the point on N-Mz-My diagram is outside resistance surface means that this forces cannot be transmitted by the given section, hence there is no need to continue with the classification of the section. The case when it is exactly on the surface is practically almost impossible, but it is trivial to determine the plastic neutral axis in this case. In most of the cases, the point (NEd, MyEd, MzEd) is somewhere inside the surface. In this case, it is necessary to reach the ultimate state by scaling some (or all, or some ratio) of the internal forces. It is not specified which ratio, nor in Eurocode or similar documents. There are several different strategies, and some of the software, such are Dlubal RFEM [8] or SOFISTIK [9] explain their strategies, while majority of others do not mention the problem in their theory reference.

Neutral axis position for biaxial bending with axial force
The papers [2] [3] show the cross-section resistances due to biaxial bending with the axial force. However, they are calculated on the assumption that in the case of biaxial bending with axial force neutral axis preserve its incline independent of the intensity of the axial force.
In general, when a cross-section which is subjected to a bending moment which does not coincide with the cross-sectional main axes, the neutral axis does not coincide with the angle of the bending axis nor with one of the main axes. Moreover, the neutral axis rotates from the bending axis towards the weaker axis with an increase of the axial force. Hence, when the axial force increases retaining the same bending moment, the plastic neutral axis gains a certain movement, but also the shift of the angle towards the weaker axis of cross section.

Method
In order to design the cross-section, it is necessary to define the conditions when the cross section is in the ultimate limit state. The conditions can be set according to the ultimate strain (as it is case for reinforced concrete sections), but the problem with the steel sections defined as they are in EN 1993-1-1 [1] is that the ultimate strain for plastic analysis is not limited. That is why it is necessary to define the ultimate limit state of the cross section using static equilibrium conditions. For creating interaction surfaces, variables n, mz and my, are needed. n represents the ratio of the applied axial load and axial capacity and my and mz represent the ratios of the applied bending moments with the corresponding plastic moment capacities. In order to determine the ultimate limit state, it is necessary to keep one value, the axial force or the ratio between bending moments My and Mz fixed, and increase the other until the ultimate limit state is reached. Considering this, there are two different strategies to determine the current cross-section state. The first one is to keep the axial force fixed (point 1 in Figure 2), while the second one is to preserve the ratio between axial force and bending moment vectors (point 2 in Figure 2).

Section resistance calculation in Matrix3D
In this research, the interaction surface is calculated by keeping the axial force fixed. This means that for given axial force the horizontal section of the resistance surface is calculated, but due to the symmetry of the standardized steel cross sections, only a quarter or the half of this section needs to be found. Still, it is impossible to find the correct point on this section in the first iteration, because the angle of the neutral axis is not known apriori. There is a possibility to search for the point 1 in the Figure 2 by some sort of numerical strategy, such as halfing method, or to calculate the whole curve (section of the surface) rigorously, and then to find the point that corresponds to the desired ratio between bending moments around My and Mz. Since this calculation should be done for each loading combination, and each member, it is desirable to use simplified representation of the section in order to make the computational process less time-consuming. Hence, the usual steel cross sections are represented with the straight lines and the transition zones between web and flanges are disregarded. This does not introduce significant error since the main goal is only to determine the class of the section.
To determine the actual position of the plastic neutral axis for certain internal forces, it is chosen to firstly determine the actual interaction diagram for the cross-section subjected to biaxial bending and given intensity of the axial force.
As mentioned, the simplified line representation of crosssections was chosen to simplify the determination of the border cases of the position of the plastic neutral axis. Furthermore, the more precise method approaches the fiber beam element [4], and fiber beam-column elements are already implemented in Matrix 3D, and used for materially and geometrically nonlinear analysis [7]. Therefore, the problem consists of calculating two unknowns, the angle of the neutral axis (α) and its position along the z axis (Δ) using three static equilibrium conditions.
The algorithm for determining an actual interaction surface in MATLAB is described through the following steps: 1. The first control is to ensure that the cross-section area is sufficient to transmit the axial force, that is, is the condition NEd / fy < A met, if not, the calculation is abandoned.
2. If the condition is met, a control whether the surface of the web itself is sufficient to transmit the axial force or whether it is necessary to include the surface of the flange is performed, depending on that, the expressions for determining the limit values (α1 α2 α3) differ as do their corresponding limit positions of the plastic neutral axis along the z axis (Δ).

3.
A loop is used to variate the angle α from 0° till 90°, with an appropriate step, in an apt segment. 4. For a single value of α, its corresponding Δ is calculated analytically from the expression obtained from the condition that the difference of the tension and compression areas is sufficient to transmit the axial force. 5. For one value of α and its corresponding Δ bending moments around the main axes are calculated. These moments correlate to one NEd, one α, and its corresponding Δ, and represent one point on the interaction surface.
In the next part the limit values (α1 α2 α3) on the example of a "I" section and rectangular hollow section are explained.
The area required to transmit the axial force for "I" section NEd/fy can be found in the following intervals: NEd/fy < A, Aweb < NEd/fy < Aweb+Aflange and Aweb+Aflange < NEd/fy < A When calculating the whole interaction curve for a given normal force, it is necessary to algorithmically cover all situations that occur for all angles of the plastic neutral axis, and in parallel, to optimize the calculation time. Therefore, the first step in the algorithm is to divide the cases into those when the normal force can be transmitted by the web, when the web is not enough to transmit the normal force and when the sum of the flange and the web is not enough to transmit the normal force. This is all necessary in order to analytically write formulas for determining Δ, for different ranges of α. So, while α varies according to the adopted step, Δ is calculated analytically. This is possible because the stress diagram is very simplea block diagram. In concrete cross-sections, this is impossible to perform analytically, the position of the neutral axis, Δ, is found iteratively in an inner loop, as described in [10].
For a situation where the web can transmit a normal force in the cross-section, there are three limit values for the angle of the plastic neutral axis as shown in Figure 3. For angles less than α1, it makes no sense to calculate Δ or divide that interval, because the whole interval has the same values and they represent a segment of the curve on the My axis.
When the area required to transmit the axial force is greater than the area of the web but smaller than the sum of the areas of the web and flange, that is, when the condition Aweb<NEd/fy<Aweb+Aflange is satisfied the plastic neutral axis for the value of the angle α=0 lies in the flange. In this case the angle α variates between the angle α3 and 90°, where the angle α3 represents the angle of the plastic neutral axis when the plastic neutral axis intersects the lower flange. In this case as appose to the previous one the plastic neutral axis for the angle of 90° exits the area of the web (Figure 4).  In the last case, the sum of the areas of the flange and web is not sufficient to transmit the axial force. Then, just like in the last case, the plastic neutral axis exits the area of the web but the difference being that the plastic neutral axis is located in the upper half of the flange for the value of the angle α=0. By variating the value of the angle α from α3 till 90° the bending moments resistances MRdy and MRdz are calculated ( Figure 4).
Also by a similar procedure to the previous example, the limit values (α1 α2 α3) and limit areas for rectangular hollow section are determined and shown in Figure 5. With the difference being that in the case of a hollow cross-section, the third case, when the sum of areas of the web and two flanges is not sufficient to transmit the axial force, is not defined, because in that case, the plastic neutral axis varies in the flange for small values of the angle, and for larger values of the angle it varies in the web of the cross-section.

Results
The interaction curves for the cross-section IPE 300 made from steel S355 are presented in Fig. 6 and 7. By variating different levels of the axial load, a series of parallel sections of the interaction surface is calculated.
In Figure 6 normalized axial force n=0.4 corresponds to the case when the web area is sufficient to transmit the axial force ( Figure 3). The value for n=0.6 correspond to the case when it is needed more than the area of the web to transmit the axial force but less than the sum of the areas of the web and flange (the first and second case in Figure 4). Lastly the value n=0.8 represent the case when is needed more than the sum of the areas of the web and flange to transmit the axial force (the third and fourth cases in Figure 4). Hence, each of these 3 curves was calculated according to a different case shown in Figures 3 and 4.
The Figure 6 also shows curves calculated by an expression given in EN 1993 [1], for the analyzed section IPE 300. As can be seen from the Figure 5, Eurocodes interaction curve matches the curve obtained by the presented procedure quite well, with some small deviations. The values of bending moment resistances around the stronger axis coincide with the moment resistances calculated by an expression given in Eurocode, while the differences are observed in the bending moment resistances around the weaker axis. The deviations are consequence of the approximate representation of the cross-section with straight lines.
Obtained by the same procedure as the I section, the interaction curves for hollow section RHS 100x100x4 are shown in Figure 7.
The problem with the resistances calculated with the expressions obtained from Eurocode is that the position of the plastic neutral axis cannot be determined. Therefore, the plastic neutral axis can only be determined for an interaction surface that is obtained by the same method that is used to determine the position of the plastic neutral axis. Figure 8 presents different values of normalized bending resistances (my and mz) for different values of normalized axial force (n) and with the same inclination of the plastic neutral axis (α=84°). It is evident that this red curve in Fig. 7 in not a meridian of the resistance surface. This means that by increasing the axial load, the inclination of plastic neutral axis will not be the same for the same direction of applied bending moment.

Discussion and conclusion
The majority of challenges of implementation of the steel section design in software are revisited in this papers. Some of the more involving issues are out of the scope of the paper, such are thin-walled sections and types of buckling there are prone to, or sections loaded by bimoment. Still, even for usual rolled sections, Eurocode tacitly implies that the ultimate resistance for each loading combination should be found in order to determine the class of the cross-section. This can easily become overly time-consuming task, or inadequately solved. It is shown in the paper that for different levels of axial force, the inclination of the plastic neutral axis could not be known a priori. Hence, the algorithm for I and RHS//SHS sections is presented. Further, the approximate resistances of section resistance given in EN 1993-1-1 [1] are compared to ones obtained by this simple algorithm.