active_earth_pressure.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false
   },
   "source": [
    "<div style=\"display: flex; align-items: center;\">\n",
    "    <img src=\"https://github.com/nagelt/Teaching_Scripts/raw/9d9e29ecca4b04eaf7397938eacbf116d37ddc93/Images/TUBAF_Logo_blau.png\" width=\"500\" height=\"auto\" height=\"auto\" style=\"margin-right: 100px;\" />\n",
    "    <div>\n",
    "        <p><strong>Prof. Dr. Thomas Nagel</strong></p>\n",
    "        <p>Chair of Soil Mechanics and Foundation Engineering<br>Geotechnical Institute<br>Technische Universität Bergakademie Freiberg.</p>\n",
    "        <p><a href=\"https://tu-freiberg.de/en/soilmechanics\">https://tu-freiberg.de/en/soilmechanics</a></p>\n",
    "    </div>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Active Earth Pressure Coulomb\n",
    "\n",
    "![Erddruck.png](https://raw.githubusercontent.com/nagelt/Teaching_Scripts/master/Images/Erddruck.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Geometry \n",
    "Let the retaining wall have a vertical height $H$ and the origin be placed at the bottom of the slip surface which has length $L$. Then, the problem is characterized by the following points in the $x-z$ plane:\n",
    "\n",
    "* Origin: $O = [0,0]$\n",
    "* Top left: $A = [-H \\tan \\alpha, H]$\n",
    "* Top right: $B = [L\\cos \\vartheta, L \\sin \\vartheta]$\n",
    "\n",
    "The length of the slip surface can be found via the following trigonometric relation:\n",
    "\n",
    "\\begin{align}\n",
    "    L \\sin \\vartheta &= (L \\cos \\vartheta + H \\tan\\alpha) \\tan \\beta\n",
    "\\end{align}\n",
    "\n",
    "Thus we arrive at\n",
    "\n",
    "\\begin{align}\n",
    "    L &= \\frac{H(1 + \\tan \\alpha \\tan \\beta)}{\\sin \\vartheta - \\cos \\vartheta \\tan \\beta}\n",
    "\\end{align}\n",
    "\n",
    "The area of the triangle (i.e. the volume of the slip body per meter) can also be calculated. We use Gauß' equation for that:\n",
    "\n",
    "$$\n",
    "    A = \\frac{1}{2} \\left| (x_2 y_3 - x_3 y_2) - (x_1 y_3 - x_3 y_1) + (x_1 y_2 - x_2 y_1) \\right|\n",
    "$$\n",
    "\n",
    "Let Python do the math for us ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "#HIDDEN\n",
    "import numpy as np #numerical methods\n",
    "#import sympy as sp #symbolic operations\n",
    "import matplotlib.pyplot as plt #plotting\n",
    "import sympy as sp\n",
    "sp.init_printing(use_latex='mathjax') #makes sympy output look nice\n",
    "\n",
    "#Some plot settings\n",
    "import plot_functions.plot_settings\n",
    "%run plot_functions/earth_pressure_plots.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first define our geometrical properties:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "H,alpha,beta,theta = sp.symbols('H alpha beta vartheta')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{H \\left(\\tan{\\left(\\alpha \\right)} \\tan{\\left(\\beta \\right)} + 1\\right)}{\\sin{\\left(\\vartheta \\right)} - \\cos{\\left(\\vartheta \\right)} \\tan{\\left(\\beta \\right)}}$"
      ],
      "text/plain": [
       "       H⋅(tan(α)⋅tan(β) + 1)        \n",
       "────────────────────────────────────\n",
       "sin(vartheta) - cos(vartheta)⋅tan(β)"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L = H*(1+sp.tan(alpha)*sp.tan(beta))/(sp.sin(theta) - sp.cos(theta) * sp.tan(beta))\n",
    "L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now define the function that calculates the area of the triangle:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "def triangle_area_Gauss(A,B,C):\n",
    "    term1 = (B[0]*C[1] - C[0] * B[1])\n",
    "    term2 = (A[0]*C[1] - C[0] * A[1])\n",
    "    term3 = (A[0]*B[1] - B[0] * A[1])\n",
    "    return -(term1 - term2 + term3)/2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The points we defined above can now be inserted into this equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{H^{2} \\left(\\sin{\\left(\\vartheta \\right)} \\tan{\\left(\\alpha \\right)} + \\cos{\\left(\\vartheta \\right)}\\right) \\left(\\tan{\\left(\\alpha \\right)} \\tan{\\left(\\beta \\right)} + 1\\right)}{2 \\left(\\sin{\\left(\\vartheta \\right)} - \\cos{\\left(\\vartheta \\right)} \\tan{\\left(\\beta \\right)}\\right)}$"
      ],
      "text/plain": [
       " 2                                                           \n",
       "H ⋅(sin(vartheta)⋅tan(α) + cos(vartheta))⋅(tan(α)⋅tan(β) + 1)\n",
       "─────────────────────────────────────────────────────────────\n",
       "          2⋅(sin(vartheta) - cos(vartheta)⋅tan(β))           "
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PO = [0,0]\n",
    "PA = [- H * sp.tan(alpha), H]\n",
    "PB = [L * sp.cos(theta), L*sp.sin(theta)]\n",
    "A = triangle_area_Gauss(PO,PA,PB).simplify()\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus, the area of our triangle is given as:\n",
    "\n",
    "$$\n",
    "    A = \\frac{H^2 (\\sin \\vartheta \\tan \\alpha + \\cos \\vartheta) (\\tan \\alpha \\tan \\beta + 1)}{2(\\sin \\vartheta - \\cos \\vartheta \\tan \\beta)}\n",
    "$$\n",
    "\n",
    "Now we can write the force balance. For that purpose, we introduce the material's specific weight $\\gamma$, friction angle $\\varphi$ and the wall friction angle $\\delta_\\text{a}$ as well as the cohesion $c$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "c, delta, phi, gamma = sp.symbols('c delta_a varphi gamma')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now calculate the downward force due to gravity as well as the cohesive force along the slip surface. The active earth pressure force $E_\\text{a}$ and the reaction force $Q$ in the slip line are introduced as unknowns:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ea, Q, G, C = sp.symbols('E_a Q G C')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now write the vertical and horizontal equilibrium conditions:\n",
    "\n",
    "\\begin{align}\n",
    "    \\text{horizontal:}\\quad  0 &= C \\cos\\vartheta - Q \\sin (\\vartheta - \\varphi) + E_\\text{a} \\cos (\\alpha + \\delta_\\text{a})\n",
    "    \\\\\n",
    "    \\text{vertikal:}\\quad  0 &= C \\sin\\vartheta + Q \\cos (\\vartheta - \\varphi) + E_\\text{a} \\sin (\\alpha + \\delta_\\text{a}) - G\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fh = sp.Eq(C * sp.cos(theta) - Q * sp.sin(theta - phi) + Ea * sp.cos(alpha + delta),0)\n",
    "Fv = sp.Eq(C * sp.sin(theta) + Q * sp.cos(theta - phi) + Ea * sp.sin(alpha + delta) - G,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle C \\cos{\\left(\\vartheta \\right)} + E_{a} \\cos{\\left(\\alpha + \\delta_{a} \\right)} + Q \\sin{\\left(\\varphi - \\vartheta \\right)} = 0$"
      ],
      "text/plain": [
       "C⋅cos(vartheta) + Eₐ⋅cos(α + δₐ) + Q⋅sin(varphi - vartheta) = 0"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Fh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle C \\sin{\\left(\\vartheta \\right)} + E_{a} \\sin{\\left(\\alpha + \\delta_{a} \\right)} - G + Q \\cos{\\left(\\varphi - \\vartheta \\right)} = 0$"
      ],
      "text/plain": [
       "C⋅sin(vartheta) + Eₐ⋅sin(α + δₐ) - G + Q⋅cos(varphi - vartheta) = 0"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Fv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now solve the equation system symbolically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = sp.solve([Fh,Fv],(Ea,Q))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left\\{ E_{a} : \\frac{C \\sin{\\left(\\vartheta \\right)} \\sin{\\left(\\varphi - \\vartheta \\right)}}{- \\sin{\\left(\\alpha + \\delta_{a} \\right)} \\sin{\\left(\\varphi - \\vartheta \\right)} + \\cos{\\left(\\alpha + \\delta_{a} \\right)} \\cos{\\left(\\varphi - \\vartheta \\right)}} - \\frac{C \\cos{\\left(\\vartheta \\right)} \\cos{\\left(\\varphi - \\vartheta \\right)}}{- \\sin{\\left(\\alpha + \\delta_{a} \\right)} \\sin{\\left(\\varphi - \\vartheta \\right)} + \\cos{\\left(\\alpha + \\delta_{a} \\right)} \\cos{\\left(\\varphi - \\vartheta \\right)}} - \\frac{G \\sin{\\left(\\varphi - \\vartheta \\right)}}{- \\sin{\\left(\\alpha + \\delta_{a} \\right)} \\sin{\\left(\\varphi - \\vartheta \\right)} + \\cos{\\left(\\alpha + \\delta_{a} \\right)} \\cos{\\left(\\varphi - \\vartheta \\right)}}, \\  Q : - \\frac{C \\sin{\\left(\\vartheta \\right)} \\cos{\\left(\\alpha + \\delta_{a} \\right)}}{- \\sin{\\left(\\alpha + \\delta_{a} \\right)} \\sin{\\left(\\varphi - \\vartheta \\right)} + \\cos{\\left(\\alpha + \\delta_{a} \\right)} \\cos{\\left(\\varphi - \\vartheta \\right)}} + \\frac{C \\sin{\\left(\\alpha + \\delta_{a} \\right)} \\cos{\\left(\\vartheta \\right)}}{- \\sin{\\left(\\alpha + \\delta_{a} \\right)} \\sin{\\left(\\varphi - \\vartheta \\right)} + \\cos{\\left(\\alpha + \\delta_{a} \\right)} \\cos{\\left(\\varphi - \\vartheta \\right)}} + \\frac{G \\cos{\\left(\\alpha + \\delta_{a} \\right)}}{- \\sin{\\left(\\alpha + \\delta_{a} \\right)} \\sin{\\left(\\varphi - \\vartheta \\right)} + \\cos{\\left(\\alpha + \\delta_{a} \\right)} \\cos{\\left(\\varphi - \\vartheta \\right)}}\\right\\}$"
      ],
      "text/plain": [
       "⎧                     C⋅sin(vartheta)⋅sin(varphi - vartheta)                   ↪\n",
       "⎨Eₐ: ────────────────────────────────────────────────────────────────────────  ↪\n",
       "⎩    -sin(α + δₐ)⋅sin(varphi - vartheta) + cos(α + δₐ)⋅cos(varphi - vartheta)  ↪\n",
       "\n",
       "↪                    C⋅cos(vartheta)⋅cos(varphi - vartheta)                    ↪\n",
       "↪ - ──────────────────────────────────────────────────────────────────────── - ↪\n",
       "↪   -sin(α + δₐ)⋅sin(varphi - vartheta) + cos(α + δₐ)⋅cos(varphi - vartheta)   ↪\n",
       "\n",
       "↪                          G⋅sin(varphi - vartheta)                            ↪\n",
       "↪  ────────────────────────────────────────────────────────────────────────, Q ↪\n",
       "↪  -sin(α + δₐ)⋅sin(varphi - vartheta) + cos(α + δₐ)⋅cos(varphi - vartheta)    ↪\n",
       "\n",
       "↪                           C⋅sin(vartheta)⋅cos(α + δₐ)                        ↪\n",
       "↪ : - ──────────────────────────────────────────────────────────────────────── ↪\n",
       "↪     -sin(α + δₐ)⋅sin(varphi - vartheta) + cos(α + δₐ)⋅cos(varphi - vartheta) ↪\n",
       "\n",
       "↪                          C⋅sin(α + δₐ)⋅cos(vartheta)                         ↪\n",
       "↪  + ────────────────────────────────────────────────────────────────────────  ↪\n",
       "↪    -sin(α + δₐ)⋅sin(varphi - vartheta) + cos(α + δₐ)⋅cos(varphi - vartheta)  ↪\n",
       "\n",
       "↪                                G⋅cos(α + δₐ)                              ⎫\n",
       "↪ + ────────────────────────────────────────────────────────────────────────⎬\n",
       "↪   -sin(α + δₐ)⋅sin(varphi - vartheta) + cos(α + δₐ)⋅cos(varphi - vartheta)⎭"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{C \\cos{\\left(\\varphi \\right)} + G \\sin{\\left(\\varphi - \\vartheta \\right)}}{\\cos{\\left(\\alpha + \\delta_{a} + \\varphi - \\vartheta \\right)}}$"
      ],
      "text/plain": [
       "-(C⋅cos(varphi) + G⋅sin(varphi - vartheta)) \n",
       "────────────────────────────────────────────\n",
       "      cos(α + δₐ + varphi - vartheta)       "
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result.get(Ea).simplify()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{C \\sin{\\left(\\alpha + \\delta_{a} - \\vartheta \\right)} + G \\cos{\\left(\\alpha + \\delta_{a} \\right)}}{\\cos{\\left(\\alpha + \\delta_{a} + \\varphi - \\vartheta \\right)}}$"
      ],
      "text/plain": [
       "C⋅sin(α + δₐ - vartheta) + G⋅cos(α + δₐ)\n",
       "────────────────────────────────────────\n",
       "    cos(α + δₐ + varphi - vartheta)     "
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result.get(Q).simplify()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gravitational and cohesive forces can be determined based on material properties and geometry (here $c=0$):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{H^{2} \\gamma \\left(\\sin{\\left(\\vartheta \\right)} \\tan{\\left(\\alpha \\right)} + \\cos{\\left(\\vartheta \\right)}\\right) \\left(\\tan{\\left(\\alpha \\right)} \\tan{\\left(\\beta \\right)} + 1\\right) \\sin{\\left(\\varphi - \\vartheta \\right)}}{2 \\left(\\sin{\\left(\\vartheta \\right)} - \\cos{\\left(\\vartheta \\right)} \\tan{\\left(\\beta \\right)}\\right) \\cos{\\left(\\alpha + \\delta_{a} + \\varphi - \\vartheta \\right)}}$"
      ],
      "text/plain": [
       "  2                                                                            ↪\n",
       "-H ⋅γ⋅(sin(vartheta)⋅tan(α) + cos(vartheta))⋅(tan(α)⋅tan(β) + 1)⋅sin(varphi -  ↪\n",
       "────────────────────────────────────────────────────────────────────────────── ↪\n",
       "        2⋅(sin(vartheta) - cos(vartheta)⋅tan(β))⋅cos(α + δₐ + varphi - varthet ↪\n",
       "\n",
       "↪           \n",
       "↪ vartheta) \n",
       "↪ ──────────\n",
       "↪ a)        "
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result.get(Ea).subs(C,0).subs(G,A*gamma).simplify()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we would have to find the maximum with respect to $\\vartheta$ in order to identify the active earth pressure and the critical slip angle $\\vartheta_\\text{a}$. The equations get very complex and we decide to go for a numerical solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "dd9cceaae0404a2eb3a4232f1cb0140d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=10, description='$c$ / kN m$^{-2}$:', max=50, step=5), IntSlider(value=3…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_interactive_earth_pressure()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tasks\n",
    "\n",
    "* Compare the numbers you get to results obtained from the analytical expressions you find in the current norm.\n",
    "* Investigate the influence of the different parameters.\n",
    "* What happens at very high cohesion values? What does this mean?\n",
    "* How high can a cohesive vertical slope be roughly? Select suitable parameters to answer this question.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Erddruck.png](https://raw.githubusercontent.com/nagelt/Teaching_Scripts/master/Images/Erddruck.png)"
   ]
  }
 ],
 "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.12.7"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}