Simulating the heat equation with insulating material

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

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
$$$$
`
• This code wont run. You have "left" which is not defined. Instead you have a bottom_left. Check your bc. Some bc parts are not specified. Do you want to apply zero neumann bc on those parts ? – cfdlab Aug 15 '19 at 13:41
• Also use "near(x[0],0.0)" to test for "x==0", etc. – cfdlab Aug 15 '19 at 13:42
• sorry, fixed now. – balborian Aug 15 '19 at 17:14

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.

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.

• I need the empty half because this is a step towards modeling the material embedded in a global domain using immersed methods. – balborian Aug 15 '19 at 17:17
• well if your conductivity is non zero, you're gonna get heat transfer. If its zero then you wont. Alternatively, I guess you could make the entire domain have normal conductivity, but set the left side boundary to the constant values, and set the right most cells on the left side to the same constant values. – EMP Aug 15 '19 at 22:06
• "very low slope on the left half due to the low conductivity" seems a wrong assumption to me, as much as this translates into the expectation that you decrease the the slope by decreasing the conductivity. See my answer for a full explanation. By the way the opposite seems true: you have a vanishing slope if you have very high conductivity in the left domain. – Stefano M Aug 15 '19 at 23:05
• If you think of the poisson equation though, which this basically is aside from the discontinuity in conductivity, then a lower coefficient gives a lower gradient, iirc. – EMP Aug 16 '19 at 20:18