Simulating the heat equation with insulating material

Question

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 left side. I figured that the left portion would be untouched and with constant temperature, given its low conductivity.

$$ \nabla \cdot \rho \nabla T = 0 ~ \text{in} ~\Omega=[0,1]\times[0,1] $$ $$ \rho = \begin{cases} 1 ~ \text{if} ~x > 0.5 \\ 10^{-8} ~ \text{else} \end{cases} $$

Boundary conditions on the right half $$ \nabla T \cdot n = 5 ~\text{in} ~ \Gamma_N = \{(x,y)~ \forall ~ x \gt 0.5 ~\text{and}~ y= 1\} $$ $$ T = 0 ~\text{in} ~ \Gamma_{D_1} = \{(x,y)~ \forall ~ x \gt 0.5 ~\text{and}~ y= 0\} $$

Boundary condition on the left half $$ T = 0 ~\text{in} ~ \Gamma_{D_2} = \{(x,y)~ \forall ~ x = 0 \} $$

I am solving this problem with first order Lagrange elements. I was expecting to see my solution in the right half and then having the left half mostly constant with a value of 0 and rapidly changing close to the interface for continuity. This would make sense given that the left half is mostly an insulating material. What I am seeing instead is a smooth transition from the interface to the boundary. Is there something wrong with my mathematical implementation? Attaching the FEniCS code if it can help

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-8", 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")
left = CompiledSubDomain("x[0] <= 0.001")
meshfunc_ds = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)

TOP, LEFT, BOTTOMRIGHT = 1, 2, 3
top.mark(meshfunc_ds, TOP)
left.mark(meshfunc_ds, LEFT)
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(0.0), meshfunc_ds, LEFT)

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

Answer 1

I was expecting to see my solution in the right half and then having the left half mostly constant with a value of 0 and rapidly changing close to the interface for continuity. This would make sense given that the left half is mostly an insulating material.

This assumption is wrong. Small conductivity simply means that you have small heat flux, but by no means this implies that the temperature gradient should be small. Think of the insulating layer of a furnace: typically you have a smooth (almost linear) gradient of temperature from the hot side to the cold one, and not constant temperature across the thickness.

What I am seeing instead is a smooth transition from the interface to the boundary.

This is what I would expect. Think it in this way: if there is a heat flux across the boundary between $D_1$ (conducting) and $D_2$ (very small conductivity), this heat flux for continuity has to be compensated by an equal heat flux across the boundary of $D_2$: this is leads to a smooth gradient across the whole $D_2$.

Suppose to compute a sequence of solutions with constant conductivity $\rho_1$ in $D_1$ but decreasing conductivity $\rho_2$ in $D_2$; once $\rho_2 \ll \rho_1$ the temperature field in $D_2$ should not change much any more, only the heat flux across the interface will become smaller and smaller.

At the limit when the conductivity is null, the problems in $D_1$ and $D_2$ become only weakly coupled. The solution in $D_1$ obeys an adiabatic condition (zero flux) at the interface. In $D_2$ instead you will have a solution with a Dirichlet boundary condition on the interface, with the temperature given by the solution of the problem in domain $D_1$. By dimensional analysis you see that this solution is independent from the actual conductivity value, hence the conclusion that for $\rho_2 \ll \rho_1$ the solution in $D_2$ does not change much with $\rho_2$. (Reasoning based only on physical intuition, but should not be to difficult to prove, or to verify by numerical experimentation.)

The OP reasoning on the contrary is that for vanishing $\rho_2$ the temperature gradient in $D_2$ localises at the interface, converging to a temperature discontinuity at zero conductivity. But this type of solution seems to me impossible for the heat equation.

Answer 2

Why even bother to simulate the left half? Why not just simulate the right half with the left BC being the constant value? Seems simpler and more accurate to what you want to actually model, because currently you're just going to have a very low slope on the left half due to the low conductivity, but it won't be constant.