Displacement and Stress Function-based Linear and Quadratic Triangular Elements for Saint-Venant Torsional Problems

Torsional problems commonly arise in frame structural members subjected to unsymmetrical loading. Saint-Venant proposed a semi inverse method to develop the exact theory of torsional bars of general cross sections. However, the solution to the problem using an analytical method for a complicated cross section is cumbersome. This paper presents the adoption of the Saint-Venant theory to develop a simple finite element program based on the displacement and stress function approaches using the standard linear and quadratic triangular elements. The displacement based approach is capable of evaluating torsional rigidity and shear stress distribution of homogeneous and nonhomogeneous; isotropic, orthotropic, and anisotropic materials; in singly and multiply-connected sections. On the other hand, applications of the stress function approach are limited to the case of singly-connected isotropic sections only, due to the complexity on the boundary conditions. The results show that both approaches converge to exact solutions with high degree of accuracy.


Introduction
Torsional problems commonly arise in three dimensional structural members subjected to unsymmetrical loading or internal twisting moment exerted to satisfy local equilibrium in members' connections [1].The twisting moment causes rotational deformation of cross section; which values are highly attributed to the resistance of cross section against twisting i.e. torsional rigidity.In many cases in steel and reinforced concrete designs, torsion, if presents, can be the governing design factors among others.Therefore, it is essential to correctly evaluate the torsional rigidity of any cross sections, such that the design conforms to torsional strength requirements.
Torsional rigidity of elastic circular sections can be evaluated using conventional mechanics and elasticity since the sections along member's length remain plain and undistorted [2,3].However, in a noncircular cross sectional member, the sections along its length warp due to various twisting.Applying the torsion theory for circular to non-circular sections results in violation of the equilibrium equations and boundary conditions.Saint-Venant proposed a semiinverse approach to solve torsional problems by assuming the unknown displacements to fulfill the equilibrium equations and boundary conditions.Nevertheless, evaluating the torsional rigidity of noncircular cross section using analytical methods is very cumbersome due to complex boundary conditions.The complexity of boundary conditions is found either in complex singly-connected cross section's shape or multiply-connected cross section.A singly-connected cross section means that the region of the cross section is bounded by a closed, piecewise smooth curve.A multiply-connected cross section is bounded by several closed, piecewise smooth curves, such as a section with several holes.The finite element method (FEM) can be used to approximate the torsional rigidity of arbitrary cross sections due to its flexibility to discretize the domain, i.e. the cross section, and thus reduces the complexities of boundary conditions.Moreover, accurate FEM results can be obtained by refining the mesh or introducing higher degree of element's shape functions [4,5].Numerous available approaches such as displacement, stress, hybrid, and mix formulations have been employed in FEM-based analysis, including for torsional problems.This research aims to develop a simple FEM-based program in MATLAB to deal with elastic torsional problems of noncircular cross sections, which formulation is based on the displacement and stress approaches.Various material types, including isotropic, orthotropic, and anisotropic, are considered in the model.Furthermore, the model also considers both singly and multiply-connected cross sections.The main outputs of the model involve torsional rigidity and shear stress contour.The outputs are further compared to the results available in literature [2,6,7].

Saint-Venant Torsion Theory
Suppose a cantilever bar of the length h with general cross section is fixed on one end and is subjected to a twisting moment on the other free end as depicted in Figure 1.Saint-Venant suggested that the deformations consist of rotation and warping [2,6].Furthermore, the cross section rotates as a rigid body motion when twisted so that no in-plane strain presents.The warping is assumed to be constant along the bar axis so that normal stresses also vanish.Hence, the remaining strains are γxz and γyz, which are respectively the shear strains of the xz-plane and yzplane and thus leaving the bar under pure shear.

Displacement Function Approach
According to Saint-Venant torsion theory [6], the displacement fields in the bar (Figure 1), can be expressed as: Here, variables u, v, and w denote displacement components in the x, y, and z-directions, respectively (Figure 1 (a)), θ denotes the twist angle per unit length (Figure 1 (c)), x, y, and z are the coordinate of the material points under consideration, and ψ(x, y) represents the displacement function.
Following the small strain theory of elasticity, the non-zero stress components at a generic point in the bar cross section can be written as In these equations, τxz and τyz are the torsional shear stress components, whereas γxz and γyz are the corresponding shear strain components (Figure 1 (b)).Matrix [G] represents constitutive material matrix, the values of which depends on the type of materials, namely, isotropic, orthotropic, or anisotropic.For example, for an isotropic material G11 = G22 = G, that is, the shear modulus, and G12 = G21 = 0.The equilibrium equation of an elemental cube in the z direction is given as The governing differential equation for the Saint-Venant torsion in the domain Ω of a cross section can be obtained by substituting Equation ( 2) and (3) into Equation.(4).For the case of a bar made from homogeneous and isotropic material, the governing equation becomes [6] Meanwhile, the boundary conditions at the perimeter of cross section, Γ, is given as [6] where s is the boundary curve parameter (Figure 1 (b)).The total potential energy of the bar, neglecting the twisting at the member's end, can be written as [6]      where h is the bar length.The unknown displacement function, ψ, can be obtained by employing the principle of minimum potential energy, by setting the first variation of Πp with respect to ψ equal to zero.

Stress Function Approach
An alternative approach to solve the Saint-Venant torsional problem is to introduce a stress function, φ.
In this approach, the shear stress components are expressed in terms of the stress function as in Equation ( 8) [6].Accordingly, the governing equation within the domain Ω and the corresponding boundary conditions at the surface, Γ, can be written as in Equation ( 9).Identical to the displacement approach, the problem can be solved by targeting the first variation of the complementary potential energy, Equation (10), equal to zero [6].In the present research, this approach is applicable only to singlyconnected section with isotropic material due to limitation of boundary conditions.

Finite Element Formulation Displacement Function Approach
The assumed displacement function, ψ, is approximated as given in Equation (11).Matrix [N] in Equation ( 12) is the matrix of element shape functions, {qψ} is the vector of displacement function values at the finite element nodes, and n is the number of finite element nodes.The standard linear three-node (n = 3) and quadratic six-node (n = 6) triangular elements are employed in the present research.
The assumed displacement function, ψ, is further partially differentiated with respect to x and y to obtain the strain-warping transformation matrix [B] as defined in Equation (13).This equation is then substituted into Equation (3) to obtain Equation (14).Finally, substituting Equation (14) into Equation (7) and employing the Rayleigh-Ritz method, the element stiffness-displacement-load vector relationship is obtained as given in Equations ( 15) and ( 16).Here, [kψ] denotes the element stiffness matrix, [qψ] denotes the nodal displacement vector, and [Qψ] denotes the nodal load vector.
The relationship between the twisting moment Mt and torsional constant is given in Equation ( 17), where J is torsional constant and D is torsional rigidity.Furthermore, twisting moment Mt can be generally expressed in the integral form as in Equation ( 18).Substituting Equation (3) into Equation (18), the torsional constant can be calculated via Equation (19).Similarly, the torsional rigidity can be evaluated through Equation (20).

Stress Function Approach
Similar to the formulation of the displacement approach, the assumed stress function, , can be approximated as given in Equation ( 21).The assumed stress function is further differentiated with respect to x and y to obtain the relationship in Equation ( 22).Substituting Equation ( 22) into Equation (10) and applying the principle of minimum complementary energy field Equations ( 23) and ( 24), where [kφ] denotes the element stiffness, [qφ] denotes the vector of nodal stress function values, and [Qφ] denotes nodal load vector.
Applying similar process as in the displacement function approach based on Equations ( 17) and (18), the torsional constant J and torsional rigidity D based on stress function approach can be evaluated using Equation ( 25) and (26), respectively.
A FEM program has been developed based on the displacement and stress function approach to deal with Saint-Venant torsional problems on various cross sections.The program was verified and tested using several numerical test presented in the following section.

Results and Discussions
The results in terms of torsional rigidity as well as shear stress distribution are presented based on the types of material used, i.e. homogeneous and nonhomogeneous.Homogeneous material covers isotropic rectangular section, orthotropic elliptical section, anisotropic elliptical section, and isotropic multiplyconnected section.Meanwhile, for nonhomogeneous material, the results of a square section are presented.Furthermore, an evaluation of torsional constant of wide flange steel section is also considered.To evaluate the convergence characteristics of the developed finite elements, five levels of meshes, i.e. coarse, medium-coarse, medium, medium-fine, and fine, are used in each problem.The first generated mesh via MATLAB PDE Toolbox is referred to as coarse mesh.Subsequent meshes are referred to as medium-coarse, medium, medium-fine, and fine meshes.

Problem 1: Homogeneous Isotropic Square Section
Table 1 shows the torsional rigidity of homogeneous isotropic square section having 2 by 2-unit length.In the problem, the shear modulus G and the angle of twist θ were assumed to be unity.It can be observed that the results converge faster when higher order elements are used.To obtain under 1% accuracy, the displacement function needs medium and mediumcoarse mesh when using linear and quadratic elements, respectively.Meanwhile, the stress function needs medium-fine and medium-coarse mesh employing linear and quadratic elements, respectively.It is important to note that the results are better when compared to the previous research [8].By comparing the error introduced by displacement function approach to stress function approach, it is notable that the rate of convergence of the preceding approach is 2.67 time faster than the later approach when linear elements are occupied.However, when the quadratic elements are adopted, the rate of convergence of the stress function increases dramatically to 0.75 times the rate of displacement function.Overall, the results obtained using both approaches converge well to the exact solution when the mesh is refined.Meanwhile, Figure 2

Problem 2: Homogeneous Orthotropic Elliptical Section
A homogeneous orthotropic elliptic bar with semimajor axis of 20-unit length and semi-minor axis of 10-unit length under unity end twisting moment was analyzed.The property of orthotropic material is taken from literature [7], that is, G11 = 1 and G22 = 8.The analysis was done using only displacement function approach because in the present formulation of the stress function approach is only applicable for the cases of cross sections made from homogeneous material.
The distribution of shear stress across the section on yz-plane is depicted in Figure 3, which shows that there is an abrupt change in shear stress distribution along yz-plane when linear elements are used.However, this phenomenon vanishes as higher order elements are adopted.
The resulting torsional rigidity is presented in Table 2.It can be seen that the linear elements give the upper bound solution, while the quadratic elements result in lower bound value [6].It also shows that the Exact solution = 2.24923 [6], e = error medium-fine mesh is needed for the linear element to achieve the accuracy under 1%, while the medium mesh is already sufficient when the quadratic element is employed.Regardless the stress distribution, the results stably converge to the exact solution.

Problem 3: Homogeneous Anisotropic Elliptical Section
The identical elliptic section as in the previous problem but using a homogeneous anisotropic material is analyzed using the displacement function approach.The anisotropic material property is given as G11 = 1, G12 = G21 = 2, and G22 = 8 [7]. Figure 4 shows the shear stress contour across the section, which displays an unsmooth contour transition for yz-plane when linear elements are used.Proper shear stress distribution can be achieved when higher order elements are used.Similar phenomenon occurs for shear stress distribution for xzplane.Furthermore, obtaining highly accurate result using linear elements requires high demand in computational resources because 1% accuracy can be obtained only if fine mesh is used, as indicated in Table 4.In contrast, quadratic elements can yield such accuracy using only the medium mesh.However, even the stress contour shows discontinueties, the results show stable convergence as the meshes are refined.

Problem 4: Homogeneous Isotropic Multiplyconnected Section
Figure 5 shows the shear stress distribution of homogeneous isotropic multiply-connected section having inner radius of 1-unit length and outer radius of 3-unit length in polar coordinate.The shear modulus (G) and the angle of twist (θ) are assumed to be unity [7].It can be observed that the shear stress distribution is in agreement with the theoretical result.

Figure 5. Shear Stress Contour of Homogeneous Isotropic Multiply-connected Section
Table 4 demonstrates the high accuracy of the finite element results, that is, an accuracy below 1% can be achieved using medium-coarse meshing.The results also indicate that the program is able to yield a good convergence.Another notable finding is that the results' convergence is highly affected by the number of elements rather than the order of elements.This can be caused by the presence of inner hole, which constraints the meshing to keep the quality of elements.Furthermore, for multiply-connected section, only displacement function approach is considered.The stress function approach is currently not considered because the complexity on boundary conditions.[2], e = error Benchmark solution = 122.8(2.28%) [9], 125.6 (0.05%) [7] Problem 5: Nonhomogeneous Square Section Figure 6 shows the shear stress distribution of 2 by 2-unit length nonhomogeneous square section.The compound section is composed of two materials with relative ratio of shear modulus equals to 2.0, i.e.G1 = 1 and G2 = 2 for the right-half and the left-half of the section, respectively.The angle of twist is assumed to be unity.Observing the shear stress distribution, it is notable that the shear stress at the left perimeter is higher than at the right, which in agreement to the results in literature [10].Another key finding is the stress jump at the boundary of two distinct materials.Furthermore, Table 7 reveals that the results stably converge to the benchmark solution.

Application to Steel Profiles W36x256
Table 6 shows the torsional constant of W36x256 steel section, which was obtained using displacement and stress function approach.Shear modulus G = 29,000 ksi and twisting moment M = 150 lb-in were used in the analysis.The results reveal that the torsional constant converge to 49.8 in 4 , which differs to the benchmark solution, i.e. 53.3 in 4 .This can be because the model neglects the radial areas around the web-flange corners, which definitely reduces the torsional constant.

Meshing Requirements
Table 7 shows the mesh density rating to achieve results of 1% accuracy compared to the benchmark solutions.It is seen that to achieve such accuracy, the displacement function approach needs medium to medium-fine mesh when linear elements are used.Whilst, medium to medium-coarse mesh is sufficient when quadratic elements are used.The result's accuracy obtained from the stress function approach introduces notable error when linear elements are used.Nevertheless, the approach shows significant improve when quadratic elements are adopted.

Conclusions
A FEM program based on the displacement and stress function approaches has been developed to solve the Saint-Venant torsion problems.Both approaches were found to be reliable to solve torsional problems on various geometries and material types.The results show that both approaches yield stable numerical results in the sense that the use of a finer mesh as well as a higher order element will give a smaller error.Nonetheless, unlike the displacement approach, the stress function approach is currently limited to evaluate singly-connected isotropic section only due to the complexity on the boundary conditions.Furthermore, the use of linear elements in the displacement and stress function approaches require medium to medium-fine and medium-fine mesh size, respectively.Meanwhile, when quadratic elements are used, the stress function approach outperforms the displacement function approach by requiring only medium-coarse mesh size compared to medium to medium-coarse mesh size of the displacement's.Furthermore, the use of linear elements in orthotropic and anisotropic material causes shear stress contour jump.The phenomenon vanishes when higher degree of elements is employed.Another shear stress jump occurs at the connection of two distinct materials.Special treatment must be made to eliminate the discrepancy.Further research can be conducted to develop the program such that it can deal with multiply-connected sections using stress function approach as well as further development based on hybrid approach.

Figure 2 .
Figure 2. Shear Stress Contour of a Homogeneous Isotropic Square Section

Figure 3 .
Figure 3. Shear Stress Contour of Homogeneous Orthotropic Elliptic Section Shear stress YZ (Cart) (a) linear elements, displacement function Shear stress YZ (Cart) (b) quadratic elements, displacement function

Figure 4 .
Figure 4. Shear Stress Contour of Homogeneous Anisotropic Elliptical Section

Table 1 .
Torsional Rigidity of Homogeneous Isotropic Square Section

Table 3 .
Torsional Rigidity of Homogeneous Anisotropic Elliptical Section

Table 4 .
Torsional Rigidity of Homogeneous Isotropic Multiply-connected Section

Table 7 .
Required Mesh Size to Achieve 1% Accuracy