passive_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": [
        "# Passive Earth Pressure Coulomb\n",
        "\n",
        "![Erddruck_passiv.png](https://raw.githubusercontent.com/nagelt/Teaching_Scripts/master/Images/Erddruck_passiv.png)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This is an educational example only. The assumption of planar slip surfaces leads to an overestimation of passive earth pressures. Curved slip surfaces should be used instead."
      ]
    },
    {
      "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": 1,
      "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": 2,
      "metadata": {},
      "outputs": [],
      "source": [
        "H,alpha,beta,theta = sp.symbols('H alpha beta vartheta')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "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": 3,
          "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": 4,
      "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": 5,
      "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": 5,
          "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": 6,
      "metadata": {},
      "outputs": [],
      "source": [
        "c, delta, phi, gamma = sp.symbols('c delta_p 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 passive earth pressure force $E_\\text{p}$ and the reaction force $Q$ in the slip line are introduced as unknowns:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {},
      "outputs": [],
      "source": [
        "Ep, Q, G, C = sp.symbols('E_p 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{p} \\cos (\\alpha - \\delta_\\text{p})\n",
        "    \\\n",
        "    \\text{vertikal:}\\quad  0 &= -C \\sin\\vartheta + Q \\cos (\\vartheta + \\varphi) + E_\\text{p} \\sin (\\alpha - \\delta_\\text{p}) - G\n",
        "\\end{align}"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {},
      "outputs": [],
      "source": [
        "Fh = sp.Eq(-C * sp.cos(theta) - Q * sp.sin(theta + phi) + Ep * sp.cos(alpha - delta),0)\n",
        "Fv = sp.Eq(-C * sp.sin(theta) + Q * sp.cos(theta + phi) + Ep * sp.sin(alpha - delta) - G,0)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/latex": [
              "$\\displaystyle - C \\cos{\\left(\\vartheta \\right)} + E_{p} \\cos{\\left(\\alpha - \\delta_{p} \\right)} - Q \\sin{\\left(\\varphi + \\vartheta \\right)} = 0$"
            ],
            "text/plain": [
              "-C⋅cos(vartheta) + Eₚ⋅cos(α - δₚ) - Q⋅sin(varphi + vartheta) = 0"
            ]
          },
          "execution_count": 9,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "Fh"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/latex": [
              "$\\displaystyle - C \\sin{\\left(\\vartheta \\right)} + E_{p} \\sin{\\left(\\alpha - \\delta_{p} \\right)} - G + Q \\cos{\\left(\\varphi + \\vartheta \\right)} = 0$"
            ],
            "text/plain": [
              "-C⋅sin(vartheta) + Eₚ⋅sin(α - δₚ) - G + Q⋅cos(varphi + vartheta) = 0"
            ]
          },
          "execution_count": 10,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "Fv"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's find a symbolic solution."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {},
      "outputs": [],
      "source": [
        "result = sp.solve([Fh,Fv],(Ep,Q))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/latex": [
              "$\\displaystyle \\left\\{ E_{p} : \\frac{C \\sin{\\left(\\vartheta \\right)} \\sin{\\left(\\varphi + \\vartheta \\right)}}{\\sin{\\left(\\alpha - \\delta_{p} \\right)} \\sin{\\left(\\varphi + \\vartheta \\right)} + \\cos{\\left(\\alpha - \\delta_{p} \\right)} \\cos{\\left(\\varphi + \\vartheta \\right)}} + \\frac{C \\cos{\\left(\\vartheta \\right)} \\cos{\\left(\\varphi + \\vartheta \\right)}}{\\sin{\\left(\\alpha - \\delta_{p} \\right)} \\sin{\\left(\\varphi + \\vartheta \\right)} + \\cos{\\left(\\alpha - \\delta_{p} \\right)} \\cos{\\left(\\varphi + \\vartheta \\right)}} + \\frac{G \\sin{\\left(\\varphi + \\vartheta \\right)}}{\\sin{\\left(\\alpha - \\delta_{p} \\right)} \\sin{\\left(\\varphi + \\vartheta \\right)} + \\cos{\\left(\\alpha - \\delta_{p} \\right)} \\cos{\\left(\\varphi + \\vartheta \\right)}}, \\  Q : \\frac{C \\sin{\\left(\\vartheta \\right)} \\cos{\\left(\\alpha - \\delta_{p} \\right)}}{\\sin{\\left(\\alpha - \\delta_{p} \\right)} \\sin{\\left(\\varphi + \\vartheta \\right)} + \\cos{\\left(\\alpha - \\delta_{p} \\right)} \\cos{\\left(\\varphi + \\vartheta \\right)}} - \\frac{C \\sin{\\left(\\alpha - \\delta_{p} \\right)} \\cos{\\left(\\vartheta \\right)}}{\\sin{\\left(\\alpha - \\delta_{p} \\right)} \\sin{\\left(\\varphi + \\vartheta \\right)} + \\cos{\\left(\\alpha - \\delta_{p} \\right)} \\cos{\\left(\\varphi + \\vartheta \\right)}} + \\frac{G \\cos{\\left(\\alpha - \\delta_{p} \\right)}}{\\sin{\\left(\\alpha - \\delta_{p} \\right)} \\sin{\\left(\\varphi + \\vartheta \\right)} + \\cos{\\left(\\alpha - \\delta_{p} \\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)   sin\n",
              "\n",
              "                     G⋅sin(varphi + vartheta)                                 \n",
              "────────────────────────────────────────────────────────────────────, Q: ─────\n",
              "(α - δₚ)⋅sin(varphi + vartheta) + cos(α - δₚ)⋅cos(varphi + vartheta)     sin(α\n",
              "\n",
              "                 C⋅sin(vartheta)⋅cos(α - δₚ)                                  \n",
              "────────────────────────────────────────────────────────────────── - ─────────\n",
              " - δₚ)⋅sin(varphi + vartheta) + cos(α - δₚ)⋅cos(varphi + vartheta)   sin(α - δ\n",
              "\n",
              "             C⋅sin(α - δₚ)⋅cos(vartheta)                                      \n",
              "────────────────────────────────────────────────────────────── + ─────────────\n",
              "ₚ)⋅sin(varphi + vartheta) + cos(α - δₚ)⋅cos(varphi + vartheta)   sin(α - δₚ)⋅s\n",
              "\n",
              "                G⋅cos(α - δₚ)                             ⎫\n",
              "──────────────────────────────────────────────────────────⎬\n",
              "in(varphi + vartheta) + cos(α - δₚ)⋅cos(varphi + vartheta)⎭"
            ]
          },
          "execution_count": 12,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "result"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/latex": [
              "$\\displaystyle \\frac{C \\cos{\\left(\\varphi \\right)} + G \\sin{\\left(\\varphi + \\vartheta \\right)}}{\\cos{\\left(- \\alpha + \\delta_{p} + \\varphi + \\vartheta \\right)}}$"
            ],
            "text/plain": [
              "C⋅cos(varphi) + G⋅sin(varphi + vartheta)\n",
              "────────────────────────────────────────\n",
              "    cos(-α + δₚ + varphi + vartheta)    "
            ]
          },
          "execution_count": 13,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "result.get(Ep).simplify()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/latex": [
              "$\\displaystyle \\frac{C \\sin{\\left(- \\alpha + \\delta_{p} + \\vartheta \\right)} + G \\cos{\\left(\\alpha - \\delta_{p} \\right)}}{\\cos{\\left(- \\alpha + \\delta_{p} + \\varphi + \\vartheta \\right)}}$"
            ],
            "text/plain": [
              "C⋅sin(-α + δₚ + vartheta) + G⋅cos(α - δₚ)\n",
              "─────────────────────────────────────────\n",
              "     cos(-α + δₚ + varphi + vartheta)    "
            ]
          },
          "execution_count": 14,
          "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:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "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_{p} + \\varphi + \\vartheta \\right)}}$"
            ],
            "text/plain": [
              " 2                                                                            \n",
              "H ⋅γ⋅(sin(vartheta)⋅tan(α) + cos(vartheta))⋅(tan(α)⋅tan(β) + 1)⋅sin(varphi + v\n",
              "──────────────────────────────────────────────────────────────────────────────\n",
              "      2⋅(sin(vartheta) - cos(vartheta)⋅tan(β))⋅cos(-α + δₚ + varphi + vartheta\n",
              "\n",
              "        \n",
              "artheta)\n",
              "────────\n",
              ")       "
            ]
          },
          "execution_count": 15,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "result.get(Ep).subs(C,0).subs(G,A*gamma).simplify()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we would have to find the minimum with respect to $\\vartheta$ in order to identify the passive earth pressure and the critical slip angle $\\vartheta_\\text{p}$. The equations get very complex and we decide to go for a numerical solution."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "metadata": {},
      "outputs": [
        {
          "data": {
            "application/vnd.jupyter.widget-view+json": {
              "model_id": "93689403b7da4f18b0533dd3d4f1ab62",
              "version_major": 2,
              "version_minor": 0
            },
            "text/plain": [
              "interactive(children=(IntSlider(value=10, description=r'$c$ / kN m$^{-2}$:', max=50, step=5), IntSlider(value=3…"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "plot_interactive_earth_pressure_passive()"
      ]
    },
    {
      "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."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "![Erddruck_passiv.png](https://raw.githubusercontent.com/nagelt/Teaching_Scripts/master/Images/Erddruck_passiv.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.9.13"
    },
    "widgets": {
      "application/vnd.jupyter.widget-state+json": {
        "state": {},
        "version_major": 2,
        "version_minor": 0
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}