# Bars under Torsional Loading: A Generalized Beam Theory Approach

Authors: Evangelos J. Sapountzakis
Journal: ISRN Civil Engineering
 (2013)

---

## Abstract

In this paper both the static and dynamic analyses of the geometrically linear or nonlinear, elastic or elastoplastic nonuniform torsion problems of bars of constant or variable arbitrary cross section are presented together with the corresponding research efforts and the conclusions drawn from examined cases with great practical interest. In the presented analyses, the bar is subjected to arbitrarily distributed or concentrated twisting and warping moments along its length, while its edges are supported by the most general torsional boundary conditions. For the dynamic problems, a distributed mass model system is employed taking into account the warping inertia. The analysis of the aforementioned problems is complete by presenting the evaluation of the torsion and warping constants of the bar, its displacement field, its stress resultants together with the torsional shear stresses and the warping normal and shear stresses at any internal point of the bar. Moreover, the construction of the stiffness matrix and the corresponding nodal load vector of a bar of arbitrary cross section taking into account warping effects are presented for the development of a beam element for static and dynamic analyses. Having in mind the disadvantages of the 3D FEM solutions, the importance of the presented beamlike analyses becomes more evident.


---

## Body

## 1. Introduction

In engineering practice, we often come across the analysis of members of structures subjected to twisting moments. Curved bridges, ribbed plates subjected to eccentric loading, or columns laid out irregularly in the interior of a plate due to functional requirements are the most common examples.When the warping of the cross section of a member is not restrained, the applied twisting moment is undertaken from the Saint-Venant [1] shear stresses. In this case the angle of twist per unit length remains constant along the bar. However, in most cases arbitrary torsional boundary conditions are applied either at the edges or at any other interior point of the bar due to construction requirements. This bar under the action of general twisting loading is leaded to nonuniform torsion, while the angle of twist per unit length is no longer constant along it. The consequences of restrained warping were first presented by Marguerre [2].Although there is extended literature on the Saint-Venant uniform torsion problem for homogeneous isotropic cylindrical bars with simply or multiply connected cross sections [3–8], the extensive use of structural elements subjected to torsional loading necessitates a reliable, accurate, and general analysis of the torsion problem of bars of arbitrary cross section avoiding the restrictions of the Saint-Venant torsion theory.In the rest of this paper, both the static and dynamic analyses of the geometrically linear or nonlinear, elastic or elastoplastic nonuniform torsion problems of bars of constant or variable arbitrary cross section are presented together with the corresponding research efforts and the conclusions drawn from examined cases with great practical interest.


## 2. Linear Elastic Nonuniform Torsion of Bars

In the last decades considerable work has been done on the elastic linear problem of nonuniform torsion of bars. Especially because of the mathematical complexity of the problem, the existing analytical solutions are limited to symmetric cross sections of simple geometry, loading, and boundary conditions [9–16]. Moreover, numerical methods such as the finite element method [17] or the boundary element method [18–20] have also been used for the analysis of the nonuniform torsional problem, in the case the geometry of the cross section, the boundary conditions, or the loading are not simple. In all the aforementioned references the analysis of the nonuniform torsional problem is not complete, since the secondary shear stresses due to warping are not evaluated. Finally, Sapountzakis and Mokos in [21, 22] developed a boundary element solution for the general linear elastic nonuniform torsion problem of homogeneous or composite prismatic bars of arbitrary cross section subjected to an arbitrarily distributed or concentrated twisting moment and supported by the most general linear torsional boundary conditions. In these latter research efforts the evaluation of the secondary warping function is accomplished leading to the computation of the secondary shear stresses due to warping, while the developed procedure is a pure BEM [23], since it requires only boundary discretization.In order to formulate the aforementioned problem, let us consider a prismatic bar of lengthl (Figure 1), of constant arbitrary cross-section of area A. The homogeneous isotropic and linearly elastic material of the bar’s cross-section, with modulus of elasticity E, shear modulus G, and Poisson’s ratio ν, occupies the two dimensional multiply connected region Ω of the y, z plane and is bounded by the Γj(j=1,2,…,K) boundary curves, which are piecewise smooth; that is, they may have a finite number of corners. In Figure 1(b)CYZ is the principal coordinate system through the cross section’s centroid C, while yC,zC are its coordinates with respect to Syz system of axes through the cross section’s shear center S. The bar is subjected to the arbitrarily distributed or concentrated conservative twisting mt=mt(x) moment acting in the x direction (Figure 1(a)). This torsional loading results in a rotation with respect to the shear center S with angle of twist θx(x). The shear center S coincides with the center of twist of the cross section provided that no other axis or rotation is imposed due to construction requirements.Prismatic bar subjected to a twisting moment (a) with a cross section of arbitrary shape occupying the two dimensional regionΩ (b).
(a)
(b)Under the aforementioned loading the displacement field of the bar with respect to theSyz system of axes is given as(1a)u(x,y,z)=θx′(x)ϕS(y,z),(1b)v(x,y,z)=-zθx(x),(1c)w(x,y,z)=yθx(x),where the displacement component u constitutes the warping of the cross section, θx′(x)=dθx/dx denotes the rate of change of the angle of twist θx regarded as the torsional curvature, and ϕS(y,z) is the warping function with respect to the shear center S, which is characterized as the warping function. The index S indicates that the warping refers to the rotation axis S.Substituting the aforementioned displacement field to the linearized strain-displacement relations (infinitesimal strain tensor) of the three-dimensional elasticity, the nonvanishing strain components are written as(2a)εxx=∂u∂x=θx′′(x)·ϕS(y,z),(2b)γxy=∂u∂y+∂v∂x=θx′(x)·(∂ϕS∂y-z),(2c)γxz=∂u∂z+∂w∂x=θx′(x)·(∂ϕS∂z+y)while the resulting nonzero stress components (Cauchy stress tensor) in the region Ω employing the stress-strain relations (constitutive relations) and assuming an isotropic and homogeneous material for zero Poisson ratio are derived as(3a)σxx=E(1+ν)(1-2ν)[(1-ν)εxx+ν(εyy+εzz)]=E·θx′′(x)·ϕS(y,z),(3b)τxy=G·γxy=G·θx′(x)·(∂ϕS∂y-z),(3c)τxz=G·γxz=G·θx′(x)·(∂ϕS∂z+y).


### 2.1. Equations of Local Equilibrium

Substituting the stress tensor components in the three equilibrium equations of the three-dimensional elasticity theory ignoring body forces yields(4a)∂∂y[G·θx′(x)·(∂ϕS∂y-z)]+∂∂z[G·θx′(x)(∂ϕS∂z+y)]+∂∂x[E·θx′′(x)·ϕS]=0,(4b)G·θx′′(x)·(∂ϕS∂z+y)=0,(4c)G·θx′′(x)·(∂ϕS∂y-z)=0.The first two of these equations, as it can be observed are not satisfied and this constitutes an inconsistency of the nonuniform torsion theory of bars. Only in the case that θx′′(x)=0, that is, the rate of change of θx′(x) is constant (case of uniform Saint Venant torsion), and the first two equilibrium equations are satisfied. From (4c) it follows that
(5)∂2ϕS∂y2+∂2ϕS∂z2=-E·θx′′′(x)G·θx′(x)·ϕS.
Having in mind that the first θx′ and third θx′′′ derivatives of the angle of twist are functions of the longitudinal coordinate x, it can be observed that the solution of the partial differential equation (5) in the general case should be a function of coordinate x as well. This fact contradicts relation (1a), where the warping function ϕS(y,z) is defined as independent of the coordinate x. Thus, the left hand side of (5) is a function only of variables y, z, since ϕS=ϕS(y,z), while at the same time the right hand side is a function of variable x. To remove this inconsistency it is assumed that the arising shear stresses are decomposed to a primaryτPpart and a secondaryτS one developed due to warping. This decomposition of shear stress is justified by the consideration of equilibrium of normal stresses due to warping in an infinitesimal element of the cross section as shown in Figure 2. Thus, in a bar of arbitrary cross section subjected to nonuniform torsion, normal stresses σxxw arise, depending on the developed warping and vary along the longitudinal axis of the bar (Figure 2(a)). Considering the intersection (A) in Figure 2(b) it can be observed that the infinitesimal change of normal stresses dσxxw due to distortion can be equilibrated only by shear stresses along the intersection of (A) and which employing Cauchy’s theorem lead to the development of secondary shear stresses τS (Figure 2(c)) on the plane of the cross section. Thus, according to the previously mentioned primary shear stresses, which express the shear stresses of uniform torsion (Saint-Venant) with the difference that θx′(x) is not constant and to the normal stresses due to warping resulting from the deformation (primary ones), secondary shear stresses result so as to equilibrate the aforementioned normal stresses σxxw.Figure 2
Normal and shear stresses due to warping.Based on the above decomposition of shear stress in primary and secondary components, that is(6a)τxy=τxyP+τxyS,(6b)τxz=τxzP+τxzSthe primary and secondary components of these stresses and the normal stresses due to warping are defined according to the relations as follows:(7a)τxyP=G·θx′(x)·(∂ϕSP∂y-z),(7b)τxyS=G·∂ϕSS∂y,(8a)
τ
x
z
P
=
G
·
θ
x
′
(
x
)
·
(
∂
ϕ
S
P
∂
z
+
y
)
,
(8b)
τ
x
z
S
=
G
·
∂
ϕ
S
S
∂
z
,
(9)
σ
x
x
w
=
E
·
θ
x
′′
(
x
)
·
ϕ
S
P
(
y
,
z
)
,
where the functions ϕSP(y,z) and ϕSS(x,y,z) are called primary and secondary warping function, respectively. Substituting relations (6a)–(9) in the first equilibrium equation of the theory of elasticity ignoring body forces it follows that
(10)∂σxxw∂x+∂τxyP∂y+∂τxzP∂z+∂τxyS∂y+∂τxzS∂z=0.
Equation (10) indicates that the stress state of the bar subjected to torsional loading results from the superposition of primary τP, secondary τS shear stress and normal stresses due to warping σxxw. In order to satisfy (10) it is required that both the terms originating from the primary shear stresses as well as these from the normal and the secondary shear stresses due to warping to vanish, that is
(11)∂τxyP∂y+∂τxzP∂z=0,∂τxyS∂y+∂τxzS∂z+∂σxxw∂x=0.
Substituting (7a)–(9) into (11) the following relations are obtained:
(12)∇2ϕSP=∂2ϕSP∂y2+∂2ϕSP∂z2=0,∇2ϕSS=∂2ϕSS∂y2+∂2ϕSS∂z2=-E·θx′′′(x)G·ϕSP.To formulate theboundary conditions of the primary and secondary warping functions the shear stress components are observed on the boundary of the cross section. By examining in Figure 3 the infinitesimal surface dydz, the following relations are obtained:
(13)τxn=τxy·ny+τxz·nz,τxt=-τxy·nz+τxz·ny,
where ny=cos(y,n)=cosα=dy/dn=dz/ds and nz=sin(z,n)=sinα=dz/dn=-dy/ds are the direction cosines of the vector n_ normal to the boundary of the cross section. Substituting relations (6a) and (6b) that express the decomposition of shear stresses in primary and secondary ones into relations (13) yields(14a)τxnP=τxyP·ny+τxzP·nz,(14b)τxtP=-τxyP·nz+τxzP·ny,Figure 3
Shear stress components on the boundary of the cross section.(15a)
τ
x
n
S
=
τ
x
y
S
·
n
y
+
τ
x
z
S
·
n
z
,
(15b)
τ
x
t
S
=
-
τ
x
y
S
·
n
z
+
τ
x
z
S
·
n
y
.
Substituting the expressions of primary and secondary shear stresses given by (7a), (7b), (8a), and (8b) in (14a), (14b), (15a), and (15b) the following relations are obtained:(16a)τxnP=G·θx′(x)(∂ϕSP∂n-z·ny+y·nz),(16b)τxtP=G·θx′(x)(∂ϕSP∂t+y·ny+z·nz),(16c)τxnS=G·∂ϕSS∂n,(16d)τxtS=G·∂ϕSS∂t.Relations (16a), (16c) and (16b), (16d) give the normal and tangential to the boundary of the cross section primary and secondary shear stresses, respectively. Since the lateral surface of the bar is unloaded along the longitudinal direction, the normal shear stresses at the boundary of the cross section should vanish in order to satisfy the equilibrium requirements. To fullfil this boundary condition it is required that both the primary and the secondary shear stress components to vanish, that is(17a)τxnP=0,(17b)τxnS=0.Substituting (16a) and (16c) into (17a) and (17b) and taking into account the fact that in general Gθx′≠0, the boundary conditions of the primary and the secondary warping functions are obtained as(18a)∂ϕSP∂n=z·ny-y·nz,(18b)∂ϕSS∂n=0.Summarizing the above mentioned it is concluded that for the determination of the primaryϕSP(y,z) and the secondary ϕSS(x,y,z) warping functions the solution of the following boundary value problems is required.For the Primary
ϕ
S
P
(
y
,
z
)
Warping Function. Neumann problem of Laplace differential equation is as follows:(19a)∇2ϕSP=∂2ϕSP∂y2+∂2ϕSP∂z2=0inthedomainofcrosssectionΩ,(19b)∂ϕSP∂n=z·ny-y·nzontheboundaryofthecrosssectionΓ.For the SecondaryϕSS(x,y,z)Warping Function.Neumann problem of the Poisson differential equation is as follows:(20a)∇2ϕSS=∂2ϕSS∂y2+∂2ϕSS∂z2=-E·θx′′′(x)G·ϕSPinthedomainofcrosssectionΩ,(20b)∂ϕSS∂n=0ontheboundaryofthecrosssectionΓ.Finally, based on the above mentioned, the displacement field of the bar ((1a), (1b), and (1c)) is formed as(21a)u(x,y,z)=θx′(x)·ϕSP(y,z)+ϕSS(x,y,z),(21b)v(x,y,z)=-zθx(x),(21c)w(x,y,z)=yθx(x).It is worth here noting that in the case the originO of the coordinate system is a point of the y¯, z¯ plane other than the shear center (Figure 1(b)), the warping function with respect to this point φOP is first established from the Neumann problem (19a) and (19b) substituting φSP by φOP. Using the evaluated warping function φOP, φSP is then established using the transformation given by the following equation [2]:
(22)ϕSP(y,z)=ϕOP(y-,z-)-y-z-S+z-y-S+c-,
where y=y--y-S, z=z--z-S, y-S, z-S are the coordinates of the shear center S with respect to the arbitrary coordinate system Oy-z- (Figure 1(b)) and c- is an integration constant. The latter is given from the solution of the following linear system of equations:(23a)Sy-y-S-Sz-z-S+Ac-=-R-SP,(23b)Iyy-y-S+Iyz-z-S+Sy-c-=-Ry-P,(23c)Iyz-y-S+Izz-z-S-Sz-c-=Rz-P,where(24a)A=∫ΩdΩ,(24b)Sy-=∫Ωz-dΩ,(24c)Sz-=∫Ωy-dΩ,(24d)Iyy-=∫Ωz-2dΩ,(24e)Izz-=∫Ωy-2dΩ,(24f)Iyz-=-∫Ωy-z-dΩ,is the cross section area, the static moments of inertia relative to the axes y-, z-, the bending moments of inertia relative to the axes y-, z-, and the product of inertia, respectively, while in relations (23a), (23b), and (23c) the warping moments R-SP,Ry-P,Rz-P, are defined as(25a)R-SP=∫ΩϕOPdΩ,(25b)Ry-P=∫Ωz-ϕOPdΩ,(25c)Rz-P=∫Ωy-ϕOPdΩ.Moreover, the evaluated warping function ϕSS from the solution of the Neumann problem (20a) and (20b) contains an integration constant cS (parallel displacement of the cross section along the beam axis), which can be obtained from [24] as follows:
(26)cS=1A∫Ωϕ~SSdΩ
and the main secondary warping function ϕ~SS is given as
(27)ϕ~SS=ϕSS+cS.


### 2.2. Equations of Global Equilibrium

The already established shear stressesτxy and τxz yield components of torque, which arise through integration over the cross section. Thus, the resulting twisting moment is obtained as
(28)Mt=∫Ω(τxz·y-τxy·z)dΩ.
Introducing the approximation of decomposition of shear stress in primary and secondary components, the twisting moment of the cross section is divided into a primary component MtP originating from the primary shear stresses τP due to twisting (as in uniform torsion) and a secondary component MtS, originating from the secondary shear stresses, that is the restraint of warping (Figure 4). Thus, according to the nonuniform torsion theory in an arbitrary cross section of the bar the following relation is valid:
(29)Mt=MtP+MtS,
where after employing (6a) and (6b) and some algebra the aforementioned twisting moment components are given as(30a)MtP=∫Ω[τxyP(∂ϕSP∂y-z)+τxzP·(∂ϕSP∂z+y)]dΩ,(30b)MtS=∫Ω(-τxyS∂ϕSP∂y-τxzS∂ϕSP∂z)dΩ.Substituting (7a)–(8b) in (30a) and (30b), employing the Gauss-Green theorem, taking into account (17a), (17b), (19a), and (20a) and after some algebra the following relations for the twisting moment components are obtained:(31a)MtP=-GItθx′(x),(31b)MtS=-ECSθx′′′(x),where(32a)It=∫Ω(y2+z2+y·∂ϕSP∂z-z·∂ϕSP∂y)dΩ,(32b)CS=∫Ω(ϕSP)2dΩ.In (32a) and (32b) the quantity It denotes the torsional constant introduced by Saint-Venant, while the quantity CS denotes the warping constant. The quantity GIt denotes the torsional rigidity, while the quantity ECS denotes the warping rigidity of the cross section. It is noted that the torsional constant is independent of the position of the coordinate system, while the warping constant refers to the center of twist S.Distribution of primary (a) and secondary (b) shear stresses and resulting torsional moments.
(a)
(b)To derive the equilibrium equation for the nonuniform torsion problem of a homogeneous isotropic bar, the equilibrium of an infinitesimal partdx of the bar against twisting moments is examined and observing Figure 5 it follows that
(33)Mt+∂Mt∂xdx-Mt+mtdx=0
or after some algebra
(34)∂Mt∂x=-mt.
Substituting relations (29) and (31a) and (31b) into (34), the following fourth-order differential equation of equilibrium of a homogeneous isotropic bar subjected to nonuniform torsion is obtained:
(35)ECSd4θxdx4-GItd2θxdx2=mt
which is independent of the torsional boundary conditions.Figure 5
Stress resultants in an infinitesimal segment of an elastic bar subjected to an arbitrary concentrated or distributed torsional loadingmt=mt(x).In analogy with the bending momentsMy, Mz a new stress resultant is defined, called warping moment (or bimoment) and given by
(36)Mw=-∫ΩϕSPσxxwdΩ.
The need for the definition of this new stress resultant stems from the fact that, whereas My=Mz=N=0, there still exist normal stresses acting on the cross section; hence if a new quantity is not considered, the elastic energy due to σxx stresses will be ignored. Substituting the expression (9) of the stress component σxxw into (36), the latter is written as
(37)Mw=-Eθx′′∫Ω(ϕSP)2dΩ
or through (32b)
(38)Mw=-ECSθx′′
and therefore relation (9) can be written in the form of
(39)σxxw=-MwCSφSP.
Equation (39) shows that the quantity of warping moment encompasses the generic characteristics of a stress resultant of theory of elasticity. Schardt [25] refers to it as a “higher order stress resultant.” Combining relations (31b) and (38), the relation which correlates the secondary twisting moment and the warping moment arises as
(40)dMwdx=MtS.
Thus, the problem of nonuniform torsion of a homogeneous isotropic bar is reduced to solving the fourth-order differential equation with respect to the angle of twist θx(x) of the cross section, given by (35). The solution of this equation depends on both the torsional loading of the bar and the torsional support conditions at the ends or inside the bar. The most general linear torsional boundary conditions at the ends of the bar are described by the relations(41a)a1θx+a2Mt=a3,(41b)β1θx′+β2Mw=β3.It is worth noting that all types of conventional boundary conditions (e.g., fixed, forked support, free end, and elastic support) arise from relations (41a) and (41b) after defining appropriately functions ai, βi(i=1,2,3). For example, in the case of a torsionally fixed support the above functions take the values a1=β1=1,  a2=a3=β2=β3=0.From the examined example problems presented in [18–22] it is concluded that (i)
the magnitude of the evaluated normal stresses, due to restrained warping compared with those due to bending, necessitates the consideration of these additional normal stresses near the restrained edges;(ii)
the magnitude of the evaluated warping shear stresses due to restrained warping is remarkable and for an “accurate” analysis these additional shear stresses should not be ignored, especially near the restrained edges.


## 2.1. Equations of Local Equilibrium

Substituting the stress tensor components in the three equilibrium equations of the three-dimensional elasticity theory ignoring body forces yields(4a)∂∂y[G·θx′(x)·(∂ϕS∂y-z)]+∂∂z[G·θx′(x)(∂ϕS∂z+y)]+∂∂x[E·θx′′(x)·ϕS]=0,(4b)G·θx′′(x)·(∂ϕS∂z+y)=0,(4c)G·θx′′(x)·(∂ϕS∂y-z)=0.The first two of these equations, as it can be observed are not satisfied and this constitutes an inconsistency of the nonuniform torsion theory of bars. Only in the case that θx′′(x)=0, that is, the rate of change of θx′(x) is constant (case of uniform Saint Venant torsion), and the first two equilibrium equations are satisfied. From (4c) it follows that
(5)∂2ϕS∂y2+∂2ϕS∂z2=-E·θx′′′(x)G·θx′(x)·ϕS.
Having in mind that the first θx′ and third θx′′′ derivatives of the angle of twist are functions of the longitudinal coordinate x, it can be observed that the solution of the partial differential equation (5) in the general case should be a function of coordinate x as well. This fact contradicts relation (1a), where the warping function ϕS(y,z) is defined as independent of the coordinate x. Thus, the left hand side of (5) is a function only of variables y, z, since ϕS=ϕS(y,z), while at the same time the right hand side is a function of variable x. To remove this inconsistency it is assumed that the arising shear stresses are decomposed to a primaryτPpart and a secondaryτS one developed due to warping. This decomposition of shear stress is justified by the consideration of equilibrium of normal stresses due to warping in an infinitesimal element of the cross section as shown in Figure 2. Thus, in a bar of arbitrary cross section subjected to nonuniform torsion, normal stresses σxxw arise, depending on the developed warping and vary along the longitudinal axis of the bar (Figure 2(a)). Considering the intersection (A) in Figure 2(b) it can be observed that the infinitesimal change of normal stresses dσxxw due to distortion can be equilibrated only by shear stresses along the intersection of (A) and which employing Cauchy’s theorem lead to the development of secondary shear stresses τS (Figure 2(c)) on the plane of the cross section. Thus, according to the previously mentioned primary shear stresses, which express the shear stresses of uniform torsion (Saint-Venant) with the difference that θx′(x) is not constant and to the normal stresses due to warping resulting from the deformation (primary ones), secondary shear stresses result so as to equilibrate the aforementioned normal stresses σxxw.Figure 2
Normal and shear stresses due to warping.Based on the above decomposition of shear stress in primary and secondary components, that is(6a)τxy=τxyP+τxyS,(6b)τxz=τxzP+τxzSthe primary and secondary components of these stresses and the normal stresses due to warping are defined according to the relations as follows:(7a)τxyP=G·θx′(x)·(∂ϕSP∂y-z),(7b)τxyS=G·∂ϕSS∂y,(8a)
τ
x
z
P
=
G
·
θ
x
′
(
x
)
·
(
∂
ϕ
S
P
∂
z
+
y
)
,
(8b)
τ
x
z
S
=
G
·
∂
ϕ
S
S
∂
z
,
(9)
σ
x
x
w
=
E
·
θ
x
′′
(
x
)
·
ϕ
S
P
(
y
,
z
)
,
where the functions ϕSP(y,z) and ϕSS(x,y,z) are called primary and secondary warping function, respectively. Substituting relations (6a)–(9) in the first equilibrium equation of the theory of elasticity ignoring body forces it follows that
(10)∂σxxw∂x+∂τxyP∂y+∂τxzP∂z+∂τxyS∂y+∂τxzS∂z=0.
Equation (10) indicates that the stress state of the bar subjected to torsional loading results from the superposition of primary τP, secondary τS shear stress and normal stresses due to warping σxxw. In order to satisfy (10) it is required that both the terms originating from the primary shear stresses as well as these from the normal and the secondary shear stresses due to warping to vanish, that is
(11)∂τxyP∂y+∂τxzP∂z=0,∂τxyS∂y+∂τxzS∂z+∂σxxw∂x=0.
Substituting (7a)–(9) into (11) the following relations are obtained:
(12)∇2ϕSP=∂2ϕSP∂y2+∂2ϕSP∂z2=0,∇2ϕSS=∂2ϕSS∂y2+∂2ϕSS∂z2=-E·θx′′′(x)G·ϕSP.To formulate theboundary conditions of the primary and secondary warping functions the shear stress components are observed on the boundary of the cross section. By examining in Figure 3 the infinitesimal surface dydz, the following relations are obtained:
(13)τxn=τxy·ny+τxz·nz,τxt=-τxy·nz+τxz·ny,
where ny=cos(y,n)=cosα=dy/dn=dz/ds and nz=sin(z,n)=sinα=dz/dn=-dy/ds are the direction cosines of the vector n_ normal to the boundary of the cross section. Substituting relations (6a) and (6b) that express the decomposition of shear stresses in primary and secondary ones into relations (13) yields(14a)τxnP=τxyP·ny+τxzP·nz,(14b)τxtP=-τxyP·nz+τxzP·ny,Figure 3
Shear stress components on the boundary of the cross section.(15a)
τ
x
n
S
=
τ
x
y
S
·
n
y
+
τ
x
z
S
·
n
z
,
(15b)
τ
x
t
S
=
-
τ
x
y
S
·
n
z
+
τ
x
z
S
·
n
y
.
Substituting the expressions of primary and secondary shear stresses given by (7a), (7b), (8a), and (8b) in (14a), (14b), (15a), and (15b) the following relations are obtained:(16a)τxnP=G·θx′(x)(∂ϕSP∂n-z·ny+y·nz),(16b)τxtP=G·θx′(x)(∂ϕSP∂t+y·ny+z·nz),(16c)τxnS=G·∂ϕSS∂n,(16d)τxtS=G·∂ϕSS∂t.Relations (16a), (16c) and (16b), (16d) give the normal and tangential to the boundary of the cross section primary and secondary shear stresses, respectively. Since the lateral surface of the bar is unloaded along the longitudinal direction, the normal shear stresses at the boundary of the cross section should vanish in order to satisfy the equilibrium requirements. To fullfil this boundary condition it is required that both the primary and the secondary shear stress components to vanish, that is(17a)τxnP=0,(17b)τxnS=0.Substituting (16a) and (16c) into (17a) and (17b) and taking into account the fact that in general Gθx′≠0, the boundary conditions of the primary and the secondary warping functions are obtained as(18a)∂ϕSP∂n=z·ny-y·nz,(18b)∂ϕSS∂n=0.Summarizing the above mentioned it is concluded that for the determination of the primaryϕSP(y,z) and the secondary ϕSS(x,y,z) warping functions the solution of the following boundary value problems is required.For the Primary
ϕ
S
P
(
y
,
z
)
Warping Function. Neumann problem of Laplace differential equation is as follows:(19a)∇2ϕSP=∂2ϕSP∂y2+∂2ϕSP∂z2=0inthedomainofcrosssectionΩ,(19b)∂ϕSP∂n=z·ny-y·nzontheboundaryofthecrosssectionΓ.For the SecondaryϕSS(x,y,z)Warping Function.Neumann problem of the Poisson differential equation is as follows:(20a)∇2ϕSS=∂2ϕSS∂y2+∂2ϕSS∂z2=-E·θx′′′(x)G·ϕSPinthedomainofcrosssectionΩ,(20b)∂ϕSS∂n=0ontheboundaryofthecrosssectionΓ.Finally, based on the above mentioned, the displacement field of the bar ((1a), (1b), and (1c)) is formed as(21a)u(x,y,z)=θx′(x)·ϕSP(y,z)+ϕSS(x,y,z),(21b)v(x,y,z)=-zθx(x),(21c)w(x,y,z)=yθx(x).It is worth here noting that in the case the originO of the coordinate system is a point of the y¯, z¯ plane other than the shear center (Figure 1(b)), the warping function with respect to this point φOP is first established from the Neumann problem (19a) and (19b) substituting φSP by φOP. Using the evaluated warping function φOP, φSP is then established using the transformation given by the following equation [2]:
(22)ϕSP(y,z)=ϕOP(y-,z-)-y-z-S+z-y-S+c-,
where y=y--y-S, z=z--z-S, y-S, z-S are the coordinates of the shear center S with respect to the arbitrary coordinate system Oy-z- (Figure 1(b)) and c- is an integration constant. The latter is given from the solution of the following linear system of equations:(23a)Sy-y-S-Sz-z-S+Ac-=-R-SP,(23b)Iyy-y-S+Iyz-z-S+Sy-c-=-Ry-P,(23c)Iyz-y-S+Izz-z-S-Sz-c-=Rz-P,where(24a)A=∫ΩdΩ,(24b)Sy-=∫Ωz-dΩ,(24c)Sz-=∫Ωy-dΩ,(24d)Iyy-=∫Ωz-2dΩ,(24e)Izz-=∫Ωy-2dΩ,(24f)Iyz-=-∫Ωy-z-dΩ,is the cross section area, the static moments of inertia relative to the axes y-, z-, the bending moments of inertia relative to the axes y-, z-, and the product of inertia, respectively, while in relations (23a), (23b), and (23c) the warping moments R-SP,Ry-P,Rz-P, are defined as(25a)R-SP=∫ΩϕOPdΩ,(25b)Ry-P=∫Ωz-ϕOPdΩ,(25c)Rz-P=∫Ωy-ϕOPdΩ.Moreover, the evaluated warping function ϕSS from the solution of the Neumann problem (20a) and (20b) contains an integration constant cS (parallel displacement of the cross section along the beam axis), which can be obtained from [24] as follows:
(26)cS=1A∫Ωϕ~SSdΩ
and the main secondary warping function ϕ~SS is given as
(27)ϕ~SS=ϕSS+cS.


## 2.2. Equations of Global Equilibrium

The already established shear stressesτxy and τxz yield components of torque, which arise through integration over the cross section. Thus, the resulting twisting moment is obtained as
(28)Mt=∫Ω(τxz·y-τxy·z)dΩ.
Introducing the approximation of decomposition of shear stress in primary and secondary components, the twisting moment of the cross section is divided into a primary component MtP originating from the primary shear stresses τP due to twisting (as in uniform torsion) and a secondary component MtS, originating from the secondary shear stresses, that is the restraint of warping (Figure 4). Thus, according to the nonuniform torsion theory in an arbitrary cross section of the bar the following relation is valid:
(29)Mt=MtP+MtS,
where after employing (6a) and (6b) and some algebra the aforementioned twisting moment components are given as(30a)MtP=∫Ω[τxyP(∂ϕSP∂y-z)+τxzP·(∂ϕSP∂z+y)]dΩ,(30b)MtS=∫Ω(-τxyS∂ϕSP∂y-τxzS∂ϕSP∂z)dΩ.Substituting (7a)–(8b) in (30a) and (30b), employing the Gauss-Green theorem, taking into account (17a), (17b), (19a), and (20a) and after some algebra the following relations for the twisting moment components are obtained:(31a)MtP=-GItθx′(x),(31b)MtS=-ECSθx′′′(x),where(32a)It=∫Ω(y2+z2+y·∂ϕSP∂z-z·∂ϕSP∂y)dΩ,(32b)CS=∫Ω(ϕSP)2dΩ.In (32a) and (32b) the quantity It denotes the torsional constant introduced by Saint-Venant, while the quantity CS denotes the warping constant. The quantity GIt denotes the torsional rigidity, while the quantity ECS denotes the warping rigidity of the cross section. It is noted that the torsional constant is independent of the position of the coordinate system, while the warping constant refers to the center of twist S.Distribution of primary (a) and secondary (b) shear stresses and resulting torsional moments.
(a)
(b)To derive the equilibrium equation for the nonuniform torsion problem of a homogeneous isotropic bar, the equilibrium of an infinitesimal partdx of the bar against twisting moments is examined and observing Figure 5 it follows that
(33)Mt+∂Mt∂xdx-Mt+mtdx=0
or after some algebra
(34)∂Mt∂x=-mt.
Substituting relations (29) and (31a) and (31b) into (34), the following fourth-order differential equation of equilibrium of a homogeneous isotropic bar subjected to nonuniform torsion is obtained:
(35)ECSd4θxdx4-GItd2θxdx2=mt
which is independent of the torsional boundary conditions.Figure 5
Stress resultants in an infinitesimal segment of an elastic bar subjected to an arbitrary concentrated or distributed torsional loadingmt=mt(x).In analogy with the bending momentsMy, Mz a new stress resultant is defined, called warping moment (or bimoment) and given by
(36)Mw=-∫ΩϕSPσxxwdΩ.
The need for the definition of this new stress resultant stems from the fact that, whereas My=Mz=N=0, there still exist normal stresses acting on the cross section; hence if a new quantity is not considered, the elastic energy due to σxx stresses will be ignored. Substituting the expression (9) of the stress component σxxw into (36), the latter is written as
(37)Mw=-Eθx′′∫Ω(ϕSP)2dΩ
or through (32b)
(38)Mw=-ECSθx′′
and therefore relation (9) can be written in the form of
(39)σxxw=-MwCSφSP.
Equation (39) shows that the quantity of warping moment encompasses the generic characteristics of a stress resultant of theory of elasticity. Schardt [25] refers to it as a “higher order stress resultant.” Combining relations (31b) and (38), the relation which correlates the secondary twisting moment and the warping moment arises as
(40)dMwdx=MtS.
Thus, the problem of nonuniform torsion of a homogeneous isotropic bar is reduced to solving the fourth-order differential equation with respect to the angle of twist θx(x) of the cross section, given by (35). The solution of this equation depends on both the torsional loading of the bar and the torsional support conditions at the ends or inside the bar. The most general linear torsional boundary conditions at the ends of the bar are described by the relations(41a)a1θx+a2Mt=a3,(41b)β1θx′+β2Mw=β3.It is worth noting that all types of conventional boundary conditions (e.g., fixed, forked support, free end, and elastic support) arise from relations (41a) and (41b) after defining appropriately functions ai, βi(i=1,2,3). For example, in the case of a torsionally fixed support the above functions take the values a1=β1=1,  a2=a3=β2=β3=0.From the examined example problems presented in [18–22] it is concluded that (i)
the magnitude of the evaluated normal stresses, due to restrained warping compared with those due to bending, necessitates the consideration of these additional normal stresses near the restrained edges;(ii)
the magnitude of the evaluated warping shear stresses due to restrained warping is remarkable and for an “accurate” analysis these additional shear stresses should not be ignored, especially near the restrained edges.


## 3. Linear Elastic Nonuniform Torsion of Bars of Variable Cross Section

Long span box shaped bridges or concrete slab and beam structures of variable height are the most common examples of structures including members of variable cross section subjected to twisting moments. The extensive use of the aforementioned structural elements necessitates a rigorous analysis. When a bar of variable cross section is subjected to general twisting loading, this member due to this variation and/or due to the arbitrary torsional boundary conditions applied either at the edges or at any other interior point is leaded to nonuniform torsion and its angle of twist per unit length is not constant along its axis.Several researchers have dealt with beams of variable cross section ignoring the warping effects resulting from the corresponding restraints at the ends of the member [26, 27]. If the aforementioned structures are analyzed or designed for torsion considering only the effect of Saint Venant torsion resistance, the analysis may underestimate the torsion in the members and the design may be unconservative. On the contrary, to the authors’ knowledge relatively little work has been done on the problem of nonuniform torsion of bars of variable cross section with pioneer the work of Cywinski [28] adopting the finite difference method. Wekezer [29] after dividing the bar into segments along its longitudinal axis approximated their shell midsurface by arbitrary triangular shell elements and employed the finite element method to the linear membrane shell theory. This approximation generates inaccuracies, as the warping of the walls of the cross section cannot be taken into account. Moreover, Eisenberger [30] employed FEM upon polynomial approximation of the torsional and warping rigidities using “exact” shape functions to derive the exact stiffness coefficients. This application of shape functions leads also to inaccuracies in stress analysis of beams of variable cross section, as static and kinematic values at nodes and in the element region are computed only approximately and the element may not satisfy local and global equilibrium conditions [31]. In all of the aforementioned procedures the torsion and warping constants have been approximated adopting the thin tube theory. Finally, Sapountzakis and Mokos in [32, 33] developed a boundary element solution for the general linear elastic nonuniform torsion problem of homogeneous or composite bars of arbitrary variable cross section subjected to an arbitrarily distributed or concentrated twisting moment and supported by the most general linear torsional boundary conditions. In these latter research efforts, three boundary value problems with respect to the variable along the beam angle of twist and to the primary and secondary warping functions are formulated and solved employing a pure BEM [23] approach; that is, only boundary discretization is used. Both the variable warping and torsion constants together with the torsional primary shear stresses and the warping normal and secondary shear stresses are computed.In order to formulate the aforementioned problem, let us consider a bar of lengthl (Figure 6), of an arbitrarily shaped variable along its axis cross section. The homogeneous isotropic and linearly elastic material of the bar’s cross-section, with modulus of elasticity E, shear modulus G, and Poisson’s ratio ν, occupies the two dimensional multiply connected region Ω(x) of the y, z plane and is bounded by the Γj(j=1,2,…,K) boundary curves, which are piecewise smooth; that is, they may have a finite number of corners. In Figure 6(b)CYZ is the principal coordinate system through the cross section’s centroid C, while yC,zC are its coordinates with respect to Syz system of axes through the cross section’s shear center S. The bar is subjected to the arbitrarily distributed or concentrated conservative twisting mt=mt(x) moment acting in the x direction (Figure 6(a)).Bar subjected to a twisting moment (a) with a variable cross section of arbitrary shape occupying the two dimensional regionΩ(x) (b).
(a)
(b)Adopting the same displacement field (21a), (21b), and (21c) with that of the constant cross section case and following the same procedure presented for the constant cross section bar in Section 2, the following boundary value problem for the angle of twist θx=θx(x) is derived:
(42)ECSd4θxdx4+2EdCSdxd3θxdx3+(Ed2CSdx2-GIt)d2θxdx2-GdItdxdθxdx=mtinsidethebeam,(43a)α1θx+α2Mt=α3atthebeamendsx=0,l,(43b)β1dθxdx+β2Mw=β3atthebeamendsx=0,l,where the total twisting moment of the cross section is once again divided into a primary and a secondary component as this is stated in (29), where these components are given from(44a)MtP=GItdθxdx,(44b)MtS=-Eddx(CSd2θxdx2)with It the torsional and CS the warping constants of the cross section varying along the length of the bar and given from (32a) and (32b), while the resulting total twisting moment of the cross section is given as
(45)Mt=GItdθxdx-EdCSdxd2θxdx2-ECSd3θxdx3.
The warping moment Mw arising from the torsional curvature, similarly with the constant cross section case, is given from (38) and also functions ai, βi(i=1,2,3) are specified at the boundary of the beam forming the most general linear torsional boundary conditions for the beam problem including also the elastic support. Finally, as for the constant cross section case the determination of the primary ϕSP(y,z) and the secondary ϕSS(x,y,z) warping functions is achieved from the solution of the boundary value problems given from (19a), (19b), (20a), and (20b), respectively, while the primary, the secondary shear stress components, and the normal stresses due to warping are defined according to the relations (7a)–(9).From the examined example problems presented in [32, 33] it is concluded that (i)
the variation of the beam height, as expected, results in increment of the nonuniform beam behavior in torsional loading;(ii)
the magnitude of the evaluated normal and warping shear stresses due to restrained warping is remarkable and necessitates the consideration of these additional stresses near the restrained edges, especially for beams with cross section of low torsional rigidity.(iii)
the inaccuracy of the thin tube theory in calculating torsional and warping rigidities even for thin walled sections is remarkable;(iv)
analyzing a thin-walled cross section in torsional loading by modeling it with shell elements cannot give accurate results, as the warping of its walls cannot be taken into account.


## 4. 3D Beam Element of Constant or Variable Cross Section Including Warping Effect

As it has been already mentioned, the analysis of rectilinear or curved members of structures of arbitrary constant or variable cross section subjected to twisting moments is often encountered in engineering practice. Nevertheless, accurate analysis of these members is difficult to achieve for two reasons.According to the first reason, generally commercial programs consider six degrees of freedom at each node of a member of a space frame, ignoring in this way the warping effects due to the corresponding restraint at the ends of the member or due to variable twisting moment along the bar [26, 27, 34]. If the aforementioned structures are analyzed or designed for torsion considering only the effect of Saint Venant torsion resistance, the analysis may underestimate the torsion in the members and the design may be unconservative. Several researchers tried to overcome this inaccuracy only for constant cross section elements by developing a 14×14 member stiffness matrix including warping degrees of freedom at the ends of a member with open thin-walled section and assuming simple [35–38] or more complicated torsional boundary conditions [39, 40].According to the second reason, when the variable cross section structures are modeled by space frame elements, they are usually analyzed using Hermite beam elements characterized by shape functions consisting of Hermite interpolation functions (cubic for bending, linear for axial components of displacement and for the angle of twist) [41] or isoparametric beam elements [42]. In these elements the variation of the cross section is generally considered by setting “average” values for the cross section parameters resulting from the corresponding parameters at the start and end nodes of the element. Moreover, application of shape functions results in inaccuracies in stress analysis of beams of variable cross section, as static and kinematic values at nodes and in the element region are computed only approximately and the element may not satisfy local and global equilibrium conditions [31]. This inaccuracy may be reduced by increasing the number of integration points for the stiffness matrix assembly (at least three-point integration) or by refining the mesh of elements, increasing at the same time the effort of preparing the input parameters.Nevertheless, Sapountzakis and Mokos in [43–45] developed a boundary element solution for the construction of the 14×14 stiffness matrix and the nodal load vector of a member of arbitrary homogeneous or composite, constant or variable cross section, subjected to an arbitrarily concentrated or distributed twisting moment, taking into account warping effects. The arbitrary low rate variation of the member of variable cross section is continuous so as to assume that its shear center is independent of the member loading and boundary conditions. The developed method can take into account the variable torsional and warping rigidities along the member length. In these latter research efforts boundary value problems with respect to the variable along the bar angle of twist and to the primary warping function are formulated and solved employing a pure BEM [23] approach; that is only boundary discretization is used.In order to formulate the aforementioned problem, let us consider a prismatic bar of lengthl (Figure 7), of constant or variable arbitrary cross-section. The homogeneous isotropic and linearly elastic material of the bar’s cross-section, with modulus of elasticity E, shear modulus G, and Poisson’s ratio ν, occupies the two dimensional multiply connected region Ω of the y, z plane and is bounded by the Γj(j=1,2,…,K) boundary curves, which are piecewise smooth; that is, they may have a finite number of corners. In Figure 7(b)CYZ is the principal coordinate system through the cross section’s centroid C, while yC,zC are its coordinates with respect to Syz system of axes through the cross section’s shear center S.Beam element (a) of a cross section of arbitrary shape occupying the two dimensional regionΩ (b). Cy~z~ is the beam principal centroidal system of axes, while Syz is the corresponding one through the shear center S.
(a)
(b)In order to include the warping behavior in the study of the aforementioned element, in each node at the element ends a seventh degree of freedom is added to the well-known six DOFs of the classical three-dimensional frame element. The additional DOF is the first derivative of the angle of twistθx′=dθx/dx denoting the rate of change of the angle of twist θx, which can be regarded as the torsional curvature (Figure 8) of the cross section. Thus, the nodal displacement vector in the local coordinate system, as shown in Figure 7(a), can be written as (46){Di}T={ux~juy~juz~jθxjθy~jθz~jθxj′ux~kuy~kuz~kθxkθy~kθz~kθxk′}and the respective nodal load vector as(47){Fi}T={NjQy~jQz~jMtjMy~jMz~jMwjNkQy~kQz~kMtkMy~kMz~kMwk},where Mt is the twisting moment at the ends of the element given from (29) and consisting of the primary part MtP defined as the resultant of the primary shear stress distribution given from (31a) or (44a) and the secondary one MtS defined as the resultant of the secondary shear stress distribution due to warping given from (31b) or (44b) for the cases of constant or variable cross section, respectively, that is,(48a)Mt=GItθx′-ECSθx′′′fortheconstantcrosssectioncase,(48b)Mt=GItθx′-ECS′θx′′-ECSθx′′′forthevariablecrosssectioncase.Moreover, Mw is the bending moment due to the torsional curvature at the same sections given from (38).Figure 8
Torsional curvature of a rectangular and a hollow square cross section.The nodal displacement and load vectors given in (46) and (47) are related with the 14×14 local stiffness matrix of the spatial beam element written as(49)[ki]=[k11i000000k18i0000000k22i000k26i00k29i000k2,13i000k33i0k35i0000k3,10i0k3,12i00000kT1i00kT2i000-kT1i00kT5i00k53i0k55i0000k5,10i0k5,12i000k62i000k66i00k69i000k6,13i0000kT2i00kT3i000-kT2i00kT4ik81i000000k88i0000000k92i000k96i00k99i000k9,13i000k10,3i0k10,5i0000k10,10i0k10,12i00000-kT1i00-kT2i000kT5i00-kT5i00k12,3i0k12,5i0000k12,10i0k12,12i000k13,2i000k13,6i00k13,9i000k13,13i0000kT5i00kT4i000-kT5i00kT6i],where the klmi coefficients (l,m=1,2,3,5,6,8,9,10,12,13) of the stiffness matrix of (49) come from the well-known 12×12 classical stiffness matrix of the classical three-dimensional frame element. In the special case of a constant cross section element the kTni(n=1,2,3,4,5,6) coefficients are given as [14, 35](50a)kT1i=koλS,(50b)kT2i=kT5i=ko(C-1),(50c)kT3i=kT6i=ko(Cl-1λS),(50d)kT4i=ko(1λS-l),where(51a)λ=GItECS,(51b)S=sinhe,(51c)C=coshe,(51d)e=lλ,(51e)ko=GIt2(1-C)+eS.In the case of a variable cross section element, the evaluation of the kTni(n=1,2,3,4,5,6) coefficients presumes the solution of the following boundary value problem with respect to the angle of twist θx=θx(x) (see also (42)):
(52)ECSθx′′′′+2ECS′θx′′′+(ECS′′-GIt)θx′′-GIt′θx′=0insidetheelement,(53a)α1θx+α2Mt=α3,(53b)β1θx′+β2Mw=β3attheelementendsx=0,lfor appropriate values of the ai, βi(i=1,2,3) functions. Thereby, for the evaluation of the kT1i coefficient it is a1=a3=β1=1,  a2=β2=β3=0 at x=0 and a1=β1=1, a2=a3=β2=β3=0 at x=l, for the evaluation of the kT3i coefficient it is a1=β1=β3=1, a2=a3=β2=0 at x=0 and a1=β1=1, a2=a3=β2=β3=0 at x=l, while for the evaluation of the kT6i coefficient it is a1=β1=1, a2=a3=β2=β3=0 at x=0 and a1=β1=β3=1, a2=a3=β2=0 at x=l. Apparently, for the constant cross section case it is CS′=CS′′=It′=0 and the governing equation (52) reduces to the one given from (35). Upon the evaluation of the angle of twist θx, the coefficients kTni(n=1,2,3,4,5,6) are established from its derivatives using relations (38), (48a), and (48b).According to the nodal load vector, assuming that the span of the bar is subjected to the arbitrarily concentrated or distributed twisting momentmt=mt(x) (Figure 7(a)), the evaluation of the elements concerning the twisting and the bending moments due to the torsional curvature is accomplished using again relations (38), (48a), and (48b) employing the derivatives of the angle of twist θx, obtained from the solution of the following boundary value problem:for the constant cross section case:
(54)ECSθx′′′′-GItθx′′=mtinsidethebar,θx=θx′=0atthebarendsx=0,l,
for the variable cross section case:
(55)ECSθx′′′′+2ECS′θx′′′+(ECS′′-GIt)θx′′-GIt′θx′=mtinsidethebar,θx=θx′=0atthebarendsx=0,l.From the examined example problems presented in [43–45] it is concluded that (i)
the discrepancy of the deflections and the internal stress resultants arising from the ignorance of the warping degrees of freedom at the ends of a member and the magnitude of the normal stresses due to warping necessitate the utilization of the14×14 member stiffness matrix, especially for beams with open shaped cross sections;(ii)
the advantages of a box shaped closed cross section beam subjected to torsional loading compared with that of an open one are verified;(iii)
warping is not constant along the thickness of the cross section walls as it is assumed in thin tube theory for thin-walled beams;(iv)
comparison of the elements of the resulting stiffness matrix of a variable cross section member taking into account the derivatives of the variable torsional and warping rigidities along the member length with those obtained using a fine mesh of elements having “average” values for the cross section parameters leads to the consideration of these derivatives;(v)
the discrepancy between the uncoupled proposed procedure and the formulation that takes into account coupled displacement components is inconsiderable;(vi)
having in mind the previous conclusion and that coupled displacement components lead to dependent shear center on the member loading and boundary conditions, and the ignorance of the effect of the transverse displacement components to the longitudinal one is justified.


## 5. Elastic Linear Torsional Vibrations of Constant or Variable Cross Section Bars

In engineering practice, we often come across the analysis of structures subjected to vibratory twisting loading. The dynamic forces acting on a structure may result from one or more of different causes, such as rotating machinery, wind, asymmetric traffic loading, blast loads, or earthquake forces. The extensive use of the aforementioned structural elements necessitates a rigorous dynamic analysis.Exact torsional vibration frequencies were presented by Gorman [46] and Belvins [47] for the case of circular cross section shafts subjected to classical boundary conditions, avoiding in this way warping effects. These efforts were extended by Kameswara Rao [48] for elastically restrained edges. Torsional vibration frequencies for beams of open thin-walled sections, subjected to several combinations of classical boundary conditions, taking into account warping effects were first derived by Gere [49]. Since then approximate methods [50] for the calculation of natural frequencies including elastic torsional and warping restraints [51] and employing either discrete [52–55] or distributed mass model systems [56–60] have been presented. In all these references the considered beam is of a constant homogeneous thin-walled cross section, while its torsion and warping constants are evaluated employing the relations of the thin tube theory.Several researchers have also dealt with beams of variable cross section ignoring the warping effects resulting from the corresponding restraints at the ends of the member [26, 27]. On the contrary, to the author’s knowledge relatively little work has been done on the problem of nonuniform torsion of bars of variable cross section with pioneer the work of Cywinski [28] adopting the finite difference method. Wekezer [29] after dividing the bar into segments along its longitudinal axis approximated their shell midsurface by arbitrary triangular shell elements and employed the finite element method to the linear membrane shell theory. This approximation generates inaccuracies, as the warping of the walls of the cross section cannot be taken into account. Moreover, Eisenberger [30, 61] employed FEM upon polynomial approximation of the torsional and warping rigidities using “exact” shape functions to derive the exact stiffness coefficients. This application of shape functions results also in inaccuracies in stress analysis of beams of variable cross section, as static and kinematic values at nodes and in the element region are computed only approximately and the element may not satisfy local and global equilibrium conditions [31]. In all the aforementioned procedures the torsion and warping constants have been approximated adopting the thin tube theory. Moreover, research efforts have been presented for the corresponding problem of composite beams limited to the formulation of a displacement-based one-dimensional finite element model for the estimation of natural frequencies and corresponding mode shapes of thin-walled composite beams after approximating again the torsion and warping constants with closed form solutions [62, 63].Only in Ganapathi et al. [64] the warping function for the constant rectangular cross section of sandwich beams is determined by solving the boundary value problem for torsion such that the displacements are continuous at the interfaces of adjacent layers, while the transverse shear stress is continuous at these interfaces and vanishes at the top and bottom surfaces of the beam. Moreover, Sapountzakis in [65, 66] developed a boundary element method for the nonuniform torsional vibration problem of doubly symmetric composite bars of arbitrary constant or variable cross section, respectively. In these latter efforts the beam is subjected to an arbitrarily distributed dynamic twisting moment, while its edges are restrained by the most general linear torsional boundary conditions. A distributed mass model system is employed which leads to the formulation of three boundary value problems with respect to the variable along the beam angle of twist and to the primary and secondary warping functions. The last two problems are solved employing a pure BEM [23] approach that is only boundary discretization is used. Finally, Sapountzakis and Mokos in [67] presented the dynamic analysis of 3D beam elements restrained at their edges by the most general linear torsional, transverse, or longitudinal boundary conditions and subjected in arbitrarily distributed dynamic twisting, bending, transverse, or longitudinal loading. For the solution of this problem, a boundary element method is developed for the construction of the 14×14 stiffness matrix and the corresponding nodal load vector, of a member of an arbitrarily shaped simply or multiply connected cross section, taking into account both warping and shear deformation effects, which together with the respective mass and damping matrices lead to the formulation of the equation of motion.In order to formulate the nonuniform torsional vibration problem of doubly symmetric bars of arbitrarily shaped simply or multiply connected constant or variable cross section, let us consider the bar of lengthl of Figure 9. The homogeneous isotropic and linearly elastic material of the bar’s cross-section, with modulus of elasticity E, shear modulus G, and Poisson’s ratio ν, occupies the two dimensional multiply connected region Ω of the y, z plane and is bounded by the Γj(j=1,2,…,K) boundary curves, which are piecewise smooth; that is, they may have a finite number of corners. In Figure 9(b)CYZ is the principal coordinate system through the cross section’s centroid C, while yC, zC are its coordinates with respect to Syz system of axes through the cross section’s shear center S. The bar is subjected to the time dependent arbitrarily distributed or concentrated twisting moment  mt=mt(x,t), t≥0 acting in the x direction (Figure 9(a)).Bar subjected to a time dependent twisting moment (a) with a variable cross section of arbitrary shape occupying the two dimensional regionΩ (b).
(a)
(b)Adopting the same displacement field (21a), (21b), and (21c) with that of the static analysis of the constant cross section case, which for the dynamic problem is written as(56a)u(x,y,z,t)=θx′(x,t)·ϕSP(y,z)+ϕSS(x,y,z,t),(56b)v(x,y,z,t)=-zθx(x,t),(56c)w(x,y,z,t)=yθx(x,t)and following the same procedure presented for the constant cross section bar in Section 2, the following initial boundary value problem for the angle of twist θx=θx(x) is derived:
(57)ECSd4θxdx4+2EdCSdxd3θxdx3+(Ed2CSdx2-GIt)d2θxdx2-GdItdxdθxdx+ρIpd2θxdt2=mtinsidethebeam,(58a)α1θx+α2Mt=α3,(58b)β1dθxdx+β2Mw=β3atthebeamendsx=0,l,(59a)
θ
x
(
x
,
0
)
=
ϑ
-
(
x
)
,
(59b)
θ
˙
x
(
x
,
0
)
=
ϑ
-
˙
(
x
)
,
where θ˙x(x,t) is the first derivative of the angle of twist with respect to time; ϑ(x), ϑ-(x) are the initial angle of twist and the corresponding initial velocity of the points of the beam axis; Ip=Ip(x) is the polar moment of inertia of the cross section about the origin S of its two axes of symmetry (see Figure 9) and ρ is the mass density of the cross section material. It is worth here noting that in the constant cross section case the governing equation (57) reduces to
(60)ECSd4θxdx4-GItd2θxdx2+ρIpd2θxdt2=mtinsidethebeam.
In (58a) and (58b) the total twisting moment of the cross section is once again divided into a primary and a secondary component as this is stated in (29), where these components are given from (44a) and (44b), the resulting total twisting moment of the cross section is given from (45), while the warping moment Mw arising from the torsional curvature, similarly with the constant cross section case is given from (38). Also, the functions ai, βi(i=1,2,3) are specified at the boundary of the beam forming the most general linear torsional boundary conditions for the beam problem including also the elastic support. Moreover, as for the static constant cross section problem the determination of the primary ϕSP(y,z) and the secondary ϕSS(x,y,z,t) warping functions is achieved from the solution of the boundary value problems given from (19a), (19b), (20a), and (20b), respectively, noting that the secondary warping function in this case is time dependent. Finally, the primary, the secondary shear stress components, and the normal stresses due to warping are defined similarly with relations (7a)–(9), which for the dynamic problem are written as(61a)τxyP(x,y,z,t)=G·θx′(x,t)·(∂ϕSP∂y-z),(61b)τxyS(x,y,z,t)=G·∂ϕSS∂y,(62a)
τ
x
z
P
(
x
,
y
,
z
,
t
)
=
G
·
θ
x
′
(
x
,
t
)
·
(
∂
ϕ
S
P
∂
z
+
y
)
,
(62b)
τ
x
z
S
(
x
,
y
,
z
,
t
)
=
G
·
∂
ϕ
S
S
∂
z
,
(63)
σ
x
x
w
(
x
,
y
,
z
,
t
)
=
E
·
θ
x
′′
(
x
,
t
)
·
ϕ
S
P
(
y
,
z
)
.From the examined example problems presented in [65–67] it is concluded that (i)
the discrepancy of the dynamic response of the bar arising from the ignorance of the warping effect necessitates the inclusion of the nonuniform torsional beam behavior;(ii)
the discrepancy of the results arising from the ignorance of the warping especially in higher eigenfrequencies is remarkable;(iii)
the discrepancy in the analysis of a thin-walled cross section beam employing the BEM after calculating the torsion and warping constants adopting the thin tube theory demonstrates the importance of the proposed procedure even in thin-walled beams, since it approximates better the torsion and warping constants and takes also into account the warping of the walls of the cross section;(iv)
the discrepancy of the results arising from the ignorance of the warping degrees of freedom at the ends of a member necessitates the utilization of the14×14 member stiffness matrix, especially for beams with open shaped cross section;(v)
warping is not constant along the thickness of the cross section walls as it is assumed in thin tube theory for thin-walled beams.


## 6. Secondary Torsional Moment Deformation Effect in Elastic Linear Torsional Analysis of Bars

As it has already been mentioned in the previous sections in engineering practice we often come across the analysis of rectilinear members of structures subjected to nonuniform torsion. In this case, as it has been analyzed in Section2, the twisting moment is split into a primary and a secondary part, where the primary and the secondary torsion moments are undertaken from the Saint-Venant shear stresses (primary shear stresses) and the warping shear stresses (secondary shear stresses), respectively, while the warping normal stresses are undertaken from the warping moment (bimoment). Moreover, in engineering practice as well as in the literature very strong torsional warping is assumed to occur only for open shaped cross sections, while for closed shaped ones the warping effect is assumed to be insignificant and therefore negligible. However, this assumption is not always valid especially for bars of closed shaped cross sections and small length. Noting that in the case of direct torsion (equilibrium torsion) the arising normal and additional shear stresses due to warping are equilibrium and not constraint stresses, which means that after a crack these stresses are redistributed, the necessity of the evaluation and inclusion of these additional warping stresses in the analysis is clearly evident.Despite the wide study from both the analytical [14, 15, 68–71] and numerical [18–22, 24, 72–75] point of view of the nonuniform torsion problem of prismatic bars, relatively little work has been done on the corresponding problem considering secondary torsional moment deformation effect, that is shear deformation due to the nonuniform torsional warping. It is worth here noting that if shear deformations due to shear forces and restrained warping are considered, flexure and torsion of bars of nonsymmetrical cross section are generally coupled, even if the bar is subjected only to twisting loading, while only for bars of doubly symmetrical cross section flexure and torsion are decoupled. According to the research efforts on the analytical solution of the aforementioned nonuniform torsion problem including secondary shear deformation effect the pioneer work of Heilig [76, 77] is mentioned, in which a theoretical formulation of the problem is presented. Later, Roik and Sedlacek [78] presented an analytical solution applying the Force Method (Flexibility Method) and employing the analogy between 2nd order beam theory with tensional axial force and torsion with warping. Moreover, in the work of Schade [79] a theoretical formulation of the corresponding problem is presented including the coupled shear deformation effect due to shear forces and restrained warping. Rubin [80] using the methodology of Roik and Sedlacek [78] presented an analytical solution for continuous prismatic bars by introducing a Three-Moment Equation (similar with the Three-Moment Equation of Clapeyron for the flexure problem). Also, the theory of the shear deformations due to the restrained torsional warping has been validated by a series of torsional experimental tests on fibre reinforced plastics I-beams by Roberts and Al-Ubaidi [81]. Furthermore, from the numerical point of view intensive research works have been made over the last years to develop new beam finite elements which take into account coupled [82] or uncoupled [83] secondary shear deformation. However, in all of the aforementioned research efforts only bars of thin-walled cross section are investigated, since the torsional cross section parameters (warping constant, primary and secondary torsion constant, and shear center) as well as the shear and normal stresses are evaluated employing the Vlasov thin tube theory [68, 69]. Finally, Kraus [84] presented a FEM solution for the calculation of the secondary torsion constant of bars of hot rolled I-sections and Mokos and Sapountzakis in [85] developed a BEM solution for the nonuniform torsion of simply or multiply connected bars of doubly symmetric arbitrary constant cross section (thin and/or thick walled), taking into account secondary torsional moment deformation effect. In this last research effort, to account for secondary shear deformations, the concept of shear deformation coefficient is used leading to a secondary torsion constant, which is computed employing an effective automatic domain integration using the Advancing Front Method (AFM) [86].In order to formulate the aforementioned problem, let us consider a rectilinear bar of lengthl, of an arbitrary doubly symmetric constant cross section with modulus of elasticity E and shear modulus G, occupying the two dimensional multiply connected region Ω of the y, z plane bounded by the K+1 curves Γ0,Γ1,…,ΓK as shown in Figure 10. The bar is subjected to arbitrarily distributed and/or concentrated twisting mt=mt(x) and warping mw=mw(x) moments, while its edges are restrained by the most general torsional boundary conditions.Rectilinear bar subjected to twisting and warping moments (a) of constant doubly symmetric arbitrary cross section occupying the two dimensional regionΩ (b).
(a)
(b)In order to account for shear deformation due to restrained torsional warping, the displacementsux, uy, uz along the three axes and the rotation θx due to twist are spilt into a primary and a secondary part arising from the primary (P) and the secondary (S) torsional moment, respectively, as(64a)uy=uyP+uyS,(64b)uz=uzP+uzS,(64c)ux=uxP+uxS,(65a)
θ
x
=
θ
x
P
+
θ
x
S
,
(65b)
d
θ
x
d
x
=
d
θ
x
P
d
x
+
d
θ
x
S
d
x
,
where the (total) first derivative of the angle of twist dθx/dx denotes the rate of change of the (total) angle of twist and can be regarded as the (total) torsional curvature. The primary dθxP/dx and the secondary dθxS/dx part of the first order derivative of the angle of twist represent the primary torsional curvature (Figure 11(a)) arising from the (primary) warping normal stresses and the secondary torsional curvature (Figure 11(b)) arising from the (secondary) warping shear stresses, respectively.Primary (a) and secondary (b) torsional curvature of a hot rolled I-profile steel HEB-100 cross section.
(a)
(b)Assuming small torsional rotation, the displacement field is defined as(66a)uyP=-zθx(x),(66b)uzP=yθx(x),(66c)uxP=dθxP(x)dxφSP(y,z),(66d)uyS=uzS=uxS=0,where φSP(y,z) is the primary warping function with respect to the shear center S of the cross section of the bar. Substituting the aforementioned displacement field to the linearized strain-displacement relations (Cauchy strain tensor) and to the stress-strain relations (constitutive relations) of the three-dimensional elasticity the nonzero stress components in the region Ω can be written as(67a)σxxP=Ed2θxP(x)dx2φSP(y,z),(67b)τxyP=Gdθx(x)dx(∂φSP(y,z)∂y-z),(67c)τxzP=Gdθx(x)dx(∂φSP(y,z)∂z+y),(68a)
τ
x
y
S
=
G
∂
φ
S
S
(
x
,
y
,
z
)
∂
y
,
(68b)
τ
x
z
S
=
G
∂
φ
S
S
(
x
,
y
,
z
)
∂
z
,
where φSS(x,y,z) is the secondary warping function with respect to the shear center S of the cross section of the bar. It is worth here mentioning that, for the derivation of the aforementioned relations of shear stresses, the arising terms (dθxS/dx)(∂φSP/∂y), (dθxS/dx)(∂φSP/∂z) have been neglected.


### 6.1. Equations of Local Equilibrium

The first elasticity equation of equilibrium of the three-dimensional elasticity (equilibrium in the axial direction of the bar) with vanishing body forces is written as(69)(∂τxyP∂y+∂τxzP∂z)P+(∂σxxP∂x+∂τxyS∂y+∂τxzS∂z)S=0.
Similarly with Section 2, requiring both the primary (P) and the secondary due to warping (S) parts of (69) to vanish, as well as the corresponding ones of the traction vector on the free lateral surface of the bar, the Neumann problems of (19a) and (19b) for the primary and of (20a) and (20b) for the secondary warping functions are obtained.Additionally, applying the stress components given from (67a), (67b), (67c), (68a), and (68b) into the stress resultants given from (29), (30a), (30b), and (36) and taking into account the Neumann problems (19a), (19b), (20a), and (20b), the following relations for the nonzero stress resultants of the bar are obtained:
(70)Mw=-ECSd2θxPdx2,(71a)MtP=GItPdθxdx,(71b)MtS=-ECSd3θxPdx3,where(72a)CS=∫Ω(ϕSP)2dΩ,(72b)ItP=∫Ω(y2+z2+y·∂ϕSP∂z-z·∂ϕSP∂y)dΩare the warping and the primary torsion constants of the cross section, respectively, while ECS and GItP are the warping and the primary torsional rigidities of the cross section, respectively. It is worth noting that the relations (70), (71a), and (71b) are valid inside the bar as well as at the bar ends. Moreover, observing (70) and (71b) equation (40) can be easily obtained, while substituting (71b) into (20a) the partial Poisson type differential equation governing the secondary warping function is obtained as
(73)∇2φSS=MtSGCSφSP(y,z)inΩ.
Since the secondary warping is much smaller with regard to the primary one [24], the secondary torsional curvature dθxS/dx can be taken into account in the calculation of the angle of twist indirectly using an effective secondary torsion constant. Thus, as in Timoshenko’s beam theory for shear deformable beams [87, 88], the secondary torsional curvature can be approximately evaluated from the following relation [76, 77]:
(74)dθxSdx=MtSGItS,
where GItS is the secondary torsional rigidity of the cross section, while ItS is the secondary torsion constant of the cross section which can be written as
(75)ItS=κxSItP=1αxSItP,
where the coefficient κxS is called warping shear correction factor and αxS is the warping shear deformation coefficient. Substituting (75) into (74) the secondary torsional curvature can be written as
(76)dθxSdx=αxSMtSGItP.
From (76) it is concluded that for αxS=0 the secondary shear deformation effect is neglected, that is dθxS/dx=0 and dθxP/dx=dθx/dx. This assumption is usually employed in the case of open shaped cross sections.The evaluation of the secondary torsion constant of the cross section can be achieved by applying an energy approach [16, 78]. Thus, taking into account (68a), (68b), and (74) and equating the approximate formula for the evaluation of the secondary shear strain energy per unit length given from
(77)Uappr=12MtSdθxSdx=12(MtS)2GItS
with the exact one given from
(78)Uexact=12G∫Ω(τxyS2+τxzS2)dΩ=12G∫Ω(∂φSS∂y∂φSS∂y+∂φSS∂z∂φSS∂z)dΩ
the secondary torsion constant can be obtained as(79)ItS=(MtS)2G2×1∫Ω((∂φSS/∂y)(∂φSS/∂y)+(∂φSS/∂z)(∂φSS/∂z))dΩ.
In order to formulate this constant independently from the loading and the material properties of the bar, the values MtS=1, G=1 are employed and the secondary torsion constant is given as
(80)ItS=1∫Ω((∂φ~SS/∂y)(∂φ~SS/∂y)+(∂φ~SS/∂z)(∂φ~SS/∂z))dΩ,
where φ~SS(y,z) is the unit secondary warping function with respect to the shear center S of the cross section of the bar, while taking into account (73) and (20b) the aforementioned function φ~SS can be established by solving independently the following Neumann problem:(81a)∇2φ~SS=φSPCSinΩ,(81b)∂φ~SS∂n=0onΓ.Moreover, applying the Green identity for the warping functions φSP(y,z) and φ~SS(y,z) and taking into account the Neumann problems (19a), (19b), (81a), (81b), and (80) the secondary torsion constant is given as
(82)ItS=CSIφ,
where Iφ is a domain integral given from the relation
(83)Iφ=-∫ΩφSPφ~SSdΩ
while form (75) and (82) the following relation for the warping shear deformation coefficient is obtained:
(84)αxS=ItPCSIφ.


### 6.2. Equations of Global Equilibrium

Furthermore, in order to formulate the governing differential equation of moment equilibrium of the bar in the axial direction, the equilibrium of a small segmentdx of the bar (Figure 12) is considered leading to the relation
(85)dMtdx=-mt-dmwdx.
Taking into account (29), (85) can be written as
(86)dMtSdx+dMtPdx=-mt-dmwdx
while substituting (40) into (86) the following equilibrium equation is obtained:
(87)d2Mwdx2+dMtPdx=-mt-dmwdx.
Having in mind that both the secondary twisting moment MtS(θxP) and the bimoment Mw(θxP) (given from relations (71b) and (70), resp.) are expressed in terms of the primary angle of twist (θxP), while the primary twisting moment MtP(θx) (given from relation (71a)) is expressed in terms of the total angle of twist (θx), an effort is given in the following to formulate also the primary twisting moment MtP with respect to the primary angle of twist (θxP).Figure 12
Equilibrium of a small segmentdx of a twisted bar.Introducing an auxiliary geometric constantκ defined as
(88)κ=11+αxS
the relation (76) can be written as
(89)dθxSdx=1-κκMtSGItP
while from relation (65b) the total torsional curvature can be written as
(90)dθxdx=dθxPdx+1-κκMtSGItP
and according to (71a) and (71b), the primary twisting moment is given with respect to the primary angle of twist (θxP) as
(91)MtP=GItPdθxPdx-EC-Sd3θxPdx3(1-κ),
where C-S is a modified warping constant given as
(92)C-S=CSκ.
It is worth here noting that the auxiliary constant κ is always smaller or equal to one, that is κ≤1. Small values of the constant κ indicate that the secondary torsional moment deformation effect is important and should be considered in the analysis, while in the case of negligible secondary shear deformations κ=1. Moreover, according to (29), (71b), and (91) the (total) twisting moment Mt is given with respect to the primary angle of twist (θxP) from the relation
(93)Mt=GItPdθxPdx-EC-Sd3θxPdx3
while using (92) the moments MtS and Mw can be written in terms of the modified warping constant as
(94)Mw=-EC-Sd2θxPdx2κ,(95)MtS=-EC-Sd3θxPdx3κ.
Hence, substituting (91) and (94) into (87) the governing differential moment equation of equilibrium of the bar in the axial direction with respect to the primary angle of twist is obtained as
(96)EC-Sd4θxPdx4-GItPd2θxPdx2=+mt+dmwdxinsidethebar
while the boundary conditions are given from(97a)α1θxP+α2Mt=α3,(97b)β1dθxPdx+β2Mw=β3atthebarendsx=0,l,where the parameters ai, βi(i=1,2,3) are functions specified at the boundary of the bar and the moments Mt and Mw are given form the relations (93) and (94), respectively. It is worth here noting that the boundary conditions (97a) and (97b) are the most general linear torsional boundary conditions including also the elastic support. It is apparent that all types of the conventional torsional boundary conditions (clamped, simply supported, free, or guided edge) can be derived from these equations by specifying appropriately the functions ai and βi (e.g., for a clamped edge it is a1=β1=1, a2=a3=β2=β3=0).Furthermore, substituting (95) into (89) and differentiating with respect to x, the governing differential equation of equilibrium of the bar in the axial direction with respect to the secondary angle of twist is obtained in terms of the primary angle of twist as
(98)d2θxSdx2=-EC-SGItPd4θxPdx4(1-κ)insidethebar
while the corresponding boundary conditions are given in terms of the primary angle of twist as
(99)γ1θxS+γ2(GItPκ1-κ)dθxSdx=γ3(-EC-Sκd3θxPdx3)atthebarendsx=0,l,
where if the primary angle of twist is fixed, that is a1=1 and a2=a3=0, then γ1=1 and γ2=γ3=0, otherwise γ1=0 and γ2=γ3=1.From the examined example problems presented in [85] it is concluded that (i)
the inaccuracy of the thin tube theory in calculating the secondary torsion constant even for thin-walled sections is noteworthy;(ii)
as it is also verified by other computational methods [80, 83], for bars with open shaped cross section the secondary torsional moment deformation effect has not significant influence. However, the inclusion of this additional effect leads to more accurate results;(iii)
as it is also verified by other computational methods [80, 83], for bars with closed shaped cross sections the secondary torsional moment deformation effect has an important influence and should be considered in the analysis (otherwise the results may be completely wrong especially in the calculation of stresses).


## 6.1. Equations of Local Equilibrium

The first elasticity equation of equilibrium of the three-dimensional elasticity (equilibrium in the axial direction of the bar) with vanishing body forces is written as(69)(∂τxyP∂y+∂τxzP∂z)P+(∂σxxP∂x+∂τxyS∂y+∂τxzS∂z)S=0.
Similarly with Section 2, requiring both the primary (P) and the secondary due to warping (S) parts of (69) to vanish, as well as the corresponding ones of the traction vector on the free lateral surface of the bar, the Neumann problems of (19a) and (19b) for the primary and of (20a) and (20b) for the secondary warping functions are obtained.Additionally, applying the stress components given from (67a), (67b), (67c), (68a), and (68b) into the stress resultants given from (29), (30a), (30b), and (36) and taking into account the Neumann problems (19a), (19b), (20a), and (20b), the following relations for the nonzero stress resultants of the bar are obtained:
(70)Mw=-ECSd2θxPdx2,(71a)MtP=GItPdθxdx,(71b)MtS=-ECSd3θxPdx3,where(72a)CS=∫Ω(ϕSP)2dΩ,(72b)ItP=∫Ω(y2+z2+y·∂ϕSP∂z-z·∂ϕSP∂y)dΩare the warping and the primary torsion constants of the cross section, respectively, while ECS and GItP are the warping and the primary torsional rigidities of the cross section, respectively. It is worth noting that the relations (70), (71a), and (71b) are valid inside the bar as well as at the bar ends. Moreover, observing (70) and (71b) equation (40) can be easily obtained, while substituting (71b) into (20a) the partial Poisson type differential equation governing the secondary warping function is obtained as
(73)∇2φSS=MtSGCSφSP(y,z)inΩ.
Since the secondary warping is much smaller with regard to the primary one [24], the secondary torsional curvature dθxS/dx can be taken into account in the calculation of the angle of twist indirectly using an effective secondary torsion constant. Thus, as in Timoshenko’s beam theory for shear deformable beams [87, 88], the secondary torsional curvature can be approximately evaluated from the following relation [76, 77]:
(74)dθxSdx=MtSGItS,
where GItS is the secondary torsional rigidity of the cross section, while ItS is the secondary torsion constant of the cross section which can be written as
(75)ItS=κxSItP=1αxSItP,
where the coefficient κxS is called warping shear correction factor and αxS is the warping shear deformation coefficient. Substituting (75) into (74) the secondary torsional curvature can be written as
(76)dθxSdx=αxSMtSGItP.
From (76) it is concluded that for αxS=0 the secondary shear deformation effect is neglected, that is dθxS/dx=0 and dθxP/dx=dθx/dx. This assumption is usually employed in the case of open shaped cross sections.The evaluation of the secondary torsion constant of the cross section can be achieved by applying an energy approach [16, 78]. Thus, taking into account (68a), (68b), and (74) and equating the approximate formula for the evaluation of the secondary shear strain energy per unit length given from
(77)Uappr=12MtSdθxSdx=12(MtS)2GItS
with the exact one given from
(78)Uexact=12G∫Ω(τxyS2+τxzS2)dΩ=12G∫Ω(∂φSS∂y∂φSS∂y+∂φSS∂z∂φSS∂z)dΩ
the secondary torsion constant can be obtained as(79)ItS=(MtS)2G2×1∫Ω((∂φSS/∂y)(∂φSS/∂y)+(∂φSS/∂z)(∂φSS/∂z))dΩ.
In order to formulate this constant independently from the loading and the material properties of the bar, the values MtS=1, G=1 are employed and the secondary torsion constant is given as
(80)ItS=1∫Ω((∂φ~SS/∂y)(∂φ~SS/∂y)+(∂φ~SS/∂z)(∂φ~SS/∂z))dΩ,
where φ~SS(y,z) is the unit secondary warping function with respect to the shear center S of the cross section of the bar, while taking into account (73) and (20b) the aforementioned function φ~SS can be established by solving independently the following Neumann problem:(81a)∇2φ~SS=φSPCSinΩ,(81b)∂φ~SS∂n=0onΓ.Moreover, applying the Green identity for the warping functions φSP(y,z) and φ~SS(y,z) and taking into account the Neumann problems (19a), (19b), (81a), (81b), and (80) the secondary torsion constant is given as
(82)ItS=CSIφ,
where Iφ is a domain integral given from the relation
(83)Iφ=-∫ΩφSPφ~SSdΩ
while form (75) and (82) the following relation for the warping shear deformation coefficient is obtained:
(84)αxS=ItPCSIφ.


## 6.2. Equations of Global Equilibrium

Furthermore, in order to formulate the governing differential equation of moment equilibrium of the bar in the axial direction, the equilibrium of a small segmentdx of the bar (Figure 12) is considered leading to the relation
(85)dMtdx=-mt-dmwdx.
Taking into account (29), (85) can be written as
(86)dMtSdx+dMtPdx=-mt-dmwdx
while substituting (40) into (86) the following equilibrium equation is obtained:
(87)d2Mwdx2+dMtPdx=-mt-dmwdx.
Having in mind that both the secondary twisting moment MtS(θxP) and the bimoment Mw(θxP) (given from relations (71b) and (70), resp.) are expressed in terms of the primary angle of twist (θxP), while the primary twisting moment MtP(θx) (given from relation (71a)) is expressed in terms of the total angle of twist (θx), an effort is given in the following to formulate also the primary twisting moment MtP with respect to the primary angle of twist (θxP).Figure 12
Equilibrium of a small segmentdx of a twisted bar.Introducing an auxiliary geometric constantκ defined as
(88)κ=11+αxS
the relation (76) can be written as
(89)dθxSdx=1-κκMtSGItP
while from relation (65b) the total torsional curvature can be written as
(90)dθxdx=dθxPdx+1-κκMtSGItP
and according to (71a) and (71b), the primary twisting moment is given with respect to the primary angle of twist (θxP) as
(91)MtP=GItPdθxPdx-EC-Sd3θxPdx3(1-κ),
where C-S is a modified warping constant given as
(92)C-S=CSκ.
It is worth here noting that the auxiliary constant κ is always smaller or equal to one, that is κ≤1. Small values of the constant κ indicate that the secondary torsional moment deformation effect is important and should be considered in the analysis, while in the case of negligible secondary shear deformations κ=1. Moreover, according to (29), (71b), and (91) the (total) twisting moment Mt is given with respect to the primary angle of twist (θxP) from the relation
(93)Mt=GItPdθxPdx-EC-Sd3θxPdx3
while using (92) the moments MtS and Mw can be written in terms of the modified warping constant as
(94)Mw=-EC-Sd2θxPdx2κ,(95)MtS=-EC-Sd3θxPdx3κ.
Hence, substituting (91) and (94) into (87) the governing differential moment equation of equilibrium of the bar in the axial direction with respect to the primary angle of twist is obtained as
(96)EC-Sd4θxPdx4-GItPd2θxPdx2=+mt+dmwdxinsidethebar
while the boundary conditions are given from(97a)α1θxP+α2Mt=α3,(97b)β1dθxPdx+β2Mw=β3atthebarendsx=0,l,where the parameters ai, βi(i=1,2,3) are functions specified at the boundary of the bar and the moments Mt and Mw are given form the relations (93) and (94), respectively. It is worth here noting that the boundary conditions (97a) and (97b) are the most general linear torsional boundary conditions including also the elastic support. It is apparent that all types of the conventional torsional boundary conditions (clamped, simply supported, free, or guided edge) can be derived from these equations by specifying appropriately the functions ai and βi (e.g., for a clamped edge it is a1=β1=1, a2=a3=β2=β3=0).Furthermore, substituting (95) into (89) and differentiating with respect to x, the governing differential equation of equilibrium of the bar in the axial direction with respect to the secondary angle of twist is obtained in terms of the primary angle of twist as
(98)d2θxSdx2=-EC-SGItPd4θxPdx4(1-κ)insidethebar
while the corresponding boundary conditions are given in terms of the primary angle of twist as
(99)γ1θxS+γ2(GItPκ1-κ)dθxSdx=γ3(-EC-Sκd3θxPdx3)atthebarendsx=0,l,
where if the primary angle of twist is fixed, that is a1=1 and a2=a3=0, then γ1=1 and γ2=γ3=0, otherwise γ1=0 and γ2=γ3=1.From the examined example problems presented in [85] it is concluded that (i)
the inaccuracy of the thin tube theory in calculating the secondary torsion constant even for thin-walled sections is noteworthy;(ii)
as it is also verified by other computational methods [80, 83], for bars with open shaped cross section the secondary torsional moment deformation effect has not significant influence. However, the inclusion of this additional effect leads to more accurate results;(iii)
as it is also verified by other computational methods [80, 83], for bars with closed shaped cross sections the secondary torsional moment deformation effect has an important influence and should be considered in the analysis (otherwise the results may be completely wrong especially in the calculation of stresses).


## 7. Nonlinear Elastic Nonuniform Torsion of Bars

As it has already been mentioned in the previous sections in engineering practice we often come across the analysis of members of structures subjected to twisting moments. Besides, since thin-walled open sections have low torsional stiffness, the torsional deformations can be of such magnitudes that it is not adequate to treat the angles of cross section rotation as small. When finite twist rotation angles are considered, the elastic nonuniform torsion problem becomes nonlinear. Moreover, this problem becomes much more complicated in the case the cross section’s centroid does not coincide with its shear center (asymmetric beams), leading to the formulation of a flexural-torsional coupled problem. The extensive use of the aforementioned structural elements necessitates a reliable and accurate analysis of bars of arbitrary cross section subjected to torsional loading taking into account the geometrical nonlinearity.Though several researchers have dealt either with the linear nonuniform torsional behaviour of beams [17, 21, 22, 40] or with the nonlinear uniform torsional behaviour of doubly symmetric beams [89, 90], to the author’s knowledge very little work has been done on the corresponding nonlinear nonuniform torsional problem of arbitrary cross section beams. Ghobarah and Tso [91], Attard [92], and Attard and Somervaille [93] have presented a set of displacement relationships for a straight prismatic thin-walled open beam applicable to situations where displacements are finite, the cross section does not distort, strains are small and flexural displacements are small to moderate while cross sectional twist can be large. The presented numerical examples in these studies are concerned only with uniform torsion of either mono- or doubly symmetric cross sections. Finally, Trahair in [94] employing the finite element method and presenting examples of only doubly symmetric cross sections and Mohri et al. in [95] employing similar equations to those established by Attard in [92] and presenting examples of either doubly symmetric cross sections subjected in nonuniform torsion or buckling or postbuckling behavior of arbitrary cross section beams also analyze the nonlinear nonuniform torsional problem. Nevertheless, all of the aforementioned studies, which are the only one considering finite angles of twist in asymmetrical bars (and taking into account all of the arising nonlinear terms) are not general since they are restricted to thin-walled beams. Finally, Sapountzakis and Tsipiras in [96, 97] presented a BE solution for the elastic nonuniform torsion analysis of simply, multiply connected or composite cylindrical bars of arbitrary cross section taking into account the effect of geometric nonlinearity. The torque-rotation relationship is computed based on the finite displacement (finite rotation) theory; that is, the transverse displacement components are expressed so as to be valid for large rotations and the longitudinal normal strain includes the second-order geometrically nonlinear term often described as the “Wagner strain.” These last formulations do not stand on the assumption of a thin-walled structure and therefore the cross section’s torsional rigidity is evaluated exactly without using the so-called Saint-Venant’s torsional constant. The torsional rigidity of the cross section is evaluated directly employing the primary warping function of the cross section [21] depending on its shape. Three boundary value problems with respect to the variable along the beam axis angle of twist, to the primary and to the secondary warping functions are formulated. The first one of these problems is numerically solved employing the Analog Equation Method [98], a BEM based method, leading to a system of nonlinear equations from which the angle of twist is computed by an iterative process. The other two problems are solved using a pure BEM [23] based method. The aforementioned formulation procedure is based on the assumption of no local or lateral torsional buckling or distortion.In order to formulate the aforementioned problem, let us consider a prismatic bar of lengthl (Figure 13), of constant arbitrary cross-section of area A. The homogeneous isotropic and linearly elastic material of the bar’s cross-section, with modulus of elasticity E, shear modulus G, and Poisson’s ratio ν, occupies the two dimensional multiply connected region Ω of the y, z plane and is bounded by the Γj(j=0,1,2,…,K) boundary curves, which are piecewise smooth; that is, they may have a finite number of corners. In Figure 13(b)CYZ is the principal coordinate system through the cross section’s centroid C, while yC, zC are its coordinates with respect to Syz system of axes through the cross section’s shear center S. The bar is subjected to the combined action of the arbitrarily distributed or concentrated conservative twisting mt=mt(x) and warping mw=mw(x) moments acting in the x direction (Figure 13(a)).Prismatic bar subjected to twisting and warping moments (a) with a cross section of arbitrary shape occupying the two dimensional regionΩ(b).
(a)
(b)Under the aforementioned loading the displacement field of the bar with respect to theSyz system of axes for large twisting rotations and small bending ones is given as(100a)u=um(x)+θY(x)·(z-zC)-θZ(x)·(y-yC)+θx′(x)·ϕSP(y,z)+ϕSS(x,y,z),(100b)v=vS(x)-z·sinθx(x)-y·(1-cosθx(x)),(100c)w=wS(x)+y·sinθx(x)-z·(1-cosθx(x)),where the transverse displacement components v, w are valid for large rotations [99]; θY, θZ are the angles of rotation due to bending of the cross section with respect to its centroid; θx′(x) denotes the rate of change of the angle of twist θx regarded as the torsional curvature; vS(x), wS(x) are the transverse displacement components of the shear center S; ϕSP, ϕSS are the primary and secondary warping functions with respect to the shear center S, respectively [21]; um(x) is an “average” axial displacement of the cross section of the bar, that will be later explained.Substituting (100a), (100b), and (100c) in the nonlinear (Green) strain-displacement relations of the nonvanishing strains(101a)εxx=∂u∂x+12·[(∂u∂x)2+(∂v∂x)2+(∂w∂x)2],(101b)γxy=∂u∂y+∂v∂x+∂u∂y·∂u∂x+∂v∂y·∂v∂x+∂w∂y·∂w∂x,(101c)γxz=∂w∂x+∂u∂z+∂u∂z·∂u∂x+∂v∂z·∂v∂x+∂w∂z·∂w∂xassuming moderate large deflections ((∂u/∂x)2≪∂u/∂x, (∂u/∂x)(∂u/∂y)≪(∂v/∂x)+(∂u/∂y), (∂u/∂x)(∂u/∂z)≪(∂w/∂x)+(∂u/∂z)) and having in mind that ignoring shear deformation effect the angles of rotation due to bending θY, θZ approximate the slope between the longitudinal axis in the deformed state and in the initial vertical and horizontal planes, respectively, the following relations are obtained:(102a)θY(x)=vS′·sinθx-wS′·cosθx,(102b)θZ(x)=vS′·cosθx+wS′·sinθxand the nonvanishing strain resultants are given as(103a)εxx=um′+κY·(z-zC)-κZ·(y-yC)+θx′′·ϕSP(y,z)-θx′·(yC·θY+zC·θZ)+12[(vS′)2+(wS′)2+(y2+z2)·(θx′)2],(103b)γxy=θx′·(∂ϕSP∂y-z)+∂ϕSS∂y,(103c)γxz=θx′·(∂ϕSP∂z+y)+∂ϕSS∂z,where the curvature components κY, κZ are given from the following relations:(104a)κY(x)=vS′′(x)·sinθx-wS′′(x)·cosθx,(104b)κZ(x)=vS′′(x)·cosθx+wS′′(x)·sinθxwhile the second-order geometrically nonlinear term in the right hand side of (103a) (y2+z2)·(θx′)2/2 is often described as the “Wagner strain” [94]. It is worth here noting that in obtaining (103a) the rate of change of the secondary warping function ϕSS, that is the arising normal stress due to the secondary shear one due to warping [24], has been ignored.Considering strains to be small, employing the second Piola-Kirchhoff stress tensor and assuming an isotropic and homogeneous material for zero Poisson ratio, the stress components are defined in terms of the strain ones as(105){SxxSxySxz}=[E000G000G]{εxxγxyγxz}
or employing (103a), (103b), and (103c) as(106a)Sxx=E[um′+κY·(z-zC)-κZ·(y-yC)+θx′′·ϕSP(y,z)-θx′·(yCθY+zCθZ)]+E2[(vS′)2+(wS′)2+(y2+z2)·(θx′)2],(106b)Sxy=G·θx′·(∂ϕSP∂y-z)+G·∂ϕSS∂y,(106c)Sxz=G·θx′·(∂ϕSP∂z+y)+G·∂ϕSS∂z.


### 7.1. Equations of Local Equilibrium

Applying the stress components (106a), (106b), and (106c) in the first equation of equilibrium neglecting the body forces [100]
(107)∂Sxx∂x+∂Sxy∂y+∂Sxz∂z=0
ignoring the independent of warping terms such as E·um′′ or warping terms due to shear such as E(wS′′′(z-zC)+vS′′′(y-yC)) and requiring both the primary and the secondary due to warping parts of both (107) and the traction vector on the free surface of the bar to vanish, the Neumann problems for the primary warping function ϕSP given from (19a) and (19b) and for the secondary warping function ϕSS given from(108a)∇2ϕSS=-EG·θx′′′·ϕSP(y,z)-EG·(y2+z2)·θx′·θx′′inΩ,(108b)∂ϕSS∂n=0onΓare obtained, where ∇2=∂2/∂y2+∂2/∂z2 is the Laplace operator and ∂/∂n denotes the directional derivative normal to the boundary Γ. Thus, the primary ϕSP and the secondary ϕSS warping functions will be evaluated from the solution of the Neumann problems described by the governing equations (19a) and (108a) inside the two dimensional multiply connected region Ω, subjected to the boundary conditions (19b) and (108b) on its boundary Γ, respectively. It is worth here noting that the evaluated warping functions due to the solution of the corresponding Neumann problems contain an integration constant (parallel displacement of the cross section along the bar axis), which can be obtained from the requests that
(109)∫ΩϕSPdΩ=0,(110)∫ΩϕSSdΩ=0.
Moreover, in the case the origin O of the coordinates is a point of the plane other than the shear center S, the primary warping function with respect to this point ϕOP is first established from the Neumann problems (19a) and (19b) substituting ϕSP by ϕOP. Using the evaluated warping function ϕOP, ϕSP is then established using the transformation presented in Section 2, while the coordinates of the shear center S with respect to the Oy¯z¯ system of coordinates are established from (109) and the conditions(111a)∫ΩϕSPydΩ=0,(111b)∫ΩϕSPzdΩ=0that define point S as the point for which the axial and bending stress resultants arising from an angle of twist at this point vanish in a geometrically linear analysis. It is worth here noting that the shear center axis of an asymmetric cross section bar subjected in pure torsion is not a rectilinear one but is assigned from the displacement field uS(x), vS(x), wS(x).Based on the conditions (109) and (110) the meaning of the “average” axial displacement of the cross section of the bar can now be explained as the relation
(112)∫ΩudΩ=um(x)·A+θY(x)·(∫ΩzdΩ-zC·A)-θZ(x)·(∫ΩydΩ-yC·A)+θx′·∫ΩϕSPdΩ+∫ΩϕSSdΩ
leads to
(113)um(x)=∫ΩudΩA
that is um generally does not coincide either with the axial displacement of the cross section’s centroid uC or with that of the shear center uS.


### 7.2. Equations of Global Equilibrium

Neglecting body forces the principle of virtual work yields(114)∫V(Sxx·δεxx+Sxy·δγxy+Sxz·δγxz)dV=∫A(tx·δu+ty·δv+tz·δw)dA,
where δ(·) denotes virtual quantities, tx, ty, tz are the components of the traction vector with respect to the undeformed surface of the bar including the cross section ends, A is the surface area, and V is the volume of the bar at its initial undeformed state. Moreover, the stress resultants of the bar are given as(115a)N=∫ΩSxxdΩ,(115b)MY=∫ΩSxx(z-zC)dΩ,(115c)MZ=-∫ΩSxx(y-yC)dΩ,(115d)MtP=∫Ω[SxyP·(∂ϕSP∂y-z)+SxzP·(∂ϕSP∂z+y)]dΩ,(115e)Mw=-∫ΩSxxϕSPdΩ,where MtP is the primary twisting moment resulting from the primary shear stress distribution SxyP, SxzP and Mw is the warping moment due to the torsional curvature. Substituting (106a), (106b), and (106c) into (115a), (115b), (115c), (115d), and (115e) the stress resultants are obtained as(116a)N=EA·[um′+12·((vS′)2+(wS′)2+IPA·(θx′)2)-θx′·(yCθY+zCθZ)((vS′)2+(wS′)2+IPA·(θ′x)2)],(116b)MY=EIYY·[κY+β2(θx′)2],(116c)MZ=EIZZ·[κZ-β1(θx′)2],(116d)MtP=GIt·θx′,(116e)Mw=-ECS·[θx′′+Uw2CS·(θx′)2],where the torsion constant It and the warping constant CS with respect to the shear center S are given from (32a) and (32b), the polar moment of inertia IP with respect to the shear center S and the moments of inertia IYY, IZZ with respect to the cross section’s centroid are given as(117a)IP=∫Ωr2dΩ,(117b)IYY=∫Ω(z-zC)2dΩ,(117c)IZZ=∫Ω(y-yC)2dΩwith r=y2+z2, while the geometric constants β1, β2, and Uw are given as(118a)β1=12IZZ·∫Ωr2·(y-yC)dΩ,(118b)β2=12IYY·∫Ωr2·(z-zC)dΩ,(118c)Uw=∫ΩϕSP·r2dΩ.It is worth here noting that the aforementioned stress resultants refer to the directions of the cross section at its deformed configuration (they take into account the cross sections’ rotations), since they have been defined with respect to the second Piola-Kirchhoff stress tensor. Moreover, though the direction of the axial forceN is the tangential one to the deformed centroid axis, it is not necessarily applied to the cross section’s centroid, as it is already mentioned the “average” axial displacement um is not necessarily the centroid’s displacement uC [92]. In deriving the relations (116a), (116b), (116c), (116d), and (116e), the properties of the principal system of axes and the orthogonality conditions of the primary warping function ϕSP with respect to the shear center S have been employed.Substituting the stress components given in (106a), (106b), and (106c) and the strain resultants given in (103a), (103b), and (103c) to the principle of virtual work (114) the equation of torque equilibrium of the bar is obtained as
(119)-NyCθZθx′+NzCθYθx′-MZκY+MYκZ-ddx[MtP+12EIn(θx′)3+Ψθx′-NyCθY-NzCθZ]-d2Mwdx2=mt(x)+ddx[mw(x)]
and the corresponding boundary conditions are written as(120a)[(Mw)′+MtP+12EIn(θx′)3+Ψθx′-NyCθY-NzCθZ-Mt¯{12EIn(θx′)3}(Mw)′12]δθx=0,(120b)(-Mw+Mw¯)δθx′=0,where the “higher order torsion constant” In and the variable Ψ are given as(121a)In=IPP-IP2A-4β12IZZ-4β22IYY-Uw2CS,(121b)Ψ=NIPA-2MZβ1+2MYβ2-MwUwCSwhile Mt¯, Mw¯ in the boundary conditions (120a) and (120b) are the externally applied conservative twisting and warping moments at the bar ends and IPP=∫Ωr4dΩ (in (121a)). It is worth here noting that in deriving the equation of torque equilibrium of the bar the secondary shear stress distribution has been ignored [24]. Considering that the axial force and the bending moments vanish along the bar(122a)N=0,(122b)MY=0,(122c)MZ=0the governing equation of the nonuniform nonlinear torsional problem is written as
(123)ECSθx′′′′-32EIn2(θx′)2θx′′-GItθx′′=mt(x)+ddx[mw(x)]
subjected to the following boundary conditions at the bar ends(124a)α1θx+α2Mt=α3,(124b)β1θx′+β2Mw=β3,where In2 is defined as
(125)In2=In+Uw2CS
and Mt, Mw are the twisting and warping moments at the boundary of the bar, respectively, given as(126a)Mt=-ECSθx′′′+GItθx′+12EIn2(θx′)3,(126b)Mw=-ECS[θx′′+Uw2CS(θx′)2]while αi, βi(i=1,2,3) are functions specified at the boundary of the bar. The boundary conditions (124a) and (124b) are the most general linear torsional boundary conditions for the beam problem including also the elastic support. It is apparent that all types of the conventional torsional boundary conditions (clamped, simply supported, free, or guided edge) can be derived from these equations by specifying appropriately the functions αi and βi (e.g., for a clamped edge it is α1=β1=1, α2=α3=β2=β3=0).The solution of the boundary value problem described by (123), (124a), and (124b) for the evaluation of the angle of twist θx assumes that the warping CS, the torsion It, and the geometric Uw constants defined from (32a), (32b), and (118c) are already established. Equations (32a), (32b), and (118c) indicate that the evaluation of the aforementioned constants presumes that the primary warping function ϕSP at any interior point of the domain Ω of the cross section of the bar is evaluated after solving the boundary value problem described by (19a) and (19b). Once θx is established, the secondary warping function ϕSS at any interior point of the domain Ω of the cross section of the bar is evaluated after solving the boundary value problem described by (108a) and (108b). Subsequently the second Piola-Kirchhoff stress components are evaluated employing (106a), (106b), and (106c). Moreover, using the evaluated angle of twist θx, the transverse displacement components of the shear center vS(x), wS(x) can be established integrating twice the following relations: (127a)vS′′=-(θx′)2·(β2·sinθx-β1·cosθx),(127b)wS′′=(θx′)2·(β2·cosθx+β1·sinθx)obtained from the solution of the linear system of equations arising from the substitution of (116b) and (116c) to (122b) and (122c) after employing (104a) and (104b). Finally, using the evaluated derivatives of the transverse displacement components vS′, wS′ the “average” axial displacement of the cross section of the bar um can be established after substitution of (116a) to (122a) as
(128)um′=-12·[(vS′)2+(wS′)2+IPA·(θx′)2]+θx′·(yCθY+zCθZ)
and subsequent integration.From the examined example problems presented in [96, 97] it is concluded that(i)
geometrical nonlinearity leads to stiffening of the bar and eventually its better response against torsional loading;(ii)
as observed, restrained warping boundary conditions lead to slightly greater torsional resistance than those of free warping for both the geometrically linear and nonlinear cases;(iii)
axial shortening of bars of doubly symmetric cross section subjected to torsional loading is observed in the nonlinear analysis, while bars of asymmetrical cross section exhibit lateral displacement of their shear center axis as well;(iv)
the nonlinear Wagner torque can reach significant values locally along the bar axis while the nonlinear warping moment can be neglected in most cases;(v)
the secondary shear stress distribution and magnitude are influenced by the geometrical nonlinearity.


## 7.1. Equations of Local Equilibrium

Applying the stress components (106a), (106b), and (106c) in the first equation of equilibrium neglecting the body forces [100]
(107)∂Sxx∂x+∂Sxy∂y+∂Sxz∂z=0
ignoring the independent of warping terms such as E·um′′ or warping terms due to shear such as E(wS′′′(z-zC)+vS′′′(y-yC)) and requiring both the primary and the secondary due to warping parts of both (107) and the traction vector on the free surface of the bar to vanish, the Neumann problems for the primary warping function ϕSP given from (19a) and (19b) and for the secondary warping function ϕSS given from(108a)∇2ϕSS=-EG·θx′′′·ϕSP(y,z)-EG·(y2+z2)·θx′·θx′′inΩ,(108b)∂ϕSS∂n=0onΓare obtained, where ∇2=∂2/∂y2+∂2/∂z2 is the Laplace operator and ∂/∂n denotes the directional derivative normal to the boundary Γ. Thus, the primary ϕSP and the secondary ϕSS warping functions will be evaluated from the solution of the Neumann problems described by the governing equations (19a) and (108a) inside the two dimensional multiply connected region Ω, subjected to the boundary conditions (19b) and (108b) on its boundary Γ, respectively. It is worth here noting that the evaluated warping functions due to the solution of the corresponding Neumann problems contain an integration constant (parallel displacement of the cross section along the bar axis), which can be obtained from the requests that
(109)∫ΩϕSPdΩ=0,(110)∫ΩϕSSdΩ=0.
Moreover, in the case the origin O of the coordinates is a point of the plane other than the shear center S, the primary warping function with respect to this point ϕOP is first established from the Neumann problems (19a) and (19b) substituting ϕSP by ϕOP. Using the evaluated warping function ϕOP, ϕSP is then established using the transformation presented in Section 2, while the coordinates of the shear center S with respect to the Oy¯z¯ system of coordinates are established from (109) and the conditions(111a)∫ΩϕSPydΩ=0,(111b)∫ΩϕSPzdΩ=0that define point S as the point for which the axial and bending stress resultants arising from an angle of twist at this point vanish in a geometrically linear analysis. It is worth here noting that the shear center axis of an asymmetric cross section bar subjected in pure torsion is not a rectilinear one but is assigned from the displacement field uS(x), vS(x), wS(x).Based on the conditions (109) and (110) the meaning of the “average” axial displacement of the cross section of the bar can now be explained as the relation
(112)∫ΩudΩ=um(x)·A+θY(x)·(∫ΩzdΩ-zC·A)-θZ(x)·(∫ΩydΩ-yC·A)+θx′·∫ΩϕSPdΩ+∫ΩϕSSdΩ
leads to
(113)um(x)=∫ΩudΩA
that is um generally does not coincide either with the axial displacement of the cross section’s centroid uC or with that of the shear center uS.


## 7.2. Equations of Global Equilibrium

Neglecting body forces the principle of virtual work yields(114)∫V(Sxx·δεxx+Sxy·δγxy+Sxz·δγxz)dV=∫A(tx·δu+ty·δv+tz·δw)dA,
where δ(·) denotes virtual quantities, tx, ty, tz are the components of the traction vector with respect to the undeformed surface of the bar including the cross section ends, A is the surface area, and V is the volume of the bar at its initial undeformed state. Moreover, the stress resultants of the bar are given as(115a)N=∫ΩSxxdΩ,(115b)MY=∫ΩSxx(z-zC)dΩ,(115c)MZ=-∫ΩSxx(y-yC)dΩ,(115d)MtP=∫Ω[SxyP·(∂ϕSP∂y-z)+SxzP·(∂ϕSP∂z+y)]dΩ,(115e)Mw=-∫ΩSxxϕSPdΩ,where MtP is the primary twisting moment resulting from the primary shear stress distribution SxyP, SxzP and Mw is the warping moment due to the torsional curvature. Substituting (106a), (106b), and (106c) into (115a), (115b), (115c), (115d), and (115e) the stress resultants are obtained as(116a)N=EA·[um′+12·((vS′)2+(wS′)2+IPA·(θx′)2)-θx′·(yCθY+zCθZ)((vS′)2+(wS′)2+IPA·(θ′x)2)],(116b)MY=EIYY·[κY+β2(θx′)2],(116c)MZ=EIZZ·[κZ-β1(θx′)2],(116d)MtP=GIt·θx′,(116e)Mw=-ECS·[θx′′+Uw2CS·(θx′)2],where the torsion constant It and the warping constant CS with respect to the shear center S are given from (32a) and (32b), the polar moment of inertia IP with respect to the shear center S and the moments of inertia IYY, IZZ with respect to the cross section’s centroid are given as(117a)IP=∫Ωr2dΩ,(117b)IYY=∫Ω(z-zC)2dΩ,(117c)IZZ=∫Ω(y-yC)2dΩwith r=y2+z2, while the geometric constants β1, β2, and Uw are given as(118a)β1=12IZZ·∫Ωr2·(y-yC)dΩ,(118b)β2=12IYY·∫Ωr2·(z-zC)dΩ,(118c)Uw=∫ΩϕSP·r2dΩ.It is worth here noting that the aforementioned stress resultants refer to the directions of the cross section at its deformed configuration (they take into account the cross sections’ rotations), since they have been defined with respect to the second Piola-Kirchhoff stress tensor. Moreover, though the direction of the axial forceN is the tangential one to the deformed centroid axis, it is not necessarily applied to the cross section’s centroid, as it is already mentioned the “average” axial displacement um is not necessarily the centroid’s displacement uC [92]. In deriving the relations (116a), (116b), (116c), (116d), and (116e), the properties of the principal system of axes and the orthogonality conditions of the primary warping function ϕSP with respect to the shear center S have been employed.Substituting the stress components given in (106a), (106b), and (106c) and the strain resultants given in (103a), (103b), and (103c) to the principle of virtual work (114) the equation of torque equilibrium of the bar is obtained as
(119)-NyCθZθx′+NzCθYθx′-MZκY+MYκZ-ddx[MtP+12EIn(θx′)3+Ψθx′-NyCθY-NzCθZ]-d2Mwdx2=mt(x)+ddx[mw(x)]
and the corresponding boundary conditions are written as(120a)[(Mw)′+MtP+12EIn(θx′)3+Ψθx′-NyCθY-NzCθZ-Mt¯{12EIn(θx′)3}(Mw)′12]δθx=0,(120b)(-Mw+Mw¯)δθx′=0,where the “higher order torsion constant” In and the variable Ψ are given as(121a)In=IPP-IP2A-4β12IZZ-4β22IYY-Uw2CS,(121b)Ψ=NIPA-2MZβ1+2MYβ2-MwUwCSwhile Mt¯, Mw¯ in the boundary conditions (120a) and (120b) are the externally applied conservative twisting and warping moments at the bar ends and IPP=∫Ωr4dΩ (in (121a)). It is worth here noting that in deriving the equation of torque equilibrium of the bar the secondary shear stress distribution has been ignored [24]. Considering that the axial force and the bending moments vanish along the bar(122a)N=0,(122b)MY=0,(122c)MZ=0the governing equation of the nonuniform nonlinear torsional problem is written as
(123)ECSθx′′′′-32EIn2(θx′)2θx′′-GItθx′′=mt(x)+ddx[mw(x)]
subjected to the following boundary conditions at the bar ends(124a)α1θx+α2Mt=α3,(124b)β1θx′+β2Mw=β3,where In2 is defined as
(125)In2=In+Uw2CS
and Mt, Mw are the twisting and warping moments at the boundary of the bar, respectively, given as(126a)Mt=-ECSθx′′′+GItθx′+12EIn2(θx′)3,(126b)Mw=-ECS[θx′′+Uw2CS(θx′)2]while αi, βi(i=1,2,3) are functions specified at the boundary of the bar. The boundary conditions (124a) and (124b) are the most general linear torsional boundary conditions for the beam problem including also the elastic support. It is apparent that all types of the conventional torsional boundary conditions (clamped, simply supported, free, or guided edge) can be derived from these equations by specifying appropriately the functions αi and βi (e.g., for a clamped edge it is α1=β1=1, α2=α3=β2=β3=0).The solution of the boundary value problem described by (123), (124a), and (124b) for the evaluation of the angle of twist θx assumes that the warping CS, the torsion It, and the geometric Uw constants defined from (32a), (32b), and (118c) are already established. Equations (32a), (32b), and (118c) indicate that the evaluation of the aforementioned constants presumes that the primary warping function ϕSP at any interior point of the domain Ω of the cross section of the bar is evaluated after solving the boundary value problem described by (19a) and (19b). Once θx is established, the secondary warping function ϕSS at any interior point of the domain Ω of the cross section of the bar is evaluated after solving the boundary value problem described by (108a) and (108b). Subsequently the second Piola-Kirchhoff stress components are evaluated employing (106a), (106b), and (106c). Moreover, using the evaluated angle of twist θx, the transverse displacement components of the shear center vS(x), wS(x) can be established integrating twice the following relations: (127a)vS′′=-(θx′)2·(β2·sinθx-β1·cosθx),(127b)wS′′=(θx′)2·(β2·cosθx+β1·sinθx)obtained from the solution of the linear system of equations arising from the substitution of (116b) and (116c) to (122b) and (122c) after employing (104a) and (104b). Finally, using the evaluated derivatives of the transverse displacement components vS′, wS′ the “average” axial displacement of the cross section of the bar um can be established after substitution of (116a) to (122a) as
(128)um′=-12·[(vS′)2+(wS′)2+IPA·(θx′)2]+θx′·(yCθY+zCθZ)
and subsequent integration.From the examined example problems presented in [96, 97] it is concluded that(i)
geometrical nonlinearity leads to stiffening of the bar and eventually its better response against torsional loading;(ii)
as observed, restrained warping boundary conditions lead to slightly greater torsional resistance than those of free warping for both the geometrically linear and nonlinear cases;(iii)
axial shortening of bars of doubly symmetric cross section subjected to torsional loading is observed in the nonlinear analysis, while bars of asymmetrical cross section exhibit lateral displacement of their shear center axis as well;(iv)
the nonlinear Wagner torque can reach significant values locally along the bar axis while the nonlinear warping moment can be neglected in most cases;(v)
the secondary shear stress distribution and magnitude are influenced by the geometrical nonlinearity.


## 8. Nonlinear Elastic Torsional Vibrations of Bars

As it has already been mentioned in the previous section, since weight saving is of paramount importance, frequently used thin-walled open sections have low torsional stiffness and their torsional deformations can be of such magnitudes that it is not adequate to treat the angles of cross section rotation as small. In these cases, the study of nonlinear effects on these members becomes essential, where this nonlinearity results from retaining the nonlinear terms in the strain-displacement relations (finite displacement—small strain theory). When finite twist rotation angles are considered, the nonuniform torsional dynamic analysis of bars becomes much more complicated, leading to the formulation of coupled and nonlinear torsional and axial equilibrium equations. In this case, accounting for the axial loading and boundary conditions becomes essential to perform a rigorous dynamic analysis of the bar.During the past few years, the nonlinear nonuniform torsional dynamic analysis of bars undergoing large rotations has received a good amount of attention in the literature. More specifically, Rozmarynowski and Szymczak in [101] studied the nonlinear free torsional vibrations of axially immovable thin-walled beams with doubly symmetric open cross section, employing the finite element method. In this research effort only free vibrations are examined, the solution is provided only at points of reversal of motion (not in the time domain), no general axial, torsional, or warping boundary conditions (elastic support case) are studied, while some nonlinear terms related to the finite twisting rotations as well as the axial inertia term are ignored. Da Silva in [102, 103] presented the nonlinear flexural-torsional-extensional vibrations of Euler-Bernoulli doubly symmetric thin-walled closed cross section beams, primarily focusing on flexural vibrations and neglecting the effect of torsional warping. Pai and Nayfeh in [104–106] studied also the nonlinear flexural-torsional-extensional vibrations of metallic and composite slewing or rotating closed cross section beams, primarily focusing on flexural vibrations and neglecting again the effect of torsional warping. Di Egidio et al. in [107, 108] presented also a FEM solution to the nonlinear flexural-torsional vibrations of shear undeformable thin-walled open beams taking into account in-plane and out-of-plane warpings and neglecting warping inertia. In these papers, the torsional-extensional coupling is taken into account but the inextensionality assumption leads to the fact that the axial boundary conditions are not general. Simo and Vu-Quoc in [109] presented a FEM solution to a fully nonlinear (small or large strains, hyperelastic material) three dimensional rod model including the effects of transverse shear and torsion-warping deformation based on a geometrically exact description of the kinematics of deformation. Moreover, Pai and Nayfeh in [110] studied a geometrically exact nonlinear curved beam model for solid composite rotor blades using the concept of local engineering stress and strain measures and taking into account the in-plane and out-of-plane warpings. In all of the aforementioned research efforts the evaluation of the warping shear stress distribution is not discussed. Finally, Sapountzakis and Tsipiras in [111, 112] developed a boundary element solution for the nonuniform torsional vibration problem of bars of arbitrary doubly symmetric constant cross section taking into account the effect of geometrical nonlinearity. In these last research efforts the bar is subjected to arbitrarily distributed or concentrated conservative dynamic twisting and warping moments along its length, while its edges are supported by the most general torsional and axial boundary conditions. A distributed mass model system is employed, taking into account the warping, rotatory, and axial inertia, which after employing a variational approach leads to the formulation of a coupled nonlinear initial boundary value problem with respect to the variable along the bar angle of twist and to an “average” axial displacement of the cross section of the bar and to two boundary value problems with respect to the primary and secondary warping functions. The numerical solution of the aforementioned initial boundary value problem is performed using the Analog Equation Method [98], a BEM based method, leading to a system of nonlinear Differential-Algebraic Equations (DAEs), which is solved using an efficient time discretization scheme. The other two problems are solved using a pure BE method, requiring exclusively boundary discretization of the bar’s cross-section. Both free and forced vibrations are considered, while the arising linear system of equations related to the secondary warping function is singular and a special technique [113] is used to perform its regularization.In order to formulate the nonlinear torsional vibration problem of doubly symmetric bars of arbitrarily shaped simply or multiply connected cross section, let us consider the bar of lengthl of Figure 14. The homogeneous isotropic and linearly elastic material of the bar’s cross-section of area A, with modulus of elasticity E, shear modulus G, and Poisson’s ratio ν, occupies the two dimensional multiply connected region Ω of the y, z plane and is bounded by the Γj(j=0,1,2,…,K) boundary curves, which are piecewise smooth; that is, they may have a finite number of corners. In Figure 14(b)Syz is the principal bending coordinate system through the cross section’s shear center. The bar is subjected to the combined action of the arbitrarily distributed or concentrated time dependent conservative axial load n(x,t), t≥0, twisting mt=mt(x,t) and warping mw=mw(x,t) moments, acting in the x direction (Figure 14(a)).Bar subjected to axial and torsional time dependent loading (a) of arbitrarily shaped doubly symmetric constant cross section occupying regionΩ (b).
(a)
(b)For the analysis of the aforementioned problem, the same displacement field (100a), (100b), and (100c) with that of the nonlinear static analysis is adopted, which for the dynamic problem and for a doubly symmetric cross section is written as(129a)u(x,y,z,t)=um(x,t)+θx′(x,t)·ϕSP(y,z)+ϕSS(x,y,z,t),(129b)v(x,y,z,t)=-z·sinθx(x,t)-y·(1-cosθx(x,t)),(129c)w(x,y,z,t)=y·sinθx(x,t)-z·(1-cosθx(x,t)).


### 8.1. Equations of Local Equilibrium

Following the same procedure presented in Section7, the local equilibrium equations are established employing the principle of virtual work
(130)δWint+δWmass=δWext,
where(131a)δWint=∫V(Sxxδεxx+Sxyδγxy+Sxzδγxz)dV,(131b)δWmass=∫Vρ·(u¨δu+v¨δv+w¨δw)dV,(131c)δWext=∫F(tx·δu+ty·δv+tz·δw)dFunder a Total Lagrangian formulation. In the above equations, δ(·) denotes virtual quantities, (·) denotes differentiation with respect to time, V, F are the volume and the surface of the bar, respectively, at the initial configuration, and tx, ty, tz are the components of the traction vector with respect to the undeformed surface of the bar.Neglecting virtual terms of the secondary warping functionϕSS, following the technique presented in [100] and after some algebra [112], the following local equilibrium equation
(132)∂Sxx∂x+∂Sxy∂y+∂Sxz∂z-ρu¨=0inΩ,∀x∈(0,l)
along with its corresponding boundary condition
(133)Sxyny+Sxznz=0onΓj,∀x∈(0,l)
is obtained. Requiring both the primary and the secondary due to warping parts of (132) and (133) to vanish, employing (106a), (106b), and (106c) after taking into account the double symmetry of the cross section and ignoring the term ϕ¨SS, the following governing equation for the secondary ϕSS warping function is obtained as
(134)∇2ϕSS=(-EGum′′+ρGu¨m)-(EGθx′′′-ρGθ¨x′)·ϕSP(y,z)-EGθx′θx′′·(y2+z2)inΩ
along with its corresponding boundary condition
(135)∂ϕSS∂n=0onΓj(j=1,2,…,K)
while as expected, the governing equation related to the primary ϕSP warping function is found to be coincident with the well-known St. Venant corresponding boundary value problem given from (19a) and (19b). As it has already been mentioned for the static case, the evaluation of the warping function contains an integration constant (parallel displacement of the cross section along the bar axis), which can be obtained from (110), arising from the request that the torsional terms of the displacement field do not result in any axial forces. It is worth pointing out that any other constraints could be used, although the use of (110) decouples the governing equations of the torsional problem at the greatest extent.


### 8.2. Equations of Global Equilibrium

Taking into account the double symmetry of the cross section and substituting the stress components ((106a), (106b), and (106c)) the strain ones ((103a), (103b), and (103c)) and the displacement components ((129a), (129b), and (129c)) to the principle of virtual work (130), the governing partial differential equations of the bar are obtained after some algebra as(136a)ρA·u¨m-EA·um′′-EIP·θx′θx′′=n(x,t),(136b)ρIP·θ¨x-ρCS·θ¨x′′-GItθx′′+ECSθx′′′′-32EIPP·(θx′)2θx′′-EIP·um′θx′′-EIP·um′′θx′=mt(x,t)+∂∂x[mw(x,t)]subjected to the initial conditions (x∈(0,l))(137a)um(x,0)=u-m0(x),(137b)u˙m(x,0)=u-˙m0(x),(137c)θx(x,0)=θ-x0(x),(137d)θ˙x(x,0)=θ-˙x0(x)together with the boundary conditions at the bar ends x=0,l(138a)a1N+α2um=α3,(138b)β1Mt+β2θx=β3,(138c)β-1Mw+β-2θx′=β-3,where N, Mt, Mw are the axial force, the twisting and warping moments at the bar ends given from (116a), (126a), and (126b), respectively. The expressions of the externally applied loads appearing in the right hand side of (136a) and (136b) with respect to the first Piola-Kirchhoff stress components can be easily deduced by virtue of (131c). It is worth here noting that damping could also be included in the analysis without any special difficulty, while the geometric cross sectional properties appearing in (136a) and (136b) are also given in Section 7.The solution of the initial boundary value problem described by (136a)–(138c) for the evaluation of the unknown kinematical components um(x,t), θx(x,t) assumes that the warping CS and the torsion It constants are already established; the evaluation of which presumes that the primary warping function ϕSP at any interior point of the domain Ω of the cross section of the bar is evaluated. Once these components are established, the secondary warping function ϕSS at any interior point of the domain Ω of the cross section of the bar is evaluated after solving the boundary value problem described by (134) and (135). Subsequently the second Piola-Kirchhoff stress components are evaluated employing (106a), (106b), and (106c) completing the computation of the stress field.A significant reduction on both the set of the governing differential equations and the boundary value problem related to the secondary warping functionϕSS can be achieved by neglecting the axial inertia term ρA·u¨m of (136a), an assumption which is common among various dynamic beam formulations. Ignoring this term, a single partial differential equation along with a single unknown kinematical component (the angle of twist θx(x,t)) is obtained, which is further simplified in the case of vanishing distributed axial load along the bar. In what follows, this procedure is described in detail for the cases of axially immovable ends and constant axial load along the bar, which are of great practical interest.


### 8.3. Reduced Problems for Special Cases of Axial Boundary Conditions

For the case of axially immovable ends, it is easily proved that [112]
(139)um′′=-IPA·θx′θx′′,∀x∈(0,l)
which after subsequent integration yields
(140)um′=-12·IPA·(θx′)2+N~EA,∀x∈[0,l],
where N~ is a time-dependent tensile axial load induced by the geometrical nonlinearity given as
(141)N~=12·EIPl·∫0l(θx′)2dx.
For the case of constant along the bar axial load equations ((139) and (140)) hold by setting N~=N-(l,t), where N-(l,t) is the externally applied axial force at the bar’s right end. Substituting (139) and (140) into (136a) and (136b) the reduced initial boundary value problem is established as
(142)ρIP·θ¨x-ρCS·θ¨x′′-(GIt+IPA·N~)θx′′+ECSθx′′′′-32EIn·(θx′)2θx′′=mt(x,t)+∂∂x[mw(x,t)],
where the pertinent initial and boundary conditions are appropriately modified, while the boundary value problem related to the secondary warping function ϕSS ((134) and (135)) is accordingly modified to(143a)∇2(ϕSS)red.=EGθx′θx′′·[IPA-(y2+z2)]-(EGθx′′′-ρGθ¨x′)·ϕSP(y,z)inΩ(143b)∂(ϕSS)red.∂n=0onΓj(j=1,2,…,K).It is worth here noting that In in (142) is a nonnegative geometric cross sectional property, related to the geometrical nonlinearity, defined from (121a), which after taking into account the double symmetry of the cross section reduces to
(144)In=IPP-IP2A.From the examined example problems presented in [111, 112] it is concluded that (i)
the geometrical nonlinearity leads to coupling between the torsional and axial equilibrium equations and alters the modeshapes of vibration;(ii)
large twisting rotations increase the stiffness of the bar, leading to higher natural frequencies and better response against torsional loading;(iii)
a tensile axial force is induced in the bar for special cases of axial boundary conditions, due to the geometrical nonlinearity;(iv)
in the treated examples, the axial inertia termρA·u¨m exhibits a time varying influence on the warping displacements and shear stresses, while its effects are more pronounced at the instants at which the aforementioned quantities get local extreme values;(v)
geometrical nonlinearity, dynamic loading and axial boundary conditions influence the warping shear stress distribution and magnitude of bars undergoing twisting deformations;(vi)
warping shear stresses should not be neglected in nonuniform nonlinear torsional vibrations of bars.


## 8.1. Equations of Local Equilibrium

Following the same procedure presented in Section7, the local equilibrium equations are established employing the principle of virtual work
(130)δWint+δWmass=δWext,
where(131a)δWint=∫V(Sxxδεxx+Sxyδγxy+Sxzδγxz)dV,(131b)δWmass=∫Vρ·(u¨δu+v¨δv+w¨δw)dV,(131c)δWext=∫F(tx·δu+ty·δv+tz·δw)dFunder a Total Lagrangian formulation. In the above equations, δ(·) denotes virtual quantities, (·) denotes differentiation with respect to time, V, F are the volume and the surface of the bar, respectively, at the initial configuration, and tx, ty, tz are the components of the traction vector with respect to the undeformed surface of the bar.Neglecting virtual terms of the secondary warping functionϕSS, following the technique presented in [100] and after some algebra [112], the following local equilibrium equation
(132)∂Sxx∂x+∂Sxy∂y+∂Sxz∂z-ρu¨=0inΩ,∀x∈(0,l)
along with its corresponding boundary condition
(133)Sxyny+Sxznz=0onΓj,∀x∈(0,l)
is obtained. Requiring both the primary and the secondary due to warping parts of (132) and (133) to vanish, employing (106a), (106b), and (106c) after taking into account the double symmetry of the cross section and ignoring the term ϕ¨SS, the following governing equation for the secondary ϕSS warping function is obtained as
(134)∇2ϕSS=(-EGum′′+ρGu¨m)-(EGθx′′′-ρGθ¨x′)·ϕSP(y,z)-EGθx′θx′′·(y2+z2)inΩ
along with its corresponding boundary condition
(135)∂ϕSS∂n=0onΓj(j=1,2,…,K)
while as expected, the governing equation related to the primary ϕSP warping function is found to be coincident with the well-known St. Venant corresponding boundary value problem given from (19a) and (19b). As it has already been mentioned for the static case, the evaluation of the warping function contains an integration constant (parallel displacement of the cross section along the bar axis), which can be obtained from (110), arising from the request that the torsional terms of the displacement field do not result in any axial forces. It is worth pointing out that any other constraints could be used, although the use of (110) decouples the governing equations of the torsional problem at the greatest extent.


## 8.2. Equations of Global Equilibrium

Taking into account the double symmetry of the cross section and substituting the stress components ((106a), (106b), and (106c)) the strain ones ((103a), (103b), and (103c)) and the displacement components ((129a), (129b), and (129c)) to the principle of virtual work (130), the governing partial differential equations of the bar are obtained after some algebra as(136a)ρA·u¨m-EA·um′′-EIP·θx′θx′′=n(x,t),(136b)ρIP·θ¨x-ρCS·θ¨x′′-GItθx′′+ECSθx′′′′-32EIPP·(θx′)2θx′′-EIP·um′θx′′-EIP·um′′θx′=mt(x,t)+∂∂x[mw(x,t)]subjected to the initial conditions (x∈(0,l))(137a)um(x,0)=u-m0(x),(137b)u˙m(x,0)=u-˙m0(x),(137c)θx(x,0)=θ-x0(x),(137d)θ˙x(x,0)=θ-˙x0(x)together with the boundary conditions at the bar ends x=0,l(138a)a1N+α2um=α3,(138b)β1Mt+β2θx=β3,(138c)β-1Mw+β-2θx′=β-3,where N, Mt, Mw are the axial force, the twisting and warping moments at the bar ends given from (116a), (126a), and (126b), respectively. The expressions of the externally applied loads appearing in the right hand side of (136a) and (136b) with respect to the first Piola-Kirchhoff stress components can be easily deduced by virtue of (131c). It is worth here noting that damping could also be included in the analysis without any special difficulty, while the geometric cross sectional properties appearing in (136a) and (136b) are also given in Section 7.The solution of the initial boundary value problem described by (136a)–(138c) for the evaluation of the unknown kinematical components um(x,t), θx(x,t) assumes that the warping CS and the torsion It constants are already established; the evaluation of which presumes that the primary warping function ϕSP at any interior point of the domain Ω of the cross section of the bar is evaluated. Once these components are established, the secondary warping function ϕSS at any interior point of the domain Ω of the cross section of the bar is evaluated after solving the boundary value problem described by (134) and (135). Subsequently the second Piola-Kirchhoff stress components are evaluated employing (106a), (106b), and (106c) completing the computation of the stress field.A significant reduction on both the set of the governing differential equations and the boundary value problem related to the secondary warping functionϕSS can be achieved by neglecting the axial inertia term ρA·u¨m of (136a), an assumption which is common among various dynamic beam formulations. Ignoring this term, a single partial differential equation along with a single unknown kinematical component (the angle of twist θx(x,t)) is obtained, which is further simplified in the case of vanishing distributed axial load along the bar. In what follows, this procedure is described in detail for the cases of axially immovable ends and constant axial load along the bar, which are of great practical interest.


## 8.3. Reduced Problems for Special Cases of Axial Boundary Conditions

For the case of axially immovable ends, it is easily proved that [112]
(139)um′′=-IPA·θx′θx′′,∀x∈(0,l)
which after subsequent integration yields
(140)um′=-12·IPA·(θx′)2+N~EA,∀x∈[0,l],
where N~ is a time-dependent tensile axial load induced by the geometrical nonlinearity given as
(141)N~=12·EIPl·∫0l(θx′)2dx.
For the case of constant along the bar axial load equations ((139) and (140)) hold by setting N~=N-(l,t), where N-(l,t) is the externally applied axial force at the bar’s right end. Substituting (139) and (140) into (136a) and (136b) the reduced initial boundary value problem is established as
(142)ρIP·θ¨x-ρCS·θ¨x′′-(GIt+IPA·N~)θx′′+ECSθx′′′′-32EIn·(θx′)2θx′′=mt(x,t)+∂∂x[mw(x,t)],
where the pertinent initial and boundary conditions are appropriately modified, while the boundary value problem related to the secondary warping function ϕSS ((134) and (135)) is accordingly modified to(143a)∇2(ϕSS)red.=EGθx′θx′′·[IPA-(y2+z2)]-(EGθx′′′-ρGθ¨x′)·ϕSP(y,z)inΩ(143b)∂(ϕSS)red.∂n=0onΓj(j=1,2,…,K).It is worth here noting that In in (142) is a nonnegative geometric cross sectional property, related to the geometrical nonlinearity, defined from (121a), which after taking into account the double symmetry of the cross section reduces to
(144)In=IPP-IP2A.From the examined example problems presented in [111, 112] it is concluded that (i)
the geometrical nonlinearity leads to coupling between the torsional and axial equilibrium equations and alters the modeshapes of vibration;(ii)
large twisting rotations increase the stiffness of the bar, leading to higher natural frequencies and better response against torsional loading;(iii)
a tensile axial force is induced in the bar for special cases of axial boundary conditions, due to the geometrical nonlinearity;(iv)
in the treated examples, the axial inertia termρA·u¨m exhibits a time varying influence on the warping displacements and shear stresses, while its effects are more pronounced at the instants at which the aforementioned quantities get local extreme values;(v)
geometrical nonlinearity, dynamic loading and axial boundary conditions influence the warping shear stress distribution and magnitude of bars undergoing twisting deformations;(vi)
warping shear stresses should not be neglected in nonuniform nonlinear torsional vibrations of bars.


## 9. Inelastic Nonuniform Torsion of Bars

Design of bars and bar assemblages subjected to twisting moments based on elastic analysis is most likely to be extremely conservative not only due to significant difference between initial yield and full plastification in a cross section, but also due to the unaccounted for yet significant reserves of strength that are not mobilized in redundant members until after inelastic redistribution takes place. Thus, material nonlinearity is important for investigating the ultimate strength of a bar that resists torsional loading, while distributed plasticity models are acknowledged in the literature [114–116] to capture more rigorously material nonlinearities than cross sectional stress resultant approaches [117] or lumped plasticity idealizations [118, 119]. If inelastic effects are considered, especially through distributed plasticity formulations, then the nonuniform [120, 121] torsional problem requires a much more rigorous analysis.Apart from research efforts in which bars are idealized with computationally demanding three dimensional [122] or shell [123] elements, several researchers proposed specialized beam elements to analyze bars under inelastic nonuniform torsion. Bathe and Wiener [123] employed Hermitian and isoparametric beam elements to model thin-walled I-beams using one element for each flange and one element for the web of the cross section. Wunderlich et al. [124] employed a power series numerical technique using an Updated Lagrangian formulation to study thin-walled beams under general loading conditions. Izzuddin and Lloyd Smith [125, 126] followed an Eulerian approach employing cubic shape functions to study thin walled bars under general loading conditions. Chen and Trahair [99] and Pi and Trahair [121], using a mitre model [127] to describe the shear strain distribution over the cross section, developed finite element models for the inelastic analysis of thin-walled I-beams under torsion. Nie and Zhong [128] and Wang et al. [129] studied thin-walled beams under general loading conditions by employing an Updated Lagrangian description and the FE method.The aforementioned contributions are applicable to cross sections of special geometry (e.g., thin-walled ones). A more general approach is presented by Baba and Kajita [120] who used a two-node, four-degree-of-freedom element for the longitudinal modeling and a four-node, 12-degree-of-freedom rectangular element for the cross sectional modeling of bars of bisymmetrical section under a Total Lagrangian formulation, taking into account plasticity effects in the plane of the cross section as well. Moreover, in the more recent contributions of Battini and Pacoste [130], Nukala and White [114] and Gruttmann et al. [131], beams of arbitrarily shaped cross section under general loading conditions are analyzed. Battini and Pacoste [130] and Nukala and White [114] use an independent warping parameter associated with the intensity of warping deformations to model nonuniform torsion. This parameter is eventually equated with the angle of twist per unit length in their analyses, thus only St. Venant shear stresses are taken into account. In all of the aforementioned research efforts warping shear stresses are not taken into account, with the exceptions of Wunderlich et al. [124], Gruttmann et al. [131], Nie and Zhong [128] and Wang et al. [129] who exploit a stress distribution arising from the introduction of an independent warping parameter in the displacement field of the bar. However, since this theory is analogous to the Timoshenko beam theory of shear-bending loading conditions, it does not satisfy local equilibrium equations under inelastic or even elastic conditions (for relevant discussions see e.g., Simo et al. [132] and Minghini et al. [133]). Finally, Sapountzakis and Tsipiras in [134] presented a boundary element solution for the inelastic nonuniform torsional problem of simply or multiply connected cylindrical bars of arbitrarily shaped doubly symmetric cross section taking into account the effect of warping shear stresses. The bar is subjected to arbitrarily distributed or concentrated torsional loading along its length, while its edges are subjected to the most general torsional boundary conditions. In this last research effort a displacement based formulation is developed and inelastic redistribution is modelled through a distributed plasticity model exploiting three-dimensional material constitutive laws and numerical integration over the cross sections. An incremental-iterative solution strategy is adopted to restore global equilibrium along with an efficient iterative process to integrate the inelastic rate equations [135]. Three boundary value problems with respect to the variable along the bar axis angle of twist to the primary and to the secondary warping functions are formulated and solved employing the boundary element method [23].In order to formulate the inelastic torsional problem, let us consider a prismatic bar of lengthl (Figure 15) with an arbitrarily shaped doubly symmetric constant cross section, occupying the two dimensional multiply connected region Ω of the y, z plane bounded by the Γj(j=0,1,2,…,K) boundary curves, which are piecewise smooth; that is, they may have a finite number of corners. In Figure 15(b)Syz is an arbitrary coordinate system through the cross section’s shear center. The normal stress-strain relationship for the material is assumed to be elastic-plastic-strain hardening with initial modulus of elasticity, shear modulus, and yield stress E, G, and σY0, respectively. The bar is subjected to the combined action of arbitrarily distributed or concentrated twisting mt=mt(x) and warping mw=mw(x) moments acting in the x direction (Figure 15(a)).Prismatic bar subjected to twisting and warping moments (a) of an arbitrarily shaped doubly symmetric cross section, occupying regionΩ (b).
(a)
(b)Under the aforementioned loading, the displacement field of the bar taking into account warping shear stresses is assumed to be given as(145a)u(x,y,z)=θx′(x)ϕSP(y,z)+θx′′′(x)ϕSS(y,z),(145b)v(x,z)=-zθx(x),(145c)w(x,y)=yθx(x),where u, v, w are the axial and transverse bar displacement components with respect to the Syz system of axes; θx′(x) denotes the rate of change of the angle of twist θx(x) regarded as the torsional curvature, while θx′′′(x) is the third order derivative of the angle of twist with respect to x; ϕSP, ϕSS are the primary and secondary warping functions, respectively with respect to the shear center S.Substituting (145a), (145b), and (145c) in the well-known three-dimensional linear strain-displacement relations, the nonvanishing (total) strain resultants are obtained as(146a)εxx=θx′′ϕSP,(146b)γxy=θx′(∂ϕSP∂y-z)+θx′′′∂ϕSS∂y,(146c)γxz=θx′(∂ϕSP∂z+y)+θx′′′∂ϕSS∂z,where in (146a) the higher order contribution of the secondary warping function to the normal strain component is neglected [21]. The first terms of (146b) and (146c) are related to the primary shear stress distribution accounting for uniform torsion, while the last ones are related to the warping shear stress distribution accounting for nonuniform torsion.Considering strains to be small, employing the Cauchy stress tensor and assuming an isotropic and homogeneous material without exhibiting any damage during its plastification, the stress rates are defined in terms of the strain ones as(147){dσxxdτxydτxz}=[E*000G000G]{dεxxeldγxyeldγxzel},
where d(·) denotes infinitesimal incremental quantities over time (rates), the superscript el denotes the elastic part of the strain components and E*=E(1-v)/(1+v)(1-2v). If the plane stress hypothesis is undertaken then E*=E/(1-v2) holds [68], while E is frequently considered instead of E* (E*≈E) in beam formulations [68, 136]. This last consideration is followed here, while any other reasonable expression of E* could also be used without any difficulty.As long as the material remains elastic or elastic unloading occurs, that is,(148){dεxxdγxydγxz}T={dεxxeldγxyeldγxzel}T.
The stress rates are given with respect to the total strain ones by combining (147) and (148). If plastic flow occurs then
(149){dεxxdγxydγxz}T={dεxxeldγxyeldγxzel}T+{dεxxpldγxypldγxzpl}T,
where the superscript pl denotes the plastic part of the strain components. A Von Mises yield criterion and an associated flow rule for the material are considered. The yield condition is described with the expression
(150)f=σxx2+3(τxy2+τxz2)-σY(εeqpl)=0,
where σY is the yield stress of the material and εeqpl is the equivalent plastic strain, the rate of which is defined in [137] and is equal to dεeqpl=dλ (dλ is the proportionality factor). Moreover, the plastic modulus h is defined as h=dσY/dεeqpl or dσY=hdλ and can be estimated from a tension test as h=EtE/(E-Et) (Figure 16). According to the associated flow rule the plastic strain rates are given as
(151){dεxxpldγxypldγxzpl}T=dλ{∂f∂σxx∂f∂τxy∂f∂τxz}T.
Using the aforementioned relation linking the yield stress rate and the proportionality factor, (147), (149)–(151) and exploiting the plastic loading condition (df=0), the stress rates—total strain rates relations are resolved as
(152){dσxxdτxydτxz}=[Del-pl]{dεxxdγxydγxz},
where [Del-pl] is the elastoplastic constitutive matrix given as
(153)[Del-pl]=1c[c11sym.c21c22c31c32c33]
with(154a)c=hσe2+Eσxx2+9G(τxy2+τxz2),(154b)c11=E[hσe2+9G(τxy2+τxz2)],(154c)c21=-3EGσxxτxy,(154d)c31=-3EGσxxτxz,(154e)c22=G[hσe2+Eσxx2+9Gτxz2],(154f)c32=-9G2τxyτxz,(154g)c33=G[hσe2+Eσxx2+9Gτxy2],(154h)σe=σxx2+3(τxy2+τxz2).By setting h=0 in the above relations, the constitutive matrix presented by Baba and Kajita [120] is obtained, while if one of the shear stress components (along with the corresponding strain one) is dropped out, the constitutive relations presented by Chen and Trahair [99] are also precisely recovered.Normal stress-strain (a) and yield stress-equivalent plastic strain (b) relationships.
(a)
(b)


### 9.1. Equations of Local Equilibrium

The primary and secondary warping functions are resolved by exploiting local equilibrium considerations along the longitudinalx axis (together with the associated boundary condition at the lateral surface of the bar) under elastic conditions. Thus, ϕSP, ϕSS are evaluated by solving the boundary value problems(155a)∇2ϕSP=0inΩ,(155b)∂ϕSP∂n=zny-ynzonΓj,(156a)
∇
2
ϕ
S
S
=
-
E
G
ϕ
S
P
(
y
,
z
)
in
Ω
,
(156b)
∂
ϕ
S
S
∂
n
=
0
on
Γ
j
which are similar to those presented in Section 2 ((19a), (19b), (20a), and (20b)). It is worth noting that the difference in the governing equation of the secondary warping function is due to the definition of the axial displacement component (145a).Apparently, local equilibrium is violated when (155a)–(156b) are introduced in (152) to perform stress calculations and inelastic redistribution is modelled only through the angle of twist θx(x). Equations (155a)–(156b) exhibit less accuracy as plastic deformations throughout the cross section increase. Introducing appropriate plastic terms in (155a) and (155b) to satisfy local equilibrium has already been achieved by Wagner and Gruttmann [138] and Sapountzakis and Tsipiras [139] for the uniform torsion problem under small and large displacements, respectively, and by Baba and Kajita [120] for the geometrically nonlinear nonuniform torsional problem of bars of bisymmetrical cross section. However, none of these research efforts accounts for warping shear stresses, while it is pointed out in the literature [130] that such efforts lead to a much higher computational cost. An interesting and simple solution to the same problem is presented by Billinghurst et al. [127] and by Chen and Trahair [99] who employ a “mitre model” to distribute the primary shear stresses along the cross section in such a manner so as to resemble the shear stress distribution of bars under uniform torsion just before plastic collapse. However, this approach is limited in thin-walled bisymmetrical I-shaped cross sections, while it is expected to be less accurate in the elastic part of the bar and during the early phases of plastic loading.


### 9.2. Equations of Global Equilibrium

To establish global equilibrium equations, the principle of virtual work neglecting body forces is employed as follows:(157)∫V(σxxδεxx+τxyδγxy+τxzδγxz)dV=∫F(txδu+tyδv+tzδw)dF,
where δ(·) denotes virtual quantities, tx, ty, tz are the components of the traction vector, V is the volume, and F is the surface area of the bar including the end cross sections.In the elastic case, shear stresses in the left hand side of (12) yield four sets of terms depending on θx′δθx′, θx′δθx′′′, θx′′′δθx′, and θx′′′δθx′′′, respectively. The first term is related to primary shear stresses, the second and third terms are proved to vanish by exploiting (155a)–(156b) and the homogeneity of the constitutive matrix, while the last term, related to warping shear stresses, yields a sixth order derivative of θx in the (Euler-Lagrange) governing equation of the bar. By dropping this term, the fourth-order governing differential equation of the bar under nonuniform torsion
(158)ECSθx′′′′-GItθx′′=mt(x)+dmw(x)dx
is obtained, where It, CS are the torsion (St. Venant) and warping constants, respectively, depending exclusively on the geometry of the cross section and defined from (32a) and (32b). Till now, in the elastic range warping shear stresses do not affect global equilibrium. However, they are still evaluated through (146a)–(148), following the procedure presented in Section 2.In the inelastic case, the linearized version of (157) neglecting body forces
(159)∫V[(σxx+Δσxx)δΔεxx+(τxy+Δτxy)δΔγxy+(τxz+Δτxz)δΔγxz]dV=∫F(txδΔu+tyδΔv+tzδΔw)dF
is employed, where Δ(·) denotes incremental quantities (over time). The incremental strains are obtained employing (146a), (146b), and (146c). As in the elastic case, a sixth order governing differential equation is formulated due to the third order derivative θx′′′ appearing in (146b) and (146c). By dropping terms of fifth order or higher, a fourth-order governing differential equation is once again formulated. However, contrary to the elastic case, the inhomogeneity of the constitutive equations (152) (dependence on the longitudinal coordinate x) and the property of the fully populated constitutive matrix lead to nonvanishing terms of fourth or lower order related to warping shear stresses. Since all these terms depend on ϕSS, they are neglected in the formulation in order to reduce the computational cost of obtaining the tangent stiffness matrix of the bar. It is pointed out that this process affects only the convergence rate of the algorithm and not the equilibrium path. Moreover, the virtual terms δΔγxy, δΔγxz, and δΔu depend on the third order derivative δΔθx′′′ which will eventually lead to the requirement of an additional (third) boundary condition apart from the well-known torsional and warping ones. In order to avoid the formulation of a third-order shear deformation theory, terms of δΔθx′′′ are neglected in (159), resulting in an approximate solution of the problem at hand. Nevertheless, it is pointed out that warping shear stresses still affect the global equilibrium of the bar, without being completely neglected due to the aforementioned simplifications.The stress resultants of the bar are defined as (see also Section2)(160a)SMt=∫Ω[τxy(∂ϕSP∂y-z)+τxz(∂ϕSP∂z+y)]dΩ,(160b)SMw=∫ΩσxxϕSPdΩ,where SMt and SMw correspond to the internal twisting and warping moments of the bar, respectively, and the relevant incremental stress resultants are defined as(161a)ΔSMt=∫Ω[Δτxy(∂ϕSP∂y-z)+Δτxz(∂ϕSP∂z+y)]dΩ,(161b)ΔSMw=∫ΩΔσxxϕSPdΩ.By exploiting (145a), (145b), (145c), (146a), (146b), (146c), (160a), (160b), (161a), and (161b) and the aforementioned simplifications, the governing equation of the bar is obtained after some algebra through (159) as
(162)-dΔSMtdx+d2ΔSMwdx2︸tangentstiffnessmatrix×incrementaldisplacements=mt(x)+dmw(x)dx︸externallyappliedloading+dSMtdx-d2SMwdx2︸internalstressresultants︷out-of-balanceforces
along with its corresponding boundary conditions(163a)α1(ΔSMt-dΔSMwdx)+α2Δθx=α3-α1(SMt-dSMwdx),(163b)β1ΔSMw+β2Δθx′=β3-β1SMw,where αi,βi(i=1,2,3) are functions specified at the bar ends. The boundary conditions (163a) and (163b) are the most general ones for the problem at hand, including also the elastic support. It is apparent that all types of the conventional boundary conditions (clamped, simply supported, free, or guided edge) may be derived from (163a) and (163b) by specifying appropriately the functions αi and βi (e.g., for a clamped edge it is α2=β2=1, α1=α3=β1=β3=0). Finally, the expressions of the externally applied loading (mt, mw) with respect to the components of the traction vector (or the Cauchy stress tensor) may be easily resolved by virtue of (159) as(164a)mt(x)=∫Γ[ty(-z)+ytz]ds,(164b)mw(x)=-∫ΓtxϕSPds.From the above presentation and the examined example problems presented in [134] it is concluded that (i)
both St. Venant and warping shear stresses are evaluated from the solution of two boundary value problems formulated under elastic conditions;(ii)
warping shear stresses slightly decrease the torsional rigidity and the ultimate torsional loading that bars can undertake under nonuniform torsion;(iii)
warping shear stresses have an influence on inelastic redistribution patterns of bars under nonuniform torsion.


## 9.1. Equations of Local Equilibrium

The primary and secondary warping functions are resolved by exploiting local equilibrium considerations along the longitudinalx axis (together with the associated boundary condition at the lateral surface of the bar) under elastic conditions. Thus, ϕSP, ϕSS are evaluated by solving the boundary value problems(155a)∇2ϕSP=0inΩ,(155b)∂ϕSP∂n=zny-ynzonΓj,(156a)
∇
2
ϕ
S
S
=
-
E
G
ϕ
S
P
(
y
,
z
)
in
Ω
,
(156b)
∂
ϕ
S
S
∂
n
=
0
on
Γ
j
which are similar to those presented in Section 2 ((19a), (19b), (20a), and (20b)). It is worth noting that the difference in the governing equation of the secondary warping function is due to the definition of the axial displacement component (145a).Apparently, local equilibrium is violated when (155a)–(156b) are introduced in (152) to perform stress calculations and inelastic redistribution is modelled only through the angle of twist θx(x). Equations (155a)–(156b) exhibit less accuracy as plastic deformations throughout the cross section increase. Introducing appropriate plastic terms in (155a) and (155b) to satisfy local equilibrium has already been achieved by Wagner and Gruttmann [138] and Sapountzakis and Tsipiras [139] for the uniform torsion problem under small and large displacements, respectively, and by Baba and Kajita [120] for the geometrically nonlinear nonuniform torsional problem of bars of bisymmetrical cross section. However, none of these research efforts accounts for warping shear stresses, while it is pointed out in the literature [130] that such efforts lead to a much higher computational cost. An interesting and simple solution to the same problem is presented by Billinghurst et al. [127] and by Chen and Trahair [99] who employ a “mitre model” to distribute the primary shear stresses along the cross section in such a manner so as to resemble the shear stress distribution of bars under uniform torsion just before plastic collapse. However, this approach is limited in thin-walled bisymmetrical I-shaped cross sections, while it is expected to be less accurate in the elastic part of the bar and during the early phases of plastic loading.


## 9.2. Equations of Global Equilibrium

To establish global equilibrium equations, the principle of virtual work neglecting body forces is employed as follows:(157)∫V(σxxδεxx+τxyδγxy+τxzδγxz)dV=∫F(txδu+tyδv+tzδw)dF,
where δ(·) denotes virtual quantities, tx, ty, tz are the components of the traction vector, V is the volume, and F is the surface area of the bar including the end cross sections.In the elastic case, shear stresses in the left hand side of (12) yield four sets of terms depending on θx′δθx′, θx′δθx′′′, θx′′′δθx′, and θx′′′δθx′′′, respectively. The first term is related to primary shear stresses, the second and third terms are proved to vanish by exploiting (155a)–(156b) and the homogeneity of the constitutive matrix, while the last term, related to warping shear stresses, yields a sixth order derivative of θx in the (Euler-Lagrange) governing equation of the bar. By dropping this term, the fourth-order governing differential equation of the bar under nonuniform torsion
(158)ECSθx′′′′-GItθx′′=mt(x)+dmw(x)dx
is obtained, where It, CS are the torsion (St. Venant) and warping constants, respectively, depending exclusively on the geometry of the cross section and defined from (32a) and (32b). Till now, in the elastic range warping shear stresses do not affect global equilibrium. However, they are still evaluated through (146a)–(148), following the procedure presented in Section 2.In the inelastic case, the linearized version of (157) neglecting body forces
(159)∫V[(σxx+Δσxx)δΔεxx+(τxy+Δτxy)δΔγxy+(τxz+Δτxz)δΔγxz]dV=∫F(txδΔu+tyδΔv+tzδΔw)dF
is employed, where Δ(·) denotes incremental quantities (over time). The incremental strains are obtained employing (146a), (146b), and (146c). As in the elastic case, a sixth order governing differential equation is formulated due to the third order derivative θx′′′ appearing in (146b) and (146c). By dropping terms of fifth order or higher, a fourth-order governing differential equation is once again formulated. However, contrary to the elastic case, the inhomogeneity of the constitutive equations (152) (dependence on the longitudinal coordinate x) and the property of the fully populated constitutive matrix lead to nonvanishing terms of fourth or lower order related to warping shear stresses. Since all these terms depend on ϕSS, they are neglected in the formulation in order to reduce the computational cost of obtaining the tangent stiffness matrix of the bar. It is pointed out that this process affects only the convergence rate of the algorithm and not the equilibrium path. Moreover, the virtual terms δΔγxy, δΔγxz, and δΔu depend on the third order derivative δΔθx′′′ which will eventually lead to the requirement of an additional (third) boundary condition apart from the well-known torsional and warping ones. In order to avoid the formulation of a third-order shear deformation theory, terms of δΔθx′′′ are neglected in (159), resulting in an approximate solution of the problem at hand. Nevertheless, it is pointed out that warping shear stresses still affect the global equilibrium of the bar, without being completely neglected due to the aforementioned simplifications.The stress resultants of the bar are defined as (see also Section2)(160a)SMt=∫Ω[τxy(∂ϕSP∂y-z)+τxz(∂ϕSP∂z+y)]dΩ,(160b)SMw=∫ΩσxxϕSPdΩ,where SMt and SMw correspond to the internal twisting and warping moments of the bar, respectively, and the relevant incremental stress resultants are defined as(161a)ΔSMt=∫Ω[Δτxy(∂ϕSP∂y-z)+Δτxz(∂ϕSP∂z+y)]dΩ,(161b)ΔSMw=∫ΩΔσxxϕSPdΩ.By exploiting (145a), (145b), (145c), (146a), (146b), (146c), (160a), (160b), (161a), and (161b) and the aforementioned simplifications, the governing equation of the bar is obtained after some algebra through (159) as
(162)-dΔSMtdx+d2ΔSMwdx2︸tangentstiffnessmatrix×incrementaldisplacements=mt(x)+dmw(x)dx︸externallyappliedloading+dSMtdx-d2SMwdx2︸internalstressresultants︷out-of-balanceforces
along with its corresponding boundary conditions(163a)α1(ΔSMt-dΔSMwdx)+α2Δθx=α3-α1(SMt-dSMwdx),(163b)β1ΔSMw+β2Δθx′=β3-β1SMw,where αi,βi(i=1,2,3) are functions specified at the bar ends. The boundary conditions (163a) and (163b) are the most general ones for the problem at hand, including also the elastic support. It is apparent that all types of the conventional boundary conditions (clamped, simply supported, free, or guided edge) may be derived from (163a) and (163b) by specifying appropriately the functions αi and βi (e.g., for a clamped edge it is α2=β2=1, α1=α3=β1=β3=0). Finally, the expressions of the externally applied loading (mt, mw) with respect to the components of the traction vector (or the Cauchy stress tensor) may be easily resolved by virtue of (159) as(164a)mt(x)=∫Γ[ty(-z)+ytz]ds,(164b)mw(x)=-∫ΓtxϕSPds.From the above presentation and the examined example problems presented in [134] it is concluded that (i)
both St. Venant and warping shear stresses are evaluated from the solution of two boundary value problems formulated under elastic conditions;(ii)
warping shear stresses slightly decrease the torsional rigidity and the ultimate torsional loading that bars can undertake under nonuniform torsion;(iii)
warping shear stresses have an influence on inelastic redistribution patterns of bars under nonuniform torsion.


## 10. Inelastic Nonuniform Torsion of Bars Including Secondary Torsional Moment Deformation Effect

As it has been already mentioned in the previous section, material nonlinearity is important for investigating the ultimate strength of a bar that resists torsional loading, while distributed plasticity models are acknowledged in the literature [114–116] to capture more rigorously material nonlinearities than cross sectional stress resultant approaches [117] or lumped plasticity idealizations [118, 119]. Moreover, in the nonuniform torsional analysis of a bar except from the primary (St. Venant) shear stress distribution forming the primary torsional moment stress resultant additional normal and secondary (warping) shear stresses arise, forming the warping moment and secondary torsional moment stress resultants, respectively. In order to include warping shear stresses in the global equilibrium of the bar, that is, to account for the secondary torsional moment deformation effect (STMDE), an additional kinematical component (along with the angle of twist) is generally required [81–140], increasing the difficulty of the problem at hand.The STMDE has been shown in the literature to be significant, especially on closed shaped section bars. Massonnet [141] presented a qualitative explanation why warping shear stresses are of the same order of magnitude as primary ones in the case of closed shaped section bars. As early as 1954, Benscoter [142] analyzed the nonuniform torsional problem of multicell section bars. Since then, a significant amount of relevant contributions has appeared in the literature as well [109, 140, 143–146]. Since the topic at hand is analogous to the geometrically nonlinear Timoshenko beam theory of shear-bending loading conditions [78, 83], it does not satisfy local equilibrium equations (for relevant discussions see e.g., Simo et al. [132] and Minghini et al. [133]). This problem is alleviated by introducing torsional shear correction factor at the global level [78, 84, 85] and suitable warping shear stress distribution at the local level [80, 85, 147]. The aforementioned contributions refer to the linear elastic regime. If inelastic effects are considered, especially through distributed plasticity formulations, then the nonuniform [120, 122] torsion problem including STMDE requires a much more rigorous analysis.In the majority of the research efforts concerning the inelastic nonuniform torsional problem of bars (already presented in the previous section) torsional shear correction factor is not included in the analyses, while the employed warping shear stress distributions though included in the global equilibrium of the bar they do not satisfy local equilibrium considerations under inelastic or even elastic conditions. A different approach is undertaken in the recent contribution of Sapountzakis and Tsipiras [139], where a single (fourth-order) governing differential equation is formulated with respect to a single kinematical component (angle of twist), leading to an approximate solution of the problem at hand. However, the adopted warping shear stress distribution verifies local equilibrium equations under elastic conditions. Moreover, an alternative methodology is presented in the very recent contribution of Le Corvec and Filippou [148], where a multitude of local section warping degrees of freedom are introduced in order to model effects such as in-plane plasticity effects, constrained warping, shear lag, and so forth. This technique does not require the introduction of torsional shear correction factor; however, the reduction of the number of employed degrees of freedom requires further investigation, while an example of a thick-walled cross section beam is not presented. In the publication of Wackerfuß and Gruttmann [149], beams of thick-walled rectangular cross sections are worked out by employing a series of polynomials as global warping functions in order to capture advanced effects such as in-plane inelastic redistribution and transverse contraction. Moreover, the same researchers [150] have employed local warping functions to capture the aforementioned effects in moderately thick arbitrarily shaped cross section beams; however, bimoment loading cannot be applied at the bar. Finally, Tsipiras and Sapountzakis [151] developed a boundary element solution for the inelastic nonuniform torsional problem of simply or multiply connected cylindrical bars of arbitrarily shaped doubly symmetric cross section taking into account the effect of secondary torsional moment deformation. The bar is subjected to arbitrarily distributed or concentrated torsional loading along its length, while its edges are subjected to the most general torsional boundary conditions. A displacement based formulation is developed and inelastic redistribution is modelled through a distributed plasticity model exploiting three dimensional material constitutive laws and numerical integration over the cross sections. An incremental-iterative solution strategy is adopted to resolve the elastic and plastic part of stress resultants along with an efficient iterative process to integrate the inelastic rate equations [135]. The one dimensional primary angle of twist per unit length, a two dimensional secondary warping function, and a scalar torsional shear correction factor are employed to account for STMDE. The latter is computed employing an energy approach under elastic conditions [85].In order to formulate the aforementioned problem, let us consider the prismatic bar of Figure15 subjected to the combined action of arbitrarily distributed or concentrated twisting mt=mt(x) and warping mw=mw(x) moments acting in the x direction (Figure 15(a)). Under the aforementioned loading, the displacement field of the bar taking into account warping shear stresses is assumed to be given as(165a)u(x,y,z)=(θxP(x))′ϕSP(y,z),(165b)v(x,z)=-zθx(x),(165c)w(x,y)=yθx(x),where u, v, w are the axial and transverse bar displacement components with respect to the Syz system of axes; θx is the (total) angle of twist; (θxP)′ is the primary angle of twist per unit length which is in general not equal to the angle of twist per unit length θx′ (see also (65a); ϕSP is the primary warping function with respect to the shear center S.Substituting (165a), (165b), and (165c) in the well-known three-dimensional linear strain-displacement relations, the nonvanishing compatible (total) strain resultants are obtained as(166a)εxx=(θxP)′′ϕSP,(166b)γxy=θx′(∂ϕSP∂y-z)︸primary+((θxP)′-θx′)∂ϕSP∂y︸secondary,(166c)γxz=θx′(∂ϕSP∂z+y)︸primary+((θxP)′-θx′)∂ϕSP∂z︸secondary,where in (166b) and (166c), the first terms correspond to St. Venant shear strains (γxyP, γxzP) and the last terms to warping shear strains (γxyS, γxzS). In order to formulate global equilibrium equations at the elastic regime that include torsional shear correction factor, the above relations are proposed to be corrected as(167a)εxx=(θxP)′′ϕSP,(167b)γxy=θx′(∂ϕSP∂y-z)︸primary+((θxP)′-θx′)κx∂ϕSP∂y︸secondary,(167c)γxz=θx′(∂ϕSP∂z+y)︸primary+((θxP)′-θx′)κx∂ϕSP∂z︸secondary,where κx is the torsional shear correction factor, while the specific form of the correction (κx) is justified from the use of an energy approach to evaluate κx [151]. The elastic shear stress distribution arising from (167a), (167b), and (167c) yields equilibrium equations that are corrected at the global level; however, it still violates the longitudinal local equilibrium equation. Thus a two dimensional secondary warping function ϕSS(y,z) [151] is introduced as(168a)εxx=(θxP)′′ϕSP,(168b)γxy=θx′(∂ϕSP∂y-z)︸primary+((θxP)′-θx′)(-ItSCS)∂ϕSS∂y︸secondary,(168c)γxz=θx′(∂ϕSP∂z+y)︸primary+((θxP)′-θx′)(-ItSCS)∂ϕSS∂z︸secondary,where ItS and CS are the secondary torsion constant and warping constant, respectively (see also Section 6, [151]). Both the primary and the secondary warping functions are determined by formulating boundary value problems based on the exploitation of the longitudinal local equilibrium equation and the associated boundary condition as presented in previous sections. Formulation arising from (166a), (166b), and (166c) corresponds to the use of constant strain distribution in Timoshenko beam theory of shear-bending loading conditions without employing a shear correction factor and is denoted as model A in this work. Formulations arising from (167a), (167b), (167c), (168a), (168b), and (168c) correspond to the use of constant and parabolic strain distribution in Timoshenko beam theory and are denoted as models B, C, respectively. It is pointed out that the motivation to perform the above corrections on the strain field of the bar (and not on the definition of the secondary torsional moment stress resultant) stems from the fact that most inelastic material state determination algorithms are strain-driven ones, while the adopted methodology conforms to the one presented in [152].Following the procedure presented in Section9, the application of the principle of virtual work neglecting body forces after some algebraic manipulations leads to the following global equilibrium equations of the bar:(169a)dSMtPdx+dSMtSdx=-mt,(169b)dSMwdx+SMtS=mwalong with their corresponding boundary conditions(170a)α1SMttot+α2θx=α3,(170b)β1SMw+β2(θxP)′=β3,where SMttot=SMtP+SMtS is the total twisting moment at the bar ends, while the stress resultants are naturally defined as(171a)SMw=∫Ω(σxxϕSP)dΩ,(171b)SMtP=∫Ω[τxy(∂ϕSP∂y-z)+τxz(∂ϕSP∂z+y)]dΩ,(171c)SMtS=ItSCS∫Ω(τxy∂ϕSS∂y+τxz∂ϕSS∂z)dΩ,where SMw, SMtP, and SMtS correspond to warping moment, primary, and secondary twisting moments, respectively. Equation (171c) gives the expression of SMtS according to model C (strain equations (168a), (168b), and (168c)), while the corresponding ones for models A, B are presented in [151]. In (170a) and (170b) αi,βi(i=1,2,3) are functions specified at the bar ends, while the given boundary conditions are the most general ones for the problem at hand, including also the elastic support. Finally, the expressions of the externally applied loading (mt,mw) with respect to the components of the traction vector (or the Cauchy stress tensor) may be easily resolved as(172a)mt(x)=∫Γ[ty(-z)+ytz]ds,(172b)mw(x)=-∫ΓtxϕSPds.Since an incremental-iterative approach is adopted for the problem at hand, the incremental version of (169a) and (169b) is written down as(173a)dΔSMtPdx+dΔSMtSdx=-Δmt,(173b)dΔSMwdx+ΔSMtS=Δmw,where Δ(·) denotes incremental quantities (over time), while the incremental stress resultants are given by virtue of (171a), (171b), (171c), (147), and (149) as(174a)ΔSMw=ECS(ΔθxP)′′+ΔSMwpl,(174b)ΔSMtP=GItPΔθx′+ΔSMtPpl,(174c)ΔSMtS=-GItS((ΔθxP)′-Δθx′)+ΔSMtSpl,where ItP, CS are the primary torsion (St. Venant) and warping constants, defined from (72b) and (72a), respectively, while ItS is the secondary torsion constant given as
(175)ItS=κx∫Ω[(∂ϕSP∂y)2+(∂ϕSP∂z)2]dΩ(modelC)
and ΔSMwpl, ΔSMtPpl, ΔSMtSpl are plastic quantities defined as(176a)ΔSMwpl=-E∫Ω(ΔεxxplϕSP)dΩ,(176b)ΔSMtPpl=-G∫Ω[Δγxypl(∂ϕSP∂y-z)+Δγxzpl(∂ϕSP∂z+y)]dΩ,(176c)ΔSMtSpl=-GItSCS∫Ω(Δγxypl∂ϕSS∂y+Δγxzpl∂ϕSS∂z)dΩ(modelC)while the expressions of ΔSMtSpl and ItS for models A, B and the details of the derivation of expression (174c) are presented in [151]. Substituting (174a), (174b), and (174c) in (173a), and (173b) and forming the incremental version of the boundary conditions (170a) and (170b), the following boundary value problem is obtained:(177a)-G(ItP+ItS)Δθx′′+GItS(ΔθxP)′′=Δmt+dΔSMtPpldx+dΔSMtSpldxinsidethebar(177b)-ECS(ΔθxP)′′′+GItS((ΔθxP)′-Δθx′)=-Δmw+dΔSMwpldx+ΔSMtSplinsidethebar,(178a)
α
1
Δ
S
M
t
tot
+
α
2
Δ
θ
x
=
Δ
α
3
at
the
bar
ends
x
=
0
,
l
,
(178b)
β
1
Δ
S
M
w
+
β
2
(
Δ
θ
x
P
)
′
=
Δ
β
3
at
the
bar
ends
x
=
0
,
l
.
In order to avoid the third order derivative of ΔθxP appearing in (177b), a one dimensional independent warping parameter ηx is set equal to the primary angle of twist per unit length, thus the boundary value problem of the above equations is reformulated as(179a)-G(ItP+ItS)Δθx′′+GItSΔηx′=Δmt+dΔSMtPpldx+dΔSMtSpldxinsidethebar(179b)-ECSΔηx′′+GItS(Δηx-Δθx′)=-Δmw+dΔSMwpldx+ΔSMtSplinsidethebar,(180a)
α
1
Δ
S
M
t
tot
+
α
2
Δ
θ
x
=
Δ
α
3
at
the
bar
ends
x
=
0
,
l
,
(180b)
β
1
Δ
S
M
w
+
β
2
Δ
η
x
=
Δ
β
3
at
the
bar
ends
x
=
0
,
l
.
Dropping the plastic quantities of the above equations, the boundary value problem of the examined problem under elastic conditions is formulated.From the above presentation and the examined example problems presented in [151] it is concluded that (i)
both St. Venant and warping shear stresses are evaluated from the solution of two boundary value problems formulated under elastic conditions;(ii)
STMD effect is negligible in inelastic nonuniform torsional analysis of open thin-walled cross section bars;(iii)
STMD effect decreases torsional rigidity of closed shaped section bars under nonuniform torsion;(iv)
modification of compatible strains (model A) with the inclusion of torsional shear correction factor (model B) and secondary warping function (model C) decreases torsional rigidity and results in better modelling of bars under inelastic nonuniform torsion, especially in predicting cross sectional plasticity patterns and stress distributions.


## 11. Concluding Remarks

In this paper both the static and dynamic analyses of the geometrically linear or nonlinear, elastic or elastoplastic nonuniform torsion problems of bars of constant or variable arbitrary cross section are presented together with the corresponding research efforts and the conclusions drawn from examined cases with great practical interest. In the presented analyses, the bar is subjected to arbitrarily distributed or concentrated twisting and warping moments along its length, while its edges are supported by the most general torsional boundary conditions. The analysis of the aforementioned problems is complete by presenting the evaluation of the torsion and warping constants of the bar, its displacement field, its stress resultants together with the torsional shear stresses and the warping normal and shear stresses at any internal point of the bar.Having in mind the disadvantages of the 3D FEM solutions and more specifically the difficulties(i)
in support modeling,(ii)
in discretizing a complex structure,(iii)
in discretizing a structure including thin-walled members (shear locking and membrane locking),(iv)
in the increased number of degrees of freedom leading to severe or unrealistic computational time,(v)
in the reduced oversight of the 3D FEM solution compared with that of the beamlike structures employing stress resultants,and the fact that the use of shell elements cannot give accurate results since the warping of the walls of a cross section cannot be taken into account (midline model), the importance of the presented beamlike analyses becomes more evident.



---

*Source: 916581-2013-03-21.xml*