Bars under Torsional Loading: A Generalized Beam Theory Approach

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 andwarping constants of the bar, its displacement field, its stress resultants togetherwith the torsional shear stresses and thewarping 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.


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][4][5][6][7][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.

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][10][11][12][13][14][15][16].Moreover, numerical methods such as the finite element method [17] or the boundary element method [18][19][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 length  (Figure 1), of constant arbitrary cross-section of area .The homogeneous isotropic and linearly elastic material of the bar's cross-section, with modulus of elasticity , shear modulus , and Poisson's ratio ], occupies the two dimensional multiply connected region Ω of the ,  plane and is bounded by the Γ  ( = 1, 2, . . ., ) boundary curves, which are piecewise smooth; that is, they may have a finite number of corners.In Figure 1(b)  is the principal coordinate system through the cross section's centroid , while   ,   are its coordinates with respect to  system of axes through the cross section's shear center .The bar is subjected to the arbitrarily distributed or concentrated conservative twisting   =   () moment acting in the  direction (Figure 1(a)).This torsional loading results in a rotation with respect to the shear center  with angle of twist   ().The shear center  coincides with the center of twist of the cross section provided that no other axis or rotation is imposed due to construction requirements.
Under the aforementioned loading the displacement field of the bar with respect to the  system of axes is given as  (, , ) =    ()   (, ) , V (, , ) = −  () , (, , ) =   () , where the displacement component  constitutes the warping of the cross section,    () =   / denotes the rate of change of the angle of twist   regarded as the torsional curvature, and   (, ) is the warping function with respect to the shear center , which is characterized as the warping function.The index  indicates that the warping refers to the rotation axis .
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 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 ⋅    () ⋅ ( ⋅    () ⋅ ( 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    () = 0, that is, the rate of change of    () is constant (case of uniform Saint Venant torsion), and the first two equilibrium equations are satisfied.From (4c) it follows that Having in mind that the first    and third    derivatives of the angle of twist are functions of the longitudinal coordinate , it can be observed that the solution of the partial differential equation (5) in the general case should be a function of coordinate  as well.This fact contradicts relation (1a), where the warping function   (, ) is defined as independent of the coordinate .Thus, the left hand side of ( 5) is a function only of variables , , since   =   (, ), while at the same time the right hand side is a function of variable .To remove this inconsistency it is assumed that the arising shear stresses are decomposed to a primary   part and a secondary   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    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    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   (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    () 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    .Based on the above decomposition of shear stress in primary and secondary components, that is =    +    (6b) the primary and secondary components of these stresses and the normal stresses due to warping are defined according to the relations as follows: =  ⋅    () ⋅ ( =  ⋅    () ⋅    (, ) , where the functions    (, ) and    (, , ) 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 Equation (10) indicates that the stress state of the bar subjected to torsional loading results from the superposition of primary   , secondary   shear stress and normal stresses due to warping    .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 Substituting (7a)-( 9) into (11) the following relations are obtained: To formulate the boundary 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  , the following relations are obtained: where   = cos(, ) = cos  = / = / and   = sin(, ) = sin  = / = −/ 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 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: 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 Substituting (16a) and (16c) into (17a) and (17b) and taking into account the fact that in general    ̸ = 0, the boundary conditions of the primary and the secondary warping functions are obtained as Summarizing the above mentioned it is concluded that for the determination of the primary    (, ) and the secondary    (, , ) warping functions the solution of the following boundary value problems is required.
For the Primary    (, ) Warping Function.Neumann problem of Laplace differential equation is as follows: on the boundary of the cross section Γ.

(19b)
For the Secondary    (, , ) Warping Function.Neumann problem of the Poisson differential equation is as follows: = 0 on the boundary of the cross section Γ.

(20b)
Finally, based on the above mentioned, the displacement field of the bar ((1a), (1b), and (1c)) is formed as It is worth here noting that in the case the origin  of the coordinate system is a point of the ,  plane other than the shear center (Figure 1  by    .Using the evaluated warping function    ,    is then established using the transformation given by the following equation [2]: where  =  −   ,  =  −   ,   ,   are the coordinates of the shear center  with respect to the arbitrary coordinate system  (Figure 1(b)) and  is an integration constant.The latter is given from the solution of the following linear system of equations: where is the cross section area, the static moments of inertia relative to the axes , , the bending moments of inertia relative to the axes , , and the product of inertia, respectively, while in relations (23a), (23b), and (23c) the warping moments    ,    ,    , are defined as Moreover, the evaluated warping function    from the solution of the Neumann problem (20a) and (20b) contains an integration constant   (parallel displacement of the cross section along the beam axis), which can be obtained from [24] as follows: and the main secondary warping function φ  is given as

Equations of Global Equilibrium.
The already established shear stresses   and   yield components of torque, which arise through integration over the cross section.Thus, the resulting twisting moment is obtained as 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    originating from the primary shear stresses   due to twisting (as in uniform torsion) and a secondary component    , 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: where after employing (6a) and (6b) and some algebra the aforementioned twisting moment components are given as 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: where In (32a) and (32b) the quantity   denotes the torsional constant introduced by Saint-Venant, while the quantity   denotes the warping constant.The quantity   denotes the torsional rigidity, while the quantity   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 .
To derive the equilibrium equation for the nonuniform torsion problem of a homogeneous isotropic bar, the equilibrium of an infinitesimal part  of the bar against twisting moments is examined and observing Figure 5 it follows that or after some algebra 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: which is independent of the torsional boundary conditions.
In analogy with the bending moments   ,   a new stress resultant is defined, called warping moment (or bimoment) and given by The need for the definition of this new stress resultant stems from the fact that, whereas   =   =  = 0, there still exist normal stresses acting on the cross section; hence if a new quantity is not considered, the elastic energy due to   stresses will be ignored.Substituting the expression (9) of the stress component    into (36), the latter is written as or through (32b) and therefore relation ( 9) can be written in the form of 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 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   () 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 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   ,   ( = 1, 2, 3).For example, in the case of a torsionally fixed support the above functions take the values From the examined example problems presented in [18][19][20][21][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.

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 length  (Figure 6), of an arbitrarily shaped variable along its axis cross section.The homogeneous isotropic and linearly elastic material of the bar's crosssection, with modulus of elasticity , shear modulus , and Poisson's ratio ], occupies the two dimensional multiply connected region Ω() of the ,  plane and is bounded by the Γ  ( = 1, 2, . . ., ) boundary curves, which are piecewise smooth; that is, they may have a finite number of corners.In Figure 6(b)  is the principal coordinate system through the cross section's centroid , while   ,   are its coordinates with respect to  system of axes through the cross section's shear center .The bar is subjected to the arbitrarily distributed or concentrated conservative twisting   =   () moment acting in the  direction (Figure 6(a)).
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   =   () is derived: 1   +  2   =  3 at the beam ends  = 0, , (43a) 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 with   the torsional and   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 The warping moment   arising from the torsional curvature, similarly with the constant cross section case, is given from (38) and also functions   ,   ( = 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    (, ) and the secondary    (, , ) 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.

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.Figure 6: Bar subjected to a twisting moment (a) with a variable cross section of arbitrary shape occupying the two dimensional region Ω() (b).

ISRN Civil Engineering
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][36][37][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][44][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 length  (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 , shear modulus , and Poisson's ratio ], occupies the two dimensional multiply connected region Ω of the ,  plane and is bounded by the Γ  ( = 1, 2, . . ., ) boundary curves, which are piecewise smooth; that is, they may have a finite number of corners.In Figure 7(b)  is the principal coordinate system through the cross section's centroid , while   ,   are its coordinates with respect to  system of axes through the cross section's shear center .
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    =   / denoting the rate of change of the angle of twist   , 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 and the respective nodal load vector as where   is the twisting moment at the ends of the element given from (29) and consisting of the primary part    defined as the resultant of the primary shear stress distribution given from (31a) or (44a) and the secondary one    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, for the constant cross section case, (48a) for the variable cross section case.
(48b) Moreover,   is the bending moment due to the torsional curvature at the same sections given from (38).
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 where the    coefficients (,  = 1, 2, ) 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    ( = 1, 2, 3, 4, 5, 6) coefficients are given as [14,35] where In the case of a variable cross section element, the evaluation of the    ( = 1, 2, 3, 4, 5, 6) coefficients presumes the solution of the following boundary value problem with respect to the angle of twist   =   () (see also (42)): for appropriate values of the   ,   ( = 1, 2, 3) functions.Thereby, for the evaluation of the Apparently, for the constant cross section case it is    =    =    = 0 and the governing equation ( 52) reduces to the one given from (35).Upon the evaluation of the angle of twist   , the coefficients    ( = 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 moment   =   () (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   , obtained from the solution of the following boundary value problem: for the constant cross section case: for the variable cross section case: From the examined example problems presented in [43][44][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 the 14 × 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.

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][53][54][55] or distributed mass model systems [56][57][58][59][60] have been presented.In all these references the considered beam is of a constant homogeneous thinwalled 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 length  of Figure 9.The homogeneous isotropic and linearly elastic material of the bar's cross-section, with modulus of elasticity , shear modulus , and Poisson's ratio ], occupies the two dimensional multiply connected region Ω of the ,  plane and is bounded by the Γ  ( = 1, 2, . . ., ) boundary curves, which are piecewise smooth; that is, they may have a finite number of corners.In Figure 9(b)  is the principal coordinate system through the cross section's centroid , while   ,   are its coordinates with respect to  system of axes through the cross section's shear center .The bar is subjected to the time dependent arbitrarily distributed or concentrated twisting moment   =   (, ),  ≥ 0 acting in the  direction (Figure 9(a)).
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) (, , , ) =   (, ) 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   =   () is derived: 1    +  2   =  3 at the beam ends  = 0, , (58b) where θ  (, ) is the first derivative of the angle of twist with respect to time; (), () are the initial angle of twist and the corresponding initial velocity of the points of the beam axis;   =   () is the polar moment of inertia of the cross section about the origin  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 ISRN Civil Engineering  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   arising from the torsional curvature, similarly with the constant cross section case is given from (38).Also, the functions   ,   ( = 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    (, ) and the secondary    (, , , ) 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 From the examined example problems presented in [65][66][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 the 14 × 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.

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 Section 2, 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 length , of an arbitrary doubly symmetric constant cross section with modulus of elasticity  and shear modulus , occupying the two dimensional multiply connected region Ω of the ,  plane bounded by the  + 1 curves Γ 0 , Γ 1 , . . ., Γ  as shown in Figure 10.The bar is subjected to arbitrarily distributed and/or concentrated twisting   =   () and warping   =   () moments, while its edges are restrained by the most general torsional boundary conditions.
In order to account for shear deformation due to restrained torsional warping, the displacements   ,   ,   along the three axes and the rotation   due to twist are spilt into a primary and a secondary part arising from the primary () and the secondary () torsional moment, respectively, as where the (total) first derivative of the angle of twist   / denotes the rate of change of the (total) angle of twist and can be regarded as the (total) torsional curvature.The primary    / and the secondary    / 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.
Assuming small torsional rotation, the displacement field is defined as where    (, ) is the primary warping function with respect to the shear center  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 threedimensional elasticity the nonzero stress components in the region Ω can be written as where    (, , ) is the secondary warping function with respect to the shear center  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 (   /)(   /), (   /)(   /) have been neglected.

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 Similarly with Section 2, requiring both the primary () and the secondary due to warping () 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: where are the warping and the primary torsion constants of the cross section, respectively, while   and    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 Since the secondary warping is much smaller with regard to the primary one [24], the secondary torsional curvature    / 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]: where    is the secondary torsional rigidity of the cross section, while    is the secondary torsion constant of the cross section which can be written as where the coefficient    is called warping shear correction factor and    is the warping shear deformation coefficient.Substituting (75) into (74) the secondary torsional curvature can be written as From ( 76) it is concluded that for    = 0 the secondary shear deformation effect is neglected, that is    / = 0 and    / =   /.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 with the exact one given from the secondary torsion constant can be obtained as In order to formulate this constant independently from the loading and the material properties of the bar, the values    = 1,  = 1 are employed and the secondary torsion constant is given as where φ  (, ) is the unit secondary warping function with respect to the shear center  of the cross section of the bar, while taking into account ( 73) and (20b) the aforementioned function φ can be established by solving independently the following Neumann problem: Moreover, applying the Green identity for the warping functions    (, ) and φ  (, ) and taking into account the Neumann problems (19a), (19b), (81a), (81b), and (80) the secondary torsion constant is given as where   is a domain integral given from the relation while form ( 75) and ( 82) the following relation for the warping shear deformation coefficient is obtained: 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 segment  of the bar (Figure 12) is considered leading to the relation Taking into account ( 29), ( 85) can be written as while substituting (40) into (86) the following equilibrium equation is obtained: Having in mind that both the secondary twisting moment    (   ) and the bimoment   (   ) (given from relations (71b) and (70), resp.) are expressed in terms of the primary angle of twist (   ), while the primary twisting moment    (  ) (given from relation (71a)) is expressed in terms of the total angle of twist (  ), an effort is given in the following to formulate also the primary twisting moment    with respect to the primary angle of twist (   ).Introducing an auxiliary geometric constant  defined as (88) the relation ( 76) can be written as while from relation (65b) the total torsional curvature can be written as and according to (71a) and (71b), the primary twisting moment is given with respect to the primary angle of twist (   ) as where   is a modified warping constant given as 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   is given with respect to the primary angle of twist (   ) from the relation while using (92) the moments    and   can be written in terms of the modified warping constant as 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 while the boundary conditions are given from where the parameters   ,   ( = 1, 2, 3) are functions specified at the boundary of the bar and the moments   and   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   and   (e.g., for a clamped edge it is Furthermore, substituting ( 95) into ( 89) and differentiating with respect to , 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 while the corresponding boundary conditions are given in terms of the primary angle of twist as where if the primary angle of twist is fixed, that is  1 = 1 and  2 =  3 = 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).

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 thinwalled 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 length  (Figure 13), of constant arbitrary cross-section of area .The homogeneous isotropic and linearly elastic material of the bar's cross-section, with modulus of elasticity , shear modulus , and Poisson's ratio ], occupies the two dimensional multiply connected region Ω of the ,  plane and is bounded by the Γ  ( = 0, 1, 2, . . ., ) boundary curves, which are piecewise smooth; that is, they may have a finite number of corners.In Figure 13(b)  is the principal coordinate system through the cross section's centroid , while   ,   are its coordinates with respect to  system of axes through the cross section's shear center .The bar is subjected to the combined action of the arbitrarily distributed or concentrated conservative twisting   =   () and warping   =   () moments acting in the  direction (Figure 13(a)).
Under the aforementioned loading the displacement field of the bar with respect to the  system of axes for large twisting rotations and small bending ones is given as where the transverse displacement components V,  are valid for large rotations [99];   ,   are the angles of rotation due to bending of the cross section with respect to its centroid;    () denotes the rate of change of the angle of twist   regarded as the torsional curvature; V  (),   () are the transverse displacement components of the shear center ;    ,    are the primary and secondary warping functions with respect to the shear center , respectively [21];   () 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 and the nonvanishing strain resultants are given as where the curvature components   ,   are given from the following relations: while the second-order geometrically nonlinear term in the right hand side of (103a) ( 2 +  2 ) ⋅ (   ) 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    , 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 or employing (103a), (103b), and (103c) as ignoring the independent of warping terms such as  ⋅    or warping terms due to shear such as (   ( −   ) + V   ( −   )) 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    given from (19a) and (19b) and for the secondary warping function    given from are obtained, where ∇ 2 =  2 / 2 +  2 / 2 is the Laplace operator and / denotes the directional derivative normal to the boundary Γ.Thus, the primary    and the secondary    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 Moreover, in the case the origin  of the coordinates is a point of the plane other than the shear center , the primary warping function with respect to this point    is first established from the Neumann problems (19a) and (19b) substituting    by    .Using the evaluated warping function    ,    is then established using the transformation presented in Section 2, while the coordinates of the shear center  with respect to the  system of coordinates are established from (109) and the conditions that define point  as the point for which the axial and bending stress resultants arising from an angle of twist at this ISRN Civil Engineering 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   (), V  (),   ().
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 leads to that is   generally does not coincide either with the axial displacement of the cross section's centroid   or with that of the shear center   .

Equations of Global Equilibrium. Neglecting body forces the principle of virtual work yields
where (⋅) denotes virtual quantities,   ,   ,   are the components of the traction vector with respect to the undeformed surface of the bar including the cross section ends,  is the surface area, and  is the volume of the bar at its initial undeformed state.Moreover, the stress resultants of the bar are given as where    is the primary twisting moment resulting from the primary shear stress distribution    ,    and   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 where the torsion constant   and the warping constant   with respect to the shear center  are given from (32a) and (32b), the polar moment of inertia   with respect to the shear center  and the moments of inertia   ,   with respect to the cross section's centroid are given as with  = √ 2 +  2 , while the geometric constants  1 ,  2 , and   are given as 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 force  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   is not necessarily the centroid's displacement   [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    with respect to the shear center  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 and the corresponding boundary conditions are written as where the "higher order torsion constant"   and the variable Ψ are given as while   ,   in the boundary conditions (120a) and (120b) are the externally applied conservative twisting and warping moments at the bar ends and   = ∫ Ω  4 Ω (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 the governing equation of the nonuniform nonlinear torsional problem is written as subjected to the following boundary conditions at the bar ends where  2 is defined as and   ,   are the twisting and warping moments at the boundary of the bar, respectively, given as while   ,   ( = 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   and   (e.g., for a clamped edge it is The solution of the boundary value problem described by ( 123), (124a), and (124b) for the evaluation of the angle of twist   assumes that the warping   , the torsion   , and the geometric   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    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   is established, the secondary warping function    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   , the transverse displacement components of the shear center V  (),   () can be established integrating twice the following relations: 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 V   ,    the "average" axial displacement of the cross ISRN Civil Engineering section of the bar   can be established after substitution of (116a) to (122a) as 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.

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][105][106] studied also the nonlinear flexuraltorsional-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 flexuraltorsional 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 length  of Figure 14.The homogeneous isotropic and linearly elastic material of the bar's cross-section of area , with modulus of elasticity , shear modulus , and Poisson's ratio ], occupies the two dimensional multiply connected region Ω of the ,  plane and is bounded by the Γ  ( = 0, 1, 2, . . ., ) boundary curves, which are piecewise smooth; that is, they may have a finite number of corners.In Figure 14(b)  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 (, ),  ≥ 0, twisting   =   (, ) and warping   =   (, ) moments, acting in the  direction (Figure 14(a)).
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 where under a Total Lagrangian formulation.In the above equations, (⋅) denotes virtual quantities, (⋅) denotes differentiation with respect to time, ,  are the volume and the surface of the bar, respectively, at the initial configuration, and   ,   ,   are the components of the traction vector with respect to the undeformed surface of the bar.
Neglecting virtual terms of the secondary warping function    , following the technique presented in [100] and after some algebra [112], the following local equilibrium equation along with its corresponding boundary condition 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 φ   , the following governing equation for the secondary    warping function is obtained as along with its corresponding boundary condition while as expected, the governing equation related to the primary    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.

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 subjected to the initial conditions ( ∈ (0, )) together with the boundary conditions at the bar ends  = 0, where ,   ,   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   (, ),   (, ) assumes that the warping   and the torsion   constants are already established; the evaluation of which presumes that the primary warping function    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    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    can be achieved by neglecting the axial inertia term  ⋅ ü  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   (, )) 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.

Reduced Problems for Special Cases of Axial Boundary
Conditions.For the case of axially immovable ends, it is easily proved that [112] which after subsequent integration yields where Ñ is a time-dependent tensile axial load induced by the geometrical nonlinearity given as For the case of constant along the bar axial load equations (( 139) and ( 140)) hold by setting Ñ = (, ), where (, ) 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 where the pertinent initial and boundary conditions are appropriately modified, while the boundary value problem related to the secondary warping function    (( 134) and ( 135)) is accordingly modified to It is worth here noting that   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 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  ⋅ ü  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.

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][115][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-offreedom 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]  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 length  (Figure 15) with an arbitrarily shaped doubly symmetric constant cross section, occupying the two dimensional multiply connected region Ω of the ,  plane bounded by the Γ  ( = 0, 1, 2, . . ., ) boundary curves, which are piecewise smooth; that is, they may have a finite number of corners.In Figure 15(b)  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 , G, and  0 , respectively.The bar is subjected to the combined action of arbitrarily distributed or concentrated twisting   =   () and warping   =   () moments acting in the  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  (, , ) =    ()    (, ) +    ()    (, ) , (145a) (, ) =   () , where , V,  are the axial and transverse bar displacement components with respect to the  system of axes;    () denotes the rate of change of the angle of twist   () regarded as the torsional curvature, while    () is the third order derivative of the angle of twist with respect to ;    ,    are the primary and secondary warping functions, respectively with respect to the shear center .
Substituting (145a), (145b), and (145c) in the well-known three-dimensional linear strain-displacement relations, the nonvanishing (total) strain resultants are obtained as 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 where (⋅) denotes infinitesimal incremental quantities over time (rates), the superscript el denotes the elastic part of the strain components and If the plane stress hypothesis is undertaken then  * = /(1 − V 2 ) holds [68], while  is frequently considered instead of  * ( * ≈ ) in beam formulations [68,136].This last consideration is followed here, while any other reasonable expression of  * could also be used without any difficulty.
As long as the material remains elastic or elastic unloading occurs, that is, The stress rates are given with respect to the total strain ones by combining ( 147) and (148).If plastic flow occurs then 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 where   is the yield stress of the material and  pl eq is the equivalent plastic strain, the rate of which is defined in [137] and is equal to  pl eq =  ( is the proportionality factor).Moreover, the plastic modulus ℎ is defined as ℎ =   / pl eq or   = ℎ and can be estimated from a tension test as ℎ =   /( −   ) (Figure 16).According to the associated flow rule the plastic strain rates are given as Using the aforementioned relation linking the yield stress rate and the proportionality factor, ( 147), ( 149)-( 151) and exploiting the plastic loading condition ( = 0), the stress rates-total strain rates relations are resolved as where [ el-pl ] is the elastoplastic constitutive matrix given as

ISRN Civil Engineering
By setting ℎ = 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.

Equations of Local
Equilibrium.The primary and secondary warping functions are resolved by exploiting local equilibrium considerations along the longitudinal  axis (together with the associated boundary condition at the lateral surface of the bar) under elastic conditions.Thus,    ,    are evaluated by solving the boundary value problems 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   ().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.

Equations of Global Equilibrium.
To establish global equilibrium equations, the principle of virtual work neglecting body forces is employed as follows: where (⋅) denotes virtual quantities,   ,   ,   are the components of the traction vector,  is the volume, and  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       ,       ,       , and       , 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   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 is obtained, where   ,   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 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    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 ) 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    , 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 Δ  , Δ  , and Δ depend on the third order derivative Δ   which will eventually lead to the requirement of an additional (third) boundary condition apart from the wellknown torsional and warping ones.In order to avoid the formulation of a third-order shear deformation theory, terms of Δ   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 Section 2) where   and   correspond to the internal twisting and warping moments of the bar, respectively, and the relevant incremental stress resultants are defined as 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 along with its corresponding boundary conditions where   ,   ( = 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   and   (e.g., for a clamped edge it is Finally, the expressions of the externally applied loading (  ,   ) with respect to the components of the traction vector (or the Cauchy stress tensor) may be easily resolved by virtue of (159) as 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.

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][115][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 , 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][144][145][146].Since the topic at hand is analogous to the geometrically nonlinear Timoshenko beam theory of shear-bending loading conditions [78,83], ISRN Civil Engineering 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 thickwalled 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 Figure 15 subjected to the combined action of arbitrarily distributed or concentrated twisting   =   () and warping   =   () moments acting in the  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  (, , ) = (   ())     (, ) , V (, ) = −  () , (, ) =   () , where , V,  are the axial and transverse bar displacement components with respect to the  system of axes;   is the (total) angle of twist; (   )  is the primary angle of twist per unit length which is in general not equal to the angle of twist per unit length    (see also (65a);    is the primary warping function with respect to the shear center S.
Substituting (165a), (165b), and (165c) in the wellknown three-dimensional linear strain-displacement relations, the nonvanishing compatible (total) strain resultants are obtained as where in (166b) and (166c), the first terms correspond to St. Venant shear strains (   ,    ) and the last terms to warping shear strains (   ,    ).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 where   is the torsional shear correction factor, while the specific form of the correction (√  ) is justified from the use of an energy approach to evaluate   [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    (, ) [151] is introduced as where    and   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 Section 9, 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: +    =   (169b) where   ,    , and    correspond to warping moment, primary, and secondary twisting moments, respectively.Equation (171c) gives the expression of    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)   ,   ( = 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 (  ,   ) with respect to the components of the traction vector (or the Cauchy stress tensor) may be easily resolved as Since an incremental-iterative approach is adopted for the problem at hand, the incremental version of (169a) and (169b) is written down as where Δ(⋅) denotes incremental quantities (over time), while the incremental stress resultants are given by virtue of (171a), (171b), (171c), (147), and (149) as 1 Δ  +  2 (Δ   )  = Δ 3 at the bar ends  = 0, .
In order to avoid the third order derivative of Δ   appearing in (177b), a one dimensional independent warping parameter   is set equal to the primary angle of twist per unit length, thus the boundary value problem of the above equations is reformulated as

− 𝐺 (𝐼
1 Δ  +  2 Δ  = Δ 3 at the bar ends  = 0, .(180b) 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.

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.

Figure 1 :
Figure 1: Prismatic bar subjected to a twisting moment (a) with a cross section of arbitrary shape occupying the two dimensional region Ω (b).
(b)), the warping function with respect to this point    is first established from the Neumann problem (19a) and (19b) substituting

Figure 4 :Figure 5 :
Figure 4: Distribution of primary (a) and secondary (b) shear stresses and resulting torsional moments.

Figure 8 :
Figure 8: Torsional curvature of a rectangular and a hollow square cross section.

Figure 9 :
Figure 9: Bar subjected to a time dependent twisting moment (a) with a variable cross section of arbitrary shape occupying the two dimensional region Ω (b).

Figure 10 :
Figure 10: Rectilinear bar subjected to twisting and warping moments (a) of constant doubly symmetric arbitrary cross section occupying the two dimensional region Ω (b).

Figure 11 :
Figure 11: Primary (a) and secondary (b) torsional curvature of a hot rolled I-profile steel HEB-100 cross section.

Figure 12 :
Figure 12: Equilibrium of a small segment  of a twisted bar.

Figure 13 :
Figure 13: Prismatic bar subjected to twisting and warping moments (a) with a cross section of arbitrary shape occupying the two dimensional region Ω (b).

Figure 14 :
Figure 14: Bar subjected to axial and torsional time dependent loading (a) of arbitrarily shaped doubly symmetric constant cross section occupying region Ω (b).

Figure 15 :
Figure 15: Prismatic bar subjected to twisting and warping moments (a) of an arbitrarily shaped doubly symmetric cross section, occupying region Ω (b).