Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-22 02:40:51 +00:00
Code Issues Packages Projects Releases Wiki Activity

elasticity tests

Browse Source
This commit is contained in:
Jukka Aho
2016-05-24 08:01:41 +03:00
parent 2aa0362812
commit f7a0fd23c8
6 changed files with 229 additions and 21 deletions
Download Patch File Download Diff File Expand all files Collapse all files
test/test_elasticity_2d_linear_with_surface_load.jl
+23 -15
View File
@@ -8,15 +8,22 @@ using JuliaFEM.Test
@testset "test 2d linear elasticity with surface load" begin
meshfile = "/geometry/2d_block/BLOCK_1elem.med"
mesh = parse_aster_med_file(Pkg.dir("JuliaFEM")*meshfile)
# field problem
body = Problem(Elasticity, "BLOCK", 2)
body.properties.formulation = :plane_stress
body_elements = aster_create_elements(mesh, :BLOCK, :QU4)
update!(body_elements, "youngs modulus", 900.0)
update!(body_elements, "poissons ratio", 0.25)
trac_elements = aster_create_elements(mesh, :TOP, :SE2)
update!(trac_elements, "displacement traction force 2", -100.0)
push!(body, body_elements..., trac_elements...)
block = Problem(Elasticity, "BLOCK", 2)
block.properties.formulation = :plane_stress
block.properties.finite_strain = false
elements = aster_create_elements(mesh, :BLOCK, :QU4)
update!(elements, "youngs modulus", 288.0)
update!(elements, "poissons ratio", 1/3)
update!(elements, "displacement load 2", 576.0)
push!(block, elements...)
traction = aster_create_elements(mesh, :TOP, :SE2)
update!(traction, "displacement traction force 2", 288.0)
push!(block, traction...)
# boundary conditions
bc_sym = Problem(Dirichlet, "symmetry bc", 2, "displacement")
bc_elements_left = aster_create_elements(mesh, :LEFT, :SE2)
@@ -24,14 +31,15 @@ using JuliaFEM.Test
update!(bc_elements_left, "displacement 1", 0.0)
update!(bc_elements_bottom, "displacement 2", 0.0)
push!(bc_sym, bc_elements_left..., bc_elements_bottom...)
solver = Solver("solve block problem")
body.properties.finite_strain = false
push!(solver, body, bc_sym)
push!(solver, block, bc_sym)
call(solver)
f = -100.0
E = 900.0
nu = 0.25
u3_expected = f/E*[-nu, 1]
u3 = reshape(body.assembly.u, 2, 4)[:,3]
f = 288.0
g = 576.0
E = 288.0
nu = 1/3
u3_expected = f/E*[-nu, 1] + g/(2*E)*[-nu, 1]
u3 = reshape(block.assembly.u, 2, 4)[:,3]
@test isapprox(u3, u3_expected)
end
test/test_elasticity_3d_linear_with_surface_load.jl
+7 -6
View File
@@ -18,9 +18,10 @@ using JuliaFEM.Test
element1 = Element(Hex8, [1, 2, 3, 4, 5, 6, 7, 8])
element2 = Element(Quad4, [5, 6, 7, 8])
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])
update!([element1], "youngs modulus", 288.0)
update!([element1], "poissons ratio", 1/3)
update!([element2], "displacement traction force 3", 288.0)
update!([element1], "displacement load 3", 576.0)
elasticity_problem = Problem(Elasticity, "solve continuum block", 3)
elasticity_problem.properties.finite_strain = false
@@ -31,9 +32,9 @@ using JuliaFEM.Test
symxz = Element(Quad4, [1, 2, 6, 5])
symyz = Element(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
symxz["displacement 2"] = 0.0
symxy["displacement 3"] = 0.0
boundary_problem = Problem(Dirichlet, "symmetry boundary conditions", 3, "displacement")
push!(boundary_problem, symxy, symxz, symyz)
@@ -44,6 +45,6 @@ using JuliaFEM.Test
disp = element1("displacement", [1.0, 1.0, 1.0], 0.0)
info("displacement at tip: $disp")
u_expected = -100.0/900.0 * [-0.25, -0.25, 1.0]
u_expected = 2.0 * [-1/3, -1/3, 1.0]
@test isapprox(disp, u_expected)
end
test/test_elasticity_3d_unit_block.jl
+99
View File
@@ -0,0 +1,99 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Test
function get_model(fn, vol, sur)
meshfile = Pkg.dir("JuliaFEM")*"/geometry/3d_blocks/BLOCK.med"
mesh = parse_aster_med_file(meshfile, fn)
block = Problem(Elasticity, fn, 3)
block.properties.finite_strain = false
elements = aster_create_elements(mesh, :BLOCK, vol)
update!(elements, "youngs modulus", 288.0)
update!(elements, "poissons ratio", 1/3)
update!(elements, "displacement load 3", 576.0)
push!(block, elements...)
traction = aster_create_elements(mesh, :LOAD, sur)
update!(traction, "displacement traction force 3", 288.0)
push!(block, traction...)
bc = Problem(Dirichlet, "symmetry boundary condition", 3, "displacement")
# bc.properties.formulation = :incremental
symyz = aster_create_elements(mesh, :SYMYZ, sur)
symxz = aster_create_elements(mesh, :SYMXZ, sur)
symxy = aster_create_elements(mesh, :SYMXY, sur)
update!(symyz, "displacement 1", 0.0)
update!(symxz, "displacement 2", 0.0)
update!(symxy, "displacement 3", 0.0)
push!(bc, symyz..., symxz..., symxy...)
return block, bc, elements, traction, symyz, symxz, symxy
end
function calc_size(elements, dim)
A = 0.0
for element in elements
Ael = 0.0
size(element, 1) == dim || continue
for (w, xi) in get_integration_points(element)
detJ = element(xi, 0.0, Val{:detJ})
Ael += w*detJ
end
for (i, X) in enumerate(element["geometry"](0.0))
info("$i : $X")
end
info("Area / volume: $Ael")
A += Ael
end
return A
end
#=
@testset "test 3d block hex8" begin
block, bc = get_model("BLOCK_HEX8", :HE8, :QU4)
V = calc_size(block.elements, 3)
info("volume of block: $V")
A = calc_size(bc.elements, 2)
info("area of boundary condition: $A")
@test isapprox(V, 1.0)
@test isapprox(A, 3.0)
solver = Solver("solver block problem")
solver.is_linear_system = true
push!(solver, block, bc)
call(solver)
max_u = maximum(abs(block.assembly.u))
info("max |u| = $max_u")
@test isapprox(max_u, 2.0)
end
=#
@testset "test 3d block TET4" begin
block, bc, elements, traction, symyz, symxz, symxy = get_model("BLOCK_TET4", :TE4, :TR3)
# block, bc = get_model("BLOCK_TET10", :T10, :TR6)
V = calc_size(block.elements, 3)
info("volume of block: $V")
A = calc_size(bc.elements, 2)
info("area of boundary condition: $A")
At = calc_size(traction, 2)
info("area of load surface: $At")
@test isapprox(V, 1.0)
@test isapprox(At, 1.0)
@test isapprox(A, 3.0)
solver = Solver("solver block problem")
#solver.is_linear_system = true
push!(solver, block, bc)
call(solver)
max_u = maximum(block.assembly.u)
nu = round(Int, length(block.assembly.u)/3)
u = reshape(block.assembly.u, 3, nu)
f = reshape(full(block.assembly.f), 3, nu)
dump(round(u', 5))
dump(round(f', 5))
info("max |u| = $max_u")
@test isapprox(max_u, 2.0)
end
test/test_elasticity_tet10_stiffness_matrix.jl
+33
View File
@@ -0,0 +1,33 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Test
@testset "test tet10 stiffness matrix" begin
el = Element(Tet10)
el["youngs modulus"] = 480.0
el["poissons ratio"] = 1/3
x1 = [2.0, 3.0, 4.0]
x2 = [6.0, 3.0, 2.0]
x3 = [2.0, 5.0, 1.0]
x4 = [4.0, 3.0, 6.0]
x5 = 0.5*(x1+x2)
x6 = 0.5*(x2+x3)
x7 = 0.5*(x3+x1)
x8 = 0.5*(x1+x4)
x9 = 0.5*(x2+x4)
x10 = 0.5*(x3+x4)
el["geometry"] = Vector{Float64}[x1, x2, x3, x4, x5, x6, x7, x8, x9, x10]
pr = Problem(Elasticity, "tet10", 3)
Kt, f = assemble(pr, el, 0.0)
eigs = real(eigvals(Kt))
eigs_expected = [8809.45, 4936.01, 2880.56, 2491.66, 2004.85,
1632.49, 1264.32, 1212.42, 817.905,
745.755, 651.034, 517.441, 255.1, 210.955,
195.832, 104.008, 72.7562, 64.4376, 53.8515,
23.8417, 16.6354, 9.54682, 6.93361, 2.22099,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
@test isapprox(eigs, eigs_expected; atol=1.0e-2)
end
test/test_elasticity_tet4_stiffness_matrix.jl
+33
View File
@@ -0,0 +1,33 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Test
@testset "test tet4 stiffness matrix" begin
el = Element(Tet4)
el["youngs modulus"] = 96.0
el["poissons ratio"] = 1/3
x1 = [2.0, 3.0, 4.0]
x2 = [6.0, 3.0, 2.0]
x3 = [2.0, 5.0, 1.0]
x4 = [4.0, 3.0, 6.0]
el["geometry"] = Vector{Float64}[x1, x2, x3, x4]
pr = Problem(Elasticity, "tet4", 3)
Kt, f = assemble(pr, el, 0.0)
Kt_expected = [
149.0 108.0 24.0 -1.0 6.0 12.0 -54.0 -48.0 0.0 -94.0 -66.0 -36.0
108.0 344.0 54.0 -24.0 104.0 42.0 -24.0 -216.0 -12.0 -60.0 -232.0 -84.0
24.0 54.0 113.0 0.0 30.0 35.0 0.0 -24.0 -54.0 -24.0 -60.0 -94.0
-1.0 -24.0 0.0 29.0 -18.0 -12.0 -18.0 24.0 0.0 -10.0 18.0 12.0
6.0 104.0 30.0 -18.0 44.0 18.0 12.0 -72.0 -12.0 0.0 -76.0 -36.0
12.0 42.0 35.0 -12.0 18.0 29.0 0.0 -24.0 -18.0 0.0 -36.0 -46.0
-54.0 -24.0 0.0 -18.0 12.0 0.0 36.0 0.0 0.0 36.0 12.0 0.0
-48.0 -216.0 -24.0 24.0 -72.0 -24.0 0.0 144.0 0.0 24.0 144.0 48.0
0.0 -12.0 -54.0 0.0 -12.0 -18.0 0.0 0.0 36.0 0.0 24.0 36.0
-94.0 -60.0 -24.0 -10.0 0.0 0.0 36.0 24.0 0.0 68.0 36.0 24.0
-66.0 -232.0 -60.0 18.0 -76.0 -36.0 12.0 144.0 24.0 36.0 164.0 72.0
-36.0 -84.0 -94.0 12.0 -36.0 -46.0 0.0 48.0 36.0 24.0 72.0 104.0]
@test isapprox(Kt, Kt_expected)
end
test/test_elasticity_tet4_volume_load.jl
+34
View File
@@ -0,0 +1,34 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Test
@testset "test tet4 element under volume load" begin
X1 = [2.0, 3.0, 4.0]
X2 = [6.0, 3.0, 2.0]
X3 = [2.0, 5.0, 1.0]
X4 = [4.0, 3.0, 6.0]
e1 = Element(Tet4, [1, 2, 3, 4])
e1["geometry"] = Node[X1, X2, X3, X4]
e1["youngs modulus"] = 96.0
e1["poissons ratio"] = 1/3
e1["displacement load 3"] = 784.0/110.0
e2 = Element(Tri3, [1, 2, 3])
e2["geometry"] = Node[X2, X1, X3]
e2["displacement 1"] = 0.0
e2["displacement 2"] = 0.0
e2["displacement 3"] = 0.0
p1 = Problem(Elasticity, "tetra", 3)
p2 = Problem(Dirichlet, "bc", 3, "displacement")
p1.properties.finite_strain = false
push!(p1, e1)
push!(p2, e2)
s = Solver()
push!(s, p1, p2)
call(s)
u_4 = p1.assembly.u[10:end]
u_expected = [-3.0/220.0, -9.0/220.0, 1.0/10.0]
@test isapprox(u_4, u_expected)
end