{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimation of stiffness properties by classical laminate theory (CLT)\n",
    "\n",
    "This python notebook goes through the steps of determining the stiffness of a fibre composite with a given layup.\n",
    "\n",
    "Initially we initiate the packages necessary for the code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Definition of input parameters\n",
    "\n",
    "Then, we will define the necessary input parameters for the calculations. Note that these stiffnesses are for each layer/ply in their own local coordinate system where the 1-direction is in the direction of the fibres and the 2-direction is in the direction transverse to the fibres. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "layup = np.array([45,-45,0])                     # Orientation of the fibres for each layer/ply of the laminate\n",
    "\n",
    "E1 = np.full((layup.shape[0],1), 40.)[:,0]       # Axial stiffness of each ply [GPa]\n",
    "#E1 = np.array([45.6,45.6,45.6])                 # (Use this instead if layer properties are not the same for all layers)\n",
    "\n",
    "E2 = np.full((layup.shape[0],1), 9.8)[:,0]       # transverse stiffness of each ply [GPa]\n",
    "#E2 = np.array([16.2,16.2.16,2])                 # (Use this instead if layer properties are not the same for all layers)\n",
    "        \n",
    "v12 = np.full((layup.shape[0],1), 0.3)[:,0]      # Poisson's ratio\n",
    "#v12 = np.array([0.278,0.278,0.278])             # (Use this instead if layer properties are not the same for all layers)\n",
    "v21 = (v12*E2)/E1                                # (Since the compliance matrix is symmetric, v12/E1=v21/E2)\n",
    "\n",
    "G12 = np.full((layup.shape[0],1), 2.8)[:,0]      # Shear modulus [GPa]\n",
    "#G12 = np.array([5.83,5.83,5.83])                # (Use this instead if layer properties are not the same for all layers)\n",
    "\n",
    "h0 = 0.1\n",
    "h = np.array([4*h0,4*h0,2*h0])                   #height of each ply [mm]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Estimate the stiffness contribution of each layer of the laminate\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hook's Law for a thin orthotropic layer (a ply) can be written as:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\left( \\begin{array}{c}\n",
    "\\sigma_1  \\\\\n",
    "\\sigma_2  \\\\\n",
    "\\tau_{12}  \\end{array} \\right) =\n",
    "\\frac{1}{1-\\nu_{12}\\nu_{21}}\n",
    "\\left[ \\begin{array}{ccc}\n",
    "E_1 & \\nu_{21}E_1 & 0\\\\\n",
    "\\nu_{12}E_2 & E_2 & 0 \\\\\n",
    "0 & 0 & G_{12}(1-\\nu_{12}\\nu_{21}) \\end{array}  \\right] \n",
    "\\left( \\begin{array}{c}\n",
    "\\varepsilon_1 \\\\\n",
    "\\varepsilon_2  \\\\\n",
    "\\gamma_{12}  \\end{array} \\right)\n",
    "\\end{equation*} \n",
    "\n",
    "or in a shorther form\n",
    "\n",
    "\\begin{equation*}\n",
    "\\sigma_l=Q_l\\varepsilon_l\n",
    "\\end{equation*} \n",
    "\n",
    "where the $l$ stands for lamina or local, and Q is the local stiffness matrix of the ply. \n",
    "\n",
    "To estimate the stiffness contribution of each layer of the composite, we calculate the local stiffness matrix for each ply and then transform them into the global coordinates. To do so, we start by calculating the local stiffness matrices for each layer of the composite. This is done by looping through each layer of the laminate as was defined by the \"layup\" parameter above, calculating the local stiffness matrices, and transforming them into global coordinates using the transformation matrix:\n",
    "\n",
    "\\begin{equation*}\n",
    "T=\\left[ \\begin{array}{ccc}\n",
    "c^2 & s^2 & -2sc\\\\\n",
    "s^2 & c^2 & 2sc \\\\\n",
    "sc & -sc & c^2-s^2 \\end{array}  \\right] \n",
    "\\end{equation*} \n",
    "\n",
    "where s and c refer to $\\sin\\theta$ and $\\cos\\theta$, respectively. The global stiffness contribution from each layer can be found from:\n",
    "\n",
    "\\begin{equation*}\n",
    "Q=TQ_lT^T\n",
    "\\end{equation*} \n",
    "\n",
    "Using this approach, Q is calculated for each layer in the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ply no. 1: theta=45, E1l=40.0  E2l=9.8\n",
      "[[ 17.03385654  11.43385654   7.7202311 ]\n",
      " [ 11.43385654  17.03385654   7.7202311 ]\n",
      " [  7.7202311    7.7202311   11.22756787]]\n",
      "Ply no. 2: theta=-45, E1l=40.0  E2l=9.8\n",
      "[[ 17.03385654  11.43385654  -7.7202311 ]\n",
      " [ 11.43385654  17.03385654  -7.7202311 ]\n",
      " [ -7.7202311   -7.7202311   11.22756787]]\n",
      "Ply no. 3: theta=0, E1l=40.0  E2l=9.8\n",
      "[[ 40.9018866    3.00628867   0.        ]\n",
      " [  3.00628867  10.02096222   0.        ]\n",
      " [  0.           0.           2.8       ]]\n"
     ]
    }
   ],
   "source": [
    "# Initiating Q as a list that can contain the stiffness matrices transformed into the global coordinates for each layer\n",
    "Q=[]\n",
    "# Looping through each layer of the laminate\n",
    "for i in range(layup.shape[0]):\n",
    "    theta = layup[i]*np.pi/180 # Current ply angle changed into radians\n",
    "    # assigning the current material properties to temporary variables (could also be used directly):\n",
    "    E1l=E1[i]\n",
    "    E2l=E2[i]\n",
    "    v12l=v12[i]\n",
    "    v21l=v21[i]\n",
    "    G12l=G12[i]\n",
    "    print('Ply no. ' + str(i+1) + ': theta='+ str(layup[i]) +', E1l=' + str(E1l), ' E2l=' + str(E2l))\n",
    "    # Establishing current local stiffness matrix, Ql:\n",
    "    Ql = 1/(1-v12l*v21l)*np.array([[E1l,v21l*E1l,0],\\\n",
    "                                   [v12l*E2l,E2l,0],\\\n",
    "                                   [0,0,G12l*(1-v12l*v21l)]])\n",
    "    # Transformation matrix:\n",
    "    T=np.array([[np.cos(theta)**2,np.sin(theta)**2,-2*np.sin(theta)*np.cos(theta)], \\\n",
    "                [np.sin(theta)**2,np.cos(theta)**2,2*np.sin(theta)*np.cos(theta)],\\\n",
    "                [np.sin(theta)*np.cos(theta),-np.sin(theta)*np.cos(theta),np.cos(theta)**2-np.sin(theta)**2]])\n",
    "    # Adding the current stiffness matrix in the global coordinate system to the Q-list variable:\n",
    "    Q.append(np.dot(np.dot(T,Ql),np.transpose(T)))\n",
    "    print(Q[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the stiffness matrices calculated in the global coordinate system for each of the layer orientations, the final contribution is calculated through the following:\n",
    "\n",
    "\\begin{equation*}\n",
    "A= \\sum _{i=1}^N Q_i h_i\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A=np.zeros((3,3))\n",
    "for i in range(layup.shape[0]):\n",
    "    A=A+Q[i]*h[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Calculating the average stiffness matrix and the global properties\n",
    "The average stiffness matrix $Q^*$ can then be found from:\n",
    "$$\n",
    "Q_*=\\frac{A}{h_{tot}}\n",
    "$$\n",
    "followed by the compliance matrix from which we can obtain the global stiffness properties of the total laminate as follows:\n",
    "$$\n",
    "S^*=(Q^*)^{-1} \\rightarrow E_x=1/S_{11}, \\quad E_y = 1/S_{22}, \\quad G_{xy}=1/S_{33} \\quad and \\quad \\nu_{xy}=\\frac{-S_{21}}{S_{1,1}}\n",
    "$$\n",
    "\n",
    "This is done in the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Q* Matrix\n",
      "[[ 21.80746255   9.74834296   0.        ]\n",
      " [  9.74834296  15.63127767   0.        ]\n",
      " [  0.           0.           9.5420543 ]]\n",
      "S* Matrix\n",
      "[[ 0.06358098 -0.03965186  0.        ]\n",
      " [-0.03965186  0.08870292  0.        ]\n",
      " [ 0.          0.          0.10479924]]\n"
     ]
    }
   ],
   "source": [
    "Qstar = A/sum(h)\n",
    "Sstar = np.linalg.inv(Qstar)\n",
    "Ex = 1/Sstar[0,0]\n",
    "Ey = 1/Sstar[1,1]\n",
    "Gxy = 1/Sstar[2,2]\n",
    "vxy=-Sstar[1,0]/Sstar[0,0]\n",
    "\n",
    "# Printing out the matrices\n",
    "print('Q* Matrix')\n",
    "print(Qstar)\n",
    "print('S* Matrix')\n",
    "print(Sstar)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The total stiffness of the laminate is Ex=15.7279729194 and Ey=11.2735863413\n",
      "The shear modulus is Gxy=9.54205429725 and the poissons ratio of vxy=0.623643387729\n"
     ]
    }
   ],
   "source": [
    "print('The total stiffness of the laminate is Ex='+ str(Ex)+ ' and Ey=' + str(Ey) + '\\n' + 'The shear modulus is Gxy='+ str(Gxy) + ' and the poissons ratio of vxy='+str(vxy))"
   ]
  }
 ],
 "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}