mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-24 19:19:54 +00:00
Geometrically nonlinear formulation for 3d.
This commit is contained in:
Jukka Aho
+114
-56
@@ -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
|
||||
|
||||
|
||||
|
||||
@@ -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
|
||||
|
||||
|
||||
Reference in New Issue
Block a user