### How to simulate heat transfer between two bodies separated by a thin thermally resistive layer?

**Creation and preparation of the geometry.**

The toy problem consists of two bodies or more specifically two blocks of different size that are connected to each other, i.e. they share a common interface. The two blocks are created with gmesh (http://gmsh.info/)*.*The two blocks are marked as physical volumes and the interface, i.e. the common face of the small block inside the face of the big block is marked as the physical surface "Interface". The mesh is saved to*HeatFlow_ThermalResistance.msh*and converted with dolfin-convert.**Mathematics**.

The task is to solve the transient heat equation in its weak form

$F^{\left(i+1\right)}=\int_{\Omega}^{ }\left(Tv+\Delta t\frac{k}{c_p\rho}\nabla T\cdot\nabla v-\left(T^{\left(i\right)}+\Delta t\frac{q_V}{c_p\rho}\right)v\right)dx-\int_{\partial\Omega}^{ }k\frac{\partial T}{\partial n}vds$`F`^{(i+1)}=∫_{Ω}(`T``v`+Δ`t``k``c`_{p}`ρ`∇`T`·∇`v`−(`T`^{(i)}+Δ`t``q`_{V}`c`_{p}`ρ`)`v`)`d``x`−∫_{∂Ω}`k`∂`T`∂`n``v``d``s`

The integration over the cells is solved for each block separately to account for different physical properties of the blocks. No boundary conditions are set, corresponding to thermally insulated boundaries. The initial condition for the transient part is $T^{\left(0\right)}\left(t=0\right)=T_0$`T`^{(0)}(`t`=0)=`T`_{0}.**Python code.**The following Python code demonstrates the problem. Here, the small block contains a heat source ( $q_v^{\left(1\right)}>0$`q`_{v}^{(1)}>0 ), hence it heats up.

`from dolfin import * vtkfile = File('HeatFlow_ThermalResistance.pvd') # get the mesh for the cylinder mesh = Mesh('HeatFlow_ThermalResistance.xml') # get markers domain_markers = MeshFunction('size_t', mesh, 'HeatFlow_ThermalResistance_physical_region.xml') boundary_markers = MeshFunction('size_t', mesh, 'HeatFlow_ThermalResistance_facet_region.xml') # setup markers (dx(1)=small block, dx(2)=big block, ds(1)=interface) dx = Measure('dx', domain=mesh, subdomain_data=domain_markers) ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers) # setup material properties class small_block: q_v = 1e3 # volumetric heat source rho = 1.0 # density c_p = 1.0 # specific heat capacity k = 300.0 # thermal conductivity class big_block: q_v = 0 rho = 1.0 c_p = 1.0 k = 1000.0 # setup function space V = FunctionSpace(mesh, 'CG', 1) T, v = TrialFunction(V), TestFunction(V) # initial condition T_n = project(Constant(300), V) # time stepping t_end = 0.1 # final time num_steps = 200 # number of time steps dt = t_end / num_steps # time step size # setup form # small box F = T * v * dx(1) + \ dt * small_block.k / (small_block.c_p * small_block.rho) * \ dot(grad(T), grad(v)) * dx(1) - \ (T_n + dt / (small_block.c_p * small_block.rho) * small_block.q_v) * v * dx(1) # big box F += T * v * dx(2) + \ dt * big_block.k / (big_block.c_p * big_block.rho) * \ dot(grad(T), grad(v)) * dx(2) - \ (T_n + dt / (big_block.c_p * big_block.rho) * big_block.q_v) * v * dx(2) a, L = lhs(F), rhs(F) T = Function(V) t = 0 for nn in range(num_steps): t += dt # solve solve(a == L, T, []) # save to vtk file vtkfile << (T, t) # update previous solution T_n.assign(T)`

**Problems/Questions.**When the small block becomes heated up, no heat is transferred to the larger block as shown in the following at end of the simulation at $t=100ms$`t`=100`m``s`:

As one can see there is no heat flux from the small block to the big block. I think this makes sense, since the interface is thermally isolated. One possibility to allow a heat flow between the two blocks is to mark the two blocks as a 'fragment' in gmsh. This is done by appending the line

to the file`BooleanFragments{ Volume{2}; Delete; }{ Volume{1}; Delete; }`

*HeatFlow_ThermalResistance.geo*presented in 1. Note that one usually must change the plane of the physical surface since the number of the surface changes due to the BooleanFragments operation. Then heat flows from the smaller block to the big block.

This is where my problems begin. The blocks should be separated by an infinitesimally thin thermally resistive layer characterized by an thermal resistance $R_{th}$`R`_{th}having the unit $\frac{Km^2}{W}$`K``m`^{2}`W`.

My first idea is to remove the operation BooleanFragments from the geo file. As a result the two blocks are separated resulting in the heat distribution shown above. Formally, I then have to put boundary conditions on the interface of the two blocks $\Gamma_1$Γ_{1}:

$k^{\left(1\right)}\vec{n_1}\cdot\nabla T|_{\Gamma_1\left(-\right)}=-\frac{T_2-T_1}{R_{th}}$`k`^{(1)}→`n`_{1}·∇`T`|_{Γ1(−)}=−`T`_{2}−`T`_{1}`R`_{th}on the "left" side (-) of the interface $\Gamma_1$Γ_{1}

$k^{\left(2\right)}\vec{n_2}\cdot\nabla T|_{\Gamma_1\left(+\right)}=-\frac{T_1-T_2}{R_{th}}$`k`^{(2)}→`n`_{2}·∇`T`|_{Γ1(+)}=−`T`_{1}−`T`_{2}`R`_{th}on the "right" side (+) of the interface $\Gamma_1$Γ_{1}

Note that $T_2$`T`_{2}is the temperature on the "right" side of the interface and $T_1$`T`_{1}is the temperature at the "left" side of the interface. Obviously there is a discontinuity, but does that matter here? How can I put this kind of boundary condition into a form? I tried something like

This form is invalid and results in an`# thermal resistance R_th = 20 F = T * v * dx(1) + \ dt * small_block.k / (small_block.c_p * small_block.rho) * \ dot(grad(T), grad(v)) * dx(1) - \ (T_n + dt / (small_block.c_p * small_block.rho) * small_block.q_v) * v * dx(1) -\ (T('+') - T('-')) / R_th * v('-') * ds(1) # big box F += T * v * dx(2) + \ dt * big_block.k / (big_block.c_p * big_block.rho) * \ dot(grad(T), grad(v)) * dx(2) - \ (T_n + dt / (big_block.c_p * big_block.rho) * big_block.q_v) * v * dx(2) -\ (T('-') - T('+')) / R_th * v('+') * ds(1) `

*IndexError: index out of range*error that is caused by the T('-') and T('+') terms. Any ideas or strategies how to solve this toy problem?

`dS = Measure('dS', domain=mesh, subdomain_data=boundary_markers)`

to mark the interior surfaces and changed the ds(1) in the form to dS(1). There is still no heat flux between the bodies.

### 1 Answer

`on_boundary`

part of the `SubDomain`

class looks like for your mesh, naively using```
V = FunctionSpace(mesh, 'CG', 1)
def boundary(x, on_boundary):
return on_boundary
dbc = DirichletBC(V, Constant(1.0), boundary)
u = Function(V)
dbc.apply(u.vector())
```

and visualising the function u. The result looks like this:

Note how the boundary that's supposed to be interior also is treated as exterior boundary.

Unfortunately don't know much about gmsh so I redid your mesh in Salome and did the same `on_boundary`

check. The result looks like this:

So the interior boundary is not treated as exterior boundary. And running your code **even without jump terms **produced this result for the temperature:

So there is heat flow now. I attached the mesh I used so maybe you can look into what's wrong. In the attached mesh, the interface is also marked as 1 and the volumes as 1 and 2.

File attached: salomemesh.zip (21.68 KB)

`dS`

instead of`ds`

. The capital`dS`

is used for interior surfaces whereas`ds`

is used for the exterior boundary.