From 1668dd658ad57be3aa826a61eb12bfcaf5f03390 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olli=20V=C3=A4in=C3=B6l=C3=A4?= Date: Sat, 12 Dec 2015 15:22:31 +0200 Subject: [PATCH] implementing ideal plastic material... --- src/abaqus_reader.jl | 2 +- src/api_types.jl | 8 +- src/elasticity.jl | 45 +++++++- src/vonmises.jl | 180 ++++++++++++-------------------- test/test_von_mises_material.jl | 15 +++ 5 files changed, 126 insertions(+), 124 deletions(-) diff --git a/src/abaqus_reader.jl b/src/abaqus_reader.jl index bd118e5..5dcb625 100644 --- a/src/abaqus_reader.jl +++ b/src/abaqus_reader.jl @@ -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 diff --git a/src/api_types.jl b/src/api_types.jl index 198e43f..fc60d6d 100644 --- a/src/api_types.jl +++ b/src/api_types.jl @@ -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 diff --git a/src/elasticity.jl b/src/elasticity.jl index 75dd52c..05ddffe 100644 --- a/src/elasticity.jl +++ b/src/elasticity.jl @@ -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' diff --git a/src/vonmises.jl b/src/vonmises.jl index 1e26bd2..f032713 100644 --- a/src/vonmises.jl +++ b/src/vonmises.jl @@ -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 + diff --git a/test/test_von_mises_material.jl b/test/test_von_mises_material.jl index 872a5b8..9ece1a4 100644 --- a/test/test_von_mises_material.jl +++ b/test/test_von_mises_material.jl @@ -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(ϕ)