How to simulate a heat diffusion on a rectangular ring with FiPy?

Question

I am new to solving a PDE and experimenting with a heat diffusion on a copper body of a rectangular ring shape using FiPy.

And this is a plot of simulation result at some times.

I am using the Grid2D() for a mesh and the CellVariable.constrain() to specify boundary conditions. The green dots are centers of exterior faces where T = 273.15 + 25 (K), and blue dots are centers of interior faces where T = 273.15 + 30 (K).

Obviously, I am doing something wrong, because the temperature goes down to 0K. How should I specify boundary conditions correctly?

These are the code.

import numpy as np
import matplotlib.pyplot as plt
import fipy

def get_mask_of_rect(mesh, x, y, w, h):
    def left_id(i, j): return mesh.numberOfHorizontalFaces + i*mesh.numberOfVerticalColumns + j
    def right_id(i, j): return mesh.numberOfHorizontalFaces + i*mesh.numberOfVerticalColumns + j + 1
    def bottom_id(i, j): return i*mesh.nx + j
    def top_id(i, j): return (i+1)*mesh.nx + j
    j0, i0 = np.floor(np.array([x, y]) / [mesh.dx, mesh.dy]).astype(int)
    n, m = np.round(np.array([w, h]) / [mesh.dx, mesh.dy]).astype(int)
    mask = np.zeros_like(mesh.exteriorFaces, dtype=bool)
    for i in range(i0, i0 + n):
        mask[left_id(i, j0)] = mask[right_id(i, j0 + m-1)] = True
    for j in range(j0, j0 + m):
        mask[bottom_id(i0, j)] = mask[top_id(i0 + n-1, j)] = True
    return mask

mesh = fipy.Grid2D(Lx = 1, Ly = 1, nx = 20, ny = 20) # Grid of size 1m x 1m
k_over_c_rho = 3.98E2 / (3.85E2 * 8.96E3) # The thermal conductivity, specific heat capacity, and density of Copper in MKS
dt = 0.1 * (mesh.dx**2 + mesh.dy**2) / (4*k_over_c_rho)
T0 = 273.15 # 0 degree Celsius in Kelvin

T = fipy.CellVariable(mesh, name='T', value=T0+25)
mask_e = mesh.exteriorFaces
T.constrain(T0+25., mask_e)

mask_i = get_mask_of_rect(mesh, 0.25, 0.25, 0.5, 0.5)
T.constrain(T0+30, mask_i)

eq = fipy.TransientTerm() == fipy.DiffusionTerm(coeff=k_over_c_rho)
viewer = fipy.MatplotlibViewer(vars=[T], datamin=0, datamax=400)
plt.ioff()
viewer._plot()
plt.plot(*mesh.faceCenters[:, mask_e], '.g')
plt.plot(*mesh.faceCenters[:, mask_i], '.b')
def update():
    for _ in range(10):
        eq.solve(var=T, dt=dt)
    viewer._plot()
    plt.draw()
timer = plt.gcf().canvas.new_timer(interval=50)
timer.add_callback(update)
timer.start()

plt.show()

How did you choose your dt ? Is it according to some CFL condition ? Isn't missing missing a square root somewhere ? Maybe not... Also I quite don't understand "The green dots are centers of exterior faces where T = 273.15 + 25 (K), and blue dots are centers of interior faces where T = 273.15 + 30 (K)." You seemed to impose condition inside you domain, as blue dots are printed inside your square, is this normal ? Is the function `get_mask_of_rect` is really doing what you want ? — Thibault Cimic, Jan 17 '23 at 10:57
@Thibault Cimic, I chose dt like the official sample - such as ```timeStepDuration = 10 * 0.9 * dx**2 / (2 * D)```. And for the blue dots, it was due to my ignorance on specifying internal boundary conditions. Now it's clear from the answer of jeguyer. — relent95, Jan 18 '23 at 01:46

score 2 · Accepted Answer · answered Jan 17 '23 at 23:22

.constrain() does not work for internal faces (see the warning at the end of Applying internal “boundary” conditions).

You can achieve an internal fixed value condition using sources, however. As a first cut, try

mask_i = get_mask_of_rect(mesh, 0.25, 0.25, 0.5, 0.5)
mask_i_cell = fipy.CellVariable(mesh, value=False)
mask_i_cell[mesh.faceCellIDs[..., mask_i]] = True
largeValue = 1e6

eq = (fipy.TransientTerm() == fipy.DiffusionTerm(coeff=k_over_c_rho)
      - fipy.ImplicitSourceTerm(mask_i_cell * largeValue)
      + mask_i_cell * largeValue * (T0 + 30))

This constrains the cells on either side of the faces identified by mask_i to be at T0+30.

score 2 · Answer 2 · answered Jan 18 '23 at 04:04

I am posting my own answer for future readers. For boundary conditions for internal faces, you should use the implicit and explicit source terms on the equation, as in the jeguyer's answer.

By using source terms, you don't need to calculate a mask for faces, like this.(The get_mask_of_rect() in my question isn't required.)

T = fipy.CellVariable(mesh, name = 'T', value = T0 + 25)
mask_e = mesh.exteriorFaces
T.constrain(T0 + 25., mask_e)

mask_i_cell = (
    (0.25 < mesh.x) & (mesh.x < 0.25 + 0.5) &
    (0.25 < mesh.y) & (mesh.y < 0.25 + 0.5)
)
large_value = 1E6
eq = fipy.TransientTerm() == (
    fipy.DiffusionTerm(coeff = k_over_c_rho) -
    fipy.ImplicitSourceTerm(mask_i_cell * large_value) +
    mask_i_cell * (large_value * (T0 + 30) # explicit source
))

viewer = fipy.MatplotlibViewer(vars = [T], datamin = T0, datamax = T0+50)

How to simulate a heat diffusion on a rectangular ring with FiPy?

2 Answers2