mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-28 20:46:58 +00:00
Added 3D formulation, no convergence yet
This commit is contained in:
Olli
@@ -23,7 +23,8 @@ using JuliaFEM.Testing
|
||||
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!,
|
||||
element.dev["plasticity"] = Dict{Any, Any}("stress" => JuliaFEM.ideal_plasticity!,
|
||||
"yield_surface" => Val{:von_mises},
|
||||
"params" => Dict("yield_stress" => 175.0))
|
||||
push!(block, element)
|
||||
|
||||
|
||||
@@ -0,0 +1,90 @@
|
||||
# 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])
|
||||
element2 = Element(Quad4, [5, 6, 7, 8])
|
||||
update!([element1, element2], "geometry", nodes)
|
||||
update!([element1], "youngs modulus", 288.0)
|
||||
update!([element1], "poissons ratio", 1/3)
|
||||
update!([element2], "displacement traction force 3", 288.0)
|
||||
update!([element1], "displacement load 3", 576.0)
|
||||
|
||||
element1.dev["plasticity"] = Dict{Any, Any}("stress" => JuliaFEM.ideal_plasticity!,
|
||||
"yield_surface" => Val{:von_mises},
|
||||
"params" => Dict("yield_stress" => 570.0))
|
||||
|
||||
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, element2)
|
||||
|
||||
symxy = Element(Quad4, [1, 2, 3, 4])
|
||||
symxz = Element(Quad4, [1, 2, 6, 5])
|
||||
symyz = Element(Quad4, [1, 4, 8, 5])
|
||||
update!([symxy, symxz, symyz], "geometry", nodes)
|
||||
symyz["displacement 1"] = 0.0
|
||||
symxz["displacement 2"] = 0.0
|
||||
symxy["displacement 3"] = 0.0
|
||||
boundary_problem = Problem(Dirichlet, "symmetry boundary conditions", 3, "displacement")
|
||||
push!(boundary_problem, symxy, symxz, symyz)
|
||||
|
||||
solver = NonlinearSolver("solve block problem")
|
||||
push!(solver, elasticity_problem, boundary_problem)
|
||||
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
|
||||
+139
-136
@@ -6,133 +6,129 @@ using JuliaFEM.Testing
|
||||
#using JuliaFEM.MaterialModels: stiffnessTensorPlaneStress
|
||||
|
||||
|
||||
# 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_3D_basic()
|
||||
|
||||
#function test_von_mises_planestress_basic()
|
||||
steps = 1000
|
||||
strain_max = 0.003
|
||||
num_cycles = 3
|
||||
E = 200.0e3
|
||||
nu = 0.3
|
||||
ν = 0.3
|
||||
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))
|
||||
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
|
||||
|
||||
info("Starting calculation")
|
||||
tic()
|
||||
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
|
||||
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 =================== #
|
||||
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")
|
||||
# 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")
|
||||
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.004
|
||||
@@ -169,20 +165,27 @@ using JuliaFEM.Testing
|
||||
info("Starting calculation")
|
||||
tic()
|
||||
|
||||
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 = 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
|
||||
JuliaFEM.plastic_von_mises!(stress, dstrain, C, params, Dtan)
|
||||
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()
|
||||
|
||||
@@ -216,11 +219,11 @@ using JuliaFEM.Testing
|
||||
push!(y_vals, s22)
|
||||
end
|
||||
|
||||
#plot(x_vals, y_vals)
|
||||
plot(x_vals, y_vals)
|
||||
plot(ee, ss)
|
||||
show()
|
||||
# end
|
||||
end
|
||||
|
||||
# test_von_mises_3D_basic()
|
||||
test_von_mises_3D_basic()
|
||||
|
||||
#test_von_mises_planestress_basic()
|
||||
# test_von_mises_planestress_basic()
|
||||
|
||||
Reference in New Issue
Block a user