Numerical Analysis Project

Numerical Method of Increasing the Critical Buckling Load for

Straight Beam-Type Elements with Variable Cross-Sections

Under Supervision of
Dr. Amany

1
Teams Members Names:

Names Sec Code

Abdallah Hesham Abdelsalam 3 9220483
Muhab Muhammed Mabrouk 6 9220881
Muhanad Hany Hamed 6 9220884
Muhanad Yasser Ibrahim 6 9220885
Youssef Abdallah Saad 6 9220992
Muhab Sayed Muhammed 6 9220879
Youssef Muhammed Muhammed 6 9211442

Youssef Essam Youssef 6 9211419

Youssef Mahmoud Said 6 9221019
Youssef Mahmoud Metamed 6 9221021
Youssef Mahmoud Attia 6 9221020

2
 Table of content:

Abstract…………………………………………………….4
Introduction…………………………………………………5
Literature review…………………………………………….7
Mathematical Modelling…………………………………….14
Work method…………………………………………………37
Conclusion……………………………………………………72
References…………………………………………………….75

3
Abstract

Structural expressionism resembles the use of slender structural elements, in particular

beam-type elements. To satisfy structural, functional, and also architectural requirements
a comprehensive structural analysis must be performed. The main issue of this study is
the buckling analysis of beam-type elements, concerning Cavalieri’s principle. The
present study is divided into two separate sections. The first part is a theoretical study, in
which a variable cross-section beam-type element is modeled.

The stability analysis is performed by an indirect variational method and the stiffness of
the support connections is also introduced. The numerical simulation highlights 6 cases
defined by the restraints of the support connections. The case study follows the
modification of the critical buckling load of the variable cross-section beam-type
element. Prior to the case study, a novel verification method is proposed to achieve a
realistic cross-section for the beam-type element.

The study revealed that with ideal characteristics of the stiffness coefficients of the
restrains significant increase of the critical buckling load is obtained, and further if an
actual situation is considered with finite values of the stiffness of the restrains, the
variable cross-section for the beam-type element is a recommended and rational choice to
make, to eliminate stability issues.

4
Introduction
General Aspects
The structural design process assumes achieving a structure that satisfies structural,
functional, and architectural requirements, the famous triad formulated by Vitruvius in his
multi-volume work De architectura in the 1st century BC. Due to the number of factors, one
should consider parametric correlation rather than individual parameters, to perform an optimal
structural design process, which in addition is an iterative one.

To quantify on a qualitative and quantitative scale the parameter correlation principle one
may choose to evaluate the carbon footprint of a structure with the Life-cycle Assessment
(LCA) technique .Through this technique, one may have an overview of the most important
stages during the life cycle of a structure, to take the most optimal solutions during the design
phase. The design and structural conformation phase have an insignificant impact on the carbon
footprint of a structure, on a quantitative scale, in comparison with the exploitation phase or the
material production phase .The structural design phase predefines the carbon footprint, on a
qualitative scale, in the means that if a structure is completed further modifications would imply
significant financial and time resources.

In the structural design phase, engineers tend to point out the most significant factors which
in the end will define the final structure: strength, stability, and stiffness. The development of
modern technologies and high-strength construction materials, allows architects and engineers
to us a smaller number of elements, therefore the presence in particular of slender beam-type
elements is more frequent.

The main goal is to avoid a structural collapse due to the loss of local or global stability, or
by extensive deformation of structural elements. Both of these undesired failure types are
closely related to structural topology, in other words, one must study the relations between
different kinds of loads and the shape of the elements.

5
 In order to achieve a comprehensive stability analysis and to compute the critical
combination of the loading parameter, a gravity and a tip force parameter may be considered.

The main issue which is a debate for more than three centuries is the study of stability, in
particular, the stability analysis of beam-type elements and frames, a bread-and-butter problem
for structural engineers according to Z.P. Bažant. One of the first pioneers of the study of elastic
buckling was Leonhard Euler who studied the strength and stability of centrally loaded
structures, under their weight. Euler’s analytical formula denotes that failure would occur due
to stresses induced by sidewise bending only, which is a particular case of long columns .Since
then the main preoccupation of engineers’ is the study of beam-type elements with material and
geometric nonlinearities. Euler’s studies have been followed by several significant theoretical
contributions to the understanding of elastic buckling, but the fathers of modern engineering
mechanics are considered to be S. Timoshenko and J.M. Gere who managed to write one of the
most comprehensive guides to the elastic stability of structures .

Centrally loaded elements in a classical interpretation are considered to be beams and

columns. Complex structures made up of beam-type elements with variable flexural rigidity are
preferably used, in addition to the architectural and functional aspects, because it provides a
better distribution of the strength and load. According to Euler’s formula, elastic buckling has a
significant role in the design of elements with great height and large spans in the case of
columns and of beams. Altogether beams do not carry a high amount of axial loads, but in the
case of tensegrity type structures, one of the main concerns of the definition of structural
topology is to reduce the length of isolated compressed elements .The present research offers a
different kind of approach for a specific situation, through a simple and high-precision method.
The main purpose of the present study is to combine stability analysis with Cavalieri’s
principle, in order to increase the critical buckling load, by modelling the geometric shape of
the beam-type element.

6
 Literature Review

-The Report includes 6 different literature reviews:

(1) Arablis and Beskos

In recent years many analytical and numerical modells have been developed for the study of
elastic buckling of structures. The feasible approach of elastic buckling has been proved to be
the numerical modells. The majority of cases do not consider the variation of material and
geometric properties. Arablis and Beskos proposed a new numerical method for the linear
elastic stability analysis of plane structures consisting of beam-type elements. Throughout the
analysis, the beams were considered with rectangular cross-sections with constant width and
variable depth.

The geometry is modeled through linearly tapered functions with constant width. The numerical
modell implies the finite element method to compute the exact flexural and axial stiffness
matrix and quasi-consistent mass and geometric matrices [1]. Das et al. proposed for the
geometric shape of the beam-type elements linearly taper, exponential and parabolic functions,
for which two types of cross-sections are adopted, solid circular and rectangular cross-sections.
In the case of rectangular cross-sections, a constant width b is maintained.

The types of structural elements that are subjected to significant axial loads are the columns,
which collect the weight of the structure and transfer it to the foundations. The columns have
become more slender, due to architectural and design requirements [3]. Topological modelling
of the beam-type element is a main issue regardless of the application field: from stability
analysis to the study of dynamic behavior with plastic deformations [2].

7
 In the case of slender columns with significant height, the establishment of the critical buckling
load is an aspect of reference in the design phase. Gradually changing cross-sections ensures
the decrease of the mass of the columns, which can be modelled throughout the moment of
inertia of the cross-section [4].

(2) Von Kármán theory

Many considerations on the numerical analysis of initial post-buckling behavior in plates and
beams have been taken. One of which is an exceedingly simple examination of the von Kármán
theory for the buckling analysis of thin plates is carried out. It is shown that the von Kármán
assumption, based on the interaction between axial and flexural deformation, is intrinsically less
“robust” in capturing the target phenomenology than the usual treatment of Euler inextensible
strut.[5]

Bearing of beams can be analyzed by means of a highly idealized model of an inextensible,

shear undeformable rod, the analysis of buckling and initial post-buckling behavior of
compressed thin plates usually takes into account the membrane strains and deformations.

This is generally done by means of the von Kármán theory [6], essentially based on an a priori
assumption on the non-linear strain–displacement relations that neglect the quadratic terms in
the derivatives of the in-plane displacements. von Kármán model provides a relatively easy and
effective way of analyzing the buckling of thin plates; but, to what extent is it accurate? This
model has been systematically employed both by means of a system of two simultaneous non-
linear differential equations and Rayleigh–Ritz energy procedures.

8
It is unquestionable from an engineering standpoint that the von Kármán theory is very
effective in capturing the phenomenology of the plate post-buckling problem and that for over a
century it has constituted a valuable tool in this extremely important topic in structural
mechanics since plates are possibly unique in their extensive use as load-carrying structural
components up to and into the post-buckling range. This theory has been also largely employed
to explain the ability of thin plates loaded in edge compression to sustain loads well above the
classical critical load in presence of abrupt changes in waveform after initial buckling on
account of non-linear coupling of buckling modes at simultaneous or near simultaneous critical
loads and these findings have been confirmed by extensive experiments.

(3) Ruocco model

Several studies investigate variable cross-section beam-type elements, in order to reduce the
weight of the structure and increase strength, or to satisfy architectural requirements. The
optimization process, through the structures own weight, is minimized and the critical buckling
load is maximized, which is the main issue regarding variable cross-section beam-type elements
.[7]

Ruocco et al.[8] proposed an optimization problem to define the exact geometric shape of
inhomogeneous columns with elastic restraints against buckling, loaded simultaneously with
concentrated and distributed compressive loads which include self-weight. In this model, the
Euler column is used, which is discretized by adopting the Hencky bar-chain model, in order to
solve the differential equations which are governing the buckling phenomena. The critical
buckling load, through this approach, represents the lowest eigenvalue of the resulting system of
algebraic equations [8]. The Hencky bar-chain model offers a physical interpretation of the
numerical finite difference method, as well as the possibility of optimization and buckling
analysis of non-uniform beam-type elements .[9]

9
(4) Leroy Gardner's study

Recent advancements in the construction industry have introduced hot-rolled and cold-formed
structural steel tubular members with elliptical cross-sections. However, there is a lack of
comprehensive understanding of their structural behavior and stability, leading to a dearth of
.design guidance

Leroy Gardner's study [10] focuses on the elastic buckling response of elliptical hollow sections
(EHS) under compression, revealing their behavior to be intermediate between circular hollow
sections and flat plates. This transition depends on the aspect ratio and relative thickness of the
.section

Through numerical and analytical investigations, the study proposes formulas to accurately
predict the elastic buckling stress of elliptical tubes [11],[12], addressing shortcomings in
existing expressions. Length effects are also explored. These findings contribute to the
development of slenderness parameters and effective section properties for slender elliptical
tubes, aiding in the classification of EHS and enhancing their structural design.

(5) Deformation theories

Over the years, the buckling behaviour of laminated composite plates has been studied
extensively by many researchers. The first order shear deformation theory (FSDT) and finite
strip are the methods that have been used considerably. Zhong and Gu (2007) presented an
exact solution for buckling of simply supported symmetrical cross-ply composite rectangular
plates under a linearly varying edge load. It was developed based on the first-order shear
deformation theory for moderately thick laminated plates. The results were verified using the
computer code ABAQUS. Ovesy et al. (2013) presented an exact finite strip for the buckling
analysis of moderately thick symmetrically cross-ply laminated composite plates and plate
structures by using First-order Shear Deformation Theory (FSDT) [13]. Ghannadpour and
Ovesy (2009) conducted theoretical developments of an exact finite strip for the buckling
analysis of symmetrically laminated composite plates and plate structures. The so-called exact
finite strip was developed based on the concept that it was effectively a plate and the von-
Kármán's equilibrium equation was solved exactly to obtain the general form of out-of-plane
buckling deflection mode for the corresponding strip.

10
The higher order shear deformation theory (HSDT) [14] has been used in various researches to
analyse buckling and bending of composite laminates. Fazzolari et al. (2013) developed an
exact dynamic stiffness element based on higher order shear deformation theory and extensive
use of symbolic algebra for the first time to carry out a buckling analysis of composite plate
assemblies. Gaiotti and Rizzo (2012) dealt with strength of composite sandwich panels
typically used in pleasure boats and naval industry. The focus was on the influence of
fabrication processes onto the mechanical behaviour, namely onto buckling strength. Both,
experimental and numerical analyses were carried out on rectangular specimens, whose
dimensions and geometry were conceived to induce buckling collapse, being this one a typical
failure mode under compressive loading for thin skins sandwich. Kheirikhah et al. (2012)
improved a high-order theory for biaxial buckling analysis of sandwich plates with soft
orthotropic core.

Third-order plate theory was used for face sheets and quadratic and cubic functions were
assumed for transverse and in-plane displacements of the core, respectively. Dash and Singh
(2012) investigated the buckling and post-buckling of laminated composite plates using higher
order shear deformation theory associated with Green–Lagrange non-linear strain–displacement
relationships. Kharghani and Guedes Soares (2015) investigated the effects of different
parameters of the composite laminate and delaminated area analytically under buckling load
based on a developed Layerwise Higher Order Shear Deformation Theory (LHSDT) and
compared the results with the three-dimensional Finite Element Analysis (FEA). Also, they
developed LHSDT for plates containing rectangular embedded delaminations undergoing non-
negligible lateral deflection effects (Kharghani and Guedes Soares, 2016). They investigated
the deflection of the composite laminates and their maximum tensile strain in bending
conditions too, using different equivalent single layer and layerwise theories based on
polynomial shape functions. Moreover, they compared the results with the three-dimensional
finite element analysis and the experimental investigation (Kharghani and Guedes Soares,
2018).

11
(6) Buckling and postbuckling of thin elastic beams
Buckling of deformable structures is quite pervasive geometric instability phenomenon and thus
it represents one of the limit states to consider in the design of structures. Specifically, buckling
in beams arises when, due to the annihilation of the effective stiffness caused by compressive
loads beyond a critical value, considerable out-of-plane bending displacements appear. The
literature on buckling and postbuckling of beams is wide and, at the same time, also scattered.
This review paper is an attempt to shed light onto the main findings especially in the context of
analytical results by recognizing the main groups of beams, based upon their boundary
conditions, which share the same models and phenomenology.

The literature presents plenty of papers that deal with buckling and postbuckling of columns. A
closer investigation of the literature would classify these papers into two main groups based on
the boundary conditions: columns with movable ends and columns with immovable ends. In
other words, the column ends may be either axially unrestrained, movable, or axially restrained,
immovable. In the first group, movable ends, one end approaches the other upon increasing the
compressive axial load during the postbuckling course of deformation.

This type of modeling enables the study of postbuckling when large deformation is involved.
No exact solution exists for the postbuckling response except for special cases where the elliptic
integrals are used, as will be shown in the following sections (Timoshenko and Gere, 1961;
Humer and Irschik, 2011a; Humer, 2013)[15]. If the deformation is relatively large, some
approximations can be applied on the bending curvature using the Taylor-series expansion and
hence an approximate solution for the postbuckling response can be found (Emam and
Lacarbonara, 2021).

12
The second group presents modeling of buckling and postbuckling while the column ends are
immovable. It is here important to notice that if the beam ends are fixed, i.e., axially restrained, increasing
the applied axial load has no effect on the buckling of the beam as this additional load will be carried by
the support reaction. In other words, one would question how it is possible to have load-deflection
relationship in the postbuckling domain in this case. In fact, the buckling problem of beams with
immovable ends is different in principle from its counterpart of beams with movable ends. The proper
way to think of buckling of beams with immovable ends is to consider for example beams under a pre-
stress, such as the case of micro-fabricated structures.

After the microfabrication process is over, and due to the residual stresses that are devoted to the
difference in thermal expansion, the beam and the substrate will be under compression. If the beam is
released from the substrate, while its ends are kept immovable, it may buckle. The buckling amplitude in
this case depends on the amount of the pre-stress. When the beam buckles while its ends are immovable,
it stretches and hence the potential energy in this case is associated with both bending and midplane
stretching. It is interesting to notice that the applied axial load does no work during the postbuckling as
the end shortening vanishes. Generally speaking, problems classified in this category involve pre-stressed
columns that are subjected to compressive loads while enforced to keep a fixed distance apart between
ends. [16].

13
Mathematical Modelling
(We have collected this data from the previous reviews of literature discussed to make our own
developed mathematical modelling)

Modelling Assumptions:
In this study the following assumptions have been made: the material, of which the column
is made, is linearly elastic; the hypothesis of small displacements and deformations; the axial
force Q is constant and is computed through a first-order analysis;
-force Q is not influenced by the bending moments, it is a conservative force, which remains
unchanged during bending; for the analysis the Euler-Bernoulli beam model is used, which does
not take into account shear deformation and rotational bending effects as the Timoshenko-
Ehrenfest beam theory; stiffness matrix assembling process neglects axial deformations; the
beam-type elements are not subjected to lateral load, external loads are applied only in the
joints; the column ax is perfectly straight, there are not taken into account any initial
imperfections; no local buckling at any cross-section along the column length is allowed.

General Aspects
A major issue in the design of beam-type elements with significant axial loads is to
improve the stability of the element and avoid their failure through buckling, which
phenomenon threatens the structural integrity of the whole structural assembly. The
phenomenon is more relevant in the case of statically determined structural assemblies.
The first step in the modelling of beam-type elements to improve their stability is to
analyze Euler’s critical load derived in 1757, which is presented in Expression(1) .
π2·E·Imin
𝐹𝑐𝑟,= 𝑙𝑓2
(1)

In Expression (1) E represents Young’s modulus, 𝐼𝑚𝑖𝑛 denotes the minimum value of the
moment of inertia of the cross-section, 𝐿 𝑓 and represents the buckling length of the beam-type
element subjected to compression forces, with respect to the imposed kinematic constraints to
the joints.

14
 Euler’s formula highlights two important aspects unaware of the kinematic restrains of the
joints, the length of the element and the moment of inertia of the cross-section. The factor of
direct proportionality with the critical load is the moment of inertia. This design parameter gain
importance when the buckling modes are studied. The deflection shape associated with the
primary buckling mode emphasizes potential critical areas of the element, which by similarity,
can be scribed as potential areas of appearance of plastic hinges. In this manner it is required to
amplify the value of the moment of inertia, therefore a variable cross-section beam-type
element is obtained. In order to avoid sudden changes of the elements cross-section, which can
be identified as stress concentration points, it is proposed the modelling of the cross-section of
the element through a continuous function.

Modelling the cross-section of the element is performed with respect to Cavalieri’s

principle. According to Cavalieri’s principle, the cross-section area of every section will remain
unchanged during the modelling process, but the moment of inertia will change in accordance
with the primary buckling mode. As a consequence of the mentioned principle, the elastic
resistance of every cross-section is identical. Modelling a beam-type element with a variable
cross-section follows the increase of the critical elastic buckling load, in an ideal case
exceeding the elastic resistance of the cross-section, thus the buckling phenomenon can be
avoided.

Modelling of the Moment Inertia

Modelling the moment of inertia of the cross-section presumes to assume a variation law.
The variation law can be formed through interpolation functions between reference values in
form of piecewise, multilinear, or polynomial functions. For the first two types of function,
results are satisfying, although in the points where the function is not continuous the iterative
process may be divergent, and the precision decreases by two orders of magnitudes. In the
particular case of polynomial functions, the passing through the critical values is continuous,
and the desired precision is in close correlation with the polynomial function degree.
The modelling of the moment of inertia in this study regardless of the kinematic constraints
of the joints is performed with cubic functions, a piecewise polynomial formula known as the
spline function. The natural boundary conditions are shown in Figure 1. Be f a function defined
on [𝑎 𝑏] which is interpolated in m + 1 nodes on n sub-intervals through a cubic interpolant S(x)
defined in Figure 1.

15
Figure 1. Cubic interpolant 𝑆𝑗(𝑥), defined on sub-interval [𝑥𝑗 𝑥𝑗+1]

Cubic interpolant S(x), is a cubic polynomial piecewise function expressed in a general

form as in Expression (2)
(𝑥)=⋃𝑚−1
𝑗=0 𝑆(𝑥) (2)

In Expression (2), 𝑆(𝑥) is defined on sub-interval [𝑥𝑗 𝑥𝑗+1],∀𝑥𝑗 ∈0,𝑚−1¯. The boundary
conditions for every sub-interval [𝑥𝑗 𝑥𝑗+1] are presented in Relations (3a)–(3c).
𝑆𝑗(𝑥𝑗 )=𝑓(𝑥𝑗 ); 𝑆𝑗(𝑥𝑗+1)=𝑓(𝑥𝑗+1 );∀𝑗 ∈0,𝑚−1 −“Left”&“Right”Interp. (3a)
𝑆𝑗′(𝑥𝑗+1 )=𝑆𝑗+1′(𝑥𝑗+1 );∀𝑗 ∈0,𝑚−2 -slope−match (3b)
𝑆𝑗′′(𝑥𝑗+1 )=𝑆𝑗+1′′(𝑥𝑗+1 );∀𝑗 ∈0,𝑚−2 -curvature−match (3c)
The degree of the cubic interpolant 𝑆𝑗(𝑥) is 3, 𝑑𝑒𝑔(𝑆𝑗(𝑥))=3, and the general expression
for 𝑆𝑗(𝑥) ∀𝑥𝑗 ∈0,𝑚−1 is presented in Expression (4) .

𝑆(𝑥)=𝑎𝑗+𝑏𝑗(𝑥−𝑥𝑗)+𝑐𝑗(𝑥−𝑥𝑗)2+𝑑𝑗(𝑥−𝑥𝑗)3 (4)

According to Expression (4), it is required to fix the values for

coefficient 𝑎𝑗, 𝑏𝑗, 𝑐𝑗 and 𝑑𝑗with respect to the boundary conditions expressed in Relation (3a)–
(3c).
One of the cubic interpolant’s properties is to minimize the definite integral, presented in
Expression (5), for every g function defined on an interval [𝑎 𝑏], which approximates
function f through nods 𝑥 .
𝑋𝑛
𝐼=∫𝑋0 (f″(x))2 dx (5)

16
 One can introduce notation ℎ𝑗=(𝑥𝑗+1−𝑥𝑗), with respect to the general form presented in
Expression (4), thus the boundary conditions from Relation (3a)–(3c) can be defined as
recursive sequences, presented in Relation (6a)–(6c).
𝑆 j+1(𝑥𝑗+1)=𝑆𝑗(𝑥𝑗+1 ) →𝑎𝑗+1=𝑎𝑗+𝑏𝑗h𝑗+𝑐𝑗h𝑗2+𝑑𝑗h𝑗3 (6a)
𝑆(𝑥𝑗+1)=𝑆𝑗+1(𝑥𝑗+1) →𝑏𝑗+1=𝑏𝑗+2𝑐𝑗h𝑗+3𝑑𝑗h𝑗2 (6b)
𝑆𝑗′(𝑥𝑗+1 )=𝑆𝑗+1′′(𝑥𝑗+1 ) →𝑐j+1=2𝑐𝑗+6𝑑𝑗h𝑗 (6c)
The numerical value of coefficient 𝑎𝑗 results from the function evaluation of f at the input
value of 𝑥𝑗, denoted with 𝑓(𝑥𝑗 ) as for 𝑎𝑗+1 results from the function evaluation of f at the
input value of 𝑥𝑗+1, denoted with 𝑓(𝑥𝑗+1 ) . For the recursive sequences from Relations (6a)–
(6c) it is required the cubic interpolants coefficients as expressions which depends only on two
specific arguments, thus coefficients 𝑎𝑗 and 𝑐𝑗 are chosen. Coefficient 𝑑𝑗, can be expressed
from Relation (6c), presented in Expression (7).
𝐶𝑗+1−𝐶𝑗
𝑑𝑗= 3ℎ𝑗
(7)

Substituting Expression (7) in the recursive Sequence (6a,b), Expression (8a,b) can be
expressed .
ℎ𝑗 2
𝑎𝑗+1=𝑎𝑗+𝑏𝑗h𝑗+ (2Cj + Cj+1) (8a)
3

𝑏𝑗+1=𝑏𝑗+h(𝑐𝑗+𝑐𝑗+1) (8b)
Substituting Expression (8a) in Expression (8b) it is possible to express coefficient 𝑏𝑗 with
arguments 𝑎𝑗 and 𝑐𝑗, as in Expression (9) .
1 3
𝑏𝑗=hj(𝑎𝑗+1−𝑎𝑗)− hj (2𝑐𝑗+𝑐𝑗+1) (9)

The recursive Sequences (8b) and (9) can be expressed also for sub-interval [𝑥𝑗−1 𝑥𝑗],
therefor Expression (8b) can be written as in Equation (10) [21].
3 3
ℎ𝑗−1𝑐𝑗−1+2(ℎ𝑗−1+ℎ𝑗)𝑗+ℎ𝑗𝑐𝑗+1=hj (𝑎𝑗+1−𝑎𝑗)− hj−1 (𝑎𝑗−𝑎𝑗−1) ∀𝑗= 1,𝑚−1

Equation (10) represents m - 1 equations, and in order to solve a linear system

of equations, two additional relations are required in the form of the natural boundary

conditions expressed in Relation (11)

(eq.11a) 𝑠0′′ (𝑥0 ) = 2𝐶0 = 0 → 𝐶0 = 0

17
 ′′ (
(eq.11b) 𝑠𝑚 𝑥𝑚 ) = 2𝐶𝑚 = 0 → 𝐶𝑚 = 0

Therefore, the linear system of equations, whit respect to Equation (10) and Condition

(11a,b), can be expressed as a matrix-vector Equation (12)

(eq.12) [A]{c}={a}

In Equation (12), [A] is the (m + 1) _ (m + 1) coefficient matrix of the system whose

implicit from is presented in Expression (13), fcg is a vector made up of the m + 1 unknowns,

presented in Expression (14a), and fag is the vector made up of the m + 1 right-hand sides

of the equations, presented in Expression (14b).

18
(eq.14a) {𝑐}𝑇 = {𝑐0𝑐1 …𝑐𝑚−1 𝑐𝑚}
3(𝑎2 −𝑎1 ) 3(𝑎1 −𝑎0 ) 3(𝑎𝑚 −𝑎𝑚−1 ) 3(𝑎𝑚−1 −𝑎𝑚−2 )
(eq.14b) {𝑎}𝑇 = {0 − … − 0}
ℎ1 ℎ0 ℎ𝑚−1 ℎ𝑚−2

Therefore, the solution of System (12) specifies the numerical values for coefficient

cj. The cubic interpolant’s coefficients defined on n sub-intervals, are presented in Expressions

(15a)–(15c)

(eq.15a) {𝑎𝑗 }; ∀𝑗 = 0… m−1; 𝑎𝑗 = 𝑓(𝑥𝑗 ); 𝑎𝑗+1 = 𝑓(𝑥𝑗+1 ); 𝑥𝑗 = 0, 𝑚 − 1;

1 ℎ𝑗
(eq.15b) {𝑏𝑗 }; ∀𝑗 = 0… m−1; 𝑏𝑗 = (𝑎𝑗+1 − 𝑎𝑗 ) − (2𝑐𝑗 + 𝑐𝑗+1);
ℎ𝑗 3

𝑐𝑗+1 −𝑐𝑗
(eq.15c) {𝑑𝑗 }; ∀𝑗 = 0… m−1; 𝑑𝑗 = ;
3ℎ𝑗

The matrix-based formulation of the stability problem assumes the definition of the

displacement field, which represents a function approximation problem. The task is to find a

function that closely matches the displacement field. In the case of a bi-dimensional stability

problem, for the beam-type element presented in Figure 2, the approximation problem

assumes the approximation of the displacement field through a polynomial function of

degree n - 1, presented in Expression (16).

𝑛−1
(eq.16) 𝑣 (𝑥 ) = 𝛼1 + 𝛼2𝑥 + 𝛼3𝑥 2 + ⋯ + 𝛼𝑘 𝑥 𝑘+1 + ⋯ 𝛼𝑛 𝑥 𝑛+1 = ∑𝑘=1 𝛼𝑘 𝑥 𝑘+1 = [𝐹 ]{𝛼 }

19
In Expression (16), v(x) is the polynomial function which describes the displacement

field, n 􀀀 1 = deg(v(x)) is the degree of v(x), fFg = n1 x x2 . . . xk+1 . . . xn􀀀1 O is a vector

with the terms of the polynomial function divided by the terms’ coefficient;

fag = fa1, a2, . . . , ak, . . . , angT a vector with the polynomial functions coefficient. For

joints i and j of the beam-type element, the kinematic boundary conditions can be written

as in Relation (17).

(eq.17)

𝑣(0) = 𝑣𝑖
𝑥 = 0{ 𝑑𝑣(𝑥 )
𝜙(0) = | 𝑥 = 0 = 𝜙𝑖,
𝑑 (𝑥 )
𝑣 (𝑙) = 𝑣𝑗
𝑥 = 𝑙{ 𝑑𝑣(𝑥 )
𝜙(𝑙) = | 𝑥 = 1 = 𝜙𝑗,
𝑑 (𝑥 )

(eq.18)

{𝑎𝑒 }𝑇 = {𝑣𝑖 =𝑎1 𝜙𝑖 = 𝑎2 𝑣𝑗 = 𝑎3 𝜙𝑗 = 𝑎4 … 𝑎𝑛 }

(eq.19)

{𝐹𝑒 }𝑇 = {𝐹𝑦𝑖 = 𝑓1 𝑀𝑖 = 𝑓2 𝐹𝑦𝑗 = 𝑓3 𝑀𝑗 = 𝑓4 𝑓5 . . . 𝑓𝑛}

20
Expression (18) highlights the nodal fa1:4 g and non-nodal fa5:n g parameters. Based

on Expressions (18) and (19) the limit displacement boundary conditions are expressed
inEquation (20).

(eq.20)

𝑣 (0) = 𝛼1 + 𝛼2 ⋅ 0 + 𝛼3 ⋅ 02 + ⋯ + 𝛼𝑘 ⋅ 0𝑘+1 + ⋯ + 𝛼𝑛 ⋅ 0𝑛−1 = 𝑣𝑖

𝜙(0) = 𝛼2 + 2𝛼3 ⋅ 0 + ⋯ + (𝑘 + 1). 𝛼𝑘 . 0𝑘 + ⋯ (𝑛 − 1). 𝛼𝑛 . 0𝑛−2 = 𝛷𝑖
𝑣(l) = 𝛼1 + 𝛼2 ⋅ 𝑖 + 𝛼3 ⋅ 𝑙2 + ⋯ + 𝛼𝑘 ⋅ 𝑙𝑘+1 + ⋯ + 𝛼𝑛 ⋅ 𝑙𝑛−1 = 𝑣𝑗
𝜙(𝑙) = 𝛼2 + 2𝛼3 ⋅ 𝑙 + ⋯ + (𝑘 + 1). 𝛼𝑘 . 𝑙𝑘 + ⋯ (𝑛 − 1). 𝛼𝑛 . 𝑙𝑛−2 = 𝛷𝑗
𝛼5 = 𝑎5
𝛼6 = 𝑎6
…
{ 𝛼𝑛 = 𝑎𝑛

The linear system of Equations presented in (20), is figured as a matrix equation in

Equation (21).
[𝐻]{α}={𝑎𝑒}
(21)
In Equation (21), the matrix [𝐻] represents the coefficient matrix of the limit displacement
boundary condition system. By expressing the vector of the polynomial function coefficients and
by introducing it in Expression (16), the interpolation function matrix [𝑁] can be defined, as in
Expression (22) and the polynomial function which describes the displacement field can be
defined as in Expression (23).
𝑣(𝑥)=[𝐹][𝐻]-1{𝑎𝑒}=[𝑁]{𝑎𝑒}
(22)
𝑣(𝑥)=𝑁1(𝑥)𝑎1+𝑁2(𝑥)𝑎2+…+𝑁k(𝑥)𝑎k+…+𝑁n(𝑥)𝑎n=[𝑁]{𝑎𝑒}
(23)

21
 In Expression (23), [𝑁]=[𝑁1(𝑥) 𝑁2(𝑥)…𝑁𝑘(𝑥)…𝑁𝑛(𝑥) ]1≤𝑘≤𝑛 represents the nodal
interpolation function matrix, which links the displacements of the beam-type element with the
nodal displacements, [𝑎𝑒]=[𝑎1 𝑎2…𝑎𝑘…𝑎𝑛 ]T 1≤𝑘≤𝑛 is the vector of nodal and non-nodal

displacements. Given Expression (23), Figure 3 is highlights the significance of non-nodal

parameters/displacements through 4 interpolation functions of a degree n − 1 = 9 polynomial
function.

Figure 3. Interpolation functions Nk(x), for k = [1 3 7 10]

Given Expression (23) and the geometric meaning of the differential relations between
displacements, presented in Figure 4, the slope of the deflected beam axis φ(𝑥), expressed in
Expression (24), the strain ε(𝑥)ε, expressed in Expression (25), and the reaction forces induced
in the beam-type element σ(𝑥). expressed in Expression (26), can be defined as a function of
nodal and non-nodal displacements {𝑎𝑒}.

22
 Figure 4. Geometric representation of the differential relations between displacements,
where dx—is the unit length of the beam-type element; dv—the unitary displacement; Q—the
axial force.
dv(x) d[N]
φ(𝑥) = = {𝑎e} = [𝑆]{𝑎e}
dx dx

(24)
In Expression (24), [𝑆] represents a matrix that links the slope of the deflected beam
axis φ(𝑥) with the nodal and non-nodal displacement vector {𝑎𝑒}
−dφ(x) −d[S]
ε(𝑥)=−χ(𝑥)= = {𝑎𝑒} = [𝐵]{𝑎𝑒}
dx dx

(25)
In Expression (25), [𝐵] represents the strain displacement matrix, which can define the
curvature χ(𝑥) depending on the nodal and non-nodal displacements {𝑎𝑒}.
σ(𝑥)=𝑀(𝑥)=𝐸·𝐼·χ(𝑥)=𝐸·𝐼·[𝐵]{𝑎𝑒}
(26)
Expression (26) is based on the validity of Hooke’s law and the hypothesis of small
displacements and deformations.
For the beam-type element, presented in Figure 2, the equilibrium is expressed through the
principle of virtual work. According to the mentioned principle, the virtual displacements or the
virtual forces are infinitesimal and consistent with the physical and elastic constraints of the
beam-type element. 𝑣Therefore for virtual displacements and the virtual forces, the same
relations can be expressed as for the real displacements and forces, presented in Expression (27).
𝑣̅ (𝑥)= [𝑁]·{𝑎̅𝑒}

̅ (𝑥)= [𝑆]·{𝑎̅𝑒}
φ

ε̅ (𝑥)= χ̅ (𝑥) = [𝐵]·{𝑎̅𝑒}

̅ (𝑥)=𝐸·𝐼·[𝐵]·{𝑎̅𝑒}
̅(𝑥)=M
σ

(27)
In Expression (27), {𝑎̅𝑒} is a vector of the nodal and non-nodal virtual displacements. The
mathematical expression of the principle of virtual work is presented in Relation (28) .

23
 𝐿̅𝑒𝑥𝑡=𝐿̅𝑒𝑓
(28)
In Relation (28), 𝐿̅𝑒𝑥𝑡denotes the virtual mechanical work of the external forces,
and 𝐿̅𝑒𝑓 represents the virtual mechanical work of the internal efforts. The mathematical
expression of the virtual mechanical work of the external forces is presented in Expression (29).
𝐿̅𝑒𝑥𝑡=𝐿̅𝑒𝑥𝑡({𝐹𝑒})+𝐿̅𝑒𝑥𝑡(𝑄)
(29)
In Expression (29), 𝐿̅𝑒𝑥𝑡({𝐹𝑒}) represents the virtual mechanical work of the external forces
applied at the joints, 𝐿̅𝑒𝑥𝑡(𝑄) denoted the virtual mechanical work of the axial force Q, presented
in Figure 4. The virtual mechanical work of the external forces applied at the joints, 𝐿̅𝑒𝑥𝑡(𝑄) can
be expressed depending on the virtual nodal and non-nodal displacements, presented in
Expression (30).
𝐿̅𝑒𝑥({𝐹𝑒})={𝑎̅𝑒}T{𝐹𝑒}
(30)

The virtual mechanical work of the axial force Q, in the case of an infinitesimal beam-type
element of length dx, presented in Figure 4, is defined in Expression (31).

(eq.31) Lext (Q) = ∫ dLext (Q)

In Expression (31), the elementary virtual mechanical work of the axial force Q, can be
expressed as in Expression (32).

(eq.32) dLext (Q) = (dv (x)/dx)T Qdv = -{ae}T [S]T Q [S] {ae}dx

By substituting Expression (32) in Expression (31) one can obtain the expression of the virtual
mechanical work of the axial force Q, as a function of the nodal and non-nodal displacements,
which is shown in Expression (33).
𝑙
(eq.33) Lext (Q) =-{ae}T Q ( ∫0 [S]T[S]dx) {ae} = -{ae}T Q [Kge] {ae}

In Expression (33), [𝐾𝑔𝑒] represents the geometric stiffness matrix, which depends on the
geometric characteristics of the beam-type element. The mathematical expression of the virtual
mechanical work of the internal efforts is presented in Expression (34).

24
(eq.34) Lef = Lef_internal +Lef_s.joint

In Expression (34), Lef_internal represents the virtual mechanical work of the internal efforts,
which is expressed in Expression (35), and Lef_s.joint represents the virtual mechanical work of
the reaction forces of the support connections, presented in Expression (36).
𝑙 𝑙
(eq.35) Lef_internal = ∫0 {𝜀 }T {𝜎 }dx = {ae}T ( ∫0 [B]T .E .I .[B] dx) {ae} = {ae}T [Ke] {ae}

In Expression (35), [𝑘𝑒] represents the material stiffness matrix of the beam-type element.

(eq.36) Lef_s.joint = {v (0)}T r1{v(0)} + {𝜑 (0)}T r2{𝜑 (0)} + {v(l)T r3{v(l)} +{𝜑 (l)}T r4{𝜑 (l)}

In Expression (36), 𝑟𝑖=1:4 denotes the value of the stiffness of the elastic support connections
with respect to the kinematic boundary conditions, presented in Figure 5.

Figure 5. Representation of the beam-type element with elastic support connections, where L—
is the length of the beam-type element; x—marks the distance to a current cross-section; X—is
the local reference system axes along the axes of the beam-type element; y—is the local
reference system axes perpendicular to the beam-type element axes; I—the moment of inertia
of the cross-section in a general form; I(x)—the moment of inertia of cross-section at section x;
Q—is the axial force; r1—denotes the stiffness for a hindered transversal displacement for joint
i; r2—denotes the stiffness for a hindered rotation for joint i; r3—denotes the stiffness for a
hindered transversal displacement for joint j; r4—denotes the stiffness for a hindered rotation
for joint j.

Expressions (35) and (36) are quadratic forms with nodal and non-nodal displacements as
variables, therefore the stiffness coefficients can be obtained by the algebraic sum of the
quadratic form’s coefficients, with respect to the connection indexes. Based on Expressions

25
(35) and (36) the virtual mechanical work of the internal efforts can be written as in Expression
(37).

Regarding Relationships (35) and (36), the procedure for assembling the material stiffness
matrix, Expression (35) can be written as in Expression (37).

(eq.37) Lef = {ae}T[ke] {ae}

In Expression (37), [𝐾𝑒]denotes the assembled material stiffness matrix. Based on Relation (28)
and Expression (29) Equation (38) can be expressed.

(eq.38) Lef – Lext (Q) = Lext ({Fe})

Equation (39) can be obtained by substituting Expressions (30), (33), and (37) in Equation (38).

(eq.39) {ae}T ([Ke] -Q [Kge]) {ae} = {ae}T {Fe}

Equation (39) must be fulfilled for every virtual displacement {𝑎̲𝑒}≠ {0}, therefore Equation
(40) can be written.

(eq.40) ([Ke]- Q [kge]) {ae} ={Fe}

A compact form of Equation (40) is presented in Equation (41).

[𝐾]{𝑎𝑒}={𝐹𝑒}
(41)
In Equation (41), [𝐾] represents the elastic stiffness matrix of a beam-type element.

Formulation of the Stability Analysis

Based on Equation (41), internal forces {𝐹𝑒}, can be obtained with respect to the nodal and
non-nodal displacements {𝑎𝑒} and the axial force Q. By expressing the nodal and non-nodal
displacement vector from Equation (41), Expression (42) is obtained .
{𝑎𝑒}=[𝐾]−1{𝐹𝑒}=1|[𝐾]|[𝐾]∗{𝐹𝑒}
(42)

26
 In Expression (42), [𝐾]* represents the adjunct matrix of the elastic stiffness matrix, whose
diagonal entries are the determinant of the elastic stiffness matrix [𝐾]. Based on Expression
(42), the failure of a beam-type element occurs when there is an infinitesimal increase in the
axial force 𝑄 the displacements tend to infinity. Therefore, in the case of beam-type elements,
the loss of stability occurs when Equation (43) is satisfied, called an eigenvalue problem .
|[𝐾]|=|[𝐾𝑒]−𝑄[𝐾𝑔𝑒]|=0
(43)
Equation (43), denotes the characteristic equation of the eigenvalue problem, which leads
to a homogeneous system of equations, presented in Equation (44) .
([𝐾𝑒]+[𝐾𝑔𝑒]){𝑣}={0}
(44)
The solution of the eigenvalue problem and of Equation (44), represents vector {λ}λ with
the proper eigenvalues λ𝑖=1 ̅̅̅̅̅
1: n and the corresponding eigenvectors in a matrix {𝑣𝑖}. The
minimum value of {λ}λ represents λ1λ1 which denotes the critical buckling load, and the
corresponding eigenvector {𝑣1} the deflection shape.

Indirect Variational Method for the Study of the Conditions and Nature of
Equilibrium

The potential energy is defined as a function of a function, that is, π=π(𝑣(𝑥))π=π, a

functional presented in Expression (45). The problem consists of determining the condition of
the minimum functional ππ of one function v of one variable x .
π[𝑣(𝑥)]=ₒ¹∫ ϕ(𝑥,𝑣,𝑣′,𝑣″)𝑑𝑥
(45)
In Expression (45), v(x) is presumed to be continuous, with continuous first four derivates
and also appropriate boundary conditions. The mathematical form of the total potential energy
is presented in Expression (46) .
π=𝑈−𝑉
(46)
In Expression (46), U denotes the strain energy, V represents the potential of the external
forces or the mechanical work of external force applied in joints. In this study, U is defined in
Expression (47) through two energetic quantities.
𝑈=𝐿𝑒𝑓−𝐿𝑒𝑥𝑡(𝑄)

27
 (47)
It is assumed that the equilibrium position corresponds to an extreme value of the total
potential energy, a condition presented in Equation (48).
δπ={0}
(48)
Equation (48) represents the condition for the Euler equation, with the proper boundary
conditions. In mechanics, Equation (48) is generally equivalent to the equilibrium condition.
For the study of the nature of equilibrium, there are considered two stages of deformation, for
the intensity λ of the axial force: stage 1, which is defined by displacements 𝑣(𝑥) and total
potential energy π(𝑣(𝑥),λ); stage 2, which corresponds to displacements 𝑣(𝑥) + δ𝑣(𝑥) and total
potential energy π(𝑣(𝑥)+ δ𝑣(𝑥),𝜆 ).
The variation of the total potential energy by passing from stage 1 to stage 2 can be
expressed in Expression (49) .
Δπ=π(𝑣(𝑥)+δ𝑣(𝑥),)−π(𝑣(𝑥),𝜆)
(49)
In Expression (49), ΔπΔπ is a continuous n-time differentiable function defined on the
interval [0 L], which is expanded into a Taylor series, presented in Expression (50).
Δπ=δπ+δ2π+δ3π+…=δπ+δ²π+𝑇𝑠(π) (50)
In Expression (50), δπ set the first-order derivative of the total potential energy, δ²π denotes
the second-order derivative of the total potential energy and 𝑇(π) represents the higher-order
derivative of the total potential energy, terms which are neglected in the modelling process .
The beam-type element in this study has a single degree of freedom v(x), which defines the
system’s state. The terms obtained after the expansion into a Taylor series of ΔπΔπ, are
presented in Expressions (51a)–(51c) .

1 ∂π(v,λ)
δπ = ∑i ∂vi (51a)
1! ∂vi

1 ∂2π(v,λ)
δ2π = ∑i∑j ∂vi ∂vj (51b)
2! ∂vi ∂vj

28
 1 ∂3π(v,λ)
δ3π = ∑i∑j∑k ∂vi ∂vj ∂vk (51c)
3! ∂vi ∂vj ∂vk

With respect to the condition of equilibrium, presented in Equation (48), and to the fact that the
potential of the external forces, for the case of constant loads, is expressed as 𝑉=∑1≤𝑖≤𝑛𝐹𝑖·𝑣𝑖
where 𝐹𝑖 represents a force associated with displacement 𝑣𝑖, which is a linear displacement
function, whose high-order derivative is canceled, the variation of the total potential energy
between stage 1 and 2 is equivalent to the second-order derivative of the strain energy,
presented in Expression (52) .

1 ∂2U(v,λ)
Δπ = 12 · δ2U = ∑i∑j ∂vi ∂vj (52)
2! ∂vi ∂vj

Expression (52) allows the study of the nature of equilibrium, as presented in Table 1.

Table 1. Study of the nature of equilibrium according to the sign of δ2 .

case The State of A Short Description of

Equilibrium in the State of
Mathematical Terms Equilibrium
1. 𝑠𝑖𝑔𝑛(δ2𝑈)<0
In stable equilibrium, the
second-order derivative

29
 of the strain energy is
positive-definite, with
respect to the Lagrange-
Dirichlet theorem.

2. 𝑠𝑖𝑔𝑛(δ2𝑈)>0 In unstable equilibrium,

the second-order
derivative of the strain
energy is not positive-
definite, according to the
Liapunov stability
theorem.
3. δ2𝑈=0 Neutral equilibrium,
corresponds to a limit
state immediately before
the loss of equilibrium.

In the case of neutral equilibrium, it is necessary to determine and study the higher-order
derivative of the total potential energy according to Expression (50). The nature of equilibrium
is indicated by the fourth-order derivative of the total potential energy, which has to be positive-
definite, positive definiteness of ΔπΔπ guarantees stability .
Based on the first theorem of Castigliano the second-order derivative of the total potential
energy with respect to the displacements represents the stiffness coefficient 𝑘𝑖𝑗, which indicates
a force at i due to unit displacement at j. According to this remark, Expression (52) can be
expressed as in Expression (53) .

1
Δπ = ∑i∑jkij ∂vi ∂vj (53)
2!

Expression (53) is a quadratic from in displacements, which in a matrix form is presented in

Expression (54).

1
Δπ=2{𝑣}𝑇[𝐾]{𝑣} (54)

30
 In Expression (54), {𝑣}1≤𝑖≤𝑛 is the displacement vector, [𝐾] is the elastic stiffness matrix
of the studied structural element. Based on the condition that the variation of the total potential
energy at the limit equilibrium state must be zero, respectively [𝐾] is obtained through the
assembly of the material stiffness matrix and the geometric stiffness matrix, the equation of
stability can be determined by an energetic method, which assumes an eigenvalue problem .
Continuous structures are analyzed through an indirect variational method, in which the
structure is not discretized. The differential equations are obtained from the minimizing
condition for the potential energy .In contrast to this approach, the stability analysis, in which
the structure is discrete or discretized and the expression of the second variation of potential
energy is reduced to a quadratic form, is the direct variational approach .

Positive Definiteness of Δπ
The positive definiteness of Δ𝜋 guarantees stability. The total potential energy is quadratic from
in displacements, and the elastic stiffness matrix [K] determines the quadratic form. Matrix [K]
is a real symmetric matrix, which is positive-definite if and only if the eigenvalues of matrix
[K] are positive .

The elastic stiffness matrix [K] is a real, Hermitian matrix and, according to Sylvester’s
criterion, the necessary and sufficient condition for the quadratic from ΔπΔπ to be positive-
definite if the leading principle minors are positive, meaning Δ𝑘>0, for ∀𝑘=1:𝑛¯, where n
represents the dimension of the real vector space of the elastic stiffness
matrix, 𝑑𝑖𝑚([𝐾])=𝑑𝑒𝑔(𝑣(𝑥))+1 .

Boundary Conditions

According to the Lagrange-Dirichlet theorem, the equilibrium is stable if the total potential
energy admits a local minimum. The condition that the potential energy must be strict minimum
determines the differential equation which describes the problem and the admissible forms of
the boundary conditions. The kinematic boundary conditions in mechanics or the essential
boundary conditions in mathematics can be expressed in Relation (55) .

31
 𝑣 0
{ }={ } (55)
Ǿ 0

Equation (56a,b) can be used to express static boundary conditions in mechanics or natural
boundary conditions in mathematics .

(𝐸𝐼𝑣″)′−𝑃𝑣′=𝑉=0 (56a)

𝐸𝐼𝑣″=𝑀=0 (56b)

The approximate deflection shape v(x), or trial function, must always satisfy all the kinetic
boundary conditions but does not have to satisfy the static boundary conditions. The static
boundary conditions and the differential conditions of equilibrium are a result of the
minimization of ππ itself. If the static boundary conditions are satisfied, v(x) approximates the
exact buckling shape, and the resulting critical load 𝑃𝑐𝑟1 approximates 𝑃𝑐𝑟 more accurately.

Rayleigh Quotient

The strain energy is a quadratic form in displacements and the neutral equilibrium state assumes
that the second-order variation of the total potential energy is zero, conditions presented
through a mathematical from in Equation (57) .

32
 1
δ2U = {v(x)}T[K]{v(x)} = {0} (57)
2

Based on Equation (48) and Expression (23), Equation (57) can be expressed as in Equation
(58).
{𝑎𝑒}([𝐾𝑒]−𝑅𝑄[𝐾𝑔𝑒]){𝑎𝑒}={0} (58)

In Equation (58) 𝑅𝑄 represents the axial force multiplicator, which is expressed in

Expression (59).
{𝑎𝑒}𝑇[𝐾𝑒]{𝑎𝑒}
𝑅𝑄 = (59)
{𝑎𝑒}𝑇[𝐾𝑔𝑒]{𝑎𝑒}

Expression (59) represents the Rayleigh quotient, which approximates the axial force, value
dependent on the displacement field v(x). In order to fix the value of RQ, an iterative process
must be launched which has as its starting value a certain displacement vector, a vector which
readjusts after every iteration, with local convergence conditions. The defection shape at the
moment of loss of equilibrium is not known, therefore the presented method is an approximative
one.
According to the Lagrange-Dirichlet condition, the equilibrium is stable if P <𝑃𝑐𝑟 for every
displacement v(x) compatible with the system’s kinematic constrains. At the limit state, the
moment before the loss of stability, the critical buckling force
represents 𝑃𝑐𝑟=𝑚𝑖𝑛({𝑅𝑄}1≤𝑖≤𝑖𝑡), where it defines the number of iterations.

In Expression (59), the numerator represents the mechanical work of the internal efforts 𝐿𝑒𝑓
and the denominator the unitary mechanical work of the axial force, in this case, RQ, 𝐿𝑒𝑥𝑡(𝑅𝑄)
Expression (59) for a beam-type element can be expressed as in Expression (60) .

1 𝐿
𝐿𝑒𝑓 ·∫ 𝐸·𝐼(𝑥)· 𝑣(𝑥)′′ 2𝑑𝑥
2 0
𝑅𝑄(𝑣(𝑥)) = = 1 𝐿 (60)
𝐿𝑒𝑥𝑡(𝑅𝑄) .∫ v(x)′𝟐dx
2 0

33
 In Expression (60), the Rayleigh quotient is defined as an analytic expression concerning
the displacement field v(x), and function I(x) is a continuous function that describes the
modification of the moment of inertia of the beam-type element .

Timoshenko Quotient

In the case of statically determined structures with the Timoshenko quotient, one may
obtain a better approximation of the critical buckling load. The method assumes the appraisal of
the strain energy on the basis of the bending moment M, which appears in the beam-type
element due to deformation. The value of the bending moment is approximated with the
product −𝑄·(𝑥) instead of the curvature, which assumes an approximative displacement field .
In the case of a structure with linear elastic behavior, the strain energy is a quadratic form,
and if (𝑥)={0}it corresponds to the equilibrium position, then the modification of the total
potential energy at the limit of equilibrium coincides with the potential energy. The
modification of the total potential energy at the critical state must be equal to zero, a condition
which can be expressed as in Equation (61).

𝐿 (Q·v(x))2 𝐿𝑄
Δπ = 𝐿𝑒𝑓−𝐿𝑒𝑥(𝑄) = ∫0 𝑑𝑥 - ∫0 . v(x)′ 2dx == Q2 · 𝐿𝑒𝑓−𝑄·𝐿𝑒𝑥𝑡(𝑄) =0 (61)
2.𝐸.𝐼 (𝑥) 2

By expressing the axial force Q from Equation (61), Expression (62) is obtained .

1 𝐿
𝐿𝑒𝑥𝑡(𝑄) · ∫ 𝑣(𝑥)′ 2𝑑𝑥
2 0
𝑇𝑄(𝑣(𝑥)) = = 1 𝐿 (62)
𝐿𝑒𝑓 ·∫ 𝑣(𝑥)2𝐸·𝐼(𝑥)𝑑𝑥
2 0

34
 Work Method:
General Aspects
In this study, it is assumed that the beam-type element is part of a new structural
assembly or part of a structure that requires structural rehabilitation. In the case of a
new structure, a main structural design criterion is to use a minimum number of
elements and a minimum amount of material. Therefore, the design or modelling of
the beam-type element assumes an iterative process of optimization, through which
the main goal is to fix the geometric characteristics of a beam-type element, of length
L, subjected to axial force, with different kinematic constraints at the joints. In this
study loads are applied only at the joints, therefore the state of strain is easier to
control. The final failure criteria is a failure due to excessive stresses rather than
buckling.

The modelling process assumes three interconnected iterative cycles, presented in

Figure 6, which assumes obtaining a geometric configuration for the beam-type
element which is not sensible to buckling.

The first iterative cycle has the purpose of fixing the maximum critical buckling load
which can be achieved for a beam-type element with length L, by defining a variation
function for the moment of inertia. To define the variation function of the moment of
inertia, five coefficients are introduced which are associated with five equidistant
points on the element, which indicate the variation of the moment of inertia as against
the reference cross- section. The reference cross-section’s geometric and inertial
characteristics are computed based on the imposed failure condition. At these points,
the quantity of material is not modified, only its distribution on the cross-section.
Between these reference points, cubic interpolation is performed in order to obtain a
continuous function to express the variation of the moment of inertia.

35
The second iterative cycle starts after fixing the values for the five coefficients, which
describe the variation of the moment of inertia, and the geometric properties of the
reference cross-section. In the case of the 5 critical sections, it is important that the
class of the sections be less than 4, in order to avoid the local verification process as
described in SREN 1993-1-6. The third iterative cycle assumes verifying the failure
condition, which presumes comparing the critical buckling load with the design
resistance of the cross-section for uniform compression.
By undergoing the presented three iterative cycles, it is possible, with respect to the
modelling assumptions, to fix the geometric characteristics for the beam-type element,
with length L and variable cross-section, in order to achieve a favorable failure
condition and avoid the loss of stability of the element

Figure 6. Flow chart of the work method

36
.

The abbreviations in Figure 6 are related to a series of subprograms that incorporate

the theoretical aspects from Section 3, as follows.

The subprogram entitled model_section_var sets the cubic interpolation function

that describes the variation of the moment of inertia, as presented in Section 3. The
input arguments are the length of the element L and the cross-section variation
coefficients vector [i1 i2 i3 i4 i5]. The output for this program is the symbolic function
Ii, which describes the modification of the moment of inertia.

The subprogram entitled sbc establishes the critical axial load, identified as the
Rayleigh quotient, based on the variation of the cross-section, the material type, and
the kinematic restraints imposed on the beam-type element. The input arguments are
the length of the element L, the Young modulus E, the symbolic function Ii, and the
stiffness vector for the joints. The outputs are the critical buckling force RQ and the
deflection shape of the element through the displacement vector fp, which are the
solution to the eigenvalue problem.

The subprogram entitled check_section is a verification program that calculates the

geometric and inertial characteristics of the 5 critical sections and indicates if the class
of the cross-section is greater than 3, based on the geometric characteristics of the
reference cross section. The input arguments are the cross-section variation
coefficients vector, the outside diameter of the reference cross-section D0 and the
thickness of the wall of the reference cross- section t0. The output of the program is a
logical one. If the answer is TRUE then the variable cross-section beam-type element
has no issues. When the answer is FALSE then modifications are required regarding
the geometric entries of the reference cross-section or the variation coefficients of the
moment of inertia.

37
The subprogram entitled failure_condition checks if the imposed failure condition is
satisfied, regarding the failure due to the resistance of the cross-section for uniform
compression rather than buckling. According to Hooke’s law, conditions can be
expressed in terms of the area of the cross-section. The input arguments are the area of
the reference cross- section A0, the steel-yield strength fy, and the critical buckling load
RQ. The outputs are the elastic resistance of the cross-section Fc,Rd and the critical
buckling load Fb,Rd.

Aspects Regarding the Input Data

Aspects Regarding the Length of the Beam-Type Element

The starting point of the modelling process is fixing the length of the beam-type
element, L. The length of the element is a numerical value that remains constant
during the iteration cycles and the kinematic constraints imposed on the joints.

In order to obtain a beam-type element with optimal shape from a topological point of
view, the criteria regarding the variation of the cross-section are a necessary but not
sufficient condition. The additional criteria are in close relation to Euler’s critical load,
presented in Expression (1). In Expression (1), the second important factor which
influences the critical buckling load is the length of the element, a numerical value
that must be correlated with the geometric characteristics of the cross-section. Given
the length of the element, a factor of inverse proportionality, there are three possible
situations to note.

The possible situations are related to the effective geometric length of the element L
and the critical length of the element 𝐿𝑐𝑟𝐿𝑐𝑟, which is an ideal length for which the
failure condition is fulfilled and the modelling of a variable cross-section makes sense.

38
In case I, the length of the element is less than the ideal length, or 𝐿<𝐿𝑐𝑟𝐿<𝐿𝑐𝑟. In this
case, the failure will occur by exceeding the value of the design resistance of the cross-
section for uniform compression. The critical buckling force is net superior. In the
present case, the modelling of a variable cross-section makes no sense, given the
computational effort, the stiffness of the joints, and the technological terms. In this
situation, the constant cross-section is an ideal solution.

In case II, the length of the element is quasi-equal to the ideal length, or 𝐿=𝐿𝑐𝑟𝐿=𝐿𝑐𝑟.
This case is considered the ideal case. The failure will occur by exceeding the value
of the design resistance of the cross-section for uniform compression. The critical
buckling force is superior. Modelling a variable cross-section is necessary in order to
achieve a superior value for the critical buckling load beyond the design resistance of
the cross-section for uniform compression.

In case III, the length of the element is greater than the ideal length, or 𝐿>𝐿𝑐𝑟𝐿>𝐿𝑐𝑟.
The failure will occur due to the loss of stability. The material is not used rationally.

For the cases for which the length of the element is a numerical value in the
neighborhood of 𝐿𝑐𝑟𝐿𝑐𝑟, there are several optimization methods, one of which is
proposed in this study is to use a minimum quantity of material that has an optimum
distribution on the cross-section

39
Material Type and Section Profile
Slender structural elements subjected to a significant amount of compression
forces appear to be characteristic of steel structural assemblies. In the case of steel
structures, the beam-type elements are modeled according to the main stress in the
cross-section. In the case of high axial forces, the most frequently used and
recommended profiles are circular tubular profiles, with s-planes of sectional
symmetry, with the inertial properties are equal in every direction.

Given the boundary conditions for the beam-type element, circular tubular profiles
are not sensible to flexural buckling. The recommended steel circular tubular
profiles are presented in EN10219. In this study, these profiles will be considered
as profiles for the reference cross-sections. The profile types presented in EN10219
are manufactured for constant cross-section elements. In this study, the reference
cross- section id is presented in Expression (63), which indicates the geometric
properties and material type for the cross-section.

CHS Dxt − L − S355JO − EN1021 (63)

In Expression (63) CHS is the abbreviation for the circular hole section, the type of
profile chosen, D represents the external diameter in [mm], t the thickness of the
profile’s wall in [mm], L represents the length of the beam-type element in [mm],
S355JO is the grade of the steel, which carries a minimum yield strength of 355
[N/mm2], JO indicates that the steel has undergone longitudinal Charpy V-Notch
impact testing at 27[J] at 0 [°C], EN10219 is the reference standard for the profile
type].

The Stiffness of the Support Connections

In a theoretical study for the hindered displacements in a support connection, infinite
value is associated with the stiffness of the connection. In the MATLAB
programming language, the concept of infinity is associated with a numerical value
of 101000 .To quantify the effects of real support connections which allow small
displacements there are imposed finite numerical values for the stiffnesses of the
hindered displacements. In the case of other not constrained connections, their
stiffness is not zero.

40
 In this study, the numerical values imposed on the stiffnesses of the
support connections are purely theoretical and they do not result from
measurements on a physical model. In the theoretical model, the support
connections are considered absolutely rigid or absolutely flexible. For any
intermediate values of the rigidity of the elastic support, a flexible joint is
obtained .The modelling of elastic support connections resumes fixing the
value of the stiffness for every connection type. The numerical values are
computed based on the stiffness properties of the cross-section at the joints
for the beam-type elements. For example, in the case of a beam-type
element that is clamped-clamped the expression for the support connection
stiffnesses, if one considers only the effect of bending, is presented in
Expression (64) .

12 ∗ 𝐸 ∗ 𝐼 6 ∗ 𝐸 ∗ 𝐼 12 ∗ 𝐸 ∗ 𝐼 6 ∗ 𝐸 ∗ 𝐼 sn
sn sn sn
𝑟=( ) (64)
3 2 3 𝐿2
𝐿 𝐿 𝐿
∗ ∗ ∗

41
In Expression (64), E represents young’s modulus, L represents the length of the
element, 𝑠𝑛 denotes the moment of inertia at the joint cross-section. In Expression

(64) for the first joint, it corresponds to the first and second values, the stiffness
for a hindered transversal displacement and a hindered rotation of the joint, and
for the second joint corresponds to the third and fourth value, the stiffness for a
hindered transversal displacement and a hindered rotation of the joint. For every
beam-type element, it is possible to fix a reference cross-section, which is
considered the geometric invariant of the iterative process .

Aspects Regarding the Variation of the Cross-Section

The cubic interpolant describes the variation of the cross-section of the element on m
sub-intervals. It is necessary to fix the values for the interpolation function at the
limit of the sub-intervals, meaning m+1 numerical values, indicated through
coefficients

𝑖𝟏 ≤ 𝒌 ≤ 𝒎 + 𝟏 = [1: 2] .

Due to rational considerations, the maximum modification for the moment of

inertia of the cross-section is 100 [%] concerning the reference cross-section. The
final values for the coefficients are fixed at the end of the iterative process. The
modification of the moment of inertia assumes the modification of the outside

diameter of the cross-section 𝐷𝟏 ≤ 𝒌 and the wall thickness of the cross-

section 𝑡𝟏 ≤ 𝒌 ≤ 𝒎 + 𝟏. ≤ 𝒎+𝟏

In this study, the interpolation function, which describes the variation of the
moment of inertia, is modeled on 4 sub-intervals. Therefore, it is necessary to fix the
value at five reference points, according to Expression (65).

𝑖𝒗 = [𝑖𝟏 𝑖𝟐 𝑖𝟑 𝑖𝟒 𝑖𝟓] 𝑖 ∈ [1: 2] (65)

42
Aspects Regarding Subprogram Check Section
The potential points where it is recommended to modify the value of the
moment of inertia are indicated by the defection shape of the compressed
element. At the limits of m + 1 sub-intervals, the values of the moment of
inertia get local maximum values. Subprogram Check section determines
the outside diameter of the cross-section 𝐷𝟏 ≤

𝒌 ≤ and the wall thickness of the cross-section 𝑡𝟏 as a function of the

≤ 𝒌 ≤ 𝒎 + 𝟏.
𝒎+𝟏

reference cross-section. At the limit of a sub-interval, m Expression (66)

can be written.

𝐼𝟏 ≤ 𝒌 ≤ 𝒎 + 𝟏 = 𝑖𝟏 ≤ 𝒌 ≤ 𝒎 + 𝟏 ∗ 𝐼𝑜 (66)

In Expression (66), 𝐼0 denotes the numerical value of the reference cross-

section, which is fixed for the constant cross-section beam-type element
according to the imposed failure condition. The numerical value of
parameter i is correlated with the amplitude of the deflection shape with
respect to the kinematic boundary conditions.
Modelling the cross-section is based on Cavalieri’s principle. The
area of the cross- section is equal for every cross-section, presented in
Relation (67).

𝐴𝟏 ≤ 𝒌 ≤ 𝒎 + 𝟏 = 𝐴𝟎 = 𝑐𝑜𝑛𝑠𝑡. (67)

43
 𝐴𝟏 ≤ 𝒌 ≤ 𝒎 + 𝟏 = 𝐴𝟎 = 𝑐𝑜𝑛𝑠𝑡. (67)

In Relation (67), 𝐴0 represents the area of the reference cross-section, and it is

highlighted the proposed modelling principle of the redistribution of material in the
cross-section, highlighted also in Figure7. The numerical value for every cross-section
is fixed for the constant cross- section beam-type element according to the imposed
failure condition.

Figure 7. Variable cross-section beam-type element, where L—is the

length of the beam- type element; x—marks the distance to a current cross-
section; X—is the local reference system axes along the axes of the beam-type
element; y—is the local reference system axes perpendicular to the beam-type
element axes; I—the moment of inertia of the cross-section in a general form;
I(x)—the moment of inertia of the cross-section at section x; Dk—the outside
diameter of the transversal cross-section k; dk—the inside diameter of the
transversal cross- section k; tk—the thickness of the wall of transversal cross-
section k; 0—the center of cross- section k-k.

In Figure 7 it is shown that a longitudinal and transversal cross-section of a variable

cross-section beam-type element. Cross-section 1-1 corresponds to a section at the first
joint, and cross-section 2-2 represents a current cross-section where a considerable
modification of the moment of inertia is recorded. For cross-section 1-1 and cross-
section 2-2 one may define the geometric and inertial characteristics presented in
Expression (68a)–(68d).

44
 𝜋
𝐴𝟏 = ∗ (𝐷𝟏2 − 𝑑𝟏2 ) (68a)
4
𝜋
𝐴𝟐 = 4 ∗ (𝐷𝟐2 − 𝑑𝟐2 ) (68b)

𝜋
𝐼𝟏 = ∗ (𝐷𝟏2 − 𝑑𝟏2 ) (68c)
64
𝜋
𝐼𝟐 = 64 ∗ (𝐷𝟐2 − 𝑑𝟐2 ) (68d)

In Expression (67), 𝐴𝑘 represents the area of the cross-section; 𝐼𝑘 denotes the moment of
inertia of the cross-section; 𝐷𝑘 represents the outside diameter and 𝑑𝑘 represents the
inside diameter of the cross-section, for ∀𝑘=1:2. By substituting Expressions (68a)–
(68d) into Expression (67) the system of Equations (69) is obtained, which can be
expressed as in Equation (70) due to a computational artifice.

By substituting Equation (2) in Equation (1) from the system (70) and with Equation
(2) from the System (71), is obtained.

If Equation (1) is gathered with Equation (2), from the system (71), then Equation
(72) can be written, from which the outside diameter of cross-section 2-2 can be
obtained as a function of the geometric characteristics of cross-section 1-1,
concerning the modification of the moment of inertia, presented in Expression (73).

45
if Equation (2) is subtracted from Equation (1), from the system (71), Equation (74)
can be written, from which the inside diameter of cross-section 2-2 can be obtained as
a function of the geometric characteristics of cross-section 1-1, concerning the
modification of the moment of inertia, presented in Expression (75).

2·𝑑22 =(𝑖2−1)·𝐷12+(𝑖2+1)·𝑑12 (74)

(𝑖2 −1).𝐷12 +(𝑖2 +1).𝑑12

𝑑2 = √ (75)
2

If a ratio, p, is accepted between the outside and inside diameter of the cross-
section for cross-section 1-1 Relation (76) can be written.
= 𝑝 ≥1 (76) 𝐷𝑑11

It is admitted the same ratio, p, between the outside and inside diameter of the
cross-section for corss-section 1-1 Relation (77) can be written.
= 𝑝≥1 (77) 𝐷𝑑22

Expressions (78) and (79) are obtained by squaring Relation (76) and substituting the outside
diameter for cross-section 1-1 in Expressions (73) and (75).
𝑖 .(𝑝2+1)+(𝑝2 −1)
𝑑2 = 𝑑1 . √2 (78)
2

𝑖 .(𝑝2+1)−(𝑝2 −1)
𝑑2 = 𝑑1 . √2 (79)
2

If Relation (77) is squared and Expressions (78) and (79) are substituted, Equation
(80) isobtained.
−𝑖2·(𝑝2+1)·(𝑝2−1)+(𝑝2+1)·(𝑝2−1) = 0 (80)

46
The solution to Equation (80), with respect 𝑝≠1, highlights that the assumed ratio for
both of the cross-sections is valid if the beam-type element has a constant cross-section.
According to the reductio ad absurdum principlee, it is not possible to associate the
same ratio for both of the cross- sections, therefore Relations (81) and (82) are written.

= 𝑝≥1 (82) 𝐷𝑑2

2

= 𝑝 ≥1 (81) 𝐷𝑑11

At this phase, the task is to determine ratio p2 as a function of ratio p1. Therefore by
dividing Expression (78) with Expression (79) and based on Relation (82) the
expression of ratio p2 can be expressed as in Expression (83).

𝑖2 .(𝑝12+1)+(𝑝12 −1)
𝑝2 = √ 2 2 (83)
𝑖 2 .(𝑝1 +1)−(𝑝1 −1)

47
To be able to control the modified cross-sections, the thickness of the profile wall of
cross- section 2-2 is expressed as a function of the thickness of the profile wall of cross-
section 1-1. Cross-section 1-1′s thickness of the profile wall can be expressed as in
Equation (84).

𝐷1−𝑑1=2·𝑡1 (84)

From Equation (84) the inside diameter of

section1–1, 𝑑1, can be expressed, which is
introduced in Expressions (78) and (79) one may obtain Expressions (85) and (86).

2 2
𝐷2 =
2.𝑡1
√𝑖2 .(𝑝1 +1)+(𝑝1 −1) (85)
(𝑝1 −1) 2

2 2
𝑑2 =
2.𝑡1
√𝑖2 .(𝑝1 +1)−(𝑝1 −1) (86)
(𝑝1 −1) 2

If Expression (86) is substracted from Expression (85) the expression for the thickness
of the profile’s wall for cross-section 2-2 is obtained, as a function of the thickness of
cross-section 1-1′s profile wall and ratio p1, as in Expression (87).

𝑡1 𝑖 .(𝑝12 +1)+(𝑝12 −1) 𝑖 .(𝑝12 +1)−(𝑝12−1)

𝑡2 = .√ 2 −√2 (87)
(𝑝12 −1) 2 2
Based on Expressions (85) and (87), the new sections can be controlled
preliminary during

the modelling phase, in order to define the class of the modified cross-section. By
modifying the moment of inertia of the cross-section and by keeping constant the area,
based on Cavalieri’s principle, the thickness of the profile’s wall will thin out. The
iterative process must keep in mind this phenomenon, to avoid slender sections and to
keep the class of the cross-section under 4, as to avoid further verification as presented
in SREN 1993-1-6.

48
By representing the buckling stress as function of slenderness as in Figure 8, the
target zone in this study is the shaded zone. The classification of the cross-section was
performed as presented in SREN1993-1-1.

Figure 8. Buckling stress as a function of slenderness, where λ-denotes the

slenderness of the beam-type element; fy - the steel-yield strength; σ𝐸σ𝐸-the
buckling stress; σλσ𝜆—the buckling stress as a function of the slenderness; E-
Young’s modulus.

Aspects Regarding the Failure Condition

The optimization process assumes defining the configuration of the beam-type

element through the moment of inertia of the cross-section in order to achieve a
failure due to exceeding the design resistance of the cross-section under uniform
compression rather than the design buckling resistance of the compression member. A
limit state is characterized by appropriate values of the design resistances, which is a
clear objective of the optimization process. The mentioned condition is expressed in
a mathematical form as in Relation (88)

1.05·𝐹𝑐Rd= FbRd (88)

49
In Relation (88), the symbols as in EN1993-1-1, 𝐹𝑐,𝐹𝑐,𝑅𝑑 represent the design
resistance of the cross-section under uniform compression and 𝐹𝑏,𝑅𝑑𝐹𝑏,𝑅𝑑
represent the design buckling resistance of the compression member. The margin of
error is considered 5 [%], and the empirical value, is quantified in Relation (88) by
coefficient 1.05

Estimating the Error of the Iterative Process

The convergence criteria of the iterative process is to compute a numerical value for
the critical buckling load which exceeds the value of the resistance of the cross-section
under uniform compression, by modelling a variable cross-section for the beam-type
element. To highlight the precision of the proposed method, a parameter is introduced
to verify if the displacement field defined based on the defection shape of the beam-
type element through the proposed method describes neutral equilibrium at critical
load.
The definition of the control parameter, Error, is conditioned by the assumption
that at the state of neutral equilibrium at critical load, the variation of the total
potential energy iszero. The displacement field which describes the equilibrium
state is approximated through a polynomial function. Therefore the variation of the
total potential energy is not equal to zero. To quantify this amount of energy, it is
related to the algebraic sum of the strain energy and the energy stored in the elastic
support connection due to deformation. The mathematical form of the control
parameter is presented in Expression (89).
𝐿𝑒𝑓 −𝑄.𝐿𝑒𝑥𝑡 (𝑄)
𝐸𝑟𝑟𝑜𝑟 = (89)
𝐿𝑒𝑓 +𝐿𝑟𝑒𝑎𝑧𝑒𝑚

50
Determining the Degree of the Polynomial Function v(x)
The displacement field is approximated through a polynomial function of degree n − 1.

To fix the degree of the polynomial function, an analysis is performed on an axially

loaded Euler column, a joint-supported statically determined system, for which the
Rayleigh quotient and the Timoshenko quotient are calculated. This analysis is
performed based on the first iteration cycle, with input values for the length of the
element, Young’s modulus, and the moment of inertia of the reference cross-section
unitary.
Based on the two methods of approximating the critical buckling loads, the analysis
assumes the variation of the polynomial function’s degree until the numerical values
for thecritical buckling load are equal. The analysis is performed on an Euler column
with a constant cross-section with results presented in Table 2 and on an Euler
column with variable cross- section, with the coefficients which describe the
variation of the moment of
inertia 𝑖𝑣=[1 1.9 2 1.9 1]𝑖𝑣=1 1.9 2 1.9 1, results presented in Table 3.

Table 2. The results of the analysis for the constant cross-section Euler
column

51
the case of the constant cross-section Euler column, the maximum error is 0.519%. The
numerical values in Table 2 are presented in a graphical form in Figure 9.

Figure 9. The variation of the critical buckling force for the constant cross-section
Euler column.

Based on the data from Table 2, for the constant cross-section Euler column, it is
sufficient to model the displacement field through a polynomial function of degree
eight.
In the case of the variable cross-section beam-type element, the maximum error is 0.573%.

The numerical values in Table 3 are presented in a graphical form in Figure 10.

Figure 10. The variation of the critical buckling force for a

variable cross-section Euler column.

52
Aspects Regarding the Modelling Phase concerning the Rayleigh Quotient
The approximation of the critical buckling load, through the Rayleigh
quotient, presumes an iterative process through it is assumed a particular
deflection shape, then after h iteration cycles, the critical buckling load is
approximated concerning the modelling assumptions. The number of
iteration cycles h is established based on the number of deflection shapes
taken into account, through the moment of inertia of the cross-section.

In this study, it is assumed a particular variation of the moment of

inertia of the cross- section, approximates the deflection shape of the
beam-type element. By modifying the
coefficients which describe the variation of the moment of inertia it is
obtained a new defection shape, which is more appropriate to the real
deflection shape, therefore the Rayleigh quotient decreases to a value more
appropriate to the critical buckling load.
With the modification of the shape of the beam-type element, a new
displacement field is obtained, and also the elastic stiffness matrix is
changing therefore the critical buckling force increases. This is the basis of
the first iterative cycle.

53
5. Case Study/Parametric Study/Computational Examples
The computational examples assume the study of statically determinate and
indeterminate beam- type elements, six situations with ideal support connections,
and also the consequence of the elastic support connections.

The starting point of every individual study is fixing the „reference” beam-type
element, with a constant cross-section. The given element is modeled to achieve a
critical buckling load superior to the resistance of the cross-section at uniform
compression, explicitly to avoid the loss of stability.

The modelling of the beam-type element assumes knowing the kinematic boundary conditions.
The kinematic boundary conditions in the case of ideal support connections can be
defined as a zero displacement field. Regardless of the imposed kinematic boundary
conditions the beam-type elements are considered with an equal geometric length of
6000 [mm], Young’s modulus of 2.1·105 [N/mm2], and yield strength of the steel, fy,
with the value of 355 [N/mm2].

The hinged-hinged beam-type element

The hinged-hinged beam-type element is presented in Figure 11a, which highlights the
kinematic boundary conditions. The type of profile imposed on the reference cross-
section is CHS163.8x8–6000– S355JO–EN10219, respectively values for the stiffness
of the support connections are [inf 0 inf 0]. The results of the iterative optimization
phase are presented in Table 4 and the graphical representation of the data from Table
4 is presented in Figure 12.

54
Figure 11. The hinged-hinged beam-type element: (a) Schematic
representation of the kinematic boundary conditions, where v—denotes the
transversal displacement; Q—the axial force; (b) The variable cross-section
beam-type element concering the imposed boundary conditions; (c)
Longitudinal section of the modeled beam-type element which highlights the
variation of the cross-section.

Figure 12. The critical buckling force variation as

against the resistance of the cross-section to uniform
compression based on the data from Table 4.

55
In Table 4 as a function of the variation of the moment of inertia the resistance of
the cross- section for uniform compression, Fc,Rd, and the buckling resistance,
Fb,Rd, is computed and in the lastcolumn is presented the estimated error of the
iteration cycle. The modeled hinged-hinged beam-type element is presented in
Figure 11b,c.In the case of the hinged-hinged beam-type element through the
variation of the cross-section, without material addition, a significant increase of
93.23% for the critical buckling load is obtained and ensuring a proper failure
condition as presented in Figure 12.

The clamped-free beam-type element

The clamped-free beam-type element is presented in Figure 13a, which
highlights the kinematic boundary conditions. The type of profile imposed
on the reference cross-section is CHS323.9x12.5– 6000–S355JO–
EN10219, respectively values for the stiffness of the support connections
are [inf inf 0 0]. The results of the iterative optimization phase are presented
in Table 4 and the graphical representation of the data from Table 5 is
presented in Figure 14

56
Figure 13. The clamped-free beam-type element: (a) Schematic
representation of the kinematic boundary conditions, where v—denotes
the transversal displacement; Q—the axial force; (b) The variable cross-
section beam-type element concerning the imposed boundary conditions;
(c) Longitudinal section of the modeled beam-type element which
highlights the variation of the cross- section.

57
Figure 14. The critical buckling force variation against the resistance of the cross-
section to uniform compression based on the data from Table 5.

.
In Table 5 as a function of the variation of the moment of inertia the resistance of
the cross- section for uniform compression, Fc,Rd, and the buckling resistance,
Fb,Rd, is computed and in the last column is presented the estimated error of the
iteration cycle. The modeled clamped-free beam-type element is presented in
Figure 13b,c.

In the case of the clamped-free beam-type element through the variation of the
cross-section, without material addition, a significant increase of 95.79% for the
critical buckling load is obtained and ensuring a proper failure condition as
presented in Figure 14.

The guided-hinged beam-type element

The guided-hinged beam-type element is presented in Figure 15a, which highlights
the kinematic boundary conditions. The type of profile imposed on the reference
cross-section is CHS323.9x12.5– 6000–S355JO–EN10219, respectively values for
the stiffness of the support connections are [0 inf inf 0]. The results of the iterative
optimization phase are presented in Table 6 and the graphical representation of the
data from Table 6 is presented in Figure 16.

58
Figure 15. The guided-hinged beam-type element: (a) Schematic
representation of the kinematic boundary conditions, where v—denotes
the transversal displacement; Q—the axial force; (b) The variable cross-
section beam-type element concerning the imposed boundary conditions;
(c) Longitudinal section of the modeled beam-type element which
highlights the variation of the cross- section.

Figure 16. The critical buckling force variation as

against the resistance of the cross-section to uniform
compression based on the data from Table 6.
59
In Table 6 as a function of the variation of the moment of inertia the resistance of
the cross- section for uniform compression, Fc,Rd, and the buckling resistance,
Fb,Rd, is computed and in the lastcolumn is presented the estimated error of the
iteration cycle. The modeled guided-hinged beam-type element is presented in
Figure 15b,c.
In the case of the guided-hinged beam-type element through the variation of the
cross-section, without material addition, a significant increase of 94.81% for the
critical buckling load is obtained and ensuring a proper failure condition as
presented in Figure 16.

The clamped-guided beam-type element

The clamped-guided beam-type element is presented in Figure 17a, which
highlights the kinematic boundary conditions. The type of profile imposed on the
reference cross-section is CHS193.7x8–6000–S355JO–EN10219, respectively
values for the stiffness of the support connections are [inf inf 0 inf]. The results of
the iterative optimization phase are presented in Table 7 and the graphical
representation of the data from Table 7 is presented in Figure 18.

60
Figure 17. The clamped-guided beam-type element: (a) Schematic
representation of the kinematic boundary conditions, where v—denotes
the transversal displacement; Q—the axial force; (b) The variable cross-
section beam-type element concerning the imposed boundary conditions;
(c) Longitudinal section of the modeled beam-type element which
highlights the variation of the cross- section.

Figure 18. The critical buckling force variation as against the resistance of the
cross-section to uniform compression based on the data from Table 7.

61
In Table 7 as a function of the variation of the moment of inertia the resistance of
the cross- section for uniform compression, Fc,Rd, and the buckling resistance,
Fb,Rd, is computed and in the last column is presented the estimated error of the
iteration cycle. The modeled clamped-guided beam-type element is presented in
Figure 17b,c. In the case of the clamped-guided beam-type element through the
variation of the cross-section, without material addition, a significant increase of
46.91% for the critical buckling load is obtained and ensuring a proper failure
condition as presented in Figure 18.

The hinged-clamped beam-type element

The hinged-clamped beam-type element is presented in Figure 19a, which highlights the
kinematic boundary conditions. The type of profile imposed on the reference cross-section is
CHS127x5–6000–S355JO–EN10219, respectively values for the stiffness of the support
connections are [inf 0 inf inf]. The results of the iterative optimization phase are presented in
Table 8 and the graphical representation of the data from Table 8 is presented in Figure 20.

62
Figure 19. The hinged-clamped beam-type element: (a) Schematic

representation of the kinematic boundary conditions, where v—denotes

the transversal displacement; Q—the axial force; (b) The variable cross-
section beam-type element concerning the imposed boundary conditions;
(c) Longitudinal section of the modeled beam-type element which
highlights the variation of the cross- section.

Figure 20. The critical buckling force variation as against the

resistance of the cross-section to uniform compression based on the
data from Table 7.
63
Figure 20. The critical buckling force variation as against the resistance of the cross-
section to uniform compression based on the data from Table 7.

In Table 8 as a function of the variation of the moment of inertia the resistance of

the cross- section for uniform compression, Fc,Rd, and the buckling resistance,
Fb,Rd, is computed and in the last column is presented the estimated error of the
iteration cycle. The modeled hinged-clamped beam-type element is presented in
Figure 19b,c.

In the case of the hinged-clamped beam-type element through the variation of the
cross-section, without material addition, a significant increase of 62.13% for the
critical buckling load is obtained and ensuring a proper failure condition as
presented in Figure 20.

The clamped-clamped beam-type element

The clamped-clamped beam-type element is presented in Figure 21a, which
highlights the kinematic boundary conditions. The type of profile imposed on the
reference cross-section is CHS101.6x6–6000–S355JO–EN10219, respectively
values for the stiffness of the support connections are [inf inf inf inf]. The results of
the iterative optimization phase are presented in Table 9 and the graphical
representation of the data from Table 9 is presented in Figure 22.

64
Figure 21. The clamped-clamped beam-type element: (a) Schematic
representation of the kinematic boundary conditions, where v—denotes the
transversal displacement; Q—the axial force; (b) The variable cross-section
beam-type element concerning the imposed boundary conditions; (c)
Longitudinal section of the modeled beam-type element which highlights
the variation of the cross- section.

Figure 22. The critical buckling force variation as against the resistance of
the cross-section to uniform compression based on the data from Table 9.

65
In Table 9 as a function of the variation of the moment of inertia the resistance of
the cross- section for uniform compression, Fc,Rd, and the buckling resistance,
Fb,Rd, is computed and in the last column is presented the estimated error of the
iteration cycle. The modeled clamped-clamped beam- type element is presented in
Figure 21b,c.In the case of the clamped-clamped beam-type element through the
variation of the cross-section, without material addition, a significant increase of
28.25% for the critical buckling load is obtained and ensuring a proper failure
condition as presented in Figure 22.

The clamped-clamped beam-type element

The beam-type element with elastic supports is presented in Figure 23a, which
highlights the kinematic boundary conditions. The type of profile imposed to the
reference cross-section is CHS101.6x6–6000–S355JO–EN10219, respectively values
for the stiffness of the support connections are [0.002 [kN/m] 28.935 [kNm/m] 0.002
[kN/m] 28.935 [kNm/m]]·104 established according to Expression (64). The results of
the iterative optimization phase are presented in Table 10 for ideal and elastic
restraints.

66
Figure 23. The beam-type element with elastic supports: (a) Schematic
representation of the kinematic boundary conditions, where L—is the length
of the beam-type element; Q—is the axial force; r1— denotes the stiffness
for a hindered transversal displacement for joint i; r2—denotes the stiffness
for a hindered rotation for joint i; r3—denotes the stiffness for a hindered
transversal displacement for
joint j; r4—denotes the stiffness for a hindered rotation for joint j; (b) The
variable cross-section beam- type element concerning the imposed boundary
conditions; (c) Longitudinal section of the modeled beam-type element
which highlights the variation of the cross-section.

67
In Table 10, the symbol fin represents the finite numerical value for the corresponding stiffness
of the elastic support, determined according to Expression (64). The symbol inf represents the
ideal beam-type element support connections stiffness. In Table 10 it is highlighted that
considering finite numerical values for the stiffness of the support connections the critical
buckling load decreases radically, for example in case 1b the buckling force decreases by
84.07% and in case 2b buckling force decreases with 55.31%.

To avoid such sudden changes, which influence negatively the failure condition,
another type of cross-section is chosen, with superior geometric and inertial
characteristics. In Table 11 the data from modelling a beam-type element with elastic
support is presented with the reference cross-section CHS139.7x6–6000–S355JO–
EN10219, and the numerical values for the stiffness of the support connections
[0.007 [kN/m] 78.996 [kNm/m] 0.007 [kN/m] 78.996 [kNm/m]]·104 determined
according to Expression (64).

In Table 11 as a function of the variation of the moment of inertia the resistance of

the cross- section for uniform compression, Fc,Rd, and the buckling resistance,
Fb,Rd, is computed and in the last column is presented the estimated error of the
iteration cycle.

In the case of the beam-type element with elastic supports through the modification of
the cross- section, without material addition, a significant increase of 1725.05% for
the critical buckling load is obtained as against the new reference cross-section beam-
type element, and it is ensured a propper failure condition as presented in Figure 24.
Figure 23b,c are presented the beam-type element with elastic supports. In Table 12 a
comparison is presented of the critical buckling forces without and with elastic
support, with the variation coefficients of the moment of inertial extracted from Table
10.

68
Figure 24. The critical buckling force variation as against the resistance of
the cross-section to uniform compression with respect to the data from Table
11.
.

Case 1b from Table 12 highlights the fact that by adding elastic supports
the failure occurs due to the loss of stability and variation modelling a
variable cross-section for the beam-type element is not
appropriate and the reference cross-section with superior geometric and inertial
properties is required.
Case 2b from Table 12 highlights that by adding elastic supports for a reference
cross-section
with superior geometric and inertial properties the failure occurs due to

69
exceeding the resistance of the cross-section under uniform compression and variation
modelling a variable cross-section for the beam-type element is necessary. Adopting
elastic support instead of ideal support connections at least a 55% decrease is
noticed in case of the critical buckling force.

The explanation for the radical decrease of the critical buckling force in the case of elastic
support connections may be found in an analysis of the displacement field of a
constant cross-section beam-type element loaded with axial forces and bending
moments at the joints. In the case of the hinged-hinged beam-type element the
displacements associated with the bending moment M, are proportional with M·L2
divided by the stiffness coefficient E·I, expressed in Relation (90a). The
displacements associated with the axial force Q, are proportional with Q·L divided
by the stiffness coefficient E·A, presented in Relation (90b).

The quantities expressed in Relation (90a,b) have the same order of magnitude if
conditions from Relation (91a,b) are satisfied, with respect to Brachmann-Landau
notation .

70
In Relation (96) is highlighted the influence of strain associated with the bending moment on
the displacement field of the beam-type element, with in this case is practically insignificant. In
this manner by adding elastic support connections with finite values for the stiffnesses,
computed based on the Expression (64), regardless of the imposed boundary conditions, the
beam-type element’s behavior is practically identical to the Euler column.

The method through the ideal support connection is replaced with elastic support connection,
with finite numerical values for the stiffness coefficients computed based on the
deformability of the structural assembly, must consider the fact that the numerical
value of the critical buckling load can be half as computed. To better understand the
phenomena, one can make an analogy with a planer truss with rigid nodes and
structural assemblies with high statical indeterminacy, due to the axial strain of the
structural elements, in the joints of the beam-type elements bending moments will
occur as a function of the relative rotation of the trusses, an effect which has been
proven to be a secondary one.

71
Conclusions
Aspects Referring to the Precision of the Proposed Method
The first important aspect which needs to be highlighted is the high precision of
the proposed method, which approximates the critical buckling load for the case of
a beam-type element with a variable cross-section. The high precision is in the first
place to the adopted methodology, more precisely by adopting a cubic interpolant
to express the variation of the cross-section and by modelling the displacement
field through a polynomial function with a degree greater than five. Secondly, it is
possible to achieve such precision due to the numerical properties of an irrational
number, which is closely related to Expression (1), Euler’s critical buckling load.
In Expression (1) the irrational number ππ appears which decimal representation
never ends, respectively its decimals do not enter a permanently repeating pattern.
This property makes it possible to estimate with high accuracy the variation of the
total potential energy at neutral equilibrium at critical load.

The Influence of the Statical Determinacy of the Structural Element

The case study represents a numerical but especially analytical analysis, which
is a different kind of approach for stability analysis. It assumes linking two
theories, namely Cavalieri’s principle and elastic buckling theory, in order to
perform a comprehensive stability analysis in the elastic domain to increase the
elastic buckling load for the beam-type element, through geometric modelling of
the cross-section for the beam-type elements, with specific boundary conditions.

The case study ranges from statically determinate to statically indeterminate

structural elements, concerning the imposed kinematic boundary conditions. In the
case of statically indeterminate structural assemblies, it is not possible to establish
a unique numerical value for the stress ratios of the elements, due to the fact that
modifying the axial stiffness of an element influences the stress redistribution in
the whole structure. In the case of statically determined structural assemblies, it is
possible to compute a unique value for the unitary efforts in the structural
elements. The case study has proven that in the case of statically determinate
structural elements have achieved the highest increase of the critical buckling load.
The final result in therm of the achieved increase rates is presented in Figure 25.

72
 A very important result is highlighted in Figure 25 for the beam-type element
with elastic supports. The elastic support connections modify significantly the
behavior of the structural element under axial loads and influence excessively the
critical buckling load. In the case of constant cross-sections increasing the area of
the cross-section is not a rational solution and a favorable failure condition is
practically impossible to achieve. Modelling a variable cross-section beam-type
element is a practical solution to increase the critical buckling load above the value
of the resistance of the cross-section subjected to uniform compression. A correct
modelling of the shape of the beam-type element, through the moment of inertia,
may point to a considerable increase of the critical buckling load, up to 1725% in
the present case as against the constant cross-section, with respect to Cavalieri’s
principle. A very important aspect is that by imposing elastic support connections
instead of ideal ones the critical buckling load is influenced and a significant
decrease is observed, up to at least 55% in the studied case. This remark is very
important for structural engineers when they establish the geometric and inertial
characteristics of a structural element subjected to considerable axial load.

A certain fact is that statically determinate structures have lack of structural

redundancy, but the state of stress and strain can be easily controlled through the
whole elastic behavior domain of the beam-type element. By adopting beam-type
elements with variable cross-sections it is possible to fix significant numerical
values for the length of the elements, with respect to the secondary effects of the
rigid joints.

73
 Based on the presented aspects the adopted beam-type elements with variable
cross-sections are recommended in the case of high-span planar or spatial trusses.
The mentioned structural assemblies must be statically determinate, made up of a
minimum number of elements, with significant length, and using high-quality
materials.

Aspects Regarding the Common Undesired Failure Conditions

The failure of a structural assembly and in essence structural elements must
take place with a high amount of warning signs which enables the evacuation of
persons and material goods from the potentially affected area. In this manner, there
must be warning signs, easily to be identified. Due to the state of stress and strain
in the structural elements, the elements failure can occur without any warning
signs. The undesired failure types are considered failure through compression,
shear force, torsion, and buckling.

In some particular situations, due to the topology of the structural assembly,

the structural elements are subjected to a state of stress and strain due to other
types of loads than has been expected at the structural design phase. In the case of
beam-type elements, compression forces can lead to unexpected failure, due to the
joints which are considered potential tension concentrators. The failure through
shear forces and torsion are in close relation with the state of stress and strain in
the structural element, respectively to avoid them one can take actions on a local
scale regarding the design of the structural element. The measures include adopting
appropriate cross-sections and quantities of materials. The loss of stability due to
extensive axial loads is high in the case of slender beam-type elements. To avoid
the present phenomena the critical buckling load must be increased above the value
of the elastic resistance of the cross-section. This is possible by modelling beam-
type elements with variable cross-sections and by designing statically determinate
structural assemblies for which the state of stress and strain can be computed
uniquely and using high-quality materials.

74
References:
1. Karabalis, D.L.; Beskos, D.E. Static, dynamic and stability analysis of
structures composed of tapered beams. Comput. Struct. 1983, 16, 731–748.
[Google Scholar] [CrossRef]
2. Das, D.; Sahoo, P.; Saha, K. Dynamic analysis of non-uniform taper bars in
post-elastic regime under body force loading. Appl. Math. Model. 2009, 33,
4163–4183. Available
online: https://www.sciencedirect.com/science/article/pii/S0307904X090
00584 (accessed on 4 January 2023). [CrossRef]
3. Yilmaz, Y.; Girgin, Z.; Evran, S. Bucling Analyses of Axially Functionally
Graded Nonuniform Columns with Elastic Restraint Using a Localized
Differential Quadrature Method. Math. Probl. Eng. 2013, 2013, 1–12.
[Google Scholar] [CrossRef] [Green Version]
4. Abdel-Lateef, T.H.; Dabaon, M.A.; Abdel-Moez, O.M.; Salama, M.I.
Buckling Loads of Columns with Gradually changing Cross-Section
subjected to Combined Axial Loading. In Proceedings of the Fourth
Alexandria International Conference on Structural and Geotechnical
Engineering, Alexandria, Egypt, 2–4 April 2001. [Google Scholar]
5. Guarracino, F., 2007. Considerations on the numerical analysis of initial
post-buckling behavior in plates and beams. Thin-Walled Structures, 45(10-
11), pp.845-848. ]Google Scholar[
6. Rhodes, J., 2002. Buckling of thin plates and members—and early work on
rectangular tubes. Thin-Walled Structures, 40(2), pp.87-108. ]Google
Scholar[

7. Saraçaoğlu, M.H.; Uzun, S. Stability analysis of columns with variable cross-

sections. J. Struct. Eng. Appl. Mech. 2020, 3, 169–179. [Google Scholar]
[CrossRef]
8. Ruocco, E.; Wang, C.M.; Zhang, H.; Challamel, N. An approximate model
for optimizing Bernoulli columns against buckling. Eng. Struct. 2017, 141,
316–327. [Google Scholar] [CrossRef] [Green Version]
9. Ruocco, E.; Zhang, H.; Wang, C.M. Hencky bar-chain model for buckling
analysis of non-uniform columns. Structures 2016, 6, 73–84. Available
online: https://www.sciencedirect.com/science/article/abs/pii/S23520124
16000187 (accessed on 4 January 2023). [CrossRef]

75
10. Ruiz-Teran, A.M. and Gardner, L., 2008. Elastic buckling of elliptical
tubes. Thin-Walled Structures, 46(11), pp.1304-1318.
a. Find online at:
(https://www.sciencedirect.com/science/article/pii/S02638231080005
42)

11. Wang, C.M., Xiang, Y., Kitipornchai, S. and Liew, K.M., 1994. Buckling
solutions for Mindlin plates of various shapes. Engineering
Structures, 16(2), pp.119-127.
12. Find online at:
(https://www.sciencedirect.com/science/article/pii/014102969490037X

13. N. Kharghani et al.

https://www.sciencedirect.com/science/article/pii/S0997753817309488
14. M. Darvizeh et al.
https://www.sciencedirect.com/science/article/pii/S0263822303001338

15. https://www.sciencedirect.com/science/article/pii/S0065215608700309

16. https://www.sciencedirect.com/science/article/pii/S0263822311002832

76