Merge 5945bc5 into bcf1446

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 = [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 = [element("temperature", time)...]
+            fe -= Ke*T
         end
         gdofs = get_gdofs(problem, element)
+        add!(assembly.K, gdofs, gdofs, Ke)
         add!(assembly.f, gdofs, fe)
     end
 

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

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])

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)

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)

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"