Optimal preliminary design of variable section beams criterion

The present paper discusses about optimal shape solution for a non-prismatic planar beam. The proposed model is based on the standard Timoshenko kinematics hypothesis (i.e., planar cross-section remains planar in consequence of a deformation, but it is able to rotate with respect to the beam center-line). The analytical solution for this type of beam is thus used to obtain deformations and stresses of the beam, under different constraints, when load is assumed as the sum of a generic external variable vertical one and the self-weight. The solution is obtained by numerical integration of the beam equation and constraints are posed both on deflection and maximum stress under the hypothesis of an ideal material. The section variability is, thus, described assuming a rectangular cross section with constant base and variable height which can be described in general with a trigonometric series. Other types of empty functions could also be analyzed in order to find the best strategy to get the optimal solution. Optimization is thus performed by minimizing the beam volume considering the effects of non-prismatic geometry on the beam behavior. Finally, several analytical and numerical solutions are compared with results existing in literature, evaluating the solutions’ sensibility to some key parameters like beam span, material density, maximum allowable stress and load distribution. In conclusion, the study finds a critical threshold in terms of emptying function beyond which it is not possible to neglect the arch effect and the curvature of the actual axis for every different case study described in this work. In order to achieve this goal, the relevance of beam span, emptying function level and maximum allowable stress are investigated.


Introduction
The adoption of variable section beams in the Civil Engineering field is really widespread. One can refer e.g. to the bridge field where it is useful to adopt a non-prismatic beam in an optical of optimization of the use of material. As a matter of fact, it is convenient to increase the cross section in the area where the bending moment demand is greater [8]. The non-prismatic beams were always adopted in the Architectural field, not only as a result of an optimization process, but rather in order to find aesthetic requirements and establish a strong relationship among the form, the structural requirements and also looking to the harmony with the place in which the structure is located. In [11], the author studied and critically analyzed the historical aspects of the different type of optimization procedures in form findings for variable section beams but, nowadays, very complex geometry can be easily treated in a parametric way with valuable technologies and powerful software such as Grasshopper as showed in [15]. The adoption of variable cross section beam involves also the knowledge of the mechanical behaviour of all available structural materials, e.g. in many applications tapered beam are mainly adopted and if they are for example in timber one can refer for instance to [9].
In the classical study of the beam, it is assumed to be valid the Navier hypothesis which states that the sections remain planar also in the deformed configuration. The classical Euler-Bernoulli beam is based on the hypothesis to neglect the shear deformation contribution whereas the Timoshenko beam takes into account it, resulting in a more complete theory [5,14]. In literature, many articles analyze the problem of variable cross section beam. For example in [2,3,14], authors analyze the general compatibility and equilibrium of non-prismatic beams proposing an approximated Timoshenko-like model solving a system of coupled ODE. Many other approaches were developed in literature, e.g. to treat also curved beams one can refer for instance to [6] for the isogeometric analysis formulated by Hughes to analyze with FEM curved beam based on the adoption of B-spline interpolation. A simpler, and usually adopted approach, is the finite differences methods which allow to take into account non-homogeneous and non-prismatic beam approximating the derivative of the functions on a discretized mesh of the structure such as in [1,16].
In this preliminary study, the main idea is to optimize the lateral shape of a prismatic beam composed of an ideal material which assumes a non-prismatic behavior because of the intervention of an emptying function which defines a new height profile. Ideally neglecting the conventional arch beam behaviour, in this study a sort of "dog-bone" beam shape is analyzed in which the barycentric beam axis remains always straight. Another study conducted on a variable section beam with the aforementioned geometry is also presented in [7] where a solution with a power series method is presented. In particular, the objectives of this study are both to optimize the shape minimizing the self-weight and, at the same time, to identify the range of validity of this simplified approach neglecting the conventional arch effect. For this reason, a final comparison is made between the dog-bone beams deployed in Matlab with the optimal real arch geometry beam modelled with FEM using the professional software Midas Gen. Substantially, the main aim of this last comparisons is to demonstrate that for different emptying values, regardless the span length, the arch effect produced by the curvature of the actual axis in the optimal shape profile it is negligible and it is also possible to adopt this simplified dog-bone approach for moderate relative emptying values.
The meta-heuristic approaches represent nowadays very interesting methods in structural optimization fields, besides the advantage of no prior requirements on differentiability conditions of the objective function (OF) and the constraints [13]. In particular, the genetic algorithms (GA) belong to the branch of evolutionary algorithms which are based on Darwin's theory of survival of the fittest. This algorithm was originally proposed by John Holland in 1975 at the University of Michigan, and it was extended later mainly by his student David E. Goldberg [10]. It takes inspiration from a crude imitation of what happens in Nature in which species whose individuals are better evolve to the next generation (individuals with the better fitness in term of OF) and the core is based directly on the genetic mechanism of transmission of the parents' features to offspring chromosomes adopting specific mathematical operators: crossover operator, mutation operator and selection operator. Due to the GA random nature, unfortunately, there are no mathematical proofs of its convergence, but numerically studies demonstrated that they are able to succeed also dealing with highly nonlinear, non-convex and discontinuous domains. Furthermore, since the GA require only the OF evaluation and not require any information about the gradient, its computational effort will be lower compared with mathematical programming techniques such as NLP methods [12]. For further readings about the GA and in general to metaheuristic algorithms, one can refer to e.g. [10].
In Sect. 2 the analytical model for the non-prismatic beam is exposed also with an insight in a dimensionless form of the differential problem. Subsequently, in the Sect. 3 the numerical adopted approach is discussed. Moreover, in Sect. 4, a comparison between different emptying imposed values and a comparison between optimal Matlab solutions with the FEM validation models are discussed. Finally, in at the end of the Sect. 4 and then in Conclusions, a comparison between some features for different emptying and span length imposed value are investigated and discussed in order to understand if the dog-bone beam model exhibits a good behaviour in comparison with respect to the reference FEM model and if this solution can be assumed as a preliminary criteria of preliminary design.

Analytical beam model
In order to solve the problem of variable section beam one has to refer to the equation of elastic line for non-prismatic beam as stated in [8], which it can be written in general as the following form: where x is the longitudinal coordinate of the beam axis, y is the beam deflection, E is the elastic modulus, J(x) represents the moment of inertia variable along the x coordinate and q(x) is the distributed load which comprises both Research Article self-weight and live load applied on the non-prismatic beam element. The Eq. (1), once written in abbreviated form and suppressing the dependencies from x, simplifies as: By developing the second order derivative, we obtain: Finally, deriving one more time and rearranging the equation's terms from the maximum order derivative to the minimum one, it yields to: where, under the assumption to adopt a rectangular cross section shape with constant base b and height h(x) variable along the x coordinate, the inertia terms becomes: are respectively the Moment of Inertia and its high order derivative, in which the first and second order derivatives of the cross sectional height, h � (x) and h �� (x) , appears. It has now introduced the beam section depth function (x) , which varies with depth beam section where h 0 represents the initial beam section depth. At this point, starting from a prismatic beam, one have to consider now to apply an emptying function from the bottom side which define a new height profile but maintaining h 0 at the fixed joints. Therefore (x) identifies a emptying function which is the sum of many harmonic as the following one: in which the A 0 term denotes the max amplitude of the emptying function and L denote the beam span lenght. Therefore, the (x) function can be written as the sum of N harmonics: The amplitude coefficients A i directly represent the maximum emptying value for each term of the sum when the trigonometric function is equal to unity for a certain value of x and, therefore, it is also possible to denote them as h i . The following relationships between the emptying function (x) and the beam section depth function h(x) are derived: Considering the new equation of the beam section depth function and its higher order derivative, it is possible to obtain Substituting the (5), (6), (7), (11), (12), (13) in the (4) it is possible to obtain the following forth order differential problem: Focusing now on the emptying function which follow a generic trigonometric series as proposed in (9), in order to obtain the explicit formulas for the inertia and its derivative, it becomes Research Article SN Applied Sciences (2021) 3:745 | https://doi.org/10.1007/s42452-021-04702-5 The term q(x) take into account two contributions: the dead load per unit length which is the product of the material density and the cross section of the beam S(x) and the variable load q 0 assumed to be constant with the length which represents the applied live load. This fact leads the forth order differential Eq. (4) to assume the following form In order to numerically solve the forth order differential equation, it is necessary to set the equation as a system of four equation of the first order introducing the vector z Considering the first derivative of z � = f (z, x) , the vector f can be written as below In the following section it is exposed another method to solve numerically the fourth order equation as suggested in [4].
In the hypothesis of a beam section fully defined by the Equation J(x) = 2 S 2 (x) and deriving two times the ratios J � ∕J and J �� ∕J , with the aim of parameterize the system problem,as we show follow: Finally, it is possible to obtain which can be solved for the only parameter h(x), making explicit the parameter J(x):

Optimal design criterion
The mathematical model exposed in the previous section has been implemented in a Matlab script. In order to solve the optimization problem the algorithm adopted is the Genetic Algorithm (GA) with the native function of Matlab. The analysis were conducted on a single span beam x ∈ [0, L] with fixed ends which is characterized with a rectangular cross section with constant base b and variable height h ∈ (0, h 0 ) along the element. Formalizing the optimization problem, the design vector is represented by the maximum amplitude of the emptying trigonometric functions i.e. the design parameter is only one if one considers to adopt only a sine function as the total emptying function. As a matter of fact, for the imposed boundary condition on this beam, the two cases that are of a valuable technical interest for the form finding are the one sine emptying function ( N = 1 lobe) and the two sine emptying function ( N = 3 lobes), as also depicted in Fig. 1.
In order to reduce the material consumption and then the related costs, the OF to be minimized is equal to the total weight of the beam W tot = ⋅ V . Dropping the material weight density which is assumed to be constant, the OF only becomes the total volume V which, for a general trigonometric series, can be calculated as in the following The optimization problem require to solve the structural analysis performed with the presented analytical model and the search space is further reduced due to the presence of structural constraints to be satisfied in terms of maximum allowable deflection in midspan, conventionally assumed as lim = L∕250 , and in terms of maximum allowable stress at the extreme fibers of each cross section evaluated with the Von Mises ideal stress criterion The id is assumed to be the yielding stress for an ideal material which exhibit the same behaviour both in tension and in compression. Setting z as the vertical local axis for each rectangular cross section, the elastic normal tension is acting with its maximum intensity at the outer fibers of the cross section accordingly to the Navier's formula, (23) whereas the tangential stress reach the maximum intensity in the center of each cross section accordingly the Jourawsky's parabolic profile, The analytical model exposed in the previous section is numerically solved with Matlab adopting the internal solver bvp4c. As suggested by [4], the general fourth order Eq. in (4) can be considered equivalent to the following system of four first order ordinary differential equations which can be rewritten in vectorial notation as s e t t i n g y 1 = y(x), y 2 = (x), y 3 = M(x) and y 4 = V (x) , where is the rotation of the cross section. In this way, the variability of the cross section is taken implicitly into account only by J(x) and, due to the fact that all the equation are coupled, this is taken into account in the entire system avoiding to explicitly solve the fourth order equation depending by the derivative of the inertia. Moreover, in this way the solutions of the system directly represent shear, moment, rotation and deflection curves.
At the end of this section, it is necessary to stress that in this preliminary study, all the equations before exposed are valid for a variable section beam with a straight axis. In reality, the beam model showed in Fig. 1 has a curvilinear axis which involved also the arch effect and it is necessary also to take into account the axial force produced by the pre-existing curvature of the system. Therefore, in this preliminary study, the results are referred to a ideally straight axis beam with variable cross section similarly to a dog-bone shape beam, as shown in Fig. 2.
In the next section, some numerical examples and the respectively main results are presented and some comparisons are discussed on the validity of approximating the real system to a sort of dog-bone variable section beam.

Numerical examples
In order to make some first comparisons, an ideal beam is analyzed with different emptying ratios and making a comparison with a full-beam without any emptying ratio. The analysis are conducted assuming a fixed ends cross section with b = 0.50 m and h = 1.5 m and taking into account the material parameters of an ideal concrete C20/25 ( E = 29962 MPa ). The span length is set to L = 20 m and three different emptying ratio are chosen for analyzing the one lobe emptying function with 1∕4, 1∕2 and 3∕4 of h as depicted in Fig. 3. Taking into account a uniformly distributed live load of q = 10 kN∕m and both considering the self-weight assuming = 25 kN∕m 3 , the results in terms of deflections, rotations, bending moments and shear diagrams are shown in Fig. 4.
As it is possible to notice, increasing the emptying ratio the beam becomes more compliant and the deformation parameters such as deflection and rotation increase. On the contrary, the shear from a linear function along the beam becomes curvilinear accordingly to the variation of the self-weight along the beam and the vertical reactions decreases because of the diminishing of the self-weight. The bending moment follows substantially the flexural rigidity along the beam. In fact, from the full-beam case where one can find the usual values of qL 2 ∕12 and qL 2 ∕24 at the fixed end and in midspan respectively, the bending moment tends to get a higher decrease in the midspan and a slight decrease at the fixed ends. This justify the aim of the optimization process to try to find a better material distribution along the beam element in order to reduce the actions but paying attention also to the deformability requirements.
Another comparison it is made with a two lobes emptying function. This time it is necessary to choose two values of the amplitude of the sine functions which are combined together, paying particular attention to not choose a combination which violates the physical domain in any point ( h ∈ [0, 1.5] m ). The amplitude for the first sine are set as before, whilst for the second sine function they are chosen as 1∕8, 3∕16 and 1∕4 of h and the resulting beam profiles are shown in Fig. 5. As before, the value of the actions decrease with the increasing of the emptying ratio and the bending moment this time looks like more regular because the decrease in midspan is lower than before. In fact, the two lobes emptying function try to remove material where the bending moment diagram of a full-beam approaches to zero. This results in a bending moment distribution which approximately reminds the bending moment diagram of a full fixed beam. Focusing on the compliance, this time the beam seems more compliant in particular there are two peaks in the rotation diagram and the deflection diagram is less sharp than before, as depicted in Fig. 6. As before mentioned, the purpose of this work is to investigate not only the ideal behavior of a dog-bone beam but especially appreciate what are the parameters or response ranges above which these particular modeling techniques deviate from the desired arch beam ones.
To achieve this goal, the behavior of a dog-bone beam for three different emptying functions, obtained with GA optimization processes, and three different imposed length values is compared with an arch beam FEM model as depicted in Tables 2 and 3.
Different intermediate emptying values are chosen because they depend on the slenderness of the beam. As we expected, the maximum emptying ratios for beams with length of 20 and 30 m reach respectively 1 and 0.76 m whilst the same maximum emptying ratio is achieved to 0.49 m for 40 m beam span as described in Table 1.  should have been defined based on the minimum emptying value obtained for each case studies. This strategy would lead to discarding larger emptying percentages obtained by slenderness beam. In this work, the bending moment at fixed joints and midspan and deflection are the only parameters used to realize these comparisons because they represent the features more sensitive to the different levels of emptying ratio.
It is worth noting that the deviation between two solutions, for both the emptying function studied, is negligible when a constant rectangular cross section is defined, or in other words, when the emptying ratio is zero. For increasing emptying values imposed, the deviation between a reference FEM modeling, which takes into account the arch effect, and the dog-bone Matlab model with a straight axis grows. In the one lobe emptying function case study, the Matlab solution bending moment and deflection tends to be overestimated compared to the more accurate FEM solution. If the maximum emptying ratio would be achieved in a dog-bone beam with 20 m of length, the fixed joint bending moment would have been larger than 50% respect the same one obtained with FEM modeling, whereas the midspan bending moment would have even grown to 2 (as one can see in Table 2).
Generally, the beam with 40 m length tends not to be affected by the model approximation neither even when investigated emptying levels reach approximately 60% of the total height. In this case, the error with respect to the reference FEM solution is about 20%. The one lobe emptying function scenario is confirmed by observing the optimization results for the case of three lobes emptying function. Referring now to the bending moment at fixed joint, the obtained trends are similar to the case of the one lobe emptying function with the exception of the beam with 20 m length. This latter case shows a deviation from the correct solution at fixed joints of about 10% whilst it is negligible at mid span. Although in the one lobe emptying function case the 20 m beam exhibits the worst deviation with respect to the FEM results, in the three lobes emptying function case it, now, represents the wellbehaved solution as illustrated in Fig. 9. As a matter of fact, it reaches the greatest emptying ratio but also exhibits the minimum error with respect to the reference FEM solution.
Moreover, for case study of a beam with a 20 m span length with three lobes emptying function, both displacement and even more for the bending moment at midspan result underestimated by the dog-bone model. In this case, as depicted in Fig. 8, the deviation from the desired arch beam is due to static considerations about normal stresses inclination and emptying function shape: bending moment induced by normal stresses tend to increase the total bending moment at midspan rather than decrease it as in the case of one lobe emptying function. The dog-bone beam model, which represents a variable cross section beam with a straight axis, does not take into account this contribution and this explains its underestimated features. Considering the 40 m span length case with three lobes, the obtained optimal result shows negligible deviations from the more accurate FEM solution regardless of the emptying ratio due to the little optimal emptying obtained value, as illustrated in Fig. 9.
In conclusion, as illustrated in Fig. 7, it is worth noting that the three lobes emptying function, with the increase of slenderness, achieves the emptying value much lower than the ones obtained using a one lobe emptying function. As a matter of fact, in the three lobes case with 40 m length, the maximum emptying percentage stopped on about 20% of the total beam height.

Conclusions
In this work, the performances of two different emptying functions (one lobe and three lobes sine functions) are tested in order to evaluate which strategy gives the best results in the proposed preliminary simplified optimization criteria. This study also aims to investigate the actual threshold values beyond which a more complex modeling of a beam with a curvilinear axis is strongly required with respect to the proposed simplified dog-bone beam model. In particular, the simplified preliminary optimization process is performed for different case studies with progressively increasing span length.
In the case of one lobe emptying function, it possible to reasonably assume that the percentages of critical emptying value beyond which the proposed dogbone model significantly deviates from the FEM solution model are about 60% of the beam total height. As a matter of fact, until this value, the error is limited within a certain engineering tolerance. A different scenario must be described in the case of the three lobes emptying function in which, also considering the optimal emptying value, the deviations of the dog-bone beam with respect to the FEM model still remain negligible. Therefore, in this range, the dog-bone beam tends to almost correctly approximate the behaviour of a non-prismatic beam with a curvilinear axis and variable cross section. Overall the analyzed emptying technique, the adoption of the using three lobes emptying function seems to The left column illustrates the optimal results adopting one lobe sine emptying function, whereas the right column contains the optimal results obtained with three lobes sine emptying function Table 1 Best optimal results obtained from the GA algorithm with a population size 20 and maximum iterations set to 50. After 10 runs, the best optimal results have been identified be the best of the one within which authors tested in this work despite it represent the hardest and the most expensive technological solution. Moreover, the different levels of emptying percentage, when slenderness increases, reached by the two emptying functions investigated is attributable both to the beam shapes and also to the adopted static scheme. A fully restrained beam tends to concentrate efforts at fixed joints in fact the one lobe emptying function, especially for beams with a high value of slenderness, tends to reduce the height section at the midspan whilst an increase of the height section occurs at the midspan with the three lobes   optimization process. This geometry results in obtaining a combination of more functions, which induced a reversal of the curvature (as depicted in Figure 7). This curvature inversion became less evident for beams with less values of slenderness as in the case of beams with 30 and 40 m of length. The optimal shape for this specific type of beam appears with an important plateau (similarly to the polycentric arch shape) localized in the central part of the beam caused by a curvature inversion: negative bending moments are greater than the positive ones at midspan and the GA algorithm is not able to perform an optimization including the curvature inversion, increasing the height section at midspan and, hence, the total weight of the beam in order to respect the geometrical and stress constraints imposed. As future developments, it will be interesting to study the behavior of beams with different slenderness ratios but realized with a material in which the traction response results different to the compression one as in the reality in order to overcome the limitations as mentioned before.
In conclusion, an optimization dimensional process is performed in order to investigate the range within which some initial hypothesis and simplified assumptions do not affect the real behavior of the structure. In this way, it is noted that variable cross section beams optimized with the three lobes emptying function regardless of their length can be used as a starting point to implement a more accurate analysis in the future.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.