diff --git a/src/api_functions.jl b/src/api_functions.jl index 7f5c4c3..a80380b 100644 --- a/src/api_functions.jl +++ b/src/api_functions.jl @@ -5,7 +5,6 @@ function add_boundary_condition!(case::Simulation, bc::NeumannBC) push!(case.neumann_boundary_conditions, bc) end - """ Add node to model and renumber for output """ diff --git a/src/api_types.jl b/src/api_types.jl index 8770f3c..198e43f 100644 --- a/src/api_types.jl +++ b/src/api_types.jl @@ -1,6 +1,7 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + type NeumannBC set_name :: ASCIIString value :: Any diff --git a/src/elasticity.jl b/src/elasticity.jl index a767b43..75dd52c 100644 --- a/src/elasticity.jl +++ b/src/elasticity.jl @@ -4,6 +4,7 @@ # Elasticity problems abstract ElasticityProblem <: AbstractProblem +abstract ElasticPlasticProblem <: AbstractProblem function get_unknown_field_name{P<:ElasticityProblem}(::Type{P}) return "displacement" @@ -13,8 +14,20 @@ 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 + function ElasticityProblem(dim::Int=3, elements=[]) - return Problem{PlaneStressElasticityProblem}(dim, elements) + return Problem{ElasticityProblem}(dim, elements) +end + +function ElasticPlasticProblem(dim::Int=3, elements=[]) + return Problem{ElasticPlasticProblem}(dim, elements) end abstract PlaneStressElasticityProblem <: ElasticityProblem @@ -104,3 +117,60 @@ function get_residual_vector{P<:ElasticityProblem}(problem::Problem{P}, element: return vec(r) end +""" + +""" +function get_residual_vector{P<:ElasticPlasticProblem}(problem::Problem{P}, element::Element, ip::IntegrationPoint, time::Number; variation=nothing) + + + # u = element("displacement", ip, time, variation) + + r = zeros(Float64, problem.dim, length(element)) + + # internal forces + if haskey(element, "youngs modulus") && haskey(element, "poissons ratio") + u = element("displacement", time, variation) + grad = element(ip, time, Val{:grad}) +# gradu = element("displacement", ip, time, Val{:grad}, variation) + gradu = grad*u + + F = I + gradu # deformation gradient + + 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 + + + E = 1/2*(F'*F - I) # large strain + #E = 1/2*(gradu + gradu') # finite strain + S = lambda*trace(E)*I + 2*mu*E + + #J = det(element, ip, time) + #T = J^-1*F*S*F' + #ip["cauchy stress"] = T + #ip["gl strain"] = E + + r += F*S*grad + end + + # external forces - volume load + if haskey(element, "displacement load") + basis = element(ip, time) + b = element("displacement load", ip, time) + r -= b*basis + end + + # external forces - surface traction force + if haskey(element, "displacement traction force") + basis = element(ip, time) + T = element("displacement traction force", ip, time) + r -= T*basis + end + + return vec(r) +end diff --git a/src/vonmises.jl b/src/vonmises.jl index ef82140..1e26bd2 100644 --- a/src/vonmises.jl +++ b/src/vonmises.jl @@ -1,6 +1,7 @@ using ForwardDiff using NLsolve + function outer_prod(a, b) out = zeros(3,3,3,3) for i=1:3 diff --git a/test/test_von_mises_material.jl b/test/test_von_mises_material.jl index 45de507..872a5b8 100644 --- a/test/test_von_mises_material.jl +++ b/test/test_von_mises_material.jl @@ -13,20 +13,20 @@ function test_von_mises_basic() ν = 0.3 C = stiffnessTensor(E, ν) - ϵ_tot = zeros(Float64, (steps, 6)) - ϵ_tot2 = zeros(Float64, (steps, 6)) - ϵ_tot3 = zeros(Float64, (steps, 6)) + 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 - ϵ_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)) + 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)) - ϵ_last = zeros(Float64, (6)) - ϵᵖ = zeros(Float64, (6)) - σ = zeros(Float64, (6, 1)) - σy = 200.0 + strain_last = zeros(Float64, (6)) + strain_p = zeros(Float64, (6)) + stress = zeros(Float64, (6, 1)) + stress_y = 200.0 ss = Float64[] ee = Float64[] @@ -47,19 +47,19 @@ function test_von_mises_basic() a[3, 2] = b[4] end - mat = State(C, σy, zeros(Float64, 6), zeros(Float64, 6)) + mat = State(C, stress_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]) + 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.σ) + fill_tensor(eig_stress, mat.stress) eig_vals[i, :] = sort(eigvals(eig_stress)) end toc()