含线性吸附与降解的一维MIM模型隐式有限差分MATLAB实现咨询

AI 小北

2026-6-1

多孔介质溶质运移MIM模型全隐式有限差分实现问题

模型背景与控制方程

针对含线性吸附与一阶降解的一维流动-不动水体模型（MIM），控制方程如下：

流动区域

$$\theta_m R_m \frac{\partial C_m}{\partial t} + \theta_{im} R_{im} \frac{\partial C_{im}}{\partial t} = \theta_m D_m \frac{\partial^2 C_m}{\partial x^2} - q \frac{\partial C_m}{\partial x} - \mu_m \theta_m C_m - \mu_{im} \theta_{im} C_{im}$$

不动区域

$$\theta_{im} R_{im} \frac{\partial C_{im}}{\partial t} = \alpha (C_m - C_{im}) - \mu_{im} \theta_{im} C_{im}$$

参数说明

$C_m$/$C_{im}$：流动区/不动区溶质浓度
$R_m$/$R_{im}$：流动区/不动区延迟因子
$\theta_m$/$θ_{im}$：流动区/不动区体积含水量
$\mu_m$/$μ_{im}$：流动区/不动区一阶降解速率
$\alpha$：流动-不动区质量传递系数

目标为实现全隐式有限差分格式（时间向后欧拉、弥散中心差分、对流迎风格式），保证无条件数值稳定性。

核心疑问

是否应先求解当前时间步的不动区方程，再代入流动区方程实现解耦？
还是构建更大的块对角稀疏矩阵，同步求解$t^{n+1}$时刻的$C_m$与$C_{im}$？若是，如何在MATLAB中用spdiags或sparse高效构建该耦合系统的索引结构？

解决方案建议

1. 解耦法：先求不动区浓度，再代入流动区方程

这种方法利用不动区方程的线性特性，先单独解出每个网格点的$C_{im}^{n+1}$，再代入流动区方程将其转化为标准对流弥散方程（ADE），复用你已熟悉的三对角矩阵求解逻辑，代码复杂度低，适合快速验证。

不动区方程离散与整理

采用向后欧拉离散时间项：
$$\theta_{im} R_{im} \frac{C_{im,i}^{n+1} - C_{im,i}^n}{\Delta t} = \alpha (C_{m,i}^{n+1} - C_{im,i}^{n+1}) - \mu_{im} \theta_{im} C_{im,i}^{n+1}$$

整理为$C_{im}^{{n+1}$关于$C_{m}}{n+1}$的线性表达式：
$$C_{im,i}^{n+1} = \frac{\theta_{im} R_{im} C_{im,i}^n + \alpha \Delta t C_{m,i}^{n+1}}{\theta_{im} R_{im} + \Delta t (\alpha + \mu_{im} \theta_{im})}$$

记为$C_{im,i}^{n+1} = a_i + b_i C_{m,i}^{n+1}$，其中：

$a_i = \frac{\theta_{im} R_{im} C_{im,i}^n}{DEN_i}$
$b_i = \frac{\alpha \Delta t}{DEN_i}$
$DEN_i = \theta_{im} R_{im} + \Delta t (\alpha + \mu_{im} \theta_{im})$

代入流动区方程

将$C_{im,i}^{{n+1}$的表达式代入流动区离散方程，即可得到仅含$C_{m,i}}{n+1}$的标准ADE形式，直接用你熟悉的三对角稀疏矩阵求解即可。

2. 耦合矩阵法：构建块三对角稀疏矩阵同步求解

适合后续扩展非线性模型的场景，矩阵为2x2块三对角结构，每个块对应单个网格点的$C_m$与$C_{im}$耦合关系。未知向量维度为$2N \times 1$，顺序为$[C_{m,1}, C_{im,1}, C_{m,2}, C_{im,2}, ..., C_{m,N}, C_{im,N}]^T$。

MATLAB构建核心示例

% 基础参数定义
N = 100;          % 网格数
dx = 1;           % 空间步长
dt = 0.1;         % 时间步长
theta_m = 0.3;    % 流动区含水量
theta_im = 0.1;   % 不动区含水量
R_m = 1.2;        % 流动区延迟因子
R_im = 1.0;       % 不动区延迟因子
D_m = 0.5;        % 弥散系数
q = 0.2;          % 达西流速
mu_m = 0.01;      % 流动区降解速率
mu_im = 0.005;    % 不动区降解速率
alpha = 0.02;     % 质量传递系数
C0 = 1.0;         % 左边界浓度

% 初始化上一步浓度（示例值）
C_m_old = zeros(N, 1);
C_im_old = zeros(N, 1);

% 稀疏矩阵构建准备
rows = []; cols = []; vals = [];
b = zeros(2*N, 1);

% 处理内部网格点（i=2到N-1）
for i = 2:N-1
    idx_m = 2*(i-1)+1;   % 当前网格C_m的全局索引
    idx_im = 2*(i-1)+2;  % 当前网格C_im的全局索引
    idx_m_left = idx_m - 2; % 左侧网格C_m的全局索引
    idx_m_right = idx_m + 2;% 右侧网格C_m的全局索引

    % 1. 不动区方程离散
    den_im = theta_im*R_im + dt*(alpha + mu_im*theta_im);
    rows = [rows, idx_im, idx_im];
    cols = [cols, idx_im, idx_m];
    vals = [vals, den_im, -alpha*dt];
    b(idx_im) = theta_im*R_im*C_im_old(i);

    % 2. 流动区方程离散（全隐式+迎风格式）
    % 时间项系数
    coeff_time_m = theta_m*R_m/dt;
    coeff_time_im = theta_im*R_im/dt;
    % 弥散项系数（中心差分）
    coeff_diff_center = 2*theta_m*D_m/(dx^2);
    coeff_diff_left = -theta_m*D_m/(dx^2);
    coeff_diff_right = -theta_m*D_m/(dx^2);
    % 对流项系数（迎风格式，q>0时取上游权重）
    coeff_conv_center = q/dx;
    coeff_conv_left = -q/dx;
    % 降解项系数
    coeff_decay_m = mu_m*theta_m;
    coeff_decay_im = mu_im*theta_im;

    % 流动区方程左边系数
    rows = [rows, idx_m, idx_m, idx_m, idx_m];
    cols = [cols, idx_m, idx_im, idx_m_left, idx_m_right];
    vals = [vals, coeff_time_m + coeff_diff_center + coeff_conv_center + coeff_decay_m, ...
            coeff_time_im + coeff_decay_im, ...
            coeff_diff_left + coeff_conv_left, ...
            coeff_diff_right];
    % 流动区方程右端项
    b(idx_m) = coeff_time_m*C_m_old(i) + coeff_time_im*C_im_old(i);
end

% 处理左边界（第一类边界C_m=C0）
rows = [rows, 1]; cols = [cols, 1]; vals = [vals, 1];
b(1) = C0;
% 左边界不动区方程
idx_im = 2;
den_im = theta_im*R_im + dt*(alpha + mu_im*theta_im);
rows = [rows, idx_im, idx_im];
cols = [cols, idx_im, 1];
vals = [vals, den_im, -alpha*dt];
b(idx_im) = theta_im*R_im*C_im_old(1);

% 处理右边界（零梯度Neumann边界）
idx_m = 2*(N-1)+1;
rows = [rows, idx_m, idx_m]; cols = [cols, idx_m, idx_m-2];
vals = [vals, 1, -1];
b(idx_m) = 0;
% 右边界不动区方程
idx_im = 2*(N-1)+2;
den_im = theta_im*R_im + dt*(alpha + mu_im*theta_im);
rows = [rows, idx_im, idx_im]; cols = [cols, idx_im, idx_m];
vals = [vals, den_im, -alpha*dt];
b(idx_im) = theta_im*R_im*C_im_old(N);

% 构建稀疏矩阵并求解
A = sparse(rows, cols, vals, 2*N, 2*N);
sol = A\b;

% 提取当前时间步浓度
C_m_new = sol(1:2:end);
C_im_new = sol(2:2:end);

两种方法对比

方法	优点	缺点
解耦法	代码简单、复用ADE逻辑、计算量小	仅适用于线性耦合场景，扩展受限
耦合矩阵法	结构清晰、支持非线性扩展	代码复杂度稍高，需处理块矩阵索引

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

含线性吸附与降解的一维MIM模型隐式有限差分MATLAB实现咨询

模型背景与控制方程 #

流动区域 #

不动区域 #