1 Introduction

In this paper, we study the torsion (twisting rotation) of prismatic bars about their longi- tudinal axis. The torsion problem has an important aspect in the theory of elasticity. It traced back to de Coulomb in 1784, who investigated the torsion of a circular cylinder^{[1, 2]}. He assumed that a plane section, normal to the bar axis, remained a plane during the torsion. Afterwards in 1824, Navier investigated the torsion of a straight bar with a non-circular section and assumed that the plane section still remained a plane during the torsion. This led to self-contradictory results^[1], and it was explained how this assumption contradicted the assumption that the lat- eral surface of the bar under torsion was free from external forces. Finally, de Saint-Venant observed that in the absence of axial symmetry, the plane section might not remain a plane. By introducing a warping function, the governing equation of the torsion of the alternative formulation was derived^{[1, 2]}.

There is extensive literature on the Saint-Venant torsion problems for bars including those with a rectangular cross-section and with other shapes of cross-section^{[3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27]}. Because of the space limitation, here we focus our review on literature related to the bars with a rectangular cross-section and the papers within the rectangular section besides other shapes of section^{[13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27]}.

In the beginning, the existing analytical solutions are reviewed. The torsion of a prismatic bar of rectangular cross-section subjected to an axial torque acting on one of its ends, the other end being fixed was analyzed by Allam and Ghaleb^[13]. The bar was made of a viscoelastic material reinforced by elastic layers parallel to the longer sides of the cross-section. Using a generalization of Ilyushin’s approximation method, the stresses were determined. The problem solution was carried out for the specific case where the applied torque was a step function of the time.

The solution to the Saint-Venant torsion problem of an orthotropic and inhomogeneous straight bar with the rectangular cross-section was investigated by Xu et al.^[14]. The problem was solved for two cases. In the first case of study, material properties of the bar varied exponentially in one direction, and the problem was solved analytically. In the second case, i.e., for the case of arbitrary variation of the material properties, a layerwise model was considered in which the functionally graded material section was modeled by multiple layers with constant material properties in each layer. The warping function, stress components, and torsional rigidity were also analytically expressed in terms of Fourier series for both the cases.

The Saint-Venant torsion problem of non-homogeneous bars with solid and hollow cross- sections was treated by Ecsedi^[15]. The shear modulus of the non-homogeneous bar was a given function which matched the Prandtl’s stress function of the cylindrical bar assuming that it is made of the homogeneous material. The problem was examined for different cases consisting of the hollow elliptical cylinder, solid equilateral triangular cross-section, thin-walled tube, rectangular cross-section, and narrow rectangular cross-section.

Mindlin^[16] solved the Saint-Venant torsion problem of elastic bars by expansion of the displacements in terms of a double power series of the transverse coordinates. The coefficients of the series were determined by solving the linear algebraic equations and deduced from a variational principle. The solution was examined for bars with elliptic, equilateral triangular, and rectangular cross-sections. The number of the coefficients depended on the shape of cross- section (for example, for the elliptic and equilateral triangular sections, it was finite, and for the rectangular cross-section, it was infinite).

Recently, the analysis of torsion of piezoelectric rods was the subject of several studies. The torsion analysis of laminated orthotropic rods with the constant cross-section and coated with thin integrated piezoelectric actuator layers was accomplished by Zehetner^[17]. The rod was subjected to a pair of constant external torques acting in opposite directions at two ends. The distribution of eigenstrains that was able to completely compensate the angle of twist caused by external torsional moments was derived. By introducing an additional warping function, the Saint-Venant’s theory of torsion for laminated orthotropic rods was extended in the presence of eigenstrains. The actuating torsional moment as a function of the eigenstrains and the additional warping function were calculated. In the case of rectangular cross-section (consisting of a homogenous isotropic layer and two piezoelectric layers), an analytical solution for the actuating moment and the additional warping function were derived.

Using the theory of elasticity, exact analytical and numerical solutions for the problem of piezoelectric rods under the static torsion were found by Maleki and Naei^[18]. The range of valid region for the assumption of linear distribution of electric potential through the cross-section and shape-effects on the induced piezoelectric deformation were studied. The problem was examined for the rectangular cross-section. In another study, Maleki et al.^[19] presented the exact solution for torsion of multilayer piezoelectric bars with the rectangular cross-section. In this study, multilayer piezoelectric materials with different properties were considered whilst the axis of torsion and the polarization axes of each layer made different orientations with respect to each other.

There is a subgroup of torsion problems with emphasis on numerical methods which are reviewed below.

The Saint-Venant torsion problem of a beam with the arbitrary and uniform cross-section was studied by Barone et al.^[20]. Among the numerical methods of analyzing the Laplace equation, some of them discretize the whole domain, such as the boundary difference method (BDM) and the finite element method (FEM). There are other methods using the complex analysis including the complex polynomial method (CPM)^{[21, 22]}, the complex variable boundary element method (CVBEM)^{[23, 24, 25]}, and the line element-less method (LEM)^[26]. In the paper by Barone et al.^[20], a comparison between the above-mentioned three boundary element methods (BEMs) for torsion problems was performed. The CPM is based on the expansion of a complex potential in the Taylor series, in which the unknown coefficients are computed by means of collocation points on the boundary. This method was limited to simply connected domains. The CVBEM is based on Cauchy’s integral formula of the theory of functions of complex variables. The Cauchy integral expresses the value of an analytic function at a point in a domain by the integral of the function along the boundary of the domain. By letting the point approach the boundary, one obtains an integral equation. The CVBEM benefits with Cauchy’s integral formula that returns the solution of Laplace equation when mixed boundary conditions on both real and imaginary parts of the complex potential are known. The CVBEM method is the generalization of Cauchy’s integral formula into the boundary integral equation method. The method consists of discretizing the boundary into segments called the complex boundary elements and replacing the analytic function on the boundary by an interpolation function which can then be integrated analytically. Both the CPM and the CVBEM require the satisfaction of local conditions on the boundary of the analyzed domain. The LEM introduces the expansion of the double-ended Laurent series involving harmonic polynomials, proposing an element-free weak form procedure, by imposing that the square of the net flux of the shear stress across the border is minimized with respect to the series coefficients. The LEM does not require any discretization, neither of the boundary nor of the domain.

The approximate solution of the torsion problem using the LEM was proposed by Di Paola et al.^[26]. In this method, using the theory of analytic complex functions, a new potential function in terms of shear stresses was proposed and expanded in the double-ended Laurent series involving harmonic polynomials. The equilibrium and compatibility equations remain fulfilled in the whole domain as those in the BEM theory. In the BEM, the shape functions are the fundamental functions, but in the LEM, the shape functions are harmonic polynomials that satisfy the field equations in the whole domain. By proposing an element-free weak form procedure and imposing the square of the net flux across the border to be minimum with respect to the coefficients of expansion, the method is implemented. Applying the LEM leads to a system of linear algebraic equations including symmetric and positive-definite matrices. As mentioned before, the LEM does not require any discretization, neither of the boundary nor of the domain, therefore, it is a truly meshless method. For a simply connected cross-section, only the regular part of the Laurent series may be used, while for multiply-connected domains, the principal part gives major contributions to the shear stress distribution. The method is evaluated for a square and a circular-hollow cross-sections.

Using the complex variable BEM, the torsion problem of composite shafts consisting of a cylindrical matrix surrounding a finite number of inclusions was solved by Shams-Ahmadi and Chou^[27]. In this method, the problem was reduced to the solution of a singular integral equation in terms of an analytic function of a complex variable using the Cauchy integral. The ensuing integral equation was then solved numerically by discretizing the boundaries into segments, namely, complex boundary elements and replacing the analytic function on the boundaries by an interpolation function. Numerical examples were given for a square shaft with a circular inclusion, as well as for an elliptic shaft with two elliptic inclusions.

Using the LEM method, approximate solutions of the Saint-Venant torsion problem for or- thotropic beams with simply and multiply connected cross-sections were obtained by Santoro^[28]. A suitable transformation of coordinates and a complex potential function in the form of double- ended Laurent series, which is analytic over the transformed domain, were utilised for the prob- lem solution. An element-free weak form procedure was adopted to minimize the square of the net flux of the shear stress across the border with respect to the series coefficients. Finally, the complete shear stress field was obtained.

The non-uniform torsion of multi-material composite bars with the constant arbitrary shape cross-section was the subject of the study by Sapountzakis^[29]. These different materials were perfectly bonded together, and the bar was subjected to an arbitrarily concentrated or twisting torque. The bar ends conditions were specified in terms of twisting and torsional moments and restricted to general linear torsional boundary conditions, including all types of the conventional boundary conditions (clamped, simply supported, free or guided edge). The Saint-Venant torsional shear stresses and the warping normal stresses were also computed. The two ensuing boundary value problems in terms of the bar angle of twist and the warping functions with respect to shear center were solved utilizing the BEM approach. Both the variable warping and torsion constants were computed by using a Gaussian integration over domains of arbitrary shape.

The torsion of non-uniform composite bars of arbitrary constant cross-section was studied by Sapountzakis and Mokos^[30]. Using the BEM, the composite bar consisted of a matrix engulfed by a finite number of inclusions with different Young’s and shear moduli. The matrix and inclusions were perfectly bonded together. An arbitrarily concentrated or distributed twisting moment was applied to the bar where the bar ends were restricted to general linear torsional boundary conditions. The Saint-Venant torsional shear stresses, the warping normal and shear stresses were computed. The problem was investigated for some problems including a cantilever beam of a composite material with three types as follows: the elliptic cross-section with two symmetrically placed elliptic inclusions, the square cross-section with a circular inclusion, and the I-beam totally encased in the rectangular concrete.

In the last part of this review, we cite investigations dealing with other concepts related to the torsion of bars such as plastic torsion, elastic-plastic torsion, and transient creep bounds.

A rectangular bar composed of two dissimilar materials with different yield limits and sub- jected to the torsion was analysed by Rychlewski^[31]. The bar consisted of two dissimilar rect- angular cross-sections attached to each other on one side. The plastic jump non-homogeneity including change of the yield limit was investigated. The yield condition was represented by equating the sum of squares of the nonzero shear stress components to a constant. The problem was solved, and the boundary of the plastic region or yield limit was derived, using a stress function and applying the traction free condition on the outer boundary of the cross-section and the continuity condition of the normal shear stress on the joint side of the two dissimilar rectangles. The problem of determining the stresses satisfying the equilibrium equation and the yield condition was therefore equivalent to the construction of a continuous stress function with the constant slope, but different in the two dissimilar rectangular cross-sections. Depending on the dimensions of the dissimilar rectangles and the ratio of yield limits of those, the problem of limit load of the bar had nine different solutions.

Using the cellular analogy, the Saint-Venant torsion problem of multi-celllular solid sections was solved by Johnson^[32]. In this method, the rectangular portion of the solid section was equalized with the thin walled rectangular cells. The entire section was sub-divided into several cells, and constant shear flows were presumed to exist in each cell. The shear flows in the cells were determined by the enforcement of geometric compatibility on the full set of cells. The method thus provided a tool for the determination of shear flows and stresses in the case of elastic analysis. The method was also convenient for elastic-plastic analysis. For a elastic-pure plastic material, the shear flow became constant when the shear stress in a cell has reached its yield value. A rectangular section and an I-section were exemplified in this context.

The transient primary and secondary creep analysis in thin rectangular strips subjected to the torsion was conducted by Edelstein^[33]. The total strain was considered to be sum of elastic strain and creep strain, whereby the elastic strain was related to stress components via constitutive relations. Also, the time rate of creep strain was stated as a function of the deviatoric stress components, the effective stress, and the effective creep strain. Assuming thin strip and using the Saint-Venant torsion theory, the displacement components were introduced in terms of the time dependent on the specific angle of twist and the time independent on the warping function. According to the above relation and applying the boundary condition, the basic creep equations were found. Individual integral equations were derived for both the primary and secondary creeps. These equations were in terms of applied moment and the only non-zero shear stress component. Finally, the limiting creep rate was obtained. The similarity between these equations and others arising in the previous work was utilized to obtain bounds, monotonicity and convexity of the stress profile and uniform approximations.

Based on the review above, it is found that the analytic solution of the Saint-Venant torsion problem of coated bars with the rectangular cross-section and with effective coating layers has not been the subject of previous studies. An effective coating layer here means that a coating with given material properties and thicknesses will create a maximum von Mises stress on the cross-section lower than the allowable value so that coating material constitutive behaviour would be in the linear elastic range. Moreover, such a solution can be useful for many important applications. In this article, the Saint-Venant torsion theory is used to derive the governing equation of the problem in terms of a warping function. It is then solved by means of the finite Fourier cosine transform (see Subsection 2.1). Next, the torsional rigidity of the coated bar is calculated (see Subsection 2.2). The corresponding results (see Section 2) are validated against existing solutions of an uncoated isotropic bar with the rectangular cross-section from literature. Numerical results are presented in Section 3, and these are compared to those of an uncoated isotropic bar. Concluding remarks are presented in Section 4.

2 Problem formulation 2.1 Solution of governing equation

We consider a dissimilar bar of rectangular cross-section with a finite length a and a width h in x- and y-directions, respectively, as shown in Fig. 1. The x-direction is located on the lower edge of the rectangle. Similarly, the y-direction is situated on the left edge of it. The cross- section consists of three isotropic sub-rectangular regions with finite heights h₁, h − (h₁ + h₂), and h₂, which are attached together along the horizontal direction (x-direction). The lower and upper rectangular planes play the role of coating layers. By the Saint-Venant torsion theory, the displacement components are given as follows^[34]:

The out of plane displacement component w(x, y) under the Saint-Venant torsion is of interest. Consequently, the two nonzero stress components in the z-direction are

where G_l (l = 1, 2, 3) are the orthotropic shear moduli of the lower, middle, and upper layers, respectively. The equilibrium equation in the absence of body force σ_{ij, j} = 0, in view of Eq. (1), reduces to

Applying the finite Fourier cosine transform to the above equation leads to

where   = π/a, and the indices 1, 2, and 3 are used to refer to the lower, middle, and upper rectangular regions, respectively. The traction-free condition on the lateral surface of the bar requires that By virtue of Eqs. (4) and (5), we arrive at Using the inverse finite Fourier cosine transform where in light of Eq. (5), we have Φ_l(0, y) = A_l0 + B_l0y (l = 1, 2, 3). The following expansions can be easily written as

Using the above expansions and Eqs. (2) and (7), the stress components are written as follows:

where . The continuity conditions in the interface of rectangular planes are as follows:

Using Eq. (1), the conditions (5) and (10) are applied to Eqs. (7) and (9), leading to

where α = G₁/G₂, and α = G₃/G₂. By solving the above equations, we arrive at where , in which H = h − (h₁ + h₂), and

Substituting the above coefficients into the stress components (9) leads to

where δ = 2 /a θ. For h₁ = h₂, we have . Therefore, in view of the relations (13), it is easy to show that for h₁ = h₂, we have

For h₁ = h₂, the stress fields of Eq. (14) are symmetric with respect to the parameters α and α. To show this symmetry, we consider two types of bars. In the first type, the shear moduli of the lower and the upper coating layers of the bar are G₁ and G₃, respectively. In the second type, the shear moduli of the lower and the upper coating layers of the bar are G₃ and G₁, respectively. According to relations (14) and (15), the stress components in the lower layer σ_xz(x, y) and σ_yz(x, y), 0 ≤ y ≤ h₁ of the first type are equal to the stress components in the upper layer of the second type at the corresponding points, i.e., σ_xz(x, h−y) and σ_yz(x, h − y), h − h₁ = h₁ + H ≤ y ≤ h.

For the special case, i.e., the case that the whole cross-section of the bar is made from a sim- ilar isotropic material, we set α = α = 1. Subsequently, we arrive at , and the stress component (13) reduces to σ_xz(x, y) = δG₂

In this case, for the isotropic bar (G = 1), the stress fields (16) reduce to identical results of Ref. ^[35].

2.2 Calculation of torsional rigidity

Substituting the stress relations (2) into = Dθ, the torsional rigidity of the bar can be calculated by using the following formula:

where M = Dθ is the applied torque, and

Noting Eq. (7) and using 2 we arrive at

Substituting Eq. (19) into Eq. (17), we have

For the special case α = α = 1, the torsional rigidity of the bar reduces to^[35]

3 Numerical results

For the verification purpose, we consider an isotropic bar with the rectangular cross-section and without any coating layer. The shear modulus of the bar is G₂, and its dimensions are a and h. The values of the parameters can be taken from the table on the page 277 of Ref. ^[34]. These values are very close to those obtained by the relationship given here for the special case of α = α = 1, which verifies the presented formulation of the paper.

The initial stage for our study is the state of stress induced by applying the axial torque M. We consider a bar with rectangular cross-section with a = 2h. The thicknesses of the upper and lower coating layers are h₁ = h₂ = 0.2h. The material of the coating layers is aluminium with the shear modulus G₀ = 26 GPa, while the substrate is made of steel with the shear modulus G₀ = 80 GPa. The contour of the dimensionless stress component σ′yz = Dσ_yz/(G₂Mh) obtained from the analytical solutions is compared with that predicted by the FEM (using ANSYS software) in Fig. 2. Continuity of the dimensionless stress component σ′yz at the interfaces between the coating layers and the bar, i.e., y = h₁ and y = h − h₂, can be observed. The results obtained from ANSYS which are for a mid-plane of the rectangular beam are in good agreement with those obtained from this study.

Also, the contour of the dimensionless stress component σ′x z = Dσ_xz/(G₂Mh) is presented in Fig. 3. Discontinuities or jumps in the dimensionless stress component σ′xz at y = h₁ and y = h−h₂ can be noticed. The maximum values of the above dimensionless stress components for a coated bar with h₁ = h₂ = 0.2h are max(σ′xz) = 0.581 5 and max(σ′yz) = 0.539 0, respectively. These values for an isotropic bar without coating (h₁ = h₂ = 0) are max(σ′xz) = 0.930 4 and max(σ′yz) = 0.723 2, respectively. Hence, the existence of the softer coated layers reduces the maximum values of dimensionless stress components.

The main objective is to find a criterion to utilise the coating layers properly by controlling the level of maximum stress in the cross-section. To this end, we use von Mises stress in the cross-section. According to the von Mises yield criterion, an isotropic and ductile metal, which is subjected to a complex loading condition, will be obtained when the von Mises stress σ_v is equal to metal’s yield strength σ_y. The von Mises stress in the Saint-Venant torsion is . Since in torsion problems we deal only with shear stresses, and the shear yield stress τ_y is approximately half of the tensile yield strength σ_y, i.e., τ_y = 0.5σ_y, therefore, here, we use the equivalent shear von Mises stress as 

We choose the shear moduli of the coating layers G₁ and G₃ or α = G₁/G₂ and α = G₃/G₂ such that max(τ_e) ≤ τ_all where all is the allowable shear stress level in the cross-section. We introduce a dimensionless von Mises stress as τ′_e = D_τe/(G₂Mh). The variation of α with α, for different values of the maximum dimensionless von Mises stress, _{τe, max}′ = max |τ_e′ | and also for h₁ = h₂ = 0.1h, can be seen in Fig. 4. By choosing α and finding the corresponding α from Fig. 4, it is possible to control the level of the maximum stress in the cross-section, which can serve as an important criterion for the structural integrity of coated bars. As seen for thin coating layers (as h₁, h₂ 6 0.1h), each segment of graph is almost linear. For example, in the cases of τ_{e, max}′ = 0.9 and τ_{e, max}′ = 0.95, the main part of each curve appears to be the lines with negative slopes as α = 1.055 4−0.097 6α and α = 1.122 3−0.102 4α, respectively. Therefore, we expect a constant value for the maximum dimensionless von Mises stress τ_{e, max}′ provided that α and α satisfy the relevant line equation. We now furnish another graph to demonstrate the aforementioned feature of Fig. 4 by comparing the effects of coating stiffness on the maximum dimensionless von Mises stress of the section. The variation of the maximum dimensionless von Mises stress with α for different line equations is shown in Fig. 5. It is noteworthy that for the two above-mentioned line equations, the maximum dimensionless von Mises stress τ_{e, max}′ does not change with the variation of α. Again, the differences between the analytical solution obtained in this study and the numerical solution obtained using ANSYS are negligible.

Analogous to Fig. 4, the variation of α versus α for different values of max |Dτ_e/(G₂Mh)| and h₁ = h₂ = 0.2h is plotted in Fig. 6. By choosing the stiffness of the lower coating layer, the parameter α is assigned, and according to the thickness of layers and the desired maximum dimensionless von Mises stress, the graphs similar to Fig. 4 and Fig. 6 are used to determine the corresponding parameter α. Also, the graphs similar to Fig. 4 and Fig. 6 are depicted in Fig. 7. Viewing the Fig. 4, Fig. 6, and Fig. 7, we find that by increasing the thickness of the coating layers, the graphs parts are not linear, and the variation of α in terms of α can not be approximated by a line equation. Also, the symmetry of the stress field with respect to parameters α and α can be noticed from Fig. 4, Fig. 6, and Fig. 7.

To evaluate the effect of the thicknesses of the layers, h₁, h₂ and the ratios α, α on the maximum von Mises stress, we consider a specific case with h₁ = h₂ and α = α. The variation of α versus h₁/h for different values of max |Dτ_e/(G₂Mh)| is shown in Fig. 8. As it can be seen, increasing the thicknesses of the layer has a less effect on the required stiffness of the layers to achieve the desired values of max |Dτ_e/(G₂Mh)|. However, to achieve low values of max |Dτ_e/(G₂Mh)|, weaker coating layers, i.e., lower values of α, are needed.

4 Conclusions

In this paper, an analytical solution based on the Saint-Venant torsion theory has been de- veloped to investigate the effects of coating layers on the maximum stress in the cross-section of rectangular bars. The ensuing results show the continuity of the dimensionless stress component σ′_yz and the axial displacement at the interface of the bar and its coating layer.

The results are summarized as follows:

A coating layer, which is softer (less stiff) compared to the substrate bar, reduces the maximum values of dimensionless stress components in the cross-section.

Parameters α = G₁/G₂ and α = G₃/G₂ representing the ratio of the material shear modulus of each coating layer with respect to that of substrate bar, can be suitably selected to vary the maximum stress in the bar cross-section.

For thin coating layers (h₁, h₂ 6 0.1h), there is a specific range of variation of parameter α, within which α can be varied linearly with α that will keep the maximum von Mises stress in the cross-section almost constant.