extension_due_to_weight_partially_completed.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we introduced the Galerkin method for deriving element equations, we worked through the example of a rod hanging under its own weight defined by the following differential equation:\n",
    "\n",
    "![](./images/hanging_rod_problem.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After applying the Galerkin method, we obtained the following element matrix equation:\n",
    "\n",
    "![](./images/hanging_rod_element_equation.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this program, we'll write a complete finite element solver that will let us look at the displacement and stress convergance properties as we vary the number of elements in the mesh. These results will be compared to the exact analytical solution for displacement $\\delta(x)=\\frac{q_0}{AE}(-\\frac{x^2}{2} + L x)$ and stress $\\sigma(x)=\\frac{q_0}{A}(-x + L)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the constants needed for the problem\n",
    "L = 1.0  # Total length of rod\n",
    "A = 0.2  # Cross sectional area of rod\n",
    "E = 10   # Elastic modulus of rod\n",
    "q_0 = 5.0*A   # Downward force per unit length of rod"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the number of elements to use in the mesh\n",
    "M = 2  # Number of elements\n",
    "N = M+1  # Number of nodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the nodes and the connectivity matrix\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the element matrices\n",
    "k_elements = []\n",
    "f_elements = []\n",
    "for current_elem in range(M):\n",
    "    node1 = nodes[connectivity[current_elem,0]]\n",
    "    node2 = nodes[connectivity[current_elem,1]]\n",
    "    elem_length = node2-node1\n",
    "    k_elements.append( E*A*np.array([[1/elem_length, -1/elem_length],[-1/elem_length, 1/elem_length]]) )\n",
    "    \n",
    "    f_elements.append( q_0*np.array([elem_length/2.0, elem_length/2.0]) )\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assemble the matrices\n",
    "K = np.zeros((N,N))\n",
    "F = np.zeros((N,1))\n",
    "\n",
    "for elem_num, elem_k in enumerate(k_elements):\n",
    "    for i in range(2):\n",
    "        row = connectivity[elem_num,i]\n",
    "        F[row] += f_elements[elem_num][i]\n",
    "        for j in range(2):\n",
    "            col = connectivity[elem_num,j]\n",
    "            K[row,col] += elem_k[i,j]\n",
    "        \n",
    "(K,F)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply the zero displacement boundary condition at the first node\n",
    "K_with_bc = np.copy(K) # make a copy of the original K, we need it to calculate the reaction force\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solve for the nodal displacements\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solve for the reaction force at node 1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the stress for each of the elements\n",
    "element_stress = []\n",
    "\n",
    "for node1,node2 in connectivity:\n",
    "    elongation = (delta[node2]-delta[node1])/(nodes[node2]-nodes[node1]) # strain = (delta L) / L\n",
    "    element_stress.append(E*elongation) # stress = E*strain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Displacement plot\n",
    "plt.plot(nodes, delta, 'b', label='FE Solution')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('displacement')\n",
    "plt.title('FE vs. Exact Displacement')\n",
    "x_values = np.linspace(0,L)\n",
    "plt.plot(x_values, (1/(A*E)) * (-((q_0*(x_values**2))/2.0) + q_0*L*x_values), 'g--', label='Exact Solution')\n",
    "plt.legend(loc='best')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Stress plot\n",
    "for (elem_num, (node1,node2)) in enumerate(connectivity):\n",
    "    if elem_num == 0:\n",
    "        plt.plot([nodes[node1],nodes[node2]], [element_stress[elem_num], element_stress[elem_num]],'b',label='FE Solution')\n",
    "    else:\n",
    "        plt.plot([nodes[node1],nodes[node2]], [element_stress[elem_num], element_stress[elem_num]],'b')\n",
    "\n",
    "plt.plot(x_values, ((q_0*L)/(A))*(1-x_values/L), 'g--', label='Exact Solution')\n",
    "\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('stress')\n",
    "plt.title('FE vs. Exact Stress')\n",
    "plt.legend(loc='best')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}