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

implementing ideal plastic material...

Browse Source
This commit is contained in:
Olli Väinölä
2015-12-12 15:22:31 +02:00
parent 2560a8791b
commit 1668dd658a
5 changed files with 126 additions and 124 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/abaqus_reader.jl
+1 -1
View File
@@ -15,7 +15,7 @@ element_has_nodes(::Type{Val{:C3D20}}) = 20
element_has_nodes(::Type{Val{:C3D20E}}) = 20
element_has_nodes(::Type{Val{:S3}}) = 3
element_has_type( ::Type{Val{:S3}}) = :Seg3
element_has_type( ::Type{Val{:S3}}) = :Tri3
element_has_nodes(::Type{Val{:STRI65}}) = 6
element_has_type(::Type{Val{:STRI65}}) = :Tri6
src/api_types.jl
+4 -4
View File
@@ -22,14 +22,14 @@ typealias HeatFluxBC NeumannBC
"""
type Material
name :: ASCIIString
scalar_data :: Dict{ASCIIString, Float64}
scalar_data :: Dict{ASCIIString, Any}
end
#Material(name, data) = Material(name, Dict(data))
Material(name) = Material(name, Dict{ASCIIString, Float64}())
Material() = Material("", Dict{ASCIIString, Float64}())
Material(name) = Material(name, Dict{ASCIIString, Any}())
Material() = Material("", Dict{ASCIIString, Any}())
function Base.setindex!{T <: AbstractString }(material::Material, val::Real, name::T)
function Base.setindex!{T <: AbstractString }(material::Material, val, name::T)
material.scalar_data[name] = val
end
src/elasticity.jl
+40 -5
View File
@@ -1,6 +1,8 @@
# 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 ElasticityProblem <: AbstractProblem
@@ -129,6 +131,22 @@ function get_residual_vector{P<:ElasticPlasticProblem}(problem::Problem{P}, elem
# internal forces
if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")
if !haskey(element, "integration points")
last_stress = zeros(3,3)
last_strain = zeros(3,3)
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
# last_ip = get_last_ip(problem, element, ip, time)
# stress_base = last_ip("stress")
u = element("displacement", time, variation)
grad = element(ip, time, Val{:grad})
# gradu = element("displacement", ip, time, Val{:grad}, variation)
@@ -139,16 +157,33 @@ function get_residual_vector{P<:ElasticPlasticProblem}(problem::Problem{P}, elem
young = element("youngs modulus", ip, time)
poisson = element("poissons ratio", ip, time)
C = stiffnessTensor(young, poisson)
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
E = 1/2*(F'*F - I) # large strain
#E = 1/2*(gradu + gradu') # finite strain
S = lambda*trace(E)*I + 2*mu*E
stress_y = element("yield stress", time).data
#E = 1/2*(F'*F - I) # large strain
E = 1/2*(gradu + gradu') # finite strain (total)
dstrain = E - last_strain
material_model = element("material model", time)
s = last_stress
de = ForwardDiff.get_value(dstrain)
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]]
#println("stress: ", s_v)
#println("de: ", de_)
#println(C)
#println("yield stress: ", stress_y)
plastic_multiplier = calculate_stress!(de_, s_v, C, stress_y, Val{:vonMises})
# dep = lambda * dfds(s)
# upate_material_parameters!(...)
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 = C * (E - dep)
#S = lambda*trace(E)*I + 2*mu*E
#J = det(element, ip, time)
#T = J^-1*F*S*F'
src/vonmises.jl
+66 -114
View File
@@ -1,80 +1,5 @@
using ForwardDiff
using NLsolve
function outer_prod(a, b)
out = zeros(3,3,3,3)
for i=1:3
for j=1:3
for k=1:3
for l=1:3
out[i, j, k, l] = a[i, j] * b[k, l]
end
end
end
end
out
end
function double_contr(a, b)
out = zeros(3, 3)
for i=1:3
for j=1:3
for k=1:3
for l=1:3
out[i, j] += a[i,j,k,l] * b[k, l]
end
end
end
end
out
end
"""
Symmetric fourth order identity tensor
Definition can be found from:
http://www.ce.berkeley.edu/~sanjay/ce231mse211/symidentity.pdf
"""
function identity_tensor_symm_4th_order()
my_kron(i,j) = i == j ? 1 : 0
II = zeros(Float64, (3, 3, 3, 3))
for i=1:3
for j=1:3
for k=1:3
for l=1:3
v1 = my_kron(i, k)
v2 = my_kron(j, l)
v3 = my_kron(i, l)
v4 = my_kron(j, k)
II[i,j,k,l] = 0.5 * (v1*v2 + v3*v4)
end
end
end
end
II
end
"""
Fourth order stiffness tensor
C = λ * I ⊗ I + 2 * μ * II
Definition: https://en.wikipedia.org/wiki/Hooke's_law
Literature from tensors and vectors
# http://www.iith.ac.in/~ashok/Maths_Lectures/Tutorial/VectTensColMat.pdf
# https://en.wikipedia.org/wiki/Tensor_product
# http://www.math.psu.edu/yzheng/m597k/m597kL11.pdf
"""
function stiffnessTensor(youngs_modulus, poissons_ratio, ::Type{Val{:isotropic}})
E = youngs_modulus
v = poissons_ratio
I = eye(3)
II = identity_tensor_symm_4th_order()
mu = E/(2*(1+v))
lambda = E*v/((1+v)*(1-2*v))
return lambda * outer_prod(I, I) + 2 * mu * II
end
# using NLsolve
"""
Create a isotropic Hooke material matrix C
@@ -109,20 +34,11 @@ end
type State
C :: Array{Float64, 2}
σ_y :: Float64
σ :: Array{Float64, 1}
ϵ :: Array{Float64, 1}
stress_y :: Float64
stress :: Array{Float64, 1}
strain :: Array{Float64, 1}
end
# using vectors with double contradiction
# http://www-2.unipv.it/compmech/teaching/available/const_mod/const_mod_mat-review_notation.pdf
M = [1 0 0 0 0 0;
0 1 0 0 0 0;
0 0 1 0 0 0;
0 0 0 2 0 0;
0 0 0 0 2 0;
0 0 0 0 0 2;]
"""
Equivalent tensile stress.
@@ -138,9 +54,13 @@ Returns
-------
Float
"""
function σₑ(σ)
s = σ[1:6] - 1/3 * sum([σ[1], σ[2], σ[3]]) * [1 1 1 0 0 0]'
return sqrt(3/2 * s' * M * s)[1]
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
@@ -160,8 +80,8 @@ Returns
-------
Float
"""
function vonMisesYield(σ, k)
σₑ(σ) - k
function vonMisesYield(stress, stress_y)
stress_eq(stress) - stress_y
end
"""
@@ -190,23 +110,22 @@ Returns
-------
Array{Float64, 7}, return values for solver
"""
function vonMisesRoot(params, dϵ, C, σ_y, σ_begin)
function vonMisesRoot(params, dstrain, C, stress_y, stress_base)
# Creating wrapper for gradient
yield_wrap(pars) = vonMisesYield(pars, σ_y)
dfdσ = ForwardDiff.gradient(yield_wrap)
vm_wrap(stress_) = vonMisesYield(stress_, stress_y)
dfds = ForwardDiff.gradient(vm_wrap)
# Stress rate
dσ = params[1:6]
σ_tot = [vec(σ_begin); 0.0] + params
# Stress rate and total strain
dstress = params[1:6]
stress_tot = vec(stress_base) + params[1:6]
# Calculating plastic strain rate
dϵp = params[end] * dfdσ(σ_tot)
dstrain_p = params[end] * dfds(stress_tot)
# Calculating equations
function_1 = dσ - C * (dϵ - dϵp[1:6])
function_2 = yield_wrap(σ_tot)
function_1 = dstress - C * (dstrain - dstrain_p)
function_2 = vm_wrap(stress_tot)
[vec(function_1); function_2]
end
@@ -233,26 +152,59 @@ Returns
Tuple
Plastic strain rate dϵᵖ and new stress vector σ
"""
function calculate_stress!(dϵ, mat::State, ::Type{Val{:vonMises}})
σ = mat.σ
function calculate_stress!(dstrain, mat::State, ::Type{Val{:vonMises}})
stress = mat.stress
C = mat.C
σ_y = mat.σ_y
stress_y = mat.stress_y
# Test stress
σ_tria = σ + C * dϵ
stress_tria = stress + C * dstrain
# Calculating and checking for yield
yield = vonMisesYield(σ_tria, σ_y)
if yield > 0
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 = [vec(σ_tria - σ); 0.1]
f(σ_) = vonMisesRoot(σ_, dϵ, C, σ_y, σ)
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.σ += result[1:6]
else
mat.σ = vec(σ_tria)
mat.stress += result[1:6]
end
end
function calculate_stress!(dstrain, stress, C, stress_y, ::Type{Val{:vonMises}})
# 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
test/test_von_mises_material.jl
+15
View File
@@ -51,6 +51,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
@@ -62,6 +63,20 @@ function test_von_mises_basic()
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(ϕ)