Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-07 19:56:23 +00:00
Code Issues Packages Projects Releases Wiki Activity

test stress extrapolation to nodes

Browse Source
This commit is contained in:
Jukka Aho
2016-12-01 11:14:32 +02:00
parent b3a14c99f4
commit e4a65c0e4e
1 changed files with 29 additions and 12 deletions
Download Patch File Download Diff File Expand all files Collapse all files
test/test_elasticity_hollow_sphere_with_surface_pressure.jl
+29 -12
View File
@@ -71,9 +71,11 @@ test subjects:
- surface pressure load in curved surface
- verification of elements wedge6 and wedge15
from Code Aster:
N38 -8.85861895037377E-01 -3.46944695195361E-18 -3.46944695195361E-18
from Code Aster, for linear model
N38 -8.85861895037377E-01 -3.46944695195361E-18 -3.46944695195361E-18
coords of N38 = (1.0, 0.0, 0.0)
N49 -4.50934684240566E-01 -4.46908405333021E-01 -5.92329377847994E-01
coords of N49 = (0.521002, 0.515482, 0.68032) # quite middle of surface
Analytical solution
uᵣ(r) = b³p/(2Er²(a³-b³)) * (a³(ν+1) + r³(-4ν+2)), where
@@ -84,7 +86,7 @@ if a=0.9, b=1.0, ν=1/3, E = 24580 and p = 7317 equation yields
http://mms2.ensmp.fr/emms_paris/plasticite3D/exercices/eSpherePress.pdf
=#
function test_wedge_sphere(model)
function test_wedge_sphere(model, u_CA, S_CA)
mesh_file = Pkg.dir("JuliaFEM") * "/test/testdata/primitives.med"
mesh = aster_read_mesh(mesh_file, model)
body = Problem(Elasticity, "hollow sphere 1/8 model", 3)
@@ -114,17 +116,25 @@ function test_wedge_sphere(model)
rtol = norm(um - 0.9) / max(norm(um), 0.9) * 100.0
info("mean umag = $um, std umag = $us, rtol = $rtol")
@test rtol < 1.5 # percents
# for linear model
#=
u_CA = [-8.85861895037377E-01, -3.46944695195361E-18, -3.46944695195361E-18]
rtol = norm(u[38] - u_CA) / max(norm(u[38]), norm(u_CA)) * 100.0
info("rel diff to CA = $rtol %")
@test isapprox(u[38], u_CA)
=#
info("Verifying displacement against Code Aster solution.. ")
for nid in keys(u_CA)
rtol = norm(u[nid] - u_CA[nid]) / max(norm(u[nid]), norm(u_CA[nid])) * 100.0
info("Node id $nid, rel diff to CA = $rtol %")
@test isapprox(u[nid], u_CA[nid])
end
S = lsq_fit(body.elements, get_stress, 0.0)
info("Verifying stress against Code Aster solution.. ")
for nid in keys(S_CA)
rtol = norm(S[nid] - S_CA[nid]) / max(norm(S[nid]), norm(S_CA[nid])) * 100.0
info("Node id $nid, S=$(S[nid]), S_CA=$(S_CA[nid]), rel diff to CA = $rtol %")
#@test isapprox(S[nid], S_CA[nid])
# http://code-aster.org/doc/default/en/man_r/r3/r3.06.03.pdf
@test rtol < 10.0 # percents
end
# Calculate principal stresses in nodes
Sp = Dict()
for (nid, s) in S
@@ -152,10 +162,17 @@ function test_wedge_sphere(model)
principal_stresses = sort(eigvals(stress_tensor))
println("princial stresses = $principal_stresses")
end
end
@testset """1/8 hollow sphere with surface load""" begin
test_wedge_sphere("HOLLOWSPHERE8_WEDGE6")
u_CA = Dict()
u_CA[38] = [-8.85861895037377E-01, -3.46944695195361E-18, -3.46944695195361E-18]
u_CA[49] = [-4.50934684240566E-01, -4.46908405333021E-01, -5.92329377847994E-01]
S_CA = Dict()
S_CA[38] = [-1.94504819940510E+03, -3.47014655438479E+04, -3.40876135857114E+04, 1.70379805227866E+03, 2.16707698388864E+03, 4.32009983342141E-12]
S_CA[49] = [-2.36067010168718E+04, -2.33820775692753E+04, -1.80606705820104E+04, 8.85783922092890E+03, 1.12122431934777E+04, 1.14035796957343E+04]
test_wedge_sphere("HOLLOWSPHERE8_WEDGE6", u_CA, S_CA)
#test_wedge_sphere("HOLLOWSPHERE8_WEDGE15")
end