Heat equation with insulating material

[miguel]

My plan is to solve the heat equation in the right half portion of the domain, while having the left half completely isolated with constant temperature. To do so, I model the left half with a very low conductivity. I then apply a heat flux to the top side of the right half and a Dirichlet BC to the bottom side. For the left portion, I simply apply a Dirichlet BC to the bottom side. I figured that the left portion would be untouched and with constant temperature, given its low conductivity, but it is not the case when running the code below:

from dolfin import *

mesh = UnitSquareMesh(100, 100)

V = FunctionSpace(mesh, 'CG', 1)
t, w = TrialFunction(V), TestFunction(V)

Rho = FunctionSpace(mesh, 'DG', 0)
rho = Function(Rho)
rho.interpolate(Expression("(x[0] > 0.5) + 1e-12", domain=mesh, degree=1))
File("test_rho.pvd") << rho

a = inner(rho*grad(t), grad(w))*dx

top = CompiledSubDomain("x[1] > 1 - 0.01 && x[0] >= 0.5")
bottom_right = CompiledSubDomain("x[0] >= 0.5 && x[1] < 0.01")
bottom_left = CompiledSubDomain("x[0] <= 0.5 && x[1] < 0.01")
meshfunc_ds = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)

TOP, BOTTOMLEFT, BOTTOMRIGHT = 1, 2, 3
top.mark(meshfunc_ds, TOP)
bottom_left.mark(meshfunc_ds, BOTTOMLEFT)
bottom_right.mark(meshfunc_ds, BOTTOMRIGHT)

File("test_measures.pvd") << meshfunc_ds

ds = Measure("ds")(subdomain_data=meshfunc_ds)
L = Constant(5.0)*w*ds(TOP)

bc1 = DirichletBC(V, Constant(0.0), meshfunc_ds, BOTTOMRIGHT)
bc2 = DirichletBC(V, Constant(-5.0), meshfunc_ds, BOTTOMLEFT)

t_sol = Function(V)
solve(a==L, t_sol, bcs=[bc1, bc2])
File("test_temperature.pvd") << t_sol

I must not be understanding something about the heat equation. How could I obtain a constant temperature in the left half?

[SantiagoOrtiz]

Miguel,

I understand that within this context, the source therm (f = Constant(5.0)) accounts for an external energetic disturbance (positive for a heat source and negative for heat removal). Indeed, the FEniCS Tutorial (Chapter 4) provides multiple guidelines to solve a similar problem as the one you come along with, defining a parameter related to heat conductivity over different subdomains.

On the other hand, the Dirichlet BCs account for the wall temperatures. I couldn’t code two distinct Dirichlet BC (DBC) at the bottom (left DBC: x[0]<=5 and right DBC: x[0]>=5), Instead I set a wall temperature of 0.0 (bc = DirichletBC(V, Constant(0.0), BoundBottom)) over the whole boundary.

Another FEniCS user may give an idea of the correct way to split the boundary and stablish two different DBCs…

from dolfin import UnitSquareMesh, FunctionSpace, TrialFunction, \
                   TestFunction, Constant, dot, grad, dx, Function, solve, \
                   SubDomain, DirichletBC, near, XDMFFile, Expression


mesh = UnitSquareMesh(100, 100)
V = FunctionSpace(mesh, 'P', 1)

tol = 1E-14
#"""
k_0 = 2.5
#k_0 = 5.0
k_1 = 0.5
kappa = Expression('x[0] <= 0.5 - tol ? k_0 : k_1', degree=0,
                   tol=tol, k_0=k_0, k_1=k_1)
#"""
#kappa = k_0 = k_1 = 2.5

u = TrialFunction(V)
v = TestFunction(V)
f = Constant(5.0)
a = kappa*dot(grad(u), grad(v))*dx
L = f*v*dx


class BoundaryBottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0, tol)


BoundBottom = BoundaryBottom()
bc = DirichletBC(V, Constant(0.0), BoundBottom)

xdmffile_U = XDMFFile('./Materials.xdmf')

# Compute solution
u = Function(V)
solve(a == L, u, bc)
xdmffile_U.write(u)

kappa = k_0 = k_1 = 2.5

k_0 = 2.5 and k_1 = 0.5

k_0 = 5.0 and k_1 = 0.5

The results evidence that by switching each Kappa parameter (k_1 & k_0), the thermal map change. If you want the left side to approximate to a constant temperature, you should keep increasing K_0.

I hope you find it helpful.

P.D.: In order to format your python code

At the beginning of your code write: ```python
And at the end of it, write again three (3) `

All the best, Santiago