Skip to content
Permalink
master
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?

MSE250/stress_principal_rotations.ipynb

Go to file
 
 
Cannot retrieve contributors at this time
339 lines (339 sloc) 9.03 KB
Raw Blame
Open in GitHub Desktop
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# MSE 250 Notebook for Stress analysis, rotation, and principcal stress calculations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, we import two packages that will help us: `numpy` for working with arrays (tensors) and linear algebra, and `math` for working with mathematical functions."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import math"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Define stress tensor and transformation matrix\n",
"Here we will define two functions that allow us to easily create and manipulate stress tensors and rotation matricies.\n",
"\n",
"The function: \n",
"\n",
"$$s(\\sigma_{11}, \\sigma_{22}, \\sigma_{33}, \\sigma_{12}, \\sigma_{13}, \\sigma_{23})$$\n",
"\n",
"returns a symmetric 3x3 numpy array with the $\\sigma_{ij}$ components properly placed.\n",
"\n",
"The function:\n",
"\n",
"$$Q(\\phi_1, \\Phi, \\phi_2[, deg=True])$$\n",
"\n",
"returns a 3x3 transformation matrix using the Bunge convention. The inputs are $\\phi_1$ - clockwise rotation about z axis [+/-180deg], $\\Phi$ - clockwise rotation about x' [0-180deg], and $\\phi-2$ - clockwise rotation along z'' [+/-180deg]. These three angle represent the Euler angles, and are by default input as degrees. The optional input deg can be set to False to allow for Euler angles in units of radians.\n",
"\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 95,
"metadata": {},
"outputs": [],
"source": [
"def s(s11,s22,s33,s12,s13,s23):\n",
" stress = np.array([[s11, s12, s13],\n",
" [s12, s22, s23],\n",
" [s13, s23, s33]])\n",
" return stress\n",
"\n",
"def Q(phi1,Phi,phi2,deg=True):\n",
" if deg == True:\n",
" phi1 = phi1*math.pi/180\n",
" Phi = Phi*math.pi/180\n",
" phi2 = phi2*math.pi/180\n",
" aphi2 = np.array([[ np.cos(phi2),np.sin(phi2),0],\n",
" [-np.sin(phi2),np.cos(phi2),0],\n",
" [ 0, 0,1]])\n",
"\n",
" aPhi = np.array([[1, 0, 0],\n",
" [0, np.cos(Phi),np.sin(Phi)],\n",
" [0,-np.sin(Phi),np.cos(Phi)]])\n",
"\n",
" aphi1 = np.array([[ np.cos(phi1),np.sin(phi1),0],\n",
" [-np.sin(phi1),np.cos(phi1),0],\n",
" [ 0, 0,1]])\n",
"\n",
" Q = np.dot(aphi2,np.dot(aPhi,aphi1))\n",
" return Q"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 1) 45deg coordinate transformation about z axis\n",
"First we will begin with a biaxial stress state:"
]
},
{
"cell_type": "code",
"execution_count": 172,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[10 0 0]\n",
" [ 0 -5 0]\n",
" [ 0 0 0]]\n"
]
}
],
"source": [
"s_1 = s(10,-5,0,0,0,0)\n",
"print(s_1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, let's transform the coordinates about the z axis 45 degrees. Since we are already in the principal orientation, this transformatio should put us in the orietation for maximum shear. The formula we use for transformation is:\n",
"\n",
"$$\\sigma'=Q \\cdot \\sigma \\cdot Q^T$$\n",
"$$\\sigma_{mn}'=Q_{mi}Q_{nj}\\sigma_{ij}$$\n",
"\n",
"We will use numpy's dot method to evaluate the first equation directly:"
]
},
{
"cell_type": "code",
"execution_count": 173,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 2.5 -7.5 0. ]\n",
" [-7.5 2.5 0. ]\n",
" [ 0. 0. 0. ]]\n"
]
}
],
"source": [
"s_1p = np.dot(Q(0,0,45),np.dot(s_1,np.transpose(Q(0,0,45))))\n",
"print(s_1p.round(decimals=4))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 2) Find the principal stresses and orientation, and maximum shear stress\n",
"We can find the eigenvalues (principal stresses) and eigenvectors (principal directions) of the newly-transformed stress matrix. We use the numpy.linalg.eig() method to find these, and it returns the principal stress list and principal direction matrix."
]
},
{
"cell_type": "code",
"execution_count": 174,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Principal Stresses\n",
"s_1 = 10.0 \n",
"s_2 = -5.0 \n",
"s_3 = 0.0\n",
"\n",
"Principal Directions\n",
"v_1 = [ 0.70710678 -0.70710678 0. ]\n",
"v_2 = [0.70710678 0.70710678 0. ]\n",
"v_3 = [0. 0. 1.]\n"
]
}
],
"source": [
"p_s,p_d = np.linalg.eig(s_1p)\n",
"\n",
"print('Principal Stresses')\n",
"print('s_1 = %s \\ns_2 = %s \\ns_3 = %s' % tuple(p_s))\n",
"print('')\n",
"print('Principal Directions')\n",
"print('v_1 = '+str(p_d[:,0]))\n",
"print('v_2 = '+str(p_d[:,1]))\n",
"print('v_3 = '+str(p_d[:,2]))"
]
},
{
"cell_type": "code",
"execution_count": 175,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"tau_max = 7.5\n"
]
}
],
"source": [
"t_max = (max(p_s)-min(p_s))/2.\n",
"print('tau_max = {}'.format(t_max))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 3) Use cosines to transform to the principal coordinate system\n",
"So far we have rotated 45 degrees around the z axis from the xyz coorindate system to the new x'y'z' system. We then calculated the principal stresses and principal directions from this new coordinate frame. Let's now transform to a new coordinate system of the principal directions (back to our original coordinate system) using the law of cosines to find our new transformation matrix.\n",
"\n",
"The Einstein notation for this transformation matrix is:\n",
"\n",
"$$Q_{ij}=P_{mi}I_{mj}$$\n",
"\n",
"where P is the new coordinate system with columns as the coordinate axes, and I is the current coordinate system with columns as the coordinate axes. In our case, the new coordinate system is our principal directions that we just calculated, and I is the identify matrix, because we are currently in a regular Cartesian coordinate frame."
]
},
{
"cell_type": "code",
"execution_count": 176,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[1. 0. 0.]\n",
" [0. 1. 0.]\n",
" [0. 0. 1.]]\n"
]
}
],
"source": [
"#Current coordinate frame, columns are axes\n",
"I=np.identity(3)\n",
"print(I)"
]
},
{
"cell_type": "code",
"execution_count": 177,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 0.70710678 0.70710678 0. ]\n",
" [-0.70710678 0.70710678 0. ]\n",
" [ 0. 0. 1. ]]\n"
]
}
],
"source": [
"#New coordinate frame, columns are axes\n",
"print(p_d)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's find the new transformation matrix using Einstein summation:"
]
},
{
"cell_type": "code",
"execution_count": 178,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 0.70710678 -0.70710678 0. ]\n",
" [ 0.70710678 0.70710678 0. ]\n",
" [ 0. 0. 1. ]]\n"
]
}
],
"source": [
"Qp = np.einsum('mi,mj->ij',p_d,I)\n",
"print(Qp)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's rotate back:"
]
},
{
"cell_type": "code",
"execution_count": 179,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[10. -0. 0.]\n",
" [-0. -5. 0.]\n",
" [ 0. 0. 0.]]\n"
]
}
],
"source": [
"s_1pp = np.dot(Qp,np.dot(s_1p,np.transpose(Qp)))\n",
"print(s_1pp.round(decimals=4))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 4
}