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

separated test from von mises and added MaterialModels namespace

Browse Source
This commit is contained in:
Olli Väinölä
2015-12-10 15:48:39 +02:00
parent 417dc0618d
commit e54c44cdcf
3 changed files with 164 additions and 100 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/JuliaFEM.jl
+4 -1
View File
@@ -19,7 +19,7 @@ include("api.jl")
end
module Preprocess
include("abaqus_reader.jl")
include("abaqus_reader.jl")
include("aster_reader.jl")
end
@@ -32,6 +32,9 @@ module Test
include("test.jl")
end
module MaterialModels
include("vonmises.jl")
end
module Interfaces
include("interfaces.jl")
src/vonMises.jl
+37 -99
View File
@@ -1,13 +1,17 @@
# imports
using ForwardDiff
using NLsolve
using PyPlot
"""
Create a isotropic Hooke material matrix C
More information: # http://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_isotropic.cfm
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
@@ -19,6 +23,22 @@ Returns
-------
Array{Float64, (6,6)}
"""
#=
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
A[i,j,k,l] = my_kron(i, j) * my_kron(k, l)
end
end
end
end
=#
#function outer_product(a, b)
#
#end
function hookeStiffnessTensor(E, ν)
a = 1 - ν
b = 1 - 2*ν
@@ -33,10 +53,12 @@ function hookeStiffnessTensor(E, ν)
end
# Pick material values
E = 200.0e3
ν = 0.3
C = hookeStiffnessTensor(E, ν)
#E = 200.0e3
#mu = 0.3
#C = hookeStiffnessTensor(E, ν)
#lambda =
#nu =
#C = λ * I ⊗ I + 2 * μ * II
type State
C :: Array{Float64, 2}
@@ -121,11 +143,11 @@ Returns
-------
Array{Float64, 7}, return values for solver
"""
function G(params, dϵ, C, k, σ_begin)
function vonMisesRoot(params, dϵ, C, σ_y, σ_begin)
# Creating wrapper for gradient
yield(pars) = vonMisesYield(pars, k)
dfdσ = ForwardDiff.gradient(yield)
yield_wrap(pars) = vonMisesYield(pars, σ_y)
dfdσ = ForwardDiff.gradient(yield_wrap)
# Stress rate
dσ = params[1:6]
@@ -137,14 +159,14 @@ function G(params, dϵ, C, k, σ_begin)
# Calculating equations
function_1 = dσ - C * (dϵ - dϵp[1:6])
function_2 = yield(σ_tot)
function_2 = yield_wrap(σ_tot)
[vec(function_1); function_2]
end
"""
Function which calculates the stress. Also handles if any yielding happens
Stress for ideal plastic von Mises material model
Parameters
----------
@@ -164,110 +186,26 @@ Returns
Tuple
Plastic strain rate dϵᵖ and new stress vector σ
"""
function calculate_stress!(dϵ, mat::State)
function calculate_stress!(dϵ, mat::State, ::Type{Val{:vonMises}})
σ = mat.σ
C = mat.C
σ_y = mat.σ_y
# Test stress
σ_tria = σ + C * dϵ
# Calculating yield
# Calculating and checking for yield
yield = vonMisesYield(σ_tria, σ_y)
if yield > 0
# Yielding happened
# Creating functions for newton: xₙ₊₁ = xₙ - df⁻¹ * f and initial values
initial_guess = [vec(σ_tria - σ); 0.1]
f(σ_) = G(σ_, dϵ, C, σ_y, σ)
f(σ_) = vonMisesRoot(σ_, dϵ, C, σ_y, σ)
df = ForwardDiff.jacobian(f)
# Calculating root
result = nlsolve(not_in_place(f, df), initial_guess).zero
mat.σ += result[1:6]
else
mat.σ = vec(σ_tria)
end
end
steps = 1000
strain_max = 0.003
num_cycles = 3
ϵ_tot = zeros(Float64, (steps, 6))
ϵ_tot2 = zeros(Float64, (steps, 6))
ϵ_tot3 = zeros(Float64, (steps, 6))
# Adding only strain in x-axis and counting for the poisson effect
ϵ_tot[:, 1] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps))
ϵ_tot[:, 2] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν
ϵ_tot[:, 3] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν
ϵ_tot[:, 4] = strain_max / 10 * sin(2 * pi * linspace(0, num_cycles, steps))
ϵ_last = zeros(Float64, (6))
ϵᵖ = zeros(Float64, (6))
σ = zeros(Float64, (6, 1))
σ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, σy, zeros(Float64, 6), zeros(Float64, 6))
info("Starting calculation")
for i=1:steps
ϵ_new = reshape(ϵ_tot[i, :, :], (6, 1))
dϵ = ϵ_new - mat.ϵ
calculate_stress!(dϵ, mat)
mat.ϵ += vec(dϵ)
push!(ss, mat.σ[1])
push!(ee, mat.ϵ[1])
fill_tensor(eig_stress, mat.σ)
eig_vals[i, :] = sort(eigvals(eig_stress))
end
# ================ 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.2:(2*pi+0.1)]
base_vec = [1 1 1] / sqrt(3)
for i=-7:7
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
#n(54.735 * pi / 180, 45 * pi/180)
info("Calculation finished")
#PyPlot.plot(ee, ss)
plot3D(eig_vals[:, 1], eig_vals[:, 2], eig_vals[:, 3], color="red")
PyPlot.title("Stress-Strain curve")
PyPlot.xlabel("Strain")
PyPlot.ylabel("Stress")
PyPlot.grid()
# PyPlot.plot(ee, ss)
PyPlot.show()
test/test_von_mises_material.jl
+123
View File
@@ -0,0 +1,123 @@
module VonMisesTests
using PyPlot
macro R_str(s)
s
end
using JuliaFEM.MaterialModels: hookeStiffnessTensor, calculate_stress!, State
function test_von_mises_basic()
steps = 1000
strain_max = 0.003
num_cycles = 3
E = 200.0e3
nu = 0.3
ν = 0.3
C = hookeStiffnessTensor(E, ν)
ϵ_tot = zeros(Float64, (steps, 6))
ϵ_tot2 = zeros(Float64, (steps, 6))
ϵ_tot3 = zeros(Float64, (steps, 6))
# Adding only strain in x-axis and counting for the poisson effect
ϵ_tot[:, 1] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps))
ϵ_tot[:, 2] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν
ϵ_tot[:, 3] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν
ϵ_tot[:, 4] = strain_max / 10 * sin(2 * pi * linspace(0, num_cycles, steps))
ϵ_last = zeros(Float64, (6))
ϵᵖ = zeros(Float64, (6))
σ = zeros(Float64, (6, 1))
σ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, σy, zeros(Float64, 6), zeros(Float64, 6))
info("Starting calculation")
tic()
for i=1:steps
ϵ_new = reshape(ϵ_tot[i, :, :], (6, 1))
dϵ = ϵ_new - mat.ϵ
calculate_stress!(dϵ, mat, Val{:vonMises})
mat.ϵ += vec(dϵ)
push!(ss, mat.σ[1])
push!(ee, mat.ϵ[1])
fill_tensor(eig_stress, mat.σ)
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
test_von_mises_basic()
end