Bars under Torsional Loading: A Generalized Beam Theory Approach

Abstract

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

1. Introduction

In engineering practice, we often come across the analysis of members of structures subjected to twisting moments. Curved bridges, ribbed plates subjected to eccentric loading, or columns laid out irregularly in the interior of a plate due to functional requirements are the most common examples.

When the warping of the cross section of a member is not restrained, the applied twisting moment is undertaken from the Saint-Venant [1] shear stresses. In this case the angle of twist per unit length remains constant along the bar. However, in most cases arbitrary torsional boundary conditions are applied either at the edges or at any other interior point of the bar due to construction requirements. This bar under the action of general twisting loading is leaded to nonuniform torsion, while the angle of twist per unit length is no longer constant along it. The consequences of restrained warping were first presented by Marguerre [2].

Although there is extended literature on the Saint-Venant uniform torsion problem for homogeneous isotropic cylindrical bars with simply or multiply connected cross sections [3–8], the extensive use of structural elements subjected to torsional loading necessitates a reliable, accurate, and general analysis of the torsion problem of bars of arbitrary cross section avoiding the restrictions of the Saint-Venant torsion theory.

In the rest of this paper, both the static and dynamic analyses of the geometrically linear or nonlinear, elastic or elastoplastic nonuniform torsion problems of bars of constant or variable arbitrary cross section are presented together with the corresponding research efforts and the conclusions drawn from examined cases with great practical interest.

2. Linear Elastic Nonuniform Torsion of Bars

In the last decades considerable work has been done on the elastic linear problem of nonuniform torsion of bars. Especially because of the mathematical complexity of the problem, the existing analytical solutions are limited to symmetric cross sections of simple geometry, loading, and boundary conditions [9–16]. Moreover, numerical methods such as the finite element method [17] or the boundary element method [18–20] have also been used for the analysis of the nonuniform torsional problem, in the case the geometry of the cross section, the boundary conditions, or the loading are not simple. In all the aforementioned references the analysis of the nonuniform torsional problem is not complete, since the secondary shear stresses due to warping are not evaluated. Finally, Sapountzakis and Mokos in [21, 22] developed a boundary element solution for the general linear elastic nonuniform torsion problem of homogeneous or composite prismatic bars of arbitrary cross section subjected to an arbitrarily distributed or concentrated twisting moment and supported by the most general linear torsional boundary conditions. In these latter research efforts the evaluation of the secondary warping function is accomplished leading to the computation of the secondary shear stresses due to warping, while the developed procedure is a pure BEM [23], since it requires only boundary discretization.

In order to formulate the aforementioned problem, let us consider a prismatic bar of 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.

(a)

(b)

Figure 1 

Prismatic bar subjected to a twisting moment (a) with a cross section of arbitrary shape occupying the two dimensional region Ω (b).

Under the aforementioned loading the displacement field of the bar with respect to the 𝑆𝑦𝑧 system of axes is given as𝑢𝑥,𝑦,𝑧=𝜃^′_𝑥(𝑥)𝜙_𝑆,𝑣𝑦,𝑧𝑥,𝑦,𝑧=−𝑧𝜃_𝑥(𝑥),𝑤𝑥,𝑦,𝑧=𝑦𝜃_𝑥(𝑥), (1a)(1b)(1c)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𝜀_𝑥𝑥=𝜕𝑢𝜕𝑥=𝜃_𝑥^′′(𝑥)⋅𝜙_𝑆,𝛾𝑦,𝑧_𝑥𝑦=𝜕𝑢+𝜕𝑦𝜕𝑣𝜕𝑥=𝜃^′_𝑥(𝑥)⋅𝜕𝜙_𝑆,𝛾𝜕𝑦−𝑧_𝑥𝑧=𝜕𝑢+𝜕𝑧𝜕𝑤𝜕𝑥=𝜃^′_𝑥(𝑥)⋅𝜕𝜙_𝑆𝜕𝑧+𝑦 (2a)(2b)(2c)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𝜎_𝑥𝑥=𝐸())𝜀1+𝜈)(1−2𝜈(1−𝜈_𝑥𝑥+𝜈𝜀_𝑦𝑦+𝜀_𝑧𝑧=𝐸⋅𝜃_𝑥^′′(𝑥)⋅𝜙_𝑆,𝜏𝑦,𝑧_𝑥𝑦=𝐺⋅𝛾_𝑥𝑦=𝐺⋅𝜃^′_𝑥(𝑥)⋅𝜕𝜙_𝑆,𝜏𝜕𝑦−𝑧_𝑥𝑧=𝐺⋅𝛾_𝑥𝑧=𝐺⋅𝜃^′_𝑥(𝑥)⋅𝜕𝜙_𝑆.𝜕𝑧+𝑦 (3a)(3b)(3c)

2.1. Equations of Local Equilibrium

Substituting the stress tensor components in the three equilibrium equations of the three-dimensional elasticity theory ignoring body forces yields𝜕𝜕𝑦𝐺⋅𝜃^′_𝑥(𝑥)⋅𝜕𝜙_𝑆+𝜕𝜕𝑦−𝑧𝜕𝑧𝐺⋅𝜃^′_𝑥(𝑥)𝜕𝜙_𝑆+𝜕𝜕𝑧+𝑦𝜕𝑥𝐸⋅𝜃_𝑥^′′(𝑥)⋅𝜙_𝑆=0,𝐺⋅𝜃_𝑥^′′(𝑥)⋅𝜕𝜙_𝑆𝜕𝑧+𝑦=0,𝐺⋅𝜃_𝑥^′′(𝑥)⋅𝜕𝜙_𝑆𝜕𝑦−𝑧=0. (4a)(4b)(4c)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 𝜕²𝜙_𝑆𝜕𝑦²+𝜕²𝜙_𝑆𝜕𝑧²=−𝐸⋅𝜃_𝑥^′′′(𝑥)𝐺⋅𝜃^′_𝑥(𝑥)⋅𝜙_𝑆. (5) 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 𝜎^𝑤_𝑥𝑥.

Figure 2 

Normal and shear stresses due to warping.

Based on the above decomposition of shear stress in primary and secondary components, that is𝜏_𝑥𝑦=𝜏^𝑃_𝑥𝑦+𝜏^𝑆_𝑥𝑦,𝜏_𝑥𝑧=𝜏^𝑃_𝑥𝑧+𝜏^𝑆_𝑥𝑧 (6a)(6b)the primary and secondary components of these stresses and the normal stresses due to warping are defined according to the relations as follows:𝜏^𝑃_𝑥𝑦=𝐺⋅𝜃^′_𝑥(𝑥)⋅𝜕𝜙^𝑃_𝑆,𝜏𝜕𝑦−𝑧_𝑆𝑥𝑦=𝐺⋅𝜕𝜙^𝑆_𝑆,𝜕𝑦 (7a)(7b)

𝜏^𝑃_𝑥𝑧=𝐺⋅𝜃^′_𝑥(𝑥)⋅𝜕𝜙^𝑃_𝑆,𝜏𝜕𝑧+𝑦_𝑆𝑥𝑧=𝐺⋅𝜕𝜙^𝑆_𝑆,𝜕𝑧 (8a)(8b)𝜎^𝑤_𝑥𝑥=𝐸⋅𝜃_𝑥^′′(𝑥)⋅𝜙^𝑃_𝑆𝑦,𝑧, (9) 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 𝜕𝜎^𝑤_𝑥𝑥+𝜕𝑥𝜕𝜏^𝑃_𝑥𝑦+𝜕𝑦𝜕𝜏^𝑃_𝑥𝑧+𝜕𝑧𝜕𝜏^𝑆_𝑥𝑦+𝜕𝑦𝜕𝜏^𝑆_𝑥𝑧𝜕𝑧=0. (10) 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 𝜕𝜏^𝑃_𝑥𝑦+𝜕𝑦𝜕𝜏^𝑃_𝑥𝑧𝜕𝑧=0,𝜕𝜏^𝑆_𝑥𝑦+𝜕𝑦𝜕𝜏^𝑆_𝑥𝑧+𝜕𝑧𝜕𝜎^𝑤_𝑥𝑥𝜕𝑥=0. (11) Substituting (7a)–(9) into (11) the following relations are obtained: ∇²𝜙^𝑃_𝑆=𝜕²𝜙^𝑃_𝑆𝜕𝑦²+𝜕²𝜙^𝑃_𝑆𝜕𝑧²∇=0,₂𝜙^𝑆_𝑆=𝜕²𝜙^𝑆_𝑆𝜕𝑦²+𝜕²𝜙^𝑆_𝑆𝜕𝑧²=−𝐸⋅𝜃_𝑥^′′′(𝑥)𝐺⋅𝜙^𝑃_𝑆. (12)

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: 𝜏_𝑥𝑛=𝜏_𝑥𝑦⋅𝑛_𝑦+𝜏_𝑥𝑧⋅𝑛_𝑧,𝜏_𝑥𝑡=−𝜏_𝑥𝑦⋅𝑛_𝑧+𝜏_𝑥𝑧⋅𝑛_𝑦, (13) where 𝑛_𝑦=cos(𝑦,𝑛)=cos𝛼=𝑑𝑦/𝑑𝑛=𝑑𝑧/𝑑𝑠 and 𝑛_𝑧=sin(𝑧,𝑛)=sin𝛼=𝑑𝑧/𝑑𝑛=−𝑑𝑦/𝑑𝑠 are the direction cosines of the vector 𝐧 normal to the boundary of the cross section. Substituting relations (6a) and (6b) that express the decomposition of shear stresses in primary and secondary ones into relations (13) yields𝜏^𝑃_𝑥𝑛=𝜏^𝑃_𝑥𝑦⋅𝑛_𝑦+𝜏^𝑃_𝑥𝑧⋅𝑛_𝑧,𝜏^𝑃_𝑥𝑡=−𝜏^𝑃_𝑥𝑦⋅𝑛_𝑧+𝜏^𝑃_𝑥𝑧⋅𝑛_𝑦, (14a)(14b)

Figure 3 

Shear stress components on the boundary of the cross section.

𝜏^𝑆_𝑥𝑛=𝜏^𝑆_𝑥𝑦⋅𝑛_𝑦+𝜏^𝑆_𝑥𝑧⋅𝑛_𝑧,𝜏^𝑆_𝑥𝑡=−𝜏^𝑆_𝑥𝑦⋅𝑛_𝑧+𝜏^𝑆_𝑥𝑧⋅𝑛_𝑦. (15a)(15b) Substituting the expressions of primary and secondary shear stresses given by (7a), (7b), (8a), and (8b) in (14a), (14b), (15a), and (15b) the following relations are obtained:𝜏^𝑃_𝑥𝑛=𝐺⋅𝜃^′_𝑥(𝑥)𝜕𝜙^𝑃_𝑆𝜕𝑛−𝑧⋅𝑛_𝑦+𝑦⋅𝑛_𝑧,𝜏^𝑃_𝑥𝑡=𝐺⋅𝜃^′_𝑥(𝑥)𝜕𝜙^𝑃_𝑆𝜕𝑡+𝑦⋅𝑛_𝑦+𝑧⋅𝑛_𝑧,𝜏^𝑆_𝑥𝑛=𝐺⋅𝜕𝜙^𝑆_𝑆,𝜏𝜕𝑛_𝑆𝑥𝑡=𝐺⋅𝜕𝜙^𝑆_𝑆.𝜕𝑡 (16a)(16b)(16c)(16d)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𝜏^𝑃_𝑥𝑛𝜏=0,_𝑆𝑥𝑛=0. (17a)(17b)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𝜕𝜙^𝑃_𝑆𝜕𝑛=𝑧⋅𝑛_𝑦−𝑦⋅𝑛_𝑧,𝜕𝜙^𝑆_𝑆𝜕𝑛=0. (18a)(18b)

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:∇²𝜙^𝑃_𝑆=𝜕²𝜙^𝑃_𝑆𝜕𝑦²+𝜕²𝜙^𝑃_𝑆𝜕𝑧²=0inthedomainofcrosssectionΩ,𝜕𝜙^𝑃_𝑆𝜕𝑛=𝑧⋅𝑛_𝑦−𝑦⋅𝑛_𝑧ontheboundaryofthecrosssectionΓ. (19a)(19b)For the Secondary 𝜙^𝑆_𝑆(𝑥,𝑦,𝑧) Warping Function. Neumann problem of the Poisson differential equation is as follows:∇²𝜙^𝑆_𝑆=𝜕²𝜙^𝑆_𝑆𝜕𝑦²+𝜕²𝜙^𝑆_𝑆𝜕𝑧²=−𝐸⋅𝜃_𝑥^′′′(𝑥)𝐺⋅𝜙^𝑃_𝑆inthedomainofcrosssectionΩ,𝜕𝜙^𝑆_𝑆𝜕𝑛=0ontheboundaryofthecrosssectionΓ. (20a)(20b)

Finally, based on the above mentioned, the displacement field of the bar ((1a), (1b), and (1c)) is formed as𝑢𝑥,𝑦,𝑧=𝜃^′_𝑥(𝑥)⋅𝜙^𝑃_𝑆𝑦,𝑧+𝜙^𝑆_𝑆,𝑣𝑥,𝑦,𝑧𝑥,𝑦,𝑧=−𝑧𝜃_𝑥(𝑥),𝑤𝑥,𝑦,𝑧=𝑦𝜃_𝑥(𝑥). (21a)(21b)(21c)

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(b)), the warping function with respect to this point 𝜑^𝑃_𝑂 is first established from the Neumann problem (19a) and (19b) substituting 𝜑^𝑃_𝑆 by 𝜑^𝑃_𝑂. Using the evaluated warping function 𝜑^𝑃_𝑂, 𝜑^𝑃_𝑆 is then established using the transformation given by the following equation [2]: 𝜙^𝑃_𝑆𝑦,𝑧=𝜙^𝑃_𝑂𝑦,𝑧−𝑦𝑧^𝑆+𝑧𝑦^𝑆+𝑐, (22) 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:𝑆_𝑦𝑦^𝑆−𝑆_𝑧𝑧^𝑆+𝐴𝑐=−𝑅^𝑃_𝑆,𝐼_𝑦𝑦𝑦^𝑆+𝐼_𝑦𝑧𝑧^𝑆+𝑆_𝑦𝑐=−𝑅^𝑃_𝑦,𝐼_𝑦𝑧𝑦^𝑆+𝐼_𝑧𝑧𝑧^𝑆−𝑆_𝑧𝑐=𝑅^𝑃_𝑧, (23a)(23b)(23c)where𝐴=_Ω𝑆𝑑Ω,_𝑦=_Ω𝑆𝑧𝑑Ω,_𝑧=_Ω𝐼𝑦𝑑Ω,_𝑦𝑦=_Ω𝑧²𝐼𝑑Ω,_𝑧𝑧=_Ω𝑦²𝐼𝑑Ω,_𝑦𝑧=−_Ω𝑦𝑧𝑑Ω, (24a)(24b)(24c)(24d)(24e)(24f)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𝑅^𝑃_𝑆=_Ω𝜙^𝑃_𝑂𝑅𝑑Ω,_𝑃𝑦=_Ω𝑧𝜙^𝑃_𝑂𝑅𝑑Ω,_𝑃𝑧=_Ω𝑦𝜙^𝑃_𝑂𝑑Ω. (25a)(25b)(25c)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: 𝑐^𝑆=1𝐴_Ω𝜙^𝑆_𝑆𝑑Ω (26) and the main secondary warping function 𝜙^𝑆_𝑆 is given as 𝜙^𝑆_𝑆=𝜙^𝑆_𝑆+𝑐^𝑆. (27)

2.2. 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 𝑀_𝑡=_Ω𝜏_𝑥𝑧⋅𝑦−𝜏_𝑥𝑦⋅𝑧𝑑Ω. (28) 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: 𝑀_𝑡=𝑀^𝑃_𝑡+𝑀^𝑆_𝑡, (29) where after employing (6a) and (6b) and some algebra the aforementioned twisting moment components are given as𝑀^𝑃_𝑡=_Ω𝜏^𝑃_𝑥𝑦𝜕𝜙^𝑃_𝑆𝜕𝑦−𝑧+𝜏^𝑃_𝑥𝑧⋅𝜕𝜙^𝑃_𝑆𝑀𝜕𝑧+𝑦𝑑Ω,_𝑆𝑡=_Ω−𝜏^𝑆_𝑥𝑦𝜕𝜙^𝑃_𝑆𝜕𝑦−𝜏^𝑆_𝑥𝑧𝜕𝜙^𝑃_𝑆𝜕𝑧𝑑Ω. (30a)(30b)Substituting (7a)–(8b) in (30a) and (30b), employing the Gauss-Green theorem, taking into account (17a), (17b), (19a), and (20a) and after some algebra the following relations for the twisting moment components are obtained:𝑀^𝑃_𝑡=−𝐺𝐼_𝑡𝜃^′_𝑥(𝑥),𝑀^𝑆_𝑡=−𝐸𝐶_𝑆𝜃_𝑥^′′′(𝑥), (31a)(31b)where𝐼_𝑡=_Ω𝑦²+𝑧²+𝑦⋅𝜕𝜙^𝑃_𝑆𝜕𝑧−𝑧⋅𝜕𝜙^𝑃_𝑆𝐶𝜕𝑦𝑑Ω,_𝑆=_Ω𝜙^𝑃_𝑆²𝑑Ω. (32a)(32b)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 𝑆.

(a)

(b)

Figure 4 

Distribution of primary (a) and secondary (b) shear stresses and resulting torsional moments.

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 𝑀_𝑡+𝜕𝑀_𝑡𝜕𝑥𝑑𝑥−𝑀_𝑡+𝑚_𝑡𝑑𝑥=0 (33) or after some algebra 𝜕𝑀_𝑡𝜕𝑥=−𝑚_𝑡. (34) Substituting relations (29) and (31a) and (31b) into (34), the following fourth-order differential equation of equilibrium of a homogeneous isotropic bar subjected to nonuniform torsion is obtained: 𝐸𝐶_𝑆𝑑⁴𝜃_𝑥𝑑𝑥⁴−𝐺𝐼_𝑡𝑑²𝜃_𝑥𝑑𝑥²=𝑚_𝑡 (35) which is independent of the torsional boundary conditions.

Figure 5 

Stress resultants in an infinitesimal segment of an elastic bar subjected to an arbitrary concentrated or distributed torsional loading 𝑚𝑡=𝑚𝑡(𝑥).

In analogy with the bending moments 𝑀_𝑦, 𝑀_𝑧 a new stress resultant is defined, called warping moment (or bimoment) and given by 𝑀_𝑤=−_Ω𝜙^𝑃_𝑆𝜎^𝑤_𝑥𝑥𝑑Ω. (36) 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 𝑀_𝑤=−𝐸𝜃_𝑥^′′_Ω𝜙^𝑃_𝑆²𝑑Ω (37) or through (32b) 𝑀_𝑤=−𝐸𝐶_𝑆𝜃_𝑥^′′ (38) and therefore relation (9) can be written in the form of 𝜎^𝑤_𝑥𝑥𝑀=−^𝑤𝐶_𝑆𝜑^𝑃_𝑆. (39) Equation (39) shows that the quantity of warping moment encompasses the generic characteristics of a stress resultant of theory of elasticity. Schardt [25] refers to it as a “higher order stress resultant.” Combining relations (31b) and (38), the relation which correlates the secondary twisting moment and the warping moment arises as 𝑑𝑀_𝑤𝑑𝑥=𝑀^𝑆_𝑡. (40) 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𝑎₁𝜃_𝑥+𝑎₂𝑀_𝑡=𝑎₃,𝛽₁𝜃^′_𝑥+𝛽₂𝑀_𝑤=𝛽₃. (41a)(41b)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 𝑎₁=𝛽₁=1,  𝑎₂=𝑎₃=𝛽₂=𝛽₃=0.

From the examined example problems presented in [18–22] it is concluded that (i)the magnitude of the evaluated normal stresses, due to restrained warping compared with those due to bending, necessitates the consideration of these additional normal stresses near the restrained edges;(ii)the magnitude of the evaluated warping shear stresses due to restrained warping is remarkable and for an “accurate” analysis these additional shear stresses should not be ignored, especially near the restrained edges.

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

Long span box shaped bridges or concrete slab and beam structures of variable height are the most common examples of structures including members of variable cross section subjected to twisting moments. The extensive use of the aforementioned structural elements necessitates a rigorous analysis. When a bar of variable cross section is subjected to general twisting loading, this member due to this variation and/or due to the arbitrary torsional boundary conditions applied either at the edges or at any other interior point is leaded to nonuniform torsion and its angle of twist per unit length is not constant along its axis.

Several researchers have dealt with beams of variable cross section ignoring the warping effects resulting from the corresponding restraints at the ends of the member [26, 27]. If the aforementioned structures are analyzed or designed for torsion considering only the effect of Saint Venant torsion resistance, the analysis may underestimate the torsion in the members and the design may be unconservative. On the contrary, to the authors’ knowledge relatively little work has been done on the problem of nonuniform torsion of bars of variable cross section with pioneer the work of Cywinski [28] adopting the finite difference method. Wekezer [29] after dividing the bar into segments along its longitudinal axis approximated their shell midsurface by arbitrary triangular shell elements and employed the finite element method to the linear membrane shell theory. This approximation generates inaccuracies, as the warping of the walls of the cross section cannot be taken into account. Moreover, Eisenberger [30] employed FEM upon polynomial approximation of the torsional and warping rigidities using “exact” shape functions to derive the exact stiffness coefficients. This application of shape functions leads also to inaccuracies in stress analysis of beams of variable cross section, as static and kinematic values at nodes and in the element region are computed only approximately and the element may not satisfy local and global equilibrium conditions [31]. In all of the aforementioned procedures the torsion and warping constants have been approximated adopting the thin tube theory. Finally, Sapountzakis and Mokos in [32, 33] developed a boundary element solution for the general linear elastic nonuniform torsion problem of homogeneous or composite bars of arbitrary variable cross section subjected to an arbitrarily distributed or concentrated twisting moment and supported by the most general linear torsional boundary conditions. In these latter research efforts, three boundary value problems with respect to the variable along the beam angle of twist and to the primary and secondary warping functions are formulated and solved employing a pure BEM [23] approach; that is, only boundary discretization is used. Both the variable warping and torsion constants together with the torsional primary shear stresses and the warping normal and secondary shear stresses are computed.

In order to formulate the aforementioned problem, let us consider a bar of length 𝑙 (Figure 6), of an arbitrarily shaped variable along its axis cross section. The homogeneous isotropic and linearly elastic material of the bar’s cross-section, with modulus of elasticity 𝐸, 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)).

(a)

(b)

Figure 6 

Bar subjected to a twisting moment (a) with a variable cross section of arbitrary shape occupying the two dimensional region Ω(𝑥) (b).

Adopting the same displacement field (21a), (21b), and (21c) with that of the constant cross section case and following the same procedure presented for the constant cross section bar in Section 2, the following boundary value problem for the angle of twist 𝜃_𝑥=𝜃_𝑥(𝑥) is derived: 𝐸𝐶_𝑆𝑑⁴𝜃_𝑥𝑑𝑥⁴+2𝐸𝑑𝐶_𝑆𝑑𝑑𝑥³𝜃_𝑥𝑑𝑥³+𝐸𝑑²𝐶_𝑆𝑑𝑥²−𝐺𝐼_𝑡𝑑²𝜃_𝑥𝑑𝑥²−𝐺𝑑𝐼_𝑡𝑑𝑥𝑑𝜃_𝑥𝑑𝑥=𝑚_𝑡insidethebeam, (42)𝛼₁𝜃_𝑥+𝛼₂𝑀_𝑡=𝛼₃𝛽atthebeamends𝑥=0,𝑙,₁𝑑𝜃_𝑥𝑑𝑥+𝛽₂𝑀_𝑤=𝛽₃atthebeamends𝑥=0,𝑙, (43a)(43b)where the total twisting moment of the cross section is once again divided into a primary and a secondary component as this is stated in (29), where these components are given from𝑀^𝑃_𝑡=𝐺𝐼_𝑡𝑑𝜃_𝑥,𝑀𝑑𝑥_𝑆𝑡𝑑=−𝐸𝐶𝑑𝑥^𝑆𝑑²𝜃_𝑥𝑑𝑥² (44a)(44b)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 𝑀_𝑡=𝐺𝐼_𝑡𝑑𝜃_𝑥𝑑𝑥−𝐸𝑑𝐶_𝑆𝑑𝑑𝑥²𝜃_𝑥𝑑𝑥²−𝐸𝐶_𝑆𝑑³𝜃_𝑥𝑑𝑥³. (45) 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.

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

As it has been already mentioned, the analysis of rectilinear or curved members of structures of arbitrary constant or variable cross section subjected to twisting moments is often encountered in engineering practice. Nevertheless, accurate analysis of these members is difficult to achieve for two reasons.

According to the first reason, generally commercial programs consider six degrees of freedom at each node of a member of a space frame, ignoring in this way the warping effects due to the corresponding restraint at the ends of the member or due to variable twisting moment along the bar [26, 27, 34]. If the aforementioned structures are analyzed or designed for torsion considering only the effect of Saint Venant torsion resistance, the analysis may underestimate the torsion in the members and the design may be unconservative. Several researchers tried to overcome this inaccuracy only for constant cross section elements by developing a 14×14 member stiffness matrix including warping degrees of freedom at the ends of a member with open thin-walled section and assuming simple [35–38] or more complicated torsional boundary conditions [39, 40].

According to the second reason, when the variable cross section structures are modeled by space frame elements, they are usually analyzed using Hermite beam elements characterized by shape functions consisting of Hermite interpolation functions (cubic for bending, linear for axial components of displacement and for the angle of twist) [41] or isoparametric beam elements [42]. In these elements the variation of the cross section is generally considered by setting “average” values for the cross section parameters resulting from the corresponding parameters at the start and end nodes of the element. Moreover, application of shape functions results in inaccuracies in stress analysis of beams of variable cross section, as static and kinematic values at nodes and in the element region are computed only approximately and the element may not satisfy local and global equilibrium conditions [31]. This inaccuracy may be reduced by increasing the number of integration points for the stiffness matrix assembly (at least three-point integration) or by refining the mesh of elements, increasing at the same time the effort of preparing the input parameters.

Nevertheless, Sapountzakis and Mokos in [43–45] developed a boundary element solution for the construction of the 14×14 stiffness matrix and the nodal load vector of a member of arbitrary homogeneous or composite, constant or variable cross section, subjected to an arbitrarily concentrated or distributed twisting moment, taking into account warping effects. The arbitrary low rate variation of the member of variable cross section is continuous so as to assume that its shear center is independent of the member loading and boundary conditions. The developed method can take into account the variable torsional and warping rigidities along the member length. In these latter research efforts boundary value problems with respect to the variable along the bar angle of twist and to the primary warping function are formulated and solved employing a pure BEM [23] approach; that is only boundary discretization is used.

In order to formulate the aforementioned problem, let us consider a prismatic bar of 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 𝑆.

(a)

(b)

Figure 7 

Beam element (a) of a cross section of arbitrary shape occupying the two dimensional region Ω (b). 𝐶𝑦𝑧 is the beam principal centroidal system of axes, while 𝑆𝑦𝑧 is the corresponding one through the 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 𝐷^𝑖𝑇=𝑢_𝑥𝑗𝑢_𝑦𝑗𝑢_𝑧𝑗𝜃_𝑥𝑗𝜃_𝑦𝑗𝜃_𝑧𝑗𝜃^′_𝑥𝑗𝑢_𝑥𝑘𝑢_𝑦𝑘𝑢_𝑧𝑘𝜃_𝑥𝑘𝜃_𝑦𝑘𝜃_𝑧𝑘𝜃^′_𝑥𝑘 (46)and the respective nodal load vector as𝐹^𝑖𝑇=𝑁_𝑗𝑄_𝑦𝑗𝑄_𝑧𝑗𝑀_𝑡𝑗𝑀_𝑦𝑗𝑀_𝑧𝑗𝑀_𝑤𝑗𝑁_𝑘𝑄_𝑦𝑘𝑄_𝑧𝑘𝑀_𝑡𝑘𝑀_𝑦𝑘𝑀_𝑧𝑘𝑀_𝑤𝑘, (47)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,𝑀_𝑡=𝐺𝐼_𝑡𝜃^′_𝑥−𝐸𝐶_𝑆𝜃_𝑥^′′′𝑀fortheconstantcrosssectioncase,_𝑡=𝐺𝐼_𝑡𝜃^′_𝑥−𝐸𝐶^′_𝑆𝜃_𝑥^′′−𝐸𝐶_𝑆𝜃_𝑥^′′′forthevariablecrosssectioncase. (48a)(48b)Moreover, 𝑀_𝑤 is the bending moment due to the torsional curvature at the same sections given from (38).

Figure 8 

Torsional curvature of a rectangular and a hollow square cross section.

The nodal displacement and load vectors given in (46) and (47) are related with the 14×14 local stiffness matrix of the spatial beam element written as𝑘^𝑖=𝑘^𝑖₁₁000000𝑘^𝑖₁₈0000000𝑘^𝑖₂₂000𝑘^𝑖₂₆00𝑘^𝑖₂₉000𝑘^𝑖_2,13000𝑘^𝑖₃₃0𝑘^𝑖₃₅0000𝑘^𝑖_3,100𝑘^𝑖_3,1200000𝐤^𝐢_𝐓𝟏00𝐤^𝐢_𝐓𝟐000−𝐤^𝐢_𝐓𝟏00𝐤^𝐢_𝐓𝟓00𝑘^𝑖₅₃0𝑘^𝑖₅₅0000𝑘^𝑖_5,100𝑘^𝑖_5,12000𝑘^𝑖₆₂000𝑘^𝑖₆₆00𝑘^𝑖₆₉000𝑘^𝑖_6,130000𝐤^𝐢_𝐓𝟐00𝐤^𝐢_𝐓𝟑000−𝐤^𝐢_𝐓𝟐00𝐤^𝐢_𝐓𝟒𝑘^𝑖₈₁000000𝑘^𝑖₈₈0000000𝑘^𝑖₉₂000𝑘^𝑖₉₆00𝑘^𝑖₉₉000𝑘^𝑖_9,13000𝑘^𝑖_10,30𝑘^𝑖_10,50000𝑘^𝑖_10,100𝑘^𝑖_10,1200000−𝐤^𝐢_𝐓𝟏00−𝐤^𝐢_𝐓𝟐000𝐤^𝐢_𝐓𝟓00−𝐤^𝐢_𝐓𝟓00𝑘^𝑖_12,30𝑘^𝑖_12,50000𝑘^𝑖_12,100𝑘^𝑖_12,12000𝑘^𝑖_13,2000𝑘^𝑖_13,600𝑘^𝑖_13,9000𝑘^𝑖_13,130000𝐤^𝐢_𝐓𝟓00𝐤^𝐢_𝐓𝟒000−𝐤^𝐢_𝐓𝟓00𝐤^𝐢_𝐓𝟔, (49)where the 𝑘^𝑖_𝑙𝑚 coefficients (𝑙,𝑚=1,2,3,5,6,8,9,10,12,13) of the stiffness matrix of (49) come from the well-known 12×12 classical stiffness matrix of the classical three-dimensional frame element. In the special case of a constant cross section element the 𝑘^𝑖_𝑇𝑛(𝑛=1,2,3,4,5,6) coefficients are given as [14, 35]𝑘^𝑖_𝑇1=𝑘_𝑜𝑘𝜆𝑆,_𝑖𝑇2=𝑘^𝑖_𝑇5=𝑘_𝑜(),𝑘𝐶−1_𝑖𝑇3=𝑘^𝑖_𝑇6=𝑘_𝑜1𝐶𝑙−𝜆𝑆,𝑘^𝑖_𝑇4=𝑘_𝑜1𝜆𝑆−𝑙, (50a)(50b)(50c)(50d)where𝜆=𝐺𝐼_𝑡𝐸𝐶_𝑆,𝑘𝑆=sinh𝑒,𝐶=cosh𝑒,𝑒=𝑙𝜆,_𝑜=𝐺𝐼_𝑡2().1−𝐶+𝑒𝑆 (51a)(51b)(51c)(51d)(51e)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)): 𝐸𝐶_𝑆𝜃_𝑥^′′′′+2𝐸𝐶^′_𝑆𝜃_𝑥^′′′+𝐸𝐶_𝑆^′′−𝐺𝐼_𝑡𝜃_𝑥^′′−𝐺𝐼^′_𝑡𝜃^′_𝑥=0insidetheelement, (52)𝛼₁𝜃_𝑥+𝛼₂𝑀_𝑡=𝛼₃,𝛽₁𝜃^′_𝑥+𝛽₂𝑀_𝑤=𝛽₃attheelementends𝑥=0,𝑙 (53a)(53b)for appropriate values of the 𝑎_𝑖, 𝛽_𝑖(𝑖=1,2,3) functions. Thereby, for the evaluation of the 𝑘^𝑖_𝑇1 coefficient it is 𝑎₁=𝑎₃=𝛽₁=1,  𝑎₂=𝛽₂=𝛽₃=0 at 𝑥=0 and 𝑎₁=𝛽₁=1, 𝑎₂=𝑎₃=𝛽₂=𝛽₃=0 at 𝑥=𝑙, for the evaluation of the 𝑘^𝑖_𝑇3 coefficient it is 𝑎₁=𝛽₁=𝛽₃=1, 𝑎₂=𝑎₃=𝛽₂=0 at 𝑥=0 and 𝑎₁=𝛽₁=1, 𝑎₂=𝑎₃=𝛽₂=𝛽₃=0 at 𝑥=𝑙, while for the evaluation of the 𝑘^𝑖_𝑇6 coefficient it is 𝑎₁=𝛽₁=1, 𝑎₂=𝑎₃=𝛽₂=𝛽₃=0 at 𝑥=0 and 𝑎₁=𝛽₁=𝛽₃=1, 𝑎₂=𝑎₃=𝛽₂=0 at 𝑥=𝑙. 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: 𝐸𝐶_𝑆𝜃_𝑥^′′′′−𝐺𝐼_𝑡𝜃_𝑥^′′=𝑚_𝑡𝜃insidethebar,_𝑥=𝜃^′_𝑥=0atthebarends𝑥=0,𝑙, (54) for the variable cross section case: 𝐸𝐶_𝑆𝜃_𝑥^′′′′+2𝐸𝐶^′_𝑆𝜃_𝑥^′′′+𝐸𝐶_𝑆^′′−𝐺𝐼_𝑡𝜃_𝑥^′′−𝐺𝐼^′_𝑡𝜃^′_𝑥=𝑚_𝑡𝜃insidethebar,_𝑥=𝜃^′_𝑥=0atthebarends𝑥=0,𝑙. (55)

From the examined example problems presented in [43–45] it is concluded that (i)the discrepancy of the deflections and the internal stress resultants arising from the ignorance of the warping degrees of freedom at the ends of a member and the magnitude of the normal stresses due to warping necessitate the utilization of 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.

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

In engineering practice, we often come across the analysis of structures subjected to vibratory twisting loading. The dynamic forces acting on a structure may result from one or more of different causes, such as rotating machinery, wind, asymmetric traffic loading, blast loads, or earthquake forces. The extensive use of the aforementioned structural elements necessitates a rigorous dynamic analysis.

Exact torsional vibration frequencies were presented by Gorman [46] and Belvins [47] for the case of circular cross section shafts subjected to classical boundary conditions, avoiding in this way warping effects. These efforts were extended by Kameswara Rao [48] for elastically restrained edges. Torsional vibration frequencies for beams of open thin-walled sections, subjected to several combinations of classical boundary conditions, taking into account warping effects were first derived by Gere [49]. Since then approximate methods [50] for the calculation of natural frequencies including elastic torsional and warping restraints [51] and employing either discrete [52–55] or distributed mass model systems [56–60] have been presented. In all these references the considered beam is of a constant homogeneous thin-walled cross section, while its torsion and warping constants are evaluated employing the relations of the thin tube theory.

Several researchers have also dealt with beams of variable cross section ignoring the warping effects resulting from the corresponding restraints at the ends of the member [26, 27]. On the contrary, to the author’s knowledge relatively little work has been done on the problem of nonuniform torsion of bars of variable cross section with pioneer the work of Cywinski [28] adopting the finite difference method. Wekezer [29] after dividing the bar into segments along its longitudinal axis approximated their shell midsurface by arbitrary triangular shell elements and employed the finite element method to the linear membrane shell theory. This approximation generates inaccuracies, as the warping of the walls of the cross section cannot be taken into account. Moreover, Eisenberger [30, 61] employed FEM upon polynomial approximation of the torsional and warping rigidities using “exact” shape functions to derive the exact stiffness coefficients. This application of shape functions results also in inaccuracies in stress analysis of beams of variable cross section, as static and kinematic values at nodes and in the element region are computed only approximately and the element may not satisfy local and global equilibrium conditions [31]. In all the aforementioned procedures the torsion and warping constants have been approximated adopting the thin tube theory. Moreover, research efforts have been presented for the corresponding problem of composite beams limited to the formulation of a displacement-based one-dimensional finite element model for the estimation of natural frequencies and corresponding mode shapes of thin-walled composite beams after approximating again the torsion and warping constants with closed form solutions [62, 63].

Only in Ganapathi et al. [64] the warping function for the constant rectangular cross section of sandwich beams is determined by solving the boundary value problem for torsion such that the displacements are continuous at the interfaces of adjacent layers, while the transverse shear stress is continuous at these interfaces and vanishes at the top and bottom surfaces of the beam. Moreover, Sapountzakis in [65, 66] developed a boundary element method for the nonuniform torsional vibration problem of doubly symmetric composite bars of arbitrary constant or variable cross section, respectively. In these latter efforts the beam is subjected to an arbitrarily distributed dynamic twisting moment, while its edges are restrained by the most general linear torsional boundary conditions. A distributed mass model system is employed which leads to the formulation of three boundary value problems with respect to the variable along the beam angle of twist and to the primary and secondary warping functions. The last two problems are solved employing a pure BEM [23] approach that is only boundary discretization is used. Finally, Sapountzakis and Mokos in [67] presented the dynamic analysis of 3D beam elements restrained at their edges by the most general linear torsional, transverse, or longitudinal boundary conditions and subjected in arbitrarily distributed dynamic twisting, bending, transverse, or longitudinal loading. For the solution of this problem, a boundary element method is developed for the construction of the 14×14 stiffness matrix and the corresponding nodal load vector, of a member of an arbitrarily shaped simply or multiply connected cross section, taking into account both warping and shear deformation effects, which together with the respective mass and damping matrices lead to the formulation of the equation of motion.

In order to formulate the nonuniform torsional vibration problem of doubly symmetric bars of arbitrarily shaped simply or multiply connected constant or variable cross section, let us consider the bar of 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)).

(a)

(b)

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).

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)(56b)(56c)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: 𝐸𝐶_𝑆𝑑⁴𝜃_𝑥𝑑𝑥⁴+2𝐸𝑑𝐶_𝑆𝑑𝑑𝑥³𝜃_𝑥𝑑𝑥³+𝐸𝑑²𝐶_𝑆𝑑𝑥²−𝐺𝐼_𝑡𝑑²𝜃_𝑥𝑑𝑥²−𝐺𝑑𝐼_𝑡𝑑𝑥𝑑𝜃_𝑥𝑑𝑥+𝜌𝐼_𝑝𝑑²𝜃_𝑥𝑑𝑡²=𝑚_𝑡insidethebeam, (57)𝛼₁𝜃_𝑥+𝛼₂𝑀_𝑡=𝛼₃,𝛽₁𝑑𝜃_𝑥𝑑𝑥+𝛽₂𝑀_𝑤=𝛽₃atthebeamends𝑥=0,𝑙, (58a)(58b)

𝜃_𝑥(𝑥,0)=𝜗(𝑥),𝜃̇_𝑥(𝑥,0)=𝜗̇(𝑥), (59a)(59b) 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 𝐸𝐶_𝑆𝑑⁴𝜃_𝑥𝑑𝑥⁴−𝐺𝐼_𝑡𝑑²𝜃_𝑥𝑑𝑥²+𝜌𝐼_𝑝𝑑²𝜃_𝑥𝑑𝑡²=𝑚_𝑡insidethebeam. (60) 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𝜏^𝑃_𝑥𝑦𝑥,𝑦,𝑧,𝑡=𝐺⋅𝜃^′_𝑥(𝑥,𝑡)⋅𝜕𝜙^𝑃_𝑆,𝜏𝜕𝑦−𝑧_𝑆𝑥𝑦𝑥,𝑦,𝑧,𝑡=𝐺⋅𝜕𝜙^𝑆_𝑆,𝜕𝑦 (61a)(61b)

𝜏^𝑃_𝑥𝑧𝑥,𝑦,𝑧,𝑡=𝐺⋅𝜃^′_𝑥(𝑥,𝑡)⋅𝜕𝜙^𝑃_𝑆,𝜏𝜕𝑧+𝑦_𝑆𝑥𝑧𝑥,𝑦,𝑧,𝑡=𝐺⋅𝜕𝜙^𝑆_𝑆,𝜕𝑧 (62a)(62b)𝜎^𝑤_𝑥𝑥𝑥,𝑦,𝑧,𝑡=𝐸⋅𝜃_𝑥^′′(𝑥,𝑡)⋅𝜙^𝑃_𝑆𝑦,𝑧. (63)

From the examined example problems presented in [65–67] it is concluded that (i)the discrepancy of the dynamic response of the bar arising from the ignorance of the warping effect necessitates the inclusion of the nonuniform torsional beam behavior;(ii)the discrepancy of the results arising from the ignorance of the warping especially in higher eigenfrequencies is remarkable;(iii)the discrepancy in the analysis of a thin-walled cross section beam employing the BEM after calculating the torsion and warping constants adopting the thin tube theory demonstrates the importance of the proposed procedure even in thin-walled beams, since it approximates better the torsion and warping constants and takes also into account the warping of the walls of the cross section;(iv)the discrepancy of the results arising from the ignorance of the warping degrees of freedom at the ends of a member necessitates the utilization of 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.

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

As it has already been mentioned in the previous sections in engineering practice we often come across the analysis of rectilinear members of structures subjected to nonuniform torsion. In this case, as it has been analyzed in 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 Γ₀,Γ₁,…,Γ_𝐾 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.

(a)

(b)

Figure 10 

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

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𝑢_𝑦=𝑢^𝑃_𝑦+𝑢^𝑆_𝑦,𝑢_𝑧=𝑢^𝑃_𝑧+𝑢^𝑆_𝑧,𝑢_𝑥=𝑢^𝑃_𝑥+𝑢^𝑆_𝑥, (64a)(64b)(64c)

𝜃_𝑥=𝜃^𝑃_𝑥+𝜃^𝑆_𝑥,𝑑𝜃_𝑥=𝑑𝑥𝑑𝜃^𝑃_𝑥+𝑑𝑥𝑑𝜃^𝑆_𝑥,𝑑𝑥 (65a)(65b) 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.

(a)

(b)

Figure 11 

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

Assuming small torsional rotation, the displacement field is defined as𝑢^𝑃_𝑦=−𝑧𝜃_𝑥(𝑥),𝑢^𝑃_𝑧=𝑦𝜃_𝑥(𝑥),𝑢^𝑃_𝑥=𝑑𝜃^𝑃_𝑥(𝑥)𝜑𝑑𝑥^𝑃𝑆,𝑢𝑦,𝑧_𝑆𝑦=𝑢^𝑆_𝑧=𝑢^𝑆_𝑥=0, (66a)(66b)(66c)(66d)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 three-dimensional elasticity the nonzero stress components in the region Ω can be written as𝜎^𝑃_𝑥𝑥𝑑=𝐸²𝜃^𝑃_𝑥(𝑥)𝑑𝑥²𝜑^𝑃_𝑆,𝜏𝑦,𝑧_𝑃𝑥𝑦=𝐺𝑑𝜃_𝑥(𝑥)𝑑𝑥𝜕𝜑^𝑃_𝑆𝑦,𝑧,𝜏𝜕𝑦−𝑧_𝑃𝑥𝑧=𝐺𝑑𝜃_𝑥(𝑥)𝑑𝑥𝜕𝜑^𝑃_𝑆𝑦,𝑧,𝜕𝑧+𝑦 (67a)(67b)(67c)

𝜏^𝑆_𝑥𝑦=𝐺𝜕𝜑^𝑆_𝑆𝑥,𝑦,𝑧,𝜏𝜕𝑦_𝑆𝑥𝑧=𝐺𝜕𝜑^𝑆_𝑆𝑥,𝑦,𝑧,𝜕𝑧 (68a)(68b) 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.

6.1. Equations of Local Equilibrium

The first elasticity equation of equilibrium of the three-dimensional elasticity (equilibrium in the axial direction of the bar) with vanishing body forces is written as 𝜕𝜏^𝑃_𝑥𝑦+𝜕𝑦𝜕𝜏^𝑃_𝑥𝑧𝜕𝑧^𝑃+𝜕𝜎^𝑃_𝑥𝑥+𝜕𝑥𝜕𝜏^𝑆_𝑥𝑦+𝜕𝑦𝜕𝜏^𝑆_𝑥𝑧𝜕𝑧^𝑆=0. (69) 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: 𝑀_𝑤=−𝐸𝐶_𝑆𝑑²𝜃^𝑃_𝑥𝑑𝑥², (70)𝑀^𝑃_𝑡=𝐺𝐼^𝑃_𝑡𝑑𝜃_𝑥,𝑀𝑑𝑥_𝑆𝑡=−𝐸𝐶_𝑆𝑑³𝜃^𝑃_𝑥𝑑𝑥³, (71a)(71b)where𝐶_𝑆=_Ω𝜙^𝑃_𝑆²𝐼𝑑Ω,_𝑃𝑡=_Ω𝑦²+𝑧²+𝑦⋅𝜕𝜙^𝑃_𝑆𝜕𝑧−𝑧⋅𝜕𝜙^𝑃_𝑆𝜕𝑦𝑑Ω (72a)(72b)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 ∇²𝜑^𝑆_𝑆=𝑀^𝑆_𝑡𝐺𝐶_𝑆𝜑^𝑃_𝑆𝑦,𝑧inΩ. (73) 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]: 𝑑𝜃^𝑆_𝑥=𝑀𝑑𝑥𝐺𝐼^𝑆𝑡𝑆_𝑡, (74) 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 𝐼^𝑆_𝑡=𝜅^𝑆_𝑥𝐼^𝑃_𝑡=1𝛼^𝑆_𝑥𝐼^𝑃_𝑡, (75) 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 𝑑𝜃^𝑆_𝑥𝑑𝑥=𝛼^𝑆_𝑥𝑀^𝑆_𝑡𝐺𝐼^𝑃_𝑡. (76) 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 𝑈_appr=12𝑀^𝑆_𝑡𝑑𝜃^𝑆_𝑥=1𝑑𝑥2𝑀^𝑆_𝑡²𝐺𝐼^𝑆_𝑡 (77) with the exact one given from 𝑈_exact=12𝐺_Ω𝜏^𝑆2_𝑥𝑦+𝜏^𝑆2_𝑥𝑧=1𝑑Ω2𝐺_Ω𝜕𝜑^𝑆_𝑆𝜕𝑦𝜕𝜑^𝑆_𝑆+𝜕𝑦𝜕𝜑^𝑆_𝑆𝜕𝑧𝜕𝜑^𝑆_𝑆𝜕𝑧𝑑Ω (78) the secondary torsion constant can be obtained as𝐼^𝑆_𝑡=𝑀^𝑆_𝑡²𝐺²×1∫_Ω𝜕𝜑^𝑆_𝑆/𝜕𝑦𝜕𝜑^𝑆_𝑆/𝜕𝑦+𝜕𝜑^𝑆_𝑆/𝜕𝑧𝜕𝜑^𝑆_𝑆./𝜕𝑧𝑑Ω (79) 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 𝐼^𝑆_𝑡=1∫_Ω𝜕𝜑^𝑆_𝑆/𝜕𝑦𝜕𝜑^𝑆_𝑆/𝜕𝑦+𝜕𝜑^𝑆_𝑆/𝜕𝑧𝜕𝜑^𝑆_𝑆,/𝜕𝑧𝑑Ω (80) 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:∇²𝜑^𝑆_𝑆=𝜑^𝑃_𝑆𝐶_𝑆𝜕𝜑inΩ,_𝑆𝑆𝜕𝑛=0onΓ. (81a)(81b)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 𝐼^𝑆_𝑡=𝐶_𝑆𝐼_𝜑, (82) where 𝐼_𝜑 is a domain integral given from the relation 𝐼_𝜑=−_Ω𝜑^𝑃_𝑆𝜑^𝑆_𝑆𝑑Ω (83) while form (75) and (82) the following relation for the warping shear deformation coefficient is obtained: 𝛼^𝑆_𝑥=𝐼^𝑃_𝑡𝐶_𝑆𝐼_𝜑. (84)

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 𝑑𝑀_𝑡𝑑𝑥=−𝑚_𝑡−𝑑𝑚_𝑤.𝑑𝑥 (85) Taking into account (29), (85) can be written as 𝑑𝑀^𝑆_𝑡+𝑑𝑥𝑑𝑀^𝑃_𝑡𝑑𝑥=−𝑚_𝑡−𝑑𝑚_𝑤𝑑𝑥 (86) while substituting (40) into (86) the following equilibrium equation is obtained: 𝑑²𝑀_𝑤𝑑𝑥²+𝑑𝑀^𝑃_𝑡𝑑𝑥=−𝑚_𝑡−𝑑𝑚_𝑤.𝑑𝑥 (87) 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 (𝜃^𝑃_𝑥).

Figure 12 

Equilibrium of a small segment 𝑑𝑥 of a twisted bar.

Introducing an auxiliary geometric constant 𝜅 defined as 1𝜅=1+𝛼^𝑆_𝑥 (88) the relation (76) can be written as 𝑑𝜃^𝑆_𝑥=𝑑𝑥1−𝜅𝜅𝑀^𝑆_𝑡𝐺𝐼^𝑃_𝑡 (89) while from relation (65b) the total torsional curvature can be written as 𝑑𝜃_𝑥=𝑑𝑥𝑑𝜃^𝑃_𝑥+𝑑𝑥1−𝜅𝜅𝑀^𝑆_𝑡𝐺𝐼^𝑃_𝑡 (90) and according to (71a) and (71b), the primary twisting moment is given with respect to the primary angle of twist (𝜃^𝑃_𝑥) as 𝑀^𝑃_𝑡=𝐺𝐼^𝑃_𝑡𝑑𝜃^𝑃_𝑥𝑑𝑥−𝐸𝐶_𝑆𝑑³𝜃^𝑃_𝑥𝑑𝑥³(1−𝜅), (91) where 𝐶_𝑆 is a modified warping constant given as 𝐶_𝑆=𝐶_𝑆𝜅. (92) 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 𝑀_𝑡=𝐺𝐼^𝑃_𝑡𝑑𝜃^𝑃_𝑥𝑑𝑥−𝐸𝐶_𝑆𝑑³𝜃^𝑃_𝑥𝑑𝑥³ (93) while using (92) the moments 𝑀^𝑆_𝑡 and 𝑀_𝑤 can be written in terms of the modified warping constant as 𝑀_𝑤=−𝐸𝐶_𝑆𝑑²𝜃^𝑃_𝑥𝑑𝑥²𝑀𝜅,_𝑆𝑡=−𝐸𝐶_𝑆𝑑³𝜃^𝑃_𝑥𝑑𝑥³𝜅. (94)(95) 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 𝐸𝐶_𝑆𝑑⁴𝜃^𝑃_𝑥𝑑𝑥⁴−𝐺𝐼^𝑃_𝑡𝑑²𝜃^𝑃_𝑥𝑑𝑥²=+𝑚_𝑡+𝑑𝑚_𝑤𝑑𝑥insidethebar (96) while the boundary conditions are given from𝛼₁𝜃^𝑃_𝑥+𝛼₂𝑀_𝑡=𝛼₃,𝛽₁𝑑𝜃^𝑃_𝑥𝑑𝑥+𝛽₂𝑀_𝑤=𝛽₃atthebarends𝑥=0,𝑙, (97a)(97b)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 𝑎₁=𝛽₁=1, 𝑎₂=𝑎₃=𝛽₂=𝛽₃=0).

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 𝑑²𝜃^𝑆_𝑥𝑑𝑥²𝐸=−𝐶_𝑆𝐺𝐼^𝑃_𝑡𝑑⁴𝜃^𝑃_𝑥𝑑𝑥⁴(1−𝜅)insidethebar (98) while the corresponding boundary conditions are given in terms of the primary angle of twist as 𝛾₁𝜃^𝑆_𝑥+𝛾₂𝐺𝐼^𝑃_𝑡𝜅1−𝜅𝑑𝜃^𝑆_𝑥𝑑𝑥=𝛾₃−𝐸𝐶_𝑆𝜅𝑑³𝜃^𝑃_𝑥𝑑𝑥³atthebarends𝑥=0,𝑙, (99) where if the primary angle of twist is fixed, that is 𝑎₁=1 and 𝑎₂=𝑎₃=0, then 𝛾₁=1 and 𝛾₂=𝛾₃=0, otherwise 𝛾₁=0 and 𝛾₂=𝛾₃=1.

From the examined example problems presented in [85] it is concluded that (i)the inaccuracy of the thin tube theory in calculating the secondary torsion constant even for thin-walled sections is noteworthy;(ii)as it is also verified by other computational methods [80, 83], for bars with open shaped cross section the secondary torsional moment deformation effect has not significant influence. However, the inclusion of this additional effect leads to more accurate results;(iii)as it is also verified by other computational methods [80, 83], for bars with closed shaped cross sections the secondary torsional moment deformation effect has an important influence and should be considered in the analysis (otherwise the results may be completely wrong especially in the calculation of stresses).

7. Nonlinear Elastic Nonuniform Torsion of Bars

As it has already been mentioned in the previous sections in engineering practice we often come across the analysis of members of structures subjected to twisting moments. Besides, since thin-walled open sections have low torsional stiffness, the torsional deformations can be of such magnitudes that it is not adequate to treat the angles of cross section rotation as small. When finite twist rotation angles are considered, the elastic nonuniform torsion problem becomes nonlinear. Moreover, this problem becomes much more complicated in the case the cross section’s centroid does not coincide with its shear center (asymmetric beams), leading to the formulation of a flexural-torsional coupled problem. The extensive use of the aforementioned structural elements necessitates a reliable and accurate analysis of bars of arbitrary cross section subjected to torsional loading taking into account the geometrical nonlinearity.

Though several researchers have dealt either with the linear nonuniform torsional behaviour of beams [17, 21, 22, 40] or with the nonlinear uniform torsional behaviour of doubly symmetric beams [89, 90], to the author’s knowledge very little work has been done on the corresponding nonlinear nonuniform torsional problem of arbitrary cross section beams. Ghobarah and Tso [91], Attard [92], and Attard and Somervaille [93] have presented a set of displacement relationships for a straight prismatic thin-walled open beam applicable to situations where displacements are finite, the cross section does not distort, strains are small and flexural displacements are small to moderate while cross sectional twist can be large. The presented numerical examples in these studies are concerned only with uniform torsion of either mono- or doubly symmetric cross sections. Finally, Trahair in [94] employing the finite element method and presenting examples of only doubly symmetric cross sections and Mohri et al. in [95] employing similar equations to those established by Attard in [92] and presenting examples of either doubly symmetric cross sections subjected in nonuniform torsion or buckling or postbuckling behavior of arbitrary cross section beams also analyze the nonlinear nonuniform torsional problem. Nevertheless, all of the aforementioned studies, which are the only one considering finite angles of twist in asymmetrical bars (and taking into account all of the arising nonlinear terms) are not general since they are restricted to thin-walled beams. Finally, Sapountzakis and Tsipiras in [96, 97] presented a BE solution for the elastic nonuniform torsion analysis of simply, multiply connected or composite cylindrical bars of arbitrary cross section taking into account the effect of geometric nonlinearity. The torque-rotation relationship is computed based on the finite displacement (finite rotation) theory; that is, the transverse displacement components are expressed so as to be valid for large rotations and the longitudinal normal strain includes the second-order geometrically nonlinear term often described as the “Wagner strain.” These last formulations do not stand on the assumption of a thin-walled structure and therefore the cross section’s torsional rigidity is evaluated exactly without using the so-called Saint-Venant’s torsional constant. The torsional rigidity of the cross section is evaluated directly employing the primary warping function of the cross section [21] depending on its shape. Three boundary value problems with respect to the variable along the beam axis angle of twist, to the primary and to the secondary warping functions are formulated. The first one of these problems is numerically solved employing the Analog Equation Method [98], a BEM based method, leading to a system of nonlinear equations from which the angle of twist is computed by an iterative process. The other two problems are solved using a pure BEM [23] based method. The aforementioned formulation procedure is based on the assumption of no local or lateral torsional buckling or distortion.

In order to formulate the aforementioned problem, let us consider a prismatic bar of 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)).

(a)

(b)

Figure 13 

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

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𝑢=𝑢_𝑚(𝑥)+𝜃_𝑌(𝑥)⋅𝑧−𝑧_𝐶−𝜃_𝑍(𝑥)⋅𝑦−𝑦_𝐶+𝜃^′_𝑥(𝑥)⋅𝜙^𝑃_𝑆𝑦,𝑧+𝜙^𝑆_𝑆𝑥,𝑦,𝑧,𝑣=𝑣_𝑆(𝑥)−𝑧⋅sin𝜃_𝑥(𝑥)−𝑦⋅1−cos𝜃_𝑥(𝑥),𝑤=𝑤_𝑆(𝑥)+𝑦⋅sin𝜃_𝑥(𝑥)−𝑧⋅1−cos𝜃_𝑥(𝑥), (100a)(100b)(100c)where the transverse displacement components 𝑣, 𝑤 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; 𝑣_𝑆(𝑥), 𝑤_𝑆(𝑥) 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𝜀_𝑥𝑥=𝜕𝑢+1𝜕𝑥2⋅𝜕𝑢𝜕𝑥²+𝜕𝑣𝜕𝑥²+𝜕𝑤𝜕𝑥²,𝛾_𝑥𝑦=𝜕𝑢+𝜕𝑦𝜕𝑣+𝜕𝑥𝜕𝑢⋅𝜕𝑦𝜕𝑢+𝜕𝑥𝜕𝑣⋅𝜕𝑦𝜕𝑣+𝜕𝑥𝜕𝑤⋅𝜕𝑦𝜕𝑤,𝛾𝜕𝑥_𝑥𝑧=𝜕𝑤+𝜕𝑥𝜕𝑢+𝜕𝑧𝜕𝑢⋅𝜕𝑧𝜕𝑢+𝜕𝑥𝜕𝑣⋅𝜕𝑧𝜕𝑣+𝜕𝑥𝜕𝑤⋅𝜕𝑧𝜕𝑤𝜕𝑥 (101a)(101b)(101c)assuming moderate large deflections ((𝜕𝑢/𝜕𝑥)²≪𝜕𝑢/𝜕𝑥, (𝜕𝑢/𝜕𝑥)(𝜕𝑢/𝜕𝑦)≪(𝜕𝑣/𝜕𝑥)+(𝜕𝑢/𝜕𝑦), (𝜕𝑢/𝜕𝑥)(𝜕𝑢/𝜕𝑧)≪(𝜕𝑤/𝜕𝑥)+(𝜕𝑢/𝜕𝑧)) and having in mind that ignoring shear deformation effect the angles of rotation due to bending 𝜃_𝑌, 𝜃_𝑍 approximate the slope between the longitudinal axis in the deformed state and in the initial vertical and horizontal planes, respectively, the following relations are obtained:𝜃_𝑌(𝑥)=𝑣^′_𝑆⋅sin𝜃_𝑥−𝑤^′_𝑆⋅cos𝜃_𝑥,𝜃_𝑍(𝑥)=𝑣^′_𝑆⋅cos𝜃_𝑥+𝑤^′_𝑆⋅sin𝜃_𝑥 (102a)(102b)and the nonvanishing strain resultants are given as𝜀_𝑥𝑥=𝑢^′_𝑚+𝜅_𝑌⋅𝑧−𝑧_𝐶−𝜅_𝑍⋅𝑦−𝑦_𝐶+𝜃_𝑥^′′⋅𝜙^𝑃_𝑆𝑦,𝑧−𝜃^′_𝑥⋅𝑦_𝐶⋅𝜃_𝑌+𝑧_𝐶⋅𝜃_𝑍+12𝑣^′_𝑆²+𝑤^′_𝑆²+𝑦²+𝑧²⋅𝜃^′_𝑥²,𝛾_𝑥𝑦=𝜃^′_𝑥⋅𝜕𝜙^𝑃_𝑆+𝜕𝑦−𝑧𝜕𝜙^𝑆_𝑆,𝛾𝜕𝑦_𝑥𝑧=𝜃^′_𝑥⋅𝜕𝜙^𝑃_𝑆+𝜕𝑧+𝑦𝜕𝜙^𝑆_𝑆,𝜕𝑧 (103a)(103b)(103c)where the curvature components 𝜅_𝑌, 𝜅_𝑍 are given from the following relations:𝜅_𝑌(𝑥)=𝑣_𝑆^′′(𝑥)⋅sin𝜃_𝑥−𝑤_𝑆^′′(𝑥)⋅cos𝜃_𝑥,𝜅_𝑍(𝑥)=𝑣_𝑆^′′(𝑥)⋅cos𝜃_𝑥+𝑤_𝑆^′′(𝑥)⋅sin𝜃_𝑥 (104a)(104b)while the second-order geometrically nonlinear term in the right hand side of (103a) (𝑦²+𝑧²)⋅(𝜃^′_𝑥)²/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 𝑆_𝑥𝑥𝑆_𝑥𝑦𝑆_𝑥𝑧=𝐸00𝜀0𝐺000𝐺^𝑥𝑥𝛾_𝑥𝑦𝛾_𝑥𝑧 (105) or employing (103a), (103b), and (103c) as𝑆_𝑥𝑥=𝐸𝑢^′_𝑚+𝜅_𝑌⋅𝑧−𝑧_𝐶−𝜅_𝑍⋅𝑦−𝑦_𝐶+𝜃_𝑥^′′⋅𝜙^𝑃_𝑆𝑦,𝑧−𝜃^′_𝑥⋅𝑦_𝐶𝜃_𝑌+𝑧_𝐶𝜃_𝑍+𝐸2𝑣^′_𝑆²+𝑤^′_𝑆²+𝑦²+𝑧²⋅𝜃^′_𝑥²,𝑆_𝑥𝑦=𝐺⋅𝜃^′_𝑥⋅𝜕𝜙^𝑃_𝑆𝜕𝑦−𝑧+𝐺⋅𝜕𝜙^𝑆_𝑆,𝑆𝜕𝑦_𝑥𝑧=𝐺⋅𝜃^′_𝑥⋅𝜕𝜙^𝑃_𝑆𝜕𝑧+𝑦+𝐺⋅𝜕𝜙^𝑆_𝑆.𝜕𝑧 (106a)(106b)(106c)

7.1. Equations of Local Equilibrium

Applying the stress components (106a), (106b), and (106c) in the first equation of equilibrium neglecting the body forces [100] 𝜕𝑆_𝑥𝑥+𝜕𝑥𝜕𝑆_𝑥𝑦+𝜕𝑦𝜕𝑆_𝑥𝑧𝜕𝑧=0 (107) ignoring the independent of warping terms such as 𝐸⋅𝑢_𝑚^′′ or warping terms due to shear such as 𝐸(𝑤_𝑆^′′′(𝑧−𝑧_𝐶)+𝑣_𝑆^′′′(𝑦−𝑦_𝐶)) 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∇²𝜙^𝑆_𝑆𝐸=−𝐺⋅𝜃_𝑥^′′′⋅𝜙^𝑃_𝑆−𝐸𝑦,𝑧𝐺⋅𝑦²+𝑧²⋅𝜃^′_𝑥⋅𝜃_𝑥^′′inΩ,𝜕𝜙^𝑆_𝑆𝜕𝑛=0onΓ (108a)(108b)are obtained, where ∇²=𝜕²/𝜕𝑦²+𝜕²/𝜕𝑧² 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 _Ω𝜙^𝑃_𝑆𝑑Ω=0,_Ω𝜙^𝑆_𝑆𝑑Ω=0. (109)(110) 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_Ω𝜙^𝑃_𝑆𝑦𝑑Ω=0,_Ω𝜙^𝑃_𝑆𝑧𝑑Ω=0 (111a)(111b)that define point 𝑆 as the point for which the axial and bending stress resultants arising from an angle of twist at this point vanish in a geometrically linear analysis. It is worth here noting that the shear center axis of an asymmetric cross section bar subjected in pure torsion is not a rectilinear one but is assigned from the displacement field 𝑢_𝑆(𝑥), 𝑣_𝑆(𝑥), 𝑤_𝑆(𝑥).

Based on the conditions (109) and (110) the meaning of the “average” axial displacement of the cross section of the bar can now be explained as the relation _Ω𝑢𝑑Ω=𝑢_𝑚(𝑥)⋅𝐴+𝜃_𝑌(𝑥)⋅_Ω𝑧𝑑Ω−𝑧_𝐶⋅𝐴−𝜃_𝑍(𝑥)⋅_Ω𝑦𝑑Ω−𝑦_𝐶⋅𝐴+𝜃^′_𝑥⋅_Ω𝜙^𝑃_𝑆+𝑑Ω_Ω𝜙^𝑆_𝑆𝑑Ω (112) leads to 𝑢_𝑚(𝑥)=∫_Ω𝑢𝑑Ω𝐴 (113) that is 𝑢_𝑚 generally does not coincide either with the axial displacement of the cross section’s centroid 𝑢_𝐶 or with that of the shear center 𝑢_𝑆.

7.2. Equations of Global Equilibrium

Neglecting body forces the principle of virtual work yields _𝑉𝑆_𝑥𝑥⋅𝛿𝜀_𝑥𝑥+𝑆_𝑥𝑦⋅𝛿𝛾_𝑥𝑦+𝑆_𝑥𝑧⋅𝛿𝛾_𝑥𝑧=𝑑𝑉_𝐴𝑡_𝑥⋅𝛿𝑢+𝑡_𝑦⋅𝛿𝑣+𝑡_𝑧⋅𝛿𝑤𝑑𝐴, (114) 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𝑁=_Ω𝑆_𝑥𝑥𝑀𝑑Ω,_𝑌=_Ω𝑆_𝑥𝑥𝑧−𝑧_𝐶𝑀𝑑Ω,_𝑍=−_Ω𝑆_𝑥𝑥𝑦−𝑦_𝐶𝑀𝑑Ω,_𝑃𝑡=_Ω𝑆^𝑃_𝑥𝑦⋅𝜕𝜙^𝑃_𝑆𝜕𝑦−𝑧+𝑆^𝑃_𝑥𝑧⋅𝜕𝜙^𝑃_𝑆𝑀𝜕𝑧+𝑦𝑑Ω,_𝑤=−_Ω𝑆_𝑥𝑥𝜙^𝑃_𝑆𝑑Ω, (115a)(115b)(115c)(115d)(115e)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𝑁=𝐸𝐴⋅𝑢^′_𝑚+12⋅𝑣^′_𝑆²+𝑤^′_𝑆²+𝐼_𝑃𝐴⋅𝜃^′_𝑥²−𝜃^′_𝑥⋅𝑦_𝐶𝜃_𝑌+𝑧_𝐶𝜃_𝑍,𝑀_𝑌=𝐸𝐼_𝑌𝑌⋅𝜅_𝑌+𝛽₂𝜃^′_𝑥²,𝑀_𝑍=𝐸𝐼_𝑍𝑍⋅𝜅_𝑍−𝛽₁𝜃^′_𝑥²,𝑀^𝑃_𝑡=𝐺𝐼_𝑡⋅𝜃^′_𝑥,𝑀_𝑤=−𝐸𝐶_𝑆⋅𝜃_𝑥^′′+𝑈_𝑤2𝐶_𝑆⋅𝜃^′_𝑥², (116a)(116b)(116c)(116d)(116e)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𝐼_𝑃=_Ω𝑟²𝐼𝑑Ω,_𝑌𝑌=_Ω𝑧−𝑧_𝐶²𝐼𝑑Ω,_𝑍𝑍=_Ω𝑦−𝑦_𝐶²𝑑Ω (117a)(117b)(117c)with 𝑟=𝑦²+𝑧², while the geometric constants 𝛽₁, 𝛽₂, and 𝑈_𝑤 are given as𝛽₁=12𝐼_𝑍𝑍⋅_Ω𝑟²⋅𝑦−𝑦_𝐶𝛽𝑑Ω,₂=12𝐼_𝑌𝑌⋅_Ω𝑟²⋅𝑧−𝑧_𝐶𝑈𝑑Ω,_𝑤=_Ω𝜙^𝑃_𝑆⋅𝑟²𝑑Ω. (118a)(118b)(118c)

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 −𝑁𝑦_𝐶𝜃_𝑍𝜃^′_𝑥+𝑁𝑧_𝐶𝜃_𝑌𝜃^′_𝑥−𝑀_𝑍𝜅_𝑌+𝑀_𝑌𝜅_𝑍−𝑑𝑀𝑑𝑥^𝑃𝑡+12𝐸𝐼_𝑛𝜃^′_𝑥³+Ψ𝜃^′_𝑥−𝑁𝑦_𝐶𝜃_𝑌−𝑁𝑧_𝐶𝜃_𝑍−𝑑²𝑀_𝑤𝑑𝑥²=𝑚_𝑡(𝑥)+𝑑𝑚𝑑𝑥^𝑤(𝑥) (119) and the corresponding boundary conditions are written as𝑀_𝑤^′+𝑀^𝑃_𝑡+12𝐸𝐼_𝑛𝜃^′_𝑥³+Ψ𝜃^′_𝑥−𝑁𝑦_𝐶𝜃_𝑌−𝑁𝑧_𝐶𝜃_𝑍−𝑀_𝑡𝛿𝜃_𝑥=0,−𝑀_𝑤+𝑀_𝑤𝛿𝜃^′_𝑥=0, (120a)(120b)where the “higher order torsion constant” 𝐼_𝑛 and the variable Ψ are given as𝐼_𝑛=𝐼_𝑃𝑃−𝐼²_𝑃𝐴−4𝛽²₁𝐼_𝑍𝑍−4𝛽²₂𝐼_𝑌𝑌−𝑈²_𝑤𝐶_𝑆,𝐼Ψ=𝑁^𝑃𝐴−2𝑀_𝑍𝛽₁+2𝑀_𝑌𝛽₂−𝑀_𝑤𝑈_𝑤𝐶_𝑆 (121a)(121b)while 𝑀_𝑡, 𝑀_𝑤 in the boundary conditions (120a) and (120b) are the externally applied conservative twisting and warping moments at the bar ends and 𝐼_𝑃𝑃=∫_Ω𝑟⁴𝑑Ω (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𝑀𝑁=0,_𝑌𝑀=0,_𝑍=0 (122a)(122b)(122c)the governing equation of the nonuniform nonlinear torsional problem is written as 𝐸𝐶_𝑆𝜃_𝑥^′′′′−32𝐸𝐼_𝑛2𝜃^′_𝑥²𝜃_𝑥^′′−𝐺𝐼_𝑡𝜃_𝑥^′′=𝑚_𝑡(𝑥)+𝑑𝑚𝑑𝑥^𝑤(𝑥) (123) subjected to the following boundary conditions at the bar ends𝛼₁𝜃_𝑥+𝛼₂𝑀_𝑡=𝛼₃,𝛽₁𝜃^′_𝑥+𝛽₂𝑀_𝑤=𝛽₃, (124a)(124b)where 𝐼_𝑛2 is defined as 𝐼_𝑛2=𝐼_𝑛+𝑈²_𝑤𝐶_𝑆 (125) and 𝑀_𝑡, 𝑀_𝑤 are the twisting and warping moments at the boundary of the bar, respectively, given as𝑀_𝑡=−𝐸𝐶_𝑆𝜃_𝑥^′′′+𝐺𝐼_𝑡𝜃^′_𝑥+12𝐸𝐼_𝑛2𝜃^′_𝑥³,𝑀_𝑤=−𝐸𝐶_𝑆𝜃_𝑥^′′+𝑈_𝑤2𝐶_𝑆𝜃^′_𝑥² (126a)(126b)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 𝛼₁=𝛽₁=1, 𝛼₂=𝛼₃=𝛽₂=𝛽₃=0).

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 𝑣_𝑆(𝑥), 𝑤_𝑆(𝑥) can be established integrating twice the following relations: 𝑣_𝑆^′′=−𝜃^′_𝑥²⋅𝛽₂⋅sin𝜃_𝑥−𝛽₁⋅cos𝜃_𝑥,𝑤_𝑆^′′=𝜃^′_𝑥²⋅𝛽₂⋅cos𝜃_𝑥+𝛽₁⋅sin𝜃_𝑥 (127a)(127b)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 𝑣^′_𝑆, 𝑤^′_𝑆 the “average” axial displacement of the cross section of the bar 𝑢_𝑚 can be established after substitution of (116a) to (122a) as 𝑢^′_𝑚1=−2⋅𝑣^′_𝑆²+𝑤^′_𝑆²+𝐼_𝑃𝐴⋅𝜃^′_𝑥²+𝜃^′_𝑥⋅𝑦_𝐶𝜃_𝑌+𝑧_𝐶𝜃_𝑍 (128) and subsequent integration.

From the examined example problems presented in [96, 97] it is concluded that(i)geometrical nonlinearity leads to stiffening of the bar and eventually its better response against torsional loading;(ii)as observed, restrained warping boundary conditions lead to slightly greater torsional resistance than those of free warping for both the geometrically linear and nonlinear cases;(iii)axial shortening of bars of doubly symmetric cross section subjected to torsional loading is observed in the nonlinear analysis, while bars of asymmetrical cross section exhibit lateral displacement of their shear center axis as well;(iv)the nonlinear Wagner torque can reach significant values locally along the bar axis while the nonlinear warping moment can be neglected in most cases;(v)the secondary shear stress distribution and magnitude are influenced by the geometrical nonlinearity.

8. Nonlinear Elastic Torsional Vibrations of Bars

As it has already been mentioned in the previous section, since weight saving is of paramount importance, frequently used thin-walled open sections have low torsional stiffness and their torsional deformations can be of such magnitudes that it is not adequate to treat the angles of cross section rotation as small. In these cases, the study of nonlinear effects on these members becomes essential, where this nonlinearity results from retaining the nonlinear terms in the strain-displacement relations (finite displacement—small strain theory). When finite twist rotation angles are considered, the nonuniform torsional dynamic analysis of bars becomes much more complicated, leading to the formulation of coupled and nonlinear torsional and axial equilibrium equations. In this case, accounting for the axial loading and boundary conditions becomes essential to perform a rigorous dynamic analysis of the bar.

During the past few years, the nonlinear nonuniform torsional dynamic analysis of bars undergoing large rotations has received a good amount of attention in the literature. More specifically, Rozmarynowski and Szymczak in [101] studied the nonlinear free torsional vibrations of axially immovable thin-walled beams with doubly symmetric open cross section, employing the finite element method. In this research effort only free vibrations are examined, the solution is provided only at points of reversal of motion (not in the time domain), no general axial, torsional, or warping boundary conditions (elastic support case) are studied, while some nonlinear terms related to the finite twisting rotations as well as the axial inertia term are ignored. Da Silva in [102, 103] presented the nonlinear flexural-torsional-extensional vibrations of Euler-Bernoulli doubly symmetric thin-walled closed cross section beams, primarily focusing on flexural vibrations and neglecting the effect of torsional warping. Pai and Nayfeh in [104–106] studied also the nonlinear flexural-torsional-extensional vibrations of metallic and composite slewing or rotating closed cross section beams, primarily focusing on flexural vibrations and neglecting again the effect of torsional warping. Di Egidio et al. in [107, 108] presented also a FEM solution to the nonlinear flexural-torsional vibrations of shear undeformable thin-walled open beams taking into account in-plane and out-of-plane warpings and neglecting warping inertia. In these papers, the torsional-extensional coupling is taken into account but the inextensionality assumption leads to the fact that the axial boundary conditions are not general. Simo and Vu-Quoc in [109] presented a FEM solution to a fully nonlinear (small or large strains, hyperelastic material) three dimensional rod model including the effects of transverse shear and torsion-warping deformation based on a geometrically exact description of the kinematics of deformation. Moreover, Pai and Nayfeh in [110] studied a geometrically exact nonlinear curved beam model for solid composite rotor blades using the concept of local engineering stress and strain measures and taking into account the in-plane and out-of-plane warpings. In all of the aforementioned research efforts the evaluation of the warping shear stress distribution is not discussed. Finally, Sapountzakis and Tsipiras in [111, 112] developed a boundary element solution for the nonuniform torsional vibration problem of bars of arbitrary doubly symmetric constant cross section taking into account the effect of geometrical nonlinearity. In these last research efforts the bar is subjected to arbitrarily distributed or concentrated conservative dynamic twisting and warping moments along its length, while its edges are supported by the most general torsional and axial boundary conditions. A distributed mass model system is employed, taking into account the warping, rotatory, and axial inertia, which after employing a variational approach leads to the formulation of a coupled nonlinear initial boundary value problem with respect to the variable along the bar angle of twist and to an “average” axial displacement of the cross section of the bar and to two boundary value problems with respect to the primary and secondary warping functions. The numerical solution of the aforementioned initial boundary value problem is performed using the Analog Equation Method [98], a BEM based method, leading to a system of nonlinear Differential-Algebraic Equations (DAEs), which is solved using an efficient time discretization scheme. The other two problems are solved using a pure BE method, requiring exclusively boundary discretization of the bar’s cross-section. Both free and forced vibrations are considered, while the arising linear system of equations related to the secondary warping function is singular and a special technique [113] is used to perform its regularization.

In order to formulate the nonlinear torsional vibration problem of doubly symmetric bars of arbitrarily shaped simply or multiply connected cross section, let us consider the bar of 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)).

(a)

(b)

Figure 14 

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

For the analysis of the aforementioned problem, the same displacement field (100a), (100b), and (100c) with that of the nonlinear static analysis is adopted, which for the dynamic problem and for a doubly symmetric cross section is written as𝑢𝑥,𝑦,𝑧,𝑡=𝑢_𝑚(𝑥,𝑡)+𝜃^′_𝑥(𝑥,𝑡)⋅𝜙^𝑃_𝑆𝑦,𝑧+𝜙^𝑆_𝑆,𝑣𝑥,𝑦,𝑧,𝑡𝑥,𝑦,𝑧,𝑡=−𝑧⋅sin𝜃_𝑥(𝑥,𝑡)−𝑦⋅1−cos𝜃_𝑥(),𝑤𝑥,𝑡𝑥,𝑦,𝑧,𝑡=𝑦⋅sin𝜃_𝑥(𝑥,𝑡)−𝑧⋅1−cos𝜃_𝑥(𝑥,𝑡). (129a)(129b)(129c)

8.1. Equations of Local Equilibrium

Following the same procedure presented in Section 7, the local equilibrium equations are established employing the principle of virtual work 𝛿𝑊_int+𝛿𝑊_mass=𝛿𝑊_ext, (130) where𝛿𝑊_int=_𝑉𝑆_𝑥𝑥𝛿𝜀_𝑥𝑥+𝑆_𝑥𝑦𝛿𝛾_𝑥𝑦+𝑆_𝑥𝑧𝛿𝛾_𝑥𝑧𝑑𝑉,𝛿𝑊_mass=_𝑉𝜌⋅(𝑢̈𝛿𝑢+𝑣𝛿̈𝑣+𝑤̈𝛿𝑤)𝑑𝑉,𝛿𝑊_ext=_𝐹𝑡_𝑥⋅𝛿𝑢+𝑡_𝑦⋅𝛿𝑣+𝑡_𝑧⋅𝛿𝑤𝑑𝐹 (131a)(131b)(131c)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 𝜕𝑆_𝑥𝑥+𝜕𝑥𝜕𝑆_𝑥𝑦+𝜕𝑦𝜕𝑆_𝑥𝑧̈()𝜕𝑧−𝜌𝑢=0inΩ,∀𝑥∈0,𝑙 (132) along with its corresponding boundary condition 𝑆_𝑥𝑦𝑛_𝑦+𝑆_𝑥𝑧𝑛_𝑧=0onΓ_𝑗,∀𝑥∈(0,𝑙) (133) 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 ∇²𝜙^𝑆_𝑆=−𝐸𝐺𝑢_𝑚^′′+𝜌𝐺𝑢̈_𝑚−𝐸𝐺𝜃_𝑥^′′′−𝜌𝐺𝜃̈^′_𝑥⋅𝜙^𝑃_𝑆−𝐸𝑦,𝑧𝐺𝜃^′_𝑥𝜃_𝑥^′′⋅𝑦²+𝑧²inΩ (134) along with its corresponding boundary condition 𝜕𝜙^𝑆_𝑆𝜕𝑛=0onΓ_𝑗𝑗=1,2,…,𝐾 (135) 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.

8.2. Equations of Global Equilibrium

Taking into account the double symmetry of the cross section and substituting the stress components ((106a), (106b), and (106c)) the strain ones ((103a), (103b), and (103c)) and the displacement components ((129a), (129b), and (129c)) to the principle of virtual work (130), the governing partial differential equations of the bar are obtained after some algebra as𝜌𝐴⋅𝑢̈_𝑚−𝐸𝐴⋅𝑢_𝑚^′′−𝐸𝐼_𝑃⋅𝜃^′_𝑥𝜃_𝑥^′′=𝑛(𝑥,𝑡),𝜌𝐼_𝑃⋅𝜃̈_𝑥−𝜌𝐶_𝑆⋅𝜃̈_𝑥^′′−𝐺𝐼_𝑡𝜃_𝑥^′′+𝐸𝐶_𝑆𝜃_𝑥^′′′′−32𝐸𝐼_𝑃𝑃⋅𝜃^′_𝑥²𝜃_𝑥^′′−𝐸𝐼_𝑃⋅𝑢^′_𝑚𝜃_𝑥^′′−𝐸𝐼_𝑃⋅𝑢_𝑚^′′𝜃^′_𝑥=𝑚_𝑡()+𝜕𝑥,𝑡𝑚𝜕𝑥^𝑤(𝑥,𝑡) (136a)(136b)subjected to the initial conditions (𝑥∈(0,𝑙))𝑢_𝑚(𝑥,0)=𝑢_𝑚0(𝑥),𝑢̇_𝑚(𝑥,0)=𝑢̇_𝑚0(𝑥),𝜃_𝑥(𝑥,0)=𝜃_𝑥0(𝑥),𝜃̇_𝑥(𝑥,0)=𝜃̇_𝑥0(𝑥) (137a)(137b)(137c)(137d)together with the boundary conditions at the bar ends 𝑥=0,𝑙𝑎₁𝑁+𝛼₂𝑢_𝑚=𝛼₃,𝛽₁𝑀_𝑡+𝛽₂𝜃_𝑥=𝛽₃,𝛽₁𝑀_𝑤+𝛽₂𝜃^′_𝑥=𝛽₃, (138a)(138b)(138c)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.

8.3. Reduced Problems for Special Cases of Axial Boundary Conditions

For the case of axially immovable ends, it is easily proved that [112] 𝑢_𝑚^′′𝐼=−^𝑃𝐴⋅𝜃^′_𝑥𝜃_𝑥^′′,∀𝑥∈(0,𝑙) (139) which after subsequent integration yields 𝑢^′_𝑚1=−2⋅𝐼_𝑃𝐴⋅𝜃^′_𝑥²+𝑁[],𝐸𝐴,∀𝑥∈0,𝑙 (140) where 𝑁 is a time-dependent tensile axial load induced by the geometrical nonlinearity given as 1𝑁=2⋅𝐸𝐼_𝑃𝑙⋅^𝑙₀𝜃^′_𝑥²𝑑𝑥. (141) 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 𝜌𝐼_𝑃⋅𝜃̈_𝑥−𝜌𝐶_𝑆⋅𝜃̈_𝑥^′′−𝐺𝐼_𝑡+𝐼_𝑃𝐴⋅𝑁𝜃_𝑥^′′+𝐸𝐶_𝑆𝜃_𝑥^′′′′−32𝐸𝐼_𝑛⋅𝜃^′_𝑥²𝜃_𝑥^′′=𝑚_𝑡()+𝜕𝑥,𝑡𝑚𝜕𝑥^𝑤(𝑥,𝑡), (142) 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∇²𝜙^𝑆_𝑆red.=𝐸𝐺𝜃^′_𝑥𝜃_𝑥^′′⋅𝐼_𝑃𝐴−𝑦²+𝑧²−𝐸𝐺𝜃_𝑥^′′′−𝜌𝐺𝜃̈^′_𝑥⋅𝜙^𝑃_𝑆𝜕𝜙𝑦,𝑧inΩ_𝑆𝑆red.𝜕𝑛=0onΓ_𝑗𝑗=1,2,…,𝐾. (143a)(143b)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 𝐼_𝑛=𝐼_𝑃𝑃−𝐼²_𝑃𝐴. (144)

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.

9. Inelastic Nonuniform Torsion of Bars

Design of bars and bar assemblages subjected to twisting moments based on elastic analysis is most likely to be extremely conservative not only due to significant difference between initial yield and full plastification in a cross section, but also due to the unaccounted for yet significant reserves of strength that are not mobilized in redundant members until after inelastic redistribution takes place. Thus, material nonlinearity is important for investigating the ultimate strength of a bar that resists torsional loading, while distributed plasticity models are acknowledged in the literature [114–116] to capture more rigorously material nonlinearities than cross sectional stress resultant approaches [117] or lumped plasticity idealizations [118, 119]. If inelastic effects are considered, especially through distributed plasticity formulations, then the nonuniform [120, 121] torsional problem requires a much more rigorous analysis.

Apart from research efforts in which bars are idealized with computationally demanding three dimensional [122] or shell [123] elements, several researchers proposed specialized beam elements to analyze bars under inelastic nonuniform torsion. Bathe and Wiener [123] employed Hermitian and isoparametric beam elements to model thin-walled I-beams using one element for each flange and one element for the web of the cross section. Wunderlich et al. [124] employed a power series numerical technique using an Updated Lagrangian formulation to study thin-walled beams under general loading conditions. Izzuddin and Lloyd Smith [125, 126] followed an Eulerian approach employing cubic shape functions to study thin walled bars under general loading conditions. Chen and Trahair [99] and Pi and Trahair [121], using a mitre model [127] to describe the shear strain distribution over the cross section, developed finite element models for the inelastic analysis of thin-walled I-beams under torsion. Nie and Zhong [128] and Wang et al. [129] studied thin-walled beams under general loading conditions by employing an Updated Lagrangian description and the FE method.

The aforementioned contributions are applicable to cross sections of special geometry (e.g., thin-walled ones). A more general approach is presented by Baba and Kajita [120] who used a two-node, four-degree-of-freedom element for the longitudinal modeling and a four-node, 12-degree-of-freedom rectangular element for the cross sectional modeling of bars of bisymmetrical section under a Total Lagrangian formulation, taking into account plasticity effects in the plane of the cross section as well. Moreover, in the more recent contributions of Battini and Pacoste [130], Nukala and White [114] and Gruttmann et al. [131], beams of arbitrarily shaped cross section under general loading conditions are analyzed. Battini and Pacoste [130] and Nukala and White [114] use an independent warping parameter associated with the intensity of warping deformations to model nonuniform torsion. This parameter is eventually equated with the angle of twist per unit length in their analyses, thus only St. Venant shear stresses are taken into account. In all of the aforementioned research efforts warping shear stresses are not taken into account, with the exceptions of Wunderlich et al. [124], Gruttmann et al. [131], Nie and Zhong [128] and Wang et al. [129] who exploit a stress distribution arising from the introduction of an independent warping parameter in the displacement field of the bar. However, since this theory is analogous to the Timoshenko beam theory of shear-bending loading conditions, it does not satisfy local equilibrium equations under inelastic or even elastic conditions (for relevant discussions see e.g., Simo et al. [132] and Minghini et al. [133]). Finally, Sapountzakis and Tsipiras in [134] presented a boundary element solution for the inelastic nonuniform torsional problem of simply or multiply connected cylindrical bars of arbitrarily shaped doubly symmetric cross section taking into account the effect of warping shear stresses. The bar is subjected to arbitrarily distributed or concentrated torsional loading along its length, while its edges are subjected to the most general torsional boundary conditions. In this last research effort a displacement based formulation is developed and inelastic redistribution is modelled through a distributed plasticity model exploiting three-dimensional material constitutive laws and numerical integration over the cross sections. An incremental-iterative solution strategy is adopted to restore global equilibrium along with an efficient iterative process to integrate the inelastic rate equations [135]. Three boundary value problems with respect to the variable along the bar axis angle of twist to the primary and to the secondary warping functions are formulated and solved employing the boundary element method [23].

In order to formulate the inelastic torsional problem, let us consider a prismatic bar of 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)).

(a)

(b)

Figure 15 

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

Under the aforementioned loading, the displacement field of the bar taking into account warping shear stresses is assumed to be given as𝑢𝑥,𝑦,𝑧=𝜃^′_𝑥(𝑥)𝜙^𝑃_𝑆𝑦,𝑧+𝜃_𝑥^′′′(𝑥)𝜙^𝑆_𝑆,𝑣()𝑦,𝑧𝑥,𝑧=−𝑧𝜃_𝑥(𝑥),𝑤𝑥,𝑦=𝑦𝜃_𝑥(𝑥), (145a)(145b)(145c)where 𝑢, 𝑣, 𝑤 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𝜀_𝑥𝑥=𝜃_𝑥^′′𝜙^𝑃_𝑆,𝛾_𝑥𝑦=𝜃^′_𝑥𝜕𝜙^𝑃_𝑆𝜕𝑦−𝑧+𝜃_𝑥^′′′𝜕𝜙^𝑆_𝑆,𝛾𝜕𝑦_𝑥𝑧=𝜃^′_𝑥𝜕𝜙^𝑃_𝑆𝜕𝑧+𝑦+𝜃_𝑥^′′′𝜕𝜙^𝑆_𝑆,𝜕𝑧 (146a)(146b)(146c)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 𝑑𝜎_𝑥𝑥𝑑𝜏_𝑥𝑦𝑑𝜏_𝑥𝑧=𝐸^∗000𝐺000𝐺𝑑𝜀^el_𝑥𝑥𝑑𝛾^el_𝑥𝑦𝑑𝛾^el_𝑥𝑧, (147) where 𝑑(⋅) denotes infinitesimal incremental quantities over time (rates), the superscript el denotes the elastic part of the strain components and 𝐸^∗=𝐸(1−𝑣)/(1+𝑣)(1−2𝑣). If the plane stress hypothesis is undertaken then 𝐸^∗=𝐸/(1−𝑣²) 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, 𝑑𝜀_𝑥𝑥𝑑𝛾_𝑥𝑦𝑑𝛾_𝑥𝑧^𝑇=𝑑𝜀^el_𝑥𝑥𝑑𝛾^el_𝑥𝑦𝑑𝛾^el_𝑥𝑧^𝑇. (148) The stress rates are given with respect to the total strain ones by combining (147) and (148). If plastic flow occurs then 𝑑𝜀_𝑥𝑥𝑑𝛾_𝑥𝑦𝑑𝛾_𝑥𝑧^𝑇=𝑑𝜀^el_𝑥𝑥𝑑𝛾^el_𝑥𝑦𝑑𝛾^el_𝑥𝑧^𝑇+𝑑𝜀^pl_𝑥𝑥𝑑𝛾^pl_𝑥𝑦𝑑𝛾^pl_𝑥𝑧^𝑇, (149) 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 𝑓=𝜎²_𝑥𝑥+3𝜏²_𝑥𝑦+𝜏²_𝑥𝑧−𝜎_𝑌𝜀^pl_eq=0, (150) 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 𝑑𝜀^pl_𝑥𝑥𝑑𝛾^pl_𝑥𝑦𝑑𝛾^pl_𝑥𝑧^𝑇=𝑑𝜆𝜕𝑓𝜕𝜎_𝑥𝑥𝜕𝑓𝜕𝜏_𝑥𝑦𝜕𝑓𝜕𝜏_𝑥𝑧^𝑇. (151) 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 𝑑𝜎_𝑥𝑥𝑑𝜏_𝑥𝑦𝑑𝜏_𝑥𝑧=𝐷^el-pl𝑑𝜀_𝑥𝑥𝑑𝛾_𝑥𝑦𝑑𝛾_𝑥𝑧, (152) where [𝐷^el-pl] is the elastoplastic constitutive matrix given as 𝐷^el-pl=1𝑐𝑐₁₁𝑐sym.₂₁𝑐₂₂𝑐₃₁𝑐₃₂𝑐₃₃ (153) with𝑐=ℎ𝜎²_𝑒+𝐸𝜎²_𝑥𝑥+9𝐺𝜏²_𝑥𝑦+𝜏²_𝑥𝑧𝑐₁₁,=𝐸ℎ𝜎²_𝑒+9𝐺𝜏²_𝑥𝑦+𝜏²_𝑥𝑧𝑐₂₁,=−3𝐸𝐺𝜎_𝑥𝑥𝜏_𝑥𝑦,𝑐₃₁=−3𝐸𝐺𝜎_𝑥𝑥𝜏_𝑥𝑧,𝑐₂₂=𝐺ℎ𝜎²_𝑒+𝐸𝜎²_𝑥𝑥+9𝐺𝜏²_𝑥𝑧𝑐₃₂,=−9𝐺²𝜏_𝑥𝑦𝜏_𝑥𝑧,𝑐₃₃=𝐺ℎ𝜎²_𝑒+𝐸𝜎²_𝑥𝑥+9𝐺𝜏²_𝑥𝑦,𝜎_𝑒=𝜎²_𝑥𝑥+3𝜏²_𝑥𝑦+𝜏²_𝑥𝑧. (154a)(154b)(154c)(154d)(154e)(154f)(154g)(154h)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.

(a)

(b)

Figure 16 

Normal stress-strain (a) and yield stress-equivalent plastic strain (b) relationships.

9.1. 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∇²𝜙^𝑃_𝑆=0inΩ,𝜕𝜙^𝑃_𝑆𝜕𝑛=𝑧𝑛_𝑦−𝑦𝑛_𝑧onΓ_𝑗, (155a)(155b)

∇²𝜙^𝑆_𝑆𝐸=−𝐺𝜙^𝑃_𝑆𝑦,𝑧inΩ,𝜕𝜙^𝑆_𝑆𝜕𝑛=0onΓ_𝑗 (156a)(156b) 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.

9.2. Equations of Global Equilibrium

To establish global equilibrium equations, the principle of virtual work neglecting body forces is employed as follows: _𝑉𝜎_𝑥𝑥𝛿𝜀_𝑥𝑥+𝜏_𝑥𝑦𝛿𝛾_𝑥𝑦+𝜏_𝑥𝑧𝛿𝛾_𝑥𝑧=𝑑𝑉_𝐹𝑡_𝑥𝛿𝑢+𝑡_𝑦𝛿𝑣+𝑡_𝑧𝛿𝑤𝑑𝐹, (157) 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 𝐸𝐶_𝑆𝜃_𝑥^′′′′−𝐺𝐼_𝑡𝜃_𝑥^′′=𝑚_𝑡(𝑥)+𝑑𝑚_𝑤(𝑥)𝑑𝑥 (158) 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 _𝑉𝜎_𝑥𝑥+Δ𝜎_𝑥𝑥𝛿Δ𝜀_𝑥𝑥+𝜏_𝑥𝑦+Δ𝜏_𝑥𝑦𝛿Δ𝛾_𝑥𝑦+𝜏_𝑥𝑧+Δ𝜏_𝑥𝑧𝛿Δ𝛾_𝑥𝑧=𝑑𝑉_𝐹𝑡_𝑥𝛿Δ𝑢+𝑡_𝑦𝛿Δ𝑣+𝑡_𝑧𝛿Δ𝑤𝑑𝐹 (159) 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 well-known 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)𝑆𝑀_𝑡=_Ω𝜏_𝑥𝑦𝜕𝜙^𝑃_𝑆𝜕𝑦−𝑧+𝜏_𝑥𝑧𝜕𝜙^𝑃_𝑆𝜕𝑧+𝑦𝑑Ω,𝑆𝑀_𝑤=_Ω𝜎_𝑥𝑥𝜙^𝑃_𝑆𝑑Ω, (160a)(160b)where 𝑆𝑀_𝑡 and 𝑆𝑀_𝑤 correspond to the internal twisting and warping moments of the bar, respectively, and the relevant incremental stress resultants are defined asΔ𝑆𝑀_𝑡=_ΩΔ𝜏_𝑥𝑦𝜕𝜙^𝑃_𝑆𝜕𝑦−𝑧+Δ𝜏_𝑥𝑧𝜕𝜙^𝑃_𝑆𝜕𝑧+𝑦𝑑Ω,Δ𝑆𝑀_𝑤=_ΩΔ𝜎_𝑥𝑥𝜙^𝑃_𝑆𝑑Ω. (161a)(161b)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 −𝑑Δ𝑆𝑀_𝑡+𝑑𝑑𝑥²Δ𝑆𝑀_𝑤𝑑𝑥²_{tangentstiffnessmatrix×incrementaldisplacements}=𝑚_𝑡(𝑥)+^{out-of-balanceforces}𝑑𝑚_𝑤(𝑥)𝑑𝑥_{externallyappliedloading}+𝑑𝑆𝑀_𝑡−𝑑𝑑𝑥²𝑆𝑀_𝑤𝑑𝑥²_{internalstressresultants} (162) along with its corresponding boundary conditions𝛼₁Δ𝑆𝑀_𝑡−𝑑Δ𝑆𝑀_𝑤𝑑𝑥+𝛼₂Δ𝜃_𝑥=𝛼₃−𝛼₁𝑆𝑀_𝑡−𝑑𝑆𝑀_𝑤,𝛽𝑑𝑥₁Δ𝑆𝑀_𝑤+𝛽₂Δ𝜃^′_𝑥=𝛽₃−𝛽₁𝑆𝑀_𝑤, (163a)(163b)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 𝛼₂=𝛽₂=1, 𝛼₁=𝛼₃=𝛽₁=𝛽₃=0). 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𝑚_𝑡(𝑥)=_Γ𝑡_𝑦(−𝑧)+𝑦𝑡_𝑧𝑚𝑑𝑠,_𝑤(𝑥)=−_Γ𝑡_𝑥𝜙^𝑃_𝑆𝑑𝑠. (164a)(164b)

From the above presentation and the examined example problems presented in [134] it is concluded that (i)both St. Venant and warping shear stresses are evaluated from the solution of two boundary value problems formulated under elastic conditions;(ii)warping shear stresses slightly decrease the torsional rigidity and the ultimate torsional loading that bars can undertake under nonuniform torsion;(iii)warping shear stresses have an influence on inelastic redistribution patterns of bars under nonuniform torsion.

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

As it has been already mentioned in the previous section, material nonlinearity is important for investigating the ultimate strength of a bar that resists torsional loading, while distributed plasticity models are acknowledged in the literature [114–116] to capture more rigorously material nonlinearities than cross sectional stress resultant approaches [117] or lumped plasticity idealizations [118, 119]. Moreover, in the nonuniform torsional analysis of a bar except from the primary (St. Venant) shear stress distribution forming the primary torsional moment stress resultant additional normal and secondary (warping) shear stresses arise, forming the warping moment and secondary torsional moment stress resultants, respectively. In order to include warping shear stresses in the global equilibrium of the bar, that is, to account for the secondary torsional moment deformation effect (STMDE), an additional kinematical component (along with the angle of twist) is generally required [81–140], increasing the difficulty of the problem at hand.

The STMDE has been shown in the literature to be significant, especially on closed shaped section bars. Massonnet [141] presented a qualitative explanation why warping shear stresses are of the same order of magnitude as primary ones in the case of closed shaped section bars. As early as 1954, Benscoter [142] analyzed the nonuniform torsional problem of multicell section bars. Since then, a significant amount of relevant contributions has appeared in the literature as well [109, 140, 143–146]. Since the topic at hand is analogous to the geometrically nonlinear Timoshenko beam theory of shear-bending loading conditions [78, 83], it does not satisfy local equilibrium equations (for relevant discussions see e.g., Simo et al. [132] and Minghini et al. [133]). This problem is alleviated by introducing torsional shear correction factor at the global level [78, 84, 85] and suitable warping shear stress distribution at the local level [80, 85, 147]. The aforementioned contributions refer to the linear elastic regime. If inelastic effects are considered, especially through distributed plasticity formulations, then the nonuniform [120, 122] torsion problem including STMDE requires a much more rigorous analysis.

In the majority of the research efforts concerning the inelastic nonuniform torsional problem of bars (already presented in the previous section) torsional shear correction factor is not included in the analyses, while the employed warping shear stress distributions though included in the global equilibrium of the bar they do not satisfy local equilibrium considerations under inelastic or even elastic conditions. A different approach is undertaken in the recent contribution of Sapountzakis and Tsipiras [139], where a single (fourth-order) governing differential equation is formulated with respect to a single kinematical component (angle of twist), leading to an approximate solution of the problem at hand. However, the adopted warping shear stress distribution verifies local equilibrium equations under elastic conditions. Moreover, an alternative methodology is presented in the very recent contribution of Le Corvec and Filippou [148], where a multitude of local section warping degrees of freedom are introduced in order to model effects such as in-plane plasticity effects, constrained warping, shear lag, and so forth. This technique does not require the introduction of torsional shear correction factor; however, the reduction of the number of employed degrees of freedom requires further investigation, while an example of a thick-walled cross section beam is not presented. In the publication of Wackerfuß and Gruttmann [149], beams of thick-walled rectangular cross sections are worked out by employing a series of polynomials as global warping functions in order to capture advanced effects such as in-plane inelastic redistribution and transverse contraction. Moreover, the same researchers [150] have employed local warping functions to capture the aforementioned effects in moderately thick arbitrarily shaped cross section beams; however, bimoment loading cannot be applied at the bar. Finally, Tsipiras and Sapountzakis [151] developed a boundary element solution for the inelastic nonuniform torsional problem of simply or multiply connected cylindrical bars of arbitrarily shaped doubly symmetric cross section taking into account the effect of secondary torsional moment deformation. The bar is subjected to arbitrarily distributed or concentrated torsional loading along its length, while its edges are subjected to the most general torsional boundary conditions. A displacement based formulation is developed and inelastic redistribution is modelled through a distributed plasticity model exploiting three dimensional material constitutive laws and numerical integration over the cross sections. An incremental-iterative solution strategy is adopted to resolve the elastic and plastic part of stress resultants along with an efficient iterative process to integrate the inelastic rate equations [135]. The one dimensional primary angle of twist per unit length, a two dimensional secondary warping function, and a scalar torsional shear correction factor are employed to account for STMDE. The latter is computed employing an energy approach under elastic conditions [85].

In order to formulate the aforementioned problem, let us consider the prismatic bar of 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𝑢𝑥,𝑦,𝑧=𝜃^𝑃_𝑥(𝑥)^′𝜙^𝑃_𝑆,𝑣()𝑦,𝑧𝑥,𝑧=−𝑧𝜃_𝑥(𝑥),𝑤𝑥,𝑦=𝑦𝜃_𝑥(𝑥), (165a)(165b)(165c)where 𝑢, 𝑣, 𝑤 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 well-known three-dimensional linear strain-displacement relations, the nonvanishing compatible (total) strain resultants are obtained as𝜀_𝑥𝑥=𝜃^𝑃_𝑥^′′𝜙^𝑃_𝑆,𝛾_𝑥𝑦=𝜃^′_𝑥𝜕𝜙^𝑃_𝑆𝜕𝑦−𝑧_primary+𝜃^𝑃_𝑥^′−𝜃^′_𝑥𝜕𝜙^𝑃_𝑆𝜕𝑦_secondary,𝛾_𝑥𝑧=𝜃^′_𝑥𝜕𝜙^𝑃_𝑆𝜕𝑧+𝑦_primary+𝜃^𝑃_𝑥^′−𝜃^′_𝑥𝜕𝜙^𝑃_𝑆𝜕𝑧_secondary, (166a)(166b)(166c)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𝜀_𝑥𝑥=𝜃^𝑃_𝑥^′′𝜙^𝑃_𝑆,𝛾_𝑥𝑦=𝜃^′_𝑥𝜕𝜙^𝑃_𝑆𝜕𝑦−𝑧_primary+𝜃^𝑃_𝑥^′−𝜃^′_𝑥√𝜅_𝑥𝜕𝜙^𝑃_𝑆𝜕𝑦_secondary,𝛾_𝑥𝑧=𝜃^′_𝑥𝜕𝜙^𝑃_𝑆𝜕𝑧+𝑦_primary+𝜃^𝑃_𝑥^′−𝜃^′_𝑥√𝜅_𝑥𝜕𝜙^𝑃_𝑆𝜕𝑧_secondary, (167a)(167b)(167c)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𝜀_𝑥𝑥=𝜃^𝑃_𝑥^′′𝜙^𝑃_𝑆,𝛾_𝑥𝑦=𝜃^′_𝑥𝜕𝜙^𝑃_𝑆𝜕𝑦−𝑧_primary+𝜃^𝑃_𝑥^′−𝜃^′_𝑥−𝐼^𝑆_𝑡𝐶_𝑆𝜕𝜙^𝑆_𝑆𝜕𝑦_secondary,𝛾_𝑥𝑧=𝜃^′_𝑥𝜕𝜙^𝑃_𝑆𝜕𝑧+𝑦_primary+𝜃^𝑃_𝑥^′−𝜃^′_𝑥−𝐼^𝑆_𝑡𝐶_𝑆𝜕𝜙^𝑆_𝑆𝜕𝑧_secondary, (168a)(168b)(168c)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:𝑑𝑆𝑀^𝑃_𝑡+𝑑𝑥𝑑𝑆𝑀^𝑆_𝑡𝑑𝑥=−𝑚_𝑡,𝑑𝑆𝑀_𝑤𝑑𝑥+𝑆𝑀^𝑆_𝑡=𝑚_𝑤 (169a)(169b)along with their corresponding boundary conditions𝛼₁𝑆𝑀^tot_𝑡+𝛼₂𝜃_𝑥=𝛼₃,𝛽₁𝑆𝑀_𝑤+𝛽₂𝜃^𝑃_𝑥^′=𝛽₃, (170a)(170b)where 𝑆𝑀^tot_𝑡=𝑆𝑀^𝑃_𝑡+𝑆𝑀^𝑆_𝑡 is the total twisting moment at the bar ends, while the stress resultants are naturally defined as𝑆𝑀_𝑤=_Ω𝜎_𝑥𝑥𝜙^𝑃_𝑆𝑑Ω,𝑆𝑀^𝑃_𝑡=_Ω𝜏_𝑥𝑦𝜕𝜙^𝑃_𝑆𝜕𝑦−𝑧+𝜏_𝑥𝑧𝜕𝜙^𝑃_𝑆𝜕𝑧+𝑦𝑑Ω,𝑆𝑀^𝑆_𝑡=𝐼^𝑆_𝑡𝐶_𝑆Ω𝜏_𝑥𝑦𝜕𝜙^𝑆_𝑆𝜕𝑦+𝜏_𝑥𝑧𝜕𝜙^𝑆_𝑆𝜕𝑧𝑑Ω, (171a)(171b)(171c)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𝑚_𝑡(𝑥)=_Γ𝑡_𝑦(−𝑧)+𝑦𝑡_𝑧𝑚𝑑𝑠,_𝑤(𝑥)=−_Γ𝑡_𝑥𝜙^𝑃_𝑆𝑑𝑠. (172a)(172b)

Since an incremental-iterative approach is adopted for the problem at hand, the incremental version of (169a) and (169b) is written down as𝑑Δ𝑆𝑀^𝑃_𝑡+𝑑𝑥𝑑Δ𝑆𝑀^𝑆_𝑡𝑑𝑥=−Δ𝑚_𝑡,𝑑Δ𝑆𝑀_𝑤𝑑𝑥+Δ𝑆𝑀^𝑆_𝑡=Δ𝑚_𝑤, (173a)(173b)where Δ(⋅) denotes incremental quantities (over time), while the incremental stress resultants are given by virtue of (171a), (171b), (171c), (147), and (149) asΔ𝑆𝑀_𝑤=𝐸𝐶_𝑆Δ𝜃^𝑃_𝑥^′′+Δ𝑆𝑀^pl_𝑤,Δ𝑆𝑀^𝑃_𝑡=𝐺𝐼^𝑃_𝑡Δ𝜃^′_𝑥+Δ𝑆𝑀^𝑃pl_𝑡,Δ𝑆𝑀^𝑆_𝑡=−𝐺𝐼^𝑆_𝑡Δ𝜃^𝑃_𝑥^′−Δ𝜃^′_𝑥+Δ𝑆𝑀^𝑆pl_𝑡, (174a)(174b)(174c)where 𝐼^𝑃_𝑡, 𝐶_𝑆 are the primary torsion (St. Venant) and warping constants, defined from (72b) and (72a), respectively, while 𝐼^𝑆_𝑡 is the secondary torsion constant given as 𝐼^𝑆_𝑡=𝜅_𝑥Ω𝜕𝜙^𝑃_𝑆𝜕𝑦²+𝜕𝜙^𝑃_𝑆𝜕𝑧²𝑑Ω(modelC) (175) and Δ𝑆𝑀^pl_𝑤, Δ𝑆𝑀^𝑃pl_𝑡, Δ𝑆𝑀^𝑆pl_𝑡 are plastic quantities defined asΔ𝑆𝑀^pl_𝑤=−𝐸_ΩΔ𝜀^pl_𝑥𝑥𝜙^𝑃_𝑆𝑑Ω,Δ𝑆𝑀^𝑃pl_𝑡=−𝐺_ΩΔ𝛾^pl_𝑥𝑦𝜕𝜙^𝑃_𝑆𝜕𝑦−𝑧+Δ𝛾^pl_𝑥𝑧𝜕𝜙^𝑃_𝑆𝜕𝑧+𝑦𝑑Ω,Δ𝑆𝑀^𝑆pl_𝑡𝐼=−𝐺^𝑆𝑡𝐶_𝑆ΩΔ𝛾^pl_𝑥𝑦𝜕𝜙^𝑆_𝑆𝜕𝑦+Δ𝛾^pl_𝑥𝑧𝜕𝜙^𝑆_𝑆()𝜕𝑧𝑑ΩmodelC (176a)(176b)(176c)while the expressions of Δ𝑆𝑀^𝑆_𝑡^pl and 𝐼^𝑆_𝑡 for models A, B and the details of the derivation of expression (174c) are presented in [151]. Substituting (174a), (174b), and (174c) in (173a), and (173b) and forming the incremental version of the boundary conditions (170a) and (170b), the following boundary value problem is obtained:−𝐺𝐼^𝑃_𝑡+𝐼^𝑆_𝑡Δ𝜃_𝑥^′′+𝐺𝐼^𝑆_𝑡Δ𝜃^𝑃_𝑥^′′=Δ𝑚_𝑡+𝑑Δ𝑆𝑀^𝑃pl_𝑡+𝑑𝑥𝑑Δ𝑆𝑀^𝑆pl_𝑡𝑑𝑥insidethebar−𝐸𝐶_𝑆Δ𝜃^𝑃_𝑥^′′′+𝐺𝐼^𝑆_𝑡Δ𝜃^𝑃_𝑥^′−Δ𝜃^′_𝑥=−Δ𝑚_𝑤+𝑑Δ𝑆𝑀^pl_𝑤𝑑𝑥+Δ𝑆𝑀^𝑆pl_𝑡insidethebar, (177a)(177b)

𝛼₁Δ𝑆𝑀^tot_𝑡+𝛼₂Δ𝜃_𝑥=Δ𝛼₃𝛽atthebarends𝑥=0,𝑙,₁Δ𝑆𝑀_𝑤+𝛽₂Δ𝜃^𝑃_𝑥^′=Δ𝛽₃atthebarends𝑥=0,𝑙. (178a)(178b) 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−𝐺𝐼^𝑃_𝑡+𝐼^𝑆_𝑡Δ𝜃_𝑥^′′+𝐺𝐼^𝑆_𝑡Δ𝜂^′_𝑥=Δ𝑚_𝑡+𝑑Δ𝑆𝑀^𝑃pl_𝑡+𝑑𝑥𝑑Δ𝑆𝑀^𝑆pl_𝑡𝑑𝑥insidethebar−𝐸𝐶_𝑆Δ𝜂_𝑥^′′+𝐺𝐼^𝑆_𝑡Δ𝜂_𝑥−Δ𝜃^′_𝑥=−Δ𝑚_𝑤+𝑑Δ𝑆𝑀^pl_𝑤𝑑𝑥+Δ𝑆𝑀^𝑆pl_𝑡insidethebar, (179a)(179b)

𝛼₁Δ𝑆𝑀^tot_𝑡+𝛼₂Δ𝜃_𝑥=Δ𝛼₃𝛽atthebarends𝑥=0,𝑙,₁Δ𝑆𝑀_𝑤+𝛽₂Δ𝜂_𝑥=Δ𝛽₃atthebarends𝑥=0,𝑙. (180a)(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.

11. Concluding Remarks

In this paper both the static and dynamic analyses of the geometrically linear or nonlinear, elastic or elastoplastic nonuniform torsion problems of bars of constant or variable arbitrary cross section are presented together with the corresponding research efforts and the conclusions drawn from examined cases with great practical interest. In the presented analyses, the bar is subjected to arbitrarily distributed or concentrated twisting and warping moments along its length, while its edges are supported by the most general torsional boundary conditions. The analysis of the aforementioned problems is complete by presenting the evaluation of the torsion and warping constants of the bar, its displacement field, its stress resultants together with the torsional shear stresses and the warping normal and shear stresses at any internal point of the bar.

Having in mind the disadvantages of the 3D FEM solutions and more specifically the difficulties(i)in support modeling,(ii)in discretizing a complex structure,(iii)in discretizing a structure including thin-walled members (shear locking and membrane locking),(iv)in the increased number of degrees of freedom leading to severe or unrealistic computational time,(v)in the reduced oversight of the 3D FEM solution compared with that of the beamlike structures employing stress resultants,and the fact that the use of shell elements cannot give accurate results since the warping of the walls of a cross section cannot be taken into account (midline model), the importance of the presented beamlike analyses becomes more evident.