Torsional shear stress with arbitrary cross-sections in homogeneous isotropic elastic material using finite element method

Article no. 30 THE CIVIL ENGINEERING JOURNAL 2-2021 DOI 10.14311/CEJ.2021.02.0030 409 TORSIONAL SHEAR STRESS WITH ARBITRARY CROSSSECTIONS IN HOMOGENEOUS ISOTROPIC ELASTIC MATERIAL USING FINITE ELEMENT METHOD Dang-Bao Tran1,2 1. VSB–Technical University of Ostrava, Faculty of Civil Engineering, Department of Structures, Ludvíka Podéště 1875/17, 708 00 Ostrava, Czech Republic; [email protected] 2. Thu Dau Mot University, Faculty of Architecture, Department of Civil Engineering, Tran Van On 06, 75000 Binh Duong Province, Vietnam;[email protected] ABSTRACT Determining the shear stress of a structural element caused by torsion is a vital problem. The analytical solution of the Saint-Venant torsion is only suitable for simple cross-sections. The numerical methods to evaluate the shear stress due to torsion of complicated cross-sections is indispensable. Many scientists have studied the torsion problem with various numerical methods. This paper aims to present an efficient finite element method for assessing the shear stress with arbitrary cross-sections in homogeneous isotropic elastic material due to torsion. MATLAB is the language for programming the numerical method. The validation examples were performed to show the reliability and efficiency of the author’s numerical method. KEYWORDS Torsional shear stress, Isoparametric eight-noded quadrilateral elements, Saint-Venant torsion, Warping function INTRODUCTION Torsion in the structure occurs due to asymmetrical loads, either by geometric dimensions or by interconnections between members. In many cases, torsion can be governing design factors. It is, therefore, essential to accurately determine the shear stress caused by torsion. Saint-Venant analysed the torsion problem using the semi-inverse method, assuming an unknown displacement to satisfy the equilibrium equations and boundary conditions. Prandtl then introduced the stress function of the Saint-Venant torsion and the method of membrane analogy [1]. However, determining the shear stress due to torsion by the analytic solution is complicated for a complex single connected domain or multiply connected domain. A single connected domain is a domain where the crosssection is bounded by a closed. The multiply connected domain is the domain where the crosssection is bordered by several closed. Nowadays, structural elements such as beams and columns are more and more complicated in their shapes. Hence the use of numerical methods to determine shear stress due to torsion is indispensable. Various numerical methods to assess torsional shear stress have been performed by many researchers [2-22]. Ely, J. F., and Zienkiewicz, O. C. [2] first solved Poisson’s equation of Prandtl’s stress function using the finite difference method and investigated the rectangular section with and without holes. Herrmann, L. R. [3] utilized the finite element method to calculate the warping function of the torsion of irregular sectional shapes. Based on the Hellinger–Reissner principle, Xiao, Q. Z., et al. [4] developed a 4-node element with four stress parameters to determine the shear Article no. 30 THE CIVIL ENGINEERING JOURNAL 2-2021 DOI 10.14311/CEJ.2021.02.0030 410 stress for the polygonal section. Gruttmann, F. et al. [5] used the finite element method to evaluate shear stress using the warping function, which is more convenient than using Prandtl’s stress function when considering a multiply connected domain. Gruttmann’s numerical method has been implemented into an enhanced version of the program FEAP [6], which used isoparametric fournoded quadrilateral elements. Fialko, S.Y. et al. [7] developed a numerical method using constant triangular elements to solve the Saint Venant problem of torsion, and torsionless bending of prismatic bars is realized in the SCAD software [8]. Jog, S. et al. [9] presented a finite-element formulation for Saint Venant torsion with a multiply-connected domain in an anisotropic material. Recently, Beheshti, A. [10] derived a finite element from strain-gradient elasticity for torsion of prismatic bars with minimal dimensions. Besides the finite difference and finite element method, many authors have applied other numerical methods such as the boundary element method (BEM) [11-16], line element-less method (LEM) [17-18], a null-field integral technique [19], and finite-volume method [20-21] to analyse the torsion problem. Katsikadelis, J. T. et al. [11] examined the torsion of general composite bars by BEM while Gaspari, D. et al. [12] tackled orthotropic beams with a polygonal cross-section. Barone, G. et al. [13] implemented the Complex Variable Boundary Element Method to examine the torsion problem in the single connected domain. Lee, J.W. et al. [14] obtained a new BEM from the general Cauchy integral formula derived from the Borel–Pompeiu formula to analyse the torsion problem. Paradiso, M. et al. [15] achieved an efficient method to determine the warping function parameter with the general cross-section. Chen, K. H. et al. [16] introduced a new error estimation technique in BEM to optimise the torsion problem with a multiply-connected domain. Di Paola, M. et al. [17] proposed LEM to deal with shear stress in torsion problem with isotropic material and arbitrary crosssection. Santoro, R. [18] handled the Saint Venant torsion problem for orthotropic beams with a general cross-section by LEM. Chen, J-T. et al. [19] introduced the null-field integral technique to analyse the torsion problem of circular cross-sections with round holes. Chen, H. et al. [20-21] developed a new finite-volume based on Bansal and Pindera’s work to investigate Saint Venant’s torsion problems of homogeneous and composite prismatic bars with multiply connected domain. In summary, all of the works show the feasibility and effectiveness in academia. Currently, to the author’s knowledge, Gruttmann’s method [5] has been developed into the FEAP program [6 ] of the University of California, Berkeley, and Allplan Bridge [22]. Fialko’s method [7] is similar to Gruttmann, developed as a module in the SCAD commercial software [8]. It means the methods of Gruttmann and Fialko are practical. However, due to the use of the four-noded quadrilateral element in FEAP and the constant strain triangle element in SCAD, to achieve high accuracy, it is necessary to mesh very smoothly, which affects the calculation speed. So, this research aims to develop a new numerical method (NMB) based on the work of Gruttmann by using the isoparametric eight-noded quadrilateral element. The validation examples were performed to show the reliability and efficiency of NMB. FINITE ELEMENT METHOD PROCEDURE Figure 1 shows the arbitrary cross-section of prismatic beam in homogeneous isotropic elastic material, the longitudinal axis is the x-axis, the cross-section denoted  is in yz plane. The multiply connected domain  is bounded by 1 2 1 , ,..., , nn     . S is the centre of gravity. On 1 2 1 , ,..., , nn     the right-handed orthogonal basis system is defined with tangent vector t and outward normal vector [ , ]T yz nnn . With t the orientation of the associated coordinate s is uniquely defined. The displacement field x y z [u ,u ,u ]T u is expressed as [1] T x u   , yx uz   , , zx uy   (1) Article no. 30 THE CIVIL ENGINEERING JOURNAL 2-2021 DOI 10.14311/CEJ.2021.02.0030 411 where x  : torsion angle, x d dx    , ( , ) Tyz  : warping function for torsion. Here, the constraint of the warping for torsion is 0. TdA     (2) Fig. 1 - Cross-section of a prismatic bar. The shear stresses are given as T xy Gz y         , . T xz Gy z         (3) The polar second moment of area can be written as . TT T I y y z z dA zy                         (4) The governing differential equation for the Saint-Venant torsion is expressed as 22 22 0. TT in yz       (5) Meanwhile, the boundary condition is given as ( 1,2,..., ), TT y z y z i n n n z n y on i n yz         (6) where y dz nds  , . z dy nds  (7) The governing differential equation (5) is transformed into the weak form by using the Galerkin’s method as below ( , ) ( ) 0, TT Tyz G dA n z n y ds y y z z                           (8) with test function 1()H   . Article no. 30 THE CIVIL ENGINEERING JOURNAL 2-2021 DOI 10.14311/CEJ.2021.02.0030 412 The equation (8) is approximated by using the finite element method. Gruttmann’s method used isoparametric four-noded quadrilateral elements. To improve the computation speed, NMB used isoparametric eight-noded quadrilateral elements. The [ , ]T yzx , the unknown function T  and the test function  are interpolated in the local coordinate system as follows   8 1 , hii i N    xx ,   8 1 ( ) , T h T ii i N       ,   8 1 ( ) , , hii i N       (9) where h denotes the approximate solution of the finite element method,   , i N  denotes the shape function of the element. Figure 2 shows the isoparametric eight-noded quadrilateral element used in NMB. The shape functions of this element can be described as follows [23, 24] 12 34 22 56 2 78 11 ( , ) (1 )(1 )( 1 ), ( , ) (1 )(1 )( 1 ), 44 11 ( , ) (1 )(1 )( 1 ), ( , ) (1 )(1 )( 1 ), 44 11 ( , ) (1 )(1 ), ( , ) (1 )(1 ), 22 11 ( , ) (1 )(1 ), ( , ) (1 )(1 22 NN NN NN NN x h x h x h x h x h x h x h x h x h x h x h x h x h x h x h x h x h x h x h x = - - - - - = + - - + - = + - - + + = - + - + + = - - = + - = - + = - - 2).h (10) Fig 2 -. Isoparametric eight-noded quadrilateral element. Inserting the derivatives of () Th  and h  into the equation (8) reads as 88 11 1 ( ) 0. numel e T e ij j i ij e KF      (11) Here, denotes the assembly operator with numel the total number finite elements. The stiffness part e ij K to the nodes i and j and the right hand e i F read 11 11 11 , ( , ) ( , ) ( , ) ( , ) ( , ) . e j j j j ei i i i ij e Q Pjj eii ij p q p q p q p q p q p q pq N N N N N N N N K dA d d y y z z y y z z NN NN K w w y y z z                                                              J J (12) Article no. 30 THE CIVIL ENGINEERING JOURNAL 2-2021 DOI 10.14311/CEJ.2021.02.0030 413 11 11 11 , ( , ) ( , ) ( , ) . e ei i i i ie Q P eii i p q p q p q p q pq N N N N F z y dA z y d d y z y z NN F w w z y yz                                              J J (13) where J denoted as Jacobian matrix defined as , yz yz Jxx hh éů ¶¶ ęú ęú ¶¶ ęú =ęú ¶¶ ęú ęú ¶¶ ëű (14) p w and q w are the weights and p  and q  are the integration points of the Gaussian integration technique. NMB uses 3 x 3 Gauss quadrature derived from the 1D case, where the quadrature points are located at 3/5 , 0 and 3/ 5 , and the corresponding weights are equal to 5/9, 8/9, and 5/9, respectively (see [23, 24]). The value T i  of one arbitrary nodal point i has to be value 0. VALIDATION EXAMPLES The objective of this section is to demonstrate the reliability and effectiveness of NMB. For this purpose, four validation examples derived from [1, 5, 9] were examined, and their results are compared with those analyzed by NMB implemented in MATLAB R2015a. A bar of square cross-section subjected to the torsion moment 1 T M kN.m with the length of the edge 2 m is analysed in the first example. Jog, C.S. et al. [9] investigated this problem using 16 isoparametric nine-noded quadrilateral elements. To achieve convergence values of the shear stress and the polar second moment of area, the square cross-section is divided into 16, 64, 256 elements with uniform mesh by NMB, respectively. Figure 3 shows the discretization of the square cross-section with 16 eight-noded quadrilateral elements. The points, A, B, C, D, correspond with the coordinates (0, 1), (2, 1), (0, 0), (0, 2), respectively. The analytical results of maximum shear stress and polar second moment of area obtained from [1, 9] are max 0.592   kPa, 2.24923 T I m4, respectively. Table 1 shows the comparison of the results of the maximum shear stress, the polar second moment of area between the methods. The results of NMB-16 elements and Jog, C.S. et al. [9] are the same. When the square cross-section is refined into 256 elements, NMB is in good agreement with analytical solutions. Tab. 1 - The maximum shear stress max  and the polar second moment of area T I of square crosssection Factors Analytical solutions Jog, C.S. et al. [9] NMB16 elements NMB256 elements () max a  [kPa] 0.592 0.6173 0.619174 0.601461 ()b T I [m4] 2.24923 2.2519 2.25187 2.24925 Error(a), (%) - 4.273 4.590 1.598 Error(b), (%) - 0.1187 0.1173 0.000 Figure 4 depicts the distribution of shear stress xz  for square cross-section, which decreases the magnitude gradually from the midpoint of the edge to the centre of gravity along the y-axis [1]. Article no. 30 THE CIVIL ENGINEERING JOURNAL 2-2021 DOI 10.14311/CEJ.2021.02.0030 414 Figure 5 and 6 present the variation of the shear stress xz  along the AB and CD segment, respectively. With 16 elements, NMB and Jog, C.S. et al. [9] cannot capture precisely the shear stress distribution along the CD segment. When the mesh is fine enough, in this case, 256 elements, the results of NMB and theory harmonize very well. Fig. 3 - Discretization of the square crosssection with 16 elements. Fig. 4 - Distributed shear stresses xz  for square cross-section. Fig. 5 - The variation of shear stress xz  along the AB segment. Fig. 6 - The variation of shear stress xz  along the CD segment. The second example is an equilateral triangle subjected to the torsion moment 1 T M kN.m with the height of the triangle 0.2 m. The analytical results of maximum shear stress max  and polar second moment of area T I are 1.62380 MPa and 6158.40 cm4, respectively [1]. The triangular crosssection was discretized 4, 8, 14, 37, 57, 109, 658 eight-noded quadrilateral elements with nonuniform mesh to obtain the convergence results. Article no. 30 THE CIVIL ENGINEERING JOURNAL 2-2021 DOI 10.14311/CEJ.2021.02.0030 415 Tab. 2 - The maximum shear stress max  and the polar second moment of area T I of triangular crosssection Factors Analytical solution NMB658 elements Error, (%) max  [MPa] 1.62380 1.6247 0.055 T I [cm4] 6158.40 6158.40 0.000 It is clear from Table 2 that the results of the maximum shear stress, the polar second moment of area obtained from NMB are in good agreement with the theoretical solution [1]. Figure 7 shows the triangular cross-section divided by 14 eight-noded quadrilateral elements. The points, E, F, correspond with the coordinates (0, 0.2 3 ), (0.2, 0.2 3 ), respectively. Figure 8 depicts the stress distribution xy  for triangular cross-section, where magnitude values do not exist along with segment EF [1]. Figure 9 presents the stress distribution xz  for triangular cross-section. It can observe from Figure 10 that the result of the variation of shear stress xz  along the EF segments of NMB and analytical method is very well matched. Fig. 7 - Discretization of the triangular crosssection with 14 elements. Fig. 8 - Distributed shear stresses xy  for triangular cross-section. Fig. 9 - Distributed shear stresses xz  for triangular cross-section. Fig. 10 - The variation of shear stress xz  along the EF segment with NMB-109 element. Article no. 30 THE CIVIL ENGINEERING JOURNAL 2-2021 DOI 10.14311/CEJ.2021.02.0030 416 The third example concerns a bar of square cross-section with hole applied to the torsion moment 1 T M kN.m. Figure 11 shows the geometrical dimensions and the meshing to 96 isoparametric eight-noded quadrilateral elements by NMB of the square cross-section with hole. The points, G, H, I, J, correspond with the coordinates (0, 1), (1.6, 1), (0, 0), (0, 2), respectively. Jog, C.S. et al. [9] analysed this problem using 96, 1536, 6144 isoparametric nine-noded quadrilateral elements, respectively. Figure 12 depicts the variation of polar second moment of area corresponding to 96, 1536, 6144 discretization elements investigated by NMB. In Table 3, the convergence of the polar second moment of area value by NMB is presented as a comparison with that obtained from Jog, C.S. et al. [9]. Figure 13 and 14 present the stress distribution xy  , xz  for square cross-section with hole, respectively. Figure 15 and 16 plot the variation of shear stress xz  along the segments, GH, IJ, respectively, evaluated with the two methods, Jog, C.S. et al. [9], and NMB. From the figures and table mentioned above, the reliability of the NMB is once more verified. Fig. 11 - The dimension and the discretization to 96 elements of the square cross-section with hole. Fig. 12 - The variation of polar second moment of area corresponding to discretization elements. Tab. 3 - The polar second moment of area of square cross-section with hole Factors Jog, C.S. et al.- 6144 elements [9] NMB6144 elements Error, (%) T I [m4] 1.707 1.70718 0.01054 Article no. 30 THE CIVIL ENGINEERING JOURNAL 2-2021 DOI 10.14311/CEJ.2021.02.0030 417 Fig. 13 - Distributed shear stresses xy  for the square cross-section with hole. Fig. 14 - Distributed shear stresses xz  for the square cross-section with hole. Fig. 15 - The variation of shear stress xz  along the GH segment. Fig. 16 - The variation of shear stress xz  along the IJ segment. According to Figure 17, the bridge cross-section is an example for multiply connected domains was considered the fourth example. With this example, the polar second moment of area of the bridge cross-sections neglected the cantilever part obtained by the analytical solution is 40.0 m4 [5]. A comparison with a theoretical explanation is not possible. The comparison between NMB and FEAP [5,6], implemented Gruttmann’s numerical method by using isoparametric four-noded quadrilateral element, was performed on the convergence speed. The bridge cross-section was divided into 50000 nodes with a uniform mesh by FEAP to obtain the convergence of the results. In NMB, the bridge cross-section is meshed with 80, 406, 1892, 4634, 11414 nodes to investigate the convergence result. Figure 18 shows the bridge cross-section divided into 80 nodes by NMB. Figure 19 depicts the convergence of the polar second moment of area obtained by NMB with 4634 nodes.