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

4 months ago by
I have a problem with a toy problem. I am going to describe my complete workflow here in the hope that anyone finds a major flaw or has useful ideas for improvement. All files needed for a working example can be downloaded here (File attached: (26.9 KB)).

  1. 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 ( 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.
  2. 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)=Ω(Tv+Δtkcpρ T·v(T(i)+ΔtqVcpρ )v)dx∂ΩkTn vds   
    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)=T0.
  3. Python code.
    The following Python code demonstrates the problem. Here, the small block contains a heat source ( $q_v^{\left(1\right)}>0$qv(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
  4. 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=100ms:

    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

    BooleanFragments{ Volume{2}; Delete; }{ Volume{1}; Delete; }
    to the file 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}$Rth having the unit $\frac{Km^2}{W}$Km2W  .
    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)n1·T|Γ1()=T2T1Rth    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)n2·T|Γ1(+)=T1T2Rth    on the "right" side (+) of the interface $\Gamma_1$Γ1
    Note that  $T_2$T2 is the temperature on the "right" side of the interface and  $T_1$T1 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
    # 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)
    This form is invalid and results in an 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?
Community: FEniCS Project
To me, your interface jump condition looks reasonable. I am guessing that you should use dS instead of ds. The capital dS is used for interior surfaces whereas ds is used for the exterior boundary.
written 4 months ago by klunkean  
Thank you very much for your answer. I added the line
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.

written 3 months ago by Christian Vorholt  

1 Answer

3 months ago by
So I looked a bit into your mesh and I looked at what the 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)

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: (21.68 KB)

Please login to add an answer/comment or follow this question.

Similar posts:
Search »