mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
Use package HeatTransfer.jl for heat problems (#194)
Heat transfer analysis is moved to its own package where the development continues. Two small modifications are needed for test files: - Instead of `problem.properties.formulation`, we have two separate problems, `PlaneHeat` for two-dimensional problems and `Heat` for three-dimensional problems. - Unnecessary prefixing of field names is changed. For example, now we simply have only "thermal conductivity" and not prefixed "temperature thermal conductivity".
This commit is contained in:
Jukka Aho
+3
-2
@@ -33,14 +33,15 @@ end
|
||||
|
||||
using AbaqusReader
|
||||
using AsterReader
|
||||
|
||||
@reexport using HeatTransfer
|
||||
include("problems_elasticity.jl")
|
||||
export Elasticity
|
||||
include("materials_plasticity.jl")
|
||||
export plastic_von_mises
|
||||
include("problems_dirichlet.jl")
|
||||
export Dirichlet
|
||||
include("problems_heat.jl")
|
||||
export Heat
|
||||
|
||||
export assemble!, postprocess!
|
||||
### Mortar methods ###
|
||||
include("problems_mortar.jl")
|
||||
|
||||
@@ -1,224 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
""" Heat equations.
|
||||
|
||||
Field equation is:
|
||||
|
||||
ρc∂u/∂t = ∇⋅(k∇u) + f
|
||||
|
||||
Weak form is: find u∈U such that ∀v in V
|
||||
|
||||
∫k∇u∇v dx = ∫fv dx + ∫gv ds,
|
||||
|
||||
where
|
||||
|
||||
k = temperature thermal conductivity defined on volume elements
|
||||
f = temperature load defined on volume elements
|
||||
g = temperature flux defined on boundary elements
|
||||
|
||||
Parameters
|
||||
----------
|
||||
temperature thermal conductivity
|
||||
temperature load
|
||||
temperature flux
|
||||
thermal conductivity
|
||||
heat source
|
||||
heat flux
|
||||
heat transfer coefficient
|
||||
external temperature
|
||||
|
||||
|
||||
Formulations
|
||||
------------
|
||||
1D, 2D, 3D
|
||||
|
||||
References
|
||||
----------
|
||||
https://en.wikipedia.org/wiki/Heat_equation
|
||||
https://en.wikipedia.org/wiki/Heat_capacity
|
||||
https://en.wikipedia.org/wiki/Heat_flux
|
||||
https://en.wikipedia.org/wiki/Thermal_conduction
|
||||
https://en.wikipedia.org/wiki/Thermal_conductivity
|
||||
https://en.wikipedia.org/wiki/Thermal_diffusivity
|
||||
https://en.wikipedia.org/wiki/Volumetric_heat_capacity
|
||||
"""
|
||||
type Heat <: FieldProblem
|
||||
formulation :: AbstractString
|
||||
store_fields :: Vector{Symbol}
|
||||
end
|
||||
|
||||
function Heat()
|
||||
return Heat("3D", [])
|
||||
end
|
||||
|
||||
function get_unknown_field_name(problem::Problem{Heat})
|
||||
return "temperature"
|
||||
end
|
||||
|
||||
function assemble!(assembly::Assembly, problem::Problem{Heat}, element::Element, time::Float64)
|
||||
formulation = Val{Symbol(problem.properties.formulation)}
|
||||
assemble!(assembly, problem, element, time, formulation)
|
||||
end
|
||||
|
||||
# 3d heat problems
|
||||
|
||||
function assemble!{E}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("3D")}})
|
||||
info("Unknown element type $E for 3d heat problem!")
|
||||
end
|
||||
|
||||
const Heat3DVolumeElements = Union{Tet4, Tet10, Pyr5, Hex8, Hex20, Hex27}
|
||||
const Heat3DSurfaceElements = Union{Tri3,Tri6,Quad4,Quad8,Quad9}
|
||||
|
||||
function assemble!{E<:Heat3DVolumeElements}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("3D")}})
|
||||
gdofs = get_gdofs(problem, element)
|
||||
field_name = get_unknown_field_name(problem)
|
||||
nnodes = length(element)
|
||||
K = zeros(nnodes, nnodes)
|
||||
fq = zeros(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, "$field_name thermal conductivity")
|
||||
dN = element(ip, time, Val{:Grad})
|
||||
k = element("$field_name thermal conductivity", ip, time)
|
||||
K += w*k*dN'*dN
|
||||
end
|
||||
if haskey(element, "thermal conductivity")
|
||||
dN = element(ip, time, Val{:Grad})
|
||||
k = element("thermal conductivity", ip, time)
|
||||
K += w*k*dN'*dN
|
||||
end
|
||||
if haskey(element, "$field_name load")
|
||||
f = element("$field_name load", ip, time)
|
||||
fq += w*N'*f
|
||||
end
|
||||
if haskey(element, "heat source")
|
||||
f = element("heat source", ip, time)
|
||||
fq += w*N'*f
|
||||
end
|
||||
end
|
||||
T = [interpolate(element[field_name], time)...]
|
||||
fq -= K*T
|
||||
add!(assembly.K, gdofs, gdofs, K)
|
||||
add!(assembly.f, gdofs, fq)
|
||||
end
|
||||
|
||||
function assemble!{E<:Heat3DSurfaceElements}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("3D")}})
|
||||
gdofs = get_gdofs(problem, element)
|
||||
field_name = get_unknown_field_name(problem)
|
||||
nnodes = length(element)
|
||||
K = zeros(nnodes, nnodes)
|
||||
fq = zeros(nnodes)
|
||||
for ip in get_integration_points(element, 2)
|
||||
detJ = element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ
|
||||
N = element(ip, time)
|
||||
if haskey(element, "$field_name flux")
|
||||
q = element("$field_name flux", ip, time)
|
||||
fq += w*N'*q
|
||||
end
|
||||
if haskey(element, "heat flux")
|
||||
q = element("heat flux", ip, time)
|
||||
fq += w*N'*q
|
||||
end
|
||||
if haskey(element, "$field_name heat transfer coefficient") && haskey(element, "$field_name external temperature")
|
||||
h = element("$field_name heat transfer coefficient", ip, time)
|
||||
Tu = element("$field_name external temperature", ip, time)
|
||||
K += w*h*N'*N
|
||||
fq += w*N'*h*Tu
|
||||
end
|
||||
if haskey(element, "heat transfer coefficient") && haskey(element, "external temperature")
|
||||
h = element("heat transfer coefficient", ip, time)
|
||||
Tu = element("external temperature", ip, time)
|
||||
K += w*h*N'*N
|
||||
fq += w*N'*h*Tu
|
||||
end
|
||||
end
|
||||
T = collect(element(field_name, time))
|
||||
fq -= K*T
|
||||
add!(assembly.K, gdofs, gdofs, K)
|
||||
add!(assembly.f, gdofs, fq)
|
||||
end
|
||||
|
||||
# 2d heat problems
|
||||
|
||||
function assemble!{E}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("2D")}})
|
||||
info("Unknown element type $E for 2d heat problem!")
|
||||
end
|
||||
|
||||
const Heat2DVolumeElements = Union{Tri3,Tri6,Quad4}
|
||||
const Heat2DSurfaceElements = Union{Seg2,Seg3}
|
||||
|
||||
function assemble!{E<:Heat2DVolumeElements}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("2D")}})
|
||||
gdofs = get_gdofs(problem, element)
|
||||
field_name = get_unknown_field_name(problem)
|
||||
nnodes = length(element)
|
||||
K = zeros(nnodes, nnodes)
|
||||
fq = zeros(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, "$field_name thermal conductivity")
|
||||
dN = element(ip, time, Val{:Grad})
|
||||
k = element("$field_name thermal conductivity", ip, time)
|
||||
K += w*k*dN'*dN
|
||||
end
|
||||
if haskey(element, "thermal conductivity")
|
||||
dN = element(ip, time, Val{:Grad})
|
||||
k = element("thermal conductivity", ip, time)
|
||||
K += w*k*dN'*dN
|
||||
end
|
||||
if haskey(element, "$field_name load")
|
||||
f = element("$field_name load", ip, time)
|
||||
fq += w*N'*f
|
||||
end
|
||||
if haskey(element, "heat source")
|
||||
f = element("heat source", ip, time)
|
||||
fq += w*N'*f
|
||||
end
|
||||
end
|
||||
T = collect(element(field_name, time))
|
||||
fq -= K*T
|
||||
add!(assembly.K, gdofs, gdofs, K)
|
||||
add!(assembly.f, gdofs, fq)
|
||||
end
|
||||
|
||||
function assemble!{E<:Heat2DSurfaceElements}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("2D")}})
|
||||
gdofs = get_gdofs(problem, element)
|
||||
field_name = get_unknown_field_name(problem)
|
||||
nnodes = length(element)
|
||||
K = zeros(nnodes, nnodes)
|
||||
fq = zeros(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, "$field_name flux")
|
||||
g = element("$field_name flux", ip, time)
|
||||
fq += w*N'*g
|
||||
end
|
||||
if haskey(element, "heat flux")
|
||||
g = element("heat flux", ip, time)
|
||||
fq += w*N'*g
|
||||
end
|
||||
if haskey(element, "$field_name heat transfer coefficient") && haskey(element, "$field_name external temperature")
|
||||
h = element("$field_name heat transfer coefficient", ip, time)
|
||||
Tu = element("$field_name external temperature", ip, time)
|
||||
K += w*h*N'*N
|
||||
fq += w*N'*h*Tu
|
||||
end
|
||||
if haskey(element, "heat transfer coefficient") && haskey(element, "external temperature")
|
||||
h = element("heat transfer coefficient", ip, time)
|
||||
Tu = element("external temperature", ip, time)
|
||||
K += w*h*N'*N
|
||||
fq += w*N'*h*Tu
|
||||
end
|
||||
end
|
||||
T = collect(element(field_name, time))
|
||||
fq -= K*T
|
||||
add!(assembly.K, gdofs, gdofs, K)
|
||||
add!(assembly.f, gdofs, fq)
|
||||
end
|
||||
+4
-2
@@ -15,8 +15,10 @@ const to = TimerOutput()
|
||||
for fn in test_files
|
||||
info("----- Running tests from file $fn -----")
|
||||
t0 = time()
|
||||
timeit(to, fn) do
|
||||
include(fn)
|
||||
@testset "$fn" begin
|
||||
timeit(to, fn) do
|
||||
include(fn)
|
||||
end
|
||||
end
|
||||
dt = round(time() - t0, 2)
|
||||
info("----- Testing file $fn completed in $dt seconds -----")
|
||||
|
||||
@@ -5,43 +5,39 @@ using JuliaFEM
|
||||
using JuliaFEM: assemble_mass_matrix!, add_elements!
|
||||
using Base.Test
|
||||
|
||||
@testset "test tet10 mass matrix" begin
|
||||
X = Dict(
|
||||
1 => [2.0, 3.0, 4.0],
|
||||
2 => [6.0, 3.0, 2.0],
|
||||
3 => [2.0, 5.0, 1.0],
|
||||
4 => [4.0, 3.0, 6.0])
|
||||
X[5] = 1/2*(X[1] + X[2])
|
||||
X[6] = 1/2*(X[2] + X[3])
|
||||
X[7] = 1/2*(X[3] + X[1])
|
||||
X[8] = 1/2*(X[1] + X[4])
|
||||
X[9] = 1/2*(X[2] + X[4])
|
||||
X[10] = 1/2*(X[3] + X[4])
|
||||
X = Dict(
|
||||
1 => [2.0, 3.0, 4.0],
|
||||
2 => [6.0, 3.0, 2.0],
|
||||
3 => [2.0, 5.0, 1.0],
|
||||
4 => [4.0, 3.0, 6.0])
|
||||
X[5] = 1/2*(X[1] + X[2])
|
||||
X[6] = 1/2*(X[2] + X[3])
|
||||
X[7] = 1/2*(X[3] + X[1])
|
||||
X[8] = 1/2*(X[1] + X[4])
|
||||
X[9] = 1/2*(X[2] + X[4])
|
||||
X[10] = 1/2*(X[3] + X[4])
|
||||
|
||||
element = Element(Tet10, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
|
||||
update!(element, "youngs modulus", 480.0)
|
||||
update!(element, "poissons ratio", 1/3)
|
||||
update!(element, "geometry", X)
|
||||
update!(element, "density", 105.0)
|
||||
element = Element(Tet10, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
|
||||
update!(element, "youngs modulus", 480.0)
|
||||
update!(element, "poissons ratio", 1/3)
|
||||
update!(element, "geometry", X)
|
||||
update!(element, "density", 105.0)
|
||||
|
||||
body = Problem(Heat, "TET", 1)
|
||||
add_elements!(body, [element])
|
||||
assemble!(body, 0.0)
|
||||
assemble_mass_matrix!(body, 0.0)
|
||||
#println(body.assembly.K)
|
||||
M = full(body.assembly.M)
|
||||
body = Problem(Heat, "TET", 1)
|
||||
add_elements!(body, [element])
|
||||
assemble_mass_matrix!(body, 0.0)
|
||||
M = full(body.assembly.M)
|
||||
|
||||
M_expected = [
|
||||
6 1 1 1 -4 -6 -4 -4 -6 -6
|
||||
1 6 1 1 -4 -4 -6 -6 -4 -6
|
||||
1 1 6 1 -6 -4 -4 -6 -6 -4
|
||||
1 1 1 6 -6 -6 -6 -4 -4 -4
|
||||
-4 -4 -6 -6 32 16 16 16 16 8
|
||||
-6 -4 -4 -6 16 32 16 8 16 16
|
||||
-4 -6 -4 -6 16 16 32 16 8 16
|
||||
-4 -6 -6 -4 16 8 16 32 16 16
|
||||
-6 -4 -6 -4 16 16 8 16 32 16
|
||||
-6 -6 -4 -4 8 16 16 16 16 32]
|
||||
M_expected = [
|
||||
6 1 1 1 -4 -6 -4 -4 -6 -6
|
||||
1 6 1 1 -4 -4 -6 -6 -4 -6
|
||||
1 1 6 1 -6 -4 -4 -6 -6 -4
|
||||
1 1 1 6 -6 -6 -6 -4 -4 -4
|
||||
-4 -4 -6 -6 32 16 16 16 16 8
|
||||
-6 -4 -4 -6 16 32 16 8 16 16
|
||||
-4 -6 -4 -6 16 16 32 16 8 16
|
||||
-4 -6 -6 -4 16 8 16 32 16 16
|
||||
-6 -4 -6 -4 16 16 8 16 32 16
|
||||
-6 -6 -4 -4 8 16 16 16 16 32]
|
||||
|
||||
@test isapprox(M, M_expected)
|
||||
end
|
||||
@test isapprox(M, M_expected)
|
||||
|
||||
+46
-98
@@ -14,18 +14,18 @@ using JuliaFEM.Postprocess
|
||||
face = Problem(Heat, "face 4", 1)
|
||||
fixed = Problem(Dirichlet, "fixed face 3", 1, "temperature")
|
||||
prob.elements = create_elements(mesh, "TET")
|
||||
update!(prob, "temperature thermal conductivity", 50.0)
|
||||
face.elements = create_elements(mesh, "FACE4")
|
||||
update!(face, "temperature external temperature", 20.0)
|
||||
update!(face, "temperature heat transfer coefficient", 60.0)
|
||||
update!(prob, "thermal conductivity", 50.0)
|
||||
face.elements = create_elements(mesh, "FACE4")
|
||||
update!(face, "external temperature", 20.0)
|
||||
update!(face, "heat transfer coefficient", 60.0)
|
||||
fixed.elements = create_elements(mesh, "FACE2")
|
||||
info("# of elements in fixed set: $(length(fixed))")
|
||||
update!(fixed, "temperature 1", 0.0)
|
||||
solver = LinearSolver(prob, face, fixed)
|
||||
solver()
|
||||
T = prob.assembly.u
|
||||
info("Solution: $T")
|
||||
T_expected = [ # using code aster
|
||||
Temp = prob.assembly.u
|
||||
info("Solution: $Temp")
|
||||
Temp_expected = [ # using code aster
|
||||
1.45606533688540E+01
|
||||
0.0
|
||||
0.0
|
||||
@@ -36,10 +36,10 @@ using JuliaFEM.Postprocess
|
||||
1.05228712963739E+01
|
||||
0.0
|
||||
0.0]
|
||||
info("Expected: $T_expected")
|
||||
rtol = norm(T-T_expected)/max(norm(T), norm(T_expected))
|
||||
info("Expected: $Temp_expected")
|
||||
rtol = norm(Temp-Temp_expected)/max(norm(Temp), norm(Temp_expected))
|
||||
info("rtol = $rtol")
|
||||
@test isapprox(T, T_expected; rtol=1.0e-6)
|
||||
@test isapprox(Temp, Temp_expected; rtol=1.0e-6)
|
||||
end
|
||||
|
||||
@testset "2d heat problem (one element)" begin
|
||||
@@ -54,19 +54,18 @@ end
|
||||
el1 = Element(Quad4, [1, 2, 3, 4])
|
||||
|
||||
update!(el1, "geometry", X)
|
||||
update!(el1, "temperature thermal conductivity", 6.0)
|
||||
update!(el1, "temperature load", 12.0)
|
||||
update!(el1, "thermal conductivity", 6.0)
|
||||
update!(el1, "heat source", 12.0)
|
||||
|
||||
# define boundary element for flux
|
||||
el2 = Element(Seg2, [1, 2])
|
||||
update!(el2, "geometry", X)
|
||||
# linear ramp from 0 -> 6 in time 0 -> 1
|
||||
update!(el2, "temperature flux", 0.0 => 0.0)
|
||||
update!(el2, "temperature flux", 1.0 => 6.0)
|
||||
update!(el2, "heat flux", 0.0 => 0.0)
|
||||
update!(el2, "heat flux", 1.0 => 6.0)
|
||||
|
||||
# define heat problem and push elements to problem
|
||||
problem = Problem(Heat, "one element heat problem", 1)
|
||||
problem.properties.formulation = "2D"
|
||||
problem = Problem(PlaneHeat, "one element heat problem", 1)
|
||||
push!(problem, el1, el2)
|
||||
|
||||
# Set constant source f=12 with k=6. Accurate solution is
|
||||
@@ -94,57 +93,6 @@ end
|
||||
@test isapprox(A[free_dofs, free_dofs] \ b[free_dofs], [2.0, 2.0])
|
||||
end
|
||||
|
||||
#=
|
||||
@testset "test 1d heat problem" begin
|
||||
|
||||
function T_acc(x)
|
||||
# accurate solution
|
||||
a = 0.01
|
||||
L = 0.20
|
||||
k = 50.0
|
||||
Tᵤ = 20.0
|
||||
h = 10.0
|
||||
P = 4*a
|
||||
A = a^2
|
||||
α = h
|
||||
β = sqrt((h*P)/(k*A))
|
||||
T̂ = 100.0
|
||||
C = [1.0 1.0; (α+k*β)*exp(β*L) (α-k*β)*exp(-β*L)] \ [T̂-Tᵤ, 0.0]
|
||||
return dot(C, [exp(β*x), exp(-β*x)]) + Tᵤ
|
||||
end
|
||||
|
||||
X = Dict{Int, Vector{Float64}}(
|
||||
1 => [0.0, 0.0, 0.0],
|
||||
2 => [0.1, 0.0, 0.0],
|
||||
3 => [0.2, 0.0, 0.0])
|
||||
e1 = Element(Seg2, [1, 2])
|
||||
e2 = Element(Seg2, [2, 3])
|
||||
e3 = Element(Poi1, [3])
|
||||
|
||||
p1 = Problem(Heat, "1d heat problem", 1)
|
||||
p1.properties.formulation = "1D"
|
||||
push!(p1, e1, e2, e3)
|
||||
update!(p1, "geometry", X)
|
||||
a = 0.010
|
||||
update!(p1, "cross-section area", a^2)
|
||||
update!(p1, "cross-section perimeter", 4*a)
|
||||
update!(p1, "temperature thermal conductivity", 50.0) # k [W/(m∘C)]
|
||||
update!(p1, "temperature heat transfer coefficient", 10.0) # h [W/(m²∘C)]
|
||||
update!(p1, "temperature external temperature", 20.0)
|
||||
|
||||
p2 = Problem(Dirichlet, "left boundary", 1, "temperature")
|
||||
e3 = Element(Poi1, [1])
|
||||
update!(e3, "geometry", X)
|
||||
update!(e3, "temperature 1", 100.0)
|
||||
push!(p2, e3)
|
||||
|
||||
solver = LinearSolver(p1, p2)
|
||||
solver()
|
||||
T_min = minimum(p1.assembly.u)
|
||||
@test isapprox(T_max, T_acc(0.2); rtol=4.5e-2)
|
||||
end
|
||||
=#
|
||||
|
||||
@testset "compare simple 3d heat problem to code aster solution" begin
|
||||
fn = @__DIR__() * "/testdata/rod_short.med"
|
||||
mesh = aster_read_mesh(fn, "Hex8")
|
||||
@@ -158,17 +106,17 @@ end
|
||||
face4 = create_elements(mesh, "FACE4")
|
||||
face5 = create_elements(mesh, "FACE5")
|
||||
face6 = create_elements(mesh, "FACE6")
|
||||
update!(rod, "temperature thermal conductivity", 50.0)
|
||||
update!(face2, "temperature external temperature", 20.0)
|
||||
update!(face2, "temperature heat transfer coefficient", 60.0)
|
||||
update!(face3, "temperature external temperature", 30.0)
|
||||
update!(face3, "temperature heat transfer coefficient", 50.0)
|
||||
update!(face4, "temperature external temperature", 40.0)
|
||||
update!(face4, "temperature heat transfer coefficient", 40.0)
|
||||
update!(face5, "temperature external temperature", 50.0)
|
||||
update!(face5, "temperature heat transfer coefficient", 30.0)
|
||||
update!(face6, "temperature external temperature", 60.0)
|
||||
update!(face6, "temperature heat transfer coefficient", 20.0)
|
||||
update!(rod, "thermal conductivity", 50.0)
|
||||
update!(face2, "external temperature", 20.0)
|
||||
update!(face2, "heat transfer coefficient", 60.0)
|
||||
update!(face3, "external temperature", 30.0)
|
||||
update!(face3, "heat transfer coefficient", 50.0)
|
||||
update!(face4, "external temperature", 40.0)
|
||||
update!(face4, "heat transfer coefficient", 40.0)
|
||||
update!(face5, "external temperature", 50.0)
|
||||
update!(face5, "heat transfer coefficient", 30.0)
|
||||
update!(face6, "external temperature", 60.0)
|
||||
update!(face6, "heat transfer coefficient", 20.0)
|
||||
push!(p1, rod, face2, face3, face4, face5, face6)
|
||||
|
||||
p2 = Problem(Dirichlet, "left support T=100", 1, "temperature")
|
||||
@@ -197,7 +145,7 @@ end
|
||||
6 => [1.74615686370955E+04, -3.73021966895897E+02, -1.38038931136833E+02],
|
||||
7 => [1.74479150854902E+04, -9.99509347888065E+01, -3.70268090662933E+01],
|
||||
8 => [1.74479150854901E+04, -3.73021966895874E+02, -1.38185932677561E+02])
|
||||
FLUX_NOEU = Dict{Int64, Vector{Float64}}(
|
||||
FLUX_NOEU = Dict{Int64, Vector{Float64}}(
|
||||
1 => [1.74596669275930E+04, 7.55555618070503E-11, 3.68594044175552E-12],
|
||||
2 => [1.74684148339734E+04, 1.10418341137120E-11, 3.48876483258209E-12],
|
||||
3 => [1.74360055518019E+04, 7.91828824731056E-11, 1.95399252334028E-13],
|
||||
@@ -207,10 +155,10 @@ end
|
||||
7 => [1.74360055518019E+04, -4.73227515822118E+02, -1.75280965396335E+02],
|
||||
8 => [1.74447696000717E+04, -4.72904678032179E+02, -1.75280965396335E+02])
|
||||
|
||||
T = p1("temperature", 0.0)
|
||||
Temp = p1("temperature", 0.0)
|
||||
|
||||
for j in sort(collect(keys(T)))
|
||||
T1 = T[j][1]
|
||||
for j in sort(collect(keys(Temp)))
|
||||
T1 = Temp[j][1]
|
||||
T2 = TEMP[j]
|
||||
rtol = norm(T1-T2)/max(T1,T2)*100.0
|
||||
@printf "node %i temp, JF: %e, CA: %e, rtol: %10.6f %%\n" j T1 T2 rtol
|
||||
@@ -227,9 +175,9 @@ end
|
||||
p2 = Problem(Dirichlet, "left support T=100", 1, "temperature")
|
||||
p1.elements = create_elements(mesh, "ROD", "FACE2")
|
||||
p2.elements = create_elements(mesh, "FACE1")
|
||||
update!(p1, "temperature thermal conductivity", 100.0)
|
||||
update!(p1, "temperature external temperature", 0.0)
|
||||
update!(p1, "temperature heat transfer coefficient", 1000.0)
|
||||
update!(p1, "thermal conductivity", 100.0)
|
||||
update!(p1, "external temperature", 0.0)
|
||||
update!(p1, "heat transfer coefficient", 1000.0)
|
||||
update!(p2, "temperature 1", 100.0)
|
||||
solver = LinearSolver(p1, p2)
|
||||
solver()
|
||||
@@ -264,17 +212,17 @@ end
|
||||
face4 = create_elements(mesh, "FACE4")
|
||||
face5 = create_elements(mesh, "FACE5")
|
||||
face6 = create_elements(mesh, "FACE6")
|
||||
update!(rod, "temperature thermal conductivity", 50.0)
|
||||
update!(face2, "temperature external temperature", 20.0)
|
||||
update!(face2, "temperature heat transfer coefficient", 60.0)
|
||||
update!(face3, "temperature external temperature", 30.0)
|
||||
update!(face3, "temperature heat transfer coefficient", 50.0)
|
||||
update!(face4, "temperature external temperature", 40.0)
|
||||
update!(face4, "temperature heat transfer coefficient", 40.0)
|
||||
update!(face5, "temperature external temperature", 50.0)
|
||||
update!(face5, "temperature heat transfer coefficient", 30.0)
|
||||
update!(face6, "temperature external temperature", 60.0)
|
||||
update!(face6, "temperature heat transfer coefficient", 20.0)
|
||||
update!(rod, "thermal conductivity", 50.0)
|
||||
update!(face2, "external temperature", 20.0)
|
||||
update!(face2, "heat transfer coefficient", 60.0)
|
||||
update!(face3, "external temperature", 30.0)
|
||||
update!(face3, "heat transfer coefficient", 50.0)
|
||||
update!(face4, "external temperature", 40.0)
|
||||
update!(face4, "heat transfer coefficient", 40.0)
|
||||
update!(face5, "external temperature", 50.0)
|
||||
update!(face5, "heat transfer coefficient", 30.0)
|
||||
update!(face6, "external temperature", 60.0)
|
||||
update!(face6, "heat transfer coefficient", 20.0)
|
||||
push!(p1, rod, face2, face3, face4, face5, face6)
|
||||
p2 = Problem(Dirichlet, "left support T=100", 1, "temperature")
|
||||
p2.elements = create_elements(mesh, "FACE1")
|
||||
@@ -294,8 +242,8 @@ end
|
||||
models = ["Tet4", "Hex8", "Hex20", "Hex27", "Tet10"]
|
||||
|
||||
for model in models
|
||||
T = calc_3d_heat_model(model)
|
||||
T_min = minimum(T)
|
||||
Temp = calc_3d_heat_model(model)
|
||||
T_min = minimum(Temp)
|
||||
T_ca = CA_sol[model]
|
||||
rtol = norm(T_min-T_ca)/max(T_min,T_ca)*100.0
|
||||
@printf "%-10s : T_min = % g, T_ca = % g, rtol = %g %%\n" model T_min T_ca rtol
|
||||
|
||||
+35
-38
@@ -6,45 +6,42 @@ using JuliaFEM.Preprocess
|
||||
using JuliaFEM.Postprocess
|
||||
using JuliaFEM.Testing
|
||||
|
||||
@testset "3d rod" begin
|
||||
mesh = aster_read_mesh(@__DIR__()*"/testdata/primitives.med", "CYLINDER_20_TET4")
|
||||
problem = Problem(Heat, "rod of length 20", 1)
|
||||
problem.elements = create_elements(mesh, "CYLINDER")
|
||||
update!(problem, "temperature thermal conductivity", 200.0)
|
||||
outer = Problem(Heat, "outer surface", 1)
|
||||
outer.elements = create_elements(mesh, "FACE2", "FACE3")
|
||||
update!(outer, "temperature external temperature", 20.0)
|
||||
update!(outer, "temperature heat transfer coefficient", 1.0)
|
||||
#midline = Problem(Heat, "midline of rod", 1)
|
||||
#midline.elements = create_elements(mesh, "INNER_LINE")
|
||||
boundary = Problem(Dirichlet, "homogeneous dirichlet boundary", 1, "temperature")
|
||||
boundary.elements = create_elements(mesh, "FACE1")
|
||||
update!(boundary, "temperature 1", 100.0)
|
||||
#solver = LinearSolver(problem, outer, boundary, midline)
|
||||
solver = LinearSolver(problem, outer, boundary)
|
||||
solver()
|
||||
mesh = aster_read_mesh(@__DIR__()*"/testdata/primitives.med", "CYLINDER_20_TET4")
|
||||
problem = Problem(Heat, "rod of length 20", 1)
|
||||
problem.elements = create_elements(mesh, "CYLINDER")
|
||||
update!(problem, "thermal conductivity", 200.0)
|
||||
outer = Problem(Heat, "outer surface", 1)
|
||||
outer.elements = create_elements(mesh, "FACE2", "FACE3")
|
||||
update!(outer, "external temperature", 20.0)
|
||||
update!(outer, "heat transfer coefficient", 1.0)
|
||||
#midline = Problem(Heat, "midline of rod", 1)
|
||||
#midline.elements = create_elements(mesh, "INNER_LINE")
|
||||
boundary = Problem(Dirichlet, "homogeneous dirichlet boundary", 1, "temperature")
|
||||
boundary.elements = create_elements(mesh, "FACE1")
|
||||
update!(boundary, "temperature 1", 100.0)
|
||||
#solver = LinearSolver(problem, outer, boundary, midline)
|
||||
solver = LinearSolver(problem, outer, boundary)
|
||||
solver()
|
||||
|
||||
L = 20
|
||||
k = 200.0
|
||||
Tu = 20.0
|
||||
h = 1.0
|
||||
P = 2*pi
|
||||
A = pi
|
||||
α = h
|
||||
β = sqrt((h*P)/(k*A))
|
||||
T0 = 100.0
|
||||
C = [1.0 1.0; (α+k*β)*exp(β*L) (α-k*β)*exp(-β*L)] \ [T0-Tu, 0]
|
||||
T(x) = dot(C, [exp(β*x), exp(-β*x)]) + Tu
|
||||
# Analytical solution
|
||||
L = 20
|
||||
k = 200.0
|
||||
Tu = 20.0
|
||||
h = 1.0
|
||||
P = 2*pi
|
||||
A = pi
|
||||
α = h
|
||||
β = sqrt((h*P)/(k*A))
|
||||
T0 = 100.0
|
||||
C = [1.0 1.0; (α+k*β)*exp(β*L) (α-k*β)*exp(-β*L)] \ [T0-Tu, 0]
|
||||
|
||||
T_diff = []
|
||||
for x in linspace(0, 20)
|
||||
T_FEM = problem("temperature", [x, 0.0, 0.0])[1]
|
||||
T_ACC = T(x)
|
||||
push!(T_diff, norm(T_FEM - T_ACC))
|
||||
info("x = $x, T_FEM = $T_FEM, T_ACC = $T_ACC")
|
||||
end
|
||||
info("mean diff = ", mean(T_diff))
|
||||
# mean diff = 1.14
|
||||
@test mean(T_diff) < 1.2
|
||||
T_diff = []
|
||||
for x in linspace(0, 20)
|
||||
T_FEM = problem("temperature", [x, 0.0, 0.0])[1]
|
||||
T_ACC = dot(C, [exp(β*x), exp(-β*x)]) + Tu
|
||||
push!(T_diff, norm(T_FEM - T_ACC))
|
||||
info("x = $x, T_FEM = $T_FEM, T_ACC = $T_ACC")
|
||||
end
|
||||
info("mean diff = ", mean(T_diff))
|
||||
|
||||
@test mean(T_diff) < 1.2 # mean diff = 1.14
|
||||
|
||||
+1
-2
@@ -13,9 +13,8 @@ using JuliaFEM.Testing
|
||||
mesh_file = @__DIR__()*"/testdata/primitives.med"
|
||||
mesh = aster_read_mesh(mesh_file, "UNITSQUARE_6X4")
|
||||
|
||||
field = Problem(Heat, "unit square, 6x4 triangular mesh", 1)
|
||||
field = Problem(PlaneHeat, "unit square, 6x4 triangular mesh", 1)
|
||||
field.elements = create_elements(mesh, "UNITSQUARE")
|
||||
field.properties.formulation = "2D"
|
||||
update!(field, "thermal conductivity", 1.0)
|
||||
update!(field, "heat source", -6.0)
|
||||
|
||||
|
||||
@@ -19,7 +19,7 @@ Results are calculated using Code Aster for comparison.
|
||||
rings = Problem(Heat, "RINGS", 1)
|
||||
# rings.elements = create_elements(mesh; element_type=:Tet4)
|
||||
rings.elements = create_elements(mesh, "RING1", "RING2")
|
||||
update!(rings.elements, "temperature thermal conductivity", 1.0)
|
||||
update!(rings.elements, "thermal conductivity", 1.0)
|
||||
bc_inner = Problem(Dirichlet, "INNER SURFACE", 1, "temperature")
|
||||
bc_inner.elements = create_elements(mesh, "RING1_INNER")
|
||||
bc_outer = Problem(Dirichlet, "OUTER SURFACE", 1, "temperature")
|
||||
@@ -49,4 +49,3 @@ Results are calculated using Code Aster for comparison.
|
||||
end
|
||||
@test passed
|
||||
end
|
||||
|
||||
|
||||
+9
-11
@@ -2,7 +2,7 @@
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using JuliaFEM
|
||||
using JuliaFEM.Testing
|
||||
using Base.Test
|
||||
|
||||
@testset "two increments, linear solver" begin
|
||||
X = Dict{Int, Vector{Float64}}(
|
||||
@@ -12,11 +12,10 @@ using JuliaFEM.Testing
|
||||
4 => [0.0,1.0])
|
||||
element = Element(Quad4, [1, 2, 3, 4])
|
||||
update!(element, "geometry", X)
|
||||
update!(element, "temperature thermal conductivity", 6.0)
|
||||
update!(element, "temperature load", 0.0 => 12.0)
|
||||
update!(element, "temperature load", 1.0 => 24.0)
|
||||
problem = Problem(Heat, "one element heat problem", 1)
|
||||
problem.properties.formulation = "2D"
|
||||
update!(element, "thermal conductivity", 6.0)
|
||||
update!(element, "heat source", 0.0 => 12.0)
|
||||
update!(element, "heat source", 1.0 => 24.0)
|
||||
problem = Problem(PlaneHeat, "one element heat problem", 1)
|
||||
push!(problem, element)
|
||||
boundary_element = Element(Seg2, [1, 2])
|
||||
update!(boundary_element, "geometry", X)
|
||||
@@ -51,11 +50,10 @@ end
|
||||
4 => [0.0,1.0])
|
||||
element = Element(Quad4, [1, 2, 3, 4])
|
||||
update!(element, "geometry", X)
|
||||
update!(element, "temperature thermal conductivity", 6.0)
|
||||
update!(element, "temperature load", 0.0 => 12.0)
|
||||
update!(element, "temperature load", 1.0 => 24.0)
|
||||
problem = Problem(Heat, "one element heat problem", 1)
|
||||
problem.properties.formulation = "2D"
|
||||
update!(element, "thermal conductivity", 6.0)
|
||||
update!(element, "heat source", 0.0 => 12.0)
|
||||
update!(element, "heat source", 1.0 => 24.0)
|
||||
problem = Problem(PlaneHeat, "one element heat problem", 1)
|
||||
push!(problem, element)
|
||||
boundary_element = Element(Seg2, [1, 2])
|
||||
update!(boundary_element, "geometry", X)
|
||||
|
||||
@@ -69,10 +69,9 @@ end
|
||||
el4 = Element(Seg2, [7, 8])
|
||||
update!([el1, el2, el3, el4], "geometry", X)
|
||||
update!([el1, el2], "density", 6.0)
|
||||
update!([el1, el2], "temperature thermal conductivity", 36.0)
|
||||
update!([el1, el2], "thermal conductivity", 36.0)
|
||||
update!([el3, el4], "temperature 1", 0.0)
|
||||
p1 = Problem(Heat, "combined body", 1)
|
||||
p1.properties.formulation = "2D"
|
||||
p1 = Problem(PlaneHeat, "combined body", 1)
|
||||
p2 = Problem(Dirichlet, "fixed ends", 1, "temperature")
|
||||
push!(p1, el1, el2)
|
||||
push!(p2, el3, el4)
|
||||
@@ -100,13 +99,11 @@ end
|
||||
el6 = Element(Seg2, [5, 6])
|
||||
update!([el1, el2, el3, el4, el5, el6], "geometry", X)
|
||||
update!([el1, el2], "density", 6.0)
|
||||
update!([el1, el2], "temperature thermal conductivity", 36.0)
|
||||
update!([el1, el2], "thermal conductivity", 36.0)
|
||||
update!([el3, el4], "temperature 1", 0.0)
|
||||
update!(el5, "master elements", [el6])
|
||||
p1 = Problem(Heat, "body 1", 1)
|
||||
p2 = Problem(Heat, "body 2", 1)
|
||||
p1.properties.formulation = "2D"
|
||||
p2.properties.formulation = "2D"
|
||||
p1 = Problem(PlaneHeat, "body 1", 1)
|
||||
p2 = Problem(PlaneHeat, "body 2", 1)
|
||||
p3 = Problem(Dirichlet, "fixed ends", 1, "temperature")
|
||||
p4 = Problem(Mortar, "interface between bodies", 1, "temperature")
|
||||
p4.properties.dimension = 1
|
||||
|
||||
@@ -60,15 +60,13 @@ end
|
||||
meshfile = @__DIR__() * "/testdata/block_2d.med"
|
||||
mesh = aster_read_mesh(meshfile)
|
||||
|
||||
upper = Problem(Heat, "upper", 1)
|
||||
upper.properties.formulation = "2D"
|
||||
upper = Problem(PlaneHeat, "upper", 1)
|
||||
upper.elements = create_elements(mesh, "UPPER")
|
||||
update!(upper.elements, "temperature thermal conductivity", 1.0)
|
||||
update!(upper.elements, "thermal conductivity", 1.0)
|
||||
|
||||
lower = Problem(Heat, "lower", 1)
|
||||
lower.properties.formulation = "2D"
|
||||
lower = Problem(PlaneHeat, "lower", 1)
|
||||
lower.elements = create_elements(mesh, "LOWER")
|
||||
update!(lower.elements, "temperature thermal conductivity", 1.0)
|
||||
update!(lower.elements, "thermal conductivity", 1.0)
|
||||
|
||||
bc_upper = Problem(Dirichlet, "upper boundary", 1, "temperature")
|
||||
bc_upper.elements = create_elements(mesh, "UPPER_TOP")
|
||||
@@ -171,7 +169,7 @@ end
|
||||
|
||||
solver = Solver(Linear)
|
||||
push!(solver, upper, lower, bc_upper, bc_lower, interface, bc_corner)
|
||||
|
||||
|
||||
solver()
|
||||
slave_elements = get_slave_elements(interface)
|
||||
node_ids, la = get_nodal_vector(slave_elements, "lambda", 0.0)
|
||||
|
||||
@@ -21,11 +21,11 @@ This is conforming mesh so result should match to the conforming situation.
|
||||
|
||||
ring1 = Problem(Heat, "RING1", 1)
|
||||
ring1.elements = create_elements(mesh, "RING1")
|
||||
update!(ring1.elements, "temperature thermal conductivity", 1.0)
|
||||
update!(ring1.elements, "thermal conductivity", 1.0)
|
||||
|
||||
ring2 = Problem(Heat, "RING2", 1)
|
||||
ring2.elements = create_elements(mesh, "RING2")
|
||||
update!(ring2.elements, "temperature thermal conductivity", 1.0)
|
||||
update!(ring2.elements, "thermal conductivity", 1.0)
|
||||
|
||||
bc_inner = Problem(Dirichlet, "INNER SURFACE", 1, "temperature")
|
||||
bc_inner.elements = create_elements(mesh, "RING1_INNER")
|
||||
@@ -60,4 +60,3 @@ This is conforming mesh so result should match to the conforming situation.
|
||||
end
|
||||
@test passed
|
||||
end
|
||||
|
||||
|
||||
@@ -16,11 +16,11 @@ tet10_meshfile = "test_problems_mortar_3d/tet10.inp"
|
||||
|
||||
upper = Problem(Heat, "UPPER", 1)
|
||||
upper.elements = create_elements(mesh, "UPPER")
|
||||
update!(upper, "temperature thermal conductivity", 1.0)
|
||||
update!(upper, "thermal conductivity", 1.0)
|
||||
|
||||
lower = Problem(Heat, "LOWER", 1)
|
||||
lower.elements = create_elements(mesh, "LOWER")
|
||||
update!(lower, "temperature thermal conductivity", 1.0)
|
||||
update!(lower, "thermal conductivity", 1.0)
|
||||
|
||||
bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature")
|
||||
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
|
||||
@@ -41,7 +41,7 @@ tet10_meshfile = "test_problems_mortar_3d/tet10.inp"
|
||||
add_results_writer!(solver, Xdmf("sl_lin_temp_results"; overwrite=true))
|
||||
|
||||
solver()
|
||||
|
||||
|
||||
node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0)
|
||||
T = [t[1] for t in temperature]
|
||||
minT = minimum(T)
|
||||
@@ -68,11 +68,11 @@ end
|
||||
|
||||
upper = Problem(Heat, "UPPER", 1)
|
||||
upper.elements = create_elements(mesh, "UPPER")
|
||||
update!(upper, "temperature thermal conductivity", 1.0)
|
||||
update!(upper, "thermal conductivity", 1.0)
|
||||
|
||||
lower = Problem(Heat, "LOWER", 1)
|
||||
lower.elements = create_elements(mesh, "LOWER")
|
||||
update!(lower, "temperature thermal conductivity", 1.0)
|
||||
update!(lower, "thermal conductivity", 1.0)
|
||||
|
||||
bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature")
|
||||
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
|
||||
@@ -95,7 +95,7 @@ end
|
||||
add_results_writer!(solver, Xdmf("dl_lin_temp_results"; overwrite=true))
|
||||
|
||||
solver()
|
||||
|
||||
|
||||
node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0)
|
||||
T = [t[1] for t in temperature]
|
||||
minT = minimum(T)
|
||||
@@ -111,11 +111,11 @@ end
|
||||
|
||||
upper = Problem(Heat, "UPPER", 1)
|
||||
upper.elements = create_elements(mesh, "UPPER")
|
||||
update!(upper, "temperature thermal conductivity", 1.0)
|
||||
update!(upper, "thermal conductivity", 1.0)
|
||||
|
||||
lower = Problem(Heat, "LOWER", 1)
|
||||
lower.elements = create_elements(mesh, "LOWER")
|
||||
update!(lower, "temperature thermal conductivity", 1.0)
|
||||
update!(lower, "thermal conductivity", 1.0)
|
||||
|
||||
bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature")
|
||||
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
|
||||
@@ -130,7 +130,7 @@ end
|
||||
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
|
||||
update!(interface_slave_elements, "master elements", interface_master_elements)
|
||||
interface.elements = [interface_master_elements; interface_slave_elements]
|
||||
|
||||
|
||||
interface.properties.linear_surface_elements = false
|
||||
interface.properties.split_quadratic_slave_elements = false
|
||||
interface.properties.split_quadratic_master_elements = false
|
||||
@@ -140,10 +140,10 @@ end
|
||||
solver = LinearSolver(upper, lower, bc_upper, bc_lower, interface)
|
||||
add_results_writer!(solver, Xdmf("sl_quad_temp_results"; overwrite=true))
|
||||
solver()
|
||||
|
||||
|
||||
node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0)
|
||||
T = [t[1] for t in temperature]
|
||||
|
||||
|
||||
minT = minimum(T)
|
||||
maxT = maximum(T)
|
||||
stdT = std(T)
|
||||
@@ -158,11 +158,11 @@ end
|
||||
|
||||
upper = Problem(Heat, "UPPER", 1)
|
||||
upper.elements = create_elements(mesh, "UPPER")
|
||||
update!(upper, "temperature thermal conductivity", 1.0)
|
||||
update!(upper, "thermal conductivity", 1.0)
|
||||
|
||||
lower = Problem(Heat, "LOWER", 1)
|
||||
lower.elements = create_elements(mesh, "LOWER")
|
||||
update!(lower, "temperature thermal conductivity", 1.0)
|
||||
update!(lower, "thermal conductivity", 1.0)
|
||||
|
||||
bc_upper = Problem(Dirichlet, "UPPER_TOP", 1, "temperature")
|
||||
bc_upper.elements = create_surface_elements(mesh, "UPPER_TOP")
|
||||
@@ -177,7 +177,7 @@ end
|
||||
interface_master_elements = create_surface_elements(mesh, "UPPER_TO_LOWER")
|
||||
update!(interface_slave_elements, "master elements", interface_master_elements)
|
||||
interface.elements = [interface_master_elements; interface_slave_elements]
|
||||
|
||||
|
||||
interface.properties.linear_surface_elements = false
|
||||
interface.properties.split_quadratic_slave_elements = false
|
||||
interface.properties.split_quadratic_master_elements = false
|
||||
@@ -187,7 +187,7 @@ end
|
||||
solver = LinearSolver(upper, lower, bc_upper, bc_lower, interface)
|
||||
add_results_writer!(solver, Xdmf("dl_quad_temp_results"; overwrite=true))
|
||||
solver()
|
||||
|
||||
|
||||
node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0)
|
||||
#node_ids, temperature = get_nodal_vector(interface_slave_elements, "temperature", 0.0)
|
||||
#=
|
||||
@@ -458,4 +458,3 @@ end
|
||||
info("tet10 block: max(abs(u3)) = $maxabsu3, std(abs(u3)) = $stdabsu3")
|
||||
@test isapprox(stdabsu3, 0.0; atol=1.0e-6)
|
||||
end
|
||||
|
||||
|
||||
Reference in New Issue
Block a user