Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-08-06 04:21:33 +00:00
Code Issues Packages Projects Releases Wiki Activity
Files
12c0f60af9ca4a576eca56e2e49570a62adca68b
JuliaFEM.jl/notebooks/2015-09-24-Ideal plastic Von Mises material.ipynb
T

396 lines
41 KiB
Plaintext
Raw Normal View History

Added notebook for ideal plastic Von Mises material
2015-09-24 15:03:02 +03:00
{
"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",
Update 2015-09-24-Ideal plastic Von Mises material.ipynb
2015-09-24 15:19:35 +03:00
"Function for NLsolve. Inside this function are the equations which we want to find root.\n",
Added notebook for ideal plastic Von Mises material
2015-09-24 15:03:02 +03:00
"Ψ 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
}
Reference in New Issue Copy Permalink