mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-02 01:06:23 +00:00
Merge branch 'plastic_material'
This commit is contained in:
Olli
@@ -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
|
||||
|
||||
|
||||
+1
-1
@@ -71,7 +71,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
|
||||
|
||||
"""
|
||||
|
||||
@@ -0,0 +1,119 @@
|
||||
using ForwardDiff
|
||||
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)
|
||||
converged = false
|
||||
for i=1:max_iter
|
||||
dx = -df(x) \ f(x)
|
||||
x += dx
|
||||
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
|
||||
"""
|
||||
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{:type_2d}})
|
||||
s1, s2, t12 = stress
|
||||
# Calculating principal stresses
|
||||
# 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
|
||||
|
||||
"""
|
||||
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, yield_surface_, type_)
|
||||
|
||||
# Creating wrapper for gradient
|
||||
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:end-1]
|
||||
stress_tot = stress_base + dstress
|
||||
|
||||
# Calculating plastic strain rate
|
||||
dstrain_p = params[end] * dfds(stress_tot)
|
||||
|
||||
# Calculating equations
|
||||
function_1 = dstress - D * (dstrain - dstrain_p)
|
||||
function_2 = vm_wrap(stress_tot)
|
||||
[vec(function_1); function_2]
|
||||
end
|
||||
|
||||
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
|
||||
stress_y = params["yield_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_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, yield_surface_, type_)
|
||||
df = x -> ForwardDiff.jacobian(f, x)
|
||||
|
||||
# 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:end-1]
|
||||
plastic_multiplier = results[end]
|
||||
|
||||
# Updating stress
|
||||
stress_new[:] = stress_last + dstress
|
||||
|
||||
# Calculating plastic strain
|
||||
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(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
|
||||
+90
-11
@@ -71,6 +71,35 @@ 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 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
|
||||
end
|
||||
|
||||
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
|
||||
end
|
||||
|
||||
""" Elasticity equations for 2d cases. """
|
||||
function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time, ::Type{Val{:plane}})
|
||||
@@ -83,6 +112,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)
|
||||
|
||||
@@ -90,7 +120,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})
|
||||
@@ -129,15 +158,45 @@ 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
|
||||
|
||||
# calculate stress
|
||||
stress_vec = D * ([1.0, 1.0, 2.0] .* strain_vec)
|
||||
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"]
|
||||
|
||||
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_, time, dt, Val{:type_2d})
|
||||
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 +204,8 @@ 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)
|
||||
@@ -383,7 +443,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)
|
||||
@@ -450,7 +509,28 @@ 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)
|
||||
|
||||
element_keys = get_keys(element)
|
||||
|
||||
if "plasticity" in element_keys
|
||||
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)
|
||||
@@ -461,8 +541,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
|
||||
|
||||
-304
@@ -1,304 +0,0 @@
|
||||
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
|
||||
|
||||
# Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values
|
||||
function find_root!(f, df, x; max_iter=50, norm_acc=1e-10)
|
||||
converged = false
|
||||
for i=1:max_iter
|
||||
dx = df(x) \ -f(x)
|
||||
x += dx
|
||||
norm(dx) < norm_acc && (converged = true; break)
|
||||
end
|
||||
converged || error("no convergence!")
|
||||
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
|
||||
|
||||
##################################################################################
|
||||
# ----- AFTER THIS POINT: VON MISES : PLANE STRESS IMPLEMENTATION ----- #
|
||||
##################################################################################
|
||||
|
||||
"""
|
||||
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)
|
||||
s1, s2, t12 = stress
|
||||
# Calculating principal stresses
|
||||
# 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
|
||||
|
||||
# https://andriandriyana.files.wordpress.com/2008/03/yield_criteria.pdf
|
||||
function vonMisesYieldPlaneStress(stress, stress_y)
|
||||
stress_eq_plane_stress(stress) - stress_y
|
||||
end
|
||||
|
||||
function vonMisesRootPlaneStress(params, dstrain, C, stress_y, stress_base)
|
||||
|
||||
# Creating wrapper for gradient
|
||||
vm_wrap(stress_) = vonMisesYieldPlaneStress(stress_, stress_y)
|
||||
dfds = ForwardDiff.gradient(vm_wrap)
|
||||
|
||||
# Stress rate and total strain
|
||||
dstress = params[1:3]
|
||||
stress_tot = vec(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_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}})
|
||||
# Test stress
|
||||
dstress = C * dstrain
|
||||
stress_tria = stress + dstress
|
||||
|
||||
# Calculating and checking for yield
|
||||
yield = vonMisesYieldPlaneStress(stress_tria, stress_y)
|
||||
if isless(yield, 0.0)
|
||||
return dstress, zeros(3)
|
||||
else
|
||||
info("yielded")
|
||||
# Yielding happened
|
||||
# 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)
|
||||
# Calculating root
|
||||
results = find_root!(f, df, x)
|
||||
dstress = results[1:3]
|
||||
stress_tot = 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
|
||||
end
|
||||
end
|
||||
@@ -70,4 +70,3 @@ using JuliaFEM.Testing
|
||||
end
|
||||
=#
|
||||
end
|
||||
|
||||
|
||||
@@ -0,0 +1,59 @@
|
||||
# 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)
|
||||
|
||||
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
|
||||
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")
|
||||
solver.time = 1.0
|
||||
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
|
||||
@@ -0,0 +1,95 @@
|
||||
# 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])
|
||||
update!([element1], "geometry", nodes)
|
||||
update!([element1], "youngs modulus", 200e3)
|
||||
update!([element1], "poissons ratio", 0.3)
|
||||
|
||||
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
|
||||
push!(elasticity_problem, element1)
|
||||
|
||||
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, 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)
|
||||
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)
|
||||
#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
|
||||
@@ -1,5 +1,6 @@
|
||||
|
||||
#using PyPlot
|
||||
using PyPlot
|
||||
using JuliaFEM
|
||||
using JuliaFEM.Testing
|
||||
#using JuliaFEM.MaterialModels: stiffnessTensor, calculate_stress, State
|
||||
#using JuliaFEM.MaterialModels: stiffnessTensorPlaneStress
|
||||
@@ -13,7 +14,14 @@ function test_von_mises_3D_basic()
|
||||
E = 200.0e3
|
||||
nu = 0.3
|
||||
ν = 0.3
|
||||
C = stiffnessTensor(E, ν)
|
||||
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))
|
||||
@@ -49,34 +57,23 @@ function test_von_mises_3D_basic()
|
||||
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)
|
||||
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
|
||||
calculate_stress!(dstrain, stress, C, stress_y, Val{:vonMises})
|
||||
strain = vec(strain_new)
|
||||
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)
|
||||
fill_tensor(eig_stress, stress_new)
|
||||
eig_vals[i, :] = sort(eigvals(eig_stress))
|
||||
stress_last[:] = stress_new[:]
|
||||
end
|
||||
|
||||
toc()
|
||||
@@ -119,8 +116,8 @@ function test_von_mises_3D_basic()
|
||||
|
||||
|
||||
info("Calculation finished")
|
||||
#PyPlot.plot(ee, ss)
|
||||
#=
|
||||
# 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")
|
||||
@@ -128,18 +125,21 @@ function test_von_mises_3D_basic()
|
||||
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
|
||||
strain_max = 0.004
|
||||
num_cycles = 1.
|
||||
E = 200000.
|
||||
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 +151,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,50 +161,39 @@ 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)
|
||||
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 = 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
|
||||
stress_inc, lambda = calculate_stress(dstrain,
|
||||
stress,
|
||||
C,
|
||||
stress_y,
|
||||
Val{:vonMises},
|
||||
Val{:PlaneStressElasticPlasticProblem})
|
||||
stress += stress_inc
|
||||
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()
|
||||
|
||||
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 +218,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()
|
||||
=#
|
||||
|
||||
plot(x_vals, y_vals)
|
||||
plot(ee, ss)
|
||||
show()
|
||||
end
|
||||
|
||||
# test_von_mises_3D_basic()
|
||||
|
||||
#test_von_mises_planestress_basic()
|
||||
test_von_mises_3D_basic()
|
||||
|
||||
# test_von_mises_planestress_basic()
|
||||
|
||||
Reference in New Issue
Block a user