Search code examples
pythonsimulationfipy

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

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. enter image description here

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()
Solution

  • .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.