From 8a47ee53cf2d5282c807e529bd694bb54f35fc07 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Olli=20V=C3=A4in=C3=B6l=C3=A4?= Date: Fri, 18 Dec 2015 10:04:16 +0200 Subject: [PATCH] Adding plastic material with linear elasticity --- src/elasticity.jl | 1 - src/elasticplastic.jl | 16 +++------- src/linear_elasticity.jl | 55 ++++++++++++++++++++++++++++++++++ src/vonmises.jl | 8 +++-- test/test_linear_elasticity.jl | 50 +++++++++++++++++++++++++++++-- 5 files changed, 112 insertions(+), 18 deletions(-) diff --git a/src/elasticity.jl b/src/elasticity.jl index 8275fb0..646c346 100644 --- a/src/elasticity.jl +++ b/src/elasticity.jl @@ -81,7 +81,6 @@ function get_residual_vector{P<:ElasticityProblem}(problem::Problem{P}, element: # strain E = 1/2*(F'*F - I) - # stress S = lambda*trace(E)*I + 2*mu*E diff --git a/src/elasticplastic.jl b/src/elasticplastic.jl index a919df0..c75df3b 100644 --- a/src/elasticplastic.jl +++ b/src/elasticplastic.jl @@ -58,16 +58,14 @@ function get_residual_vector{P<:PlaneStressElasticPlasticProblem}(problem::Probl F = I + gradu E = 1/2*(F'*F - I) + #E = 1/2*(gradu + gradu') # finite strain (total) + # 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 + de_v = [dstrain[1,1], dstrain[2,2], dstrain[1,2]] material_model = element("material model", time) s = last_stress de = copy(ForwardDiff.get_value(dstrain)) @@ -90,14 +88,8 @@ function get_residual_vector{P<:PlaneStressElasticPlasticProblem}(problem::Probl 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 + s_v += C * (de_v - dep) info("--: ", ForwardDiff.get_value(s_v)) # stress if P == PlaneStressElasticPlasticProblem diff --git a/src/linear_elasticity.jl b/src/linear_elasticity.jl index 96d9d14..a57245b 100644 --- a/src/linear_elasticity.jl +++ b/src/linear_elasticity.jl @@ -45,6 +45,13 @@ function assemble!{E<:CG, P<:LinearElasticityProblem}(assembly::Assembly, proble B[6, 3*(i-1)+1] = dN[3,i] B[6, 3*(i-1)+3] = dN[1,i] end + # L = b * B' + # D = 0.5 * (L' + L) + # F = ... + # E = 0.5 * (F'*F - I) + # de = E - E_last + # S = vonMisesStress(de, stress) + # K = B' * S * J * w add!(assembly.stiffness_matrix, gdofs, gdofs, w*B'*C*B*det(J)) end if haskey(element, "displacement load") @@ -114,3 +121,51 @@ function assemble!{E<:CG, P<:PlaneStressLinearElasticityProblem}(assembly::Assem end end +############################### +# Plastic material # +############################### +include("vonmises.jl") +abstract PlaneStressLinearElasticPlasticProblem <: LinearElasticityProblem + +function PlaneStressLinearElasticPlasticProblem(name="plane stress linear elasticity", dim::Int=2, elements=[]) + return Problem{PlaneStressLinearElasticPlasticProblem}(name, dim, elements) +end + +""" Elasticity equations, plane stress. """ +function assemble!{E<:CG, P<:PlaneStressLinearElasticPlasticProblem}(assembly::Assembly, problem::Problem{P}, element::Element{E}, time::Real) + + gdofs = get_gdofs(element, problem.dim) + ndim, nnodes = size(E) + B = zeros(3, 2*nnodes) + for ip in get_integration_points(element) + w = ip.weight + J = get_jacobian(element, ip, time) + N = element(ip, time) + if haskey(element, "youngs modulus") && haskey(element, "poissons ratio") + nu = element("poissons ratio", ip, time) + E_ = element("youngs modulus", ip, time) + C = E_/(1.0 - nu^2) .* [ + 1.0 nu 0.0 + nu 1.0 0.0 + 0.0 0.0 (1.0-nu)/2.0] + dN = element(ip, time, Val{:grad}) + fill!(B, 0.0) + for i=1:size(dN, 2) + B[1, 2*(i-1)+1] = dN[1,i] + B[2, 2*(i-1)+2] = dN[2,i] + B[3, 2*(i-1)+1] = dN[2,i] + B[3, 2*(i-1)+2] = dN[1,i] + end + add!(assembly.stiffness_matrix, gdofs, gdofs, w*B'*C*B*det(J)) + end + if haskey(element, "displacement load") + b = element("displacement load", ip, time) + add!(assembly.force_vector, gdofs, w*N'*b*det(J)) + end + if haskey(element, "displacement traction force") + T = element("displacement traction force", ip, time) + L = w*T*N*norm(J) + add!(assembly.force_vector, gdofs, vec(L)) + end + end +end diff --git a/src/vonmises.jl b/src/vonmises.jl index 653bd67..ecc477f 100644 --- a/src/vonmises.jl +++ b/src/vonmises.jl @@ -291,12 +291,14 @@ function calculate_stress(dstrain, stress, C, stress_y, df = ForwardDiff.jacobian(f) # Calculating root results = find_root!(f, df, x) - stress_tot = stress + results[1:3] + dstress = results[1:3] + stress_tot = stress + dstress plastic_multiplier = results[end] vm_wrap(stress_) = vonMisesYieldPlaneStress(stress_, stress_y) dfds = ForwardDiff.gradient(vm_wrap) - dep = plastic_multiplier * dfds(stress_tot) + dep = plastic_multiplier * dfds(vec(stress_tot)) info("II ", stress_tot) - return results[1:3], dep + info(vm_wrap(stress_tot)) + return dstress, dep end end diff --git a/test/test_linear_elasticity.jl b/test/test_linear_elasticity.jl index 3eb13f6..2047173 100644 --- a/test/test_linear_elasticity.jl +++ b/test/test_linear_elasticity.jl @@ -10,7 +10,7 @@ using JuliaFEM.Core: Seg2, Quad4, Hex8, LinearElasticityProblem, get_connectivit assemble, PlaneStressLinearElasticityProblem, DirichletProblem, LinearSolver using JuliaFEM.Preprocess: aster_parse_nodes - +using JuliaFEM.Core: PlaneStressLinearElasticPlasticProblem function test_plane_stress_linear_elasticity_with_surface_load() nodes = Dict{Int64, Vector{Float64}}( @@ -54,8 +54,54 @@ function test_plane_stress_linear_elasticity_with_surface_load() # 2015-10-22-plane-stress/cplan_linear_traction_force.* @test isapprox(u[:,3], [2.77777777777778E-03, -1.11111111111111E-02]) end + #test_plane_stress_linear_elasticity_with_surface_load() +function test_plane_stress_linear_elasticplastic_with_surface_load() + nodes = Dict{Int64, Vector{Float64}}( + 1 => [0.0, 0.0], + 2 => [1.0, 0.0], + 3 => [1.0, 1.0], + 4 => [0.0, 1.0]) + + function set_geometry!(element, nodes) + element["geometry"] = Vector{Float64}[nodes[i] for i in get_connectivity(element)] + end + element1 = Quad4([1, 2, 3, 4]) + set_geometry!(element1, nodes) + element1["youngs modulus"] = 9000.0 + element1["poissons ratio"] = 0.25 + + element2 = Seg2([3, 4]) + set_geometry!(element2, nodes) + element2["displacement traction force"] = Vector{Float64}[[0.0, -100.0] for i=1:2] + + # problem = PlaneStressLinearElasticityProblem() + problem = PlaneStressLinearElasticPlasticProblem() + push!(problem, element1) + push!(problem, element2) + + free_dofs = Int64[3, 5, 6, 8] + + ass = assemble(problem, 0.0) + f = full(ass.force_vector) + K = full(ass.stiffness_matrix) +# info("initial force vector") +# dump(reshape(f, 2, 4)) +# info("initial stiffness matrix") +# dump(round(Int, K)[free_dofs, free_dofs]) + + u = zeros(2, 4) + u[free_dofs] = K[free_dofs, free_dofs] \ f[free_dofs] + + info("result vector") + dump(u) + # verified using Code Aster. + # 2015-10-22-plane-stress/cplan_linear_traction_force.* + @test isapprox(u[:,3], [2.77777777777778E-03, -1.11111111111111E-02]) +end +# test_plane_stress_linear_elasticplastic_with_surface_load() + function test_continuum_elasticity_with_surface_load() nodes = JuliaFEM.Preprocess.aster_parse_nodes(""" @@ -141,7 +187,7 @@ function test_continuum_elasticity_with_surface_load() # [1/36, 1/36, -1/9] @test isapprox(u, [2.77777777777778E-02, 2.77777777777778E-02, -1.11111111111111E-01]) end -#test_continuum_elasticity_with_surface_load() +# test_continuum_elasticity_with_surface_load() end