Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-08-06 04:21:33 +00:00
Code Issues Packages Projects Releases Wiki Activity

initial dev for ideal plastic material

Browse Source
This commit is contained in:
Olli
2016-10-02 17:34:05 +03:00
parent c79b4d5514
commit 9749778420
7 changed files with 274 additions and 235 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/JuliaFEM.jl
+3
View File
@@ -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
src/elements.jl
+4 -3
View File
@@ -9,14 +9,15 @@ type Element{E<:AbstractElement}
integration_points :: Vector{IP}
fields :: Dict{AbstractString, Field}
properties :: E
dev :: Dict{Any, Any}
end
function Element{E<:AbstractElement}(::Type{E}, id::Int64, connectivity=[])
return Element{E}(id, connectivity, [], Dict(), E())
return Element{E}(id, connectivity, [], Dict(), E(), Dict{Any, Any}())
end
function Element{E<:AbstractElement}(::Type{E}, connectivity=[])
return Element{E}(-1, connectivity, [], Dict(), E())
return Element{E}(-1, connectivity, [], Dict(), E(), Dict{Any, Any}())
end
function getindex(element::Element, field_name::AbstractString)
@@ -71,7 +72,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
"""
src/vonmises.jl → src/materials_plasticity.jl
+42 -65
View File
@@ -1,45 +1,18 @@
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
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-10)
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
dx = -df(x) \ f(x)
norm(dx) < norm_acc && (converged = true; break)
x += dx
end
converged || error("no convergence!")
x
return x
end
type State
@@ -227,15 +200,6 @@ end
"""
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)
@@ -244,6 +208,7 @@ function stress_eq_plane_stress(stress)
# 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
@@ -252,53 +217,65 @@ function vonMisesYieldPlaneStress(stress, stress_y)
stress_eq_plane_stress(stress) - stress_y
end
function vonMisesRootPlaneStress(params, dstrain, C, stress_y, stress_base)
function vonMisesRootPlaneStress(params, dstrain, D, stress_y, stress_base)
# Creating wrapper for gradient
vm_wrap(stress_) = vonMisesYieldPlaneStress(stress_, stress_y)
dfds = ForwardDiff.gradient(vm_wrap)
dfds = x -> ForwardDiff.gradient(vm_wrap, x)
# Stress rate and total strain
dstress = params[1:3]
stress_tot = vec(stress_base) + params[1:3]
stress_tot = 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_1 = dstress - D * (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}})
# http://homes.civil.aau.dk/lda/continuum/plast.pdf
function plastic_von_mises!(stress, dstrain_vec, D, params, Dtan)
# Test stress
dstress = C * dstrain
dstress = vec(D * dstrain_vec)
stress_tria = stress + dstress
stress_y = params["yield_stress"]
# Calculating and checking for yield
yield = vonMisesYieldPlaneStress(stress_tria, stress_y)
if isless(yield, 0.0)
return dstress, zeros(3)
stress[:] = stress_tria[:]
Dtan[:,:] = D[:,:]
else
info("yielded")
# Yielding happened
# Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values
# 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)
f = stress_ -> vonMisesRootPlaneStress(stress_, dstrain_vec, D, stress_y, stress)
df = x -> ForwardDiff.jacobian(f, x)
# Calculating root
results = find_root!(f, df, x)
results = nlsolve(not_in_place(f), x).zero
dstress = results[1:3]
stress_tot = stress + dstress
stress_new = 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
f_ = stress_ -> vonMisesYieldPlaneStress(stress_, stress_y)
dfds_ = x -> ForwardDiff.gradient(f_, x)
dep = plastic_multiplier * dfds_(vec(stress_new))
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[:,:]
end
end
src/problems_elasticity.jl
+25 -2
View File
@@ -71,6 +71,14 @@ 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}})
if !(ip_id in keys(params))
params[ip_id] = Dict{Any, Any}()
params[ip_id]["last_stress"] = [0.0,0.0,0.0]
params[ip_id]["last_strain"] = [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}})
@@ -137,7 +145,22 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity},
error("unknown plane formulation: $(props.formulation)")
end
# calculate stress
stress_vec = D * ([1.0, 1.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"]
(stress_last, strain_last) = get_internal_params(element.dev, ip.id, Val{:planestress})
dstrain_vec = strain_vec - strain_last
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
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 +168,7 @@ 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)
test/test_elasticity_2d_nonhomogeneous_boundary_conditions.jl
-1
View File
@@ -70,4 +70,3 @@ using JuliaFEM.Testing
end
=#
end
test/test_elasticplastic_2d_nonhomogenious_boundary_conditions.jl
+53
View File
@@ -0,0 +1,53 @@
# 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)
element.dev["plasticity"] = Dict{Any, Any}("stress" => JuliaFEM.plastic_von_mises!,
"params" => Dict("yield_stress" => 175.0))
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")
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
test/test_von_mises_material.jl
+147 -164
View File
@@ -1,145 +1,149 @@
#using PyPlot
using PyPlot
using JuliaFEM
using JuliaFEM.Testing
#using JuliaFEM.MaterialModels: stiffnessTensor, calculate_stress, State
#using JuliaFEM.MaterialModels: stiffnessTensorPlaneStress
function test_von_mises_3D_basic()
# 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_planestress_basic()
steps = 1000
strain_max = 0.003
num_cycles = 3
E = 200.0e3
nu = 0.3
strain_max = 0.004
num_cycles = 1.
E = 200000.
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_planestress_basic()
steps = 1000
strain_max = 0.003
num_cycles = 5
E = 200.0e3
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 +155,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,35 +165,18 @@ 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)
strain = zeros(Float64, 3)
params = Dict("yield_stress" => stress_y)
Dtan = C
for i=1:steps
strain_new = reshape(strain_tot[i, :, :], (3, 1))
dstrain = strain_new - strain
stress_inc, lambda = calculate_stress(dstrain,
stress,
C,
stress_y,
Val{:vonMises},
Val{:PlaneStressElasticPlasticProblem})
stress += stress_inc
JuliaFEM.plastic_von_mises!(stress, dstrain, C, params, Dtan)
strain = vec(strain_new)
s1, s2, t12 = stress
se1 = (s1 + s2)/2 + sqrt(((s1 - s2)/2)^2 + t12^2)
@@ -197,14 +184,13 @@ function test_von_mises_planestress_basic()
push!(ss, se1)
push!(ee, se2)
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 +215,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()
=#
end
#plot(x_vals, y_vals)
plot(ee, ss)
show()
# end
# test_von_mises_3D_basic()
#test_von_mises_planestress_basic()