diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index 0aecff0..9a543d3 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -66,8 +66,8 @@ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::El add!(assembly.f, gdofs, f) end -const Elasticity2DSurfaceElements = Union{Poi1,Seg2,Seg3} -const Elasticity2DVolumeElements = Union{Tri3,Tri6,Quad4,Quad8,Quad9} +include("problems_elasticity_2d.jl") + const Elasticity3DSurfaceElements = Union{Poi1,Tri3,Tri6,Quad4,Quad8,Quad9} const Elasticity3DVolumeElements = Union{Tet4, Pyr5, Wedge6, Wedge15, Hex8, Tet10, Hex20, Hex27} @@ -100,210 +100,6 @@ function get_keys(element) map(x -> all_keys[x], idx) end -#function initialize_internal_params!(params, ip_id, ::Type{Val{:type_3d}}) -# if !(ip_id in keys(params)) -# params[ip_id] = Dict{Any, Any}() -# params[ip_id]["last_stress"] = [0.0,0.0,0.0,0.0,0.0,0.0] -# params[ip_id]["last_strain"] = [0.0,0.0,0.0,0.0,0.0,0.0] -# end -#end - -""" Elasticity equations for 2d cases. """ -function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time, ::Type{Val{:plane}}) - - props = problem.properties - dim = get_unknown_field_dimension(problem) - nnodes = length(element) - BL = zeros(3, dim*nnodes) - BNL = zeros(4, dim*nnodes) - Km = zeros(dim*nnodes, dim*nnodes) - Kg = zeros(dim*nnodes, dim*nnodes) - f = zeros(dim*nnodes) - Dtan = zeros(3,3) - - for ip in get_integration_points(element) - - detJ = element(ip, time, Val{:detJ}) - w = ip.weight*detJ - N = element(ip, time) - dN = element(ip, time, Val{:Grad}) - # kinematics - - gradu = element("displacement", ip, time, Val{:Grad}) - fill!(BL, 0.0) - - if props.finite_strain - strain = 1/2*(gradu + gradu' + gradu'*gradu) - F = eye(dim) + gradu - for i=1:size(dN, 2) - BL[1, 2*(i-1)+1] += F[1,1]*dN[1,i] - BL[1, 2*(i-1)+2] += F[2,1]*dN[1,i] - BL[2, 2*(i-1)+1] += F[1,2]*dN[2,i] - BL[2, 2*(i-1)+2] += F[2,2]*dN[2,i] - BL[3, 2*(i-1)+1] += F[1,1]*dN[2,i] + F[1,2]*dN[1,i] - BL[3, 2*(i-1)+2] += F[2,1]*dN[2,i] + F[2,2]*dN[1,i] - end - else # linearized strain - strain = 1/2*(gradu + gradu') - F = eye(dim) - for i=1:size(dN, 2) - BL[1, 2*(i-1)+1] = dN[1,i] - BL[2, 2*(i-1)+2] = dN[2,i] - BL[3, 2*(i-1)+1] = dN[2,i] - BL[3, 2*(i-1)+2] = dN[1,i] - end - end - - strain_vec = [strain[1,1]; strain[2,2]; strain[1,2]] - - # calculate stress - E = element("youngs modulus", ip, time) - nu = element("poissons ratio", ip, time) - if props.formulation == :plane_stress - D = E/(1.0 - nu^2) .* [ - 1.0 nu 0.0 - nu 1.0 0.0 - 0.0 0.0 (1.0-nu)/2.0] - elseif props.formulation == :plane_strain - D = E/((1.0+nu)*(1.0-2.0*nu)) .* [ - 1.0-nu nu 0.0 - nu 1.0-nu 0.0 - 0.0 0.0 (1.0-2.0*nu)/2.0] - else - error("unknown plane formulation: $(props.formulation)") - end - - # calculate stress - element_keys = get_keys(element) - - if "plasticity" in element_keys - plastic_def = element("plasticity")[ip.id] - - calculate_stress! = plastic_def["type"] - yield_surface_ = plastic_def["yield_surface"] - params = plastic_def["params"] - - initialize_internal_params!(params, ip, Val{:type_2d}) - - if time == 0.0 - error("Given step time = $(time). Please select time > 0.0") - end - - t_last = ip("prev_time", time) - update!(ip, "prev_time", time => t_last) - - dt = time - t_last - - stress_last = ip("stress", t_last) - strain_last = ip("strain", t_last) - - dstrain_vec = strain_vec - strain_last - stress_vec = [0.0, 0.0, 0.0] - pstrain = zeros(3) - calculate_stress!(stress_vec, stress_last, dstrain_vec, pstrain, D, params, Dtan, yield_surface_, time, dt, Val{:type_2d}) - else - stress_vec = D * ([1.0, 1.0, 2.0] .* strain_vec) - Dtan[:,:] = D[:,:] - end - - :strain in props.store_fields && update!(ip, "strain", time => strain_vec) - :stress in props.store_fields && update!(ip, "stress", time => stress_vec) - :stress11 in props.store_fields && update!(ip, "stress11", time => stress_vec[1]) - :stress22 in props.store_fields && update!(ip, "stress22", time => stress_vec[2]) - :stress12 in props.store_fields && update!(ip, "stress12", time => stress_vec[3]) - - Km += w*BL'*Dtan*BL - - - # stress = [stress_vec[1] stress_vec[3]; stress_vec[3] stress_vec[2]] - # cauchy_stress = F'*stress*F/det(F) - # cauchy_stress = [cauchy_stress[1,1]; cauchy_stress[2,2]; cauchy_stress[1,2]] - # update!(ip, "cauchy stress", time => cauchy_stress) - - # material stiffness end - - if props.geometric_stiffness - # take geometric stiffness into account - - fill!(BNL, 0.0) - for i=1:size(dN, 2) - BNL[1, 2*(i-1)+1] = dN[1,i] - BNL[2, 2*(i-1)+1] = dN[2,i] - BNL[3, 2*(i-1)+2] = dN[1,i] - BNL[4, 2*(i-1)+2] = dN[2,i] - end - - S2 = zeros(2*dim, 2*dim) - S2[1,1] = stress_vec[1] - S2[2,2] = stress_vec[2] - S2[1,2] = S2[2,1] = stress_vec[3] - S2[3:4,3:4] = S2[1:2,1:2] - - Kg += w*BNL'*S2*BNL # geometric stiffness - - end - - # rhs, internal and external load - - f -= w*BL'*stress_vec - - if haskey(element, "displacement load") - b = element("displacement load", ip, time) - f += w*vec(b*N) - end - - for i=1:dim - if haskey(element, "displacement load $i") - b = element("displacement load $i", ip, time) - f[i:dim:end] += w*vec(b*N) - end - end - - end - - return Km, Kg, f -end - -function assemble{El<:Elasticity2DSurfaceElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:plane}}) - - props = problem.properties - dim = get_unknown_field_dimension(problem) - nnodes = length(element) - Km = zeros(dim*nnodes, dim*nnodes) - Kg = zeros(dim*nnodes, dim*nnodes) - f = zeros(dim*nnodes) - - for ip in get_integration_points(element) - - detJ = element(ip, time, Val{:detJ}) - w = ip.weight*detJ - N = element(ip, time) - - if haskey(element, "displacement traction force") - T = element("displacement traction force", ip, time) - f += w*vec(T*N) - end - - for i=1:dim - # traction force for ith component - if haskey(element, "displacement traction force $i") - T = element("displacement traction force $i", ip, time) - f[i:dim:end] += w*vec(T*N) - end - end - - if haskey(element, "nt displacement traction force") - # traction force given in normal-tangential direction - T = element("nt displacement traction force", ip, time) - Q = element("normal-tangential coordinates", ip, time) - f += w*vec(Q'*T*N) - end - - end - - return Km, Kg, f -end - """ Elasticity equations, 3d nonlinear. """ function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum}}) props = problem.properties diff --git a/src/problems_elasticity_2d.jl b/src/problems_elasticity_2d.jl new file mode 100644 index 0000000..3eb8559 --- /dev/null +++ b/src/problems_elasticity_2d.jl @@ -0,0 +1,219 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +const Elasticity2DSurfaceElements = Union{Poi1,Seg2,Seg3} +const Elasticity2DVolumeElements = Union{Tri3,Tri6,Quad4,Quad8,Quad9} + +function assemble!(assembly::Assembly, problem::Problem{Elasticity}, elements::Vector{Element}, time::Real, ::Type{Val{:plane_stress}}) + for element in elements + gdofs = get_gdofs(problem, element) + Km, Kg, f = assemble(problem, element, time, Val{:plane}) + add!(assembly.K, gdofs, gdofs, Km) + add!(assembly.Kg, gdofs, gdofs, Kg) + add!(assembly.f, gdofs, f) + end +end + +function assemble!(assembly::Assembly, problem::Problem{Elasticity}, elements::Vector{Element}, time::Real, ::Type{Val{:plane_strain}}) + for element in elements + gdofs = get_gdofs(problem, element) + Km, Kg, f = assemble(problem, element, time, Val{:plane}) + add!(assembly.K, gdofs, gdofs, Km) + add!(assembly.Kg, gdofs, gdofs, Kg) + add!(assembly.f, gdofs, f) + end +end + +""" Plane elasticity equations (plane stress, plane strain). """ +function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time, ::Type{Val{:plane}}) + + props = problem.properties + dim = get_unknown_field_dimension(problem) + nnodes = length(element) + BL = zeros(3, dim*nnodes) + BNL = zeros(4, dim*nnodes) + Km = zeros(dim*nnodes, dim*nnodes) + Kg = zeros(dim*nnodes, dim*nnodes) + f = zeros(dim*nnodes) + Dtan = zeros(3,3) + + for ip in get_integration_points(element) + + detJ = element(ip, time, Val{:detJ}) + w = ip.weight*detJ + N = element(ip, time) + dN = element(ip, time, Val{:Grad}) + # kinematics + + gradu = element("displacement", ip, time, Val{:Grad}) + fill!(BL, 0.0) + + if props.finite_strain + strain = 1/2*(gradu + gradu' + gradu'*gradu) + F = eye(dim) + gradu + for i=1:size(dN, 2) + BL[1, 2*(i-1)+1] += F[1,1]*dN[1,i] + BL[1, 2*(i-1)+2] += F[2,1]*dN[1,i] + BL[2, 2*(i-1)+1] += F[1,2]*dN[2,i] + BL[2, 2*(i-1)+2] += F[2,2]*dN[2,i] + BL[3, 2*(i-1)+1] += F[1,1]*dN[2,i] + F[1,2]*dN[1,i] + BL[3, 2*(i-1)+2] += F[2,1]*dN[2,i] + F[2,2]*dN[1,i] + end + else # linearized strain + strain = 1/2*(gradu + gradu') + F = eye(dim) + for i=1:size(dN, 2) + BL[1, 2*(i-1)+1] = dN[1,i] + BL[2, 2*(i-1)+2] = dN[2,i] + BL[3, 2*(i-1)+1] = dN[2,i] + BL[3, 2*(i-1)+2] = dN[1,i] + end + end + + strain_vec = [strain[1,1]; strain[2,2]; strain[1,2]] + + # calculate stress + E = element("youngs modulus", ip, time) + nu = element("poissons ratio", ip, time) + if props.formulation == :plane_stress + D = E/(1.0 - nu^2) .* [ + 1.0 nu 0.0 + nu 1.0 0.0 + 0.0 0.0 (1.0-nu)/2.0] + elseif props.formulation == :plane_strain + D = E/((1.0+nu)*(1.0-2.0*nu)) .* [ + 1.0-nu nu 0.0 + nu 1.0-nu 0.0 + 0.0 0.0 (1.0-2.0*nu)/2.0] + else + error("unknown plane formulation: $(props.formulation)") + end + + # calculate stress + + if haskey(element, "plasticity") + plastic_def = element("plasticity")[ip.id] + + calculate_stress! = plastic_def["type"] + yield_surface_ = plastic_def["yield_surface"] + params = plastic_def["params"] + + initialize_internal_params!(params, ip, Val{:type_2d}) + + if time == 0.0 + error("Given step time = $(time). Please select time > 0.0") + end + + t_last = ip("prev_time", time) + update!(ip, "prev_time", time => t_last) + dt = time - t_last + + stress_last = ip("stress", t_last) + strain_last = ip("strain", t_last) + + dstrain_vec = strain_vec - strain_last + stress_vec = [0.0, 0.0, 0.0] + pstrain = zeros(3) + calculate_stress!(stress_vec, stress_last, dstrain_vec, pstrain, D, params, Dtan, yield_surface_, time, dt, Val{:type_2d}) + else + stress_vec = D * ([1.0, 1.0, 2.0] .* strain_vec) + Dtan[:,:] = D[:,:] + end + + :strain in props.store_fields && update!(ip, "strain", time => strain_vec) + :stress in props.store_fields && update!(ip, "stress", time => stress_vec) + :stress11 in props.store_fields && update!(ip, "stress11", time => stress_vec[1]) + :stress22 in props.store_fields && update!(ip, "stress22", time => stress_vec[2]) + :stress12 in props.store_fields && update!(ip, "stress12", time => stress_vec[3]) + + Km += w*BL'*Dtan*BL + + + # stress = [stress_vec[1] stress_vec[3]; stress_vec[3] stress_vec[2]] + # cauchy_stress = F'*stress*F/det(F) + # cauchy_stress = [cauchy_stress[1,1]; cauchy_stress[2,2]; cauchy_stress[1,2]] + # update!(ip, "cauchy stress", time => cauchy_stress) + + # material stiffness end + + if props.geometric_stiffness + # take geometric stiffness into account + + fill!(BNL, 0.0) + for i=1:size(dN, 2) + BNL[1, 2*(i-1)+1] = dN[1,i] + BNL[2, 2*(i-1)+1] = dN[2,i] + BNL[3, 2*(i-1)+2] = dN[1,i] + BNL[4, 2*(i-1)+2] = dN[2,i] + end + + S2 = zeros(2*dim, 2*dim) + S2[1,1] = stress_vec[1] + S2[2,2] = stress_vec[2] + S2[1,2] = S2[2,1] = stress_vec[3] + S2[3:4,3:4] = S2[1:2,1:2] + + Kg += w*BNL'*S2*BNL # geometric stiffness + + end + + # rhs, internal and external load + + f -= w*BL'*stress_vec + + if haskey(element, "displacement load") + b = element("displacement load", ip, time) + f += w*vec(b*N) + end + + for i=1:dim + if haskey(element, "displacement load $i") + b = element("displacement load $i", ip, time) + f[i:dim:end] += w*vec(b*N) + end + end + + end + + return Km, Kg, f +end + +function assemble{El<:Elasticity2DSurfaceElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:plane}}) + + props = problem.properties + dim = get_unknown_field_dimension(problem) + nnodes = length(element) + Km = zeros(dim*nnodes, dim*nnodes) + Kg = zeros(dim*nnodes, dim*nnodes) + f = zeros(dim*nnodes) + + for ip in get_integration_points(element) + + detJ = element(ip, time, Val{:detJ}) + w = ip.weight*detJ + N = element(ip, time) + + if haskey(element, "displacement traction force") + T = element("displacement traction force", ip, time) + f += w*vec(T*N) + end + + for i=1:dim + # traction force for ith component + if haskey(element, "displacement traction force $i") + T = element("displacement traction force $i", ip, time) + f[i:dim:end] += w*vec(T*N) + end + end + + if haskey(element, "nt displacement traction force") + # traction force given in normal-tangential direction + T = element("nt displacement traction force", ip, time) + Q = element("normal-tangential coordinates", ip, time) + f += w*vec(Q'*T*N) + end + + end + + return Km, Kg, f +end