Slope_Stability_Control.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e2c03424-604d-4629-bfba-e3f2b785fe86",
   "metadata": {},
   "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": "code",
   "execution_count": 37,
   "id": "f6aef413-a2bd-48c2-896e-eccfcfe5f1bc",
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6066ade8-8b08-4513-ae33-ce735a79b8c3",
   "metadata": {},
   "source": [
    "## Slope failure from a material point stability perspective\n",
    "\n",
    "*Based on: Buscarnera and Di Prisco (2011)*\n",
    "\n",
    "**Disclaimer: the following is not finished or checked, yet. Handle with care.**\n",
    "\n",
    "We depart from the assumption of an infinite slope and are interested in the stress and strain evolution in a plane parallel to the slope surface (simple shear kinematics).\n",
    "\n",
    "Furthermore, we assume the following set of equations to hold:\n",
    "\n",
    "\\begin{align}\n",
    "    \\dot{\\sigma}_\\text{n}' &= E_\\text{s} \\dot{\\epsilon}_\\text{el}\n",
    "    \\\\\n",
    "    \\dot{\\tau} &= G \\dot{\\gamma}_\\text{el}\n",
    "    \\\\\n",
    "    f &= \\tau - \\sigma_\\text{n}' \\tan \\varphi - c\n",
    "    \\\\\n",
    "    g &= \\tau - \\sigma_\\text{n}' \\tan \\psi\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf1ae0ed-1489-4613-b582-771242c641c8",
   "metadata": {},
   "source": [
    "The flow rule yields\n",
    "\n",
    "\\begin{align}\n",
    "    \\dot{\\epsilon}_\\text{pl} &= - \\lambda \\tan \\psi\n",
    "    \\\\\n",
    "    \\dot{\\gamma}_\\text{pl} &= \\lambda\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7703053b-d2d5-4ce5-9add-bae79d69b5d1",
   "metadata": {},
   "source": [
    "Consistency in the plastic state requires\n",
    "\n",
    "\\begin{align}\n",
    "    \\dot{f} = 0 = \\dot{\\tau} - \\dot{\\sigma}_\\text{n}' \\tan \\varphi\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9eab82d0-b777-4b94-9741-9ad96bfe6807",
   "metadata": {},
   "source": [
    "Inserting the stress evolution\n",
    "\n",
    "\\begin{align}\n",
    "    \\dot{\\sigma}_\\text{n}' &= E_\\text{s} \\left( \\dot{\\epsilon} + \\lambda \\tan \\psi \\right)\n",
    "    \\\\\n",
    "    \\dot{\\tau} &= G \\left( \\dot{\\gamma} - \\lambda \\right)\n",
    "\\end{align}\n",
    "\n",
    "into the consistency condition yields\n",
    "\n",
    "\\begin{align}\n",
    "    0 &= G \\left( \\dot{\\gamma} - \\lambda \\right) - E_\\text{s} \\left( \\dot{\\epsilon} + \\lambda \\tan \\psi \\right) \\tan \\varphi\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "449ef64f-2a5d-4788-b758-a6bde264545e",
   "metadata": {},
   "source": [
    "$$\n",
    "    \\lambda = \\frac{G \\dot{\\gamma} - E_\\text{s} \\tan\\varphi \\dot{\\epsilon}}{G + E_\\text{s} \\tan\\varphi \\tan\\psi}\n",
    "$$\n",
    "\n",
    "### From governing equations, over material stability to controllability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "089d901b-64bd-4737-b695-951a160ab66a",
   "metadata": {},
   "outputs": [],
   "source": [
    "lam, sign, tau, gamma, eps = sp.symbols('lambda \\dot{\\sigma}_\\mathrm{n} \\dot{\\\\tau} \\dot{\\gamma} \\dot{\\epsilon}')\n",
    "E, G, phi, psi = sp.symbols('E_\\mathrm{s} G varphi psi', positive=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "5c5c4dbb-29ae-4d86-8aaf-a059dddca27b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[ E_\\mathrm{s} \\left(\\dot{\\epsilon} + \\lambda \\tan{\\left(\\psi \\right)}\\right) = \\dot{\\sigma}_\\mathrm{n}, \\  G \\left(\\dot{\\gamma} - \\lambda\\right) = \\dot{\\tau}, \\  - E_\\mathrm{s} \\left(\\dot{\\epsilon} + \\lambda \\tan{\\left(\\psi \\right)}\\right) \\tan{\\left(\\varphi \\right)} + G \\left(\\dot{\\gamma} - \\lambda\\right) = 0\\right]$"
      ],
      "text/plain": [
       "[E_\\mathrm{s}⋅(\\dot{\\epsilon} + λ⋅tan(ψ)) = \\dot{\\sigma}_\\mathrm{n}, G⋅(\\dot{\\\n",
       "gamma} - λ) = \\dot{\\tau}, -E_\\mathrm{s}⋅(\\dot{\\epsilon} + λ⋅tan(ψ))⋅tan(varphi\n",
       ") + G⋅(\\dot{\\gamma} - λ) = 0]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Fs = sp.Eq(E * (eps + lam * sp.tan(psi)), sign)\n",
    "Ft = sp.Eq(G*(gamma - lam),tau)\n",
    "Fl = sp.Eq(G*(gamma - lam) - E*(eps+lam*sp.tan(psi))*sp.tan(phi),0)\n",
    "display([Fs,Ft,Fl])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "1857452f-6b19-48f6-a8f8-0582ddc9ab5b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left\\{ \\dot{\\sigma}_\\mathrm{n} : \\frac{E_\\mathrm{s} G \\dot{\\epsilon} + E_\\mathrm{s} G \\dot{\\gamma} \\tan{\\left(\\psi \\right)}}{E_\\mathrm{s} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G}, \\  \\dot{\\tau} : \\frac{E_\\mathrm{s} G \\dot{\\epsilon} \\tan{\\left(\\varphi \\right)} + E_\\mathrm{s} G \\dot{\\gamma} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)}}{E_\\mathrm{s} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G}, \\  \\lambda : \\frac{- E_\\mathrm{s} \\dot{\\epsilon} \\tan{\\left(\\varphi \\right)} + G \\dot{\\gamma}}{E_\\mathrm{s} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G}\\right\\}$"
      ],
      "text/plain": [
       "⎧                         E_\\mathrm{s}⋅G⋅\\dot{\\epsilon} + E_\\mathrm{s}⋅G⋅\\dot{\n",
       "⎨\\dot{\\sigma}_\\mathrm{n}: ────────────────────────────────────────────────────\n",
       "⎩                                        E_\\mathrm{s}⋅tan(ψ)⋅tan(varphi) + G  \n",
       "\n",
       "\\gamma}⋅tan(ψ)              E_\\mathrm{s}⋅G⋅\\dot{\\epsilon}⋅tan(varphi) + E_\\mat\n",
       "──────────────, \\dot{\\tau}: ──────────────────────────────────────────────────\n",
       "                                                       E_\\mathrm{s}⋅tan(ψ)⋅tan\n",
       "\n",
       "hrm{s}⋅G⋅\\dot{\\gamma}⋅tan(ψ)⋅tan(varphi)     -E_\\mathrm{s}⋅\\dot{\\epsilon}⋅tan(\n",
       "────────────────────────────────────────, λ: ─────────────────────────────────\n",
       "(varphi) + G                                            E_\\mathrm{s}⋅tan(ψ)⋅ta\n",
       "\n",
       "varphi) + G⋅\\dot{\\gamma}⎫\n",
       "────────────────────────⎬\n",
       "n(varphi) + G           ⎭"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result = sp.solve([Fs,Ft,Fl],(sign,tau,lam))\n",
    "sp.simplify(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96393f3a-773a-40d9-866b-3d2d14bfab5f",
   "metadata": {},
   "source": [
    "With this, the stress-strain relations in the elasto-plastic state reads\n",
    "\n",
    "$$\n",
    "    \\begin{pmatrix}\n",
    "    \\dot{\\sigma}_\\text{n}' \\\\ \\dot{\\tau}\n",
    "    \\end{pmatrix}\n",
    "    =\n",
    "    \\frac{E_\\text{s} G}{E_\\text{s}\\tan\\psi\\tan\\varphi + G} \\begin{pmatrix}\n",
    "    1 & \\tan\\psi\n",
    "    \\\\\n",
    "    \\tan\\varphi & \\tan\\varphi\\tan\\psi\n",
    "    \\end{pmatrix}\n",
    "    \\begin{pmatrix}\n",
    "    \\dot{\\epsilon} \\\\ \\dot{\\gamma}\n",
    "    \\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aefbe1a8-1b68-4ec1-be38-88fbb0bc92a2",
   "metadata": {},
   "source": [
    "The system matrix is the control matrix for a drained slope."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "3b1cdfbb-ca69-461d-b18d-86caa1b13b1e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{E_\\mathrm{s} G}{E_\\mathrm{s} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G} & \\frac{E_\\mathrm{s} G \\tan{\\left(\\psi \\right)}}{E_\\mathrm{s} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G}\\\\\\frac{E_\\mathrm{s} G \\tan{\\left(\\varphi \\right)}}{E_\\mathrm{s} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G} & \\frac{E_\\mathrm{s} G \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)}}{E_\\mathrm{s} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡           E_\\mathrm{s}⋅G                   E_\\mathrm{s}⋅G⋅tan(ψ)       ⎤\n",
       "⎢───────────────────────────────────  ───────────────────────────────────⎥\n",
       "⎢E_\\mathrm{s}⋅tan(ψ)⋅tan(varphi) + G  E_\\mathrm{s}⋅tan(ψ)⋅tan(varphi) + G⎥\n",
       "⎢                                                                        ⎥\n",
       "⎢     E_\\mathrm{s}⋅G⋅tan(varphi)       E_\\mathrm{s}⋅G⋅tan(ψ)⋅tan(varphi) ⎥\n",
       "⎢───────────────────────────────────  ───────────────────────────────────⎥\n",
       "⎣E_\\mathrm{s}⋅tan(ψ)⋅tan(varphi) + G  E_\\mathrm{s}⋅tan(ψ)⋅tan(varphi) + G⎦"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D_drained = sp.Matrix([[1,sp.tan(psi)],[sp.tan(phi),sp.tan(phi)*sp.tan(psi)]]) * E*G/(E*sp.tan(psi)*sp.tan(phi)+G)\n",
    "D_drained"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "4ab8089c-2bc7-4e48-947c-df39b6c9bdae",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.det(D_drained)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7443e70-28af-45ba-b47d-20400bc29a4a",
   "metadata": {},
   "source": [
    "The vanishing determinant indicates a loss of stability independent of the chosen material parameters and the tendency for the MC model to localize (elastic ideal plastic).\n",
    "\n",
    "For the undrained case, we assume isochoric deformation and thus $\\dot{\\epsilon} = 0$ becomes a controlling variable in the present case. The system matrix is found by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "a2178725-faee-42bb-aa40-0bdccf3ff968",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left\\{ \\dot{\\epsilon} : \\frac{- E_\\mathrm{s} G \\dot{\\gamma} \\tan{\\left(\\psi \\right)} + E_\\mathrm{s} \\dot{\\sigma}_\\mathrm{n} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G \\dot{\\sigma}_\\mathrm{n}}{E_\\mathrm{s} G}, \\  \\dot{\\tau} : \\dot{\\sigma}_\\mathrm{n} \\tan{\\left(\\varphi \\right)}, \\  \\lambda : \\frac{G \\dot{\\gamma} - \\dot{\\sigma}_\\mathrm{n} \\tan{\\left(\\varphi \\right)}}{G}\\right\\}$"
      ],
      "text/plain": [
       "⎧                -E_\\mathrm{s}⋅G⋅\\dot{\\gamma}⋅tan(ψ) + E_\\mathrm{s}⋅\\dot{\\sigm\n",
       "⎨\\dot{\\epsilon}: ─────────────────────────────────────────────────────────────\n",
       "⎩                                                                      E_\\math\n",
       "\n",
       "a}_\\mathrm{n}⋅tan(ψ)⋅tan(varphi) + G⋅\\dot{\\sigma}_\\mathrm{n}                  \n",
       "────────────────────────────────────────────────────────────, \\dot{\\tau}: \\dot\n",
       "rm{s}⋅G                                                                       \n",
       "\n",
       "                                    G⋅\\dot{\\gamma} - \\dot{\\sigma}_\\mathrm{n}⋅t\n",
       "{\\sigma}_\\mathrm{n}⋅tan(varphi), λ: ──────────────────────────────────────────\n",
       "                                                             G                \n",
       "\n",
       "an(varphi)⎫\n",
       "──────────⎬\n",
       "          ⎭"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result = sp.solve([Fs,Ft,Fl],(eps,tau,lam))\n",
    "sp.simplify(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f581bb48-dcca-4755-84d2-3ade8d001250",
   "metadata": {},
   "source": [
    "$$\n",
    "    \\begin{pmatrix}\n",
    "    \\dot{\\epsilon} \\\\ \\dot{\\tau}\n",
    "    \\end{pmatrix}\n",
    "    =\n",
    "    \\begin{pmatrix}\n",
    "    \\frac{G + E_\\text{s}\\tan\\psi\\tan{\\varphi}}{E_\\text{s}G} & -\\tan\\psi\n",
    "    \\\\\n",
    "    \\tan\\varphi & 0\n",
    "    \\end{pmatrix}\n",
    "    \\begin{pmatrix}\n",
    "    \\dot{\\sigma}_\\text{n}' \\\\ \\dot{\\gamma}\n",
    "    \\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "ece68a77-91da-40e8-b42e-43be1cbe8caa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{E_\\mathrm{s} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G}{E_\\mathrm{s} G} & - \\tan{\\left(\\psi \\right)}\\\\\\tan{\\left(\\varphi \\right)} & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡E_\\mathrm{s}⋅tan(ψ)⋅tan(varphi) + G         ⎤\n",
       "⎢───────────────────────────────────  -tan(ψ)⎥\n",
       "⎢           E_\\mathrm{s}⋅G                   ⎥\n",
       "⎢                                            ⎥\n",
       "⎣            tan(varphi)                 0   ⎦"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D_undrained = sp.Matrix([[(G+E*sp.tan(phi)*sp.tan(psi))/(E*G),-sp.tan(psi)],[sp.tan(phi),0]])\n",
    "D_undrained"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "376a2846-c93b-431c-9b21-91989070b4c7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)}$"
      ],
      "text/plain": [
       "tan(ψ)⋅tan(varphi)"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.det(D_undrained)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77231e6a-bfe3-4570-aaf9-5e3a2b92e3f9",
   "metadata": {},
   "source": [
    "This determinant vanishes for vanishing friction or dilatancy angles. The plausible solution in our case for a vanishing (or negative) second-order work is thus the transition to negative dilatancy angles, i.e. contractile material behaviour. For positive dilatancy the condition indicates stability.\n",
    "\n",
    "Note in passing: this determinant corresponds to $D_{22} = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc5b43cd-7394-421d-bf11-cb21e1dfb565",
   "metadata": {},
   "source": [
    "### Eigendecomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fc44aab-6a35-4a60-87a1-32ab7f68c6d2",
   "metadata": {},
   "source": [
    "Note: function returns eigenvalue, multiplicity and eigenvector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "25d0ba7c-e27a-4cb5-8e24-6f5ec6e187d3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[ \\left( 0, \\  1, \\  \\left[ \\left[\\begin{matrix}- \\tan{\\left(\\psi \\right)}\\\\1\\end{matrix}\\right]\\right]\\right), \\  \\left( \\frac{E_\\mathrm{s} G \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + E_\\mathrm{s} G}{E_\\mathrm{s} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G}, \\  1, \\  \\left[ \\left[\\begin{matrix}\\frac{1}{\\tan{\\left(\\varphi \\right)}}\\\\1\\end{matrix}\\right]\\right]\\right)\\right]$"
      ],
      "text/plain": [
       "⎡                     ⎛                                                       \n",
       "⎢⎛      ⎡⎡-tan(ψ)⎤⎤⎞  ⎜E_\\mathrm{s}⋅G⋅tan(ψ)⋅tan(varphi) + E_\\mathrm{s}⋅G     \n",
       "⎢⎜0, 1, ⎢⎢       ⎥⎥⎟, ⎜──────────────────────────────────────────────────, 1, \n",
       "⎢⎝      ⎣⎣   1   ⎦⎦⎠  ⎜       E_\\mathrm{s}⋅tan(ψ)⋅tan(varphi) + G             \n",
       "⎣                     ⎝                                                       \n",
       "\n",
       "⎡⎡     1     ⎤⎤⎞⎤\n",
       "⎢⎢───────────⎥⎥⎟⎥\n",
       "⎢⎢tan(varphi)⎥⎥⎟⎥\n",
       "⎢⎢           ⎥⎥⎟⎥\n",
       "⎣⎣     1     ⎦⎦⎠⎦"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D_drained.eigenvects()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3969f19-a5d4-4a3a-bdf6-f9b8b83e6bfc",
   "metadata": {},
   "source": [
    "Vanishing eigenstiffness associated with the dilatant shear mode."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "8c705714-3de0-4ed1-bc81-5c075baf6c44",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[ \\left( \\frac{E_\\mathrm{s} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G}{2 E_\\mathrm{s} G} - \\frac{\\sqrt{- 4 E_\\mathrm{s}^{2} G^{2} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + E_\\mathrm{s}^{2} \\tan^{2}{\\left(\\psi \\right)} \\tan^{2}{\\left(\\varphi \\right)} + 2 E_\\mathrm{s} G \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G^{2}}}{2 E_\\mathrm{s} G}, \\  1, \\  \\left[ \\left[\\begin{matrix}\\frac{\\frac{E_\\mathrm{s} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G}{2 E_\\mathrm{s} G} - \\frac{\\sqrt{- 4 E_\\mathrm{s}^{2} G^{2} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + E_\\mathrm{s}^{2} \\tan^{2}{\\left(\\psi \\right)} \\tan^{2}{\\left(\\varphi \\right)} + 2 E_\\mathrm{s} G \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G^{2}}}{2 E_\\mathrm{s} G}}{\\tan{\\left(\\varphi \\right)}}\\\\1\\end{matrix}\\right]\\right]\\right), \\  \\left( \\frac{E_\\mathrm{s} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G}{2 E_\\mathrm{s} G} + \\frac{\\sqrt{- 4 E_\\mathrm{s}^{2} G^{2} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + E_\\mathrm{s}^{2} \\tan^{2}{\\left(\\psi \\right)} \\tan^{2}{\\left(\\varphi \\right)} + 2 E_\\mathrm{s} G \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G^{2}}}{2 E_\\mathrm{s} G}, \\  1, \\  \\left[ \\left[\\begin{matrix}\\frac{\\frac{E_\\mathrm{s} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G}{2 E_\\mathrm{s} G} + \\frac{\\sqrt{- 4 E_\\mathrm{s}^{2} G^{2} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + E_\\mathrm{s}^{2} \\tan^{2}{\\left(\\psi \\right)} \\tan^{2}{\\left(\\varphi \\right)} + 2 E_\\mathrm{s} G \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G^{2}}}{2 E_\\mathrm{s} G}}{\\tan{\\left(\\varphi \\right)}}\\\\1\\end{matrix}\\right]\\right]\\right)\\right]$"
      ],
      "text/plain": [
       "⎡⎛                                                                            \n",
       "⎢⎜                                         ___________________________________\n",
       "⎢⎜                                        ╱                 2  2              \n",
       "⎢⎜E_\\mathrm{s}⋅tan(ψ)⋅tan(varphi) + G   ╲╱  - 4⋅E_\\mathrm{s} ⋅G ⋅tan(ψ)⋅tan(va\n",
       "⎢⎜─────────────────────────────────── - ──────────────────────────────────────\n",
       "⎢⎜          2⋅E_\\mathrm{s}⋅G                                                  \n",
       "⎢⎜                                                                            \n",
       "⎢⎜                                                                            \n",
       "⎣⎝                                                                            \n",
       "\n",
       "                                                                              \n",
       "______________________________________________________________________________\n",
       "                    2    2       2                                            \n",
       "rphi) + E_\\mathrm{s} ⋅tan (ψ)⋅tan (varphi) + 2⋅E_\\mathrm{s}⋅G⋅tan(ψ)⋅tan(varph\n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                2⋅E_\\mathrm{s}⋅G                                              \n",
       "                                                                              \n",
       "                                                                              \n",
       "                                                                              \n",
       "\n",
       "             ⎡⎡                                         ______________________\n",
       "________     ⎢⎢                                        ╱                 2  2 \n",
       "      2      ⎢⎢E_\\mathrm{s}⋅tan(ψ)⋅tan(varphi) + G   ╲╱  - 4⋅E_\\mathrm{s} ⋅G ⋅\n",
       "i) + G       ⎢⎢─────────────────────────────────── - ─────────────────────────\n",
       "────────, 1, ⎢⎢          2⋅E_\\mathrm{s}⋅G                                     \n",
       "             ⎢⎢───────────────────────────────────────────────────────────────\n",
       "             ⎢⎢                                                               \n",
       "             ⎢⎢                                                               \n",
       "             ⎣⎣                                                               \n",
       "\n",
       "______________________________________________________________________________\n",
       "                                 2    2       2                               \n",
       "tan(ψ)⋅tan(varphi) + E_\\mathrm{s} ⋅tan (ψ)⋅tan (varphi) + 2⋅E_\\mathrm{s}⋅G⋅tan\n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                             2⋅E_\\mathrm{s}⋅G                                 \n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "            tan(varphi)                                                       \n",
       "                                                                              \n",
       "                 1                                                            \n",
       "\n",
       "_____________________⎤⎤⎞  ⎛                                                   \n",
       "                   2 ⎥⎥⎟  ⎜                                         __________\n",
       "(ψ)⋅tan(varphi) + G  ⎥⎥⎟  ⎜                                        ╱          \n",
       "─────────────────────⎥⎥⎟  ⎜E_\\mathrm{s}⋅tan(ψ)⋅tan(varphi) + G   ╲╱  - 4⋅E_\\ma\n",
       "                     ⎥⎥⎟, ⎜─────────────────────────────────── + ─────────────\n",
       "─────────────────────⎥⎥⎟  ⎜          2⋅E_\\mathrm{s}⋅G                         \n",
       "                     ⎥⎥⎟  ⎜                                                   \n",
       "                     ⎥⎥⎟  ⎜                                                   \n",
       "                     ⎦⎦⎠  ⎝                                                   \n",
       "\n",
       "                                                                              \n",
       "______________________________________________________________________________\n",
       "       2  2                                  2    2       2                   \n",
       "thrm{s} ⋅G ⋅tan(ψ)⋅tan(varphi) + E_\\mathrm{s} ⋅tan (ψ)⋅tan (varphi) + 2⋅E_\\mat\n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                         2⋅E_\\mathrm{s}⋅G                     \n",
       "                                                                              \n",
       "                                                                              \n",
       "                                                                              \n",
       "\n",
       "                                      ⎡⎡                                      \n",
       "_________________________________     ⎢⎢                                      \n",
       "                               2      ⎢⎢E_\\mathrm{s}⋅tan(ψ)⋅tan(varphi) + G   \n",
       "hrm{s}⋅G⋅tan(ψ)⋅tan(varphi) + G       ⎢⎢─────────────────────────────────── + \n",
       "─────────────────────────────────, 1, ⎢⎢          2⋅E_\\mathrm{s}⋅G            \n",
       "                                      ⎢⎢──────────────────────────────────────\n",
       "                                      ⎢⎢                                      \n",
       "                                      ⎢⎢                                      \n",
       "                                      ⎣⎣                                      \n",
       "\n",
       "   ___________________________________________________________________________\n",
       "  ╱                 2  2                                  2    2       2      \n",
       "╲╱  - 4⋅E_\\mathrm{s} ⋅G ⋅tan(ψ)⋅tan(varphi) + E_\\mathrm{s} ⋅tan (ψ)⋅tan (varph\n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                                      2⋅E_\\mathrm{s}⋅G        \n",
       "──────────────────────────────────────────────────────────────────────────────\n",
       "                                     tan(varphi)                              \n",
       "                                                                              \n",
       "                                          1                                   \n",
       "\n",
       "______________________________________________⎤⎤⎞⎤\n",
       "                                            2 ⎥⎥⎟⎥\n",
       "i) + 2⋅E_\\mathrm{s}⋅G⋅tan(ψ)⋅tan(varphi) + G  ⎥⎥⎟⎥\n",
       "──────────────────────────────────────────────⎥⎥⎟⎥\n",
       "                                              ⎥⎥⎟⎥\n",
       "──────────────────────────────────────────────⎥⎥⎟⎥\n",
       "                                              ⎥⎥⎟⎥\n",
       "                                              ⎥⎥⎟⎥\n",
       "                                              ⎦⎦⎠⎦"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D_undrained.eigenvects()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7524f0c4-2d02-4eb8-9688-6452a4074dd9",
   "metadata": {},
   "source": [
    "Instead of the general result, let's look at the previously identified point of instability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "cc255993-6df5-4b6f-9077-b312d2930242",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[ \\left( 0, \\  1, \\  \\left[ \\left[\\begin{matrix}0\\\\1\\end{matrix}\\right]\\right]\\right), \\  \\left( \\frac{1}{E_\\mathrm{s}}, \\  1, \\  \\left[ \\left[\\begin{matrix}\\frac{1}{E_\\mathrm{s} \\tan{\\left(\\varphi \\right)}}\\\\1\\end{matrix}\\right]\\right]\\right)\\right]$"
      ],
      "text/plain": [
       "⎡               ⎛                 ⎡⎡           1            ⎤⎤⎞⎤\n",
       "⎢⎛      ⎡⎡0⎤⎤⎞  ⎜     1           ⎢⎢────────────────────────⎥⎥⎟⎥\n",
       "⎢⎜0, 1, ⎢⎢ ⎥⎥⎟, ⎜────────────, 1, ⎢⎢E_\\mathrm{s}⋅tan(varphi)⎥⎥⎟⎥\n",
       "⎢⎝      ⎣⎣1⎦⎦⎠  ⎜E_\\mathrm{s}     ⎢⎢                        ⎥⎥⎟⎥\n",
       "⎣               ⎝                 ⎣⎣           1            ⎦⎦⎠⎦"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D_undrained.subs(psi,0).eigenvects()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a309bf96-b119-4f23-bdf6-59e1591de65a",
   "metadata": {},
   "source": [
    "Eigenmode corresponds to constant normal stress shear deformation.... "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "cf13bfae-0bd7-4dba-b62a-3e090d396e10",
   "metadata": {},
   "outputs": [],
   "source": [
    "D_drained_s = (D_drained + D_drained.T)/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "6fbad504-8e78-422e-9bf3-109ff64a7133",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & \\frac{E_\\mathrm{s} G \\left(\\tan{\\left(\\psi \\right)} - \\tan{\\left(\\varphi \\right)}\\right)}{2 \\left(E_\\mathrm{s} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G\\right)}\\\\\\frac{E_\\mathrm{s} G \\left(- \\tan{\\left(\\psi \\right)} + \\tan{\\left(\\varphi \\right)}\\right)}{2 \\left(E_\\mathrm{s} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G\\right)} & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡                                          E_\\mathrm{s}⋅G⋅(tan(ψ) - tan(varphi\n",
       "⎢                   0                     ────────────────────────────────────\n",
       "⎢                                         2⋅(E_\\mathrm{s}⋅tan(ψ)⋅tan(varphi) +\n",
       "⎢                                                                             \n",
       "⎢ E_\\mathrm{s}⋅G⋅(-tan(ψ) + tan(varphi))                                      \n",
       "⎢───────────────────────────────────────                     0                \n",
       "⎣2⋅(E_\\mathrm{s}⋅tan(ψ)⋅tan(varphi) + G)                                      \n",
       "\n",
       ")) ⎤\n",
       "───⎥\n",
       " G)⎥\n",
       "   ⎥\n",
       "   ⎥\n",
       "   ⎥\n",
       "   ⎦"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.simplify(D_drained - D_drained_s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "0eec1573-2a3f-4865-9fc7-fd6b4e3e68e5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{E_\\mathrm{s}^{2} G^{2} \\left(\\tan{\\left(\\psi \\right)} - \\tan{\\left(\\varphi \\right)}\\right)^{2}}{4 \\left(E_\\mathrm{s} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G\\right)^{2}}$"
      ],
      "text/plain": [
       "             2  2                       2 \n",
       "-E_\\mathrm{s} ⋅G ⋅(tan(ψ) - tan(varphi))  \n",
       "──────────────────────────────────────────\n",
       "                                        2 \n",
       " 4⋅(E_\\mathrm{s}⋅tan(ψ)⋅tan(varphi) + G)  "
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "detDs = sp.det(D_drained_s).simplify()\n",
    "detDs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "2ac24c10-3078-4900-8d5c-ba15ac45f6d7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{E_\\mathrm{s}^{2} G^{2} \\left(\\tan{\\left(\\psi \\right)} - \\tan{\\left(\\varphi \\right)}\\right)^{2}}{4 \\left(E_\\mathrm{s} \\tan{\\left(\\psi \\right)} \\tan{\\left(\\varphi \\right)} + G\\right)^{2}} = 0$"
      ],
      "text/plain": [
       "             2  2                       2     \n",
       "-E_\\mathrm{s} ⋅G ⋅(tan(ψ) - tan(varphi))      \n",
       "────────────────────────────────────────── = 0\n",
       "                                        2     \n",
       " 4⋅(E_\\mathrm{s}⋅tan(ψ)⋅tan(varphi) + G)      "
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Fds = sp.Eq(detDs,0)\n",
    "Fds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "id": "100db6c2-4799-4861-9b38-03c30055b5d2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\operatorname{atan}{\\left(\\tan{\\left(\\varphi \\right)} \\right)}$"
      ],
      "text/plain": [
       "atan(tan(varphi))"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.solve(Fds,psi)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9daf8df4-47f5-4d6b-82f0-2945145c3e74",
   "metadata": {},
   "source": [
    "Possible stability issues due to non-associativity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66f5317e-2322-48d2-b63b-a5f4e7de8961",
   "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.9.13"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}