如何使用 FiPy 模拟矩形环上的热扩散？答案

【问题标题】：How to simulate a heat diffusion on a rectangular ring with FiPy?如何使用 FiPy 模拟矩形环上的热扩散？
【发布时间】：2023-01-18 13:41:21
【问题描述】：

我刚开始使用 FiPy 求解 PDE 并在矩形环状铜体上进行热扩散实验。

这是某些时候的模拟结果图。

我使用 Grid2D() 表示网格，使用 CellVariable.constrain() 指定边界条件。绿点是外表面的中心吨= 273.15 + 25 (钾), 蓝点是内面的中心吨= 273.15 + 30 (钾).

显然，我做错了什么，因为温度下降到 0钾.我应该如何正确指定边界条件？

这些是代码。

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

【问题讨论】：

你是如何选择你的 dt 的？是否符合某些 CFL 条件？是不是缺少某处缺少平方根？也许不是......我也很不明白“绿点是外表面的中心，其中T = 273.15 + 25（K），蓝点是内表面的中心，其中T = 273.15 + 30（K）。”您似乎在您的域内强加了条件，因为您的正方形内印有蓝点，这正常吗？ get_mask_of_rect 函数真的在做你想做的事吗？

标签： python simulation fipy

【解决方案1】：

.constrain() 不适用于内部面孔（请参阅Applying internal “boundary” conditions 末尾的警告）。

但是，您可以使用源实现 internal fixed value 条件。作为第一次切割，尝试

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

这将mask_i 标识的面部两侧的单元格限制为T0+30。

【讨论】：