Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-28 20:46:58 +00:00
Code Issues Packages Projects Releases Wiki Activity

Added 3D formulation, no convergence yet

Browse Source
This commit is contained in:
Olli
2016-10-03 08:47:20 +03:00
parent 206db3a2b3
commit 662eaf2638
5 changed files with 322 additions and 355 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/materials_plasticity.jl
+50 -203
View File
@@ -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
src/problems_elasticity.jl
+41 -15
View File
@@ -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