考虑非线性扩散系数的Richards方程1D模拟结果偏差问题求助

我正在开展Richards方程的一维土柱渗流模拟：底部（一维左侧）设置Dirichlet边界（u=0，水压为0），顶部（一维右侧）设置Neumann边界（降雨q=1E-5）。扩散系数（导水率k(u)）是依赖负压的非线性函数。我用Fipy编写了迭代更新k(u)的脚本，但计算结果和专业渗流软件SEEP/W的结果偏差很大，求解决方法。

原代码 #

import fipy as fp
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Simulation with professional software SEEP/W
seep = pd.read_csv('https://gist.githubusercontent.com/essicolo/e670981aa0171a449a554d7763157374/raw/8d06b07117f6f31a57b2156d07a2ba971a355e97/gistfile1.txt')

# Material parameters
alpha = 2.0 / 9.807  # [1/m]
n = 1.5
m = 1 - 1/n
l = 0.5
Ks = 1.0E-4  # [m/s]

# Top infiltration
rain = 5E-5

# Domain
nz = 200
Lz = 20.0
dz = Lz / nz
mesh = fp.Grid1D(nx=nz, dx=dz)

# Variables
z = fp.CellVariable(name="elevation", mesh=mesh, value=mesh.cellCenters[0])
u = fp.CellVariable(name="matrix pressure", mesh=mesh, value=-z.value)
u.setValue(u_min, where = z >= -u_min) # close to the real solution
u_old = fp.CellVariable(name="old matrix pressure", mesh=mesh, value=-u.value)
u_convergence = fp.CellVariable(name="convergence check", mesh=mesh, value=0.)
K = fp.CellVariable(mesh=mesh, value=Ks)

# Boundary conditions
u.constrain(0., mesh.facesLeft) # bottom
g = fp.FaceVariable(mesh=mesh, value=-rain)
u.faceGrad.constrain(g * mesh.faceNormals, mesh.facesRight) # top

# Solver
max_iters = 100
tolerance = 1e-10
solver = fp.LinearLUSolver()

convergence = []
iterations = []
for iter in range(max_iters):

    plt.plot(u.value, z.value, label=f'Iter {iter}')
    u_old.value = u.value

    K.setValue(Ks * (1 - (-alpha * u)**(n-1) * (1 + (-alpha * u)**n)**(-m)) * (1 + (-alpha * u)**n)**(-m*l), where=u < 0)
    K.setValue(Ks, where=u >= 0)
    
    eq = (fp.DiffusionTerm(coeff=K, var=u) + K == 0)
    eq.solve(var=u, solver=solver)

    u_convergence.value = np.abs(u.value - u_old.value)

    iterations.append(iter)
    convergence.append(np.mean(u_convergence.value))
    
    if np.all(u_convergence.value < tolerance):
        break


plt.plot(seep['u'], seep['z'], c='k', label='with SEEP/W')
plt.legend()

#plt.plot(iterations, convergence)

核心问题排查与修正方案 #

1. Richards方程形式不匹配 #

原代码构建的方程fp.DiffusionTerm(coeff=K, var=u) + K == 0存在推导误差。标准Richards方程（总水头形式）为：
$$\nabla \cdot \left(K(h)\nabla h\right) = 0$$
其中$h = u + z$（$u$为基质势，$z$为高程）。直接用基质势$u$构建方程时，需补充对流项修正高程影响，或转换为总水头变量简化方程：

# 转换为总水头h计算
h = u + z
eq = fp.DiffusionTerm(coeff=K, var=h) == 0

2. 导水率公式实现错误 #

Van Genuchten-Mualem模型的标准导水率公式中，饱和比$S_e$的部分是平方项，原代码仅为一次方，这是结果偏差的核心原因：

# 修正导水率计算
mask = u < 0
se = (1 - (-alpha * u[mask])**(n-1) * (1 + (-alpha * u[mask])**n)**(-m))**2
kr = (1 + (-alpha * u[mask])**n)**(-m*l)
K[mask] = Ks * se * kr
K[u >= 0] = Ks

3. 边界条件设置错误 #

顶部Neumann边界的通量约束逻辑错误，应直接使用FixedFluxBoundary定义降雨通量，避免梯度计算的方向混淆：

boundary_top = fp.FixedFluxBoundary(mesh.facesRight, value=rain)
# 求解时传入边界条件
eq.solve(var=u, solver=solver, boundaryConditions=[boundary_top])

4. 初始值与收敛准则问题 #

• 原代码中u_min未定义，需补充初始值截断阈值（如u_min = -10）
• 收敛准则np.all(u_convergence.value < tolerance)过于严苛，改用最大绝对误差判断更合理：

max_diff = np.max(np.abs(u.value - u_old.value))
if max_diff < 1e-6:
    break

修正后完整代码 #

import fipy as fp
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# SEEP/W模拟结果
seep = pd.read_csv('gistfile1.txt')

# 材料参数
alpha = 2.0 / 9.807  # [1/m]
n = 1.5
m = 1 - 1/n
l = 0.5
Ks = 1.0E-4  # [m/s]

# 顶部降雨强度
rain = 5E-5  # m/s

# 计算域
nz = 200
Lz = 20.0
dz = Lz / nz
mesh = fp.Grid1D(nx=nz, dx=dz)

# 变量定义
z = fp.CellVariable(name="高程", mesh=mesh, value=mesh.cellCenters[0])
u = fp.CellVariable(name="基质势", mesh=mesh, value=-z.value)
u_min = -10
u.setValue(u_min, where=z >= -u_min)
u_old = fp.CellVariable(name="旧基质势", mesh=mesh, value=u.value)
K = fp.CellVariable(mesh=mesh, value=Ks)

# 边界条件
u.constrain(0., mesh.facesLeft)  # 底部Dirichlet边界
boundary_top = fp.FixedFluxBoundary(mesh.facesRight, value=rain)

# 求解器设置
max_iters = 100
tolerance = 1e-6
solver = fp.LinearLUSolver()

convergence = []
iterations = []

for iter in range(max_iters):
    u_old.value = u.value
    
    # 修正导水率计算
    mask = u < 0
    se = (1 - (-alpha * u[mask])**(n-1) * (1 + (-alpha * u[mask])**n)**(-m))**2
    kr = (1 + (-alpha * u[mask])**n)**(-m*l)
    K[mask] = Ks * se * kr
    K[u >= 0] = Ks
    
    # 总水头形式的Richards方程
    h = u + z
    eq = fp.DiffusionTerm(coeff=K, var=h) == 0
    
    eq.solve(var=h, solver=solver, boundaryConditions=[boundary_top])
    u.value = h.value - z.value
    
    # 收敛检查
    max_diff = np.max(np.abs(u.value - u_old.value))
    convergence.append(max_diff)
    iterations.append(iter)
    
    if max_diff < tolerance:
        print(f"迭代收敛，迭代次数：{iter+1}")
        break

# 绘图对比
plt.figure(figsize=(8, 10))
plt.plot(u.value, z.value, label='Fipy模拟结果')
plt.plot(seep['u'], seep['z'], c='k', linestyle='--', label='SEEP/W模拟结果')
plt.xlabel('基质势 u (m)')
plt.ylabel('高程 z (m)')
plt.legend()
plt.grid(True)
plt.show()

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