CODEWITHSUNDEEP
pythonsimulationfipy
relent95
relent95

Reputation: 4732

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

Upvotes: 0

Views: 187

Answers (2)

relent95
relent95

Reputation: 4732

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)

Upvotes: 2

jeguyer
jeguyer

Reputation: 2484

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

Upvotes: 2

Related Questions