From d055875e94254d7090eddb1f5b629e5dee3fde10 Mon Sep 17 00:00:00 2001 From: Olli Date: Sun, 2 Oct 2016 17:34:05 +0300 Subject: [PATCH 1/6] initial dev for ideal plastic material --- src/JuliaFEM.jl | 3 + src/elements.jl | 7 +- src/{vonmises.jl => materials_plasticity.jl} | 107 +++--- src/problems_elasticity.jl | 27 +- ...y_2d_nonhomogeneous_boundary_conditions.jl | 1 - ...c_2d_nonhomogenious_boundary_conditions.jl | 53 +++ test/test_von_mises_material.jl | 311 +++++++++--------- 7 files changed, 274 insertions(+), 235 deletions(-) rename src/{vonmises.jl => materials_plasticity.jl} (76%) create mode 100644 test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 3240258..f5a25f0 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -67,6 +67,9 @@ export Problem, AbstractProblem, FieldProblem, BoundaryProblem, include("problems_elasticity.jl") export Elasticity +include("materials_plasticity.jl") +export plastic_von_mises + include("problems_dirichlet.jl") export Dirichlet diff --git a/src/elements.jl b/src/elements.jl index 59b9c0b..4183c81 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -9,14 +9,15 @@ type Element{E<:AbstractElement} integration_points :: Vector{IP} fields :: Dict{AbstractString, Field} properties :: E + dev :: Dict{Any, Any} end function Element{E<:AbstractElement}(::Type{E}, id::Int64, connectivity=[]) - return Element{E}(id, connectivity, [], Dict(), E()) + return Element{E}(id, connectivity, [], Dict(), E(), Dict{Any, Any}()) end function Element{E<:AbstractElement}(::Type{E}, connectivity=[]) - return Element{E}(-1, connectivity, [], Dict(), E()) + return Element{E}(-1, connectivity, [], Dict(), E(), Dict{Any, Any}()) end function getindex(element::Element, field_name::AbstractString) @@ -71,7 +72,7 @@ julia> el([0.0, 0.0], 0.0, 1) julia> el([0.0, 0.0], 0.0, 2) 2x8 Array{Float64,2}: - 0.25 0.0 0.25 0.0 0.25 0.0 0.25 0.0 + 0.25 0.0 0.25 0.0 0.25 0.0 0.25 0.0 0.0 0.25 0.0 0.25 0.0 0.25 0.0 0.25 """ diff --git a/src/vonmises.jl b/src/materials_plasticity.jl similarity index 76% rename from src/vonmises.jl rename to src/materials_plasticity.jl index ecc477f..257737d 100644 --- a/src/vonmises.jl +++ b/src/materials_plasticity.jl @@ -1,45 +1,18 @@ using ForwardDiff - -""" -Create a isotropic Hooke material matrix C - -More information: http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_isotropic.cfm - https://en.wikipedia.org/wiki/Hooke's_law - http://www.ce.berkeley.edu/~sanjay/ce231mse211/symidentity.pdf -Parameters ----------- - E: Float - Elastic modulus - ν: Float - Poisson constant - -Returns -------- - Array{Float64, (6,6)} -""" -function stiffnessTensor(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 +using NLsolve # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values -function find_root!(f, df, x; max_iter=50, norm_acc=1e-10) +function find_root!(f, df, x; max_iter=50, norm_acc=1e-9) converged = false for i=1:max_iter - dx = df(x) \ -f(x) - x += dx + dx = -df(x) \ f(x) norm(dx) < norm_acc && (converged = true; break) + + x += dx + end converged || error("no convergence!") - x + return x end type State @@ -227,15 +200,6 @@ end """ http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_plane_stress.cfm """ -function stiffnessTensorPlaneStress(E, ν) - a = 1 - ν^2 - b = 1 - ν - multiplier = E / a - return Float64[1 ν 0; - ν 1 0; - 0 0 b].*multiplier -end - # von mises: plane stress # https://andriandriyana.files.wordpress.com/2008/03/yield_criteria.pdf function stress_eq_plane_stress(stress) @@ -244,6 +208,7 @@ function stress_eq_plane_stress(stress) # http://www.engineersedge.com/material_science/principal_vonmises_stress__13418.htm se1 = (s1 + s2)/2 + sqrt(((s1 - s2)/2)^2 + t12^2) se2 = (s1 + s2)/2 - sqrt(((s1 - s2)/2)^2 + t12^2) + return sqrt(se1^2 -se1*se2 + se2^2) end @@ -252,53 +217,65 @@ function vonMisesYieldPlaneStress(stress, stress_y) stress_eq_plane_stress(stress) - stress_y end -function vonMisesRootPlaneStress(params, dstrain, C, stress_y, stress_base) +function vonMisesRootPlaneStress(params, dstrain, D, stress_y, stress_base) # Creating wrapper for gradient vm_wrap(stress_) = vonMisesYieldPlaneStress(stress_, stress_y) - dfds = ForwardDiff.gradient(vm_wrap) + dfds = x -> ForwardDiff.gradient(vm_wrap, x) # Stress rate and total strain dstress = params[1:3] - stress_tot = vec(stress_base) + params[1:3] + stress_tot = stress_base + params[1:3] # Calculating plastic strain rate dstrain_p = params[end] * dfds(stress_tot) # Calculating equations - function_1 = dstress - C * (dstrain - dstrain_p) + function_1 = dstress - D * (dstrain - dstrain_p) function_2 = vm_wrap(stress_tot) [vec(function_1); function_2] end -function calculate_stress(dstrain, stress, C, stress_y, - ::Type{Val{:vonMises}}, - ::Type{Val{:PlaneStressElasticPlasticProblem}}) +# http://homes.civil.aau.dk/lda/continuum/plast.pdf +function plastic_von_mises!(stress, dstrain_vec, D, params, Dtan) # Test stress - dstress = C * dstrain + dstress = vec(D * dstrain_vec) stress_tria = stress + dstress + stress_y = params["yield_stress"] # Calculating and checking for yield yield = vonMisesYieldPlaneStress(stress_tria, stress_y) if isless(yield, 0.0) - return dstress, zeros(3) + stress[:] = stress_tria[:] + Dtan[:,:] = D[:,:] else - info("yielded") - # Yielding happened - # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values + # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ \ f and initial values x = [vec(stress_tria - stress); 0.0] - f(stress_) = vonMisesRootPlaneStress(stress_, dstrain, C, stress_y, stress) - df = ForwardDiff.jacobian(f) + f = stress_ -> vonMisesRootPlaneStress(stress_, dstrain_vec, D, stress_y, stress) + df = x -> ForwardDiff.jacobian(f, x) + # Calculating root - results = find_root!(f, df, x) + results = nlsolve(not_in_place(f), x).zero dstress = results[1:3] - stress_tot = stress + dstress + + stress_new = stress + dstress plastic_multiplier = results[end] - vm_wrap(stress_) = vonMisesYieldPlaneStress(stress_, stress_y) - dfds = ForwardDiff.gradient(vm_wrap) - dep = plastic_multiplier * dfds(vec(stress_tot)) - info("II ", stress_tot) - info(vm_wrap(stress_tot)) - return dstress, dep + f_ = stress_ -> vonMisesYieldPlaneStress(stress_, stress_y) + dfds_ = x -> ForwardDiff.gradient(f_, x) + dep = plastic_multiplier * dfds_(vec(stress_new)) + + D2g = x -> ForwardDiff.hessian(f_, x) + Dc = (D^-1 + plastic_multiplier * D2g(stress_new))^-1 + dfds = dfds_(stress_new) + Dtan = Dc - (Dc * dfds * dfds' * Dc) / (dfds' * Dc * dfds)[1] + println("plastic stress") + println(stress_new) + println(Dtan) + stress[:] = stress_new[:] + # stress[:] = D * dstrain_vec + println("elastic stress") + println(stress) + println(D) + Dtan[:,:] = D[:,:] end end diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index b7dc871..2aad6ab 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -71,6 +71,14 @@ typealias Elasticity2DVolumeElements Union{Tri3, Tri6, Quad4, Quad8, Quad9} typealias Elasticity3DSurfaceElements Union{Poi1, Tri3, Tri6, Quad4, Quad8, Quad9} typealias Elasticity3DVolumeElements Union{Tet4, Wedge6, Hex8, Tet10, Hex20, Hex27} +function get_internal_params(params, ip_id, ::Type{Val{:planestress}}) + if !(ip_id in keys(params)) + params[ip_id] = Dict{Any, Any}() + params[ip_id]["last_stress"] = [0.0,0.0,0.0] + params[ip_id]["last_strain"] = [0.0,0.0,0.0] + end + return (params[ip_id]["last_stress"], params[ip_id]["last_strain"]) +end """ Elasticity equations for 2d cases. """ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time, ::Type{Val{:plane}}) @@ -137,7 +145,22 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, error("unknown plane formulation: $(props.formulation)") end # calculate stress - stress_vec = D * ([1.0, 1.0, 2.0] .* strain_vec) + if "plasticity" in keys(element.dev) + plastic_def = element.dev["plasticity"] + calculate_stress! = plastic_def["stress"] + params = plastic_def["params"] + (stress_last, strain_last) = get_internal_params(element.dev, ip.id, Val{:planestress}) + dstrain_vec = strain_vec - strain_last + Dtan = [0.0 0.0 0.0; + 0.0 0.0 0.0; + 0.0 0.0 0.0] + + calculate_stress!(stress_last, dstrain_vec, D, params, Dtan) + stress_vec = stress_last + else + stress_vec = D * ([1.0, 1.0, 2.0] .* strain_vec) + Dtan = D + end :strain in props.store_fields && update!(ip, "strain", time => strain_vec) :stress in props.store_fields && update!(ip, "stress", time => stress_vec) @@ -145,7 +168,7 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, :stress22 in props.store_fields && update!(ip, "stress22", time => stress_vec[2]) :stress12 in props.store_fields && update!(ip, "stress12", time => stress_vec[3]) - Km += w*BL'*D*BL + Km += w*BL'*Dtan*BL # stress = [stress_vec[1] stress_vec[3]; stress_vec[3] stress_vec[2]] # cauchy_stress = F'*stress*F/det(F) diff --git a/test/test_elasticity_2d_nonhomogeneous_boundary_conditions.jl b/test/test_elasticity_2d_nonhomogeneous_boundary_conditions.jl index d4d3754..268fe87 100644 --- a/test/test_elasticity_2d_nonhomogeneous_boundary_conditions.jl +++ b/test/test_elasticity_2d_nonhomogeneous_boundary_conditions.jl @@ -70,4 +70,3 @@ using JuliaFEM.Testing end =# end - diff --git a/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl b/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl new file mode 100644 index 0000000..f8a1260 --- /dev/null +++ b/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl @@ -0,0 +1,53 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +using JuliaFEM +using JuliaFEM.Preprocess +using JuliaFEM.Testing + +# @testset "2d nonlinear elasticity: test nonhomogeneous boundary conditions and stress calculation" begin + + # field problem + block = Problem(Elasticity, "BLOCK", 2) + block.properties.formulation = :plane_stress + block.properties.finite_strain = true + block.properties.geometric_stiffness = true + + nodes = Dict{Int, Vector{Float64}}( + 1 => [0.0, 0.0], + 2 => [1.0, 0.0], + 3 => [1.0, 1.0], + 4 => [0.0, 1.0]) + + element = Element(Quad4, [1, 2, 3, 4]) + update!(element, "geometry", nodes) + update!(element, "youngs modulus", 288.0) + update!(element, "poissons ratio", 1/3) + element.dev["plasticity"] = Dict{Any, Any}("stress" => JuliaFEM.plastic_von_mises!, + "params" => Dict("yield_stress" => 175.0)) + push!(block, element) + + # boundary conditions + bc = Problem(Dirichlet, "bc", 2, "displacement") + bel1 = Element(Seg2, [1, 2]) + bel2 = Element(Seg2, [3, 4]) + bel3 = Element(Seg2, [4, 1]) + update!([bel1, bel2, bel3], "geometry", nodes) + update!(bel1, "displacement 2", 0.0) + update!(bel2, "displacement 2", 0.5) + update!(bel3, "displacement 1", 0.0) + push!(bc, bel1, bel2, bel3) + + solver = NonlinearSolver("solve block problem") + push!(solver, block, bc) + solver() + + # from code aster + eps_expected = [-2.08333312468287E-01, 6.25000000000000E-01, 0.0] + sig_expected = [ 4.50685020821470E-06, 4.62857140373777E+02, 0.0] + u3_expected = [-2.36237356855269E-01, 5.00000000000000E-01] + + u3 = reshape(block.assembly.u, 2, 4)[:, 3] + info("u3 = $u3") + #@test isapprox(u3, u3_expected, atol=1.0e-5) +# end diff --git a/test/test_von_mises_material.jl b/test/test_von_mises_material.jl index 27012b9..601401b 100644 --- a/test/test_von_mises_material.jl +++ b/test/test_von_mises_material.jl @@ -1,145 +1,149 @@ -#using PyPlot +using PyPlot +using JuliaFEM using JuliaFEM.Testing #using JuliaFEM.MaterialModels: stiffnessTensor, calculate_stress, State #using JuliaFEM.MaterialModels: stiffnessTensorPlaneStress -function test_von_mises_3D_basic() +# function test_von_mises_3D_basic() +# +# steps = 1000 +# strain_max = 0.003 +# num_cycles = 3 +# E = 200.0e3 +# nu = 0.3 +# ν = 0.3 +# C = stiffnessTensor(E, ν) +# +# strain_tot = zeros(Float64, (steps, 6)) +# strain_tot2 = zeros(Float64, (steps, 6)) +# strain_tot3 = zeros(Float64, (steps, 6)) +# +# # Adding only strain in x-axis and counting for the poisson effect +# strain_tot[:, 1] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)) +# strain_tot[:, 2] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν +# strain_tot[:, 3] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν +# strain_tot[:, 4] = strain_max / 10 * sin(2 * pi * linspace(0, num_cycles, steps)) +# +# strain_last = zeros(Float64, (6)) +# strain_p = zeros(Float64, (6)) +# stress = zeros(Float64, (6, 1)) +# stress_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, stress_y, zeros(Float64, 6), zeros(Float64, 6)) +# +# info("Starting calculation") +# tic() +# #= +# for i=1:steps +# strain_new = reshape(strain_tot[i, :, :], (6, 1)) +# dstrain = strain_new - mat.strain +# calculate_stress!(dstrain, mat, Val{:vonMises}) +# mat.strain += vec(dstrain) +# push!(ss, mat.stress[1]) +# push!(ee, mat.strain[1]) +# +# fill_tensor(eig_stress, mat.stress) +# eig_vals[i, :] = sort(eigvals(eig_stress)) +# end +# =# +# stress = zeros(Float64, 6) +# strain = zeros(Float64, 6) +# for i=1:steps +# strain_new = reshape(strain_tot[i, :, :], (6, 1)) +# dstrain = strain_new - strain +# calculate_stress!(dstrain, stress, C, stress_y, Val{:vonMises}) +# strain = vec(strain_new) +# push!(ss, stress[1]) +# push!(ee, strain[1]) +# fill_tensor(eig_stress, stress) +# eig_vals[i, :] = sort(eigvals(eig_stress)) +# end +# +# toc() +# # ================ 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.15:(2*pi+0.1)] +# base_vec = [1 1 1] / sqrt(3) +# +# for i=-5:5 +# 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 +# +# tt = [w[x] + vec(base_vec) + 50 * -5 for x=1:length(w)] +# x_start = map(x->tt[x][1], collect(1:length(w)))[1:5:end] +# y_start = map(x->tt[x][2], collect(1:length(w)))[1:5:end] +# z_start = map(x->tt[x][3], collect(1:length(w)))[1:5:end] +# +# +# tt = [w[x] + vec(base_vec) + 50 * 5 for x=1:length(w)] +# x_end = map(x->tt[x][1], collect(1:length(w)))[1:5:end] +# y_end = map(x->tt[x][2], collect(1:length(w)))[1:5:end] +# z_end = map(x->tt[x][3], collect(1:length(w)))[1:5:end] +# +# for i=1:length(x_start) +# x = [x_start[i], x_end[i]] +# y = [y_start[i], y_end[i]] +# z = [z_start[i], z_end[i]] +# plot3D(x, y, z, color="blue") +# end +# +# +# info("Calculation finished") +# #PyPlot.plot(ee, ss) +# #= +# plot3D(eig_vals[:, 1], eig_vals[:, 2], eig_vals[:, 3], color="red") +# PyPlot.title("Stress path and von Mises yield surface") +# PyPlot.xlabel("Eig Stress 1") +# PyPlot.ylabel("Eig Stress 2") +# PyPlot.zlabel("Eig Stress 3") +# PyPlot.grid() +# PyPlot.show() +# =# +# end + +#function test_von_mises_planestress_basic() steps = 1000 - strain_max = 0.003 - num_cycles = 3 - E = 200.0e3 - nu = 0.3 + strain_max = 0.004 + num_cycles = 1. + E = 200000. + nu = 0.3 ν = 0.3 - C = stiffnessTensor(E, ν) - - strain_tot = zeros(Float64, (steps, 6)) - strain_tot2 = zeros(Float64, (steps, 6)) - strain_tot3 = zeros(Float64, (steps, 6)) - - # Adding only strain in x-axis and counting for the poisson effect - strain_tot[:, 1] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)) - strain_tot[:, 2] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν - strain_tot[:, 3] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν - strain_tot[:, 4] = strain_max / 10 * sin(2 * pi * linspace(0, num_cycles, steps)) - - strain_last = zeros(Float64, (6)) - strain_p = zeros(Float64, (6)) - stress = zeros(Float64, (6, 1)) - stress_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, stress_y, zeros(Float64, 6), zeros(Float64, 6)) - - info("Starting calculation") - tic() - #= - for i=1:steps - strain_new = reshape(strain_tot[i, :, :], (6, 1)) - dstrain = strain_new - mat.strain - calculate_stress!(dstrain, mat, Val{:vonMises}) - mat.strain += vec(dstrain) - push!(ss, mat.stress[1]) - push!(ee, mat.strain[1]) - - fill_tensor(eig_stress, mat.stress) - eig_vals[i, :] = sort(eigvals(eig_stress)) - end - =# - stress = zeros(Float64, 6) - strain = zeros(Float64, 6) - for i=1:steps - strain_new = reshape(strain_tot[i, :, :], (6, 1)) - dstrain = strain_new - strain - calculate_stress!(dstrain, stress, C, stress_y, Val{:vonMises}) - strain = vec(strain_new) - push!(ss, stress[1]) - push!(ee, strain[1]) - fill_tensor(eig_stress, stress) - eig_vals[i, :] = sort(eigvals(eig_stress)) - end - - toc() - # ================ 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.15:(2*pi+0.1)] - base_vec = [1 1 1] / sqrt(3) - - for i=-5:5 - 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 - - tt = [w[x] + vec(base_vec) + 50 * -5 for x=1:length(w)] - x_start = map(x->tt[x][1], collect(1:length(w)))[1:5:end] - y_start = map(x->tt[x][2], collect(1:length(w)))[1:5:end] - z_start = map(x->tt[x][3], collect(1:length(w)))[1:5:end] - - - tt = [w[x] + vec(base_vec) + 50 * 5 for x=1:length(w)] - x_end = map(x->tt[x][1], collect(1:length(w)))[1:5:end] - y_end = map(x->tt[x][2], collect(1:length(w)))[1:5:end] - z_end = map(x->tt[x][3], collect(1:length(w)))[1:5:end] - - for i=1:length(x_start) - x = [x_start[i], x_end[i]] - y = [y_start[i], y_end[i]] - z = [z_start[i], z_end[i]] - plot3D(x, y, z, color="blue") - end - - - info("Calculation finished") - #PyPlot.plot(ee, ss) - #= - plot3D(eig_vals[:, 1], eig_vals[:, 2], eig_vals[:, 3], color="red") - PyPlot.title("Stress path and von Mises yield surface") - PyPlot.xlabel("Eig Stress 1") - PyPlot.ylabel("Eig Stress 2") - PyPlot.zlabel("Eig Stress 3") - PyPlot.grid() - PyPlot.show() - =# -end - -function test_von_mises_planestress_basic() - - steps = 1000 - strain_max = 0.003 - num_cycles = 5 - E = 200.0e3 - nu = 0.3 - ν = 0.3 - C = stiffnessTensorPlaneStress(E, ν) + C = E/((1+nu)*(1-2*nu)) .* [ + 1-nu nu 0 + nu 1-nu 0 + 0 0 (1-2*nu)/2] strain_tot = zeros(Float64, (steps, 3)) @@ -151,7 +155,7 @@ function test_von_mises_planestress_basic() strain_last = zeros(Float64, (3)) strain_p = zeros(Float64, (3)) stress = zeros(Float64, (3, 1)) - stress_y = 200.0 + stress_y = 400 ss = Float64[] ee = Float64[] @@ -161,35 +165,18 @@ function test_von_mises_planestress_basic() eig_stress = zeros(Float64, (3, 3)) eig_vals = zeros(Float64, (steps, 3)) - #mat = State(C, stress_y, zeros(Float64, 6), zeros(Float64, 6)) info("Starting calculation") tic() - #= - for i=1:steps - strain_new = reshape(strain_tot[i, :, :], (6, 1)) - dstrain = strain_new - mat.strain - calculate_stress!(dstrain, mat, Val{:vonMises}) - mat.strain += vec(dstrain) - push!(ss, mat.stress[1]) - push!(ee, mat.strain[1]) - fill_tensor(eig_stress, mat.stress) - eig_vals[i, :] = sort(eigvals(eig_stress)) - end - =# stress = zeros(Float64, 3) strain = zeros(Float64, 3) + params = Dict("yield_stress" => stress_y) + Dtan = C for i=1:steps strain_new = reshape(strain_tot[i, :, :], (3, 1)) dstrain = strain_new - strain - stress_inc, lambda = calculate_stress(dstrain, - stress, - C, - stress_y, - Val{:vonMises}, - Val{:PlaneStressElasticPlasticProblem}) - stress += stress_inc + JuliaFEM.plastic_von_mises!(stress, dstrain, C, params, Dtan) strain = vec(strain_new) s1, s2, t12 = stress se1 = (s1 + s2)/2 + sqrt(((s1 - s2)/2)^2 + t12^2) @@ -197,14 +184,13 @@ function test_von_mises_planestress_basic() push!(ss, se1) push!(ee, se2) end - toc() function vm_upper(a, c) vals = f(a[1], a[2], c) vm(vals[1], vals[2], 200) end - vm(a,b) = sqrt(a^2 - a*b + b^2) - 200 + vm(a,b) = sqrt(a^2 - a*b + b^2) - stress_y f(m,c) = [600*cos(c) 600*sin(c)].*m x_vals = [] max_iter = 100 @@ -229,15 +215,12 @@ function test_von_mises_planestress_basic() push!(x_vals, s11) push!(y_vals, s22) end - #= - PyPlot.plot(x_vals, y_vals) - PyPlot.plot(ee, ss) - PyPlot.grid() - PyPlot.show() - =# -end + + #plot(x_vals, y_vals) + plot(ee, ss) + show() +# end # test_von_mises_3D_basic() #test_von_mises_planestress_basic() - From 206db3a2b3d7091ea6bb0c9867df2560288edb6e Mon Sep 17 00:00:00 2001 From: Olli Date: Sun, 2 Oct 2016 18:43:34 +0300 Subject: [PATCH 2/6] Changed functions forms --- src/materials_plasticity.jl | 415 ++++++++++++++++++------------------ src/problems_elasticity.jl | 7 +- 2 files changed, 206 insertions(+), 216 deletions(-) diff --git a/src/materials_plasticity.jl b/src/materials_plasticity.jl index 257737d..1020acb 100644 --- a/src/materials_plasticity.jl +++ b/src/materials_plasticity.jl @@ -1,208 +1,200 @@ using ForwardDiff -using NLsolve +# using NLsolve # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values -function find_root!(f, df, x; max_iter=50, norm_acc=1e-9) +function find_root!(f, df, x; max_iter=100, norm_acc=1e-9) converged = false + iter_num = 0 for i=1:max_iter dx = -df(x) \ f(x) - norm(dx) < norm_acc && (converged = true; break) - x += dx - + norm(dx) < norm_acc && (converged = true; iter_num = i; break) end - converged || error("no convergence!") + converged || error("No convergence in radial return!") return x end -type State - C :: Array{Float64, 2} - stress_y :: Float64 - stress :: Array{Float64, 1} - strain :: Array{Float64, 1} -end - -""" -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 stress_eq(stress) - stress_ten = [stress[1] stress[6] stress[5]; - stress[6] stress[2] stress[4]; - stress[5] stress[4] stress[3]] - stress_dev = stress_ten - 1/3 * trace(stress_ten) * eye(3) - s = vec(stress_dev) - return sqrt(3/2 * dot(s, s)) -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(stress, stress_y) - stress_eq(stress) - stress_y -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 vonMisesRoot(params, dstrain, C, stress_y, stress_base) - - # Creating wrapper for gradient - vm_wrap(stress_) = vonMisesYield(stress_, stress_y) - dfds = ForwardDiff.gradient(vm_wrap) - - # Stress rate and total strain - dstress = params[1:6] - stress_tot = vec(stress_base) + params[1:6] - - # Calculating plastic strain rate - dstrain_p = params[end] * dfds(stress_tot) - - # Calculating equations - function_1 = dstress - C * (dstrain - dstrain_p) - function_2 = vm_wrap(stress_tot) - [vec(function_1); function_2] -end - - - -""" -Stress for ideal plastic von Mises material model - -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!(dstrain, mat::State, ::Type{Val{:vonMises}}) - stress = mat.stress - C = mat.C - stress_y = mat.stress_y - # Test stress - stress_tria = stress + C * dstrain - - # Calculating and checking for yield - yield = vonMisesYield(stress_tria, stress_y) - if isless(yield, 0.0) - mat.stress = vec(stress_tria) - else - # Yielding happened - # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values - initial_guess = Float64[vec(stress_tria - stress); 0.1] - f(stress_) = vonMisesRoot(stress_, dstrain, C, stress_y, stress) - df = ForwardDiff.jacobian(f) - - # Calculating root - result = nlsolve(not_in_place(f, df), initial_guess).zero - mat.stress += result[1:6] - end -end - -function calculate_stress(dstrain, stress, C, stress_y, - ::Type{Val{:vonMises}}, - ::Type{Val{:ElasticPlasticProblem}}) - # Test stress - stress_tria = stress + C * dstrain - - # Calculating and checking for yield - yield = vonMisesYield(stress_tria, stress_y) - if isless(yield, 0.0) - # stress[i] = stress_tria[i] - return 0.0 - else - # Yielding happened - # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values - x = [vec(stress_tria - stress); 0.0] - f(stress_) = vonMisesRoot(stress_, dstrain, C, stress_y, stress) - df = ForwardDiff.jacobian(f) - - # Calculating root - # result = nlsolve(not_in_place(f, df), initial_guess).zero - max_iter = 10 - converged = false - for i=1:5 - dx = df(x) \ -f(x) - x += dx - # println(x) - norm(dx) < 1e-10 && (converged = true; break) - end - converged || error("no convergence!") - # stress[:] += x[1:6] - return x[end] - end -end +# """ +# 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 stress_eq(stress) +# stress_ten = [stress[1] stress[6] stress[5]; +# stress[6] stress[2] stress[4]; +# stress[5] stress[4] stress[3]] +# stress_dev = stress_ten - 1/3 * trace(stress_ten) * eye(3) +# s = vec(stress_dev) +# return sqrt(3/2 * dot(s, s)) +# 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(stress, stress_y) +# stress_eq(stress) - stress_y +# 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 vonMisesRoot(params, dstrain, C, stress_y, stress_base) +# +# # Creating wrapper for gradient +# vm_wrap(stress_) = vonMisesYield(stress_, stress_y) +# dfds = ForwardDiff.gradient(vm_wrap) +# +# # Stress rate and total strain +# dstress = params[1:6] +# stress_tot = vec(stress_base) + params[1:6] +# +# # Calculating plastic strain rate +# dstrain_p = params[end] * dfds(stress_tot) +# +# # Calculating equations +# function_1 = dstress - C * (dstrain - dstrain_p) +# function_2 = vm_wrap(stress_tot) +# [vec(function_1); function_2] +# end +# +# +# +# """ +# Stress for ideal plastic von Mises material model +# +# 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!(dstrain, mat, ::Type{Val{:vonMises}}) +# stress = mat.stress +# C = mat.C +# stress_y = mat.stress_y +# # Test stress +# stress_tria = stress + C * dstrain +# +# # Calculating and checking for yield +# yield = vonMisesYield(stress_tria, stress_y) +# if isless(yield, 0.0) +# mat.stress = vec(stress_tria) +# else +# # Yielding happened +# # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values +# initial_guess = Float64[vec(stress_tria - stress); 0.1] +# f(stress_) = vonMisesRoot(stress_, dstrain, C, stress_y, stress) +# df = ForwardDiff.jacobian(f) +# +# # Calculating root +# result = nlsolve(not_in_place(f, df), initial_guess).zero +# mat.stress += result[1:6] +# end +# end +# +# function calculate_stress(dstrain, stress, C, stress_y, +# ::Type{Val{:vonMises}}, +# ::Type{Val{:ElasticPlasticProblem}}) +# # Test stress +# stress_tria = stress + C * dstrain +# +# # Calculating and checking for yield +# yield = vonMisesYield(stress_tria, stress_y) +# if isless(yield, 0.0) +# # stress[i] = stress_tria[i] +# return 0.0 +# else +# # Yielding happened +# # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values +# x = [vec(stress_tria - stress); 0.0] +# f(stress_) = vonMisesRoot(stress_, dstrain, C, stress_y, stress) +# df = ForwardDiff.jacobian(f) +# +# # Calculating root +# # result = nlsolve(not_in_place(f, df), initial_guess).zero +# max_iter = 10 +# converged = false +# for i=1:5 +# dx = df(x) \ -f(x) +# x += dx +# # println(x) +# norm(dx) < 1e-10 && (converged = true; break) +# end +# converged || error("no convergence!") +# # stress[:] += x[1:6] +# return x[end] +# end +# end ################################################################################## # ----- AFTER THIS POINT: VON MISES : PLANE STRESS IMPLEMENTATION ----- # ################################################################################## -""" -http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_plane_stress.cfm -""" +#""" +#http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_plane_stress.cfm +#""" # von mises: plane stress # https://andriandriyana.files.wordpress.com/2008/03/yield_criteria.pdf -function stress_eq_plane_stress(stress) +function equivalent_stress(stress, ::Type{Val{:planestress}}) s1, s2, t12 = stress # Calculating principal stresses # http://www.engineersedge.com/material_science/principal_vonmises_stress__13418.htm @@ -213,14 +205,14 @@ function stress_eq_plane_stress(stress) end # https://andriandriyana.files.wordpress.com/2008/03/yield_criteria.pdf -function vonMisesYieldPlaneStress(stress, stress_y) - stress_eq_plane_stress(stress) - stress_y +function yield_function(stress, stress_y, ::Type{Val{:von_mises}}, ::Type{Val{:plane_stress}}) + equivalent_stress(stress, Val{:planestress}) - stress_y end -function vonMisesRootPlaneStress(params, dstrain, D, stress_y, stress_base) +function radial_return(params, dstrain, D, stress_y, stress_base, ::Type{Val{:von_mises}}, ::Type{Val{:plane_stress}}) # Creating wrapper for gradient - vm_wrap(stress_) = vonMisesYieldPlaneStress(stress_, stress_y) + vm_wrap(stress_) = yield_function(stress_, stress_y, Val{:von_mises}, Val{:plane_stress}) dfds = x -> ForwardDiff.gradient(vm_wrap, x) # Stress rate and total strain @@ -236,46 +228,45 @@ function vonMisesRootPlaneStress(params, dstrain, D, stress_y, stress_base) [vec(function_1); function_2] end -# http://homes.civil.aau.dk/lda/continuum/plast.pdf -function plastic_von_mises!(stress, dstrain_vec, D, params, Dtan) +function plastic_von_mises!(stress_new, stress_last, dstrain_vec, D, params, Dtan) # Test stress dstress = vec(D * dstrain_vec) - stress_tria = stress + dstress + stress_tria = stress_last + dstress stress_y = params["yield_stress"] # Calculating and checking for yield - yield = vonMisesYieldPlaneStress(stress_tria, stress_y) + yield = yield_function(stress_tria, stress_y, Val{:von_mises}, Val{:plane_stress}) + if isless(yield, 0.0) - stress[:] = stress_tria[:] + stress_new[:] = stress_tria[:] Dtan[:,:] = D[:,:] else # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ \ f and initial values - x = [vec(stress_tria - stress); 0.0] - f = stress_ -> vonMisesRootPlaneStress(stress_, dstrain_vec, D, stress_y, stress) + f = stress_ -> radial_return(stress_, dstrain_vec, D, stress_y, stress_last, Val{:von_mises}, Val{:plane_stress}) df = x -> ForwardDiff.jacobian(f, x) # Calculating root - results = nlsolve(not_in_place(f), x).zero - dstress = results[1:3] + vals = [vec(stress_tria - stress_last); 0.0] + results = find_root!(f, df, vals) - stress_new = stress + dstress + # extracting results + dstress = results[1:3] plastic_multiplier = results[end] - f_ = stress_ -> vonMisesYieldPlaneStress(stress_, stress_y) + + # Updating stress + stress_new[:] = stress_last + dstress + + # Calculating plastic strain + f_ = stress_ -> yield_function(stress_, stress_y, Val{:von_mises}, Val{:plane_stress}) dfds_ = x -> ForwardDiff.gradient(f_, x) dep = plastic_multiplier * dfds_(vec(stress_new)) + # Equations for consistent tangent matrix can be found from: + # http://homes.civil.aau.dk/lda/continuum/plast.pdf + # equations: 152 & 153 D2g = x -> ForwardDiff.hessian(f_, x) Dc = (D^-1 + plastic_multiplier * D2g(stress_new))^-1 dfds = dfds_(stress_new) - Dtan = Dc - (Dc * dfds * dfds' * Dc) / (dfds' * Dc * dfds)[1] - println("plastic stress") - println(stress_new) - println(Dtan) - stress[:] = stress_new[:] - # stress[:] = D * dstrain_vec - println("elastic stress") - println(stress) - println(D) - Dtan[:,:] = D[:,:] + Dtan[:,:] = Dc - (Dc * dfds * dfds' * Dc) / (dfds' * Dc * dfds)[1] end end diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index 2aad6ab..d27691a 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -151,12 +151,11 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, params = plastic_def["params"] (stress_last, strain_last) = get_internal_params(element.dev, ip.id, Val{:planestress}) dstrain_vec = strain_vec - strain_last + stress_vec = [0.0, 0.0, 0.0] Dtan = [0.0 0.0 0.0; 0.0 0.0 0.0; 0.0 0.0 0.0] - - calculate_stress!(stress_last, dstrain_vec, D, params, Dtan) - stress_vec = stress_last + calculate_stress!(stress_vec, stress_last, dstrain_vec, D, params, Dtan) else stress_vec = D * ([1.0, 1.0, 2.0] .* strain_vec) Dtan = D @@ -169,7 +168,7 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, :stress12 in props.store_fields && update!(ip, "stress12", time => stress_vec[3]) Km += w*BL'*Dtan*BL - + # stress = [stress_vec[1] stress_vec[3]; stress_vec[3] stress_vec[2]] # cauchy_stress = F'*stress*F/det(F) # cauchy_stress = [cauchy_stress[1,1]; cauchy_stress[2,2]; cauchy_stress[1,2]] From 662eaf263877c48bedf8962029d03580e37ab894 Mon Sep 17 00:00:00 2001 From: Olli Date: Mon, 3 Oct 2016 08:47:20 +0300 Subject: [PATCH 3/6] Added 3D formulation, no convergence yet --- src/materials_plasticity.jl | 253 ++++------------ src/problems_elasticity.jl | 56 +++- ...c_2d_nonhomogenious_boundary_conditions.jl | 3 +- ...sticplastic_3d_linear_with_surface_load.jl | 90 ++++++ test/test_von_mises_material.jl | 275 +++++++++--------- 5 files changed, 322 insertions(+), 355 deletions(-) create mode 100644 test/test_elasticplastic_3d_linear_with_surface_load.jl diff --git a/src/materials_plasticity.jl b/src/materials_plasticity.jl index 1020acb..7d984ef 100644 --- a/src/materials_plasticity.jl +++ b/src/materials_plasticity.jl @@ -1,200 +1,42 @@ using ForwardDiff -# using NLsolve +using NLsolve -# Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values -function find_root!(f, df, x; max_iter=100, norm_acc=1e-9) +""" +Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values +""" +function find_root!(f, df, x; max_iter=50, norm_acc=1e-9) converged = false - iter_num = 0 for i=1:max_iter dx = -df(x) \ f(x) x += dx - norm(dx) < norm_acc && (converged = true; iter_num = i; break) + norm(dx) < norm_acc && (converged = true; break) end converged || error("No convergence in radial return!") return x end -# """ -# 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 stress_eq(stress) -# stress_ten = [stress[1] stress[6] stress[5]; -# stress[6] stress[2] stress[4]; -# stress[5] stress[4] stress[3]] -# stress_dev = stress_ten - 1/3 * trace(stress_ten) * eye(3) -# s = vec(stress_dev) -# return sqrt(3/2 * dot(s, s)) -# 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(stress, stress_y) -# stress_eq(stress) - stress_y -# 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 vonMisesRoot(params, dstrain, C, stress_y, stress_base) -# -# # Creating wrapper for gradient -# vm_wrap(stress_) = vonMisesYield(stress_, stress_y) -# dfds = ForwardDiff.gradient(vm_wrap) -# -# # Stress rate and total strain -# dstress = params[1:6] -# stress_tot = vec(stress_base) + params[1:6] -# -# # Calculating plastic strain rate -# dstrain_p = params[end] * dfds(stress_tot) -# -# # Calculating equations -# function_1 = dstress - C * (dstrain - dstrain_p) -# function_2 = vm_wrap(stress_tot) -# [vec(function_1); function_2] -# end -# -# -# -# """ -# Stress for ideal plastic von Mises material model -# -# 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!(dstrain, mat, ::Type{Val{:vonMises}}) -# stress = mat.stress -# C = mat.C -# stress_y = mat.stress_y -# # Test stress -# stress_tria = stress + C * dstrain -# -# # Calculating and checking for yield -# yield = vonMisesYield(stress_tria, stress_y) -# if isless(yield, 0.0) -# mat.stress = vec(stress_tria) -# else -# # Yielding happened -# # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values -# initial_guess = Float64[vec(stress_tria - stress); 0.1] -# f(stress_) = vonMisesRoot(stress_, dstrain, C, stress_y, stress) -# df = ForwardDiff.jacobian(f) -# -# # Calculating root -# result = nlsolve(not_in_place(f, df), initial_guess).zero -# mat.stress += result[1:6] -# end -# end -# -# function calculate_stress(dstrain, stress, C, stress_y, -# ::Type{Val{:vonMises}}, -# ::Type{Val{:ElasticPlasticProblem}}) -# # Test stress -# stress_tria = stress + C * dstrain -# -# # Calculating and checking for yield -# yield = vonMisesYield(stress_tria, stress_y) -# if isless(yield, 0.0) -# # stress[i] = stress_tria[i] -# return 0.0 -# else -# # Yielding happened -# # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values -# x = [vec(stress_tria - stress); 0.0] -# f(stress_) = vonMisesRoot(stress_, dstrain, C, stress_y, stress) -# df = ForwardDiff.jacobian(f) -# -# # Calculating root -# # result = nlsolve(not_in_place(f, df), initial_guess).zero -# max_iter = 10 -# converged = false -# for i=1:5 -# dx = df(x) \ -f(x) -# x += dx -# # println(x) -# norm(dx) < 1e-10 && (converged = true; break) -# end -# converged || error("no convergence!") -# # stress[:] += x[1:6] -# return x[end] -# end -# end +""" +Equivalent tensile stress. -################################################################################## -# ----- AFTER THIS POINT: VON MISES : PLANE STRESS IMPLEMENTATION ----- # -################################################################################## +More info can be found from: https://en.wikipedia.org/wiki/Von_Mises_yield_criterion + Section: Reduced von Mises equation for different stress conditions +""" +function equivalent_stress(stress, ::Type{Val{:type_3d}}) + stress_ten = [stress[1] stress[6] stress[5]; + stress[6] stress[2] stress[4]; + stress[5] stress[4] stress[3]] + stress_dev = stress_ten - 1/3 * trace(stress_ten) * eye(3) + s = vec(stress_dev) + return sqrt(3/2 * dot(s, s)) +end -#""" -#http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_plane_stress.cfm -#""" -# von mises: plane stress -# https://andriandriyana.files.wordpress.com/2008/03/yield_criteria.pdf -function equivalent_stress(stress, ::Type{Val{:planestress}}) +""" +http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_plane_stress.cfm + +von mises: plane stress +https://andriandriyana.files.wordpress.com/2008/03/yield_criteria.pdf +""" +function equivalent_stress(stress, ::Type{Val{:type_2d}}) s1, s2, t12 = stress # Calculating principal stresses # http://www.engineersedge.com/material_science/principal_vonmises_stress__13418.htm @@ -204,20 +46,22 @@ function equivalent_stress(stress, ::Type{Val{:planestress}}) return sqrt(se1^2 -se1*se2 + se2^2) end -# https://andriandriyana.files.wordpress.com/2008/03/yield_criteria.pdf -function yield_function(stress, stress_y, ::Type{Val{:von_mises}}, ::Type{Val{:plane_stress}}) - equivalent_stress(stress, Val{:planestress}) - stress_y +""" +https://andriandriyana.files.wordpress.com/2008/03/yield_criteria.pdf +""" +function yield_function(stress, stress_y, ::Type{Val{:von_mises}}, type_) + equivalent_stress(stress, type_) - stress_y end -function radial_return(params, dstrain, D, stress_y, stress_base, ::Type{Val{:von_mises}}, ::Type{Val{:plane_stress}}) +function radial_return(params, dstrain, D, stress_y, stress_base, yield_surface_, type_) # Creating wrapper for gradient - vm_wrap(stress_) = yield_function(stress_, stress_y, Val{:von_mises}, Val{:plane_stress}) + vm_wrap(stress_) = yield_function(stress_, stress_y, yield_surface_, type_) dfds = x -> ForwardDiff.gradient(vm_wrap, x) # Stress rate and total strain - dstress = params[1:3] - stress_tot = stress_base + params[1:3] + dstress = params[1:end-1] + stress_tot = stress_base + dstress # Calculating plastic strain rate dstrain_p = params[end] * dfds(stress_tot) @@ -228,45 +72,48 @@ function radial_return(params, dstrain, D, stress_y, stress_base, ::Type{Val{:vo [vec(function_1); function_2] end -function plastic_von_mises!(stress_new, stress_last, dstrain_vec, D, params, Dtan) +function ideal_plasticity!(stress_new, stress_last, dstrain_vec, D, params, Dtan, yield_surface_, type_) # Test stress dstress = vec(D * dstrain_vec) - stress_tria = stress_last + dstress + stress_trial = stress_last + dstress stress_y = params["yield_stress"] - # Calculating and checking for yield - yield = yield_function(stress_tria, stress_y, Val{:von_mises}, Val{:plane_stress}) + yield_curr = x -> yield_function(x, stress_y, yield_surface_, type_) + # Calculating and checking for yield + yield = yield_curr(stress_trial) if isless(yield, 0.0) - stress_new[:] = stress_tria[:] + stress_new[:] = stress_trial[:] Dtan[:,:] = D[:,:] else # Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ \ f and initial values - f = stress_ -> radial_return(stress_, dstrain_vec, D, stress_y, stress_last, Val{:von_mises}, Val{:plane_stress}) + f = stress_ -> radial_return(stress_, dstrain_vec, D, stress_y, stress_last, yield_surface_, type_) df = x -> ForwardDiff.jacobian(f, x) - # Calculating root - vals = [vec(stress_tria - stress_last); 0.0] + # Calculating root (two options) + vals = [vec(stress_trial - stress_last); 0.0] + + #results = nlsolve(not_in_place(f), vals).zero results = find_root!(f, df, vals) # extracting results - dstress = results[1:3] + dstress = results[1:end-1] plastic_multiplier = results[end] # Updating stress stress_new[:] = stress_last + dstress # Calculating plastic strain - f_ = stress_ -> yield_function(stress_, stress_y, Val{:von_mises}, Val{:plane_stress}) - dfds_ = x -> ForwardDiff.gradient(f_, x) + dfds_ = x -> ForwardDiff.gradient(yield_curr, x) dep = plastic_multiplier * dfds_(vec(stress_new)) # Equations for consistent tangent matrix can be found from: # http://homes.civil.aau.dk/lda/continuum/plast.pdf # equations: 152 & 153 - D2g = x -> ForwardDiff.hessian(f_, x) + D2g = x -> ForwardDiff.hessian(yield_curr, x) Dc = (D^-1 + plastic_multiplier * D2g(stress_new))^-1 dfds = dfds_(stress_new) Dtan[:,:] = Dc - (Dc * dfds * dfds' * Dc) / (dfds' * Dc * dfds)[1] + end end diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index d27691a..40f32f1 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -71,7 +71,7 @@ typealias Elasticity2DVolumeElements Union{Tri3, Tri6, Quad4, Quad8, Quad9} typealias Elasticity3DSurfaceElements Union{Poi1, Tri3, Tri6, Quad4, Quad8, Quad9} typealias Elasticity3DVolumeElements Union{Tet4, Wedge6, Hex8, Tet10, Hex20, Hex27} -function get_internal_params(params, ip_id, ::Type{Val{:planestress}}) +function get_internal_params(params, ip_id, ::Type{Val{:type_2d}}) if !(ip_id in keys(params)) params[ip_id] = Dict{Any, Any}() params[ip_id]["last_stress"] = [0.0,0.0,0.0] @@ -80,6 +80,15 @@ function get_internal_params(params, ip_id, ::Type{Val{:planestress}}) return (params[ip_id]["last_stress"], params[ip_id]["last_strain"]) end +function get_internal_params(params, ip_id, ::Type{Val{:type_3d}}) + if !(ip_id in keys(params)) + params[ip_id] = Dict{Any, Any}() + params[ip_id]["last_stress"] = [0.0,0.0,0.0,0.0,0.0,0.0] + params[ip_id]["last_strain"] = [0.0,0.0,0.0,0.0,0.0,0.0] + end + return (params[ip_id]["last_stress"], params[ip_id]["last_strain"]) +end + """ Elasticity equations for 2d cases. """ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time, ::Type{Val{:plane}}) @@ -91,6 +100,7 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, Km = zeros(dim*nnodes, dim*nnodes) Kg = zeros(dim*nnodes, dim*nnodes) f = zeros(dim*nnodes) + Dtan = zeros(3,3) for ip in get_integration_points(element) @@ -137,10 +147,10 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, nu 1.0 0.0 0.0 0.0 (1.0-nu)/2.0] elseif props.formulation == :plane_strain - D = E/((1+nu)*(1-2*nu)) .* [ - 1-nu nu 0 - nu 1-nu 0 - 0 0 (1-2*nu)/2] + D = E/((1.0+nu)*(1.0-2.0*nu)) .* [ + 1.0-nu nu 0.0 + nu 1.0-nu 0.0 + 0.0 0.0 (1.0-2.0*nu)/2.0] else error("unknown plane formulation: $(props.formulation)") end @@ -148,17 +158,15 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, if "plasticity" in keys(element.dev) plastic_def = element.dev["plasticity"] calculate_stress! = plastic_def["stress"] + yield_surface_ = plastic_def["yield_surface"] params = plastic_def["params"] - (stress_last, strain_last) = get_internal_params(element.dev, ip.id, Val{:planestress}) + (stress_last, strain_last) = get_internal_params(element.dev, ip.id, Val{:type_2d}) dstrain_vec = strain_vec - strain_last stress_vec = [0.0, 0.0, 0.0] - Dtan = [0.0 0.0 0.0; - 0.0 0.0 0.0; - 0.0 0.0 0.0] - calculate_stress!(stress_vec, stress_last, dstrain_vec, D, params, Dtan) + calculate_stress!(stress_vec, stress_last, dstrain_vec, D, params, Dtan, yield_surface_, Val{:type_2d}) else stress_vec = D * ([1.0, 1.0, 2.0] .* strain_vec) - Dtan = D + Dtan[:,:] = D[:,:] end :strain in props.store_fields && update!(ip, "strain", time => strain_vec) @@ -168,7 +176,7 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, :stress12 in props.store_fields && update!(ip, "stress12", time => stress_vec[3]) Km += w*BL'*Dtan*BL - + # stress = [stress_vec[1] stress_vec[3]; stress_vec[3] stress_vec[2]] # cauchy_stress = F'*stress*F/det(F) # cauchy_stress = [cauchy_stress[1,1]; cauchy_stress[2,2]; cauchy_stress[1,2]] @@ -472,7 +480,26 @@ function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, 0.0 0.0 0.0 0.5-nu 0.0 0.0 0.0 0.0 0.0 0.0 0.5-nu 0.0 0.0 0.0 0.0 0.0 0.0 0.5-nu] - stress_vec = D * ([1.0, 1.0, 1.0, 2.0, 2.0, 2.0].*strain_vec) + + if "plasticity" in keys(element.dev) + plastic_def = element.dev["plasticity"] + calculate_stress! = plastic_def["stress"] + params = plastic_def["params"] + yield_surface_ = plastic_def["yield_surface"] + (stress_last, strain_last) = get_internal_params(element.dev, ip.id, Val{:type_3d}) + dstrain_vec = strain_vec - strain_last + stress_vec = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + Dtan = [0.0 0.0 0.0 0.0 0.0 0.0; + 0.0 0.0 0.0 0.0 0.0 0.0; + 0.0 0.0 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 0.0; + 0.0 0.0 0.0 0.0 0.0 0.0; + 0.0 0.0 0.0 0.0 0.0 0.0] + calculate_stress!(stress_vec, stress_last, dstrain_vec, D, params, Dtan, yield_surface_, Val{:type_3d}) + else + stress_vec = D * ([1.0, 1.0, 1.0, 2.0, 2.0, 2.0].*strain_vec) + Dtan = D + end :strain in props.store_fields && update!(ip, "strain", time => strain_vec) :stress in props.store_fields && update!(ip, "stress", time => stress_vec) @@ -483,8 +510,7 @@ function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, :stress23 in props.store_fields && update!(ip, "stress23", time => stress_vec[5]) :stress13 in props.store_fields && update!(ip, "stress13", time => stress_vec[6]) - Km += w*BL'*D*BL - + Km += w*BL'*Dtan*BL # material stiffness end if props.geometric_stiffness diff --git a/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl b/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl index f8a1260..60c6ba5 100644 --- a/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl +++ b/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl @@ -23,7 +23,8 @@ using JuliaFEM.Testing update!(element, "geometry", nodes) update!(element, "youngs modulus", 288.0) update!(element, "poissons ratio", 1/3) - element.dev["plasticity"] = Dict{Any, Any}("stress" => JuliaFEM.plastic_von_mises!, + element.dev["plasticity"] = Dict{Any, Any}("stress" => JuliaFEM.ideal_plasticity!, + "yield_surface" => Val{:von_mises}, "params" => Dict("yield_stress" => 175.0)) push!(block, element) diff --git a/test/test_elasticplastic_3d_linear_with_surface_load.jl b/test/test_elasticplastic_3d_linear_with_surface_load.jl new file mode 100644 index 0000000..226d88d --- /dev/null +++ b/test/test_elasticplastic_3d_linear_with_surface_load.jl @@ -0,0 +1,90 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +using JuliaFEM +using JuliaFEM.Preprocess +using JuliaFEM.Testing + +#@testset "test continuum 3d linear elasticity with surface load" begin + nodes = Dict{Int64, Node}( + 1 => [0.0, 0.0, 0.0], + 2 => [1.0, 0.0, 0.0], + 3 => [1.0, 1.0, 0.0], + 4 => [0.0, 1.0, 0.0], + 5 => [0.0, 0.0, 1.0], + 6 => [1.0, 0.0, 1.0], + 7 => [1.0, 1.0, 1.0], + 8 => [0.0, 1.0, 1.0]) + + element1 = Element(Hex8, [1, 2, 3, 4, 5, 6, 7, 8]) + element2 = Element(Quad4, [5, 6, 7, 8]) + update!([element1, element2], "geometry", nodes) + update!([element1], "youngs modulus", 288.0) + update!([element1], "poissons ratio", 1/3) + update!([element2], "displacement traction force 3", 288.0) + update!([element1], "displacement load 3", 576.0) + + element1.dev["plasticity"] = Dict{Any, Any}("stress" => JuliaFEM.ideal_plasticity!, + "yield_surface" => Val{:von_mises}, + "params" => Dict("yield_stress" => 570.0)) + + elasticity_problem = Problem(Elasticity, "solve continuum block", 3) + elasticity_problem.properties.finite_strain = false + elasticity_problem.properties.geometric_stiffness = false + push!(elasticity_problem, element1) + push!(elasticity_problem, element2) + + symxy = Element(Quad4, [1, 2, 3, 4]) + symxz = Element(Quad4, [1, 2, 6, 5]) + symyz = Element(Quad4, [1, 4, 8, 5]) + update!([symxy, symxz, symyz], "geometry", nodes) + symyz["displacement 1"] = 0.0 + symxz["displacement 2"] = 0.0 + symxy["displacement 3"] = 0.0 + boundary_problem = Problem(Dirichlet, "symmetry boundary conditions", 3, "displacement") + push!(boundary_problem, symxy, symxz, symyz) + + solver = NonlinearSolver("solve block problem") + push!(solver, elasticity_problem, boundary_problem) + solver() + + disp = element1("displacement", [1.0, 1.0, 1.0], 0.0) + info("displacement at tip: $disp") + u_expected = 2.0 * [-1/3, -1/3, 1.0] + @test isapprox(disp, u_expected) +#end + +# function solve_rod_model_elasticity(eltype) +# fn = Pkg.dir("JuliaFEM") * "/test/testdata/rod_short.med" +# mesh = aster_read_mesh(fn, eltype) +# element_sets = join(keys(mesh.element_sets), ", ") +# info("element sets: $element_sets") +# p1 = Problem(Elasticity, "rod", 3) +# p2 = Problem(Elasticity, "trac", 3) +# p3 = Problem(Dirichlet, "fixed", 3, "displacement") +# p4 = Problem(Dirichlet, "fixed", 3, "displacement") +# p5 = Problem(Dirichlet, "fixed", 3, "displacement") +# p1.elements = create_elements(mesh, "ROD") +# p2.elements = create_elements(mesh, "FACE2") +# p3.elements = create_elements(mesh, "FACE1") +# p4.elements = create_elements(mesh, "FACE3") +# p5.elements = create_elements(mesh, "FACE5") +# update!(p1, "youngs modulus", 96.0) +# update!(p1, "poissons ratio", 1/3) +# update!(p2, "displacement traction force 1", 96.0) +# update!(p3, "displacement 1", 0.0) +# update!(p4, "displacement 2", 0.0) +# update!(p5, "displacement 3", 0.0) +# solver = LinearSolver(p1, p2, p3, p4, p5) +# solver() +# u_max = maximum(p1.assembly.u) +# info("$eltype, u_max = $u_max") +# return u_max +# end +# @testset "compare 3d rod to CA solution" begin +# @test isapprox(solve_rod_model_elasticity("Tet4"), 0.2) +# @test isapprox(solve_rod_model_elasticity("Tet10"), 0.2) +# @test isapprox(solve_rod_model_elasticity("Hex8"), 0.2) +# @test isapprox(solve_rod_model_elasticity("Hex20"), 0.2) +# @test isapprox(solve_rod_model_elasticity("Hex27"), 0.2) +# end diff --git a/test/test_von_mises_material.jl b/test/test_von_mises_material.jl index 601401b..313fe79 100644 --- a/test/test_von_mises_material.jl +++ b/test/test_von_mises_material.jl @@ -6,133 +6,129 @@ using JuliaFEM.Testing #using JuliaFEM.MaterialModels: stiffnessTensorPlaneStress -# function test_von_mises_3D_basic() -# -# steps = 1000 -# strain_max = 0.003 -# num_cycles = 3 -# E = 200.0e3 -# nu = 0.3 -# ν = 0.3 -# C = stiffnessTensor(E, ν) -# -# strain_tot = zeros(Float64, (steps, 6)) -# strain_tot2 = zeros(Float64, (steps, 6)) -# strain_tot3 = zeros(Float64, (steps, 6)) -# -# # Adding only strain in x-axis and counting for the poisson effect -# strain_tot[:, 1] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)) -# strain_tot[:, 2] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν -# strain_tot[:, 3] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν -# strain_tot[:, 4] = strain_max / 10 * sin(2 * pi * linspace(0, num_cycles, steps)) -# -# strain_last = zeros(Float64, (6)) -# strain_p = zeros(Float64, (6)) -# stress = zeros(Float64, (6, 1)) -# stress_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, stress_y, zeros(Float64, 6), zeros(Float64, 6)) -# -# info("Starting calculation") -# tic() -# #= -# for i=1:steps -# strain_new = reshape(strain_tot[i, :, :], (6, 1)) -# dstrain = strain_new - mat.strain -# calculate_stress!(dstrain, mat, Val{:vonMises}) -# mat.strain += vec(dstrain) -# push!(ss, mat.stress[1]) -# push!(ee, mat.strain[1]) -# -# fill_tensor(eig_stress, mat.stress) -# eig_vals[i, :] = sort(eigvals(eig_stress)) -# end -# =# -# stress = zeros(Float64, 6) -# strain = zeros(Float64, 6) -# for i=1:steps -# strain_new = reshape(strain_tot[i, :, :], (6, 1)) -# dstrain = strain_new - strain -# calculate_stress!(dstrain, stress, C, stress_y, Val{:vonMises}) -# strain = vec(strain_new) -# push!(ss, stress[1]) -# push!(ee, strain[1]) -# fill_tensor(eig_stress, stress) -# eig_vals[i, :] = sort(eigvals(eig_stress)) -# end -# -# toc() -# # ================ 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.15:(2*pi+0.1)] -# base_vec = [1 1 1] / sqrt(3) -# -# for i=-5:5 -# 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 -# -# tt = [w[x] + vec(base_vec) + 50 * -5 for x=1:length(w)] -# x_start = map(x->tt[x][1], collect(1:length(w)))[1:5:end] -# y_start = map(x->tt[x][2], collect(1:length(w)))[1:5:end] -# z_start = map(x->tt[x][3], collect(1:length(w)))[1:5:end] -# -# -# tt = [w[x] + vec(base_vec) + 50 * 5 for x=1:length(w)] -# x_end = map(x->tt[x][1], collect(1:length(w)))[1:5:end] -# y_end = map(x->tt[x][2], collect(1:length(w)))[1:5:end] -# z_end = map(x->tt[x][3], collect(1:length(w)))[1:5:end] -# -# for i=1:length(x_start) -# x = [x_start[i], x_end[i]] -# y = [y_start[i], y_end[i]] -# z = [z_start[i], z_end[i]] -# plot3D(x, y, z, color="blue") -# end -# -# -# info("Calculation finished") -# #PyPlot.plot(ee, ss) -# #= -# plot3D(eig_vals[:, 1], eig_vals[:, 2], eig_vals[:, 3], color="red") -# PyPlot.title("Stress path and von Mises yield surface") -# PyPlot.xlabel("Eig Stress 1") -# PyPlot.ylabel("Eig Stress 2") -# PyPlot.zlabel("Eig Stress 3") -# PyPlot.grid() -# PyPlot.show() -# =# -# end +function test_von_mises_3D_basic() -#function test_von_mises_planestress_basic() + steps = 1000 + strain_max = 0.003 + num_cycles = 3 + E = 200.0e3 + nu = 0.3 + ν = 0.3 + nu = 0.3 + C = E/((1.0+nu)*(1.0-2.0*nu)) * [ + 1.0-nu nu nu 0.0 0.0 0.0 + nu 1.0-nu nu 0.0 0.0 0.0 + nu nu 1.0-nu 0.0 0.0 0.0 + 0.0 0.0 0.0 0.5-nu 0.0 0.0 + 0.0 0.0 0.0 0.0 0.5-nu 0.0 + 0.0 0.0 0.0 0.0 0.0 0.5-nu] + + strain_tot = zeros(Float64, (steps, 6)) + strain_tot2 = zeros(Float64, (steps, 6)) + strain_tot3 = zeros(Float64, (steps, 6)) + + # Adding only strain in x-axis and counting for the poisson effect + strain_tot[:, 1] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)) + strain_tot[:, 2] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν + strain_tot[:, 3] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν + strain_tot[:, 4] = strain_max / 10 * sin(2 * pi * linspace(0, num_cycles, steps)) + + strain_last = zeros(Float64, (6)) + strain_p = zeros(Float64, (6)) + stress = zeros(Float64, (6, 1)) + stress_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 + + info("Starting calculation") + tic() + params = Dict("yield_stress" => stress_y) + stress_new = zeros(Float64, 6) + stress_last = zeros(Float64, 6) + strain = zeros(Float64, 6) + Dtan = zeros(6,6) + for i=1:steps + strain_new = reshape(strain_tot[i, :, :], (6, 1)) + dstrain = strain_new - strain + JuliaFEM.plastic_von_mises!(stress_new, stress_last, dstrain, C, params, Dtan, Val{:type_3d}) + strain[:] = vec(strain_new)[:] + push!(ss, stress[1]) + push!(ee, strain[1]) + fill_tensor(eig_stress, stress_new) + eig_vals[i, :] = sort(eigvals(eig_stress)) + stress_last[:] = stress_new[:] + end + + toc() + # ================ 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.15:(2*pi+0.1)] + base_vec = [1 1 1] / sqrt(3) + + for i=-5:5 + 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 + + tt = [w[x] + vec(base_vec) + 50 * -5 for x=1:length(w)] + x_start = map(x->tt[x][1], collect(1:length(w)))[1:5:end] + y_start = map(x->tt[x][2], collect(1:length(w)))[1:5:end] + z_start = map(x->tt[x][3], collect(1:length(w)))[1:5:end] + + + tt = [w[x] + vec(base_vec) + 50 * 5 for x=1:length(w)] + x_end = map(x->tt[x][1], collect(1:length(w)))[1:5:end] + y_end = map(x->tt[x][2], collect(1:length(w)))[1:5:end] + z_end = map(x->tt[x][3], collect(1:length(w)))[1:5:end] + + for i=1:length(x_start) + x = [x_start[i], x_end[i]] + y = [y_start[i], y_end[i]] + z = [z_start[i], z_end[i]] + plot3D(x, y, z, color="blue") + end + + + info("Calculation finished") + # plot3D(ee, ss) + + plot3D(eig_vals[:, 1], eig_vals[:, 2], eig_vals[:, 3], color="red") + PyPlot.title("Stress path and von Mises yield surface") + PyPlot.xlabel("Eig Stress 1") + PyPlot.ylabel("Eig Stress 2") + PyPlot.zlabel("Eig Stress 3") + PyPlot.grid() + PyPlot.show() + +end + +function test_von_mises_planestress_basic() steps = 1000 strain_max = 0.004 @@ -169,20 +165,27 @@ using JuliaFEM.Testing info("Starting calculation") tic() - stress = zeros(Float64, 3) + stress_new = zeros(Float64, 3) + stress_last = zeros(Float64, 3) strain = zeros(Float64, 3) + strain_last = zeros(Float64, 3) params = Dict("yield_stress" => stress_y) - Dtan = C + #Dtan = C + Dtan = zeros(3,3) for i=1:steps - strain_new = reshape(strain_tot[i, :, :], (3, 1)) + println("last stress: ", round(stress_last, 2)) + strain_new = vec(strain_tot[i, :, :]) dstrain = strain_new - strain - JuliaFEM.plastic_von_mises!(stress, dstrain, C, params, Dtan) - strain = vec(strain_new) - s1, s2, t12 = stress + println("analytical stress: ", round((C * strain_new)', 2)) + + JuliaFEM.plastic_von_mises!(stress_new, stress_last, dstrain, C, params, Dtan, Val{:type_2d}) + strain[:] = vec(strain_new)[:] + s1, s2, t12 = stress_new se1 = (s1 + s2)/2 + sqrt(((s1 - s2)/2)^2 + t12^2) se2 = (s1 + s2)/2 - sqrt(((s1 - s2)/2)^2 + t12^2) push!(ss, se1) push!(ee, se2) + stress_last[:] = stress_new[:] end toc() @@ -216,11 +219,11 @@ using JuliaFEM.Testing push!(y_vals, s22) end - #plot(x_vals, y_vals) + plot(x_vals, y_vals) plot(ee, ss) show() -# end +end -# test_von_mises_3D_basic() +test_von_mises_3D_basic() -#test_von_mises_planestress_basic() +# test_von_mises_planestress_basic() From e0ca5811db1a0b76e6d6b799c0c3afee17141902 Mon Sep 17 00:00:00 2001 From: Olli Date: Tue, 4 Oct 2016 08:27:54 +0300 Subject: [PATCH 4/6] 3d formulation still not converging --- src/materials_plasticity.jl | 2 + src/problems_elasticity.jl | 3 +- ...sticplastic_3d_linear_with_surface_load.jl | 40 ++++++++++--------- 3 files changed, 24 insertions(+), 21 deletions(-) diff --git a/src/materials_plasticity.jl b/src/materials_plasticity.jl index 7d984ef..128165c 100644 --- a/src/materials_plasticity.jl +++ b/src/materials_plasticity.jl @@ -115,5 +115,7 @@ function ideal_plasticity!(stress_new, stress_last, dstrain_vec, D, params, Dtan dfds = dfds_(stress_new) Dtan[:,:] = Dc - (Dc * dfds * dfds' * Dc) / (dfds' * Dc * dfds)[1] + println("equivalent stress") + println(equivalent_stress(stress_new, type_)) end end diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index 40f32f1..c4aa901 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -108,7 +108,6 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, w = ip.weight*detJ N = element(ip, time) dN = element(ip, time, Val{:Grad}) - # kinematics gradu = element("displacement", ip, time, Val{:Grad}) @@ -413,7 +412,6 @@ end """ Elasticity equations, 3d nonlinear. """ function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum}}) - props = problem.properties dim = get_unknown_field_dimension(problem) nnodes = length(element) @@ -501,6 +499,7 @@ function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, Dtan = D end + println(round(stress_vec, 4)) :strain in props.store_fields && update!(ip, "strain", time => strain_vec) :stress in props.store_fields && update!(ip, "stress", time => stress_vec) :stress11 in props.store_fields && update!(ip, "stress11", time => stress_vec[1]) diff --git a/test/test_elasticplastic_3d_linear_with_surface_load.jl b/test/test_elasticplastic_3d_linear_with_surface_load.jl index 226d88d..010eed1 100644 --- a/test/test_elasticplastic_3d_linear_with_surface_load.jl +++ b/test/test_elasticplastic_3d_linear_with_surface_load.jl @@ -17,41 +17,43 @@ using JuliaFEM.Testing 8 => [0.0, 1.0, 1.0]) element1 = Element(Hex8, [1, 2, 3, 4, 5, 6, 7, 8]) - element2 = Element(Quad4, [5, 6, 7, 8]) - update!([element1, element2], "geometry", nodes) - update!([element1], "youngs modulus", 288.0) - update!([element1], "poissons ratio", 1/3) - update!([element2], "displacement traction force 3", 288.0) - update!([element1], "displacement load 3", 576.0) + update!([element1], "geometry", nodes) + update!([element1], "youngs modulus", 200e3) + update!([element1], "poissons ratio", 0.3) - element1.dev["plasticity"] = Dict{Any, Any}("stress" => JuliaFEM.ideal_plasticity!, + element1.dev["plastdicity"] = Dict{Any, Any}("stress" => JuliaFEM.ideal_plasticity!, "yield_surface" => Val{:von_mises}, - "params" => Dict("yield_stress" => 570.0)) + "params" => Dict("yield_stress" => 500.0)) elasticity_problem = Problem(Elasticity, "solve continuum block", 3) elasticity_problem.properties.finite_strain = false elasticity_problem.properties.geometric_stiffness = false push!(elasticity_problem, element1) - push!(elasticity_problem, element2) - symxy = Element(Quad4, [1, 2, 3, 4]) - symxz = Element(Quad4, [1, 2, 6, 5]) - symyz = Element(Quad4, [1, 4, 8, 5]) - update!([symxy, symxz, symyz], "geometry", nodes) - symyz["displacement 1"] = 0.0 - symxz["displacement 2"] = 0.0 - symxy["displacement 3"] = 0.0 + bc = Element(Quad4, [1,4,8,5]) + update!([bc], "geometry", nodes) + bc["displacement 1"] = 0.0 + bc["displacement 2"] = 0.0 + bc["displacement 3"] = 0.0 boundary_problem = Problem(Dirichlet, "symmetry boundary conditions", 3, "displacement") - push!(boundary_problem, symxy, symxz, symyz) + push!(boundary_problem, bc) + + disp = Element(Quad4, [2,3,7,6]) + update!([disp], "geometry", nodes) + disp["displacement 1"] = 0.002 + boundary_motion = Problem(Dirichlet, "displacement bc", 3, "displacement") + push!(boundary_motion, disp) solver = NonlinearSolver("solve block problem") - push!(solver, elasticity_problem, boundary_problem) + push!(solver, elasticity_problem) + push!(solver, boundary_problem) + push!(solver, boundary_motion) solver() disp = element1("displacement", [1.0, 1.0, 1.0], 0.0) info("displacement at tip: $disp") u_expected = 2.0 * [-1/3, -1/3, 1.0] - @test isapprox(disp, u_expected) + # @test isapprox(disp, u_expected) #end # function solve_rod_model_elasticity(eltype) From 35a2ac28ec2781932459fc8fbf33017c7a2102b7 Mon Sep 17 00:00:00 2001 From: Olli Date: Sun, 9 Oct 2016 13:52:21 +0300 Subject: [PATCH 5/6] added parameters to ip & element fields --- src/elements.jl | 5 +- src/materials_plasticity.jl | 4 +- src/problems_elasticity.jl | 57 ++++++++++++++----- ...c_2d_nonhomogenious_boundary_conditions.jl | 11 +++- 4 files changed, 55 insertions(+), 22 deletions(-) diff --git a/src/elements.jl b/src/elements.jl index 4183c81..5925559 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -9,15 +9,14 @@ type Element{E<:AbstractElement} integration_points :: Vector{IP} fields :: Dict{AbstractString, Field} properties :: E - dev :: Dict{Any, Any} end function Element{E<:AbstractElement}(::Type{E}, id::Int64, connectivity=[]) - return Element{E}(id, connectivity, [], Dict(), E(), Dict{Any, Any}()) + return Element{E}(id, connectivity, [], Dict(), E()) end function Element{E<:AbstractElement}(::Type{E}, connectivity=[]) - return Element{E}(-1, connectivity, [], Dict(), E(), Dict{Any, Any}()) + return Element{E}(-1, connectivity, [], Dict(), E()) end function getindex(element::Element, field_name::AbstractString) diff --git a/src/materials_plasticity.jl b/src/materials_plasticity.jl index 128165c..d1c000a 100644 --- a/src/materials_plasticity.jl +++ b/src/materials_plasticity.jl @@ -72,7 +72,7 @@ function radial_return(params, dstrain, D, stress_y, stress_base, yield_surface_ [vec(function_1); function_2] end -function ideal_plasticity!(stress_new, stress_last, dstrain_vec, D, params, Dtan, yield_surface_, type_) +function ideal_plasticity!(stress_new, stress_last, dstrain_vec, D, params, Dtan, yield_surface_, time, dt, type_) # Test stress dstress = vec(D * dstrain_vec) stress_trial = stress_last + dstress @@ -115,7 +115,5 @@ function ideal_plasticity!(stress_new, stress_last, dstrain_vec, D, params, Dtan dfds = dfds_(stress_new) Dtan[:,:] = Dc - (Dc * dfds * dfds' * Dc) / (dfds' * Dc * dfds)[1] - println("equivalent stress") - println(equivalent_stress(stress_new, type_)) end end diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index c4aa901..764125a 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -71,22 +71,34 @@ typealias Elasticity2DVolumeElements Union{Tri3, Tri6, Quad4, Quad8, Quad9} typealias Elasticity3DSurfaceElements Union{Poi1, Tri3, Tri6, Quad4, Quad8, Quad9} typealias Elasticity3DVolumeElements Union{Tet4, Wedge6, Hex8, Tet10, Hex20, Hex27} -function get_internal_params(params, ip_id, ::Type{Val{:type_2d}}) - if !(ip_id in keys(params)) - params[ip_id] = Dict{Any, Any}() - params[ip_id]["last_stress"] = [0.0,0.0,0.0] - params[ip_id]["last_strain"] = [0.0,0.0,0.0] +function initialize_internal_params!(params, ip, ::Type{Val{:type_2d}}) + param_keys = keys(params) + all_keys = ip.fields.keys + ip_fields = filter(x->isdefined(all_keys, x), collect(1:length(all_keys))) + + if !("params_initialized" in ip_fields) + for key in param_keys + update!(ip, key, 0.0 => params[key]) + end + update!(ip, "stress", 0.0 => [0.0,0.0,0.0]) + update!(ip, "strain", 0.0 => [0.0,0.0,0.0]) + update!(ip, "prev_time", 0.0 => 0.0) + update!(ip, "params_initialized", 0.0 => true) end - return (params[ip_id]["last_stress"], params[ip_id]["last_strain"]) end -function get_internal_params(params, ip_id, ::Type{Val{:type_3d}}) +function get_keys(element) + all_keys = element.fields.keys + idx = filter(x->isdefined(all_keys, x), collect(1:length(all_keys))) + map(x -> all_keys[x], idx) +end + +function initialize_internal_params!(params, ip_id, ::Type{Val{:type_3d}}) if !(ip_id in keys(params)) params[ip_id] = Dict{Any, Any}() params[ip_id]["last_stress"] = [0.0,0.0,0.0,0.0,0.0,0.0] params[ip_id]["last_strain"] = [0.0,0.0,0.0,0.0,0.0,0.0] end - return (params[ip_id]["last_stress"], params[ip_id]["last_strain"]) end """ Elasticity equations for 2d cases. """ @@ -153,16 +165,34 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, else error("unknown plane formulation: $(props.formulation)") end + # calculate stress - if "plasticity" in keys(element.dev) - plastic_def = element.dev["plasticity"] - calculate_stress! = plastic_def["stress"] + element_keys = get_keys(element) + + if "plasticity" in element_keys + plastic_def = element("plasticity")[ip.id] + + calculate_stress! = plastic_def["type"] yield_surface_ = plastic_def["yield_surface"] params = plastic_def["params"] - (stress_last, strain_last) = get_internal_params(element.dev, ip.id, Val{:type_2d}) + + initialize_internal_params!(params, ip, Val{:type_2d}) + + if time == 0.0 + error("Given step time = $(time). Please select time > 0.0") + end + + t_last = ip("prev_time", time) + update!(ip, "prev_time", time => t_last) + + dt = time - t_last + + stress_last = ip("stress", t_last) + strain_last = ip("strain", t_last) + dstrain_vec = strain_vec - strain_last stress_vec = [0.0, 0.0, 0.0] - calculate_stress!(stress_vec, stress_last, dstrain_vec, D, params, Dtan, yield_surface_, Val{:type_2d}) + calculate_stress!(stress_vec, stress_last, dstrain_vec, D, params, Dtan, yield_surface_, time, dt, Val{:type_2d}) else stress_vec = D * ([1.0, 1.0, 2.0] .* strain_vec) Dtan[:,:] = D[:,:] @@ -176,6 +206,7 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, Km += w*BL'*Dtan*BL + # stress = [stress_vec[1] stress_vec[3]; stress_vec[3] stress_vec[2]] # cauchy_stress = F'*stress*F/det(F) # cauchy_stress = [cauchy_stress[1,1]; cauchy_stress[2,2]; cauchy_stress[1,2]] diff --git a/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl b/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl index 60c6ba5..5eee961 100644 --- a/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl +++ b/test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl @@ -23,9 +23,13 @@ using JuliaFEM.Testing update!(element, "geometry", nodes) update!(element, "youngs modulus", 288.0) update!(element, "poissons ratio", 1/3) - element.dev["plasticity"] = Dict{Any, Any}("stress" => JuliaFEM.ideal_plasticity!, - "yield_surface" => Val{:von_mises}, - "params" => Dict("yield_stress" => 175.0)) + + plastic_parameters = Dict{Any, Any}("type" => JuliaFEM.ideal_plasticity!, + "yield_surface" => Val{:von_mises}, + "params" => Dict("yield_stress" => 175.0)) + to_integ_points = Dict() + map(x-> to_integ_points[x] = plastic_parameters, get_connectivity(element)) + update!(element, "plasticity", to_integ_points) push!(block, element) # boundary conditions @@ -40,6 +44,7 @@ using JuliaFEM.Testing push!(bc, bel1, bel2, bel3) solver = NonlinearSolver("solve block problem") + solver.time = 1.0 push!(solver, block, bc) solver() From 739c0cedc8f7735d7f32942bca5ae896c8bb45c4 Mon Sep 17 00:00:00 2001 From: Olli Date: Sun, 9 Oct 2016 14:51:06 +0300 Subject: [PATCH 6/6] fixed 3d, since there is no dev field anymore --- src/problems_elasticity.jl | 5 +++-- ...test_elasticplastic_3d_linear_with_surface_load.jl | 11 +++++++---- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index 764125a..1a12c5a 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -510,7 +510,9 @@ function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, 0.0 0.0 0.0 0.0 0.5-nu 0.0 0.0 0.0 0.0 0.0 0.0 0.5-nu] - if "plasticity" in keys(element.dev) + element_keys = get_keys(element) + + if "plasticity" in element_keys plastic_def = element.dev["plasticity"] calculate_stress! = plastic_def["stress"] params = plastic_def["params"] @@ -530,7 +532,6 @@ function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, Dtan = D end - println(round(stress_vec, 4)) :strain in props.store_fields && update!(ip, "strain", time => strain_vec) :stress in props.store_fields && update!(ip, "stress", time => stress_vec) :stress11 in props.store_fields && update!(ip, "stress11", time => stress_vec[1]) diff --git a/test/test_elasticplastic_3d_linear_with_surface_load.jl b/test/test_elasticplastic_3d_linear_with_surface_load.jl index 010eed1..38f3102 100644 --- a/test/test_elasticplastic_3d_linear_with_surface_load.jl +++ b/test/test_elasticplastic_3d_linear_with_surface_load.jl @@ -21,10 +21,13 @@ using JuliaFEM.Testing update!([element1], "youngs modulus", 200e3) update!([element1], "poissons ratio", 0.3) - element1.dev["plastdicity"] = Dict{Any, Any}("stress" => JuliaFEM.ideal_plasticity!, - "yield_surface" => Val{:von_mises}, - "params" => Dict("yield_stress" => 500.0)) - + plastic_parameters = Dict{Any, Any}("type" => JuliaFEM.ideal_plasticity!, + "yield_surface" => Val{:von_mises}, + "params" => Dict("yield_stress" => 175.0)) + to_integ_points = Dict() + map(x-> to_integ_points[x] = plastic_parameters, get_connectivity(element)) + update!(element, "plasticity", to_integ_points) + elasticity_problem = Problem(Elasticity, "solve continuum block", 3) elasticity_problem.properties.finite_strain = false elasticity_problem.properties.geometric_stiffness = false