如何在Gmsh中实现从圆心向边缘的扩散模拟[Python]

你当前代码出现反向扩散，完全是初始条件和边界条件设置错误导致的：

• 初始浓度场全域设为0，没有在圆心位置设置营养物质的初始高浓度，不符合“营养初始集中在圆心”的需求
• 边界判断阈值完全不合理：圆形域半径仅为1，你判断边界位置的阈值设为50，所有外边界到圆心的距离都是1，全部满足dr<50的判定条件，等于把整个圆形边缘的浓度强制固定为1，自然会出现物质从边缘往圆心扩散的反向效果。

按照营养从圆心向外扩散的需求，按以下逻辑调整即可：

1. 正确设置圆心初始高浓度
   不再把全域初始值统一设为0，先计算每个网格单元到圆心(0,0)的距离，把距离圆心小于指定小半径（和网格尺寸匹配即可，比如0.05）的核心区域浓度初始设为1，其余区域初始为0，对应营养物质初始集中在圆心供给区的设定。
2. 修正边界条件
   删掉原来无效的距离阈值判断，直接根据模拟场景设置外边界：如果要模拟营养扩散到肿瘤边缘后被代谢/带走，就给外边界设置固定浓度0；如果不考虑边缘营养流失，不用额外加约束，FiPy默认是零通量边界，营养会扩散到全域均匀分布。
3. 调整求解参数
   原代码时间步长乘了10的放大系数，容易出现数值振荡，建议去掉这个放大系数，同时适当增加模拟步数，就能完整看到扩散锋面从中心向外推进的全过程。

from fipy import CellVariable, Gmsh2D, TransientTerm, DiffusionTerm, Viewer
import numpy as np

# 网格基础参数
cellSize = 0.05
radius = 1.

# 生成圆形计算域网格
mesh = Gmsh2D('''
              cellSize = %(cellSize)g;
              radius = %(radius)g;
              Point(1) = {0, 0, 0, cellSize};
              Point(2) = {-radius, 0, 0, cellSize};
              Point(3) = {0, radius, 0, cellSize};
              Point(4) = {radius, 0, 0, cellSize};
              Point(5) = {0, -radius, 0, cellSize};
              Circle(6) = {2, 1, 3};
              Circle(7) = {3, 1, 4};
              Circle(8) = {4, 1, 5};
              Circle(9) = {5, 1, 2};
              Line Loop(10) = {6, 7, 8, 9};
              Plane Surface(11) = {10};
              ''' % locals())

# 初始化营养浓度场
phi = CellVariable(name="营养物质浓度",
                   mesh=mesh,
                   value=0.) 

# 设置圆心初始高浓度区
cell_center_dr = np.linalg.norm(mesh.cellCenters, axis=0)
phi.setValue(1., where=cell_center_dr < 0.05)

# 设置外边界条件：边缘浓度固定为0，可根据需求注释掉该行切换为零通量边界
phi.constrain(0., mesh.exteriorFaces)

# 初始化可视化窗口
viewer = None
if __name__ == '__main__':
    viewer = Viewer(vars=phi, datamin=0, datamax=1.)
    viewer.plotMesh()

# 定义纯扩散方程
D = 1.
eq = TransientTerm() == DiffusionTerm(coeff=D) 

# 调整时间步长保证数值稳定
timeStepDuration = 0.9 * cellSize**2 / (2 * D)
steps = 200

# 迭代求解
for step in range(steps):
    eq.solve(var=phi, dt=timeStepDuration)
    if viewer is not None and step % 5 == 0:
        viewer.plot()

• 初始圆心供给区的半径可以根据你模拟的实际血管尺寸调整，不用固定为0.05
• 如果要更贴合肿瘤微环境真实场景，可以在方程中添加营养消耗项、空间异质性扩散系数，对应补充ImplicitSourceTerm或修改DiffusionTerm的系数即可
• 如果需要输出定量结果，可以在迭代过程中按步提取不同半径位置的浓度值，绘制浓度随时间变化的径向分布曲线

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