JuliaFEM.jl/notebooks/2015-09-24-Ideal plastic Von Mises material.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ideal plastic Von Mises material\n",
    "\n",
    "Author(s): Olli Väinölä <olli.vainola@student.oulu.fi>\n",
    "\n",
    "In this notebook is an small tutorial, how to create a Von Mises material without any hardening. Equations are formulated into rate depended form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# imports\n",
    "using ForwardDiff\n",
    "using NLsolve\n",
    "using PyPlot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create a isotropic Hooke material."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6x6 Array{Float64,2}:\n",
       " 2.69231e5  1.15385e5  1.15385e5  0.0        0.0        0.0      \n",
       " 1.15385e5  2.69231e5  1.15385e5  0.0        0.0        0.0      \n",
       " 1.15385e5  1.15385e5  2.69231e5  0.0        0.0        0.0      \n",
       " 0.0        0.0        0.0        1.53846e5  0.0        0.0      \n",
       " 0.0        0.0        0.0        0.0        1.53846e5  0.0      \n",
       " 0.0        0.0        0.0        0.0        0.0        1.53846e5"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "\n",
    "\"\"\"\n",
    "Create a isotropic Hooke material matrix C \n",
    "\n",
    "More information: # http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_isotropic.cfm\n",
    "\n",
    "Parameters\n",
    "----------\n",
    "    E: Float\n",
    "        Elastic modulus\n",
    "    ν: Float\n",
    "        Poisson constant\n",
    "\n",
    "Returns\n",
    "-------\n",
    "    Array{Float64, (6,6)}\n",
    "\"\"\"\n",
    "function hookeStiffnessTensor(E, ν)\n",
    "    a = 1 - ν\n",
    "    b = 1 - 2*ν\n",
    "    c = 1 + ν\n",
    "    multiplier = E / (b * c)\n",
    "    return Float64[a ν ν 0 0 0;\n",
    "                   ν a ν 0 0 0;\n",
    "                   ν ν a 0 0 0;\n",
    "                   0 0 0 b 0 0;\n",
    "                   0 0 0 0 b 0;\n",
    "                   0 0 0 0 0 b].*multiplier\n",
    "end\n",
    "\n",
    "# Pick material values\n",
    "E = 200.0e3\n",
    "ν =  0.3\n",
    "C = hookeStiffnessTensor(E, ν)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "# Defining equations for the calculation\n",
    "\n",
    "Functions are defined for strain controller simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "calculate_stress (generic function with 2 methods)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "# using vectors with double contradiction\n",
    "# http://www-2.unipv.it/compmech/teaching/available/const_mod/const_mod_mat-review_notation.pdf\n",
    "M = [1 0 0 0 0 0;\n",
    "     0 1 0 0 0 0;\n",
    "     0 0 1 0 0 0;\n",
    "     0 0 0 2 0 0;\n",
    "     0 0 0 0 2 0;\n",
    "     0 0 0 0 0 2;]\n",
    "\n",
    "\"\"\"\n",
    "Equivalent tensile stress. \n",
    "\n",
    "More info can be found from: https://en.wikipedia.org/wiki/Von_Mises_yield_criterion\n",
    "    Section: Reduced von Mises equation for different stress conditions\n",
    "\n",
    "Parameters\n",
    "----------\n",
    "    σ: Array{Float64, 6}\n",
    "        Stress in Voigt notation\n",
    "\n",
    "Returns\n",
    "-------\n",
    "    Float\n",
    "\"\"\"\n",
    "function σₑ(σ)\n",
    "    s = σ[1:6] - 1/3 * sum([σ[1], σ[2], σ[3]]) * [1 1 1 0 0 0]'\n",
    "    return sqrt(3/2 * s' * M * s)[1]\n",
    "end\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "Von Mises Yield criterion\n",
    "\n",
    "More info can be found from: http://csm.mech.utah.edu/content/wp-content/uploads/2011/10/9tutorialOnJ2Plasticity.pdf\n",
    "\n",
    "Parameters\n",
    "----------\n",
    "    σ: Array{Float64, 6}\n",
    "        Stress in Voigt notation\n",
    "    k: Float64\n",
    "        Material constant, Yield limit\n",
    "\n",
    "Returns\n",
    "-------\n",
    "    Float\n",
    "\"\"\"\n",
    "function vonMisesYield(σ, k)\n",
    "    σₑ(σ) - k\n",
    "end\n",
    "\n",
    "\"\"\"\n",
    "Function for NLsolve. Inside this function are the equations which we want to find root.\n",
    "Ψ is the yield function below. Functions defined here:\n",
    "\n",
    "    dσ - C (dϵ - dλ*dΨ/dσ) = 0\n",
    "                 σₑ(σ) - k = 0\n",
    "\n",
    "Parameters\n",
    "----------\n",
    "    params: Array{Float64, 7}\n",
    "        Array containing values from solver\n",
    "    dϵ: Array{Float64, 6}\n",
    "        Strain rate vector in Voigt notation\n",
    "    C: Array{Float64, (6, 6)}\n",
    "        Material tensor\n",
    "    k: Float\n",
    "        Material constant, yield limit\n",
    "    Δt: Float\n",
    "        time increment\n",
    "    σ_begin:Array{Float64, 6}\n",
    "        Stress vector in Voigt notation\n",
    "\n",
    "Returns\n",
    "-------\n",
    "    Array{Float64, 7}, return values for solver\n",
    "\"\"\"\n",
    "function G(params, dϵ, C, k, σ_begin)\n",
    "\n",
    "    # Creating wrapper for gradient\n",
    "    yield(pars) = vonMisesYield(pars, k)\n",
    "    dfdσ = ForwardDiff.gradient(yield)\n",
    "    \n",
    "    # Stress rate\n",
    "    dσ = params[1:6]\n",
    "    \n",
    "    σ_tot = [vec(σ_begin); 0.0] + params\n",
    "    \n",
    "    # Calculating plastic strain rate\n",
    "    dϵp = params[end] * dfdσ(σ_tot)\n",
    "    \n",
    "    # Calculating equations\n",
    "    function_1 = dσ - C * (dϵ - dϵp[1:6])\n",
    "    function_2 = yield(σ_tot)\n",
    "    [vec(function_1); function_2]\n",
    "end\n",
    "\n",
    "\"\"\"\n",
    "Function which calculates the stress. Also handles if any yielding happens\n",
    "\n",
    "Parameters\n",
    "----------\n",
    "    dϵ: Array{Float64, 6}\n",
    "        Strain rate vector in Voigt notation\n",
    "    Δt: Float\n",
    "        time increment\n",
    "    σ: Array{Float64, 6}\n",
    "        Last stress vector in Voigt notation\n",
    "    C: Array{Float64, (6, 6)}\n",
    "        Material tensor\n",
    "    k: Float\n",
    "        Material constant, yield limit\n",
    "\n",
    "Returns\n",
    "-------\n",
    "    Tuple\n",
    "    Plastic strain rate dϵᵖ and new stress vector σ\n",
    "\"\"\"\n",
    "function calculate_stress(dϵ, σ, C, k)\n",
    "    # Test stress\n",
    "    σ_tria = σ + C * dϵ\n",
    "    \n",
    "    # Calculating yield\n",
    "    yield = vonMisesYield(σ_tria, k)\n",
    "\n",
    "    if yield > 1\n",
    "        # Yielding happened\n",
    "        # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values\n",
    "        initial_guess = [vec(σ_tria - σ); 0.1]\n",
    "        f(σ_)  = G(σ_, dϵ, C, k, σ)\n",
    "        df     = ForwardDiff.jacobian(f)\n",
    "        \n",
    "        # Calculating root \n",
    "        result = nlsolve(not_in_place(f, df), initial_guess).zero\n",
    "\n",
    "        σ[:]  += result[1:6]\n",
    "    else\n",
    "        σ = σ_tria\n",
    "    end\n",
    "    return σ\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Defining strain history"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "300-element Array{Float64,1}:\n",
       " -0.0        \n",
       " -5.67002e-5 \n",
       " -0.000113175\n",
       " -0.0001692  \n",
       " -0.000224554\n",
       " -0.000279014\n",
       " -0.000332367\n",
       " -0.000384399\n",
       " -0.000434904\n",
       " -0.00048368 \n",
       " -0.000530536\n",
       " -0.000575283\n",
       " -0.000617745\n",
       "  ⋮          \n",
       "  0.000575283\n",
       "  0.000530536\n",
       "  0.00048368 \n",
       "  0.000434904\n",
       "  0.000384399\n",
       "  0.000332367\n",
       "  0.000279014\n",
       "  0.000224554\n",
       "  0.0001692  \n",
       "  0.000113175\n",
       "  5.67002e-5 \n",
       "  6.61309e-19"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "steps = 300\n",
    "strain_max = 0.003\n",
    "num_cycles = 3\n",
    "\n",
    "ϵ_tot = zeros(Float64, (steps, 6))\n",
    "ϵ_tot2 = zeros(Float64, (steps, 6))\n",
    "ϵ_tot3 = zeros(Float64, (steps, 6))\n",
    "\n",
    "# Adding only strain in x-axis and counting for the poisson effect\n",
    "ϵ_tot[:, 1] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps))\n",
    "ϵ_tot[:, 2] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν\n",
    "ϵ_tot[:, 3] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation\n",
    "\n",
    "Ok, we're good to go! Now we just need to define yield limit and the main loop.\n",
    "\n",
    "This simulation is not time dependent, but since it's already defined in the equations we'll give it value 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAt8AAAI6CAYAAAD/gyT1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XlYVGX7B/DvjLKpgICoIKG+bqnlbm6VW+64K+pbmuTST1xI7bU0zXpdX5XUtMzUcjcXzARMciH3FBXSXBLNNAFTwAWRdeb8/niaGQ4zFiLMOTN8P9fl1TPPOTNzD4N2zz33eR6NJEkSiIiIiIio2GmVDoCIiIiIqKRg8k1EREREZCVMvomIiIiIrITJNxERERGRlTD5JiIiIiKyEibfRERERERWwuSbiIiIiMhKmHwTEREREVkJk28iIiIiIith8k1ERMWqWrVqqF69utJhEBGpApNvIrIZOp0Oq1atQtu2beHp6QlHR0dUqlQJDRs2xKhRoxAeHi47f+3atdBqtVi3bp1CERe/zMxMLFq0CC1atIC7uzucnJzg6+uLZs2aYfz48Th8+LDs/I8++ghardZsvjhpNBpoNBqrPR8RkZqVVjoAIqKC0Ol0CAgIQFRUFDw8PBAQEAA/Pz9kZ2fjl19+webNm/Hrr7+iZ8+eZve118Tv0aNHaNu2LWJjY+Hj44OBAweicuXKePToEeLi4vDll1/iwYMHePXVVxWN8+DBg4o+PxGRmjD5JiKbsGXLFkRFRaFRo0Y4dOgQXF1dZcczMjJw6tQpi/eVJMkaIVrdkiVLEBsbiy5duiA8PBylS8v/Sb9//z4uX75s8b7W/Jmw5YSIyIRtJ0RkE44fPw4AGD58uFniDQAuLi5o27at8Xa7du3w1ltvAQCCgoKg1WqNf27evAnA1IJx6NAhbN68GS1atEC5cuVkyeLjx48xb948NGrUCOXKlYOrqytat26Nb775xmKc69atQ+vWreHt7Q0XFxf4+/uja9eu2LZtm+y8c+fOYciQIahWrRqcnZ1RsWJFNG3aFBMnTkRubu5T/UzGjBljlngDQPny5dGyZUvj7WrVquG///0vAKB9+/ayn4nB8OHDodVqcf36dSxbtgwNGjRAmTJl0L59ewBATk4Oli9fju7du6Nq1apwdnaGl5cXOnXqhL1791qM01LPd96WoOjoaLRr1w5ubm5wd3dHQEDAEz80/J0ffvgBPXv2RMWKFeHs7Ax/f3/06dMHBw4csPi8lmi1WuNrNfi735OTJ09Cq9WiX79+T4yrbt26cHZ2xv3792XzUVFR6N69OypUqABnZ2fUrFkTU6ZMwYMHD576tROR7WDlm4hsQoUKFQAAv/76a4HODwoKgoeHB7777jv06dMHjRo1Mh5zd3eXnRsaGop9+/ahV69e6NixozH5uX//Pjp06IC4uDg0bdoUI0aMgF6vx969e/Hvf/8bFy5cwKxZs4yPM23aNMyfPx//+te/MHjwYLi7uyMxMRExMTHYsWMHAgMDAYjEu0WLFihVqhR69eqF6tWr4+HDh4iPj8eKFSswZ84ci8n0s/5MJk6ciF27duHQoUMYPnw4qlWr9sRzQ0JCcOTIEQQEBCAgIAClSpUCAKSkpOCdd95BmzZt0KVLF3h7eyMxMRHh4eHo3r07Vq1ahREjRpg93pNafyIiIvDdd9+he/fuGDNmDC5cuIA9e/YgJiYGFy9ehJeXV4Fe28yZMzFr1iy4urqiT58+eO6555CQkIDjx49j06ZN6NixY4Hi+btjln5PWrRogTp16mDPnj1ITU2Fp6en7D6nTp3Cr7/+igEDBqB8+fLG+Y8//hgff/wxvLy8jB8Yfv75ZyxatAh79uzBiRMnLH7IJCI7IBER2YDY2FjJ0dFR0mq10tChQ6WdO3dKv//++9/e5+uvv5Y0Go20bt06i8dnzpwpaTQaqVy5clJcXJzZ8TfffFPSaDTSwoULZfOZmZlS165dJa1WK7ufp6en9Nxzz0kZGRlmj5WcnGwcT5o0SdJoNNLu3bvNzrt//76k1+v/9nUZRERESBqNRnJycpKCg4OlyMhIKTEx8W/vY3jNhw4dsnjc8Jr9/Pws/nyzsrKkhIQEs/kHDx5IL7zwguTp6Wn2+qtWrSpVr15dNmd4bxwcHKSDBw/Kjk2dOlXSaDTSggUL/va1GERFRUkajUaqUaOGxdd/69Yts+d90u+ERqOR2rdvL5v7p9+TefPmSRqNRlq+fLnZseDgYEmj0UgRERHGuYMHD0oajUZq06aN9ODBA9n5a9eulTQajTRx4sS/f9FEZLPYdkJENqFRo0bYuHEjKlWqhI0bN6J///6oXr06vLy80K9fP0RERBT6sUePHo2GDRvK5lJSUrBx40Y0b94c7777ruyYk5MT5s+fD0mSsHnzZuO8RqOBg4ODrI3DwFIF19nZ2WzO3d29wBeI9ujRA0uXLoWLiwtWrFiBgIAAVKlSBT4+PnjjjTdw5MiRAj2OJVOmTEHVqlXN5h0dHeHr62s27+bmhqCgINy7dw8xMTEFfp7BgwebtXmMHj0aAAr8OMuWLQMgKtM+Pj5mx6tUqVLgeP6Opd8TABg6dKjFVpbs7Gx88803qFSpErp162ac//TTTwEAq1atgpubm+w+b775Jho2bIhNmzYVScxEpD5sOyEimzFw4ED07dsX0dHROHbsGGJjY3H06FHs2rULu3btwrBhw7B27dqnftyXXnrJbC4mJgZ6vR6A6PnNLycnBwBw6dIl49zrr7+OZcuWoV69eggMDETbtm3RsmVLszaXwYMH49NPP0WfPn0wYMAAdOzYEW3atEGNGjVk58XFxWHXrl2yOQ8PD4SEhBhvjx8/HiNHjsS+fftw4sQJxMbG4vjx49i8eTM2b96MGTNm4OOPP366Hwgs/0wMLly4gIULF+Lw4cO4ffs2MjMzZccTExML/DzNmjUzm/Pz8wMA3Lt3r0CP8dNPP0Gr1aJr164Fft7CeNLPpEqVKujYsSP27duHS5cuoW7dugCA8PBw3Lt3D5MmTZJ9IDtx4gQcHBywbds2ixe+Zmdn4+7du7h37x48PDyK58UQkWKYfBORTSldujQ6deqETp06AQD0ej3CwsLw1ltvYf369ejbty969+79VI9ZuXJls7mUlBQAIgl/UgVWo9EgPT3deHvx4sX417/+ha+//hrz58/H/PnzUbp0aXTv3h2hoaHG5Lp58+Y4cuQI5syZgx07dmDDhg0AgDp16mDmzJkYPHgwAODnn382XiBpUK1aNVnyDYiLTXv16oVevXoBEB8MVq1ahZCQEMyaNQv9+vWzWLF92p8JIBLdDh06QK/Xo2PHjujTpw/c3Nyg1WoRGxuL7777DllZWQV+nrx90AaGfnedTlegx7h//z48PDzg5ORU4OctjCf9TABxoeq+ffuwbt06zJ8/HwCMlfA333xTdm5KSgp0Ot3ffijSaDR49OgRk28iO8S2EyKyaVqtFgMHDsTEiRMBANHR0U/9GJbaPAzV6kmTJkGv11v8o9PpZCtpaLVahISEIC4uDn/++SfCwsLQt29f7N69G127dkV2drbx3JYtWyI8PBz379/HsWPHMGPGDPz555/497//bXzMN9980+w5f/vtt398PQ4ODggODsaQIUMAFG6d7Se1vsyePRuZmZn44YcfEBkZiU8++QQfffQRPvzww7+tlhen8uXL4969e2YVeEsMFWhLK8rkX40kv79rB+rbty/c3NywceNGSJKEO3fu4Pvvv0ejRo3w4osvys51d3eHp6fnE3+vDL9bzz333D++HiKyPUy+icgulCtXDoB8/WrDCh0FraDm1aJFi2faCdLb2xt9+/bF1q1b0b59e1y7dg0XLlwwO8/BwQGtWrXCxx9/bOwF3r17d6GeMz/DzySvZ/mZAMDVq1fh5eVlceOeQ4cOFeoxn1WrVq2Mq9D8E0Ml2bDcZF6nT58udAzOzs4IDAxEYmIi9u3bh82bN0On05lVvQ3xpqam4uLFi4V+PiKyXUy+icgmbNmyBfv377fYI3v79m2sWrUKAGRJoeEixxs3bjz183l7e+P111/H6dOnMXv2bGP/d17Xrl3D77//DkD06R47dszsnJycHKSmpkKj0aBMmTIAxPrclqq
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x7f6b660944d0>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "PyObject <matplotlib.text.Text object at 0x7f6b63ad4f90>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ϵ_last = zeros(Float64, (6))\n",
    "ϵᵖ = zeros(Float64, (6))\n",
    "σ = zeros(Float64, (6, 1))\n",
    "Δt = 1.0\n",
    "σy =  200.0\n",
    "ss = Float64[]\n",
    "ee = Float64[]\n",
    "\n",
    "for i=1:steps\n",
    "    dϵ = reshape(ϵ_tot[i, :, :], (6, 1)) - ϵ_last\n",
    "    σ = calculate_stress(dϵ, σ, C, σy)\n",
    "    ϵ_last += dϵ * Δt\n",
    "    push!(ss, σ[1])\n",
    "    push!(ee, ϵ_last[1])    \n",
    "end\n",
    "\n",
    "PyPlot.plot(ee, ss)\n",
    "PyPlot.title(\"Stress-Strain curve\")\n",
    "PyPlot.xlabel(\"Strain\")\n",
    "PyPlot.ylabel(\"Stress\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.1-pre",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}