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

396 lines
41 KiB
Plaintext
Raw Normal View History

{
"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 PyPlot\n",
"using ForwardDiff\n",
"using NLsolve"
]
},
{
"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": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"calculate_stress (generic function with 1 method)"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\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",
" e1 = (σ[1] - σ[2])^2\n",
" e2 = (σ[2] - σ[3])^2\n",
" e3 = (σ[3] - σ[1])^2\n",
" e4 = σ[4]^2 \n",
" e5 = σ[5]^2\n",
" e6 = σ[6]^2\n",
" return sqrt((e1 + e2 + e3 + 6 * (e4 + e5 + e6)) / 2.)\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 functions where 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, Δt, σ_begin)\n",
" # Creating wrapper for gradient\n",
" yield(pars) = vonMisesYield(pars, k)\n",
" dfdσ = ForwardDiff.gradient(yield)\n",
" \n",
" # Stress rate\n",
" dσ = (σ_begin - params[1:6]) / Δt\n",
" \n",
" # Calculating plastic strain rate\n",
" dϵp = params[end] * dfdσ(params)\n",
" \n",
" # Calculating equations\n",
" function_1 = dσ - C * (dϵ - dϵp[1:6])\n",
" function_2 = yield(params)\n",
" [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ϵ, Δt, σ, C, k)\n",
"\n",
" # Test stress\n",
" σ_tria = σ + Δt * C * dϵ\n",
" \n",
" # Calculating yield\n",
" yield = vonMisesYield(σ_tria, k)\n",
" \n",
" if yield > 1\n",
" \n",
" # Yielding happened\n",
" # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values\n",
" initial_guess = [σ[:]; 0.0]\n",
" f(σ_) = G(σ_, dϵ, C, k, Δt, σ_tria)\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",
" dλ = result[end]\n",
" \n",
" # these are not necessary to calculate for our calculation but why not?\n",
" yield_f(σ_) = vonMisesYield(σ_, k)\n",
" dfdσ_ = ForwardDiff.gradient(yield_f)\n",
" dϵᵖ = dλ * dfdσ_(σ)\n",
" else\n",
" dϵᵖ = zeros(Float64, (6, 1))\n",
" σ = σ_tria\n",
" end\n",
" return (dϵᵖ, σ)\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Defining strain history"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"linspace(-0.0,-0.0006,200)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"steps = 200\n",
"\n",
"ϵ_tot = zeros(Float64, (steps, 6))\n",
"\n",
"# Adding only strain in x-axis and counting for the poisson effect\n",
"ϵ_tot[:, 1] = linspace(0, 0.002, steps)\n",
"ϵ_tot[:, 2] = linspace(0, 0.002, steps).*-ν\n",
"ϵ_tot[:, 3] = linspace(0, 0.002, 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": 11,
"metadata": {
"collapsed": false,
"scrolled": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAtcAAAI6CAYAAADsVGQBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xt8jvXjx/H3fbNTDrPN+ayvQyjNKSH54os55TRzzlkh852kkKSJpRS1kqKcySnHWA6LkcMUKodSiS+mzJxtxnb9/rhz/9xtmLm36972ej4e9yO77s997X3d3a23a5/rc1kMwzAEAAAA4IFZzQ4AAAAAZBeUawAAAMBJKNcAAACAk1CuAQAAACehXAMAAABOQrkGAAAAnIRyDQAAADgJ5RoAAABwEso1AAAA4CSUawBAmpUtW1blypUzOwYAuCzKNQBTJCUl6dNPP1XDhg3l6+srd3d3FSlSRI8//rgGDBigNWvWOIyfPXu2rFar5syZY1LijJeQkKB33nlHderUkbe3tzw8PFS8eHHVqlVLQ4cO1bZt2xzGv/7667JarSm2ZySLxSKLxZJp3w8AsprcZgcAkPMkJSWpdevWioiIkI+Pj1q3bq2SJUsqMTFRP/30kxYuXKiff/5Zbdq0SfHa7Frsrly5ooYNG2rfvn0qVqyYOnXqpKJFi+rKlSvav3+/PvnkE128eFFPP/20qTm3bNli6vcHAFdHuQaQ6RYtWqSIiAj5+/tr69atypcvn8Pz8fHx2rNnT6qvNQwjMyJmuqlTp2rfvn1q3ry51qxZo9y5HX88X7hwQUeOHEn1tZn5njAlBADujmkhADLdt99+K0nq3bt3imItSV5eXmrYsKH963//+9/q27evJKlPnz6yWq32x4kTJyT9/xSJrVu3auHChapTp47y5s3rUAavXbumSZMmyd/fX3nz5lW+fPlUr149LV68ONWcc+bMUb169VSoUCF5eXmpdOnSCggI0JIlSxzG/fDDD+ratavKli0rT09PFS5cWDVr1lRISIhu3rx5X+/JoEGDUhRrSSpQoICefPJJ+9dly5bVG2+8IUlq1KiRw3tyS+/evWW1WnXs2DF98MEHqlatmh566CE1atRIknTjxg2Fh4erZcuWKlOmjDw9PeXn56emTZtqw4YNqeZMbc717VN2IiMj9e9//1v58+eXt7e3Wrdufce/FNzN119/rTZt2qhw4cLy9PRU6dKl1a5dO23evDnV75saq9VqP9Zb7vY52b17t6xWqzp06HDHXJUrV5anp6cuXLjgsD0iIkItW7ZUwYIF5enpqfLly2vkyJG6ePHifR87gKyNM9cAMl3BggUlST///HOaxvfp00c+Pj5atWqV2rVrJ39/f/tz3t7eDmOnTJmijRs36plnnlGTJk3s5ebChQtq3Lix9u/fr5o1a6pfv35KTk7Whg0b1K1bNx08eFChoaH2/YwePVphYWF6+OGH1aVLF3l7e+v06dOKjo7WsmXLFBQUJMlWrOvUqaNcuXLpmWeeUbly5XTp0iUdPXpU06dP15tvvplqWX7Q9yQkJEQrV67U1q1b1bt3b5UtW/aOY4cNG6aoqCi1bt1arVu3Vq5cuSRJ586d03//+1/Vr19fzZs3V6FChXT69GmtWbNGLVu21Keffqp+/fql2N+dpuasXbtWq1atUsuWLTVo0CAdPHhQX331laKjo3Xo0CH5+fml6djGjRun0NBQ5cuXT+3atVOpUqV06tQpffvtt1qwYIGaNGmSpjx3ey61z0mdOnVUqVIlffXVV4qLi5Ovr6/Da/bs2aOff/5ZgYGBKlCggH37+PHjNX78ePn5+dn/QnDgwAG98847+uqrr7Rz585U/xIJIJsyACCT7du3z3B3dzesVqvRs2dPY8WKFcYff/xx19d8/vnnhsViMebMmZPq8+PGjTMsFouRN29eY//+/Sme79Wrl2GxWIy3337bYXtCQoIREBBgWK1Wh9f5+voapUqVMuLj41PsKzY21v7n4cOHGxaLxVi9enWKcRcuXDCSk5Pvely3rF271rBYLIaHh4cxePBgY926dcbp06fv+ppbx7x169ZUn791zCVLlkz1/b1+/bpx6tSpFNsvXrxoPProo4avr2+K4y9TpoxRrlw5h223/t24ubkZW7ZscXhu1KhRhsViMSZPnnzXY7klIiLCsFgsxr/+9a9Uj//kyZMpvu+dPhMWi8Vo1KiRw7Z7fU4mTZpkWCwWIzw8PMVzgwcPNiwWi7F27Vr7ti1bthgWi8WoX7++cfHiRYfxs2fPNiwWixESEnL3gwaQrTAtBECm8/f31/z581WkSBHNnz9fHTt2VLly5eTn56cOHTpo7dq16d73wIED9fjjjztsO3funObPn6/atWtrxIgRDs95eHgoLCxMhmFo4cKF9u0Wi0Vubm4O0yxuSe0MrKenZ4pt3t7eab4As1WrVpo2bZq8vLw0ffp0tW7dWiVKlFCxYsXUo0cPRUVFpWk/qRk5cqTKlCmTYru7u7uKFy+eYnv+/PnVp08fnT9/XtHR0Wn+Pl26dEkxDWPgwIGSlOb9fPDBB5JsZ5aLFSuW4vkSJUqkOc/dpPY5kaSePXumOtUkMTFRixcvVpEiRdSiRQv79vfff1+S9Omnnyp//vwOr+nVq5cef/xxLViwwCmZAWQNTAsBYIpOnTqpffv2ioyM1I4dO7Rv3z5t375dK1eu1MqVK/Xss89q9uzZ973fJ554IsW26OhoJScnS7LNuf2nGzduSJIOHz5s39a9e3d98MEHqlKlioKCgtSwYUM9+eSTKaahdOnSRe+//77atWunwMBANWnSRPXr19e//vUvh3H79+/XypUrHbb5+Pho2LBh9q+HDh2q/v37a+PGjdq5c6f27dunb7/9VgsXLtTChQs1duxYjR8//v7eEKX+ntxy8OBBvf3229q2bZvOnDmjhIQEh+dPnz6d5u9Tq1atFNtKliwpSTp//nya9rFr1y5ZrVYFBASk+fumx53ekxIlSqhJkybauHGjDh8+rMqVK0uS1qxZo/Pnz2v48OEOf+HauXOn3NzctGTJklQvLE1MTNTZs2d1/vx5+fj4ZMzBAHAplGsApsmdO7eaNm2qpk2bSpKSk5O1fPly9e3bV3PnzlX79u3Vtm3b+9pn0aJFU2w7d+6cJFvJvtMZVIvFoqtXr9q/fu+99/Twww/r888/V1hYmMLCwpQ7d261bNlSU6ZMsZfn2rVrKyoqSm+++aaWLVumefPmSZIqVaqkcePGqUuXLpKkAwcO2C9AvKVs2bIO5VqyXcz5zDPP6JlnnpFkK/6ffvqphg0bptDQUHXo0CHVM673+55ItiLbuHFjJScnq0mTJmrXrp3y588vq9Wqffv2adWqVbp+/Xqav8/t85BvuTXfPCkpKU37uHDhgnx8fOTh4ZHm75sed3pPJNuFoBs3btScOXMUFhYmSfYz2b169XIYe+7cOSUlJd31Lz0Wi0VXrlyhXAM5BNNCALgMq9WqTp06KSQkRJIUGRl53/tIbRrGrbPNw4cPV3JycqqPpKQkh5UorFarhg0bpv379+vPP//U8uXL1b59e61evVoBAQFKTEy0j33yySe1Zs0aXbhwQTt27NDYsWP1559/qlu3bvZ99urVK8X3/P333+95PG5ubho8eLC6du0qKX3rTN9pasqECROUkJCgr7/+WuvWrdO7776r119/Xa+99tpdz3ZnpAIFCuj8+fMpzqCn5tYZ5NRWZPnnah7/dLfpOu3bt1f+/Pk1f/58GYahv/76S+vXr5e/v78ee+wxh7He3t7y9fW94+fq1merVKlS9zweANkD5RqAy8mbN68kx/Wbb61wkdYzoLerU6fOA93JsFChQmrfvr2++OILNWrUSL/99psOHjyYYpybm5vq1q2r8ePH2+firl69Ol3f859uvSe3e5D3RJJ+/fVX+fn5pXpjmq1bt6Zrnw+qbt269lVc7uXWmeBbyzHebu/evenO4OnpqaCgIJ0+fVobN27UwoULlZSUlOKs9a28cXFxOnToULq/H4DshXININMtWrRImzZtSnWO6pkzZ/Tpp59KkkPpu3UR4fHjx+/7+xUqVEjdu3fX3r17NWHCBPv869v99ttv+uOPPyTZ5snu2LE
"text/plain": [
"PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x7f31a0ba76d0>)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"(0,250)"
]
},
"execution_count": 11,
"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 = zeros(Float64, steps)\n",
"ee = zeros(Float64, steps)\n",
"for i=1:steps\n",
" ss[i] = σ[1]\n",
" ee[i] = ϵ_last[1]\n",
" dϵ = (reshape(ϵ_tot[i, :, :], (6, 1)) - ϵ_last) / Δt \n",
" dϵᵖ, σ = calculate_stress(dϵ, Δt, σ, C, σy)\n",
" ϵ_last += dϵ * Δt\n",
"end\n",
"PyPlot.plot(ee, ss)\n",
"PyPlot.title(\"Stress-Strain curve\")\n",
"PyPlot.xlabel(\"Strain\")\n",
"PyPlot.ylabel(\"Stress\")\n",
"PyPlot.ylim([0, 250])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.5.0-dev",
"language": "julia",
"name": "julia-0.5"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "0.5.0"
}
},
"nbformat": 4,
"nbformat_minor": 0
}