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

Ideal plasticity converged, both 2D and 3D

Browse Source
This commit is contained in:
Olli Väinölä
2016-10-28 15:46:04 +03:00
parent b7ddc4d4c8
commit c3bc708f01
4 changed files with 144 additions and 45 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/materials_plasticity.jl
+4 -2
View File
@@ -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_, time, dt, type_)
function ideal_plasticity!(stress_new, stress_last, dstrain_vec, pstrain, D, params, Dtan, yield_surface_, time, dt, type_)
# Test stress
dstress = vec(D * dstrain_vec)
stress_trial = stress_last + dstress
@@ -83,6 +83,7 @@ function ideal_plasticity!(stress_new, stress_last, dstrain_vec, D, params, Dtan
# Calculating and checking for yield
yield = yield_curr(stress_trial)
if isless(yield, 0.0)
stress_new[:] = stress_trial[:]
Dtan[:,:] = D[:,:]
else
@@ -103,6 +104,7 @@ function ideal_plasticity!(stress_new, stress_last, dstrain_vec, D, params, Dtan
# Updating stress
stress_new[:] = stress_last + dstress
# Calculating plastic strain
dfds_ = x -> ForwardDiff.gradient(yield_curr, x)
dep = plastic_multiplier * dfds_(vec(stress_new))
@@ -114,6 +116,6 @@ function ideal_plasticity!(stress_new, stress_last, dstrain_vec, D, params, Dtan
Dc = (D^-1 + plastic_multiplier * D2g(stress_new))^-1
dfds = dfds_(stress_new)
Dtan[:,:] = Dc - (Dc * dfds * dfds' * Dc) / (dfds' * Dc * dfds)[1]
pstrain[:] = plastic_multiplier * dfds
end
end
src/problems_elasticity.jl
+59 -24
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 initialize_internal_params!(params, ip, ::Type{Val{:type_2d}})
function initialize_internal_params!(params, ip, type_) #::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)))
@@ -80,8 +80,15 @@ function initialize_internal_params!(params, ip, ::Type{Val{:type_2d}})
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])
if type_ == Val{:type_2d}
update!(ip, "stress", 0.0 => [0.0,0.0,0.0])
update!(ip, "strain", 0.0 => [0.0,0.0,0.0])
elseif type_ == Val{:type_3d}
update!(ip, "stress", 0.0 => [0.0,0.0,0.0,0.0,0.0,0.0])
update!(ip, "strain", 0.0 => [0.0,0.0,0.0,0.0,0.0,0.0])
else
error("daa")
end
update!(ip, "prev_time", 0.0 => 0.0)
update!(ip, "params_initialized", 0.0 => true)
end
@@ -93,13 +100,13 @@ function get_keys(element)
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
#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}})
@@ -513,21 +520,48 @@ function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity},
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})
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_3d})
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, 0.0, 0.0, 0.0]
plastic_strain = [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, plastic_strain, D, params, Dtan, yield_surface_, time, dt, Val{:type_3d})
# 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
# 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
@@ -540,6 +574,7 @@ function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity},
:stress12 in props.store_fields && update!(ip, "stress12", time => stress_vec[4])
:stress23 in props.store_fields && update!(ip, "stress23", time => stress_vec[5])
:stress13 in props.store_fields && update!(ip, "stress13", time => stress_vec[6])
:plastic_strain in props.store_fields && update!(ip, "plastic_strain", time => plastic_strain)
Km += w*BL'*Dtan*BL
# material stiffness end
test/test_elasticplastic_3d_linear_with_surface_load.jl
+10 -8
View File
@@ -16,22 +16,23 @@ using JuliaFEM.Testing
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)
element = Element(Hex8, [1, 2, 3, 4, 5, 6, 7, 8])
update!([element], "geometry", nodes)
update!([element], "youngs modulus", 200e3)
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))
"params" => Dict("yield_stress" => 400.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)
push!(elasticity_problem.properties.store_fields, :plastic_strain)
push!(elasticity_problem, element)
bc = Element(Quad4, [1,4,8,5])
update!([bc], "geometry", nodes)
@@ -48,12 +49,13 @@ using JuliaFEM.Testing
push!(boundary_motion, disp)
solver = NonlinearSolver("solve block problem")
solver.time = 1.0
push!(solver, elasticity_problem)
push!(solver, boundary_problem)
push!(solver, boundary_motion)
solver()
disp = element1("displacement", [1.0, 1.0, 1.0], 0.0)
disp = element("displacement", [1.0, 1.0, 1.0], 1.0)
info("displacement at tip: $disp")
u_expected = 2.0 * [-1/3, -1/3, 1.0]
# @test isapprox(disp, u_expected)
test/test_von_mises_material.jl
+71 -11
View File
@@ -64,17 +64,17 @@ function test_von_mises_3D_basic()
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
#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 =================== #
@@ -118,6 +118,66 @@ function test_von_mises_3D_basic()
info("Calculation finished")
# plot3D(ee, ss)
# plot the surface
xx = zeros(10, 10)
yy = zeros(10, 10)
for i=1:10
for j=1:10
xx[i, j] = (i - 5) * 100
yy[i, j] = (j - 5) * 100
end
end
# calculate corresponding z
z = zeros(10, 10)
for i=1:10
for j=1:10
z[i, j] = 1
end
end
# ==================================================================
# plot the surface
plot_surface(xx, yy, z, color="blue")
stress_y = 200.0
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) - stress_y
f(m,c) = [600*cos(c) 600*sin(c)].*m
x_vals = []
max_iter = 100
y_vals = []
for i=0:0.1:(2*pi+0.3)
wf(x) = f(x, i)
t = 0.01
step = 2
merkki = -1
s11, s22 = wf(t)
ii = 0
while (abs(vm(s11, s22)) > 1e-7) && ii < max_iter
val = vm(s11, s22)
if sign(val) != merkki
merkki *= -1
step *= -0.5
end
t += step
s11, s22 = wf(t)
ii += 1
end
push!(x_vals, s11)
push!(y_vals, s22)
end
plot(x_vals, y_vals, zeros(length(y_vals)), color="yellow")
axis("equal")
# ==================================================================
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")