mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-19 01:48:47 +00:00
Von mises ideal plastic material, still couple of bugs
This commit is contained in:
Olli Väinölä
+3
-15
@@ -1,13 +1,11 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
include("vonmises.jl")
|
||||
include("elasticplastic.jl")
|
||||
|
||||
# Elasticity problems
|
||||
|
||||
abstract ElasticityProblem <: AbstractProblem
|
||||
|
||||
abstract PlaneStressElasticityProblem <: ElasticityProblem
|
||||
abstract ElasticityProblem <: AbstractProblem
|
||||
abstract PlaneStressElasticityProblem <: ElasticityProblem
|
||||
|
||||
function get_unknown_field_name{P<:ElasticityProblem}(::Type{P})
|
||||
return "displacement"
|
||||
@@ -17,14 +15,6 @@ function get_unknown_field_type{P<:ElasticityProblem}(::Type{P})
|
||||
return Vector{Float64}
|
||||
end
|
||||
|
||||
function get_unknown_field_name{P<:ElasticPlasticProblem}(::Type{P})
|
||||
return "displacement"
|
||||
end
|
||||
|
||||
function get_unknown_field_type{P<:ElasticPlasticProblem}(::Type{P})
|
||||
return Vector{Float64}
|
||||
end
|
||||
|
||||
# 3D Elasticity problems
|
||||
function ElasticityProblem(dim::Int=3, elements=[])
|
||||
return Problem{ElasticityProblem}("elasticity problem", dim, elements)
|
||||
@@ -67,8 +57,6 @@ https://en.wikipedia.org/wiki/Hooke's_law
|
||||
|
||||
"""
|
||||
function get_residual_vector{P<:ElasticityProblem}(problem::Problem{P}, element::Element, ip::IntegrationPoint, time::Number; variation=nothing)
|
||||
#function get_residual_vector{P<:ElasticPlasticProblem}(problem::Problem{P}, element::Element, ip::IntegrationPoint, time::Number; variation=nothing)
|
||||
|
||||
r = zeros(Float64, problem.dim, length(element))
|
||||
|
||||
J = get_jacobian(element, ip, time)
|
||||
|
||||
@@ -0,0 +1,247 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
include("vonmises.jl")
|
||||
|
||||
# Elasticity problems
|
||||
abstract ElasticPlasticProblem <: AbstractProblem
|
||||
abstract PlaneStressElasticPlasticProblem <: ElasticPlasticProblem
|
||||
|
||||
function get_unknown_field_name{P<:ElasticPlasticProblem}(::Type{P})
|
||||
return "displacement"
|
||||
end
|
||||
|
||||
function get_unknown_field_type{P<:ElasticPlasticProblem}(::Type{P})
|
||||
return Vector{Float64}
|
||||
end
|
||||
|
||||
# 3D Elasticity problems
|
||||
function ElasticPlasticProblem(dim::Int=3, elements=[])
|
||||
return Problem{ElasticPlasticProblem}("elasticplastic problem", dim, elements)
|
||||
end
|
||||
|
||||
# 2D Plane stress elasticity problems
|
||||
function PlaneStressElasticPlasticProblem(dim::Int=2, elements=[])
|
||||
return Problem{PlaneStressElasticPlasticProblem}("plane stress elasticplastic problem", dim, elements)
|
||||
end
|
||||
|
||||
|
||||
function get_residual_vector{P<:PlaneStressElasticPlasticProblem}(problem::Problem{P}, element::Element, ip::IntegrationPoint, time::Number; variation=nothing)
|
||||
r = zeros(Float64, problem.dim, length(element))
|
||||
|
||||
J = get_jacobian(element, ip, time)
|
||||
|
||||
# internal forces
|
||||
if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")
|
||||
if !haskey(element, "integration points")
|
||||
if P == PlaneStressElasticPlasticProblem
|
||||
last_stress = zeros(2,2)
|
||||
last_strain = zeros(2,2)
|
||||
else
|
||||
last_stress = zeros(3,3)
|
||||
last_strain = zeros(3,3)
|
||||
end
|
||||
else
|
||||
for each_ip in element("integration points", time)
|
||||
if isapprox(each_ip.xi, ip.xi)
|
||||
last_stress = ip("stress", time)
|
||||
last_strain = ip("stress", time)
|
||||
break
|
||||
end
|
||||
end
|
||||
end
|
||||
u = element("displacement", time, variation)
|
||||
grad = element(ip, time, Val{:grad})
|
||||
gradu = grad*u
|
||||
|
||||
# deformation gradient
|
||||
F = I + gradu
|
||||
E = 1/2*(F'*F - I)
|
||||
|
||||
# material
|
||||
young = element("youngs modulus", ip, time)
|
||||
poisson = element("poissons ratio", ip, time)
|
||||
stress_y = element("yield stress", time).data
|
||||
mu = young/(2*(1+poisson))
|
||||
lambda = young*poisson/((1+poisson)*(1-2*poisson))
|
||||
if P == PlaneStressElasticPlasticProblem
|
||||
lambda = 2*lambda*mu/(lambda + 2*mu) # <- correction for 2d problems
|
||||
end
|
||||
dstrain = E - last_strain
|
||||
material_model = element("material model", time)
|
||||
s = last_stress
|
||||
de = copy(ForwardDiff.get_value(dstrain))
|
||||
|
||||
if P == PlaneStressElasticPlasticProblem
|
||||
C = stiffnessTensorPlaneStress(young, poisson)
|
||||
s_v = [s[1,1], s[2,2], s[1,2]]
|
||||
de_ = [de[1,1], de[2,2], de[1,2]]
|
||||
problem_stress_type = :PlaneStressElasticPlasticProblem
|
||||
else
|
||||
C = stiffnessTensor(young, poisson)
|
||||
s_v = [s[1,1], s[2,2], s[3,3], s[2,3], s[1,3], s[1,2]]
|
||||
de_ = [de[1,1], de[2,2], de[3,3], de[2,3], de[1,3], de[1,2]]
|
||||
problem_stress_type = :ElasticPlasticProblem
|
||||
end
|
||||
dep = zeros(3)
|
||||
stress_inc, dep = calculate_stress(de_,
|
||||
s_v,
|
||||
C,
|
||||
stress_y,
|
||||
Val{:vonMises},
|
||||
Val{problem_stress_type})
|
||||
nd = [dep[1] dep[3];
|
||||
dep[3] dep[2]]
|
||||
# a = dstrain - nd
|
||||
# b = C * a
|
||||
info("%% ", dep)
|
||||
dif = dstrain - nd
|
||||
mm = [dif[1,1], dif[2,2], dif[1,2]]
|
||||
s_v += C * mm
|
||||
info("--: ", ForwardDiff.get_value(s_v))
|
||||
# stress
|
||||
if P == PlaneStressElasticPlasticProblem
|
||||
S = [s_v[1] s_v[3];
|
||||
s_v[3] s_v[2]]
|
||||
else
|
||||
S = [s_v[1] s_v[6] s_v[5];
|
||||
s_v[6] s_v[2] s_v[4];
|
||||
s_v[5] s_v[4] s_v[3]]
|
||||
end
|
||||
r += F*S*grad*det(J)
|
||||
end
|
||||
|
||||
# external forces - volume load
|
||||
if haskey(element, "displacement load")
|
||||
basis = element(ip, time)
|
||||
b = element("displacement load", ip, time)
|
||||
r -= b*basis*det(J)
|
||||
end
|
||||
|
||||
# external forces - surface traction force
|
||||
if haskey(element, "displacement traction force")
|
||||
basis = element(ip, time)
|
||||
T = element("displacement traction force", ip, time)
|
||||
JT = transpose(J)
|
||||
s = size(JT, 2) == 1 ? JT : cross(JT[:,1], JT[:,2])
|
||||
r -= T*basis*norm(s)
|
||||
end
|
||||
|
||||
return vec(r)
|
||||
end
|
||||
|
||||
|
||||
|
||||
#=
|
||||
function get_residual_vector{P<:ElasticPlasticProblem}(problem::Problem{P}, element::Element, ip::IntegrationPoint, time::Number; variation=nothing)
|
||||
r = zeros(Float64, problem.dim, length(element))
|
||||
|
||||
J = get_jacobian(element, ip, time)
|
||||
|
||||
info("_____________________")
|
||||
# internal forces
|
||||
if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")
|
||||
|
||||
if !haskey(element, "integration points")
|
||||
if P == PlaneStressElasticPlasticProblem
|
||||
last_stress = zeros(2,2)
|
||||
last_strain = zeros(2,2)
|
||||
else
|
||||
last_stress = zeros(3,3)
|
||||
last_strain = zeros(3,3)
|
||||
end
|
||||
else
|
||||
for each_ip in element("integration points", time)
|
||||
if isapprox(each_ip.xi, ip.xi)
|
||||
last_stress = ip("stress", time)
|
||||
last_strain = ip("stress", time)
|
||||
break
|
||||
end
|
||||
end
|
||||
end
|
||||
u = element("displacement", time, variation)
|
||||
grad = element(ip, time, Val{:grad})
|
||||
gradu = grad*u
|
||||
|
||||
# deformation gradient
|
||||
F = I + gradu
|
||||
|
||||
# material
|
||||
young = element("youngs modulus", ip, time)
|
||||
poisson = element("poissons ratio", ip, time)
|
||||
mu = young/(2*(1+poisson))
|
||||
lambda = young*poisson/((1+poisson)*(1-2*poisson))
|
||||
if P == PlaneStressElasticityProblem
|
||||
lambda = 2*lambda*mu/(lambda + 2*mu) # <- correction for 2d problems
|
||||
end
|
||||
|
||||
# strain
|
||||
E = 1/2*(F'*F - I)
|
||||
#E = 1/2*(gradu + gradu') # finite strain (total)
|
||||
|
||||
young = element("youngs modulus", ip, time)
|
||||
poisson = element("poissons ratio", ip, time)
|
||||
stress_y = element("yield stress", time).data
|
||||
dstrain = E - last_strain
|
||||
material_model = element("material model", time)
|
||||
s = last_stress
|
||||
de = ForwardDiff.get_value(dstrain)
|
||||
|
||||
if P == PlaneStressElasticPlasticProblem
|
||||
C = stiffnessTensorPlaneStress(young, poisson)
|
||||
s_v = [s[1,1], s[2,2], s[1,2]]
|
||||
de_ = [de[1,1], de[2,2], de[1,2]]
|
||||
problem_stress_type = :PlaneStressElasticPlasticProblem
|
||||
else
|
||||
C = stiffnessTensor(young, poisson)
|
||||
s_v = [s[1,1], s[2,2], s[3,3], s[2,3], s[1,3], s[1,2]]
|
||||
de_ = [de[1,1], de[2,2], de[3,3], de[2,3], de[1,3], de[1,2]]
|
||||
problem_stress_type = :ElasticPlasticProblem
|
||||
end
|
||||
|
||||
stress_inc, lambda = plastic_multiplier = calculate_stress(de_,
|
||||
s_v,
|
||||
C,
|
||||
stress_y,
|
||||
Val{:vonMises},
|
||||
Val{problem_stress_type})
|
||||
|
||||
# dep = lambda * dfds(s)
|
||||
# upate_material_parameters!(...)
|
||||
s_new = s_v + stress_inc
|
||||
#S = [s_v[1] s_v[6] s_v[5];
|
||||
# s_v[6] s_v[2] s_v[4];
|
||||
# s_v[5] s_v[4] s_v[3]]
|
||||
S = [s_new[1] s_new[3];
|
||||
s_new[3] s_new[2]]
|
||||
# S = C * (E - dep)
|
||||
|
||||
|
||||
info("Stress: ", vec(ForwardDiff.get_value(S)))
|
||||
# stress
|
||||
#S = lambda*trace(E)*I + 2*mu*E
|
||||
|
||||
r += F*S*grad*det(J)
|
||||
|
||||
end
|
||||
|
||||
|
||||
# external forces - volume load
|
||||
if haskey(element, "displacement load")
|
||||
basis = element(ip, time)
|
||||
b = element("displacement load", ip, time)
|
||||
r -= b*basis*det(J)
|
||||
end
|
||||
|
||||
# external forces - surface traction force
|
||||
if haskey(element, "displacement traction force")
|
||||
basis = element(ip, time)
|
||||
T = element("displacement traction force", ip, time)
|
||||
JT = transpose(J)
|
||||
s = size(JT, 2) == 1 ? JT : cross(JT[:,1], JT[:,2])
|
||||
r -= T*basis*norm(s)
|
||||
end
|
||||
|
||||
return vec(r)
|
||||
end
|
||||
=# #fff
|
||||
+96
-4
@@ -1,5 +1,4 @@
|
||||
using ForwardDiff
|
||||
# using NLsolve
|
||||
|
||||
"""
|
||||
Create a isotropic Hooke material matrix C
|
||||
@@ -31,6 +30,17 @@ function stiffnessTensor(E, ν)
|
||||
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}
|
||||
@@ -176,7 +186,9 @@ function calculate_stress!(dstrain, mat::State, ::Type{Val{:vonMises}})
|
||||
end
|
||||
end
|
||||
|
||||
function calculate_stress!(dstrain, stress, C, stress_y, ::Type{Val{:vonMises}})
|
||||
function calculate_stress(dstrain, stress, C, stress_y,
|
||||
::Type{Val{:vonMises}},
|
||||
::Type{Val{:ElasticPlasticProblem}})
|
||||
# Test stress
|
||||
stress_tria = stress + C * dstrain
|
||||
|
||||
@@ -184,7 +196,7 @@ function calculate_stress!(dstrain, stress, C, stress_y, ::Type{Val{:vonMises}})
|
||||
yield = vonMisesYield(stress_tria, stress_y)
|
||||
if isless(yield, 0.0)
|
||||
# stress[i] = stress_tria[i]
|
||||
return 0.0
|
||||
return 0.0
|
||||
else
|
||||
# Yielding happened
|
||||
# Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values
|
||||
@@ -193,7 +205,7 @@ function calculate_stress!(dstrain, stress, C, stress_y, ::Type{Val{:vonMises}})
|
||||
df = ForwardDiff.jacobian(f)
|
||||
|
||||
# Calculating root
|
||||
# result = nlsolve(not_in_place(f, df), initial_guess).zero
|
||||
# result = nlsolve(not_in_place(f, df), initial_guess).zero
|
||||
max_iter = 10
|
||||
converged = false
|
||||
for i=1:5
|
||||
@@ -208,3 +220,83 @@ function calculate_stress!(dstrain, stress, C, stress_y, ::Type{Val{:vonMises}})
|
||||
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)
|
||||
stress_tot = stress + results[1:3]
|
||||
plastic_multiplier = results[end]
|
||||
vm_wrap(stress_) = vonMisesYieldPlaneStress(stress_, stress_y)
|
||||
dfds = ForwardDiff.gradient(vm_wrap)
|
||||
dep = plastic_multiplier * dfds(stress_tot)
|
||||
info("II ", stress_tot)
|
||||
return results[1:3], dep
|
||||
end
|
||||
end
|
||||
|
||||
@@ -0,0 +1,230 @@
|
||||
|
||||
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
module DirectSolverVonMisesTests
|
||||
|
||||
using JuliaFEM.Test
|
||||
using JuliaFEM.Core: Seg2, Quad4
|
||||
using JuliaFEM.Core: PlaneStressElasticityProblem, DirichletProblem
|
||||
using JuliaFEM.Core: PlaneStressElasticPlasticProblem
|
||||
using JuliaFEM.Core: DirectSolver
|
||||
|
||||
function test_solver_multiple_dirichlet_bc()
|
||||
|
||||
N = Vector[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]
|
||||
|
||||
e1 = Quad4([1, 2, 4, 3])
|
||||
e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
|
||||
e1["youngs modulus"] = 900.0
|
||||
e1["poissons ratio"] = 0.25
|
||||
e1["yield stress"] = 100.0
|
||||
e1["material model"] = :vonMises
|
||||
b1 = Seg2([3, 4])
|
||||
b1["geometry"] = Vector[N[3], N[4]]
|
||||
b1["displacement traction force"] = (
|
||||
0.0 => Vector[[0.0, 0.0], [0.0, 0.0]],
|
||||
1.0 => Vector[[0.0, -100.0], [0.0, -100.0]])
|
||||
|
||||
#problem = PlaneStressElasticityProblem()
|
||||
problem = PlaneStressElasticPlasticProblem()
|
||||
push!(problem, e1)
|
||||
push!(problem, b1)
|
||||
|
||||
# boundary elements for dirichlet dx=0
|
||||
dx = Seg2([1, 3])
|
||||
dx["geometry"] = Vector[N[1], N[3]]
|
||||
dx["displacement 1"] = 0.0
|
||||
|
||||
# boundary elements for dirichlet dy=0
|
||||
dy = Seg2([1, 2])
|
||||
dy["geometry"] = Vector[N[1], N[2]]
|
||||
dy["displacement 2"] = 0.0
|
||||
|
||||
problem2 = DirichletProblem("displacement", 2)
|
||||
push!(problem2, dx)
|
||||
|
||||
problem3 = DirichletProblem("displacement", 2)
|
||||
push!(problem3, dy)
|
||||
|
||||
solver = DirectSolver()
|
||||
#solver.dump_matrices = true
|
||||
solver.name = "test_solver_multiple_dirichlet_bc"
|
||||
push!(solver, problem)
|
||||
push!(solver, problem2)
|
||||
push!(solver, problem3)
|
||||
|
||||
# launch solver
|
||||
#norm = solver(0.0)
|
||||
norm = solver(1.0)
|
||||
disp = e1("displacement", [1.0, 1.0], 1.0)
|
||||
info("displacement at tip: $disp")
|
||||
#@test isapprox(disp, [3.17431158889468E-02, -1.38591518927826E-01])
|
||||
|
||||
end
|
||||
test_solver_multiple_dirichlet_bc()
|
||||
|
||||
#=
|
||||
function test_direct_cholesky_with_non_homogeneous_dirichlet_conditions()
|
||||
|
||||
N = Vector[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]
|
||||
|
||||
e1 = Quad4([1, 2, 4, 3])
|
||||
e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
|
||||
e1["youngs modulus"] = 900.0
|
||||
e1["poissons ratio"] = 0.25
|
||||
|
||||
problem = PlaneStressElasticityProblem()
|
||||
push!(problem, e1)
|
||||
|
||||
# left boundary: dx=-0.1, dy=0.1
|
||||
bc1 = Seg2([1, 3])
|
||||
bc1["geometry"] = Vector[N[1], N[3]]
|
||||
bc1["displacement 1"] = -0.1
|
||||
bc1["displacement 2"] = 0.1
|
||||
|
||||
# right boundary: dx=0.2, dy=-0.2
|
||||
bc2 = Seg2([2, 4])
|
||||
bc2["geometry"] = Vector[N[2], N[4]]
|
||||
bc2["displacement 1"] = 0.2
|
||||
bc2["displacement 2"] = -0.2
|
||||
|
||||
boundary = DirichletProblem("displacement", 2)
|
||||
push!(boundary, bc1)
|
||||
push!(boundary, bc2)
|
||||
|
||||
solver = DirectSolver("test_direct_cholesky_with_non_homogeneous_dirichlet_boundary_conditions")
|
||||
push!(solver, problem)
|
||||
push!(solver, boundary)
|
||||
|
||||
# launch solver
|
||||
solver.method = :UMFPACK
|
||||
solver.dump_matrices = true
|
||||
solver.max_iterations = 1
|
||||
iters, status = solver(0.0)
|
||||
# FIXME: solver gives no convergence warning when all dofs are fixed.
|
||||
n1disp = e1("displacement", [-1.0, -1.0], 0.0)
|
||||
n2disp = e1("displacement", [ 1.0, -1.0], 0.0)
|
||||
n3disp = e1("displacement", [-1.0, 1.0], 0.0)
|
||||
n4disp = e1("displacement", [ 1.0, 1.0], 0.0)
|
||||
udisp = [n1disp n2disp n3disp n4disp]
|
||||
info("nodal disp = ", udisp)
|
||||
@test isapprox(n1disp, [-0.1, 0.1])
|
||||
@test isapprox(n3disp, [-0.1, 0.1])
|
||||
@test isapprox(n2disp, [ 0.2, -0.2])
|
||||
@test isapprox(n4disp, [ 0.2, -0.2])
|
||||
@test status == true
|
||||
end
|
||||
#test_direct_cholesky_with_non_homogeneous_dirichlet_conditions()
|
||||
|
||||
function test_solver_no_convergence()
|
||||
|
||||
N = Vector[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]
|
||||
|
||||
e1 = Quad4([1, 2, 4, 3])
|
||||
e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
|
||||
e1["youngs modulus"] = 900.0
|
||||
e1["poissons ratio"] = 0.25
|
||||
b1 = Seg2([3, 4])
|
||||
b1["geometry"] = Vector[N[3], N[4]]
|
||||
b1["displacement traction force"] = Vector[[100.0, 100.0], [100.0, 100.0]]
|
||||
|
||||
problem = PlaneStressElasticityProblem()
|
||||
push!(problem, e1)
|
||||
push!(problem, b1)
|
||||
|
||||
# boundary elements for dirichlet dx=0
|
||||
dx = Seg2([1, 3])
|
||||
dx["geometry"] = Vector[N[1], N[3]]
|
||||
dx["displacement 1"] = 0.0
|
||||
|
||||
# boundary elements for dirichlet dy=0
|
||||
dy = Seg2([1, 2])
|
||||
dy["geometry"] = Vector[N[1], N[2]]
|
||||
dy["displacement 2"] = 0.0
|
||||
|
||||
problem2 = DirichletProblem("displacement", 2)
|
||||
push!(problem2, dx)
|
||||
|
||||
problem3 = DirichletProblem("displacement", 2)
|
||||
push!(problem3, dy)
|
||||
|
||||
solver = DirectSolver()
|
||||
solver.max_iterations = 1
|
||||
push!(solver, problem)
|
||||
push!(solver, problem2)
|
||||
push!(solver, problem3)
|
||||
|
||||
# launch solver
|
||||
iterations, status = solver(0.0)
|
||||
@test status == false
|
||||
end
|
||||
|
||||
|
||||
function test_solver_multiple_bodies_multiple_dirichlet_bc()
|
||||
N = Vector[
|
||||
[0.0, 0.0], [1.0, 0.0],
|
||||
[0.0, 1.0], [1.0, 1.0],
|
||||
[0.0, 2.0], [1.0, 2.0]]
|
||||
|
||||
e1 = Quad4([1, 2, 4, 3])
|
||||
e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
|
||||
e2 = Quad4([3, 4, 6, 5])
|
||||
e2["geometry"] = Vector[N[3], N[4], N[6], N[5]]
|
||||
for el in [e1, e2]
|
||||
el["youngs modulus"] = 900.0
|
||||
el["poissons ratio"] = 0.25
|
||||
end
|
||||
b1 = Seg2([5, 6])
|
||||
b1["geometry"] = Vector[N[5], N[6]]
|
||||
b1["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]]
|
||||
|
||||
body1 = PlaneStressElasticityProblem()
|
||||
push!(body1, e1)
|
||||
|
||||
body2 = PlaneStressElasticityProblem()
|
||||
push!(body2, e2)
|
||||
push!(body2, b1)
|
||||
|
||||
# boundary elements for dirichlet dx=0
|
||||
dx1 = Seg2([1, 3])
|
||||
dx1["geometry"] = Vector[N[1], N[3]]
|
||||
dx2 = Seg2([3, 5])
|
||||
dx2["geometry"] = Vector[N[3], N[5]]
|
||||
for dx in [dx1, dx2]
|
||||
dx["displacement 1"] = 0.0
|
||||
end
|
||||
|
||||
boundary1 = DirichletProblem("displacement", 2)
|
||||
push!(boundary1, dx1)
|
||||
push!(boundary1, dx2)
|
||||
|
||||
# boundary elements for dirichlet dy=0
|
||||
dy1 = Seg2([1, 2])
|
||||
dy1["geometry"] = Vector[N[1], N[2]]
|
||||
dy1["displacement 2"] = 0.0
|
||||
|
||||
boundary2 = DirichletProblem("displacement", 2)
|
||||
push!(boundary2, dy1)
|
||||
|
||||
|
||||
solver = DirectSolver()
|
||||
push!(solver, body1)
|
||||
push!(solver, body2)
|
||||
push!(solver, boundary1)
|
||||
push!(solver, boundary2)
|
||||
|
||||
# launch solver
|
||||
norm = solver(0.0)
|
||||
|
||||
disp = e2("displacement", [1.0, 1.0], 0.0)
|
||||
info("displacement at tip: $disp")
|
||||
# code aster verification, two_elements.comm
|
||||
@test isapprox(disp, [3.17431158889468E-02, -2.77183037855653E-01])
|
||||
|
||||
end
|
||||
|
||||
#test_solver_multiple_bodies_multiple_dirichlet_bc()
|
||||
=#
|
||||
end
|
||||
@@ -1,9 +1,12 @@
|
||||
module VonMisesTests
|
||||
|
||||
using PyPlot
|
||||
using JuliaFEM.MaterialModels: stiffnessTensor, calculate_stress!, State
|
||||
using JuliaFEM.Test
|
||||
using JuliaFEM.MaterialModels: stiffnessTensor, calculate_stress, State
|
||||
using JuliaFEM.MaterialModels: stiffnessTensorPlaneStress
|
||||
|
||||
function test_von_mises_basic()
|
||||
|
||||
function test_von_mises_3D_basic()
|
||||
|
||||
steps = 1000
|
||||
strain_max = 0.003
|
||||
@@ -51,7 +54,7 @@ function test_von_mises_basic()
|
||||
|
||||
info("Starting calculation")
|
||||
tic()
|
||||
#=
|
||||
#=
|
||||
for i=1:steps
|
||||
strain_new = reshape(strain_tot[i, :, :], (6, 1))
|
||||
dstrain = strain_new - mat.strain
|
||||
@@ -76,7 +79,7 @@ function test_von_mises_basic()
|
||||
fill_tensor(eig_stress, stress)
|
||||
eig_vals[i, :] = sort(eigvals(eig_stress))
|
||||
end
|
||||
|
||||
|
||||
toc()
|
||||
# ================ Plotting =================== #
|
||||
n(θ, ϕ) = [sin(θ)*cos(ϕ)
|
||||
@@ -127,6 +130,112 @@ function test_von_mises_basic()
|
||||
PyPlot.show()
|
||||
end
|
||||
|
||||
test_von_mises_basic()
|
||||
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, ν)
|
||||
|
||||
strain_tot = zeros(Float64, (steps, 3))
|
||||
|
||||
# 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_last = zeros(Float64, (3))
|
||||
strain_p = zeros(Float64, (3))
|
||||
stress = zeros(Float64, (3, 1))
|
||||
stress_y = 200.0
|
||||
ss = Float64[]
|
||||
ee = Float64[]
|
||||
|
||||
|
||||
ss2 = Float64[]
|
||||
ee2 = Float64[]
|
||||
|
||||
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)
|
||||
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
|
||||
strain = vec(strain_new)
|
||||
s1, s2, t12 = stress
|
||||
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)
|
||||
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
|
||||
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
|
||||
PyPlot.plot(x_vals, y_vals)
|
||||
PyPlot.plot(ee, ss)
|
||||
PyPlot.grid()
|
||||
PyPlot.show()
|
||||
end
|
||||
|
||||
# test_von_mises_3D_basic()
|
||||
|
||||
test_von_mises_planestress_basic()
|
||||
|
||||
end
|
||||
|
||||
Reference in New Issue
Block a user