Principal stresses and transformation¶
Problem statement¶
Find the pricipal stresses and directions for the plane stress state
\begin{equation} \boldsymbol{\sigma}= \begin{bmatrix} \sigma_{x} & \tau_{xy} & \tau_{xz} \\ \tau_{xy} & \sigma_{y} & \tau_{yz} \\ \tau_{xz} & \tau_{yz} & \sigma_{z} \end{bmatrix}= \begin{bmatrix} 200 & -75 & 0 \\ -75 & 50 & 0 \\ 0 & 0 & 0 \end{bmatrix} \end{equation}by using the relations for stress transformation.
Solution¶
Write the stress state as a 1-dimensional array:
import numpy as np
sxyz=np.array([200,50,0,0,0,-75])
The transformation matrix for stress by rotation about the z-axis as implemented in Transformations:
from compositelib import T3Dsz
Then, a straight forward sweep through rotation angles from -90 to 90 in order to find where the shear stress is zero, or alternatively, where normal stresses show maximum or minimum:
s_xyz=(200,50,0,0,0,-75)
angles=np.linspace(-90,90,10000)
s11,s22,s12=[],[],[]
for a in angles:
s_123=np.dot(T3Dsz(a),s_xyz)
s11.append(s_123[0])
s22.append(s_123[1])
s12.append(s_123[5])
import matplotlib.pyplot as plt
%matplotlib inline
fig,ax=plt.subplots(figsize=(10,4))
ax.plot(angles,s11,color='red', label='$\sigma_1$')
ax.plot(angles,s22,color='blue', label='$\sigma_2$')
ax.plot(angles,s12,color='green', label=r'$\tau_{12}$')
ax.set_xlim(-90,90)
ax.legend(loc='best')
ax.grid(True)
plt.show()
max_s11 = max(s11)
ind_angle1 = s11.index(max_s11)
angle1=angles[ind_angle1]
print('Angle for max S11: {}'.format(angle1))
max_s22 = max(s22)
ind_angle2 = s22.index(max_s22)
angle2=angles[ind_angle2]
print('Angle for max S22: {}'.format(angle2))
Verify that the shear stresses are zeros:
T=T3Dsz(angle1)
s=np.dot(T,sxyz)
print(np.round(s,1))
T=T3Dsz(angle2)
s=np.dot(T,sxyz)
print(np.round(s,1))
Hence, the first principal stress is 231.1, the second is 18.9 and the third is zero. The direction of the first principal stress is inclined by -22.5 degrees in the x-y plane relatively to the x-axis.
Verify using numpy.linalg.eig() (see Stress)
ST=np.array ( [ [200.0, -75.0, 0.0 ],
[ -75.0, 50.0, 0.0 ],
[ 0.0, 0.0, 0.0 ] ])
L,v = np.linalg.eig(ST)
print('Eigenvalues:')
print(L)
print()
print('Eigenvectors:')
print(v)
Remember that the eigenvectors are given by the columns:
v1 = v[:,0]
v2 = v[:,1]
v3 = v[:,2]
import matplotlib.pyplot as plt
%matplotlib inline
fig,ax = plt.subplots(figsize=(4,4))
ax.set_axis_off()
ax.set_xlim(-1.1,1.1)
ax.set_ylim(-1.1,1.1)
ax.arrow(0, 0, 1, 0.0, head_width=0.1, head_length=0.1, fc='black', ec='black')
ax.arrow(0, 0, 0.0, 1, head_width=0.1, head_length=0.1, fc='black', ec='black')
ax.text(0.9, 0.1, 'x')
ax.text(0.1, 0.9, 'y')
ax.arrow(0, 0, v1[0], v1[1], head_width=0.1, head_length=0.1, fc='red', ec='red')
ax.arrow(0, 0, v2[0], v2[1], head_width=0.1, head_length=0.1, fc='blue', ec='blue')
plt.show()
