mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
Olli Väinölä
333 lines
8.7 KiB
Plaintext
333 lines
8.7 KiB
Plaintext
{
|
||
"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": null,
|
||
"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": null,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [],
|
||
"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": null,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [],
|
||
"source": [
|
||
"\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, Δt, σ_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ϵ, Δt, σ, C, k)\n",
|
||
" # Test stress\n",
|
||
" σ_tria = σ + Δt * 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, Δt, σ)\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",
|
||
" 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": null,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [],
|
||
"source": [
|
||
"steps = 10\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": null,
|
||
"metadata": {
|
||
"collapsed": false,
|
||
"scrolled": false
|
||
},
|
||
"outputs": [],
|
||
"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",
|
||
" dϵ = (reshape(ϵ_tot[i, :, :], (6, 1)) - ϵ_last) / Δt \n",
|
||
" dϵᵖ, σ = calculate_stress(dϵ, Δt, σ, C, σy)\n",
|
||
" ϵ_last += dϵ * Δt\n",
|
||
" ss[i] = σ[1]\n",
|
||
" ee[i] = ϵ_last[1]\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
|
||
}
|