diff --git a/REQUIRE b/REQUIRE index 8e2ccd7..264ed47 100644 --- a/REQUIRE +++ b/REQUIRE @@ -11,3 +11,4 @@ AsterReader FEMBasis FEMQuad Reexport +HeatTransfer diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 0a628b2..e98358f 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -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") diff --git a/src/problems_heat.jl b/src/problems_heat.jl deleted file mode 100644 index 5198e4e..0000000 --- a/src/problems_heat.jl +++ /dev/null @@ -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 diff --git a/test/runtests.jl b/test/runtests.jl index 2e5c321..bf6b4d3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -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 -----") diff --git a/test/test_elasticity_tet10_mass_matrix.jl b/test/test_elasticity_tet10_mass_matrix.jl index 7f05a5c..5a1352c 100644 --- a/test/test_elasticity_tet10_mass_matrix.jl +++ b/test/test_elasticity_tet10_mass_matrix.jl @@ -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) diff --git a/test/test_heat.jl b/test/test_heat.jl index e3345bd..7987045 100644 --- a/test/test_heat.jl +++ b/test/test_heat.jl @@ -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 diff --git a/test/test_heat_2.jl b/test/test_heat_2.jl index 8f2cdc5..93b6c23 100644 --- a/test/test_heat_2.jl +++ b/test/test_heat_2.jl @@ -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 diff --git a/test/test_heat_3.jl b/test/test_heat_3.jl index 8c332c4..2b2fd28 100644 --- a/test/test_heat_3.jl +++ b/test/test_heat_3.jl @@ -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) diff --git a/test/test_heat_3d_two_rings.jl b/test/test_heat_3d_two_rings.jl index d6efd8d..009e57c 100644 --- a/test/test_heat_3d_two_rings.jl +++ b/test/test_heat_3d_two_rings.jl @@ -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 - diff --git a/test/test_heat_4.jl b/test/test_heat_4.jl index b0cf699..a36cb39 100644 --- a/test/test_heat_4.jl +++ b/test/test_heat_4.jl @@ -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) diff --git a/test/test_modal_analysis.jl b/test/test_modal_analysis.jl index 5e05957..a414ad7 100644 --- a/test/test_modal_analysis.jl +++ b/test/test_modal_analysis.jl @@ -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 diff --git a/test/test_mortar_2d_mesh_tie.jl b/test/test_mortar_2d_mesh_tie.jl index de17751..56d1333 100644 --- a/test/test_mortar_2d_mesh_tie.jl +++ b/test/test_mortar_2d_mesh_tie.jl @@ -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) diff --git a/test/test_mortar_3d_mesh_tie_two_rings.jl b/test/test_mortar_3d_mesh_tie_two_rings.jl index 2bc81d7..ab36211 100644 --- a/test/test_mortar_3d_mesh_tie_two_rings.jl +++ b/test/test_mortar_3d_mesh_tie_two_rings.jl @@ -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 - diff --git a/test/test_problems_mortar_3d.jl b/test/test_problems_mortar_3d.jl index 73caaa2..51c868f 100644 --- a/test/test_problems_mortar_3d.jl +++ b/test/test_problems_mortar_3d.jl @@ -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 -