After static temperature field solutions, the heat flow in to solid bodies through specific boundary regions has to be measured.
I started with a simple test case: a rectangle with temperature clamping (Dirichlet boundary conditions) at two sides.
The first experiment uses 1. Kelvin at the left side and 0. Kelvin at the right side as temperature boundary conditions. After solving the heat equation, the heat flow “through” the left and right side are measured. The magnitude of the heat flows match the analytically predicted values. Both heat flows have the same magnitude and inverse signs. This is also expected because the same amount of energy goes in and out of the system.
The second experiment uses 1. Kelvin at the left side and 0. Kelvin at the bottom side as temperature boundary conditions. The magnitude of the measured heat flows differ significantly. This difference is not expected. The heat flow is measured in the normal direction of the edge. Which I suppose is the right way. What am I doing wrong?
from dolfin import *
import matplotlib.pyplot as plt
WIDTH = 1.
HEIGHT = .5
WIDTH_DIV = 20
HEIGHT_DIV = 10
mesh = RectangleMesh(Point(0., 0.), Point(WIDTH, HEIGHT), WIDTH_DIV, HEIGHT_DIV, "right/left")
plot(mesh, title=f"Mesh (right/left)")
class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0) and on_boundary
class RightBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], WIDTH) and on_boundary
class BottomBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0) and on_boundary
left_bdry = LeftBoundary()
right_bdry = RightBoundary()
bottom_bdry = BottomBoundary()
boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
LEFT, RIGHT, BOTTOM = 1, 2, 3
left_bdry.mark(boundary_parts, LEFT)
right_bdry.mark(boundary_parts, RIGHT)
bottom_bdry.mark(boundary_parts, BOTTOM)
ds = Measure("ds", subdomain_data=boundary_parts)
V = FunctionSpace(mesh, "CG", 2)
u, v = TestFunction(V), TrialFunction(V)
temp = Function(V, name="Temperature")
a = dot(grad(u), grad(v)) * dx # lhs
L = Constant(0) * u * dx # rhs
solver_parameters = {'linear_solver': 'umfpack'}
n = FacetNormal(mesh)
print("heat flow from left side to right side:")
bc_left_to_right = [
DirichletBC(V, Constant(1.), left_bdry),
DirichletBC(V, Constant(0.), right_bdry)
]
solve(a == L, temp, bc_left_to_right, solver_parameters=solver_parameters)
plt.figure()
p = plot(temp, mode="contour", title="Temperature")
plt.colorbar(p)
heat_left0 = dot(n, nabla_grad(temp)) * ds(LEFT)
heat_right = dot(n, nabla_grad(temp)) * ds(RIGHT)
print(f"measured heat left/right {assemble(heat_left0): .2f}/{assemble(heat_right): .2f}")
print("heat flow from left side to bottom side:")
bc_left_to_bottom = [
DirichletBC(V, Constant(1.), left_bdry),
DirichletBC(V, Constant(0.), bottom_bdry)
]
solve(a == L, temp, bc_left_to_bottom, solver_parameters=solver_parameters)
plt.figure()
p = plot(temp, mode="contour", title="Temperature")
plt.colorbar(p)
heat_left1 = dot(n, nabla_grad(temp)) * ds(LEFT)
heat_bottom = dot(n, nabla_grad(temp)) * ds(BOTTOM)
print(f"measured heat left/bottom {assemble(heat_left1): .2f}/{assemble(heat_bottom): .2f}")
# plt.show()
Output:
heat flow from left side to right side:
measured heat left/right 0.50/-0.50
heat flow from left side to bottom side:
measured heat left/bottom 2.63/-3.63