简介二维渗漏承压含水层流动方程的数值求解是地下水动力学教学与科研中的常见问题这份资源基于MATLAB 2019a实现采用有限差分法FDM结合高斯-赛德尔迭代解算器对描述承压含水层渗漏的泊松方程进行离散与迭代求解适合本科、硕士阶段的地下水数值模拟课程设计及入门科研参考。压缩包共4个文件包含1个可直接运行的.m主程序与3个结果展示图片代码结构简洁便于对照图像分析渗流场分布。资源包大小仅31KB轻量易用目前已有185人学习。借助该代码读者可快速掌握泊松方程在含水层渗漏问题中的建模思路、FDM网格剖分方法以及高斯-赛德尔迭代的收敛过程亦可在此基础上扩展边界条件或修改源汇项用于更复杂的承压流动分析。1. 承压含水层二维渗漏流从Laplace到Poisson的跃迁如果只考虑水平二维展布、顶底板完全隔水的承压含水层稳态水头满足Laplace方程这是地下水数值模拟课程里最熟悉的起点。可一旦含水层顶板是弱透水层、上方又有持续补给水头分布就会出现“渗漏”效应方程右端不再是零而是带上越流补给项的Poisson型方程∇²H -R/T。这个差别直接影响有限差分FDM离散格式Laplace方程的系数矩阵是严格的五对角负定结构而Poisson方程引入右端源项后虽然仍可用迭代解算器求解但边界处理和收敛行为都变了。这个项目就是拿MATLAB 2019a把网格剖分、五点差分、高斯-赛德尔迭代串成一条能跑通的链路适合本科毕业设计、硕士课程作业里需要“从公式到代码”完整走一遍的人也适合后续改用稀疏矩阵直接求解或加入SOR加速的团队做参照。2. FDM五点差分与承压含水层网格离散理论到代码之间最重要的一步是把连续方程变成离散代数系统。我一般先写物理参数和网格再给出五点差分公式最后处理边界条件顺序不能反因为边界格式取决于离散格式。2.1 含水层渗漏方程与越流因子的定义承压含水层稳态二维渗漏流的常用控制方程为∂/∂x (T ∂H/∂x) ∂/∂y (T ∂H/∂y) (k / b) (H0 - H) 0T是导水系数b是弱透水层厚度k是弱透水层垂向渗透系数H0是上覆含水层或地表水体的固定水头。当T在空间上均匀时方程整理为∂²H/∂x² ∂²H/∂y² (H0 - H) / λ²其中λ² T b/k量纲是长度平方代表渗漏影响的特征尺度。λ越小越流补给越强水头越被拉向H0λ越大方程越接近Laplace方程。很多误用来自直接把右端设成常数但越流项实际上是H的函数只有把它当作“上一轮迭代已知的源项”处理才能在泊松框架里收敛。下面的表格给出本代码默认参数的量级这些值是教研场景里常用的“中等尺度含水层”参数不是野外反演结果但足以保证迭代稳定。参数含义默认值Lx, Ly矩形区域长宽m1000, 1000nx, nyx/y方向节点数101, 101T导水系数m²/s1e-4b_prime弱透水层厚度m10k_leak弱透水层垂向渗透系数m/s1e-8H0越流源水头m50H_left, H_right左右边界水头m40, 35网格剖分我习惯用均匀网格。矩形区域长Lx、宽Lyx方向节点数nxy方向节点数ny则Δx Lx/(nx-1)Δy Ly/(ny-1)。节点编号采用行优先外层循环j从1到ny内层循环i从1到nx这样MATLAB访问连续内存更快后续用reshape和surf绘图也方便。2.2 五点差分格式与迭代式推导对内部节点(i,j)二阶中心差分格式为∂²H/∂x² ≈ (H(i1,j) - 2H(i,j) H(i-1,j)) / Δx² ∂²H/∂y² ≈ (H(i,j1) - 2H(i,j) H(i,j-1)) / Δy²代入带越流的控制方程得到离散关系(H(i1,j) - 2H(i,j) H(i-1,j)) / Δx² (H(i,j1) - 2H(i,j) H(i,j-1)) / Δy² (H0 - H(i,j)) / λ²整理成高斯-赛德尔可用的显式更新式H(i,j) [ (H(i1,j)H(i-1,j))/Δx² (H(i,j1)H(i,j-1))/Δy² H0/λ² ] / [ 2/Δx² 2/Δy² 1/λ² ]这个格式天然满足对角占优迭代不需要矩阵分解也没有选主元的问题。缺点是收敛速度随节点数增加而下降nx从51加到201迭代次数大约翻倍。实际代码中不用构造系数矩阵直接用这个表达式在网格上扫描即可。需要注意的是如果Δx远大于Δyy方向二阶项的权重更大会放大边界误差所以我倾向于让Δx和Δy尽量接近避免强迫收敛。2.3 边界条件处理虚节点法与固定水头承压含水层常见的定解条件有三类第一类Dirichlet如河流切割含水层给定水头第二类Neumann如隔水边界给定零通量第三类混合边界。这个项目只实现左右Dirichlet、上下Neumann因为最能体现渗漏项的影响也容易验证。上边界j1为零通量时∂H/∂y0中心差分需要虚节点H(i,0)H(i,2)。代入五点格式后边界节点离散式变成(H(i1,1)-2H(i,1)H(i-1,1))/Δx² 2(H(i,2)-H(i,1))/Δy² (H0 - H(i,1))/λ²也就是说Neumann边界不新增未知数而是调整y方向二阶差的系数。千万不要把边界节点当内部节点直接套原格式那等于强制给边界加了一个错误的第二类条件迭代会在边界附近震荡。下面用文本示意网格节点类型黑色代表固定水头灰色代表隔水边界。j ny ┌─────────────────────────────┐ Neumann │ · · · · · · │ │ · · · · · · │ j 1 │ · · · · · · │ └─────────────────────────────┘ Neumann i 1 i 2 ... i nx3. 高斯-赛德尔迭代解算器的MATLAB实现理论离散清楚后这一章把公式变成可运行的MATLAB函数。核心文件是leaky_semi_confined_steady_state.m封装成函数而不是脚本方便后续批量扫描参数和输出收敛历史。3.1 为什么不构造稀疏矩阵高斯-赛德尔迭代的教材写法是先把系数矩阵拆成A D - L - U然后套H_new D^{-1}(LU)H_new D^{-1}b。很多同学会先用sparse构造五对角矩阵再在循环里做矩阵向量乘法。这个思路没有错但在这个问题里系数矩阵结构非常规则只有主对角线和四个偏移量显式构造矩阵反而增加内存和出错概率。我一般直接写标量更新式。高斯-赛德尔与雅可比的关键区别是计算H(i,j)时H(i-1,j)和H(i,j-1)已经来自本轮迭代而H(i1,j)和H(i,j1)还是上一轮值。这种就地更新使内存占用只有两套网格且收敛速度大约是雅可比的两倍。扫描次序会影响收敛路径但最终收敛解唯一。如果希望更快可以改成红黑排序或SOR后面会提到。3.2 完整可运行的MATLAB 2019a函数下面的函数实现了全部流程。为了兼容MATLAB 2019a没有使用arguments块和string数组只用传统语法。输入是一个结构体A输出含水头矩阵、坐标网格和残差历史。function [H, X, Y, residual_history] leaky_semi_confined_steady_state(A) % 求解二维渗漏承压含水层稳态水头分布 % 方程d2H/dx2 d2H/dy2 (H0 - H) / (T * b_prime / k_leak) % 边界左右Dirichlet上下Neumann % 输入A为结构体至少包含以下字段 % nx, ny, Lx, Ly, T, k_leak, b_prime, H0, H_left, H_right, tol, max_iter % 可选字段omega默认1.0即高斯-赛德尔 nx A.nx; ny A.ny; Lx A.Lx; Ly A.Ly; T A.T; k_leak A.k_leak; b_prime A.b_prime; H0 A.H0; H_left A.H_left; H_right A.H_right; tol A.tol; max_iter A.max_iter; if ~isfield(A, omega) || isempty(A.omega) omega 1.0; else omega A.omega; end dx Lx / (nx - 1); dy Ly / (ny - 1); lambda2 T * b_prime / k_leak; % 越流因子单位m^2 if lambda2 0 error(lambda2必须为正数请检查T/b_prime/k_leak的量级); end % 坐标网格 x linspace(0, Lx, nx); y linspace(0, Ly, ny); [X, Y] meshgrid(x, y); % 初始水头线性插值左到右作为收敛初场 H zeros(ny, nx); for j 1:ny H(j,:) H_left (H_right - H_left) * (x / Lx); end % 左右固定水头边界 H(:,1) H_left; H(:,nx) H_right; % 系数缓存避免循环内重复计算 inv_dx2 1 / dx^2; inv_dy2 1 / dy^2; inv_lambda2 1 / lambda2; denominator 2 * inv_dx2 2 * inv_dy2 inv_lambda2; residual_history zeros(1, max_iter); H_old zeros(ny, nx); for iter 1:max_iter H_old H; % 保存上一轮水头用于计算收敛增量和残差曲线 % 上下边界Neumann用虚节点法更新 for i 2:nx-1 H_new_top ( (H(1,i1)H(1,i-1))*inv_dx2 2*H(2,i)*inv_dy2 H0*inv_lambda2 ) / denominator; H(1,i) (1 - omega) * H(1,i) omega * H_new_top; H_new_bottom ( (H(ny,i1)H(ny,i-1))*inv_dx2 2*H(ny-1,i)*inv_dy2 H0*inv_lambda2 ) / denominator; H(ny,i) (1 - omega) * H(ny,i) omega * H_new_bottom; end % 内部节点统一更新 for j 2:ny-1 for i 2:nx-1 H_gs ( (H(j,i1)H(j,i-1))*inv_dx2 (H(j1,i)H(j-1,i))*inv_dy2 H0*inv_lambda2 ) / denominator; H(j,i) (1 - omega) * H(j,i) omega * H_gs; end end % 计算收敛判据取最大绝对增量 residual_history(iter) max(max(abs(H - H_old))); if residual_history(iter) tol residual_history residual_history(1:iter); fprintf(收敛于第%d次迭代残差%.3e\n, iter, residual_history(iter)); return; end end warning(达到最大迭代次数%d未收敛到tol%.1e请检查网格步长或松弛因子, max_iter, tol); end这段代码的逻辑分四步从结构体A读出参数计算网格步长、越流因子λ²和denominator。lambda2是T*b_prime/k_leak单位m²它直接决定右端源项的强度。用linspace和meshgrid建立坐标。初始场用左右边界线性插值而不是全零这样在强越流时能显著减少前几十次迭代的冲量。边界更新和内部更新分开。上下边界使用Neumann公式内部使用五点格式。所有重复系数先算成inv_dx2、inv_dy2、inv_lambda2避免内循环做除法。每轮末尾比较H与H_old的最大绝对增量并保存到residual_history。注意这里H_old是在边界更新前保存的所以它代表上一轮的完整水头场。参数方面tol通常取1e-6到1e-8单位是米。如果只是教研演示1e-6足够如果要做定量对比建议1e-9。max_iter不要小于10000因为101×101网格下GS默认收敛大约需要几千次。3.3 调用示例与初值敏感性调用代码非常简单把参数塞进结构体后执行一次函数即可A.nx 101; A.ny 101; A.Lx 1000; A.Ly 1000; A.T 1e-4; A.k_leak 1e-8; A.b_prime 10; A.H0 50; A.H_left 45; A.H_right 35; A.tol 1e-6; A.max_iter 20000; [H, X, Y, res] leaky_semi_confined_steady_state(A); figure; contourf(X, Y, H, 20); colorbar; xlabel(x (m)); ylabel(y (m)); title(Steady-state head in leaky confined aquifer);如果看三维水头面把contourf换成surf(X,Y,H,EdgeColor,none)再加view(2)。初始场敏感性方面线性插值初值在λ²较小时能比全零初值少5%到10%的迭代次数但最终结果一致说明解的唯一性没有被破坏。4. 收敛判据、松弛因子与MATLAB排错实践这一章专门解决“为什么迭代卡住、震荡、慢得像爬”的问题。GS迭代本身简单但实际运行中一半问题出在判据和边界另一半出在参数量级不匹配。4.1 绝对残差与相对残差怎么选上面代码使用最大绝对增量max(|H_new - H_old|)。这个判据直观但有一个陷阱当水头以米计时1e-6的绝对残差有物理意义如果换成毫米计1e-6就太严苛。我一般同时计算相对残差让判据无量纲化relative_res norm(H(:) - H_old(:), 2) / norm(H(:), 2); if relative_res rel_tol % 达到相对精度要求提前退出 end用norm函数每轮只调用一次不影响主循环。工程上可以设置双条件绝对残差小于tol或相对残差小于rel_tol满足一个就退出。下表对比了三种常用判据的适用场景。判据优点缺点max abs diff实现简单、物理直观受水头绝对量级影响大L2相对残差无量纲适合不同单位体系需要额外norm计算能量范数可直接量化误差需要解析解不通用对于课程作业单用绝对残差没问题但报告里要写清为什么选这个阈值。对于我自己的项目我会把tol取1e-8rel_tol取1e-10两者取先到者。4.2 松弛因子与SOR加速的实现高斯-赛德尔是SOR在ω1时的特例。当节点数超过151×151后GS收敛速度明显下降常见做法是加松弛因子。更新公式变为H_gs ( (H(j,i1)H(j,i-1))*inv_dx2 (H(j1,i)H(j-1,i))*inv_dy2 H0*inv_lambda2 ) / denominator; H(j,i) (1 - omega) * H(j,i) omega * H_gs;ω取值通常在1.0到1.7之间。对均匀矩形网格最优ω可以通过红黑扫描求出但工程上有个简单办法先跑三组ω1.0、1.2、1.4比较达到tol的迭代次数选最快的一组。下面这段代码演示批量对比omega_list [1.0, 1.2, 1.4, 1.6]; iter_count zeros(size(omega_list)); for k 1:length(omega_list) A.omega omega_list(k); [~, ~, ~, res] leaky_semi_confined_steady_state(A); iter_count(k) length(res); end fprintf(omega%s 时迭代次数分别为 %s\n, ... mat2str(omega_list), mat2str(iter_count));在101×101网格、λ≈316m的默认参数下ω1.4通常能比ω1.0快约30%。网格越密最优ω越接近2但超过1.8后数值溢出风险陡增我一般不推荐超过1.7。提示残差历史曲线建议用semilogy绘制线性坐标下小残差会被压成一条直线无法区分停滞和缓慢下降。4.3 边界条件引发的震荡和发散最常见的振荡源是Neumann边界与Dirichlet边界衔接的角落点。在我的实现里四个角落属于Dirichlet因为H(:,1)和H(:,nx)先被赋成固定值而上下边界更新时i从2开始不会覆盖角落所以一致性没问题。如果同学把全边界都当Dirichlet赋同一个值角落会被两层循环反复覆盖最终结果取决于最后一次写入非常隐蔽。第二个陷阱是上下边界公式里的2*H(2,i)*inv_dy2。它基于均匀网格的虚节点只有在Δx和Δy各自均匀时成立。如果使用非均匀网格这个系数要按实际距离重新推导不能直接套。第三个陷阱是λ²太小。λ²小于四个网格步长平方时右端源项主导显式迭代会表现出类似“刚性”的振荡。对策是加密网格或改用SOR但最稳妥的是检查量级lambda2 4*max(dx^2, dy^2)。4.4 用残差历史定位问题把residual_history画出来能够快速区分三种情况残差持续下降但极慢网格太密或ω偏小对策是换SOR或先用粗网格试算。残差先降后升λ²过小越流强度过大对策是检查量级并加密网格。残差来回震荡边界条件设置矛盾或初始场不连续对策是打印角落和边界值。画残差曲线的代码figure; semilogy(res, o-); xlabel(迭代次数); ylabel(max |H^{k1} - H^k|); grid on;如果曲线像一条平线说明达到迭代上限如果出现周期性波动优先查边界。5. 算例验证与越流参数敏感性分析这一章给出两个可以写进课设报告的进阶内容一是解析解对照验证二是越流强度对水头形态的影响。5.1 与解析解做快速对拍最简单有效的验证是把k_leak设为0让方程退化为Laplace方程。此时上下Neumann边界下解应该只沿x方向线性变化即H(x) H_left (H_right - H_left)*x/Lx。用代码算一次与线性解析解的最大误差应小于1e-10。这一步能过滤掉八成离散错误。如果想保留越流项可以让y方向只设3个节点上下边界都改成Dirichlet并赋同一解析解然后看二维解在每一列上是否重合。一维越流解析解为H(x) H0 C1*cosh(x/λ) C2*sinh(x/λ)其中C1和C2由边界条件H(0)H_left、H(Lx)H_right确定。MATLAB里可以用\求解2×2线性方程组避免手推公式出错。5.2 越流因子λ²对水头形态的影响λ² T b/k它把导水系数和弱透水层参数缩成一个长度尺度。λ越小越流越强水头越被拉向H0λ越大水头越接近无渗漏的线性分布。下面的扫描循环可以画出中心点水头随λ的变化A.nx 81; A.ny 81; A.Lx 1000; A.Ly 1000; A.H0 50; A.H_left 45; A.H_right 35; A.tol 1e-7; A.max_iter 20000; lambda_list [50, 150, 300, 600, 1200]; center_head zeros(size(lambda_list)); figure; hold on; for k 1:length(lambda_list) A.T 1e-4; A.b_prime 10; A.k_leak A.T * A.b_prime / lambda_list(k)^2; % 由lambda反推k_leak [H, X, Y, ~] leaky_semi_confined_steady_state(A); center_head(k) H(round(A.ny/2), round(A.nx/2)); end plot(lambda_list, center_head, k-o); xlabel(lambda (m)); ylabel(中心点水头 (m));结果会显示λ从1200降到50中心水头从接近线性插值的40.1 m升高到接近48 m。这个趋势可以解释为渗漏强度增加时越流源水头H0对承压含水层的主导作用增强。写报告时可以直接把中心点曲线和两幅水头等值面截图放在一起对比比纯文字更有说服力。5.3 把脚本改造成通用求解器的小技巧最后一个实用技巧不要把所有参数都堆在函数参数列表里而是用一个结构体承载。教研场景经常要跑几十组参数对比如果函数签名是(nx, ny, Lx, Ly, T, b_prime, k_leak, H0, ...)调用一个参数就要数一遍顺序极易错位。改成结构体后批量扫描只需要在循环里更新相应字段。另一个技巧是把H_old和residual_history预分配并把残差历史作为输出返回这样后续画图、写表都不用重跑代码。调试时如果怀疑某个区域不收敛可以临时加一行if jny/2 inx/2, fprintf(%d %.6f\n, iter, H(j,i)); end观察中心点水头是否在正确区间内变化。这样做能快速区分迭代算法问题还是边界条件问题。本文还有配套的精品资源点击获取