简介LearnEDFM是一套基于MATLAB的开源数值程序采用嵌入式离散裂缝方法EDFM求解裂隙多孔介质中的流动与运移问题面向地下模拟、油藏工程及裂缝性介质研究领域的研究人员与研究生适合具备一定数值方法与MATLAB基础、希望快速上手并应用EDFM的学习者。资源包共111个文件约5.15MB以43个.m源码文件为核心配套19个png与14个jpg示意图、若干js/html/css网页文档、pdf用户手册及mat数据文件另含license与md说明结构完整便于查阅与二次开发。压缩包内提供examples示例目录运行ex1即可测试程序并附有从GitHub仓库克隆构建的说明方便读者对照源码理解裂缝建模、流动与运移耦合求解流程。目前已有504人学习下载可作为学习EDFM原理、复现算例与开展裂隙介质模拟研究的实用参考。1. 从裂隙介质到求解器为什么需要自己写 EDFM 求导代码裂隙多孔介质里的流动和运移问题商业软件能算但一旦你要改本构、换耦合方式、加自定义的传质项闭源求解器就成了黑匣子。嵌入式离散断裂方法EDFM的思路很直接不把裂缝显式剖成网格而是把它作为低维单元“嵌”进基质网格里通过非相邻连接NNC计算基质与裂缝、裂缝与裂缝之间的流量。这样做的好处是网格量不随裂缝数量爆炸坏处是你得自己推导并实现那些连接对的传导率还要处理时间项、对流项和扩散项的离散。标题里的“matlab求导代码”其实指向一个很具体的需求在 EDFM 框架下流动和运移方程离散后会产生大量对压力、饱和度、浓度等变量的偏导数也就是雅可比矩阵的元素。这些导数如果推错一个符号牛顿迭代就崩。所以一套可复现的 EDFM 数值码核心不只是组装矩阵还包括怎么用 MATLAB 把解析导数写对、怎么用数值差分做验证。这套东西适合两类人一是做裂隙介质数值模拟的研究生和工程师二是想从零理解 EDFM 而不是只会点按钮的从业者。下面按“概念—实现—避坑—进阶”的顺序拆开讲。2. EDFM 的几何降维与连接对计算从三维裂缝到一维单元2.1 为什么裂缝要降维降维后丢了什么三维基质网格里如果按真实开度剖分裂缝网格尺度会被开度通常微米到毫米级拖死而储层尺度是米到千米级跨尺度直接剖分不现实。EDFM 的做法是把裂缝当成 (n-1) 维单元三维问题里裂缝是二维面二维问题里裂缝是一维线。裂缝单元有自己的开度、渗透率和孔隙度但它不占基质网格的体积。基质和裂缝之间的流动通过“基质网格—裂缝单元”连接对来算裂缝之间的相交则通过“裂缝单元—裂缝单元”连接对来算。降维丢掉的首先是裂缝内部的横向压力梯度细节因为一个裂缝单元内假设压力均匀其次是裂缝开度变化对基质网格几何的扰动EDFM 假设裂缝嵌在基质里但不改变基质网格的几何。这些假设在裂缝开度远小于网格尺度时成立如果开度和网格尺度可比EDFM 的精度就会下降这时候要么加密网格要么换显式剖分。2.2 三类连接对的传导率公式与 MATLAB 实现EDFM 里需要计算三类连接基质—基质MM、基质—裂缝MF、裂缝—裂缝FF。MM 连接就是常规相邻网格的传导率MF 和 FF 是 EDFM 特有的。下面这段 MATLAB 代码计算一个二维问题里基质网格与裂缝单元之间的传导率假设裂缝是直线段基质是矩形网格。function T_MF calc_T_MF(mat_grid, frac_elem, k_mat, mu) % mat_grid: 基质网格结构体含 xmin,xmax,ymin,ymax,dx,dy % frac_elem: 裂缝单元结构体含 x1,y1,x2,y2, aperture, perm % k_mat: 基质渗透率标量或矩阵 % mu: 流体粘度 % 返回基质网格与裂缝单元之间的传导率矩阵 nx mat_grid.nx; ny mat_grid.ny; n_frac length(frac_elem); T_MF sparse(nx*ny, n_frac); for i 1:nx for j 1:ny mat_id (j-1)*nx i; % 基质网格中心 xc mat_grid.xmin (i-0.5)*mat_grid.dx; yc mat_grid.ymin (j-0.5)*mat_grid.dy; for f 1:n_frac % 裂缝单元中点 xf 0.5*(frac_elem(f).x1 frac_elem(f).x2); yf 0.5*(frac_elem(f).y1 frac_elem(f).y2); % 判断裂缝是否穿过该基质网格 if is_frac_in_grid(frac_elem(f), mat_grid, i, j) % 裂缝在基质网格内的长度 L calc_frac_length_in_grid(frac_elem(f), mat_grid, i, j); % 基质到裂缝的平均距离常用等效半径近似 d sqrt(mat_grid.dx * mat_grid.dy / pi); % 传导率 渗透率 * 面积 / (粘度 * 距离) % 二维问题面积用长度 L 代替开度在裂缝侧已考虑 T_MF(mat_id, f) k_mat * L / (mu * d); end end end end end这段代码里is_frac_in_grid和calc_frac_length_in_grid需要自己实现前者判断裂缝线段是否与基质网格相交后者算交线长度。传导率公式k*L/(mu*d)里L是裂缝在基质网格内的长度d是基质到裂缝的特征距离常用等效半径sqrt(dx*dy/pi)近似。注意这里没有乘裂缝开度因为二维问题里裂缝单元本身是一维的开度在裂缝单元的孔隙体积和裂缝间连接里体现。如果是三维问题L换成面积d换成等效半径(dx*dy*dz)^(1/3)的某个倍数。FF 连接的传导率计算更麻烦因为两条裂缝相交时相交点附近的流动是三维的EDFM 用“相交线长度”或“相交点”来近似。常见做法是如果两条裂缝相交在相交处各生成一个虚拟单元然后用相交线长度乘上调和平均渗透率和开度除以距离。下面是一个简化版的 FF 传导率计算。function T_FF calc_T_FF(frac_elem, mu) % 计算裂缝单元之间的传导率假设裂缝单元已离散为线段 n_frac length(frac_elem); T_FF sparse(n_frac, n_frac); for i 1:n_frac for j i1:n_frac [flag, L_int] frac_intersect(frac_elem(i), frac_elem(j)); if flag % 调和平均渗透率 k_eq 2 * frac_elem(i).perm * frac_elem(j).perm / ... (frac_elem(i).perm frac_elem(j).perm); % 等效开度 w_eq 0.5*(frac_elem(i).aperture frac_elem(j).aperture); % 相交线长度 L_int二维问题里是相交点用 1 代替 T_FF(i,j) k_eq * w_eq * L_int / mu; T_FF(j,i) T_FF(i,j); end end end endfrac_intersect返回两条线段是否相交以及相交部分的长度二维里相交是点长度取 1 或特征长度。k_eq用调和平均是因为两条裂缝串联w_eq用算术平均是简化。实际写的时候FF 传导率对结果影响很大尤其是裂缝网络密集时建议用数值实验标定。2.3 连接对列表的组装与稀疏矩阵存储算完三类传导率后要把它们组装成全局连接对列表。MATLAB 里用稀疏矩阵最方便但 EDFM 的连接对不是简单的相邻关系需要先建立单元索引映射基质单元编号 1 到 Nm裂缝单元编号 Nm1 到 NmNf。然后构造一个N x N的稀疏矩阵非零元就是传导率。下面这段代码展示组装过程。function T assemble_T(mat_grid, frac_elem, k_mat, mu) Nm mat_grid.nx * mat_grid.ny; Nf length(frac_elem); N Nm Nf; % 基质-基质连接 T_MM calc_T_MM(mat_grid, k_mat, mu); % 基质-裂缝连接 T_MF calc_T_MF(mat_grid, frac_elem, k_mat, mu); % 裂缝-裂缝连接 T_FF calc_T_FF(frac_elem, mu); % 组装 T sparse(N, N); T(1:Nm, 1:Nm) T_MM; T(1:Nm, Nm1:N) T_MF; T(Nm1:N, 1:Nm) T_MF; T(Nm1:N, Nm1:N) T_FF; end组装时注意T_MF的转置因为基质到裂缝和裂缝到基质的传导率相同但方向相反。稀疏矩阵存储后后续求解线性系统直接用 MATLAB 的反斜杠或迭代求解器。这里有个坑如果裂缝单元数量大T_FF的循环会慢可以用向量化或 KD 树加速相交判断。3. 流动与运移方程的离散从传导率到雅可比矩阵3.1 单相流动的离散与解析导数先看最简单的单相不可压缩流动。控制方程是∇·(k/μ ∇p) 0EDFM 离散后变成对每个单元的质量守恒sum_j T_ij (p_j - p_i) 0其中T_ij是连接对传导率。这个方程对压力是线性的雅可比矩阵就是传导率矩阵本身。但如果是可压缩流体或考虑重力方程就非线性了。下面以可压缩单相流为例展示残差和雅可比的计算。function [R, J] residual_single_phase(p, T, V, phi, ct, dt, p_old) % p: 当前压力向量 % T: 传导率矩阵 % V: 单元体积向量 % phi: 孔隙度向量 % ct: 总压缩系数 % dt: 时间步长 % p_old: 上一时间步压力 % 累积项 accum V .* phi .* ct .* (p - p_old) / dt; % 流动项 flux T * p; % 注意 T 的对角线是负的流出和 % 残差 R accum flux; % 雅可比累积项对 p 的导数是 V.*phi.*ct/dt % 流动项对 p 的导数是 T J spdiags(V .* phi .* ct / dt, 0, length(p), length(p)) T; end这里T的构造要保证每行对角线是负的流出总和非对角线是正的传导率。flux T * p实际上算的是sum_j T_ij p_j但T_ii是负的sum_j T_ij所以T*p的第 i 个分量是sum_j T_ij (p_j - p_i)。雅可比J就是diag(accum_deriv) T。这段代码的解析导数很简单因为方程对压力是线性的可压缩性只影响累积项且是线性的。如果考虑非线性的 PVT 关系累积项对压力的导数就要用链式法则。3.2 两相流动的雅可比饱和度与压力的耦合导数两相流动油水或气水的方程更复杂未知量是压力和饱和度。以油水两相为例忽略毛管力方程是油相∇·(k kr_o/μ_o ∇p) q_o ∂(φ S_o)/∂t 水相∇·(k kr_w/μ_w ∇p) q_w ∂(φ S_w)/∂t离散后每个单元有两个残差油相质量守恒和水相质量守恒。未知量是p和S_w。雅可比矩阵是2N x 2N的块矩阵。下面这段代码计算一个单元对另一个单元的导数展示怎么用解析法写。function [R, J] residual_two_phase(p, Sw, T, V, phi, dt, p_old, Sw_old, fluid) % fluid: 含 kr_o, kr_w, mu_o, mu_w, Bo, Bw 等 N length(p); % 油相流度 lambda_o fluid.kr_o(Sw) ./ (fluid.mu_o .* fluid.Bo); % 水相流度 lambda_w fluid.kr_w(Sw) ./ (fluid.mu_w .* fluid.Bw); % 油相传导率 T_o T .* lambda_o; % 注意这里需要按连接对平均简化写法 T_w T .* lambda_w; % 累积项 accum_o V .* phi .* (1 - Sw) / fluid.Bo / dt - ... V .* phi .* (1 - Sw_old) / fluid.Bo / dt; accum_w V .* phi .* Sw / fluid.Bw / dt - ... V .* phi .* Sw_old / fluid.Bw / dt; % 残差 R [accum_o T_o * p; accum_w T_w * p]; % 雅可比块 % dR_o/dp T_o % dR_o/dSw d(lambda_o)/dSw * T * p d(accum_o)/dSw % dR_w/dp T_w % dR_w/dSw d(lambda_w)/dSw * T * p d(accum_w)/dSw dlambda_o fluid.dkr_o(Sw) ./ (fluid.mu_o .* fluid.Bo); dlambda_w fluid.dkr_w(Sw) ./ (fluid.mu_w .* fluid.Bw); dT_o_dSw T .* dlambda_o; dT_w_dSw T .* dlambda_w; daccum_o_dSw -V .* phi / fluid.Bo / dt; daccum_w_dSw V .* phi / fluid.Bw / dt; J11 T_o; J12 spdiags(dT_o_dSw * p daccum_o_dSw, 0, N, N); J21 T_w; J22 spdiags(dT_w_dSw * p daccum_w_dSw, 0, N, N); J [J11, J12; J21, J22]; end这里的关键是dT_o_dSw的计算。T_o T .* lambda_o而lambda_o是Sw的函数所以dT_o/dSw T .* dlambda_o/dSw。但注意T是矩阵lambda_o是向量T .* lambda_o在 MATLAB 里是每列乘对应的lambda_o这要求T的列对应单元编号。实际实现时传导率矩阵的每一行对应一个单元非零元对应连接对流度要用上游权或调和平均所以T_o的构造比这里写的复杂。上面的代码是简化版假设流度按单元平均且T已经是对称的。J12和J22的对角线元素是dT_o_dSw * p daccum_o_dSw但dT_o_dSw * p实际上是sum_j dT_ij/dSw * p_j这里用spdiags只取了对角线忽略了非对角线项。严格来说dR_o_i/dSw_j对j≠i也有非零值因为T_ij依赖于Sw_i和Sw_j的平均。如果流度用上游权T_ij只依赖于上游单元的Sw那么dR_o_i/dSw_j只在j是上游单元时非零。这部分是 EDFM 两相流动里最容易写错的地方建议先用数值差分验证雅可比。3.3 运移方程的离散对流扩散与时间项运移方程通常是组分浓度或示踪剂的输运∂(φ C)/∂t ∇·(v C) - ∇·(D ∇C) qEDFM 离散后对流项用上游权扩散项用连接对传导率。下面是对流扩散方程的残差和雅可比。function [R, J] residual_transport(C, p, T, V, phi, dt, C_old, D, fluid) % C: 浓度向量 % p: 压力向量用于算速度 % D: 扩散系数矩阵或连接对扩散传导率 N length(C); % 速度用达西定律v -k/mu * grad p离散后用传导率算流量 % 这里简化假设流量 q_ij T_ij * (p_j - p_i) % 对流项用上游权 q T * p; % 每个单元的净流量 % 上游权需要逐连接对判断这里用简化写法 % 实际实现时遍历连接对如果 q_ij 0上游是 j否则是 i % 下面用向量化近似假设流度已知对流项为 sum_j max(q_ij,0)*C_i - max(-q_ij,0)*C_j % 这里省略逐连接对循环直接给残差形式 accum V .* phi .* (C - C_old) / dt; conv zeros(N,1); diff zeros(N,1); % 需要连接对列表 I, J, T_ij % 假设已有 conn_i, conn_j, conn_T for k 1:length(conn_i) i conn_i(k); j conn_j(k); q_ij conn_T(k) * (p(j) - p(i)); if q_ij 0 conv(i) conv(i) q_ij * C(i); conv(j) conv(j) - q_ij * C(i); else conv(i) conv(i) q_ij * C(j); conv(j) conv(j) - q_ij * C(j); end % 扩散项 diff(i) diff(i) D(k) * (C(j) - C(i)); diff(j) diff(j) D(k) * (C(i) - C(j)); end R accum conv - diff; % 雅可比累积项对 C 的导数是 V.*phi/dt % 对流项对 C 的导数如果 q_ij 0d(conv_i)/dC_i q_ijd(conv_j)/dC_i -q_ij % 扩散项对 C 的导数是 -D 矩阵 % 这里省略具体组装建议用数值差分验证 end运移方程的雅可比比对流项更复杂因为上游权依赖于流动方向而流动方向又依赖于压力。如果压力是隐式求解的运移方程的雅可比还要包含对压力的导数形成全耦合。常见做法是顺序求解先解压力再解运移运移方程里压力视为已知。这样雅可比只对浓度求导简单很多。但顺序求解有时间步限制CFL 条件要求dt V/(q)否则对流项会震荡。4. 避坑与排查EDFM 求导代码里最容易翻车的五个地方4.1 传导率符号写反牛顿迭代不收敛现象残差在第一次迭代后反而变大或者压力出现负值。原因T矩阵的对角线应该是负的流出总和非对角线是正的传导率。如果写成对角线正、非对角线负T*p算出来的就是sum_j T_ij (p_i - p_j)符号反了。解决组装T后检查sum(T,2)是否接近零质量守恒或者用T*ones(N,1)看是否为零向量。4.2 流度平均方式与导数不一致现象雅可比矩阵用了解析导数但残差里流度用的是算术平均导数却按上游权推导致牛顿迭代收敛慢或震荡。原因残差和雅可比必须用同一套离散格式。如果残差里T_ij的流度是0.5*(lambda_i lambda_j)雅可比里dT_ij/dSw_i就是0.5*dlambda_i/dSw_i不能写成上游权的形式。解决先确定流度平均方式再推导数。建议用上游权因为导数简单且物理上更合理。4.3 裂缝单元体积为零导致累积项奇异现象裂缝单元的累积项V.*phi.*ct/dt为零雅可比对角线出现零元线性系统奇异。原因EDFM 里裂缝单元是降维的体积是L * aperture二维或A * aperture三维如果忘记乘开度体积就是零。解决裂缝单元的孔隙体积要显式乘开度且开度不能为零。如果开度极小累积项也很小这时候可以用稳态假设忽略裂缝内的累积但要在代码里显式处理。4.4 数值差分验证时步长选错现象用数值差分验证解析雅可比时误差很大甚至符号都不对。原因差分步长太大截断误差主导步长太小舍入误差主导。解决用中心差分步长取1e-6 * max(abs(x), 1)然后比较相对误差。如果相对误差在1e-4到1e-6之间说明解析导数基本正确。下面是一个验证雅可比的代码模板。function check_jacobian(fun, x, J_analytic) % fun: 返回残差的函数句柄 % x: 当前变量 % J_analytic: 解析雅可比 N length(x); J_num zeros(N, N); for i 1:N dx 1e-6 * max(abs(x(i)), 1); x_plus x; x_plus(i) x_plus(i) dx; x_minus x; x_minus(i) x_minus(i) - dx; R_plus fun(x_plus); R_minus fun(x_minus); J_num(:,i) (R_plus - R_minus) / (2*dx); end err norm(J_analytic - J_num, fro) / norm(J_num, fro); fprintf(相对误差: %e\n, err); if err 1e-4 warning(雅可比可能写错了); end end这个模板对每个变量做中心差分然后算 Frobenius 范数下的相对误差。如果误差大于1e-4就要检查解析导数。注意fun必须返回列向量且x的扰动不能导致物理上无意义的状态比如负饱和度。4.5 时间步长与 CFL 条件现象运移方程的解出现非物理震荡浓度超过 1 或小于 0。原因显式对流项的时间步长不满足 CFL 条件。解决要么用隐式对流要么限制dt min(V./q)。在 EDFM 里裂缝单元的V很小q可能很大所以裂缝处的 CFL 条件很苛刻。常见做法是对裂缝单元用隐式基质单元用显式或者全隐式但用迎风格式。下面是一个估算最大时间步长的代码。function dt_max calc_dt_cfl(T, p, V, phi) q T * p; % 净流量 % 每个单元的流出量 q_out max(q, 0); % CFL 条件 dt_max min(V .* phi ./ max(q_out, 1e-12)); endq_out是流出量V.*phi是孔隙体积比值就是流体充满一个单元的时间。取最小值就是全局时间步长。如果dt_max太小说明网格或裂缝开度需要调整。5. 进阶技巧用自动微分与稀疏结构加速雅可比组装5.1 为什么手动推导数容易错自动微分怎么帮上忙手动推雅可比在单相流里还行两相流加运移再加毛管力导数项会爆炸。MATLAB 没有内置的自动微分AD但可以用dlarray做深度学习的那种 AD或者用第三方工具箱。更实际的做法是把残差函数写成只依赖基本运算的形式然后用数值差分验证再用符号计算工具箱Symbolic Math Toolbox生成解析导数。下面是一个用符号计算生成雅可比的例子。syms p1 p2 Sw1 Sw2 real % 定义残差 T12 1e-3; % 传导率 lambda_o1 (1-Sw1)^2; % 简化油相流度 lambda_o2 (1-Sw2)^2; T_o12 T12 * 0.5*(lambda_o1 lambda_o2); R1 T_o12 * (p2 - p1); R2 T_o12 * (p1 - p2); % 雅可比 J jacobian([R1; R2], [p1; p2; Sw1; Sw2]); % 生成 MATLAB 函数 matlabFunction(J, File, J_two_phase.m);这段代码生成一个J_two_phase.m文件里面是解析雅可比。符号计算的好处是不会推错坏处是表达式可能很长生成的文件很大。实际用的时候可以把符号计算的结果复制到主代码里或者用matlabFunction的优化选项。5.2 稀疏结构复用与颜色分组EDFM 的雅可比是稀疏的但非零元的位置在时间步之间不变只要网格和连接对不变。所以可以预先计算稀疏结构然后用sparse的nzmax预分配。MATLAB 的sparse在组装时如果不断扩容会很慢建议先算非零元个数再用spalloc分配。下面是一个预分配的例子。% 假设连接对列表 conn_i, conn_j 已知 N max([conn_i; conn_j]); nnz_est length(conn_i) * 4; % 每个连接对贡献 4 个非零元两相 J spalloc(N, N, nnz_est); % 然后逐连接对填充 J for k 1:length(conn_i) i conn_i(k); j conn_j(k); % 填充 J(i,i), J(i,j), J(j,i), J(j,j) endspalloc预分配后填充速度比动态扩容快很多。如果变量是2N个压力和饱和度雅可比是2N x 2N非零元个数大约是4 * length(conn_i) 2N累积项对角线。颜色分组是进一步加速数值差分的方法把不相关的列分到同一组一次差分算多列。MATLAB 没有内置颜色分组但可以自己实现或者用numjac函数。5.3 用numjac做雅可比验证与替换MATLAB 的numjac是专门为稀疏雅可比设计的数值差分函数它支持颜色分组和稀疏结构。用法如下% fun: 残差函数返回列向量 % t: 当前时间 % y: 当前变量 % F: 当前残差 % thresh: 差分阈值 % JPattern: 稀疏结构矩阵 [J, nfe] numjac((t,y) fun(t,y), t, y, F, 1e-6, JPattern, 0);JPattern是稀疏结构可以用连接对列表构造。numjac会自动做颜色分组减少函数调用次数。如果解析雅可比写对了numjac的结果应该和解析雅可比很接近。如果差很多就用numjac的结果替换解析雅可比先保证收敛再慢慢查解析导数。这是血泪经验不要死磕解析导数先用数值雅可比跑通再优化。5.4 一个完整的验证流程最后给一个验证流程从简单到复杂先用单相流验证传导率矩阵再用两相流验证雅可比最后用运移验证对流扩散。每一步都用numjac或中心差分检查。如果某一步误差大就退回上一步。这个流程能帮你省下大量调试时间。希望帮到你。本文还有配套的精品资源点击获取