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]]