PDF (letter size)

PDF (legal size)

Solving the torsion problem for isotropic matrial with a rectangular cross section using the FEM and FVM methods with triangular elements
MAE 207, Computational methods. UCI
Fall 2006

Nasser M. Abbasi

2006   Compiled on September 4, 2021 at 6:36pm


 1 Introduction
  1.1 Problem setup
 2 Analytical solution using Prandtl stress function
  2.1 Stress components
  2.2 Strain components
  2.3 Determining the twist angle \(\alpha \)
  2.4 Displacement calculations
 3 References

1 Introduction

We consider bar made of isotropic martial with rectangular cross section subjected to twisting torque \(T\). The following diagram illustrate the basic geometry.

Figure 1:basic geometry

Experiments show that rectangular cross sections do wrap and that cross sections do not remain plane as shown in this diagram (in the case of a circular cross section, cross section do NOT wrap).

Figure 2:cross section wrap

This is another diagram showing a bar under torsion

Figure 3:bar under torsion

1.1 Problem setup

1.1.1 What are the assumptions?

  1. The twist rate (called \(k\) in this problem) and defined as \(\frac {d\alpha }{dz}\) where \(\alpha \) is the twist angle is assumed to be constant.

  2. Cross section can wrap also in the \(z\) direction (i.e. the cross section does not have to remain in the \(xy\) plane) but if this happens, all cross sections will wrap in the \(z\) section by the same amount.

  3. Material is isotropic

1.1.2 What is the input and what is the output?

The input to the problem are the following (these are the known or given):

  1. The width \(b\) and height \(a\) of the cross section.

  2. Material Modulus of rigidity or sheer modulus \(G\) which is the ratio of the shearing stress \(\tau \) to the shearing strain \(\gamma \)

  3. The applied torque \(T\)

  4. \(J\) the torsion constant for the a rectangular cross section. For a rectangular section of dimensions \(a,b\) it is given by

    \begin {equation} J=\frac {16}{3}a^{3}b\left (1-\frac {192\ a}{b\pi ^{5}}{\displaystyle \sum \limits _{n=1,3,5\cdots }^{\infty }} \frac {1}{n^{5}}\tanh \left (\frac {n\ \pi b}{2\ a}\right ) \right ) \tag {1} \end {equation}

Hence the torsional rigidity \(GJ\) is known since \(G\) is given (material) and \(J\) is from above (geometry).

1.1.3 The output from the problem (the things we need to calculate)

  1. The stress distribution in the cross section (stress tensor field). Once this is found then using the material constitutive relation we can the strain tensor field.

  2. The angle of twist \(\alpha \) as a function of \(z\) (the length of the beam).

2 Analytical solution using Prandtl stress function

First we solve for the Prandtl stress function \(\Phi \left (x,y\right ) \) by solving the Poisson equation

\[ \nabla ^{2}\Phi \left (x,y\right ) =-2GK \]

Where \(G\) is the sheer modulus and \(k\) is the twist rate (which was assumed to be constant).

The boundary conditions (\(\Phi \left (x,y\right ) \) at any point on the edge of the cross section and at the ends of the beam) is an arbitrary constant. We take this constant to be zero. Hence at the cross section boundary we have\[ \Phi =0 \]

The analytical solution to the above equation is from book Theory of elasticity by S. P. Timoshenko and J. N. Goodier

\begin {equation} \Phi \left (x,y\right ) =\frac {32\ Gk\ a^{2}}{\pi ^{3}}{\displaystyle \sum \limits _{n=1,3,5,\cdots }^{\infty }} \frac {1}{n^{3}}\left (-1\right ) ^{\frac {\left (n-1\right ) }{2}}\left ( 1-\frac {\cosh \left (\frac {n\pi y}{2a}\right ) }{\cosh \left (\frac {n\pi b}{2a}\right ) }\right ) \cos \left (\frac {n\pi x}{2a}\right ) \tag {2} \end {equation}

where the linear twist \(k\)

\[ k=\frac {T}{GJ} \]

Hence (2) becomes

\[ \Phi \left (x,y\right ) =\frac {32\ T\ a^{2}}{J\pi ^{3}}{\displaystyle \sum \limits _{n=1,3,5,\cdots }^{\infty }} \frac {1}{n^{3}}\left (-1\right ) ^{\frac {\left (n-1\right ) }{2}}\left ( 1-\frac {\cosh \left (\frac {n\pi y}{2a}\right ) }{\cosh \left (\frac {n\pi b}{2a}\right ) }\right ) \cos \left (\frac {n\pi x}{2a}\right ) \]

Where \(J\) is given by (1)

2.1 Stress components

\begin {align*} \tau _{yz} & =-\frac {\partial \Phi }{\partial x}\\ \tau _{xz} & =\frac {\partial \Phi }{\partial y}\\ \tau _{yx} & =0\\ \sigma _{x} & =\sigma _{y}=\sigma _{z}=0 \end {align*}

Hence \[ \tau _{yz}=-\frac {\partial \Phi }{\partial x}=\frac {16Gka}{\pi ^{2}}{\displaystyle \sum \limits _{n=1,3,5,\cdots }^{\infty }} \frac {1}{n^{2}}\left (-1\right ) ^{\frac {\left (n-1\right ) }{2}}\left ( 1-\frac {\cosh \left (\frac {n\pi y}{2a}\right ) }{\cosh \left (\frac {n\pi b}{2a}\right ) }\right ) \sin \left (\frac {n\pi x}{2a}\right ) \]


\[ \tau _{xz}=\frac {\partial \Phi }{\partial y}=\frac {16\ T\ a}{J\pi ^{2}}{\displaystyle \sum \limits _{n=1,3,5,\cdots }^{\infty }} \frac {1}{n^{2}}\left (-1\right ) ^{\frac {\left (n-1\right ) }{2}}\left ( -\cos \left (\frac {n\pi x}{2a}\right ) sech\left (\frac {bn\pi }{2a}\right ) \sinh \left (\frac {n\pi y}{2a}\right ) \right ) \]

Timoshenko gives the maximum sheer stress, which is \(\sqrt {\tau _{yz}^{2}+\tau _{xz}^{2}}\) as\[ \tau _{\max }=\frac {16Gka}{\pi ^{2}}{\displaystyle \sum \limits _{n=1,3,5,\cdots }^{\infty }} \frac {1}{n^{2}}\left (1-\frac {1}{\cosh \left (\frac {n\pi b}{2a}\right ) }\right ) \]

2.2 Strain components

Given that \(E\) is Young’s modulus for the material, \(\upsilon \) is Poisson’s ratio for the material, and \(G=\frac {E}{2\left (1+\upsilon \right ) }\) we can now obtain the strain components from the constitutive equations (stress-strain equations) since we have determined the stress components from the above solution.

\begin {align*} \varepsilon _{x} & =\frac {1}{E}\left (\sigma _{x}-\upsilon \left (\sigma _{y}+\sigma _{z}\right ) \right ) =0\\ \varepsilon _{y} & =\frac {1}{E}\left (\sigma _{y}-\upsilon \left (\sigma _{x}+\sigma _{z}\right ) \right ) =0\\ \varepsilon _{z} & =\frac {1}{E}\left (\sigma _{z}-\upsilon \left (\sigma _{x}+\sigma _{y}\right ) \right ) =0\\ \gamma _{xy} & =\frac {2\left (1+\upsilon \right ) }{E}\tau _{xy}=\frac {1}{G}\tau _{xy}=0\\ \gamma _{yz} & =\frac {2\left (1+\upsilon \right ) }{E}\tau _{yz}=\frac {1}{G}\tau _{yz}\\ \gamma _{xz} & =\frac {2\left (1+\upsilon \right ) }{E}\tau _{xz}=\frac {1}{G}\tau _{xz} \end {align*}

Hence only \(\gamma _{yz}\) and \(\gamma _{xz}\) are non-zero.

2.3 Determining the twist angle \(\alpha \)

If we look at a cross section of the bar at some distance \(z\) from the end of the bar, the angle that this specific cross section has twisted due to the torque is \(\alpha \).

Figure 4:twist angles

This angle is given by the solution to the equation

\[ \frac {d\alpha \relax (z) }{dz}=k \]

But \(k\) is the linear twist and is given by \(k=\frac {T}{GJ}\) hence the above equation becomes

\[ \frac {d\alpha \relax (z) }{dz}=\frac {T}{GJ} \]


\[ \alpha \relax (z) =\frac {T}{GJ}z+C_{1} \]

Where \(C_{1}\) is the constant of integration. Assuming \(\alpha =0\) at \(z=0\) we obtain that \[ \alpha \relax (z) =\frac {T}{GJ}z \] and using the expression \(J\) given in equation (1) above we can determine \(\alpha \) for each \(z.\)

2.4 Displacement calculations

Figure 5:Displacement calculations

\begin {align*} r & =\sqrt {x^{2}+y^{2}}\\ u & =r\alpha \left (-\sin \beta \right ) \\ v & =r\alpha \left (\cos \beta \right ) \end {align*}

we see that

\begin {align*} \sin \beta & =\frac {y}{r}\\ \cos \beta & =\frac {x}{r} \end {align*}


\begin {align*} u & =-\alpha y\\ v & =\alpha x \end {align*}

Where \(\alpha =\frac {T}{GJ}z\)

3 References

  1. Mathematica Structural Mechanics help page

  2. MIT course 16.20 lecture notes. MIT open course website.

  3. Theory of elasticity by S. P. Timoshenko and J. N. Goodier. chapter 10