From 662eaf263877c48bedf8962029d03580e37ab894 Mon Sep 17 00:00:00 2001 From: Olli Date: Mon, 3 Oct 2016 08:47:20 +0300 Subject: [PATCH] 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()