Saint-Venant Torsion of Beams with Rectangular Cross-Section Using the Fourier Analysis

1. Introduction

The Saint-Venant torsion is characterized by the absence of axial stresses due to torsion and only shear stresses are developed. The correct solution to the problem of torsion of prismatic beams was given by Saint-Venant [1]; he made assumptions on the deformation of the twisted bar that could satisfy the equations of equilibrium and the boundary conditions. Prandtl [2] introduced a stress function in terms of which the shear stresses were defined. Therefore, the solution to the torsion problem consists of finding a stress function which satisfies the governing equation and the boundary conditions. Numerical methods were considered by numerous authors. Ritz [3] and Trefftz [4] used energy methods whereby the stress function was approximately determined from the minimum condition of the strain energy of the twisted beam. Tran [5] used isoparametric eight-noded quadrilateral elements in order to improve Gruttmann’s isoparametric four-noded quadrilateral elements. Fogang [6] used finite difference method for the torsional analysis of beams with arbitrary cross-sections; a grid network was considered for the cross-section and additional nodes were introduced at the boundaries. The additional nodes permitted to apply the governing equations at the boundary nodes and to satisfy the boundary conditions. In a thesis Chen [7] developed a new approach to Saint-Venant’s torsion problems based on the finite-volume method. The approach employed the displacement formulation expressed in terms of the warping function subject to Neumann-type boundary condition. Furthermore, the torsional problem was solved for rectangular cross-sections using infinite trigonometric series [8,9]. Ike et al. [10] used the double finite sine transform method to obtain closed-form mathematical solutions for the Saint-Venant torsional analysis of beams with rectangular cross-section. Ike [11] applied the single Fourier sine transform method to the stress compatibility equation, converting the partial differential equation to an easier to solve ordinary differential equation in the transformed domain. The boundary conditions were used to find the integration constants and inversion was used to find the solution in the physical domain.

In this paper, Fourier sine series was utilized to describe the constant value of the Laplacian of the Prandtl stress function, while double sine series or simple sine series was utilized to describe the Prandtl stress function. Closed- form solutions in form of infinite series were derived for shear stresses, torsional constant and warping function.

2. Materials and Methods

2.1. Elasticity Equations of Solid Beams Subjected to Saint-Venant Torsion

In this section, equations of the linear theory of elasticity of beams subjected to Saint-Venant torsion are recalled. The axis convention in the cross-section is represented in Figure 1, whereas the x-axis is the longitudinal beam’s axis. The dimensions of the beam’s cross section in y- and z-directions are denoted by a and b, respectively.

The displacements in x-, y-, and z-direction are denoted by u(x,y,z), v(x,y,z), and w(x,y,z), respectively, and the rotation of cross-section or angle of twist (positive anticlockwise) is denoted by θ(x).

The cross-section is assumed to rotate about an axis through the center of torsion T (y_T, z_T).

$$v=-\left(z-z_T\right)\theta ,\begin{array}{cc}& w=\left(y-y_T\right)\theta \end{array}$$

(1a, b)

The displacements v(x,y,z) and w(x,y,z) are related to the angle of twist θ(x) as follows

The shearing strains γ_xy and γ_xz and the shearing stresses τ_xy and τ_xz are related to the displacements as follows

$${\gamma }_{xy}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x},\begin{array}{cc}& {\gamma }_{xz}=\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\end{array}$$

(2a, b)

$${\tau }_{xy}=G\left(\frac{\partial u}{\partial y}-\left(z-z_T\right)\frac{d\theta }{dx}\right),\begin{array}{cc}& {\tau }_{xz}=G\left(\frac{\partial u}{\partial z}+\left(y-y_T\right)\frac{d\theta }{dx}\right)\end{array}$$

(3a, b)

where G is the shear modulus.

The Prandtl [2] stress function (PSF) ϕ(x,y,z) is defined in terms of stress components

$${\tau }_{xy}=\frac{\partial \varphi }{\partial z},\begin{array}{cc}& {\tau }_{xz}=-\frac{\partial \varphi }{\partial y}\end{array}$$

(4a, b)

The PSF ϕ(x,y,z) must satisfy the following differential equation (a compatibility equation)

$$\frac{{\partial }^2\varphi }{\partial y^2}+\frac{{\partial }^2\varphi }{\partial z^2}=-2G\frac{d\theta }{dx}$$

(5)

It is shown that the PSF has a constant value throughout the boundary: that value is taken zero for the rectangular solid section of the present study. The torque Mt is given by

$$\begin{array}{l}{M}_t={\displaystyle \iint \left[{\tau }_{xz}\left(y-y_T\right)-{\tau }_{xy}\left(z-z_T\right)\right]dydz}\\ \begin{array}{cc}& =2{\displaystyle \iint \varphi dydz}\end{array}\end{array}$$

(6a, b)

The torsional stiffness GIt is related to the torque Mt and the angle of twist θ(x) as follows

$$M_t=G{I}_t\frac{d\theta }{dx}$$

(7)

Combining Equations (6b) and (7) yields the torsion constant It

$$I_t=\frac{2}{G\frac{d\theta }{dx}}{\displaystyle \iint \varphi dydz}$$

(8)

The analysis is then reduced of determining the Prandtl stress function: the latter is determined using double sine series or single sine series.

2.2. Analysis Using Double Sine Series for the Prandtl Stress Function

As a recall, the dimensions of the beam’s cross section in y- and z-directions are denoted by a and b, respectively, as represented in Figure 1. Noting that the sine series hypothesis satisfies the condition that the Prandtl stress function (PSF) is zero along the boundaries of the cross-section, the PSF is described using a double sine series

$$\varphi (x,y,z)=G\frac{d\theta }{dx}{\displaystyle \sum _m{\displaystyle \sum _n{\varphi }_{mn}\sin \frac{m\pi y}{a}}}\sin \frac{n\pi z}{b}$$

(9)

The double Fourier sine series expansion of unity is given by

$$1={\displaystyle \sum _{m=1,3,5\dots }{\displaystyle \sum _{n=1,3,5\dots }\frac{16}{\pi ²mn}\sin \frac{m\pi y}{a}}}\sin \frac{n\pi z}{b}$$

(10a)

Using Equation (10a), the right-hand side of Equation (5) can be expanded into a double Fourier sine series as follows

$$-2G\frac{d\theta }{dx}=-2G\frac{d\theta }{dx}{\displaystyle \sum _{m=1,3,5\dots }{\displaystyle \sum _{n=1,3,5\dots }\frac{16}{\pi ²mn}\sin \frac{m\pi y}{a}}}\sin \frac{n\pi z}{b}$$

(10b)

Substituting Equations (9) and (10b) into (5) yields

$${\displaystyle \sum _{m=1,3,5\dots }{\displaystyle \sum _{n=1,3,5\dots }G\frac{d\theta }{dx}\left[\left(-{\left(\frac{m\pi }{a}\right)}^2-{\left(\frac{n\pi }{b}\right)}^2\right){\varphi }_{mn}+\frac{32}{\pi ²mn}\right]\sin \frac{m\pi y}{a}}}\sin \frac{n\pi z}{b}=0$$

(10c)

Since Equation (10c) must be satisfied for any value of y and z it yields

$$\begin{array}{l}\left(-{\left(\frac{m\pi }{a}\right)}^2-{\left(\frac{n\pi }{b}\right)}^2\right){\varphi }_{mn}+\frac{32}{\pi ²mn}=0\begin{array}{cc}& \end{array}\to {\varphi }_{mn}=\frac{32}{\pi ²mn\left[{\left(\frac{m\pi }{a}\right)}^2+{\left(\frac{n\pi }{b}\right)}^2\right]}\begin{array}{cc}& \to \end{array}\\ \varphi (x,y,z)=G\frac{d\theta }{dx}\times \frac{32}{{\pi }^4}{\displaystyle \sum _{m=1,3,5\dots }{\displaystyle \sum _{n=1,3,5\dots }\frac{1}{mn\left[{\left(\frac{m}{a}\right)}^2+{\left(\frac{n}{b}\right)}^2\right]}\sin \frac{m\pi y}{a}}}\sin \frac{n\pi z}{b}\end{array}$$

(11a, b)

Substituting Equation (11b) into (8) yields the torsional constant

$$\begin{array}{l}{I}_t=\frac{64}{{\pi }^4}{\displaystyle \underset{0}{\overset{b}{\int }}{\displaystyle \underset{0}{\overset{a}{\int }}{\displaystyle \sum _{m=1,3,5\dots }{\displaystyle \sum _{n=1,3,5\dots }\frac{1}{mn\left[{\left(\frac{m}{a}\right)}^2+{\left(\frac{n}{b}\right)}^2\right]}\sin \frac{m\pi y}{a}}}\sin \frac{n\pi z}{b}dydz}}\\ \begin{array}{cc}& =\end{array}\frac{256}{{\pi }^6}ab{\displaystyle \sum _{m=1,3,5\dots }{\displaystyle \sum _{n=1,3,5\dots }\frac{1}{m²n²\left[{\left(\frac{m}{a}\right)}^2+{\left(\frac{n}{b}\right)}^2\right]}}}\end{array}$$

(12a, b)

The torsional constant can be expressed as a function of the ratio a/b as follows

$$\begin{array}{l}{I}_t=k\left(a/b\right)\times b{a}^3\\ \to k\left(a/b\right)=\frac{256}{{\pi }^6}{\displaystyle \sum _{m=1,3,5\dots }{\displaystyle \sum _{n=1,3,5\dots }\frac{1}{m²n²\left[m^2+{\left(\frac{na}{b}\right)}^2\right]}}}\end{array}$$

(13a, b)

The shear stresses are then determined using Equations (4a, b), (7), and (11b)

$$\begin{array}{l}{\tau }_{xy}\left(y,z\right)=\frac{32}{{\pi }^3}\frac{a^2{M}_t}{b{I}_t}{\displaystyle \sum _{m=1,3,5\dots }{\displaystyle \sum _{n=1,3,5\dots }\frac{1}{m\left[m^2+{\left(\frac{na}{b}\right)}^2\right]}\sin \frac{m\pi y}{a}}}\cos \frac{n\pi z}{b}\\ {\tau }_{xz}\left(y,z\right)=-\frac{32}{{\pi }^3}\frac{a{M}_t}{I_t}{\displaystyle \sum _{m=1,3,5\dots }{\displaystyle \sum _{n=1,3,5\dots }\frac{1}{n\left[m^2+{\left(\frac{na}{b}\right)}^2\right]}\cos \frac{m\pi y}{a}}}\sin \frac{n\pi z}{b}\end{array}$$

(14a, b)

Observing that z_T = b/2 and that the warping is zero along the symmetry axis y = a/2, The warping u(x,y,z) is expressed using Equations (3a, b), (7), and (14a) as follows

$$\begin{array}{l}\frac{\partial u}{\partial y}=\frac{{\tau }_{xy}}{G}+\left(z-z_T\right)\frac{M_t}{G{I}_t}\to \\ u=\frac{M_t}{G{I}_t}\left(z-\frac{b}{2}\right)\left(y-\frac{a}{2}\right)-\frac{32}{{\pi }^4}\frac{a^3{M}_t}{bG{I}_t}{\displaystyle \sum _{m=1,3,5\dots }{\displaystyle \sum _{n=1,3,5\dots }\frac{1}{m^2\left[m^2+{\left(\frac{na}{b}\right)}^2\right]}\cos \frac{m\pi y}{a}}}\cos \frac{n\pi z}{b}\end{array}$$

(15a, b)

2.3. Analysis Using Single Sine Series for the Prandtl Stress Function

For simplification purpose the y-axis was shifted to the middle of the plate as shown in Figure 2.

The Prandtl stress function can be described using a single sine series, satisfying so the zero value boundary condition at y = 0 and y = a.

$$\varphi (x,y,z)=G\frac{d\theta }{dx}{\displaystyle \sum _m{\varphi }_m(z)\sin \frac{m\pi y}{a}}$$

(16)

Observing that the single Fourier sine series expansion of unity is

$$1={\displaystyle \sum _{m=1,3,5\dots }\frac{4}{\pi m}\sin \frac{m\pi y}{a}}$$

(17)

the right-hand side of Equation (5) can be expanded into a single Fourier sine series as follows

$$-2G\frac{d\theta }{dx}=-2G\frac{d\theta }{dx}{\displaystyle \sum _{m=1,3,5\dots }\frac{4}{\pi m}\sin \frac{m\pi y}{a}}$$

(18)

Substituting Equations (16) and (18) into (5) yields

$${\displaystyle \sum _{m=1,3,5\dots }G\frac{d\theta }{dx}\left[\frac{d^2{\varphi }_m(z)}{d{z}^2}-{\left(\frac{m\pi }{a}\right)}^2{\varphi }_m(z)+\frac{8}{\pi m}\right]\sin \frac{m\pi y}{a}}=0$$

(18a)

Since Equation (18a) holds for any value of y it yields the following differential equation with ϕ_m(z)

$$\frac{d^2{\varphi }_m(z)}{d{z}^2}-{\left(\frac{m\pi }{a}\right)}^2{\varphi }_m(z)+\frac{8}{\pi m}=0$$

(19)

The solution to Equation (19) that satisfies the zero value boundary conditions at z = ± b/2 is

$${\varphi }_m(z)=-\frac{8a^2}{{\pi }^3{m}^3}\frac{\cosh \frac{m\pi z}{a}}{\cosh \frac{m\pi b}{2a}}+\frac{8a^2}{{\pi }^3{m}^3}$$

(20)

Substituting Equations (16) and (20) into (8) yields the torsional constant

$$\begin{array}{l}{I}_t=2{\displaystyle \underset{0}{\overset{a}{\int }}{\displaystyle \underset{-b/2}{\overset{b/2}{\int }}{\displaystyle \sum _{m=1,3,5\dots }\left(-\frac{8a^2}{{\pi }^3{m}^3}\frac{\cosh \frac{m\pi z}{a}}{\cosh \frac{m\pi b}{2a}}+\frac{8a^2}{{\pi }^3{m}^3}\right)\sin \frac{m\pi y}{a}}dzdy}}\\ \begin{array}{cc}& =\end{array}\frac{32b{a}^3}{{\pi }^4}{\displaystyle \sum _{m=1,3,5\dots }\frac{1}{m^4}-}\frac{64a^4}{{\pi }^5}{\displaystyle \sum _{m=1,3,5\dots }\frac{1}{m^5}\tanh \frac{m\pi b}{2a}}\end{array}$$

(21a, b)

It is well known that

$${\displaystyle \sum _{m=1,3,5\dots }\frac{1}{m^4}=\frac{{\pi }^4}{96}}$$

(22)

Substituting Equation (22) into (21b) yields

$$I_t=\frac{1}{3}b{a}^3-\frac{64a^4}{{\pi }^5}{\displaystyle \sum _{m=1,3,5\dots }\frac{1}{m^5}\tanh \frac{m\pi b}{2a}}$$

(23)

Similarly to Equation (13a, b) the torsional constant can be expressed as a function of the ratio a/b as follows

$$k\left(a/b\right)=\frac{1}{3}-\frac{64}{{\pi }^5}\frac{a}{b}{\displaystyle \sum _{m=1,3,5\dots }\frac{1}{m^5}\tanh \frac{m\pi b}{2a}}$$

(24)

The shear stresses are determined using Equations (4a, b), (7), (16), and (20)

$$\begin{array}{l}{\tau }_{xy}=-\frac{8a{M}_t}{{\pi }^2{I}_t}{\displaystyle \sum _{m=1,3,5\dots }\frac{1}{m^2}\frac{\sinh \frac{m\pi z}{a}}{\cosh \frac{m\pi b}{2a}}\sin \frac{m\pi y}{a}}\\ {\tau }_{xz}=-\frac{8a{M}_t}{{\pi }^2{I}_t}{\displaystyle \sum _{m=1,3,5\dots }\frac{1}{m^2}\left(1-\frac{\cosh \frac{m\pi z}{a}}{\cosh \frac{m\pi b}{2a}}\right)\cos \frac{m\pi y}{a}}\end{array}$$

(25a, b)

Observing that z_T = 0 and that the warping is zero along the symmetry axis y = a/2, the warping function u(x,y,z) is expressed using Equations (3a, b), (7), and (25a) as follows

$$u\left(x,y,z\right)=\frac{M_t}{G{I}_t}z\left(y-\frac{a}{2}\right)+\frac{8a^2{M}_t}{{\pi }^{3}G{I}_t}{\displaystyle \sum _{m=1,3,5\dots }\frac{1}{m^3}\frac{\sinh \frac{m\pi z}{a}}{\cosh \frac{m\pi b}{2a}}\cos \frac{m\pi y}{a}}$$

(26)

3. Results and Discussion

3.1. Analysis Using Double Sine Series for the Prandtl Stress Function

• The expressions for the Prandtl stress function, the shear stresses, and the torsional constant obtained in this study are identical to those obtained by Ike [10].

3.2. Analysis Using Single Sine Series for the Prandtl Stress Function

• The expression for the torsional constant obtained in this study is identical to that obtained by Timoshenko [8] and Sadd [9] (the beam dimensions in the latters are 2A and 2B).
• The expression for the shear stress τxy obtained in this study is identical to that obtained by Sadd (9.5.11) with τxz as the corresponding shear stress. Details are presented in Appendix A.
• For a < b the maximum stress τxz will occur at y = a and z = 0, the midpoint of the longest side. The expression for max τxz obtained in this study is identical to that obtained by Sadd (9.5.13). Details are presented in Appendix B.
• The expression for the warping obtained in this study is identical to that obtained by Sadd (9.5.14). Details are presented in Appendix C.