FiPy中能否仅在外边界添加源项？1D扩散问题技术问询

问题分析 #

你尝试的两种方法都无法精准在边界添加源/汇项：

• sourceTerm.constrain(sink, mesh.facesLeft)：constrain是用于约束变量（如浓度c）的边界值，而非源项。源项sourceTerm是单元格变量，不支持面约束。
• sourceTerm[0] = sink：这是将第一个单元格内部设为源/汇，并非真正的边界位置。

FiPy中，边界的源/汇项需要通过边界专属的项来实现，因为这类源汇本质是边界上的通量贡献，而非单元格内部的体积源。

解决方案1：使用BoundaryTerm（推荐） #

BoundaryTerm专门用于处理边界上的源/汇，直接作用于指定的边界面，精准对应物理上的边界位置。

修改后的代码：

L = 20                            # 厚度
nx = 100                          # 网格点数
dx = L/nx                         # 网格间距
mesh = Grid1D(nx = nx, dx = dx)   # 1D网格

sink = -10  # 汇项，负号表示物质消耗

c = CellVariable(mesh = mesh, name = r'$c_i$')
c.constrain(100, mesh.facesRight)  # 右边界Dirichlet条件

diffusionCoefficient = CellVariable(mesh = mesh, value = 5, name = r'$D_i$')

# 左边界添加精准汇项：BoundaryTerm作用于左边界面
eq = DiffusionTerm(coeff = -diffusionCoefficient, var = c) == BoundaryTerm(
    coeff=sink, 
    var=c, 
    faces=mesh.facesLeft
)

eq.solve()

viewer = Matplotlib1DViewer(vars = c)
viewer.plot()
input('Press return to continue')

原理说明

• BoundaryTerm(coeff=sink, faces=mesh.facesLeft)直接将汇项绑定到左边界的面上，完全对应物理边界位置。
• 原方程-∇·(D∇c) = S中，边界汇项通过BoundaryTerm作为方程右端项，等价于在边界处施加了物质消耗的通量条件。

解决方案2：使用固定通量边界条件 #

边界汇项本质等价于边界处的固定通量，因此也可以通过设置面梯度约束实现：

L = 20                            # 厚度
nx = 100                          # 网格点数
dx = L/nx                         # 网格间距
mesh = Grid1D(nx = nx, dx = dx)   # 1D网格

sink = -10  # 汇项，负号表示物质消耗

c = CellVariable(mesh = mesh, name = r'$c_i$')
c.constrain(100, mesh.facesRight)  # 右边界Dirichlet条件

diffusionCoefficient = CellVariable(mesh = mesh, value = 5, name = r'$D_i$')

# 左边界设置固定通量：通量J = -D dc/dx = sink
c.faceGrad.constrain(
    sink / (-diffusionCoefficient), 
    mesh.facesLeft
)

# 内部无源性，方程右端为0
eq = DiffusionTerm(coeff = -diffusionCoefficient, var = c) == 0

eq.solve()

viewer = Matplotlib1DViewer(vars = c)
viewer.plot()
input('Press return to continue')

原理说明

• 汇项对应边界处的物质通量，通过约束浓度的面梯度faceGrad，间接实现边界处的物质消耗，效果与BoundaryTerm完全一致。

两相耦合场景适配 #

对于你提到的两相质量守恒需求，只需在另一相的网格右边界添加符号相反的BoundaryTerm即可：

# 另一相的右边界添加源项（与当前相左边界汇项符号相反）
eq_phase_j = DiffusionTerm(coeff=-D_j, var=c_j) == BoundaryTerm(coeff=-sink, var=c_j, faces=mesh_j.facesRight)

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