Nonlinear_vibrating_cantilever.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using TensorFlow backend.\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import random\n",
    "import math\n",
    "\n",
    "import numpy as np\n",
    "np.random.seed(1)\n",
    "from keras.models import Sequential\n",
    "from keras.layers import Dense\n",
    "from keras.optimizers import SGD\n",
    "#%matplotlib notebook\n",
    "\n",
    "from numpy.linalg import inv\n",
    "from numpy import linalg as LA\n",
    "from scipy import linalg\n",
    "import random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_nodes = 10\n",
    "n_dofs = n_nodes*2\n",
    "n_elements = n_nodes-1\n",
    "x = np.linspace(0, 1, n_nodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = 1000 #mm\n",
    "rho = 0.006 #kg/mm\n",
    "E = 210000 #N/mm2\n",
    "I = 0.801*10**6 #mm4\n",
    "A = 764 #mm^2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Linear element stiffnes matrix\n",
    "L = L/n_elements\n",
    "k_lin = np.matrix([[12, -6*L, -12, -6*L], [-6*L, 4*L**2, 6*L, 2*L**2], [-12, 6*L, 12, 6*L], [-6*L, 2*L**2, 6*L, 4*L**2]])*E*I/L**3\n",
    "\n",
    "#Element mass matrix\n",
    "m = np.matrix([[156, 22*L, 54, -13*L],[22*L, 4*L**2, 13*L, -3*L**2], [54, 13*L, 156, -22*L], [-13*L, -3*L**2, -22*L, 4*L**2]] )*rho*A*L/420"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Global linear stiffness matrix\n",
    "n_dofs = n_nodes*2\n",
    "K = np.zeros([n_dofs, n_dofs])\n",
    "for i in range(n_elements):\n",
    "    K[2*i:2*i+4,2*i:2*i+4] += k_lin\n",
    "    \n",
    "#Global mass matrix\n",
    "M = np.zeros([n_dofs, n_dofs])\n",
    "for i in range(n_elements):\n",
    "    M[2*i:2*i+4,2*i:2*i+4] += m\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Parameters\n",
    "gamma = 0.5\n",
    "beta = 0.25\n",
    "dt = 0.01\n",
    "n_timesteps = 1000\n",
    "C = np.zeros([n_dofs, n_dofs])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Fixed end at x=0\n",
    "K = K[2:,2:]\n",
    "M = M[2:,2:]\n",
    "C = C[2:,2:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Linear matrices\n",
    "\n",
    "A1 = M + gamma*dt*C + beta*dt**2*K\n",
    "A2 = -2*M + (1-2*gamma)*dt*C + (0.5-2*beta+gamma)*dt**2*K\n",
    "A3 = M - (1-gamma)*dt*C + (0.5+beta-gamma)*dt**2*K\n",
    "\n",
    "A1inv = inv(A1)\n",
    "Avec = [A1,A2,A3,A1inv]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Linear force vector\n",
    "def get_Fh(Fvec):\n",
    "    F0 = Fvec[0]\n",
    "    F1 = Fvec[1]\n",
    "    F2 = Fvec[2]\n",
    "    \n",
    "    Fh = dt**2*(beta*F2 + (0.5-2*beta*gamma)*F1 + (0.5+beta-gamma)*F0)\n",
    "    \n",
    "    return Fh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Linear and nonlinear solution\n",
    "def get_Qd2(Fh, Qvec, Avec):\n",
    "    A1 = Avec[0]\n",
    "    A2 = Avec[1]\n",
    "    A3 = Avec[2]\n",
    "    A1inv = Avec[3]\n",
    "    \n",
    "    Qd0 = Qvec[0]\n",
    "    Qd1 = Qvec[1]\n",
    "    Qd2 = np.matmul(A1inv, np.matmul(A2, Qd1)) - np.matmul(A1inv, np.matmul(A3, Qd0)) + dt**2*np.matmul(A1inv, Fh)\n",
    "    \n",
    "    return Q2d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Linear force vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'c' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-15-704127addf2d>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;31m#Nonlinear matrices\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      2\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m \u001b[0ma1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mM\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mgamma\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mdt\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mc\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mbeta\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mdt\u001b[0m\u001b[1;33m**\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      4\u001b[0m \u001b[0ma2\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m-\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mM\u001b[0m \u001b[1;33m+\u001b[0m \u001b[1;33m(\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'c' is not defined"
     ]
    }
   ],
   "source": [
    "#Nonlinear matrices\n",
    "\n",
    "a1 = M + gamma*dt*c + beta*dt**2*k\n",
    "a2 = -2*M + (1-2*gamma)*dt*c + (0.5-2*beta+gamma)*dt**2*k\n",
    "A3 = M - (1-gamma)*dt*c + (0.5+beta-gamma)*dt**2*k\n",
    "\n",
    "A1inv = inv(A1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Solution loop\n",
    "error = 1\n",
    "TOL = 0.1\n",
    "\n",
    "#Initial conditions\n",
    "Qd0 = np.zeros(n_dofs-2)\n",
    "Qd1 = np.zeros(n_dofs-2)\n",
    "qd0 = np.zeros(n_dofs-2)\n",
    "qd1 = np.zeros(n_dofs-2)\n",
    "\n",
    "J = \n",
    "\n",
    "#Timestep loop\n",
    "for i in range(n_timesteps):\n",
    "    \n",
    "    #Linear solution\n",
    "    \n",
    "    #Nonlinear solution\n",
    "    \n",
    "    #Error\n",
    "    \n",
    "    #Convergence loop\n",
    "    while error > TOL:\n",
    "    \n",
    "    #Q = q\n",
    "    \n",
    "    #New nonlinear soltuion\n",
    "    \n",
    "    #new error\n",
    "    \n",
    "    "
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}