From f811f4c48a67a85baeaaed9f4c18afb12c4f5300 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olli=20V=C3=A4in=C3=B6l=C3=A4?= Date: Wed, 9 Dec 2015 15:45:15 +0200 Subject: [PATCH] Added a ideal plastic von Mises with plot --- src/vonMises.jl | 273 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 273 insertions(+) create mode 100644 src/vonMises.jl diff --git a/src/vonMises.jl b/src/vonMises.jl new file mode 100644 index 0000000..c3bfd88 --- /dev/null +++ b/src/vonMises.jl @@ -0,0 +1,273 @@ +# imports +using ForwardDiff +using NLsolve +using PyPlot + +""" +Create a isotropic Hooke material matrix C + +More information: # http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_isotropic.cfm + +Parameters +---------- + E: Float + Elastic modulus + ν: Float + Poisson constant + +Returns +------- + Array{Float64, (6,6)} +""" +function hookeStiffnessTensor(E, ν) + a = 1 - ν + b = 1 - 2*ν + c = 1 + ν + multiplier = E / (b * c) + return Float64[a ν ν 0 0 0; + ν a ν 0 0 0; + ν ν a 0 0 0; + 0 0 0 b 0 0; + 0 0 0 0 b 0; + 0 0 0 0 0 b].*multiplier +end + +# Pick material values +E = 200.0e3 +ν = 0.3 +C = hookeStiffnessTensor(E, ν) + + +type State + C :: Array{Float64, 2} + σ_y :: Float64 + σ :: Array{Float64, 1} + ϵ :: Array{Float64, 1} +end + +# using vectors with double contradiction +# http://www-2.unipv.it/compmech/teaching/available/const_mod/const_mod_mat-review_notation.pdf +M = [1 0 0 0 0 0; + 0 1 0 0 0 0; + 0 0 1 0 0 0; + 0 0 0 2 0 0; + 0 0 0 0 2 0; + 0 0 0 0 0 2;] + +""" +Equivalent tensile stress. + +More info can be found from: https://en.wikipedia.org/wiki/Von_Mises_yield_criterion + Section: Reduced von Mises equation for different stress conditions + +Parameters +---------- + σ: Array{Float64, 6} + Stress in Voigt notation + +Returns +------- + Float +""" +function σₑ(σ) + s = σ[1:6] - 1/3 * sum([σ[1], σ[2], σ[3]]) * [1 1 1 0 0 0]' + return sqrt(3/2 * s' * M * s)[1] +end + + +""" +Von Mises Yield criterion + +More info can be found from: http://csm.mech.utah.edu/content/wp-content/uploads/2011/10/9tutorialOnJ2Plasticity.pdf + +Parameters +---------- + σ: Array{Float64, 6} + Stress in Voigt notation + k: Float64 + Material constant, Yield limit + +Returns +------- + Float +""" +function vonMisesYield(σ, k) + σₑ(σ) - k +end + +""" +Function for NLsolve. Inside this function are the equations which we want to find root. +Ψ is the yield function below. Functions defined here: + + dσ - C (dϵ - dλ*dΨ/dσ) = 0 + σₑ(σ) - k = 0 + +Parameters +---------- + params: Array{Float64, 7} + Array containing values from solver + dϵ: Array{Float64, 6} + Strain rate vector in Voigt notation + C: Array{Float64, (6, 6)} + Material tensor + k: Float + Material constant, yield limit + Δt: Float + time increment + σ_begin:Array{Float64, 6} + Stress vector in Voigt notation + +Returns +------- + Array{Float64, 7}, return values for solver +""" +function G(params, dϵ, C, k, σ_begin) + + # Creating wrapper for gradient + yield(pars) = vonMisesYield(pars, k) + dfdσ = ForwardDiff.gradient(yield) + + # Stress rate + dσ = params[1:6] + + σ_tot = [vec(σ_begin); 0.0] + params + + # Calculating plastic strain rate + dϵp = params[end] * dfdσ(σ_tot) + + # Calculating equations + function_1 = dσ - C * (dϵ - dϵp[1:6]) + function_2 = yield(σ_tot) + [vec(function_1); function_2] +end + + + +""" +Function which calculates the stress. Also handles if any yielding happens + +Parameters +---------- + dϵ: Array{Float64, 6} + Strain rate vector in Voigt notation + Δt: Float + time increment + σ: Array{Float64, 6} + Last stress vector in Voigt notation + C: Array{Float64, (6, 6)} + Material tensor + k: Float + Material constant, yield limit + +Returns +------- + Tuple + Plastic strain rate dϵᵖ and new stress vector σ +""" +function calculate_stress!(dϵ, mat::State) + σ = mat.σ + C = mat.C + σ_y = mat.σ_y + # Test stress + σ_tria = σ + C * dϵ + + # Calculating yield + yield = vonMisesYield(σ_tria, σ_y) + + if yield > 0 + # Yielding happened + # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values + initial_guess = [vec(σ_tria - σ); 0.1] + f(σ_) = G(σ_, dϵ, C, σ_y, σ) + df = ForwardDiff.jacobian(f) + + # Calculating root + result = nlsolve(not_in_place(f, df), initial_guess).zero + + mat.σ += result[1:6] + else + mat.σ = vec(σ_tria) + end +end + +steps = 1000 +strain_max = 0.003 +num_cycles = 3 + +ϵ_tot = zeros(Float64, (steps, 6)) +ϵ_tot2 = zeros(Float64, (steps, 6)) +ϵ_tot3 = zeros(Float64, (steps, 6)) + +# Adding only strain in x-axis and counting for the poisson effect +ϵ_tot[:, 1] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)) +ϵ_tot[:, 2] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν +ϵ_tot[:, 3] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν +ϵ_tot[:, 4] = strain_max / 10 * sin(2 * pi * linspace(0, num_cycles, steps)) + +ϵ_last = zeros(Float64, (6)) +ϵᵖ = zeros(Float64, (6)) +σ = zeros(Float64, (6, 1)) +σy = 200.0 +ss = Float64[] +ee = Float64[] + +eig_stress = zeros(Float64, (3, 3)) +eig_vals = zeros(Float64, (steps, 3)) +function fill_tensor(a, b) + a[1, 1] = b[1] + a[2, 2] = b[2] + a[3, 3] = b[3] + + a[1, 2] = b[6] + a[1, 3] = b[5] + a[2, 3] = b[4] + + a[2, 1] = b[6] + a[3, 1] = b[5] + a[3, 2] = b[4] +end +mat = State(C, σy, zeros(Float64, 6), zeros(Float64, 6)) + +info("Starting calculation") +for i=1:steps + ϵ_new = reshape(ϵ_tot[i, :, :], (6, 1)) + dϵ = ϵ_new - mat.ϵ + calculate_stress!(dϵ, mat) + mat.ϵ += vec(dϵ) + push!(ss, mat.σ[1]) + push!(ee, mat.ϵ[1]) + + fill_tensor(eig_stress, mat.σ) + eig_vals[i, :] = sort(eigvals(eig_stress)) +end + +# ================ Plotting =================== # +n(θ, ϕ) = [sin(θ)*cos(ϕ) + sin(θ)*sin(ϕ) + cos(θ)] +m(θ, ϕ, χ) = [-sin(ϕ)*cos(χ)-cos(θ)*cos(ϕ)*sin(χ) + cos(ϕ)*cos(χ)-cos(θ)*sin(ϕ)*sin(χ) + sin(θ)*sin(χ)] + +w = [sqrt(2/3) * 200 * m(54.735 * pi / 180, 45 * pi/180, x) for x=0:0.2:(2*pi+0.1)] +base_vec = [1 1 1] / sqrt(3) + +for i=-7:7 + tt = [w[x] + vec(base_vec) + 50 * i for x=1:length(w)] + x = map(x->tt[x][1], collect(1:length(w))) + y = map(x->tt[x][2], collect(1:length(w))) + z = map(x->tt[x][3], collect(1:length(w))) + plot3D(x, y, z, color="blue") +end + + +#n(54.735 * pi / 180, 45 * pi/180) +info("Calculation finished") +#PyPlot.plot(ee, ss) +plot3D(eig_vals[:, 1], eig_vals[:, 2], eig_vals[:, 3], color="red") +PyPlot.title("Stress-Strain curve") +PyPlot.xlabel("Strain") +PyPlot.ylabel("Stress") +PyPlot.grid() +# PyPlot.plot(ee, ss) +PyPlot.show()