From dfec88629ecd68c1ea299dfe5116aaf7a4ae6fee Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Wed, 2 May 2018 19:42:52 +0300 Subject: [PATCH 1/2] New boundary condition + solve increment * New keywords "heat transfer coefficient" and "external temperature" defines new boundary condition for surface (heat exchange). * Make code ready for nonlinear iteration by solving K*du = f-KT instead of KT = f. --- src/HeatTransfer.jl | 31 ++++++++++++++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/src/HeatTransfer.jl b/src/HeatTransfer.jl index 90bb1a7..5096361 100644 --- a/src/HeatTransfer.jl +++ b/src/HeatTransfer.jl @@ -14,6 +14,18 @@ - `thermal conductivity` - `heat source` - `heat flux` +- `external temperature` +- `heat transfer coefficient` + +# 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 """ module HeatTransfer @@ -50,6 +62,10 @@ function assemble_elements!(problem::Problem{P}, assembly::Assembly, fe += s * N'*f end end + if haskey(element, "temperature") + T = [interpolate(element["temperature"], time)...] + fe -= Ke*T + end gdofs = get_gdofs(problem, element) add!(assembly.K, gdofs, gdofs, Ke) add!(assembly.f, gdofs, fe) @@ -74,19 +90,32 @@ function assemble_boundary_elements!{B}(problem::Problem, assembly::Assembly, bi = BasisInfo(B) ndofs = length(bi) + Ke = zeros(ndofs, ndofs) fe = zeros(ndofs) for element in elements fill!(fe, 0.0) - for ip in get_integration_points(element) + fill!(Ke, 0.0) + for ip in get_integration_points(element, 2) J, detJ, N, dN = element_info!(bi, element, ip, time) s = ip.weight * detJ if haskey(element, "heat flux") g = element("heat flux", ip, time) fe += s * N'*g end + if haskey(element, "heat transfer coefficient") && haskey(element, "external temperature") + h = element("heat transfer coefficient", ip, time) + Tu = element("external temperature", ip, time) + Ke += s * h*N'*N + fe += s * N'*h*Tu + end + end + if haskey(element, "temperature") + T = [interpolate(element["temperature"], time)...] + fe -= Ke*T end gdofs = get_gdofs(problem, element) + add!(assembly.K, gdofs, gdofs, Ke) add!(assembly.f, gdofs, fe) end From 5945bc554d0424bd607c3dd1b4e71629b132a768 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Thu, 3 May 2018 10:36:55 +0300 Subject: [PATCH 2/2] Separate tests and add coverage back to 100 %. --- src/HeatTransfer.jl | 4 +- test/runtests.jl | 55 +++---------------- test/test_heat_exchange_bc.jl | 27 +++++++++ test/test_heat_flux_bc_2d.jl | 25 +++++++++ test/test_heat_flux_bc_3d.jl | 27 +++++++++ ...est_stiffness_matrix_and_heat_source_2d.jl | 36 ++++++++++++ 6 files changed, 126 insertions(+), 48 deletions(-) create mode 100644 test/test_heat_exchange_bc.jl create mode 100644 test/test_heat_flux_bc_2d.jl create mode 100644 test/test_heat_flux_bc_3d.jl create mode 100644 test/test_stiffness_matrix_and_heat_source_2d.jl diff --git a/src/HeatTransfer.jl b/src/HeatTransfer.jl index 5096361..7dc3c41 100644 --- a/src/HeatTransfer.jl +++ b/src/HeatTransfer.jl @@ -63,7 +63,7 @@ function assemble_elements!(problem::Problem{P}, assembly::Assembly, end end if haskey(element, "temperature") - T = [interpolate(element["temperature"], time)...] + T = [element("temperature", time)...] fe -= Ke*T end gdofs = get_gdofs(problem, element) @@ -111,7 +111,7 @@ function assemble_boundary_elements!{B}(problem::Problem, assembly::Assembly, end end if haskey(element, "temperature") - T = [interpolate(element["temperature"], time)...] + T = [element("temperature", time)...] fe -= Ke*T end gdofs = get_gdofs(problem, element) diff --git a/test/runtests.jl b/test/runtests.jl index 294405c..5cb05d4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,57 +7,20 @@ using HeatTransfer @testset "HeatTransfer.jl" begin -@testset "stiffness matrix and heat source" begin +@testset "stiffness matrix and heat source boundary condition" begin + include("test_stiffness_matrix_and_heat_source_2d.jl") +end - X = Dict(1 => [0.0,0.0], 2 => [1.0,0.0], 3 => [1.0,1.0], 4 => [0.0,1.0]) - el1 = Element(Quad4, [1, 2, 3, 4]) - update!(el1, "geometry", X) - update!(el1, "thermal conductivity", 6.0) - update!(el1, "heat source", 4.0) - assembly = Assembly() - problem = Problem(PlaneHeat, "test problem", 1) - time = 0.0 - assemble_elements!(problem, assembly, [el1], time) - K = full(assembly.K) - f = full(assembly.f) - K_expected = [ - 4.0 -1.0 -2.0 -1.0 - -1.0 4.0 -1.0 -2.0 - -2.0 -1.0 4.0 -1.0 - -1.0 -2.0 -1.0 4.0 - ] - f_expected = [1.0; 1.0; 1.0; 1.0] - @test isapprox(K, K_expected) - @test isapprox(f, f_expected) - @test get_unknown_field_name(problem) == "temperature" +@testset "heat flux boundary condition for 2d problem" begin + include("test_heat_flux_bc_2d.jl") end -@testset "heat flux, plane heat" begin - X = Dict(1 => [0.0,0.0], 2 => [1.0,0.0]) - el2 = Element(Seg2, [1, 2]) - update!(el2, "geometry", X) - update!(el2, "heat flux", 2.0) - assembly = Assembly() - problem = Problem(PlaneHeat, "test problem", 1) - time = 0.0 - assemble_elements!(problem, assembly, [el2], time) - f = full(assembly.f) - f_expected = [1.0; 1.0] - @test isapprox(f, f_expected) +@testset "heat flux boundary condition for 3d problem" begin + include("test_heat_flux_bc_3d.jl") end -@testset "heat flux, 3d heat" begin - X = Dict(1 => [0.0,0.0,0.0], 2 => [1.0,0.0,0.0], 3 => [0.0,1.0,0.0]) - el = Element(Tri3, [1, 2, 3]) - update!(el, "geometry", X) - update!(el, "heat flux", 6.0) - assembly = Assembly() - problem = Problem(Heat, "test problem", 1) - time = 0.0 - assemble_elements!(problem, assembly, [el], time) - f = full(assembly.f) - f_expected = [1.0; 1.0; 1.0] - @test isapprox(f, f_expected) +@testset "heat exchange boundary condition" begin + include("test_heat_exchange_bc.jl") end end diff --git a/test/test_heat_exchange_bc.jl b/test/test_heat_exchange_bc.jl new file mode 100644 index 0000000..3378b66 --- /dev/null +++ b/test/test_heat_exchange_bc.jl @@ -0,0 +1,27 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/HeatTransfer.jl/blob/master/LICENSE + +using HeatTransfer +using FEMBase +using FEMBase.Test + +X = Dict( + 1 => [0.0,0.0], + 2 => [1.0,0.0]) + +T = Dict( + 1 => 0.0, + 2 => 0.0) + +element = Element(Seg2, [1, 2]) +update!(element, "geometry", X) +update!(element, "temperature", 0.0 => T) +update!(element, "heat transfer coefficient", 6.0) +update!(element, "external temperature", 1.0) +problem = Problem(PlaneHeat, "test problem", 1) +add_elements!(problem, [element]) +assemble!(problem, 0.0) +K = full(problem.assembly.K) +f = full(problem.assembly.f) +@test isapprox(K, [2.0 1.0; 1.0 2.0]) +@test isapprox(f, [3.0; 3.0]) diff --git a/test/test_heat_flux_bc_2d.jl b/test/test_heat_flux_bc_2d.jl new file mode 100644 index 0000000..a561d1d --- /dev/null +++ b/test/test_heat_flux_bc_2d.jl @@ -0,0 +1,25 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/HeatTransfer.jl/blob/master/LICENSE + +using HeatTransfer +using FEMBase +using FEMBase.Test + +X = Dict( + 1 => [0.0,0.0], + 2 => [1.0,0.0]) + +T = Dict( + 1 => 0.0, + 2 => 0.0) + +element = Element(Seg2, [1, 2]) +update!(element, "geometry", X) +update!(element, "temperature", 0.0 => T) +update!(element, "heat flux", 2.0) +problem = Problem(PlaneHeat, "test problem", 1) +add_elements!(problem, [element]) +assemble!(problem, 0.0) +f = full(problem.assembly.f) +f_expected = [1.0; 1.0] +@test isapprox(f, f_expected) diff --git a/test/test_heat_flux_bc_3d.jl b/test/test_heat_flux_bc_3d.jl new file mode 100644 index 0000000..03d3d4a --- /dev/null +++ b/test/test_heat_flux_bc_3d.jl @@ -0,0 +1,27 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/HeatTransfer.jl/blob/master/LICENSE + +using HeatTransfer +using FEMBase +using FEMBase.Test + +X = Dict( + 1 => [0.0,0.0,0.0], + 2 => [1.0,0.0,0.0], + 3 => [0.0,1.0,0.0]) + +T = Dict( + 1 => 0.0, + 2 => 0.0, + 3 => 0.0) + +element = Element(Tri3, [1, 2, 3]) +update!(element, "geometry", X) +update!(element, "temperature", T) +update!(element, "heat flux", 6.0) +problem = Problem(Heat, "test problem", 1) +add_elements!(problem, [element]) +assemble!(problem, 0.0) +f = full(problem.assembly.f) +f_expected = [1.0; 1.0; 1.0] +@test isapprox(f, f_expected) diff --git a/test/test_stiffness_matrix_and_heat_source_2d.jl b/test/test_stiffness_matrix_and_heat_source_2d.jl new file mode 100644 index 0000000..b8ddf42 --- /dev/null +++ b/test/test_stiffness_matrix_and_heat_source_2d.jl @@ -0,0 +1,36 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/HeatTransfer.jl/blob/master/LICENSE + +X = Dict( + 1 => [0.0,0.0], + 2 => [1.0,0.0], + 3 => [1.0,1.0], + 4 => [0.0,1.0]) + +T = Dict( + 1 => 0.0, + 2 => 0.0, + 3 => 0.0, + 4 => 0.0) + +element = Element(Quad4, [1, 2, 3, 4]) +update!(element, "geometry", X) +update!(element, "temperature", 0.0 => T) +update!(element, "thermal conductivity", 6.0) +update!(element, "heat source", 4.0) + +problem = Problem(PlaneHeat, "test problem", 1) +add_elements!(problem, [element]) +assemble!(problem, 0.0) +K = full(problem.assembly.K) +f = full(problem.assembly.f) +K_expected = [ + 4.0 -1.0 -2.0 -1.0 + -1.0 4.0 -1.0 -2.0 + -2.0 -1.0 4.0 -1.0 + -1.0 -2.0 -1.0 4.0 + ] +f_expected = [1.0; 1.0; 1.0; 1.0] +@test isapprox(K, K_expected) +@test isapprox(f, f_expected) +@test get_unknown_field_name(problem) == "temperature"