From 04f3c0e7053d06c0a55c5ff08d7234798e3fa335 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Thu, 11 Feb 2016 16:16:15 +0200 Subject: [PATCH] Geometrically nonlinear formulation for 3d. --- src/elasticity.jl | 170 ++++++++++++++++++--------- test/test_elasticity_surface_load.jl | 104 +++++----------- 2 files changed, 145 insertions(+), 129 deletions(-) diff --git a/src/elasticity.jl b/src/elasticity.jl index 915dd88..9e82178 100644 --- a/src/elasticity.jl +++ b/src/elasticity.jl @@ -29,14 +29,14 @@ end function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time::Real) props = problem.properties + gdofs = get_gdofs(problem, element) if props.formulation == :continuum - return assemble!(assembly, problem, element, time, Val{:continuum}) + Kt, f = assemble(problem, element, time, Val{:continuum}) elseif (props.formulation == :plane_stress) || (props.formulation == :plane_strain) - gdofs = get_gdofs(problem, element) Kt, f = assemble(problem, element, time, Val{:plane}) - add!(assembly.K, gdofs, gdofs, Kt) - add!(assembly.f, gdofs, f) end + add!(assembly.K, gdofs, gdofs, Kt) + add!(assembly.f, gdofs, f) end @@ -163,71 +163,129 @@ end """ Elasticity equations, continuum formulation. """ -function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time::Real, ::Type{Val{:continuum}}) +function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum}}) + + props = problem.properties + dim = get_unknown_field_dimension(problem) + nnodes = size(element, 2) + BL = zeros(6, dim*nnodes) + BNL = zeros(9, dim*nnodes) + Kt = zeros(dim*nnodes, dim*nnodes) + f = zeros(dim*nnodes) - gdofs = get_gdofs(problem, element) - ndim, nnodes = size(element) - B = zeros(6, 3*nnodes) for ip in get_integration_points(element) - w = ip.weight J = get_jacobian(element, ip, time) + w = ip.weight*det(J) N = element(ip, time) - if haskey(element, "youngs modulus") && haskey(element, "poissons ratio") - v = element("poissons ratio", ip, time) - E_ = element("youngs modulus", ip, time) - a = 1 - v - b = 1 - 2*v - c = 1 + v - C = E_/(b*c) .* [ - a v v 0 0 0 - v a v 0 0 0 - v v a 0 0 0 - 0 0 0 b 0 0 - 0 0 0 0 b 0 - 0 0 0 0 0 b] - dN = element(ip, time, Val{:grad}) - fill!(B, 0.0) - for i=1:size(dN, 2) - B[1, 3*(i-1)+1] = dN[1,i] - B[2, 3*(i-1)+2] = dN[2,i] - B[3, 3*(i-1)+3] = dN[3,i] - B[4, 3*(i-1)+1] = dN[2,i] - B[4, 3*(i-1)+2] = dN[1,i] - B[5, 3*(i-1)+2] = dN[3,i] - B[5, 3*(i-1)+3] = dN[2,i] - 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 - Kt = w*B'*C*B*det(J) - add!(assembly.K, gdofs, gdofs, Kt) + dN = element(ip, time, Val{:grad}) + + # kinematics; calculate deformation gradient and strain + F = eye(dim) + if haskey(element, "displacement") + gradu = element("displacement", ip, time, Val{:grad}) + F += gradu end + GL = 1/2*(F'*F - I) # green-lagrange strain + + E = element("youngs modulus", ip, time) + nu = element("poissons ratio", ip, time) + + a = 1 - nu + b = 1 - 2*nu + c = 1 + nu + D = E/(b*c) .* [ + a nu nu 0 0 0 + nu a nu 0 0 0 + nu nu a 0 0 0 + 0 0 0 b 0 0 + 0 0 0 0 b 0 + 0 0 0 0 0 b] + + # # PK2 stress tensor in voigt notation + S = D*[GL[1,1]; GL[2,2]; GL[3,3]; 2*GL[2,3]; 2*GL[1,3]; 2*GL[1,2]] + + # add contributions: material and geometric stiffness + internal forces + fill!(BL, 0.0) + for i=1:size(dN, 2) + BL[1, 3*(i-1)+1] = F[1,1]*dN[1,i] + BL[1, 3*(i-1)+2] = F[2,1]*dN[1,i] + BL[1, 3*(i-1)+3] = F[3,1]*dN[1,i] + BL[2, 3*(i-1)+1] = F[1,2]*dN[2,i] + BL[2, 3*(i-1)+2] = F[2,2]*dN[2,i] + BL[2, 3*(i-1)+3] = F[3,2]*dN[2,i] + BL[3, 3*(i-1)+1] = F[1,3]*dN[3,i] + BL[3, 3*(i-1)+2] = F[2,3]*dN[3,i] + BL[3, 3*(i-1)+3] = F[3,3]*dN[3,i] + BL[4, 3*(i-1)+1] = F[1,1]*dN[2,i] + F[1,2]*dN[1,i] + BL[4, 3*(i-1)+2] = F[2,1]*dN[2,i] + F[2,2]*dN[1,i] + BL[4, 3*(i-1)+3] = F[3,1]*dN[2,i] + F[3,2]*dN[1,i] + BL[5, 3*(i-1)+1] = F[1,2]*dN[3,i] + F[1,3]*dN[2,i] + BL[5, 3*(i-1)+2] = F[2,2]*dN[3,i] + F[2,3]*dN[2,i] + BL[5, 3*(i-1)+3] = F[3,2]*dN[3,i] + F[3,3]*dN[2,i] + BL[6, 3*(i-1)+1] = F[1,3]*dN[1,i] + F[1,1]*dN[3,i] + BL[6, 3*(i-1)+2] = F[2,3]*dN[1,i] + F[2,1]*dN[3,i] + BL[6, 3*(i-1)+3] = F[3,3]*dN[1,i] + F[3,1]*dN[3,i] + end + fill!(BNL, 0.0) + for i=1:size(dN, 2) + BNL[1, 3*(i-1)+1] = dN[1,i] + BNL[2, 3*(i-1)+1] = dN[2,i] + BNL[3, 3*(i-1)+1] = dN[3,i] + BNL[4, 3*(i-1)+2] = dN[1,i] + BNL[5, 3*(i-1)+2] = dN[2,i] + BNL[6, 3*(i-1)+2] = dN[3,i] + BNL[7, 3*(i-1)+3] = dN[1,i] + BNL[8, 3*(i-1)+3] = dN[2,i] + BNL[9, 3*(i-1)+3] = dN[3,i] + end + S3 = zeros(3*dim, 3*dim) + S3[1,1] = S[1] + S3[2,2] = S[2] + S3[3,3] = S[3] + S3[2,3] = S3[3,2] = S[4] + S3[1,3] = S3[3,1] = S[5] + S3[1,2] = S3[2,1] = S[6] + S3[4:6,4:6] = S3[7:9,7:9] = S3[1:3,1:3] + + Kt += w*(BL'*D*BL + BNL'*S3*BNL) + f -= w*BL'*S + + # volume load if haskey(element, "displacement load") - b = element("displacement load", ip, time) - add!(assembly.f, gdofs, w*N'*b*det(J)) + T = element("displacement load", ip, time) + f += vec(w*T*N) end + + end + + return Kt, f +end + +""" Elasticity equations, surface traction for continuum formulation. """ +function assemble{El<:Union{Tri3, Tri6, Quad4}}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum}}) + + props = problem.properties + dim = get_unknown_field_dimension(problem) + nnodes = size(element, 2) + Kt = zeros(dim*nnodes, dim*nnodes) + f = zeros(dim*nnodes) + + for ip in get_integration_points(element) + JT = transpose(get_jacobian(element, ip, time)) + N = element(ip, time) + w = ip.weight*norm(cross(JT[:,1], JT[:,2])) if haskey(element, "displacement traction force") T = element("displacement traction force", ip, time) - JT = transpose(J) - L = w*T*N*norm(cross(JT[:,1], JT[:,2])) - add!(assembly.f, gdofs, vec(L)) + f += vec(w*T*N) end - for dim in 1:get_unknown_field_dimension(problem) - if haskey(element, "displacement traction force $dim") - T = element("displacement traction force $dim", ip, time) - ldofs = gdofs[dim:unknown_field_dimension(problem):end] - JT = transpose(J) - L = w*T*N*norm(cross(JT[:,1], JT[:,2])) - add!(assembly.f, ldofs, vec(L)) + for i in 1:dim + if haskey(element, "displacement traction force $i") + T = element("displacement traction force $i", ip, time) + f[i:dim:end] += vec(w*T*N) end end end + return Kt, f end diff --git a/test/test_elasticity_surface_load.jl b/test/test_elasticity_surface_load.jl index 9d0059a..959379c 100644 --- a/test/test_elasticity_surface_load.jl +++ b/test/test_elasticity_surface_load.jl @@ -2,7 +2,7 @@ # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md using JuliaFEM.Test -using JuliaFEM.Core: Node, update!, Quad4, Seg2, Problem, Elasticity, Solver, Dirichlet +using JuliaFEM.Core: Node, update!, Quad4, Seg2, Hex8, Problem, Elasticity, Solver, Dirichlet using JuliaFEM.Preprocess: aster_parse_nodes @testset "test 2d linear elasticity with surface load" begin @@ -40,7 +40,7 @@ using JuliaFEM.Preprocess: aster_parse_nodes # type, name, dimension, unknown_field_name boundary_problem = Problem(Dirichlet, "symmetry boundaries", 2, "displacement") push!(boundary_problem, sym13, sym23) - + solver = Solver("solve block problem") solver.is_linear_system = true # to get linear solution push!(solver, elasticity_problem) @@ -95,87 +95,45 @@ end @testset "test continuum linear elasticity with surface load" begin - nodes = JuliaFEM.Preprocess.aster_parse_nodes(""" - COOR_3D - N1 0.0 0.0 0.0 - N2 1.0 0.0 0.0 - N3 1.0 1.0 0.0 - N4 0.0 1.0 0.0 - N5 0.0 0.0 1.0 - N6 1.0 0.0 1.0 - N7 1.0 1.0 1.0 - N8 0.0 1.0 1.0 - FINSF - """) + nodes = Dict{Int64, Node}( + 1 => [0.0, 0.0, 0.0], + 2 => [1.0, 0.0, 0.0], + 3 => [1.0, 1.0, 0.0], + 4 => [0.0, 1.0, 0.0], + 5 => [0.0, 0.0, 1.0], + 6 => [1.0, 0.0, 1.0], + 7 => [1.0, 1.0, 1.0], + 8 => [0.0, 1.0, 1.0]) - function set_geometry!(element, nodes) - element["geometry"] = Vector{Float64}[nodes[i] for i in get_connectivity(element)] - end element1 = Hex8([1, 2, 3, 4, 5, 6, 7, 8]) - set_geometry!(element1, nodes) - element1["youngs modulus"] = 900.0 - element1["poissons ratio"] = 0.25 - element1["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0, 0.0] for i=1:8]) - element2 = Quad4([5, 6, 7, 8]) - set_geometry!(element2, nodes) - element2["displacement traction force"] = Vector{Float64}[[0.0, 0.0, -100.0] for i=1:4] - element2["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0, 0.0] for i=1:4]) + update!([element1, element2], "geometry", nodes) + update!([element1], "youngs modulus", 900.0) + update!([element1], "poissons ratio", 0.25) + update!([element2], "displacement traction force", Vector{Float64}[[0.0, 0.0, -100.0] for i=1:4]) - problem = ElasticityProblem() - push!(problem, element1) - push!(problem, element2) + elasticity_problem = Problem(Elasticity, "solve continuum block", 3) + push!(elasticity_problem, element1) + push!(elasticity_problem, element2) - free_dofs = zeros(Bool, 8, 3) - x = 1 - y = 2 - z = 3 - free_dofs[2, x] = true - free_dofs[3, [x, y]] = true - free_dofs[4, y] = true - free_dofs[5, z] = true - free_dofs[6, [x, z]] = true - free_dofs[7, [x, y, z]] = true - free_dofs[8, [y, z]] = true - free_dofs = find(vec(free_dofs')) - info("free dofs: $free_dofs") + symxy = Quad4([1, 2, 3, 4]) + symxz = Quad4([1, 2, 6, 5]) + symyz = Quad4([1, 4, 8, 5]) + update!([symxy, symxz, symyz], "geometry", nodes) + symxy["displacement 3"] = 0.0 + symxz["displacement 2"] = 0.0 + symyz["displacement 1"] = 0.0 + boundary_problem = Problem(Dirichlet, "symmetry boundary conditions", 3, "displacement") + push!(boundary_problem, symxy, symxz, symyz) - info("initial force vector") - ass = JuliaFEM.Core.assemble(problem, 0.0) - info(reshape(full(ass.force_vector), 3, 8)) - info("initial stiffness matrix") - dump(round(Int, full(ass.stiffness_matrix))[free_dofs, free_dofs]) - solve!(problem, free_dofs, 0.0; max_iterations=10) - -#= - dx = Quad4([1, 4, 8, 5]) - dx["displacement 1"] = 0.0 - dy = Quad4([1, 5, 6, 2]) - dy["displacement 2"] = 0.0 - dz = Quad4([1, 2, 3, 4]) - dz["displacement 3"] = 0.0 - bc = DirichletProblem("displacement", 3) - for el in [dx, dy, dz] - set_geometry!(el, nodes) - push!(bc, el) - end - solver = JuliaFEM.Core.DirectSolver() - push!(solver, problem) - push!(solver, bc) - solver.dump_matrices = true - solver.name = "3d_hex8" - solver(0.0) -=# + solver = Solver("solve 3d block") + push!(solver, elasticity_problem) + push!(solver, boundary_problem) + call(solver) disp = element1("displacement", [1.0, 1.0, 1.0], 0.0) info("displacement at tip: $disp") - info("displacement on element: ") - for (i, d) in enumerate(element1("displacement", 0.0)) - @printf "%d % f % f % f\n" [i;d]... - end # verified using Code Aster. # 2015-12-12-continuum-elasticity/vim c3d_grot_gdep_traction_force.comm @test isapprox(disp, [3.17431158889468E-02, 3.17431158889468E-02, -1.38591518927826E-01]) - #@test isapprox(disp, [2.80559539222183E-03, 2.80559539222183E-03, -1.13019918093242E-02]) end -