如何用FiPy正确模拟矩形环热扩散？边界条件设置求助

问题描述 #

使用FiPy模拟铜制矩形环热扩散，设置外边界（绿点）温度298.15K，内边界（蓝点）温度303.15K，但模拟中出现温度降至0K的异常。

问题根源 #

原代码中手动计算内边界面ID的逻辑错误，导致错误标记大量非边界面，引发边界条件冲突，最终造成求解器计算异常。FiPy的面ID顺序与手动计算逻辑不匹配，直接操作面ID极易出错。

修正方案 #

改用坐标筛选的方式生成内边界面掩码，避免手动计算面ID的错误，确保边界条件仅应用在正确的面上。

修正后完整代码 #

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

def get_inner_boundary_faces(mesh):
    # 通过面中心坐标筛选内边界面：x=0.25/0.75且y在0.25-0.75之间，或y=0.25/0.75且x在0.25-0.75之间
    x_face = mesh.faceCenters[0]
    y_face = mesh.faceCenters[1]
    left_face = (np.isclose(x_face, 0.25)) & (y_face >= 0.25) & (y_face <= 0.75)
    right_face = (np.isclose(x_face, 0.75)) & (y_face >= 0.25) & (y_face <= 0.75)
    bottom_face = (np.isclose(y_face, 0.25)) & (x_face >= 0.25) & (x_face <= 0.75)
    top_face = (np.isclose(y_face, 0.75)) & (x_face >= 0.25) & (x_face <= 0.75)
    return left_face | right_face | bottom_face | top_face

# 初始化网格与参数
mesh = fipy.Grid2D(Lx=1, Ly=1, nx=20, ny=20)  # 1m×1m网格
k_over_c_rho = 3.98E2 / (3.85E2 * 8.96E3)  # 铜的热扩散率（MKS单位）
dt = 0.1 * (mesh.dx**2 + mesh.dy**2) / (4 * k_over_c_rho)  # 稳定时间步长
T0 = 273.15  # 0℃对应的开尔文温度

# 初始化温度场
T = fipy.CellVariable(mesh, name='T', value=T0 + 25)  # 初始温度298.15K

# 设置外边界条件
mask_e = mesh.exteriorFaces
T.constrain(T0 + 25., mask_e)

# 设置内边界条件
mask_i = get_inner_boundary_faces(mesh)
T.constrain(T0 + 30, mask_i)

# 定义热扩散方程
eq = fipy.TransientTerm() == fipy.DiffusionTerm(coeff=k_over_c_rho)

# 可视化设置
viewer = fipy.MatplotlibViewer(vars=[T], datamin=290, datamax=305)  # 缩小显示范围，聚焦有效温度区间
plt.ioff()
viewer._plot()
plt.plot(*mesh.faceCenters[:, mask_e], '.g', label='外边界')
plt.plot(*mesh.faceCenters[:, mask_i], '.b', label='内边界')
plt.legend()

# 更新模拟函数
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()

关键修改说明 #

1. 内边界面生成逻辑：通过mesh.faceCenters直接获取面中心坐标，利用坐标匹配筛选内边界面，完全避免手动计算面ID的错误。
2. 可视化优化：调整datamin和datamax为合理的温度区间（290-305K），避免显示无关的0K异常值。
3. 边界条件验证：确保外边界和内边界的掩码没有重叠，不会出现同一面被多次约束的冲突。

内容的提问来源于stack exchange，提问作者relent95