1 s2.0 S0020768321000688 Main

International Journal of Solids and Structures 222–223 (2021) 111003
Contents lists available at ScienceDirect
International Journal of Solids and Structures

journal homepage: www.elsevier.com/locate/ijsolstr
A unified approach to the Timoshenko 3D beam-column element

tangent stiffness matrix considering higher-order terms in the strain
tensor and large rotations
Marcos Antonio Campos Rodrigues a,⇑, Rodrigo Bird Burgos b, Luiz Fernando Martha c
a
Espírito Santo Federal University, Department of Civil Engineering, Avenida Fernando Ferrari, 514, Goiabeiras, Vitória, ES 29075-910, Brazil
b
State University of Rio de Janeiro, Department of Structures and Foundations, Rua São Francisco Xavier, 524, Maracanã, Rio de Janeiro, RJ 20550-900, Brazil
c
Pontifical Catholic University of Rio de Janeiro, Department of Civil Engineering, Rua Marques de São Vicente, 225, Gávea, Rio de Janeiro, RJ 22451-900, Brazil
a r t i c l e i n f o a b s t r a c t
Article history: A structural geometric nonlinear analysis using the finite element method (FEM) depends on the consid-
Received 10 July 2020 eration of five aspects: the interpolation (shape) functions, the bending theory, the kinematic description,
Received in revised form 14 January 2021 the strain–displacement relations, and the nonlinear solution scheme. As the FEM provides a numerical
Accepted 16 February 2021
solution, the structure discretization has a great influence on the analysis response. However, when
Available online 04 March 2021
applying interpolation functions calculated from the homogenous solution of the differential equation
of the problem, a numerical solution closer to the analytical response of the structure is obtained, and
Keywords:
the level of discretization could be reduced, as in the case of linear analysis. Thus, to reduce this influence
Tangent stiffness matrix
Analytical interpolation functions
and allow a minimal discretization of the structure for a geometric nonlinearity problem, this work uses
Timoshenko beam theory interpolation functions obtained directly from the solution of the equilibrium differential equation of a
Nonlinear geometric analysis deformed infinitesimal element, which includes the influence of axial forces. These shape functions are
used to develop a complete tangent stiffness matrix in an updated Lagrangian formulation, which also
integrates the Timoshenko beam theory, to consider shear deformation and higher-order terms in the
strain tensor. This formulation was implemented, and its results for minimal discretization were com-
pared with those from conventional formulations, analytical solutions, and Mastan2 v3.5 software. The
results clearly show the efficiency of the developed formulation to predict the critical load of plane
and spatial structures using a minimum discretization.
Ó 2021 Elsevier Ltd. All rights reserved.
1. Introduction so-called cubic Hermitian interpolation functions (Rodrigues

et al., 2019). In this case, the formulation does not consider any
The continuous (analytical) behavior of a solid can be approxi- other approximation except those already covered in the analytical
mated by a discrete solution. Usually, the discrete response is idealization of the element behavior. This explains the fact that, in
obtained by nodal displacements, and an approximated continuous linear elastic analysis, the structure response of this type of model
solution can be found by means of interpolating (shape) functions. does not depend on the level of discretization.
However, the discrete solution using the FEM introduces simplifi- However, for geometric nonlinear or second-order analysis, in
cations in the mathematical idealization of the structure behavior which equilibrium should be considered in the deformed configu-
as the interpolation functions that define the deformed configura- ration, Hermitian interpolation functions do not represent the
tion of a structure are not compatible with the mathematical ide- analytical response of the structure. To cope with this problem,
alization of the response of a continuous medium (Martha, 2018). high-order finite elements can be used (So and Chan, 1991;
In a linear elastic analysis of frame models with beam elements Zheng and Dong, 2011; Rodrigues et al., 2016). Burgos et al.,
with constant cross-sections, interpolating functions are obtained (2005) employed classic linearization from the stability problem
from the homogeneous solution of the equilibrium differential and used additional degrees of freedom within the elements to
equation of an undeformed infinitesimal element, leading to the calculate the critical load by eigenvalue analysis.
Another way of improving a second-order analysis is to use sta-
⇑ Corresponding author at: Avenida Fernando Ferrari, 514, Nexem, Goiabeiras, bility functions (Chen and Lui, 1991; Aristizábal-Ochoa, 1997,
Vitória, ES 29075-910, Brazil. 2007, 2008, 2012). Some authors have used the consistent field
E-mail address: rodriguesma.civil@gmail.com (M.A.C. Rodrigues).
https://doi.org/10.1016/j.ijsolstr.2021.02.014
0020-7683/Ó 2021 Elsevier Ltd. All rights reserved.
Marcos Antonio Campos Rodrigues, Rodrigo Bird Burgos and Luiz Fernando Martha International Journal of Solids and Structures 222–223 (2021) 111003
Nomenclature
General geometric parameters Nh2 ; Nh5 cross-section rotation interpolation functions of a bar
x; y; z bar local axis in longitudinal (x) and transversal (y; z) for a transverse displacement.
directions. Nh3 ; Nh6 cross-section rotation interpolation functions of a bar
dx; dy; ds infinitesimal increment in the (x,y) directions or in arc for a rotation displacement.
length.
l; L bar length. Strain and stress tensor parameters
V spatial vector position of a point in space. exx normal deformation.
R spatial transformation matrix. cxy ; cxz shear distortion in the plane xy and xz.
fag; fbg point in a cross-section. c shear distortion.
h generic cross-section height. eijðtþDtÞ Green-Lagrange strain tensor in an unknown configura-
A cross-section area. tion.
v form factor that defines the effective shear area. eðtÞ
ij
Green-Lagrange strain tensor in a known configuration.
I cross-section moment of inertia in relation to an axis. Deij incremental strain.
Jp cross-section polar moment of inertia. Deij linear component of incremental strain.
k element slenderness given by:L=h Dexx axial deformation from linear component of incremen-
X generic auxiliary parameter for characterization of bar tal strain.
elements. Dexy shear distortion from linear component of incremental
strain.
Displacement field parameters Dgij nonlinear component of incremental strain.
u axial displacement of a point inside a bar in the x direc- Dgxx axial deformation from nonlinear component of incre-
tion. mental strain.
v; w transverse displacement of a point inside a bar in the y Dgxy shear distortion from nonlinear component of incre-
or z direction. mental strain.
u0 ; v 0 ; w0 axial and transversal displacements at the cross- sxx normal stress.
section center of gravity. sxy ; sxz shear stress.
v h ; v p homogeneous and particular parts of analytical solution ðtþDtÞ
Sij second Piola-Kirchoff stress tensor.
for displacements and rotations. sðtÞ
ij
Cauchy stress tensor.
vg; vl global and local solutions for displacements and rota- Dsij stress increment.
tions. E Young’s modulus.
h cross-section rotation in relation to an axis. G modulus of rigidity (shear modulus).
0 0
d1 ; d4 axial displacements at the initial and final nodes of a m Poisson’s ratio.
bar. C ijkl material constitutive tensor.
0 0
d2 ; d5 transverse displacements at the initial and final nodes of
a bar. External loads and internal forces parameters
0 0
d3 ; d6 cross-section rotation at the initial and final nodes of a q load rate of transverse force on the bar.
0 bar. P constant generic axial load acting on infinitesimal ele-
d displacement vector of a bar. ment.
fug axial displacement vector of a bar. M generic bending moment acting on a bar element.
fv g; fwg transverse displacements vectors of a bar. l; K auxiliary parameters for the differential equation devel-
ci ; fC g constant of integration. opment of the deformed infinitesimal element equilib-
½X interpolation polynomials matrix of the displacement rium.
field. N normal force.
½H interpolation polynomials matrix of the displacement V vertical component of acting force in a cross-section.
field at nodes. Q y; Q z shearing force in directions y and z.
½N interpolation functions matrix. Q shearing force.
fNu g axial displacement interpolation functions vector. Mx twisting moment.
fNv g transverse displacement interpolation functions vector My ; Mz bending moment.
on the local xy plane of a bar. M y1 ; M y2 bending moment acting at the initial and final nodes of
fNw g transverse displacement interpolation functions vector a bar in y direction.
on the local xz plane of a bar. M z1 ; Mz2 bending moment acting at the initial and final nodes of
fNhx g cross-section rotation interpolation functions vector a bar in z direction.
around the local axis x.
Nhy cross-section rotation interpolation functions vector Stiffness matrix parameters
around the local axis y. RðtþDtÞ virtual work due to external loading.
fNhz g cross-section rotation interpolation functions vector U; U NL linear and nonlinear component of virtual work expres-
around the local axis z.
sion.
Nu1 ; Nu4 axial displacement interpolation functions at the initial K g;Rotfin finite rotations influence on the geometric stiffness ma-
and final nodes of a bar. trix of a generic plane element.
Nv2 ; N v5 transverse displacement interpolation functions of a bar
for a transverse displacement.
Nv3 ; N v6 transverse displacement interpolation functions of a bar
for a rotation displacement.
2
approach and transformed the shape interpolation functions (axial placement relations. The matrix is constructed for a spatial ele-
and transverse displacements and rotation) into a power series ment by an updated Lagrangian formulation considering higher-
(Yunhua, 1998; Tang et al., 2015). The consistent field approach order terms in the strain tensor using interpolation functions
eliminates ‘‘shear locking” and ‘‘membrane locking” effects that obtained directly from the equilibrium differential equation of a
usually appear in a geometric nonlinear analysis that employs just deformed infinitesimal element, including shear deformation
one element per member. according to the Timoshenko beam theory. In addition, the element
Other works have a formulation based on the equilibrium equa- tangent stiffness matrix is adjusted to consider finite rotations. The
tion of the infinitesimal element. Some include the axial effect con- fifth important aspect of a nonlinear analysis is the incremental
sidering the equilibrium in the deformed configuration (Davis solution scheme. This paper does not address this problem, and a
et al., 1972; Nukulchai et al., 1981; Goto and Chen, 1987; Chan standard Newton-Raphson solution procedure was used.
and Gu, 2000; Balling and Lyon, 2011). Others incorporate addi- Since this formulation involves trigonometric and hyperbolic
tional effects, such as elastic foundations (Areiza-Hurtado et al., basic functions, the resulting stiffness matrix can be conveniently
2005; Aydogan, 1995; Burgos and Martha, 2013; Chiwanga and rewritten using a Taylor series expansion with more terms than
Valsangkar, 1988; Eisenberger and Yankelevsky, 1985; Morfidis, the usual formulations. In this paper, tangent matrices with up to
2007; Morfidis and Avramidis, 2002; Onu, 2000, 2008; Shirima 3 and 4 terms were developed. This formulation was implemented
and Giger, 1992; Ting and Mockry, 1984; Zhaohua and Cook, in Framoop (Martha and Parente Junior, 2002), the solver used in
1983) and interlayer slip for sandwich beams (Ha, 1993; Ftool (Martha, 1999) structural analysis program.
Girhammar and Gopu, 1993). However, these studies directly for-
mulate the stiffness coefficients of the elements without explicitly 2. Differential equilibrium relationships in beam-columns
presenting expressions for interpolation functions based on the
differential equation of the problem. This section formulates the differential equations that define
In addition to the enrichment of interpolating the field variables the analytical behavior of a deformed infinitesimal beam element
in an element, another important aspect to improve the response considering both the Euler-Bernoulli and Timoshenko beam theo-
of a structure when performing a geometric nonlinear analysis is ries. Fig. 2.1presents a deformed infinitesimal element subjected
related to the bending theory. The most commonly used bending to a distributed transverse load q and a constant axial load P.
solution for frame elements is the Euler-Bernoulli beam theory The equilibrium of the deformed infinitesimal beam element
(EBBT). However, the effects of shear deformation are important leads to Eqs. (2.1) and (2.2)
when predicting the behavior of beam-columns with a moderate
X dV ðxÞ
slenderness ratio or with a small shear-to-bending ratio, and the F y ! dV þ qðxÞdx ¼ 0 ! ¼ qðxÞ ð2:1Þ
Timoshenko beam theory (TBT) provides good results dx
(Timoshenko and Gere, 1963; Friedman and Kosmatka, 1993; X dx2
Pilkey et al., 1995; Schramm et al., 1994). For this reason, many M o ! dM ðV þ dV Þdx P:dv þ qðxÞ ¼0 ð2:2Þ
works cited above also consider this effect, such as references 2
(Burgos and Martha, 2013; Onu, 2008). where v ðxÞ is the infinitesimal element transverse displace-
Moreover, a kinematic description of motion is needed when ment, qðxÞ is the transverse distributed load, V ðxÞ is the vertical
formulating an incremental geometrically nonlinear analysis. component of the force acting on the cross-section, P is the hori-
Based on the choice of the reference configuration, the following zontal component, and MðxÞ is the bending moment.
three kinematic descriptions are commonly used in structural Considering that the element has a constant cross-section, using
mechanics to formulate the nonlinear system of equations: Total Eq. (2.2) and the approximate relation between bending moment
Lagrangian (TL), Updated Lagrangian (UL), and Corotational (CR). and curvature, M ðxÞ ¼ EI dh=dx, where hðxÞ corresponds to the
When consistently developed, the total Lagrangian and the cross-section rotation, the differential equation of the problem
updated formulation produce the same results (Mcguire et al., can be written according to expression (2.3).
2000). These formulations are developed in (Mcguire et al., 2000;
d v ð xÞ
2
Bathe and Bolourchi, 1979; Conci, 1988; Yang and Leu, 1994;
d hðxÞ
EI V ðxÞ P
Yang and Kuo, 1994; Chen, 1994; Bathe, 1996). However, in gen- dx2 dx
d v ð xÞ
3 2
eral, the Timoshenko theory is not considered, or in some cases, d hðxÞ dV ðxÞ
¼ 0 ! EI P ¼0 ð2:3Þ
the formulation yields additional degrees of freedom, which dx3 dx dx2
requires static condensation to reduce the matrix order (Bathe
and Bolourchi, 1979; Aguiar et al., 2014). Additionally, most stud-
ies that consider shear deformation employ a corotational formu-
lation, as proposed by (Battini, 2002; Crisfield, 1991; Pacoste and
Eriksson, 1995, 1997; Santana and Silveira, 2019; Silva et al., 2016).
Furthermore, the strain–displacement relation plays an impor-
tant role when developing a geometric stiffness matrix. The con-
sideration of higher-order terms in the strain tensor improves
the accuracy of the geometric nonlinear analysis (Chen, 1994;
Conci, 1988; Yang and Kuo, 1994; Yang and Leu, 1994). More
recently, the authors (Rodrigues et al., 2019) considered higher-
order terms in the strain tensor in the formulation of a Timoshenko
element using Hermitian interpolation functions and an updated
Lagrangian formulation.
The main objective of this paper, which differentiates this work
from others found in the literature, is to formulate the tangent
stiffness matrix of a frame element integrating four important
aspects that improve geometric nonlinear analysis: interpolation
functions, beam theory, kinematic description, and strain–dis- Fig. 2.1. Equilibrium of a deformed beam element.
3
The static fundamental relation in Eq. (2.1) shows that the gra- Equation (2.9) can be written using the relation
dient of the vertical component of the force acting on the cross- Q ðxÞ ¼ dMðxÞ=dx in Eq. (2.8). According to the relation between
section is equal to the acting distributed load. Therefore, substitut- bending moment and curvature M ðxÞ ¼ EIdh=dx, Eq. (2.10), which
ing this relation in Eq. (2.3), the differential equilibrium relation of relates transverse displacement v ðx) with cross-section rotation
a deformed infinitesimal element can be found, as shown in hðxÞ, is obtained.
expression (2.4).
dMðxÞ dv ðxÞ
¼ vGA: hðxÞ ð2:9Þ
d v ð xÞ
3 2
d h dx dx
EI P ¼ qðxÞ ð2:4Þ
dx3 dx2
dv ðxÞ
2
EI d h
¼ hðxÞ ð2:10Þ
2.1. Euler-Bernoulli beam theory dx vGA dx2
Finally, applying equation (2.10) in expression (2.4), the equilib-
The EBBT considers the rotation as the gradient of the trans- rium differential relation of a deformed infinitesimal element, con-
verse displacement (h ¼ dv =dx). Thus, the differential relation sidering Timoshenko beam theory, can be written according to
developed in expression (2.4) can be written according to (2.5). expressions (2.11) or (2.12).
d
4
v ð xÞ P d v ðxÞ qðxÞ
2
P
3
d h dhðxÞ
¼ ð2:5Þ EI 1 þ P ¼ qðxÞ ð2:11Þ
dx4 EI dx2 EI vGA dx3 dx
2.2. Timoshenko beam theory 3

d h P dhðxÞ qðxÞ
¼ ð2:12Þ
dx3 1 þ vGA
P
EI dx 1 þ vGA
P
EI
In TBT, shear distortion ðcÞ is constant for each cross-section (no
warping) and is considered an additional rotation of the section.
Equation (2.11), or alternatively Eq. (2.12), relates the cross-
Therefore, the section rotation and the transverse displacement
section rotation hðxÞ with the applied distributed transverse
are not associated, and both are considered independent variables,
loadqðx), considering the element axial force and the shear and
as indicated in Fig. 2.2.
bending rigidity parameters. The next section shows the solution
Therefore, for the deformed element, internal forces are calcu-
of differential Eqs. (2.10) and (2.12) for an isolated beam element.
lated using the cross-sectional rotation ðdv =dx ¼ h þ cÞ shown in
Fig. 2.3 and expressed in Eqs. (2.6) and (2.7).
3. Solutions of the differential equations
NðxÞ ¼ Pcosðh þ cÞ V ðxÞsenðh þ cÞ ! N ðxÞ
dv ðxÞ The solution of the equilibrium differential relations presented
¼ P V ðxÞ ð2:6Þ before can be obtained by the composition of a homogeneous solu-
dx
tion, v h ðxÞ, with a particular solution, v p ðxÞ, according to
Q ðxÞ ¼ Psenðh þ cÞ þ V ðxÞcosðh þ cÞ ! Q ðxÞ v ðxÞ ¼ v h ðxÞ þ v p ðxÞ.
dv ðxÞ In the direct stiffness method, the solution is obtained by
¼P þ VðxÞ ð2:7Þ the superposition of a global and a local solution, hence
dx
v ðxÞ ¼ v g ðxÞ þ v l ðxÞ. The global solution is the one obtained by
Substituting Eqs. (2.7) into (2.2), the static fundamental relation using nodal displacements and rotations as coefficients for the
between the shear force and the bending moment is obtained: interpolating functions. The local solution is a fixed-end ele-
Q ðxÞ ¼ dMðxÞ=dx. The shear force acting on the section is given ment solution, thus presenting null values for nodal displace-
by Eq. (2.8) ments and notations. Martha (2018) used the scheme in

dv ðxÞ Fig. 3.1 to explain this superposition. Nodal displacements and
Q ðxÞ ¼ vGA:cðxÞ ! Q ðxÞ ¼ vGA: hðxÞ ð2:8Þ rotations of the final solution are obtained directly from the
dx
global solution. Within an element, displacements are obtained
in which G is the material shear modulus, A is the cross-section by superposition of the local and global solutions. In a classical
area, and v is the factor that defines the cross-section effective area FEM analysis for continuous domains, generally, the local solu-
for shear. tion is not available, and the final solution is approximated
by the global solution.
If the right-side of the differential equation is null, then there is
no need for a particular solution and the homogeneous part is the
final solution, then v ðxÞ ¼ v h ðxÞ. Similarly, if there is no transverse
force within the element (q(x) = 0), the local solution is null, giving
v ðxÞ ¼ v g ðxÞ. Thus, the homogeneous and the global solutions are
equivalent: v h ðxÞ ¼ v g ðxÞ. Therefore, when adopting analytical
solutions for the behavior of the element, the displacements and
nodal rotations of the discrete model are exact, independently of
discretization. Moreover, in case a fixed-end local element solution
is available, nodal displacements and rotations of the final solution
are obtained directly from the global solution, also independently
of discretization.
These observations justify the use of interpolation functions cal-
culated from the homogenous solution. When the differential
equation is calculated from the equilibrium of an undeformed
Fig. 2.2. Shear deformation in the Timoshenko beam theory. infinitesimal element, the so-called Hermitian cubic functions cor-
respond to the homogeneous solution for EBBT, and because of
4
Fig. 2.3. Internal forces, displacements, and rotations in the Timoshenko beam theory.
Fig. 3.1. Solution composition from the direct stiffness method.
that, linear analyses can be performed without discretization. Sim- d

4
v ð xÞ 2
d v ðxÞ pffiffiffiffiffiffiffiffiffiffi
ilar cubic functions are obtained from the undeformed infinitesi- l2 ¼ 0; l ¼ P=EI ð3:1Þ
dx4 dx2
mal equilibrium for TBT. The same applies to the functions that
represent the homogenous solution of the equations obtained from v h ðxÞ ¼ c1 elx þ c2 elx þ c3 x þ c4 ð3:2Þ
the equilibrium of a deformed infinitesimal element.
where c1 , c2 ,c3 and c4 are the coefficients of an exponential
Therefore, this section seeks the homogeneous solution of the
function.
differential equations of equilibrium of a frame element in the
In EBBT, h ¼ dv =dx, and the cross-section rotation is obtained by
deformed configuration. The next sections indicates how the inter-
equation (3.3).
polation (shape) functions resulting from the solutions to these dif-
ferential equations are obtained; and Section 5 shows how these hh ðxÞ ¼ c1 lelx c2 lelx þ c3 ð3:3Þ
interpolation functions are used in the formulation of the tangent
The particular solution does not introduce any additional coef-
stiffness matrix of a frame element in its local axis system. It
ficients that need to be determined. Thus, the coefficients of the
should be noted that the homogenous solution considers approxi-
exponential function can be calculated from the boundary
mations inherent to the analytical model adopted for the behavior
conditions.
of beam-column elements, such as the hypothesis of small rota-
For a tensile force, l is a real number, and the homogeneous
tions in the equilibrium of the infinitesimal element in the
solution of the differential equation can be written by hyperbolic
deformed configuration.
functions according to (3.4) and (3.5). In the case of a compressive
force, l is a complex number, and the displacement can be written
3.1. Euler-Bernoulli beam theory
by trigonometric functions (3.6) and (3.7).
pffiffiffiffiffiffiffiffiffiffi
Considering the EBBT and qðxÞ ¼ 0, the homogeneous solution v h ðxÞ ¼ c1 sinhðlxÞ þ c2 coshðlxÞ þ c3 x þ c4 ; l ¼ P=EI ð3:4Þ
of the differential equation presented in (2.5) can be obtained from
relation (3.1) and written as given in (3.2), hh ðxÞ ¼ c1 lcoshðlxÞ þ c2 lsinhðlxÞ þ c3 ð3:5Þ
5
pffiffiffiffiffiffiffiffiffiffiffiffiffi
v h ðxÞ ¼ c1 sinðlxÞ þ c2 cosðlxÞ þ c3 x þ c4 ; l ¼ P=EI ð3:6Þ
hh ðxÞ ¼ c1 lcosðlxÞ c2 lsinðlxÞ þ c3 ð3:7Þ
3.2. Timoshenko beam theory
In TBT, the cross-section rotation is not the gradient of the

transverse displacement, and it is not possible to write just one dif-
ferential relation; thus, Eqs. (2.10) and (2.12) need to be solved.
The homogenous solution can be calculated considering qðxÞ ¼ 0;
thus, the differential equations are written as shown in (3.8) and Fig. 4.1. Deformed configuration of an isolated element.
(3.9).
4. Interpolation functions
dv ðxÞ
2
EI d h
¼ hðxÞ ð3:8Þ
dx vGA dx2 The deformed configuration of an isolated element can be
described by interpolating nodal displacements according to
3
d h dhðxÞ l pffiffiffiffiffiffiffiffiffiffi Fig. 4.1 and Eqs. (4.1) to (4.4).
K2 ¼ 0; K ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; l ¼ P=EI ð3:9Þ
dx3 dx 1þ EI
:l 2 0 0
vGA u0 ðxÞ ¼ Nu1 ðxÞd1 þ Nu4 ðxÞd4 ! u0 ðxÞ ¼ fNu ðxÞgfug ð4:1Þ
Using constantX, expression (3.10), introduced by Reddy
v 0 ðxÞ ¼ Nv2 ðxÞd2 þ Nv3 ðxÞd3 þ Nv5 ðxÞd5 þ Nv6 ðxÞd6 ! v 0 ðxÞ ¼ fNv ðxÞgfv g
0 0 0 0
(Reddy, 1997) and used in (Burgos and Martha, 2013; Martha
and Burgos, 2015; Martha and Burgos, 2014), the differential equa- ð4:2Þ
tions can be rewritten as given in (3.11) and (3.12).
hðxÞ ¼ N h2 ðxÞd2 þ N h3 ðxÞd3 þ N h5 ðxÞd5 þ N h6 ðxÞd6 ! hz ðxÞ ¼ fN hz ðxÞgfv g
0 0 0 0
EI 1 ð4:3Þ
X¼ ð3:10Þ
vGA L2
u0 ðxÞ ¼ fNu ðxÞgfug v 0 ðxÞ ¼ fNv ðxÞgfv g hz ðxÞ ¼ fNhz ðxÞgfv g
dv ðxÞ ð4:4Þ
2
d h
¼ hðxÞ XL2 2 ð3:11Þ
dx dx When interpolating functions are obtained from the homoge-
3 neous solution of the equilibrium differential equation of an
d h dhðxÞ l
K2 ¼ 0; K ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð3:12Þ infinitesimal element, the analytical behavior of the element is
dx3 dx
1 þ Xl2 L2 represented. In this work, the equilibrium of a deformed infinites-
imal element is considered taking the axial force into account.
Solving equation (3.12), the homogeneous solution for the
cross-section rotation is given by (3.13). This solution can be 4.1. Euler-Bernoulli beam theory
applied in Eq. (3.11). Thus, the transverse displacement is written
according to (3.14). The interpolation functions can be calculated based on the
Kx Kx
homogeneous solution for the EBBT presented in Eqs. (3.2) and
hh ðxÞ ¼ K c1 e c2 e þ c3 ð3:13Þ
(3.3). First, the displacement is written in matrix form according
Kx to Eq. (4.5).
v h ð xÞ ¼ 1 XL2 K2 c1 e c2 eKx þ c3 x þ c4 ð3:14Þ 8 9
v h ðxÞ ¼ c1 elx þ c2 elx þ c3 x þ c 4 > 1 >
>
c
>
v 0 ðxÞ lx
e elx x 1 < c2 =
Finally, for a tensile force, l is a real number, and the homoge- !
hðxÞ
¼
lelx lelx 1 0 > >
c3 >
¼ ½X fC g
>
: ;
neous solution of the differential equation can be written by the hh ðxÞ ¼ c1 lelx c2 lelx þ c3 c4
hyperbolic functions in (3.15) and (3.16); for a compressive force, ð4:5Þ
l is a complex number, and the displacement is written by the
The boundary conditions are obtained by evaluating the homo-
trigonometric functions in (3.17) and (3.18).
geneous solution of these displacements (v h ðxÞ and hh ðxÞ) at the
extreme nodes of the bar, as shown in expression (4.6).
hh ðxÞ ¼ K½c1 coshðKxÞ þ c2 sinhðKxÞ þ c3 ; 8 0 9 2 3
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8 9
> d20 >
> > > v 0 ð0Þ > 1 1 0 1
K ¼ l= 1 þ Xl2 L2 ð3:15Þ n 0o >< > = < = n 0o 6
l l 7
d3 hð0Þ 1 0
d ¼ ¼ ! d ¼6 7
>
> d
0
>
5 >
>
: v 0 ðLÞ >
;
4 eLl e Ll
L 1 5
>
: 0 ; >
hðLÞ l eLl
l e Ll
1 0
d
v h ð xÞ ¼
1 XL2 K2 ½c1 sinhðKxÞ þ c2 coshðKxÞ þ c3 x þ c4 ;
6
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8 9
K ¼ l= 1 þ Xl2 L2 ð3:16Þ >
> c1 >
>
>
<c >
=
2
: ¼ ½H f C g ð4:6Þ
> c3
> >
>
hh ðxÞ ¼ K½c1 cosðKxÞ c2 sinðKxÞ þ c3 ; >
: >
;
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c4
K ¼ l= 1 Xl2 L2 ð3:17Þ
Finally, interpolation functions are calculated using Eqs. (4.4)–(4.6),
resulting in relation (4.7). The results are presented in the work by
Rodrigues (2019) for exponential, hyperbolic and trigonometric
v h ð xÞ ¼
1 þ XL2 K2 ½c1 sinðKxÞ þ c2 cosðKxÞ þ c3 x þ c4 ;
functions and were implemented in file ShpFuncBeamEulerBernoulli
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
using MATLAB and C in an open source code that can be accessed
K ¼ l= 1 Xl2 L2 ð3:18Þ
and used (Rodrigues et al., 2020, 2021).
6

v 0 ð xÞ 5. Local stiffness matrix
¼ ½X ½H1 fd0 g ) ½N ¼ ½X ½H1 ð4:7Þ
hðxÞ
According to Fig. 5.1, the stiffness matrix development can be
performed by initially analyzing a two-dimensional structure (xy
4.2. Timoshenko beam theory and xz planes). Then, for spatial structures, as shown in Fig. 5.2,
the integration between the planes is coupled to the stiffness
The homogeneous solution for the TBT is presented in Eqs. matrix. Finally, to consider large rotations, the concept of finite
(3.13) and (3.14). Thus, the same development is performed to cal- rotations was added to the matrix.
culate the interpolation functions, according to expressions (4.8)
and (4.9). 5.1. Updated Lagrangian formulation

v h ð xÞ ¼ 1 XL2 K2 c1 eKx c2 eKx þ c3 x þ c 4 The stiffness matrices are calculated considering the updated
v 0 ð xÞ Lagrangian formulation, and the steps shown in this work have
! ¼ ½X fC g
hðxÞ
hh ðxÞ ¼ K c1 eKx þ c2 eKx þ c3 been presented in the work by Mcguire et al. (2000), Chen
(1994), Bathe (1996), Aguiar et al. (2014) and Rodrigues et al.
8 9 (2019).
2 3 > c1 >
>
< c >= The virtual work of the internal forces must be equal to the vir-
1 XL2 K2 eKx 1 XL2 K2 eKx x 1
½X ¼ 4 5 2
fC g ¼ ð4:8Þ
KeKx KeKx > c3 >
> > tual work of the external forces, according to (5.1),
1 0 : ;
c4 Z
ðtþDt Þ ðtþDt Þ
Sij deij dV ¼ RðtþDtÞ ð5:1Þ
8 0 9 8 9 V
> d2 > 8 9 >
> 1 XL2 K2 ðc1 c2 Þ þ c4 >
>
> v 0 ð0Þ
>
> >
> > >
>
> > > >
> >
> >
> ðtþDt Þ
< d0 >
= > < hð0Þ >
= >
< Kðc1 þ c2 Þ þ c3
>
= where Sij corresponds to the second Piola-Kirchoff stress
fd0 g ¼ 3
¼ ¼ ¼ ½HfC g
> v 0 ðLÞ
> ðtþDt Þ
> d5 >
> > > > >
0
>
> >
> > :
>
>
;
>
>
>
1 XL2 K2 c1 eKL c2 eKL þ c3 L þ c4 >
>
>
tensor, ij e corresponds to the Green-Lagrange strain tensor,
: 0 >
> ; hðLÞ >
> >
>
d6 : KL
K c1 e þ c2 e KL
þ c3 ; ðtþDt Þ
and R corresponds to the virtual work due to external loading,
which could include body, surface, inertia and damping forces
2 3 (Aguiar et al., 2014). It is important to note that the term d is
ð1 XL2 K2 Þ ðXL2 K2 1Þ 0 1 related to the application of a virtual displacement to consider
6 7
6K K 1 0 7 the virtual displacement principle, while D is related to a small dis-
½H ¼ 6
6 eLK ð1 XL2 K2 Þ
7
7 ð4:9Þ
4 eLK ðXL2 K2 1Þ L 1 5 placement increment, changing the configuration.
LK However, the equilibrium equations of an unknown configura-
Ke KeLK 1 0
tion t þ Dt must be written using a known reference configuration
Once again, equation (4.7) can be used to calculate the interpo- t. Thus, the linearized incremental equation requires small dis-
lation functions. The results are also presented in the work by placement increments according to the equations in (5.2),
Rodrigues (2019) and implemented in file ShpFuncBeamTi- ðtþDt Þ
moshenko; an open source code that can be accessed and used
Sij ¼ stij þ Dsij eðijtþDtÞ ¼ etij þ Deij ð5:2Þ
((Rodrigues et al., 2020, 2021)). where stij corresponds to the Cauchy stress tensor, Dsij is the
stress increment and Deij is the deformation increment. In the
4.3. Axial interpolation functions known reference configuration, there is no deformation so
detij = 0, and the element is only subject to rigid body motion. Thus,
For the axial displacement and the cross-section rotation about the left side of Eq. (5.1) can be rewritten as in (5.3).
the x axis, in this research, a linear interpolation was adopted Z Z
ðtþDt Þ ðtþDt Þ
according to Mcguire et al. (2000), which is given in expression Sij deij dV ¼ stij þ Dsij d Deij dV
V
(4.10). ZV Z
¼ Dsij dDeij dV þ stij dDeij dV ð5:3Þ
0 0 x x
u0 ðxÞ ¼ Nu1 ðxÞd1 þ Nu4 ðxÞd4 Nu1 ðxÞ ¼1 Nu4 ðxÞ ¼ V V
L L The Green-Lagrange strain tensor has a linear (Deij ) and a non-
ð4:10Þ
linear part (Dgij ), according to expressions (5.4) to (5.6) in the xy
Some authors, such as Tang et al. (2015) and Silva et al., (2016); plane.
employed consistent interpolation functions for the axial displace-
Deij ¼ Deij þ Dgij ð5:4Þ
ment and not just linear interpolation.
Fig. 5.1. Bending and torsion for a spatial element.
7
5.3. Local stiffness matrix equations
The development of the local stiffness matrix using cubic func-

tions can be seen in Mcguire et al. (2000) considering Euler-
Bernoulli beam theory and in Rodrigues et al. (2019) for a
Timoshenko beam element. In this research, the adopted form to
achieve the complete matrix is the same as that presented in the
cited literature; however, the interpolation function is complete,
i.e., the axial force is considered.
From the displacement field in Eq. (5.14), the linear and nonlin-
Fig. 5.2. 3-D element (Mcguire et al., 2000). ear parts of the Green-Lagrange strain tensor components in Eq.
(5.6), can be rewritten as expressions (5.15) and (5.16), respec-
tively, for the Euler-Bernoulli beam theory.
@u @u @ v
Dexx ¼ Dexy ¼ þ ð5:5Þ @u @u0 @2v @ v @u
@x @y @x exx ¼ ¼ y 2 cxy ¼ þ ¼0 ð5:15Þ
@x @x @x @x @y
" #
@v @u @u @ v @ v
2 2
1 @u
Dgxx ¼ þ Dgxy ¼ þ ð5:6Þ 2 2
@u2
þ @@xv @u 2
þ @@xv þ y2 @@xv2 y @@xv2 @u
2
2 @x @x @x @y @x @y
2 2
gxx ¼ 12 @x
¼ 12 @x @x
ð5:16Þ
The stress increment is obtained from the material constitutive @u @u @ v @ v @ 2 v @ v @ v @u
relation. The consideration of the linear approximation for the gxy ¼ þ ¼y 2
@x @y @x @y @x @x @x @x
stress and strain increment leads to the expression in Eq. (5.7).
Meanwhile, for the Timoshenko beam theory, the linear and
Dsij ¼ C ijkl Dekl ¼ C ijkl Dekl Deij ¼ Deij ð5:7Þ nonlinear parts of the strain tensor components in Eq. (5.6) can
Based on equations (5.3) and (5.7), the virtual work equation be written according to Eqs. (5.17) and (5.18), using the displace-
becomes (5.8). ment field presented in (5.13).
Z Z Z @u @u0 @hZ @ v @u @ v 0
C ijkl Dekl dDeij dV þ stij dDeij dV þ stij Dgij dV ¼ RðtþDtÞ ð5:8Þ exx ¼ ¼ y cxy ¼ þ ¼ hZ ð5:17Þ
V V V
@x @x @x @x @y @x
The first integral leads to the elastic stiffness matrix, while the 2 2 2
gxx ¼ 12 @u 2
@x
þ @@xv ¼ 12 @u 2
@x
þ @@xv þ y2 @h
@x
z
y @u @hz
@x @x
third integral leads to the geometric stiffness matrix. The second
integral represents the virtual work of forces acting on the element @u @u @ v @ v @hz @u ð5:18Þ
gxy ¼ þ ¼y hz hz
in configuration t and is usually represented at the right-hand side @x @y @x @y @x @x
of the expression. However, in this case, the rotation components
are approximated by displacement derivatives, and the finite rota-
5.3.1. Elastic matrix
tion effect is not considered (Mcguire et al., 2000).
In the xy plane, the stress vector, the constitutive matrix, the
Plane xy (EBBT)
linear and nonlinear strain vectors are given according to equation
From the displacement field in (5.14) and the linear part of the
(5.9). Thus, relations (5.10)–(5.12) can be written.
( ) ( ) strain tensor in (5.15), the first integral of the virtual work equa-

s E 0 exx gxx tion, relation (5.10), can be written according to Eq. (5.19).
s ¼ xx C¼ e¼ g¼ Z Z
sxy 0 G cxy gxy
dU ¼ exx :Edexx dV þ cxy :Gdcxy dV
ð5:9Þ
V V
Z Z Z Z Z ! ! !
L
@u @2v @u @2v
C ijkl Dekl dDeij dV ¼ exx :Edexx dV þ cxy :Gdcxy dV ð5:10Þ ¼ y 2 E d yd 2 dx dA
V V V A 0 @x @x @x @x
Z Z Z
RL R R RL
stij dDeij dV ¼ sxx dexx dV þ sxy dcxy dV ð5:11Þ ! dU ¼ @u @u
d dx
0 @x @x
d @@xv2 dx E A y2 dAþ
E
2
A
dA þ @2 v
0 @x2
V V V
RL RL R R ð5:19Þ
@ 2 v @u
Z Z Z d @ v dx E A ydA
2
0 d dx E A ydA 0 @u
@x2 @x @x @x2
sij Dgij dV ¼ sxx dgxx dV þ sxy dgxy dV ð5:12Þ R R
V V V In the beam centroidal axis, A
y2 dA ¼ Iz ; A
ydA ¼ 0, and the
expression is reduced to (5.20).
5.2. Displacement field Z Z !
L
@u @u L
@2v @2v
dU ¼ d dx EA þ d dx EIz ð5:20Þ
The displacement field of a beam element is shown in Fig. 5.3 0 @x @x 0 @x2 @x2
and defined by Eq. (5.13).
From Eq. (4.4), the displacements can be written using the
uðx; yÞ ¼ u0 ðxÞ hðxÞ:y v ðx; yÞ ¼ v 0 ðxÞ ð5:13Þ interpolation functions resulting in (5.21). The notation
@=@x ¼ ð Þ’ is adopted for simplicity.
Considering the Euler-Bernoulli beam theory, the rotation is the
Z
derivative of the transverse displacement (h ¼ dv =dx); thus, the L T
displacement field can be rewritten according to (5.14).
dU ¼ fdugT EA N0u N0u dxfug þ fdv gT
0
Z
v 0 ð xÞ L T
uðx; yÞ ¼ u0 ðxÞ :y v ðx; yÞ ¼ v 0 ðxÞ ð5:14Þ EIz N00v N00v dxfv g ð5:21Þ
dx 0
8
Fig. 5.3. Beam displacement field.
The elastic matrix is formed when cubic functions are used. The axial part was already used for the xy plane. Finally, writing
However, in this work, the interpolation functions are those calcu- the expression using the interpolation functions, equation (5.26)
lated in section 4. can be achieved:
Z
Plane xy (TBT)
L T
dU ¼ fdwgT EIy N00w N00w dxfwg ð5:26Þ
0
The same idea can be used for the Timoshenko beam theory,
however, considering now the displacement field presented in where N w corresponds to the same N v interpolation functions.
(5.13). Thus, the first integral of the virtual work equation, relation Plane xz (TBT)
(5.10), is written by Eq. (5.22) or according to expression (5.23) For the Timoshenko beam theory considering the xz plane, the
when interpolation functions are used. displacement field becomes (5.27), and the equations for the elas-
Z Z tic matrix are given by relations (5.28) and (5.29).
dU ¼ exx :Edexx dV þ cxy :Gdcxy dV ! dU
@u @u0 @hy @w @u @w0
V V exx ¼ ¼ z cxz ¼ þ ¼ hy ð5:27Þ
Z Z L @x @x @x @x @z @x
@u @hz @u @hz
¼ y E y dx dA
A 0 @x @x @x @x RL @u @u
RL @hy @hy RL @w @w

Z Z L dU ¼ d dx EA þ d dx EIy þ d @x dx GA
@v @v
0 @x @x 0 @x @x 0 @x
þ hZ G d dhZ dx dA RL
A 0 @x @x þ 0
hy dhy dx GAþ
Z Z RL RL
L
@u @u @hZ @hZ L 0
hy d @w
@x
dx GA @w
0 @x
dhy dx GA
! dU ¼ d dx EA þ d dx EIz
0 @x @x 0 @x @x ð5:28Þ
Z L Z L
@v @v n on oT
þ d dx GA þ hZ dhZ dx GA RL R L T
0 @x @x 0 dU ¼ fdwgT EIy N 0hy N 0hy dxfwg þ fdwgT 0 GA N 0w N 0w dxfwg
0
Z L Z L R L T R L T
@v @v þ fdwgT 0 GA N hy N hy dxfwg fdwgT 0 GA Nhy N0w dxfwgþ
hZ d dx GA dhZ dx GA ð5:22Þ R
0 @x 0 @x
T
fdwgT 0 GA N0w N hy dxfwg
L
R L T R L T ð5:29Þ
dU ¼ fdugT 0 EA N0u N0u dxfug þ fdv gT 0 EIz N0hz N0hz dxfv g
R L T RL
þ fdv gT 0 GA N0v N0v dxfv g þ fdv gT 0 GAfNhz gfNhz gT dxfv gþ 5.3.2. Geometric matrix
RL T RL
fdv gT 0 GAfNhz g N0v dxfv g fdv gT 0 GA N0v fNhz gT dxfv g
Plane xy (EBBT)
ð5:23Þ
The third integral from the virtual work expression (5.8) gener-
ates the usual geometric stiffness matrix. Thus, the nonlinear part
Plane xz (EBBT) of the virtual work principle, considering the complete strain ten-
sor and higher order terms, is given by Eq. (5.30).
For the plane xz behavior, the linear part of the strain tensor can
be written just by switching the displacement v with w in equation Z Z
(5.24). Thus, the first integral of the virtual work equation is given dU NL ¼ sxx dgxx dV þ sxy dgxy dV
by expression (5.25). V V
@u @u0 @2w @w @u @w0 @w0 0 0 0 !2 1 1 1

Z Z
exx ¼ ¼ z 2 cxz ¼ þ ¼ ¼0 ð5:24Þ @v @2 v A y @ v @uAdxAdA
2 2 2
@
L
1 @u
@x @x @x @x @y @x @x dU NL ¼ t xx d@ @ þ þy 2
A 0 2 @x @x @x2 @x2 @x
Z Z ! Z Z L ! !
L
@u @u L
@2w @2w @ 2 v @ v @ v @u
dU ¼ d dx EA þ d dx EIy ð5:25Þ þ t xy d y: 2 : : dx dA
@x @x @x2 @x2 A 0 @x @x @x @x
0 0
9
Z !Z
2
1 @u
L
! dU NL ¼ d dx t xx dA
2 0 @x A
Z L !Z
@v
2
1
þ d dx txx dA
2 0 @x A
0 !2 1 Z Fig. 5.4. Frame element.
Z
1 @ L @2v
þ d dxA y2 t xx dA
2 0 @x2 A
Z L 2 !Z RL R 2RLR RL 2 R
@ v @u d @@xv dx A txx dA þ 12 0 d @h
@u 2
! dU NL ¼ 12 0
d @x
dx t dA
A xx 0
þ 12 @x
z
dx A y2 txx dA
d 2 dx txx ydA RL R R L @hz R RL R
0 @x @x A 0
d @h
@x
z @u
@x
dx t ydA þ 0 d @x hz dx A txy ydA 0 dhz @u
A xx @x
dx A txy dA
Z L 2 !Z ð5:34Þ
@ v @v
þ d 2 dx t xy ydA " ! !#
@x @x Z
@v
2 2 2
0 A
1 @u L
Iz @hz
Z L Z dU NL ¼ þ Pd þP d dx
@ v @u 2 0 @x @x A @x
d dx txy dA ð5:30Þ
@x @x Z L
0 A
@hz @u @u
R R R þ Mz d Q y d hz dx ð5:35Þ
Applying relations, A txx dA ¼ P; A txx ydA ¼ Mz ; A txy dA ¼ Q y ; 0 @x @x @x
the expression becomes (5.31).
Writing Eq. (5.35) with interpolation functions leads to Eq.
2 ! 0 !2 13 (5.36).
Z
@v I z @ @ 2 v A5
2 2
1 L4 @u R L T R L T
dU NL ¼ Pd þ þP d dx dU NL ¼ fdugT 0 P N0u N0u dxfug þ fdv gT 0 P N 0v N0v dxfv g
2 0 @x @x A @x2
RL T R L T
Z L" ! # þfdv gT 0 P IAz N0hz N0hz dxfv g þ fdv gT 0 M z N0hz N0u dxfug
@ 2 v @u @ v @u R L T RL T
þ Mz d Q yd dx ð5:31Þ þfdugT 0 Mz N0u N0hz dxfv g fdv gT 0 Q y fNhz g N0u dxfug
0 @x2 @x @x @x
RL
fdugT 0 Q y N0u fNhz gT dxfv g
Once again, the displacements can be written using the interpo- ð5:36Þ
lation functions in equation (4.4), resulting in expression (5.32).
R L T R L T
dU NL ¼ fdugT 0 P N0u N0u dxfug þ fdv gT 0 P N0v N0v dxfv gþ
RL T R L T Plane xz (EBBT)
þfdv gT 0 P IAz N00v N00v dxfv g þ fdv gT 0 Mz N00v N 0u dxfugþ
R L T R L T The nonlinear parts of the Green-Lagrange strain tensor consid-
þfdugT 0 Mz N0u N00v dxfv g fdv gT 0 Q y N0v N0u dxfug ering bending in plane xz are written according to expressions
R L T (5.37) and (5.38).
fdugT 0 Q y N0u N0v dxfv g !
2 2
ð5:32Þ 1 @u @w
gxx ¼ þ
2 @x @x
Considering a constant shear force, the bending moment and 0
!2 1
the shear force equations of a planar frame, as shown in Fig. 5.4, 2 2
1 @ @u @w @2w A @ 2 w @u
can be calculated by the expressions in (5.33). ¼ þ þy 2
y 2 ð5:37Þ
2 @x @x @x 2 @x @x
ðMZ1 þ M Z2 Þx ðMZ1 þ M Z2 Þ
M Z ¼ M Z1 þ Qy ¼ ð5:33Þ
L L @u @u @w @w @ 2 w @w @w @u
gxz ¼ þ ¼y 2 ð5:38Þ
@x @y @x @y @x @x @x @x
Substituting the complete shape functions and solving the inte-
grals of the problem leads to the geometric stiffness matrix, con- With the same development made in plane xy, the third integral
sidering the Euler-Bernoulli beam theory and higher-order terms from the virtual work expression in Eq. (5.8) can be written as
in the strain tensor. (5.39) for bending in plane xz.

RL
@u 2
R RL
@w 2
R RL @2 w
2 R
Plane xy (TBT) dU NL ¼ 12 0
d @x
dx t dA
A xx
þ 12 0
d @x
dx t dA
A xx
þ 12 0
d @x2
dx A
z2 txx dAþ
RL 2 R RL 2 R RL R
0
d @@xw2 @u
@x
dx t zdA
A xx
þ 0
d @@xw2 @w
@x
dx t zdA
A xz
0
d @w
@x
@u
@x
dx t dA
A xz
When the Timoshenko beam theory is considered, the nonlinear ð5:39Þ
part of the virtual work principle in Eq. (5.8) is written by Eqs. R R R
(5.34) and (5.35) applying the same relations presented for Applying relations A txx dA ¼ P; t zdA
A xx
¼ My ; t dA
A xz
¼ Q z , the
Euler-Bernoulli beam theory. expression becomes (5.40).
Z Z 1
dU NL ¼
dU NL ¼ sxx dgxx dV þ sxy dgxy dV 2
2 0
Z ! !2 13
2 2 2
V V
4Pd
L
@w I @u@ w
þ P d@ A5dx
y
þ
Z Z ! ! ! 0 @x A @x @x2
@v
2 2 2
L
1 @u @hz @hz @u
dU NL ¼ txx d
2 @x
þ
@x
þ y2
@x
y
@x @x
dx dA Z L" ! #
A 0 @ 2 w @u @w @u
Z Z L My d þ Q zd dx
@hz @u 0 @x2 @x @x @x
þ t xy d y hz hz dx dA
@x @x
A 0
ð5:40Þ
10
The axial part was already used for the xy plane. The expression R L 0 0 T R L I n on oT
dU NL ¼ fdwgT 0
P N w N w dxfwg þ fdwgT 0 P Ay N 0hy N0hy dxfwgþ
using the interpolation functions is given according to (5.41): RL n o RL n oT
T
R L T R L I T fdwgT 0 M y N 0hy N 0u dxfug fdugT 0 M y N 0u N 0hy dxfwgþ
dU NL ¼ fdwgT 0 P N 0w N 0w dxfwg þ fdwgT 0 P Ay N00w N 00w dxfwgþ R L T R L T
RL T RL T fdwgT 0 Q z Nhy N 0u dxfug fdugT 0 Q z N 0u N hy dxfwg
fdwgT 0 M y N 00w N 0u dxfug fdugT 0 M y N0u N 00w dxfwgþ
T RL 0 0 T T RL T ð5:47Þ
fdwg 0 Q z N w N u dxfug fdug 0 Q z N 0u N 0w dxfwg
ð5:41Þ
5.3.3. Combined torsion and axial force
where N w corresponds to the same N v interpolation functions. The interaction between torsion and axial force, as shown in
Additionally, considering a constant shear force, the bending Fig. 5.5, has an important influence on the geometric stiffness
moment and the shear force equations of a planar frame in plane matrix. The transverse displacements, v and w, need to take these
xz can be calculated by the expressions in (5.42). effects into account.

M y1 þ My2 x My1 þ M y2
M y ¼ My1 þ Qz ¼ ð5:42Þ eBBT
L L
Substituting the complete shape functions that were previously Considering the Euler-Bernoulli beam theory, the displacement
calculated and solving the integrals of the problem leads to the field is written according to expression (5.48)
geometric stiffness matrix, which considers the Euler-Bernoulli
@w @v
beam theory and higher-order terms in the strain tensor. u ¼ u0 z y v ¼ v 0 zhx w ¼ w0 þ yhx ð5:48Þ
@x @x
Plane xz (TBT) However, most terms of this equation have already been
employed when planes were analyzed independently. Therefore,
Meanwhile, for this beam theory, the nonlinear parts of the it is necessary to consider only terms that correspond to the rota-
strain tensor considering bending in plane xz are written by tion about the x axis, which is given in equation (5.49). Thus, the
expressions (5.43) and (5.44). Thus, the third integral from the vir- nonlinear parts of the Green-Lagrange strain tensor can be written
tual work expression is given according to (5.45) or with (5.46). according to expressions (5.50) to (5.52).
!
1@u
2
@w
2
v ¼ zhx w ¼ yhx ð5:49Þ
gxx ¼ þ
2 @x @x !

! @ v @hx
22
1 @hx @hx @w @hx
1 @u
2
@w
2
@hy
2
@u @hy gxx ¼ z þy z þy
¼ þ þ z2 z ð5:43Þ 2 @x @x @x @x @x @x
2 @x @x @x @x @x !

1 @hx 2 @ v @hx @w @hx
¼ z2 þ y 2 z þy ð5:50Þ
@u @u @w @w @hy @u 2 @x @x @x @x @x
gxz ¼ þ ¼z hy hy ð5:44Þ
@x @z @x @z @x @x
@u @u @ v @ v @w @w
RL R RL R RL @hy
2 R gxy ¼ þ þ
dU NL ¼ 12 0
d @u 2
@x
dx t dA
A xx
þ 12 0
d @w 2
@x
dx t dA
A xx
þ 12 0
d @x
dx A
z2 txx dAþ @x @y @x @y @x @y

RL
d
@hy @u
dx
R
t zdA þ
RL
d
@hy
hy dx
R
t zdA
RL
dhy @u dx
R
t dA @u @ v @2w @v @ 2 v @ v @w @hx
0 @x @x A xx 0 @x A xz 0 @x A xz
¼ þz 2 þy 2 þ hx þ y hx ð5:51Þ
ð5:45Þ @x @x @x @x @x @x @x @x
1 @u @u @ v @ v @w @w
dU NL ¼ gxz ¼ þ þ
2 @x @z @x @z @x @z
" ! !#
Z L
@u
Iy
2
@hy

@w
2 2 @u @w 2
@ w @w @ 2 v @w @ v @hx
Pd d þ þP dx ¼ þz 2 þy 2 hx þ z hx ð5:52Þ
@x
A @x @x @x @x @x @x @x @x @x @x
0
Z L The virtual work principle is then written as (5.53) and (5.54).
@hy @u @u Z Z
My d þ Q z d hy dx L
0 @x @x @x dU NL ¼ txx dgxx dx dA
ð5:46Þ A
Z
0
Z L " ! # !
@hx @ v @hx
2
1 2 @w @hx
Equation (5.47) can be written by considering only terms that ¼ t xx d z þ y2 z þy dx dA
A 0 2 @x @x @x @x @x
have not been used and applying interpolation functions.
ð5:53Þ
Fig. 5.5. Combined torsion and axial force (Mcguire et al., 2000).
11
Z Z L Z Z L @u @u @ v @ v @w @w
dU NL ¼ t xy dgxy dx dA þ txz dgxz dx dA gxy ¼ þ þ
A 0 A 0
@x @y @x @y @x @y
@u @hy @hz @w @hx
RL R RL R RL R ¼ hz þ z hz þ y hz þ hx þ y hx ð5:62Þ
dU NL ¼ 0
d @2 w @v
@x2 @x
dx t zdA
A xy
þ 0
d @w
h
@x x
dx A txy dA þ 0 d @h x
@x x
h dx A txy ydAþ @x @x @x @x @x
RL @ 2 v @w
R RL R R R
d @@xv hx dx A txz dA þ 0 d @h
L
þ
@u @u @ v @ v @w @w
0
d @x2 @x
dx t ydA
A xz 0
x
@x x
h dx A t xz zdA
ð5:54Þ gxz ¼ þ þ
@x @z @x @z @x @z
Adopting the relations presented before for the bending @u @hy @hz @v @hx
R ¼ hy þ z hy þ y hy hx þ z hx ð5:63Þ
moment and A z2 þ y2 dA ¼ J p , equation (5.53) becomes (5.55). @x @x @x @x @x
Expression (5.56) is written using interpolation functions. Equation (5.61) is the same for both beam theories because it is
Z Z L equal to (5.50). Finally, the virtual work principle for the
@ v @hx
2
1 PJp L @hx Timoshenko beam theory only changes in Eq. (5.64) due to the
dU NL ¼ d dx My d dx
2 A 0 @x 0 @x @x nonlinear distortion.
Z L Z Z Z Z
@w @hx L L
þ ðMz Þ d dx ð5:55Þ dU NL ¼ t xy dgxy dx dA þ t xz dgxz dx dA
0 @x @x
A 0 A 0
PJ p 0 0 T
RL RL T
dU NL ¼ fdhx gT A
Nhx N hx dxfhx g fdv gT 0 M y N 0v N 0hx dxfhx g
0 dU NL ¼
RL
d
@hy R RL
hz dx A txy zdA þ 0 d @w
R RL
h dx A t xy dA þ 0 d @h x
R
h dx A t xy ydAþ
@x x @x x
R R T 0 @x
fdhx gT 0 M y N 0hx N 0v dxfv g fdwgT 0 M z N 0w N 0hx dxfhx g

L T L RL R R R R R
h dx A txz ydA 0 d @@xv hx dx A txz dA þ 0 d @h
L L
þ d @h z x
h dx A txz zdA
R T
0 @x y @x x
fdhx gT 0 M z N0hx N 0w dxfwg

L
ð5:64Þ
ð5:56Þ
Expression (5.64) is rewritten according to (5.65) or (5.66) using
For bisymmetric sections, the shear force and torsion moment interpolation functions with the same relations presented before in
are defined by (5.57). Thus, expression (5.54) can be written (5.57).
according to (5.58). RL @hy RL
dU NL ¼ hz dx ða 1ÞMx þ 0 d @w
d @x
h dx Q y þ
@x x
Z Z Z 0
ð5:65Þ
RL RL
t xy dA ¼ Q y ; t xz dA ¼ Q z ; t xz ydA ¼ aMx ; þ 0 d @hz
@x y
h dx aM x 0 d @@xv hx dx Q z
ZA A A
t xy zdA ¼ ða 1ÞM x ð5:57Þ T R L M x n 0 o RL n oT

A dU NL ¼ dhy 0 2
N hy fNhz gT dxfhz g fdhz gT 0 M2x fN hz g N 0hy dx hy
R L T T R L M x 0 T
RL RL þfdhz gT 0 M2x N 0hz N hy dx hy þ dhy N hy N hz dxfhz g
d @@xw2 @@xv dx ða 1ÞM x þ 0 d @w
2 0 2
dU NL ¼ h dx Q y
@x x
RL RL T
0
ð5:58Þ þfdwgT 0 Q y N0w fN hx gT dxfhx g þ fdhx gT 0 Q y fNhx g N0w dxfwg
RL RL R R T
þ 0 d @@xv2 @w dx aM x 0 d @@xv hx dx Q z
2
fdv gT 0 Q z N 0v fN hx gT dxfhx g fdhx gT 0 Q z fNhx g N0v dxfv g
L L
@x
ð5:66Þ
Considering Saint Venant pure torsion, M x ¼ Mx2 and a ¼ 1=2.
Then, Eq. (5.59) is written by rewriting (5.58) using interpolation
5.3.4. Finite rotations
functions.
In spatial elements, the geometric stiffness matrix needs to con-
RL 00 0 T RL T
dU NL ¼ fdwgT 0
Mx
2
N w N v dxfv g fdv gT 0 Mx
2
N 0v N 00w dxfwg sider finite rotations to satisfy the kinematic compatibility at the
RL T RL joint of angled elements (Mcguire et al., 2000; Conci, 1988). Thus,
þfdv g T
0
Mx
2
Nv 00
N0w dxfwg þ fdwg T
0
Mx
2
N0w Nv 00 T
dxfv g for a spatial vector, the rotation is given by equation (5.67), as
RL RL T shown in Fig. 5.6.
þfdwgT 0 Q y N 0w fN hx gT dxfhx g þ fdhx gT 0 Q y fN hx g N 0w dxfwg
RL RL T
V 1 ¼ RðhÞV 0 ð5:67Þ
fdv gT 0
Q z N0v fN hx gT dxfhx g fdhx gT 0 Q z fN hx g N 0v dxfv g
ð5:59Þ
The combined torsion and axial force contribution is found by
substituting the complete shape functions and solving the
integrals.
TBT
When the Timoshenko beam theory is considered, the displace-

ment field cannot be written using the transverse displacement
derivative and instead should be given by Eq. (5.60). Thus, the non-
linear part of the strain tensor that is not used before is written
according to expressions (5.61) to (5.63).
u ¼ u0 zhy yhz v ¼ v 0 zhx w ¼ w0 þ yhx ð5:60Þ
!
@ v @hx
2 2
1 @hx @hx @w @hx
gxx ¼ z þy z þy
2 @x @x @x @x @x @x

@hx @ v @hx
2
1 2 @w @hx
¼ z þ y2 z þy ð5:61Þ
2 @x @x @x @x @x Fig. 5.6. Spatial transformation between two vectors (Aguiar et al., 2014).
12
Z Z Z Z L
R¼Iþ
senðhÞ
W ðhÞ þ
1 cosðhÞ
W ðhÞ2
L
@u @ v @u
h 2 ¼ txx d dx dA þ t xy d þ dx dA
2 3h 8 9 A 0 @x A 0 @x @y
0 hz hy >
< hx >
= Z Z L
6
W ðhÞ ¼ 4 hz 0
7
hx 5; h ¼ hy @w @u
> þ txz d þ dx dA
: > ; @x @z
hy hx 0 hz A 0
Using a trigonometric series approximation given by expression Z Z Z

t xx @ hx hy
L
@ ðhx hz Þ
(5.68) and considering only terms up to the second order, the rota- ! t ij Deij ¼ d yþd z dx dA
V A 0 2 @x @x
tion matrix can be rewritten according to (5.69). Z Z L
t xy @ hy hz
h3 h5 h7 þ d hx hy þ d z dx dA
senðhÞ ¼ h þ þ A 0 2 @x
3! 5! 7! Z Z L
2 t xz @ hy hz
h h4 h6 þ dðhx hz Þ þ d y dx dA ð5:73Þ
cosðhÞ ¼ 1 þ þ ð5:68Þ A 0 2 @x
2! 4! 6!
Equation (5.74) can be written by employing the relations for
2
1 the bending moment, torsion, shear and axial force.
RðhÞ ¼ I þ W ðhÞ þ W ðhÞ ð5:69Þ
2 R R L z @ ðhx hy Þ RL M RL Q
Therefore, considering a bisymmetric section as shown in V
tij Deij ¼ 0 M 2
d @x dx þ 0 2y d @ðh@xx hz Þ dx þ 0 2y d hx hy dx
Fig. 5.7, the position vector of a point, fag, and after rotation,fbg, RL @ ðhy hz Þ RL R L x @ ðhy hz Þ
þ 0 ða1ÞM2
x
d @x dx þ 0 Q2z dðhx hz Þdx þ 0 aM 2
d @x dx
is written by equation (5.70).
" # ð5:74Þ
2
1
fbg ¼ I þ W ðhÞ þ W ðhÞ fag ð5:70Þ Expression (5.75) is obtained by adopting a ¼ 1=2 for Saint
2 Venant torsion and using integration by parts. Finally, the finite
rotation matrix is given in (5.76) and is independent of the beam
The displacement can be written as: fbg fag or fug. However,
theory considered.
W ðhÞ has already been taken into account in the previous analysis
with the principle of virtual work. Thus, the displacement can be R L M
tij Deij ¼ M z
d hx hy 0 þ 2y dðhx hz ÞjL0
written only as a function of W ðhÞ , according to expressions 2 RV 2

t Deij ¼ 12 M z1 d hx1 hy1 þ M y1 dðhx1 hz1 Þ M z2 d hx2 hy2 þ M y2 dðhx2 hz2 Þ
V ij
(5.71) and (5.72). ð5:75Þ
2 3
8 9 hy 2 þhz 2 hx hy 8 9
> hx hz
2 3
6 <u>
= 2 2 2 >0>
7< = 0 M z1 M y1 0 0 0 hx1
1 6 7
W ðhÞ2 fag ) v ¼ 6
2
fug ¼ þhz 2
hx hy
hx
hy hz
7 y 6 M z1
2 >
: > ; 4 2 2 2 5>
: > ; 6 0 0 0 0 0 7 7 y1
h
w hy hz h 2 þh 2 z 6 7
hx hz
x y 16 My1 0 0 0 0 0 7 hz1
¼ 6 7
2 2 2
K g;Rotfin
ð5:71Þ 26
6 0 0 0 0 M z1 M y2 7
7 hx2
6 7
! 4 0 0 0 Mz1 0 0 5 hy2
hx :hy :y hx :hz :z hy :hz :z h2 h2
u¼ þ v ¼ y x þ z 0 0 0 My2 0 0 hz2
2 2 2 2 2 ð5:76Þ
2
!
hy :hz :y h2 hy Therefore, considering these effects, the local complete tangent
w¼ z x þ ð5:72Þ
2 2 2 stiffness matrix can be written by adding all the components for
both bending theories. The results are presented in the work by
To consider the finite rotations, the second integral of the virtual Rodrigues (Rodrigues, 2019) and with an open source code avail-
work principle needs to be used according to expression (5.73). able in the files StfBeamEulerBernoulliT (positive axial force) and
Z Z Z L Z Z L StfBeamEulerBernoulliC (negative axial force) considering Euler-
tij Deij ¼ t xx dexx dx dA þ t xy dexy dx dA Bernoulli beam theory. For the Timoshenko beam theory, the files
V A 0 A 0
Z Z L are StfBeamTimoshenkoT (positive axial force) and StfBeamTi-
þ t xz dexz dx dA ¼ moshenkoC (negative axial force) (Rodrigues et al., 2020, 2021).
A 0 To avoid numerical instability, the local tangent stiffness matrix
can be written using a Taylor series expansion, providing more
terms for the tangent stiffness matrix than the usual formulations.
These approximations provide greater precision in the results, and
in this work, tangent matrices with up to 3 and 4 terms were
developed. It is important to note that the 2-term approximation
provides the usual elastic and tangent stiffness matrices. These
matrices are also presented in the work by Rodrigues (Rodrigues,
2019) and implemented in (Rodrigues et al., 2020, 2021) in the files
GeoStfBeamEulerBernoulli and GeoStfBeamTimoshenko.
6. Numerical applications
In this section, numerical examples were developed with no

discretization to test the formulation from this work, the complete
formulation and the Taylor series expansion with 3 and 4 terms.
The elements can be identified based on the following
Fig. 5.7. Spatial transformation between two vectors (Mcguire et al., 2000). descriptions:
13
EBBT_Large_Complete – Euler-Bernoulli beam theory, hyper- 6.1. Isolated columns

bolic and geometric functions in the tangent stiffness matrix, and
higher order terms in the strain tensor. To verify the developed tangent stiffness matrix, the buckling
EBBT_Large_4tr – Euler-Bernoulli beam theory, 4 terms in the load of columns was studied as shown in Fig. 6.1, adopting just
tangent stiffness matrix (1 elastic + 3 geometric), and higher order one element per member. Additionally, a reduced slenderness ratio
terms in the strain tensor. (k ¼ L=h) was employed because of the influence of Timoshenko
EBBT_Large_3tr – Euler-Bernoulli beam theory, 3 terms in the beam theory. The columns have a length of L = 1 m, Young’s mod-
tangent stiffness matrix (1 elastic + 2 geometric), and higher order ulus of E = 107 kN/m2, section form factor of v ¼ 1 and null Pois-
terms in the strain tensor. son’s ratio (m ¼ 0).
EBBT_Large or Small_2tr – Euler-Bernoulli beam theory, 2 Fig. 6.2 shows the equilibrium paths for the clamped column,
terms in the tangent stiffness matrix (1 elastic + 1 geometric), while Figs. 6.3 and 6.4 represent equilibrium paths for the simply
and higher order terms in the strain tensor (Large) or (Small) if this supported and fixed and simply supported columns, respectively,
influence is not considered (Yang and Leu, 1994; Yang and Kuo, with a slenderness ratio of k ¼ 10:0 for the Euler-Bernoulli beam
1994; Chen, 1994). This formulation is the most conventional. theory and k ¼ 4:0 for the Timoshenko beam theory. The numerical
TBT_Large_Complete – Timoshenko beam theory, hyperbolic results for the buckling loads for the columns are shown in Table 1.
and geometric functions in the tangent stiffness matrix, and higher Analyzing the equilibrium paths, it can be noted that the complete
order terms in the strain tensor. formulation developed in this work, where the tangent stiffness
TBT_Large_4tr – Timoshenko beam theory, 4 terms in the tan- matrix was calculated with the complete interpolation functions
gent stiffness matrix (1 elastic + 3 geometric), and higher order (EBBT_Large_complete and TBT_Large_complete), provides the cor-
terms in the strain tensor. rect prediction for the buckling loads of columns using just one ele-
TBT_Large_3tr – Timoshenko beam theory, 3 terms in the tan- ment per member for both beam theories. With no discretization,
gent stiffness matrix (1 elastic + 2 geometric), and higher order the complete formulation provides the best approximation for the
terms in strain tensor. analytical response for both beam theories. In Table 1 the results
TBT_Large or Small_2tr – Timoshenko beam theory, 2 terms in can be seen numerically and the efficiency of the complete formula-
the tangent stiffness matrix (1 elastic + 1 geometric), and higher tion is observed when compared with literature values.
order terms in the strain tensor (Large) or (Small) if this influence Additionally, it can be seen that geometric matrices obtained
is not considered (Rodrigues et al., 2019). This formulation is the from a Taylor series approximation (Large_3tr and Large_4tr) also
most conventional. provide an accurate prediction for the buckling load. The approxi-
The results were compared with the results obtained with Mas- mation using 4 terms in some cases provides the same response
tan2 v3.5. The software is able to perform a geometric nonlinear from the complete formulation, and the approximation with 3
analysis considering both beam theories while considering cubic terms follows those curves closely. Table 1 also shows the proxim-
interpolation functions and disregarding higher order terms in ity of the complete formulation and the Taylor series expansion
the strain tensor. The element tangent stiffness matrix considered and the differences with the usual formulation.
in Mastan2 v3.5 is described in the work by McGuire et al.
(Mcguire et al., 2000). In this research, Mastan2 v3.5 elements 6.2. Continuous beam-column
are labeled as follows:
EBBT_Mastan – Mastan2 v3.5 software Euler-Bernoulli beam The next example studies a continuous beam-column; this prob-
theory; lem is presented in the work by (Timoshenko and Gere (1963) and
TBT_Mastan – Mastan2 v3.5 software Timoshenko beam is shown in Fig. 6.5. The geometry of the structure, material and sec-
theory. tion proprieties are the same as in the first example with a length of
In cases that are no available analytical solution, the reference L = 1 m, Young’s modulus of E = 107 kN/m2, section form factor of
response for comparison purposes was given by a numerical solu- v ¼ 1 and null Poisson’s ratio (m ¼ 0). Fig. 6.6 shows the equilibrium
tion with discretized structure (EBBT or TBT_Large_2tr_ number of path considering a slenderness ratio of k ¼ 10:0 for the Euler-
elements in each member). Bernoulli beam theory and k ¼ 4:0 for Timsohenko beam theory.
Euler Critical Load – Analytical Euler buckling load The structures were modeled with one element in each span.
Timoshenko and Gere – Analytical Timoshenko buckling load Again, using just one element per member, the complete formu-
(Timoshenko and Gere, 1963) lation and the Taylor series expansion with 4 terms provide the
Fig. 6.1. Analyzed columns (Silva et al., 2016).
14
Fig. 6.2. Equilibrium paths for a clamped column.
Fig. 6.3. Equilibrium paths for a simply supported column.
Fig. 6.4. Equilibrium paths for a fixed and simply supported column.
correct buckling load of the structure. The expansion with just 3 has a length of L = 1 m, Young’s modulus of E = 107 kN/m2, sec-
terms gives approximate results. These results can be seen numer- tion form factor of v ¼ 5=6 and Poisson’s ratio of m ¼ 0:3. The
ically in Table 2. frame is loaded by vertical loads (P) and two small lateral dis-
The conventional geometric stiffness matrix and the Mastan2 turbing loads (H ¼ 0:001P). The structure was modeled with one
v3.5 software cannot predict the buckling load of the continuous element per member, and the response was compared with a dis-
beam-column, providing differences on the order of 50% to the cretized structure using the usual formulation. The equilibrium
expected result with a discretized structure. paths can be seen in Fig. 6.8, considering k ¼ 10:0 for the Euler-
Bernoulli beam theory and k ¼ 4:0 for the Timoshenko beam
6.3. Spatial frame theory.
Once again, the complete formulation and the Taylor series
The example presented is the first to evaluate the developed expansion with 4 terms provide the buckling load of the spatial
formulation in a spatial structure, according to Fig. 6.7. The frame frame. The expansion with 3 terms overlaps these results for the
15
Table 1
Numerical Buckling Loads for Columns.
Euler Bernoulli Beam Theory - EBBT

Discretization 1 Element Euler Critical Load
EBBT- Element Mastan Small Large Large Large Large
2tr 2tr 3tr 4tr Complete
Clamped col. (k = 10.0) 2.496 2.491 2.486 2.469 2.468 2.468 2.467
Simply supported col. (k = 10.0) 12.255 12.312 12.197 10.348 10.047 9.95 9.869
Fixed – simply sup. col. (k = 10.0) 31.651 31.630 30.798 23.583 21.779 20.541 20.142
Timoshenko Beam Theory – TBT
Discretization 1 Element (Timoshenko and Gere, 1963)
TBT - Element Mastan Small Large Large Large Large
2tr 2tr 3tr 4tr Complete
Clamped col. (k = 4.0) 2.475 2.475 2.444 2.419 2.417 2.417 2.408
Simply supported col. (k = 4.0) 13.914 13.940 12.866 10.064 9.546 9.368 8.949
Fixed - simply sup. col. (k = 4.0) 40.642 > 60 38.065 22.572 19.892 17.826 16.401
Fig. 6.5. Beam-column (Timoshenko and Gere, 1963).
Fig. 6.6. Equilibrium paths for a beam-column.
Euler-Bernoulli beam theory and provides intermediate results for 6.4. Roorda spatial frame
the Timoshenko theory. The usual formulation and Mastan2 v3.5
software provide more elevated buckling loads. Table 3 presents Using the same bar length, material and section proprieties as
these conclusions numerically. the structure presented before, another spatial frame was studied
Table 2
Numerical buckling loads for the beam-columns.
Discretization 1 Element (Timoshenko and Gere, 1963)

Element Mastan Small Large Large Large Large
2tr 2tr 3tr 4tr complete
EBBT -k = 10.0 12.240 12.310 12.240 10.453 10.145 9.953 9.869
TBT - k = 4.0 13.989 13.843 12.867 10.074 9.553 9.368 8.949
16
Fig. 6.9. Roorda spatial frame.
3 terms. These formulations give better approximations when

Fig. 6.7. Spatial frame.
compared with the usual formulations with 2 terms, which present
relevant errors. Table 4 shows the numerical results.
When considering just one element per member, Table 4 shows
(Fig. 6.9). The frame is loaded by a vertical load (P) and by two
that the conventional formulation doubles the buckling load, while
small lateral disturbing loads (H ¼ 0:001P). The structure was
the complete formulation maintains the error at a rate of 12% for
modeled with one element per member, and the response was
TBT elements and 8% for EBBT theory.
compared with a discretized structure using the usual formulation.
The equilibrium paths can be seen in Fig. 6.10 with k ¼ 10:0 for the
Euler-Bernoulli beam theory and k ¼ 4:0 for the Timoshenko beam 6.5. Spatial frame with inclined columns
theory.
There is no available analytical solution for the critical load of The frame presented in Fig. 6.11 is loaded by two vertical loads
this frame. In this case, the reference solution for comparison pur- (P) and by small lateral disturbing loads (H ¼ 0:001P) with a length
poses was the Large_2tr formulation with 4 segments in each of L = 1 m, Young’s modulus of E = 107 kN/m2, section form factor of
member. One may observe that, with no discretization, the com- v ¼ 5=6 and Poisson’s ratio of m ¼ 0:3.
plete formulation provides the best approximation for the refer- The equilibrium paths found for the structure employing one
ence response. The second-best solution is for the Taylor series element per member and both beam theories, considering
expansion with 4 terms and then the Taylor series expansion with k ¼ 4:0, are given in Fig. 6.12 and numerically presented in Table 5.
Fig. 6.8. Equilibrium paths for a spatial frame.
Table 3
Numerical buckling loads for a spatial frame.
Discretization 1 Element 4 Elements

Element Mastan Small Large Large Large Large Large
2tr 2tr 3tr 4tr Complete 2tr_4el
EBBT -k = 10.0 5.577 5.569 5.554 5.507 5.501 5.500 5.488
TBT-k = 4.0 5.068 5.115 5.030 4.977 4.837 4.934 4.927
17
Fig. 6.10. Equilibrium paths for a Roorda spatial frame.
Table 4
Numerical buckling loads for a Roorda spatial frame.

EBBT -k = 10.0 19.008 24.461 20.897 15.262 14.484 14.136 13.200
TBT - k = 4.0 21.135 21.166 19.630 13.670 12.445 11.753 10.414
form factor of v ¼ 5=6 and Poisson’s ratio of m ¼ 0:3 and is loaded

by four vertical loads (P) and by two lateral loads (H ¼ 0:01P)
according to Fig. 6.13. Fig. 6.14 shows the deformed structure illus-
trating the torsion experienced by the space frame with the scale
factor increased tenfold. Finally, the equilibrium path of the struc-
ture considering k ¼ 10 for the Euler-Bernoulli beam theory and
k ¼ 4:0 for the Timoshenko beam theory is illustrated.Fig. 6.15.
The solution obtained by the complete formulation and for the
Taylor series expansion provides the correct buckling load employ-
ing just one element per member for both beam theories, while the
usual formulation proceeds to higher values for the critical load.
This is clearly illustrated numerically in Table 6.
6.7. Asymmetric spatial frame
The last example studies an asymmetric frame with moderate

slenderness, k ¼ 6:6 for TBT, as shown in Fig. 6.16, loaded by two
vertical loads (P) and small lateral disturbing loads (H ¼ 0:001P),
with a length of L = 1 m, Young’s modulus of E = 107 kN/m2, section
form factor of v ¼ 5=6 and Poisson’s ratio of m ¼ 0:3. Fig. 6.17
shows the influence of torsion on this structure, and Fig. 6.18
Fig. 6.11. Spatial frame with inclined columns (adapted from Zugic et al. (Zugic shows the equilibrium paths.
et al., 2016). Finally, the complete formulation and the Taylor series expan-
sion provide the buckling load for the spatial frame, with the
As concluded for the other examples, the complete formulation curves overlapping. The usual formulation and Mastan software
calculates the exact buckling load of the structure. In this example, overcome this load for both beam theories. Table 7 shows these
the Taylor series expansion with 4 or even 3 terms reaches this conclusions numerically.
load as well, and their results overlap. Although the usual formula-
tion gives a good approximation with just one element, employing
the formulation developed in this work improves the results. In 7. Conclusions
addition, consideration of Timoshenko beam theory reduces the
critical buckling load for structures with small slenderness. This research developed a complete formulation of the tangent
stiffness matrices of a spatial frame element for the Euler-Bernoulli
and Timoshenko beam theories, considering interpolation func-
6.6. Spatial frame with torsion tions obtained directly from the solution of the equilibrium differ-
ential equation of a deformed infinitesimal element, which
This example explores the influence of torsion. The frame also includes the influence of axial forces. Additionally, the formulation
has a length of L = 1 m, Young’s modulus of E = 107 kN/m2, section uses higher-order terms in the strain tensor, and the geometric
18
Fig. 6.12. Equilibrium paths for the spatial frame with inclined columns.
stiffness matrix is corrected by employing finite rotations. The

complete formulation was rewritten using a Taylor series expan-
sion considering 3 and 4 terms, improving the usual tangent stiff-
ness matrices for both beam theories that involve only 2 terms (1
elastic and 1 geometric).
The formulation was evaluated in several numerical tests for
planar and spatial structures to verify the geometric nonlinear
response of the structure considering distinct bending theories
and the consideration of high-order terms for the strain. The first
set of examples employs the proposed element to analyze isolated
columns with different boundary conditions and a continuous
beam-column. Similarly, this formulation was evaluated for spatial
structures, and some examples explored the torsion effect, apply-
ing a load or geometric asymmetry.
The complete formulation presented precisely predicted results
for the critical load of the developed examples with a minimum
discretization of the structure, while the usual tangent stiffness
Fig. 6.13. Spatial frame with asymmetric loads.
matrix would require a more refined mesh to predict this behavior.
Table 5
Numerical buckling loads for the spatial frame with inclined columns.

EBBT-k = 4.0 7.070 7.179 7.157 6.952 6.942 6.942 6.919
TBT- k = 4.0 6.182 6.371 6.227 6.155 6.155 6.082 6.082
Fig. 6.14. Deformed spatial frame with torsion.
19
Fig. 6.15. Equilibrium paths for the spatial frame with torsion.
Table 6
Numerical buckling loads for the spatial frame with torsion.

EBBT- k = 10.0 7.488 7.328 7.313 7.259 7.257 7.256 7.260
TBT-k = 4.0 6.389 6.548 6.374 6.367 6.367 6.341 6.341
Fig. 6.17. Deformed asymmetric spatial frame.
Due to the numerical instability that can occur using the hyper-
bolic and trigonometric expressions of the complete formulation,
an alternative geometric stiffness matrix using a Taylor series
expansion with 4 terms is proposed to solve this problem with
accurate results. All necessary equations for this implementation
Fig. 6.16. Asymmetric spatial frame with inclined columns (adapted from Zugic are available online in open source.
et al. (Zugic et al., 2016).
Although precise buckling loads are achieved using just one ele-
ment if one is interested in following the actual equilibrium path
up to the critical load, a refined discretization should be used. This
As expected, the examples have also shown that for structures with
is suggested because one element is not able to capture local
small slenderness, the consideration of the Timoshenko beam the-
effects, such as the P-‘‘small delta”.
ory leads to lower buckling loads, and the proposed formulation
The continuation of this work will focus on monitoring the
illustrated this effect. The Timoshenko beam theory is also used
equilibrium path in the postcritical stage. For this, more improved
in cases of structures with a small shear-to-bending ratio, regard-
solution schemes for incremental geometric nonlinear analysis
less of their slenderness.
should be employed.
The reported results clearly illustrate the efficiency of the com-
plete formulation and the Taylor series expansion using 4 terms in
the solution of nonlinear geometric analyses using only one ele-
ment for planar and spatial structures. Using 3 terms also provides
better results than considering the usual geometric stiffness 8. Funding sources
matrix, which overestimates the buckling load.
The presented formulation solves geometric nonlinear prob- This work has been partially supported by Conselho Nacional de
lems and can be implemented in any structural analysis software. Desenvolvimento Científico e Tecnológico (CNPq) and FAPERJ.
20
Fig. 6.18. Equilibrium paths for an asymmetric spatial frame.
Table 7
Numerical buckling loads for asymmetric spatial frames.

EBBT- k = 10 7.584 7.606 7.586 7.499 7.486 7.483 7.452
TBT- k = 6.6 7.282 7.298 7.259 7.152 7.132 7.083 7.089
Declaration of Competing Interest Chiwanga, M., Valsangkar, A.J., 1988. Generalized beam element on two-parameter
elastic foundation. J. Struct. Eng. 114 (6), 1414–1427.
Conci, A., 1988. Analysis of reticulated steel structures with warping consideration and
The authors declare that they have no known competing finan- geometric and material nonlinearities (in Portuguese). PUC, Rio de Janeiro, Brazil.
cial interests or personal relationships that could have appeared Crisfield, M.A., 1991. Non-linear finite element analysis of solids and structures.
John Wiley & Sons Inc, NY, USA.
to influence the work reported in this paper.
Davis, R., Henshell, R.D., Warburton, G.B., 1972. A Timoshenko beam element. J.
Sound Vib. 22 (4), 475–487.
References Eisenberger, M., Yankelevsky, D.Z., 1985. Exact stiffness matrix for beams on elastic
foundation. Comput. Struct. 21 (6), 1355–1359.
Friedman, Z., Kosmatka, J.B., 1993. An improved two-node Timoshenko beam finite
Aguiar, F.L.L., Almeida, C.A., Paulino, G.H., 2014. A three-dimensional multilayered
element. Comput. Struct. 47 (3), 473–481.
pipe beam element: nonlinear analysis. Comput. Struct. 138, 142–161.
Girhammar, U.A., Gopu, V.K.A., 1993. Composite beam-columns with interlayer slip
Areiza-Hurtado, M., Vega-Posada, C., Aristizábal-Ochoa, J.D., 2005. Second-order
– exact analysis. J. Struct. Eng. 119 (4), 1265–1282.
stiffness matrix and loading vector of a beam-column with semirigid
Goto, Y., Chen, W.F., 1987. Second-order elastic analysis for frame design. J. Struct.
connections on an elastic foundation. J. Eng. Mech. 131 (7), 752–762.
Eng. 113 (7), 1501–1519.
Aristizábal-Ochoa, J.D., 1997. First- and second-order stiffness matrices and load
Ha, K.H., 1993. Stiffness matrix for exact solution of sandwich beam and frame
vector of beam-columns with semirigid connections. J. Struct. Eng. 123 (5),
systems. J. Struct. Eng. 119 (4), 1150–1167.
669–678.
Martha, L.F., 1999. Ftool: A structural educational interactive tool. In: In:
Aristizábal-Ochoa, J.D., 2007. Tension and compression stability and second-order
Proceedings of Workshop in Multimedia Computer Techniques in Engineering
analyses of three-dimensional multicolumn systems: effects of shear
Education, Institute for Structural Analysis. Technical University of Graz,
deformations. J. Eng. Mech. 133 (1), 106–116.
Austria, pp. 51–65.
Aristizábal-Ochoa, J.D., 2008. Slope-deflection equations for stability and second-
Martha, L.F., 2018. Matrix structural analysis with object oriented programming (in
order analysis of Timoshenko beam-column structures with semi-rigid
Portuguese). Elsevier, Rio de Janeiro.
connections. Eng. Struct. 30, 2517–2527.
Martha L. F., Burgos R. B., Alternative forms of considering shear deformation in
Aristizábal-Ochoa, J.D., 2012. Matrix method for stability and second-order analysis
axially loaded timoshenko beams (in Portuguese). In: XXXVI South American
of Timoshenko beam-column structures with semi-rigid connections. Eng.
Structural Engineering Journeys (ASAEE), Montevidéu, 2014.
Struct. 34, 289–302.
Martha, L.F., Burgos, R.B., 2015. Possible inconsistencies when considering shear
Aydogan, M., 1995. Stiffness-matrix formulation of beams with shear effect on
deformation in axially loaded beams (in Portuguese). In: 57th CBC (Brazilian
elastic foundation. J. Struct. Eng. 121 (9), 1265–1270.
Conference on Concrete), Bonito, Brazil.
Balling, R.J., Lyon, J.W., 2011. Second-order analysis of plane frames with one
Martha, L.F., Parente Junior, E., 2002. An object-oriented framework for finite
element per membe. J. Struct. Eng. 137 (11), 1350–1358.
element programming. in: World Congr. Comput. Mech., Vienna 01–10.
Bathe, K.J., 1996. Finite element procedures. Prentice-Hall, Englewoods Cliffs, NJ,
Mcguire, W., Gallagher, R.H., Ziemian, R.D., 2000. Matrix structural analysis. John
USA.
Wiley & Sons Inc, NY, USA.
Bathe, K.J., Bolourchi, S., 1979. Large displacement analysis of three-dimensional
Morfidis, K., 2007. Exact matrices for beams on three-parameter elastic foundation.
beam structures. Int. J. Numer. Meth. Eng. 14, 961–986.
Comput. Struct. 85, 1243–1256.
Battini, J.M., 2002. Co-rotational beam elements in instability problems. Royal
Morfidis, K., Avramidis, I.E., 2002. Formulation of a generalized beam element on a
Institute of Technology, Stockholm, Sweden.
two-parameter elastic foundation with semi-rigid connections and rigid offsets.
Burgos, R.B., Martha, L.F., 2013. Exact shape functions and tangent stiffness matrix
Comput. Struct. 80, 1919–1934.
for the buckling of beam-columns considering shear deformation. XXXIV
Nukulchai, W.K., Dayawansa, P.H., Karasudhi, P., 1981. An exact finite element
Iberian Latin American Congres on Computational Methods in Engineering.
model for deep beams. Int. J. Struct. 1 (1), 1–7.
Burgos, R.B., Silva, R.R., Martha, L.F., 2005. Evaluation of Critical Loads and Initial
Onu, G., 2000. Shear effect in beam finite element on two-parameter elastic
PostBuckling Behavior of Portal Frames (in Portuguese). In: XXXV Iberian Latin
foundation. J. Struct. Eng. 126 (9), 1104–1107.
American Congres on Computational Methods in Engineering, Guarapari, Brazil.
Onu, G., 2008. Finite elements on generalized elastic foundation in Timoshenko
Chan, S.L., Gu, J.X., 2000. Exact tangent stiffness for imperfect beam-column
beam theory. J. Eng. Mech. 134 (9), 763–776.
members. J. Struct. Eng. 126 (9), 1094–1102.
Pacoste, C., Eriksson, A., 1995. Element behavior in post-critical plane frames
Chen, D.C., 1994. Geometric nonlinear analysis of three-dimensional structures.
analysis. Comput. Methods Appl. Mech. Eng. 125, 319–343.
Cornell University, Ithaca, NY.
Pacoste, C., Eriksson, A., 1997. Beam elements in instability problems. Comput.
Chen, W.F., Lui, E.M., 1991. Stability design of steel frames. CRC Press, Boca Raton,
Methods Appl. Mech. Eng. 144, 163–197.
USA.
21
Pilkey, W.D., Kang, W., Schramm, U., 1995. New structural matrices for a beam Marcos Antonio Campos Rodrigues. Adjunct professor in
element with shear deformation. Finite Elem. Anal. Des. 19, 25–44. the Department of Civil Engineering at Federal Univer-
Reddy, J.N., 1997. On locking-free shear deformable beam finite elements. Comput. sity of Espírito Santo (UFES). B.Sc. in Civil Engineering
Methods Appl. Mech. Eng. 149, 113–132. from Federal University of Espírito Santo (UFES) in 2011.
Rodrigues, M.A.C., 2019. Integrated solutions for the formulations of the geometric M.Sc. in Aeronautical and Mechanical Engineering from
nonlinearity problem (in Portuguese). PUC, Rio de Janeiro, Brazil. the Aeronautics Institute of Technology (ITA) in 2014
Rodrigues, M.A.C., Burgos, R.B., Martha, L.F., 2019. A unified approach to the and Ph.D. in Civil Engineering from Pontifical Catholic
Timoshenko geometric stiffness matrix considering higher-order terms in the University of Rio de Janeiro (PUC-Rio) in 2019. He has
strain tensor. Latin Am. J. Solids Struct. 16, 1–22. experience in computational mechanics and structural
Rodrigues, M.A.C., Burgos, R.B., Martha, L.F., 2020. CENLG – Complete expressions
design. His main research interests include Finite Ele-
for geometric non linear analysis. GitLab Git-repository, Project ID 19532538
ment Method, Nonlinear Analysis, Design of Reinforced
https://gitlab.com/marcos.a.rodrigues/cenlg-complete-expressions-for-
geometric-non-linear-analysis. Concrete, Prestressed and Steel Structures.
Rodrigues, M.A.C., Burgos, R.B., Martha, L.F., 2021. CENLG – Complete expressions
for geometric non linear analysis. File Exchange, MathWorks. https://
www.mathworks.com/matlabcentral/fileexchange/77380-cenlg-complete-
expressions-for-geometric-non-linear-analys. Rodrigo Bird Burgos Associate Professor in the Depart-
Rodrigues, C.F., Suzuki, J.L., Bittencourt, M.L., 2016. Construction of minimum ment of Structures and Foundations at the State
energy high-order Helmholtz bases for structured elements. J. Comput. Phys. University of Rio de Janeiro (UERJ). B.Sc. (2003), M.Sc.
306, 269–290. (2005) and Ph.D. (2009) in Civil Engineering from Pon-
Santana, M.V.B., Silveira, R.A.M., 2019. Numerical fundamentals and interactive tifical Catholic University of Rio de Janeiro (PUC-Rio).
computer graphics system for the nonlinear analysis of planar frames. REM, Int. Postdoctoral researcher at Pontifical Catholic University
Eng. J. 72(2), 199–207.
of Rio de Janeiro from 2009 to 2012. He has experience
Schramm, U., Kitis, L., Kang, W., Pilkey, W.D., 1994. On the shear deformation
in computational modeling and advanced numerical
coefficient in beam theory. Finite Elem. Anal. Des. 16, 141–162.
methods for solving differential equations. His main
Shirima, L.M., Giger, M.W., 1992. Timoshenko beam element resting on two-
parameter elastic foundation. J. Eng. Mech. 118 (2), 280–295. research interests include instability and structural
Silva, J.L., Lemes, I.J.M, Silveira, R.A.M, Silva, A.R.D., 2016. Influence of beam theory dynamics, wavelet functions, interpolets, Wavelet-
on geometrically nonlinear analysis of reticulated structures (in Portuguese). In: Galerkin Method, Finite Element Method.
XXXVII Iberian Latin American Congress on Computational Methods in
Engineering, Brasilia, Brazil.
So, A.K.W., Chan, S.L., 1991. Buckling and geometrically nonlinear analysis of frames
using one element / member. J. Constr. Steel Res. 20, 271–289. Luiz Fernando Martha Associate professor of the
Tang, Y.Q., Zhou, Z.H., Chan, S.L., 2015. Nonlinear beam-column element under Department of Civil Engineering at Pontifical Catholic
consistent deformation. World Scientific Publishing Company 15 (5), 123–144. University of Rio de Janeiro. B.Sc. (1977), M.Sc. (1980) in
Timoshenko, S.P., Gere, J.M., 1963. Theory of Elastic Stability. McGraw Hill, Civil Engineering from Pontifical Catholic University of
Singapura.
Rio de Janeiro (PUC-Rio) and Ph.D. (1989) in Structural
Ting, B.Y., Mockry, E.F., 1984. Beam on elastic foundation finite element. J. Struct.
Engineering from Cornell University. Postdoctoral
Eng. 110 (10), 2324–2339.
Yang, Y.B., Kuo, S.R., 1994. Theory & analysis of nonlinear framed structures, researcher at Superior Technical Institute from Techni-
Prentice Hall. Simon & Schuster (Asia) Pte ltd, Singapura. cal University of Lisbon in 2012. Acts as coordinator of
Yang, Y.B., Leu, L.J., 1994. Non-linear stiffnesses in analysis of planar frames. research projects at Tecgraf / PUC-Rio (Tecgraf Institute
Comput. Methods Appl. Mech. Eng. 117, 233–247. for Technical and Scientific Software Development). He
Yunhua, L., 1998. Explanation and elimination of shear locking and membrane has experience in structural analysis and his main
locking with field consistence approach. Comput. Methods Appl. Mech. Eng. research interests include Computer Graphics, Geo-
162, 249–269. metric Modeling, Numerical Methods Applied to Engineering and Geology Simu-
Zhaohua, F., Cook, R.D., 1983. Beam elements on two-parameter elastic foundations. lations, Computational Fracture Mechanics, and Educational Software for
J. Eng. Mech. 109 (6), 1390–1402. Engineering Teaching.
Zheng, X., Dong, S., 2011. An eigen-based high-order expansion basis for structured
spectral elements. Journal of Computational Physics 230, 8602–8753.
Zugic, L., Brcic, S., Gopcevic, S., 2016. Computer-based analysis of spatial frames
according to second order theory. Gradevinar 68 (5), 381–398.
22

1 s2.0 S0020768321000688 Main

Uploaded by

Copyright:

Available Formats

1 s2.0 S0020768321000688 Main

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

1 s2.0 S0020768321000688 Main

Uploaded by

Copyright:

Available Formats

International Journal of Solids and Structures 222–223 (2021) 111003

Contents lists available at ScienceDirect

International Journal of Solids and Structures

A unified approach to the Timoshenko 3D beam-column element

1. Introduction so-called cubic Hermitian interpolation functions (Rodrigues

2.2. Timoshenko beam theory 3

Fig. 3.1. Solution composition from the direct stiffness method.

that, linear analyses can be performed without discretization. Sim- d

hh ðxÞ ¼ c1 lcosðlxÞ c2 lsinðlxÞ þ c3 ð3:7Þ

3.2. Timoshenko beam theory

In TBT, the cross-section rotation is not the gradient of the

Fig. 5.1. Bending and torsion for a spatial element.

5.3. Local stiffness matrix equations

The development of the local stiffness matrix using cubic func-

Fig. 5.3. Beam displacement field.

@u @u0 @2w @w @u @w0 @w0 0 0 0 !2 1 1 1

fdhx gT 0 M y N 0hx N 0v dxfv g fdwgT 0 M z N 0w N 0hx dxfhx g   

fdhx gT 0 M z N0hx N 0w dxfwg

t xy zdA ¼ ða 1ÞM x ð5:57Þ T R L M x n 0 o RL n oT

When the Timoshenko beam theory is considered, the displace-

u ¼ u0 zhy yhz v ¼ v 0 zhx w ¼ w0 þ yhx ð5:60Þ

Using a trigonometric series approximation given by expression Z Z Z    

In this section, numerical examples were developed with no

EBBT_Large_Complete – Euler-Bernoulli beam theory, hyper- 6.1. Isolated columns

Fig. 6.1. Analyzed columns (Silva et al., 2016).

Fig. 6.2. Equilibrium paths for a clamped column.

Fig. 6.3. Equilibrium paths for a simply supported column.

Euler Bernoulli Beam Theory - EBBT

Fig. 6.5. Beam-column (Timoshenko and Gere, 1963).

Fig. 6.6. Equilibrium paths for a beam-column.

Discretization 1 Element (Timoshenko and Gere, 1963)

Fig. 6.9. Roorda spatial frame.

3 terms. These formulations give better approximations when

Fig. 6.8. Equilibrium paths for a spatial frame.

Discretization 1 Element 4 Elements

Fig. 6.10. Equilibrium paths for a Roorda spatial frame.

Discretization 1 Element 4 Elements

form factor of v ¼ 5=6 and Poisson’s ratio of m ¼ 0:3 and is loaded

6.7. Asymmetric spatial frame

The last example studies an asymmetric frame with moderate

stiffness matrix is corrected by employing finite rotations. The

Discretization 1 Element 4 Elements

Fig. 6.14. Deformed spatial frame with torsion.

Discretization 1 Element 4 Elements

Fig. 6.17. Deformed asymmetric spatial frame.

Fig. 6.18. Equilibrium paths for an asymmetric spatial frame.

Discretization 1 Element 4 Elements

You might also like

fdhx gT 0 M y N 0hx N 0v dxfv g fdwgT 0 M z N 0w N 0hx dxfhx g

Using a trigonometric series approximation given by expression Z Z Z