简介基于离散元方法DEM的岩石物理建模MATLAB代码包面向岩石力学、地质工程和矿产开采等领域的研究者与学生用于模拟颗粒材料在复杂受力条件下的力学行为并分析弹性模量、剪切模量、P波/S波速度及饱和度等因素对岩石物理特性的影响可应用于地震响应预测、地基处理评估和矿层开采设计等工程场景。包内共8个.m脚本涵盖弹性模量、速度-饱和度关系、纵横波速度关联等核心计算模块各脚本既可独立运行也可组合调用便于按需研究或二次开发。RAR压缩包仅5KB体量虽小但聚焦关键算法适合作为DEM在岩石物理应用中的快速入门参考。目前已有650人学习/下载作者自编代码并欢迎留言交流建模思路读者可借此快速掌握实现逻辑并针对实际课题进行扩展。1. 先分清两个DEM别让热词带偏方向打开搜索框直接敲“DEM”前几页大概率是ArcGIS、OPENTopography、12.5米DEM下载、DSM生成DEM这类测绘关键词但在这个缩写后面加上“岩石物理”画风会立刻转变。我经常看到有同行搜“岩石物理DEM”时被数字高程模型刷屏问岩石物理为什么和地形数据扯上关系。其实两者只是缩写撞车。在岩石物理学里DEM 全称是 Differential Effective Medium中文通常翻译成“微分等效介质”和数字高程模型没有任何关系。它解决的是从岩石微观组分出发预测宏观弹性性质的问题某块砂岩主要由石英组成孔隙度15%孔隙是球形还是扁裂缝对应的纵波速度可能差出不少。DEM 就是能把“孔隙形状”这个变量带进计算的理论。这篇内容适合三类人看一是做测井岩性解释、地震储层预测的地球物理工程师需要做岩石物理正演模板二是做岩石物理实验、想把实验室数据和理论模型对上的研究人员三是对有效介质理论感兴趣、想弄明白Differential Effective Medium到底怎么算的同学。我会把DEM的物理图像、关键公式、可运行代码和踩坑经验一次讲清楚。1.1 测绘圈的数字高程模型先花一点时间把容易混淆的概念处理掉。测绘、GIS领域常说的DEMDigital Elevation Model是数字高程模型也就是包含地面高程信息的栅格数据常用来做坡度分析、流域提取、地形可视化。热搜词里“12.5米DEM下载”“ArcGIS用DEM计算河流流域面积”“OPENTopography下载DEM教程”指的都是这个。这类DEM数据和岩石物理正演毫无关系。如果目标是做岩石物理建模搜索时建议用“DEM岩石物理”“微分等效介质Berryman”“Differential Effective Medium”来定位否则翻好几页看到的都是ArcMap和镶嵌工具的操作步骤浪费大量时间。1.2 岩石物理为什么需要“微分等效介质”岩石不是均匀材料。拿最常见的中粒砂岩来说骨架矿物主要是石英粒间有孔隙孔隙里可能填充空气、水或油气。如果只知道各组分的体积分数想预测整体弹性模量最简单的思路是做体积平均。可惜弹性模量对“几何”极其敏感同一个孔隙度球形孔和硬币状裂缝对Vp的影响天差地别这是体积平均完全无法刻画的短板。Voigt上限和Reuss下限只能给出一个很宽的范围实际岩石往往落在中间Hashin-Shtrikman界限比前两者紧一些但仍不包含孔隙形状信息。DEM模型则模拟了一个更接近真实的过程把矿物当作初始宿主相按微小的增量一步一步向宿主中加入孔隙夹杂物每一步都用当下的有效弹性性质重新计算夹杂物对整体的扰动。这个“边算边更新”的过程让孔隙形状得以进入方程也让结果在孔隙度达到中等水平时依然拥有较好的物理合理性。2. DEM的物理图像与数学方程2.1 从Eshelby夹杂到应变集中因子要理解DEM先要理解Eshelby在1957年解决的一个经典问题在无限大均匀弹性介质中放入一个椭球状夹杂物远场施加均匀应变夹杂物内部的应变场仍然是均匀的并且可以由一个四阶张量把远场应变与内部应变联系起来。这个张量通常记作A称为应变集中因子它依赖于夹杂物的形状纵横比和宿主介质的弹性模量。用生活化的方式理解把一块柔软果冻里压入一个硬球球周围的果冻受力时变形在球内是均匀的如果压入的是薄片即使体积相同薄片附近应力集中的程度也完全不一样。DEM模型反复利用这一思想把“一粒一粒加孔隙”的过程拆解成许多个Eshelby问题每次都用当前有效介质作为新的宿主。2.2 微分形式的等效介质方程DEM的常见微分方程写成dC / dφ (C_i - C) · A / (1 - φ)其中C 是当前有效介质的刚度矩阵C_i 是即将加入的夹杂物刚度矩阵A 是基于当前C和夹杂物形状计算出的应变集中因子φ 是当前孔隙度。分母中的1-φ是体积守恒修正因为我们是在已有孔隙的基础上继续加入新的孔隙需要把新加入的体积比例折算到剩余固体骨架中。这个微分增量过程从纯矿物的φ0出发一直积分到目标孔隙度。一个特别容易被忽略的要点是微分方程右侧的A必须用“当前”的有效模量计算不是用初始矿物模量。如果固定使用初始宿主去计算A那得到的其实是低孔隙度近似模型比如Kuster-Toksöz一阶近似DEM的独特之处就在于每一步都让夹杂物“看见”已经被前序孔隙软化过的介质这在高孔隙度时非常关键。2.3 DEM和自洽模型、VRH界限的区别有些教材会把DEM和SCASelf-Consistent Approximation自洽近似放在一起讨论。两者都基于Eshelby夹杂解但对夹杂物所处的环境假设不同。SCA假设夹杂物嵌入的是某种未知的、最终要求的有效介质然后通过自洽方程求解DEM则相当于一步一步把夹杂物加进不断更新的宿主介质中。三种方法的对比可以整理成下面这张表模型是否考虑孔隙形状处理夹杂物相互作用主要特点适用范围Voigt/Reuss界限否无上下限范围大只做边界检查自洽近似SCA是强相互作用近似对矿物相和孔隙平等的处理可能出现多解中高孔隙度混合物微分等效介质DEM是增量式处理逐步加入孔隙物理过程清晰低到中等孔隙度岩石实际项目中我喜欢把DEM和VRH配合使用先用DEM计算随孔隙度变化的骨架模量再对照VRH边界检查结果有没有落在合理区间内。如果计算出来的干岩石体积模量低于Reuss下限基本可以判断输入的孔隙纵横比或矿物模量设置有误。3. 手写一个DEM正演球孔模型教学版理论讲再多不如跑一段代码。下面我用Python实现一个严格但只针对球形孔隙的教学版DEM。这部分代码适合用于理解主流程工程中如果涉及裂缝或软孔需要更换Eshelby张量的计算我会在后面的问题章节说明。3.1 参数设置把石头拆成基质孔隙假设我们要模拟的岩石是石英砂岩骨架矿物为石英。理想化的矿物弹性参数如下参数取值单位说明石英体积模量K_m37.0GPa矿物基质背景模量石英剪切模量G_m44.0GPa矿物基质背景模量石英密度ρ_m2.65g/cm³骨架密度孔隙流体体积模量K_f2.25GPa水的常见取值流体密度ρ_f1.00g/cm³水密度目标孔隙度范围0 ~ 0.30无量纲常见砂岩有效范围夹杂物纵横比α1.0无量纲球形孔隙这里先固定成球孔是因为球形夹杂的Eshelby张量表达式最简单能让你最快看懂DEM主循环。实际储层孔隙通常不是球形的这组参数只作为教学起点。3.2 完整Python代码代码使用6x6 Voigt记号表示刚度矩阵用SciPy的solve_ivp对微分方程做自适应积分。注释里我尽量把每一步物理含义写清楚。import numpy as np from scipy.integrate import solve_ivp import matplotlib.pyplot as plt # 石英初始模量单位为GPa Km 37.0 Gm 44.0 def stiffness6(K, G): 由体积模量K和剪切模量G构造各向同性刚度矩阵(6x6 Voigt) C np.zeros((6, 6)) C11 K 4.0 * G / 3.0 C12 K - 2.0 * G / 3.0 C[0, 0] C[1, 1] C[2, 2] C11 C[0, 1] C[1, 0] C[0, 2] C[2, 0] C[1, 2] C[2, 1] C12 C[3, 3] C[4, 4] C[5, 5] G return C def kg_from_stiffness6(C): 从6x6刚度矩阵反推K和G用于检查各向同性 K (C[0, 0] 2.0 * C[0, 1]) / 3.0 G C[3, 3] return K, G def sphere_eshelby6(nu): 球形夹杂物的Eshelby张量(6x6 Voigt记号) a (7.0 - 5.0 * nu) / (15.0 * (1.0 - nu)) b (5.0 * nu - 1.0) / (15.0 * (1.0 - nu)) c (4.0 - 5.0 * nu) / (15.0 * (1.0 - nu)) S np.array([ [a, b, b, 0, 0, 0], [b, a, b, 0, 0, 0], [b, b, a, 0, 0, 0], [0, 0, 0, c, 0, 0], [0, 0, 0, 0, c, 0], [0, 0, 0, 0, 0, c], ]) return S def dem_rhs(phi, y): DEM微分方程右端项dC/dphi (Ci-C)A/(1-phi) C6 y.reshape(6, 6) K, G kg_from_stiffness6(C6) if K 0 or G 0: return np.zeros(36) nu (3.0 * K - 2.0 * G) / (2.0 * (3.0 * K G)) S6 sphere_eshelby6(nu) Ci np.zeros((6, 6)) # 孔隙刚度为零 I6 np.eye(6) A np.linalg.inv(I6 S6 np.linalg.inv(C6) (Ci - C6)) dC6 (Ci - C6) A / (1.0 - phi) return dC6.flatten() # 对孔隙度从0到0.3积分 phi_eval np.linspace(0.0, 0.30, 60) y0 stiffness6(Km, Gm).flatten() sol solve_ivp( dem_rhs, [0.0, 0.30], y0, t_evalphi_eval, methodLSODA, rtol1e-8, atol1e-8, ) # 提取干岩石模量 Kdry np.array([kg_from_stiffness6(sol.y[:, i].reshape(6, 6))[0] for i in range(len(phi_eval))]) Gdry np.array([kg_from_stiffness6(sol.y[:, i].reshape(6, 6))[1] for i in range(len(phi_eval))]) # 用Gassmann流体替换得到饱和岩石模量 Kfl 2.25 Ks Kdry (1.0 - Kdry / Km) ** 2 / ( phi_eval / Kfl (1.0 - phi_eval) / Km - Kdry / Km ** 2 ) Gs Gdry # Gassmann方程中剪切模量不变 # 计算纵横波速度 rho_m, rho_f 2.65, 1.00 rho (1.0 - phi_eval) * rho_m phi_eval * rho_f Vp_dry np.sqrt((Kdry 4.0 / 3.0 * Gdry) * 1e9 / (rho * 1e3)) Vs_dry np.sqrt(Gdry * 1e9 / (rho * 1e3)) Vp_sat np.sqrt((Ks 4.0 / 3.0 * Gs) * 1e9 / (rho * 1e3)) Vs_sat np.sqrt(Gs * 1e9 / (rho * 1e3)) plt.figure(figsize(10, 4)) plt.subplot(1, 2, 1) plt.plot(phi_eval, Kdry, -, labelK_dry) plt.plot(phi_eval, Ks, --, labelK_sat) plt.xlabel(porosity) plt.ylabel(bulk modulus (GPa)) plt.legend() plt.subplot(1, 2, 2) plt.plot(phi_eval, Vp_dry, -, labelVp_dry) plt.plot(phi_eval, Vp_sat, --, labelVp_sat) plt.plot(phi_eval, Vs_dry, :, labelVs) plt.xlabel(porosity) plt.ylabel(velocity (m/s)) plt.legend() plt.show()运行这段代码你会看到干岩石体积模量、剪切模量随着孔隙度增加单调下降饱和岩石体积模量高于干岩石并且Vp对孔隙度的敏感程度明显高于Vs。这正是砂岩储层岩石物理模板里最常见的趋势。3.3 读结果时需要注意的一点上面代码的积分起点是纯石英终点是孔隙度30%。在这个范围内球形孔隙DEM曲线通常平滑K和G不会出现非物理突变。但如果把孔隙度推到45%以上细孔的连通效应会变得显著DEM模型的假设就开始失效曲线可能出现过度刚硬的现象。实际项目中我通常把DEM正演的适用孔隙度控制在25%以内再高需要换其他模型或做实测数据标杆。另外球孔的DEM结果会明显高于相同孔隙度、更扁孔隙的计算结果。很多初学者第一次看到球孔曲线时会觉得“砂岩孔隙度20%体积模量怎么还剩下这么多”这其实是形状假设不同导致的。孔隙纵横比才是决定软硬的灵魂参数我建议你拿到岩心或铸体薄片资料后先用图像统计孔隙形状再去设定DEM中的α值。4. 常见问题与排错速查DEM模型在工程实现中容易踩的坑远不止代码本身。这里整理一些我在实际项目中反复遇到的问题。4.1 搜索结果端到端对不上先看孔隙纵横比最典型的场景是用DEM预测Vp发现结果和实测井曲线偏差很大。我第一步不是怀疑DEM方程而是检查输入的孔隙纵横比α。同样20%孔隙度α1的球孔和α0.05的扁孔预测的纵波速度可以差出接近20%。如果手头没有实测岩心数据至少可以做一个孔隙纵横比敏感性分析把α从0.01扫描到1落在实测范围内的α区间就是你的合理选择。4.2 干岩石Gassmann和“直接加流体”的DEM结果不一样这个问题最容易让人困惑。DEM理论上可以直接把水当作夹杂物加入矿物基质从而得到饱和岩石模量也可以先加入刚度为零的孔隙得到干岩石再用Gassmann方程做流体替换。两种做法结果并不完全一致因为DEM在加流体时Eshelby夹杂解和Gassmann的液体压强假设不同。工程实践上我更推荐后者先干岩石DEM再Gassmann流体替换。这样计算流程清晰也更方便对干岩石实测数据做校准。4.3 多矿物基质怎么设置初始宿主实际岩石很少是纯石英泥质砂岩、灰质砂岩都很常见。DEM起点可以不是单矿物而是一个多矿物混合基质的等效模量。通常做法是先用Voigt-Reuss-Hill平均算矿物混合物的等效模量把它作为DEM的初始宿主。要注意的是矿物混合体本身仍假设为各向同性、均匀的如果存在定向排列的黏土矿物还需要单独考虑各向异性问题。4.4 积分失败或负模量如果孔隙度过高或孔隙纵横比太小比如极薄的裂缝DEM微分方程可能变得刚硬显式欧拉步进容易产生负模量。我在代码里用solve_ivp的LSODA方法就是为了避免这个问题。如果你自己写循环做显式积分要特别注意步长最好设置动态步长或者把最大孔隙度限制在0.35以内。最后再分享一个多年实践的体会DEM不是“放之四海而皆准”的万能模板它最大的价值是提供了一个受控的物理框架让你能把“孔隙形状”“矿物模量”“流体类型”这几个地质参数定量地投射到弹性波速度上。我用它做储层正演模拟从来不追求一次算准而是拿着与测井数据做匹配反推出该层段可能的主导孔隙形状这比单纯算一条曲线有用得多。下次看到“DEM”搜索热词记得先停下来确认一下你要加载的是地形栅格还是准备向矿物里加孔隙的微分等效介质。本文还有配套的精品资源点击获取
岩石物理DEM正演：微分等效介质模型从原理到代码实现
简介基于离散元方法DEM的岩石物理建模MATLAB代码包面向岩石力学、地质工程和矿产开采等领域的研究者与学生用于模拟颗粒材料在复杂受力条件下的力学行为并分析弹性模量、剪切模量、P波/S波速度及饱和度等因素对岩石物理特性的影响可应用于地震响应预测、地基处理评估和矿层开采设计等工程场景。包内共8个.m脚本涵盖弹性模量、速度-饱和度关系、纵横波速度关联等核心计算模块各脚本既可独立运行也可组合调用便于按需研究或二次开发。RAR压缩包仅5KB体量虽小但聚焦关键算法适合作为DEM在岩石物理应用中的快速入门参考。目前已有650人学习/下载作者自编代码并欢迎留言交流建模思路读者可借此快速掌握实现逻辑并针对实际课题进行扩展。1. 先分清两个DEM别让热词带偏方向打开搜索框直接敲“DEM”前几页大概率是ArcGIS、OPENTopography、12.5米DEM下载、DSM生成DEM这类测绘关键词但在这个缩写后面加上“岩石物理”画风会立刻转变。我经常看到有同行搜“岩石物理DEM”时被数字高程模型刷屏问岩石物理为什么和地形数据扯上关系。其实两者只是缩写撞车。在岩石物理学里DEM 全称是 Differential Effective Medium中文通常翻译成“微分等效介质”和数字高程模型没有任何关系。它解决的是从岩石微观组分出发预测宏观弹性性质的问题某块砂岩主要由石英组成孔隙度15%孔隙是球形还是扁裂缝对应的纵波速度可能差出不少。DEM 就是能把“孔隙形状”这个变量带进计算的理论。这篇内容适合三类人看一是做测井岩性解释、地震储层预测的地球物理工程师需要做岩石物理正演模板二是做岩石物理实验、想把实验室数据和理论模型对上的研究人员三是对有效介质理论感兴趣、想弄明白Differential Effective Medium到底怎么算的同学。我会把DEM的物理图像、关键公式、可运行代码和踩坑经验一次讲清楚。1.1 测绘圈的数字高程模型先花一点时间把容易混淆的概念处理掉。测绘、GIS领域常说的DEMDigital Elevation Model是数字高程模型也就是包含地面高程信息的栅格数据常用来做坡度分析、流域提取、地形可视化。热搜词里“12.5米DEM下载”“ArcGIS用DEM计算河流流域面积”“OPENTopography下载DEM教程”指的都是这个。这类DEM数据和岩石物理正演毫无关系。如果目标是做岩石物理建模搜索时建议用“DEM岩石物理”“微分等效介质Berryman”“Differential Effective Medium”来定位否则翻好几页看到的都是ArcMap和镶嵌工具的操作步骤浪费大量时间。1.2 岩石物理为什么需要“微分等效介质”岩石不是均匀材料。拿最常见的中粒砂岩来说骨架矿物主要是石英粒间有孔隙孔隙里可能填充空气、水或油气。如果只知道各组分的体积分数想预测整体弹性模量最简单的思路是做体积平均。可惜弹性模量对“几何”极其敏感同一个孔隙度球形孔和硬币状裂缝对Vp的影响天差地别这是体积平均完全无法刻画的短板。Voigt上限和Reuss下限只能给出一个很宽的范围实际岩石往往落在中间Hashin-Shtrikman界限比前两者紧一些但仍不包含孔隙形状信息。DEM模型则模拟了一个更接近真实的过程把矿物当作初始宿主相按微小的增量一步一步向宿主中加入孔隙夹杂物每一步都用当下的有效弹性性质重新计算夹杂物对整体的扰动。这个“边算边更新”的过程让孔隙形状得以进入方程也让结果在孔隙度达到中等水平时依然拥有较好的物理合理性。2. DEM的物理图像与数学方程2.1 从Eshelby夹杂到应变集中因子要理解DEM先要理解Eshelby在1957年解决的一个经典问题在无限大均匀弹性介质中放入一个椭球状夹杂物远场施加均匀应变夹杂物内部的应变场仍然是均匀的并且可以由一个四阶张量把远场应变与内部应变联系起来。这个张量通常记作A称为应变集中因子它依赖于夹杂物的形状纵横比和宿主介质的弹性模量。用生活化的方式理解把一块柔软果冻里压入一个硬球球周围的果冻受力时变形在球内是均匀的如果压入的是薄片即使体积相同薄片附近应力集中的程度也完全不一样。DEM模型反复利用这一思想把“一粒一粒加孔隙”的过程拆解成许多个Eshelby问题每次都用当前有效介质作为新的宿主。2.2 微分形式的等效介质方程DEM的常见微分方程写成dC / dφ (C_i - C) · A / (1 - φ)其中C 是当前有效介质的刚度矩阵C_i 是即将加入的夹杂物刚度矩阵A 是基于当前C和夹杂物形状计算出的应变集中因子φ 是当前孔隙度。分母中的1-φ是体积守恒修正因为我们是在已有孔隙的基础上继续加入新的孔隙需要把新加入的体积比例折算到剩余固体骨架中。这个微分增量过程从纯矿物的φ0出发一直积分到目标孔隙度。一个特别容易被忽略的要点是微分方程右侧的A必须用“当前”的有效模量计算不是用初始矿物模量。如果固定使用初始宿主去计算A那得到的其实是低孔隙度近似模型比如Kuster-Toksöz一阶近似DEM的独特之处就在于每一步都让夹杂物“看见”已经被前序孔隙软化过的介质这在高孔隙度时非常关键。2.3 DEM和自洽模型、VRH界限的区别有些教材会把DEM和SCASelf-Consistent Approximation自洽近似放在一起讨论。两者都基于Eshelby夹杂解但对夹杂物所处的环境假设不同。SCA假设夹杂物嵌入的是某种未知的、最终要求的有效介质然后通过自洽方程求解DEM则相当于一步一步把夹杂物加进不断更新的宿主介质中。三种方法的对比可以整理成下面这张表模型是否考虑孔隙形状处理夹杂物相互作用主要特点适用范围Voigt/Reuss界限否无上下限范围大只做边界检查自洽近似SCA是强相互作用近似对矿物相和孔隙平等的处理可能出现多解中高孔隙度混合物微分等效介质DEM是增量式处理逐步加入孔隙物理过程清晰低到中等孔隙度岩石实际项目中我喜欢把DEM和VRH配合使用先用DEM计算随孔隙度变化的骨架模量再对照VRH边界检查结果有没有落在合理区间内。如果计算出来的干岩石体积模量低于Reuss下限基本可以判断输入的孔隙纵横比或矿物模量设置有误。3. 手写一个DEM正演球孔模型教学版理论讲再多不如跑一段代码。下面我用Python实现一个严格但只针对球形孔隙的教学版DEM。这部分代码适合用于理解主流程工程中如果涉及裂缝或软孔需要更换Eshelby张量的计算我会在后面的问题章节说明。3.1 参数设置把石头拆成基质孔隙假设我们要模拟的岩石是石英砂岩骨架矿物为石英。理想化的矿物弹性参数如下参数取值单位说明石英体积模量K_m37.0GPa矿物基质背景模量石英剪切模量G_m44.0GPa矿物基质背景模量石英密度ρ_m2.65g/cm³骨架密度孔隙流体体积模量K_f2.25GPa水的常见取值流体密度ρ_f1.00g/cm³水密度目标孔隙度范围0 ~ 0.30无量纲常见砂岩有效范围夹杂物纵横比α1.0无量纲球形孔隙这里先固定成球孔是因为球形夹杂的Eshelby张量表达式最简单能让你最快看懂DEM主循环。实际储层孔隙通常不是球形的这组参数只作为教学起点。3.2 完整Python代码代码使用6x6 Voigt记号表示刚度矩阵用SciPy的solve_ivp对微分方程做自适应积分。注释里我尽量把每一步物理含义写清楚。import numpy as np from scipy.integrate import solve_ivp import matplotlib.pyplot as plt # 石英初始模量单位为GPa Km 37.0 Gm 44.0 def stiffness6(K, G): 由体积模量K和剪切模量G构造各向同性刚度矩阵(6x6 Voigt) C np.zeros((6, 6)) C11 K 4.0 * G / 3.0 C12 K - 2.0 * G / 3.0 C[0, 0] C[1, 1] C[2, 2] C11 C[0, 1] C[1, 0] C[0, 2] C[2, 0] C[1, 2] C[2, 1] C12 C[3, 3] C[4, 4] C[5, 5] G return C def kg_from_stiffness6(C): 从6x6刚度矩阵反推K和G用于检查各向同性 K (C[0, 0] 2.0 * C[0, 1]) / 3.0 G C[3, 3] return K, G def sphere_eshelby6(nu): 球形夹杂物的Eshelby张量(6x6 Voigt记号) a (7.0 - 5.0 * nu) / (15.0 * (1.0 - nu)) b (5.0 * nu - 1.0) / (15.0 * (1.0 - nu)) c (4.0 - 5.0 * nu) / (15.0 * (1.0 - nu)) S np.array([ [a, b, b, 0, 0, 0], [b, a, b, 0, 0, 0], [b, b, a, 0, 0, 0], [0, 0, 0, c, 0, 0], [0, 0, 0, 0, c, 0], [0, 0, 0, 0, 0, c], ]) return S def dem_rhs(phi, y): DEM微分方程右端项dC/dphi (Ci-C)A/(1-phi) C6 y.reshape(6, 6) K, G kg_from_stiffness6(C6) if K 0 or G 0: return np.zeros(36) nu (3.0 * K - 2.0 * G) / (2.0 * (3.0 * K G)) S6 sphere_eshelby6(nu) Ci np.zeros((6, 6)) # 孔隙刚度为零 I6 np.eye(6) A np.linalg.inv(I6 S6 np.linalg.inv(C6) (Ci - C6)) dC6 (Ci - C6) A / (1.0 - phi) return dC6.flatten() # 对孔隙度从0到0.3积分 phi_eval np.linspace(0.0, 0.30, 60) y0 stiffness6(Km, Gm).flatten() sol solve_ivp( dem_rhs, [0.0, 0.30], y0, t_evalphi_eval, methodLSODA, rtol1e-8, atol1e-8, ) # 提取干岩石模量 Kdry np.array([kg_from_stiffness6(sol.y[:, i].reshape(6, 6))[0] for i in range(len(phi_eval))]) Gdry np.array([kg_from_stiffness6(sol.y[:, i].reshape(6, 6))[1] for i in range(len(phi_eval))]) # 用Gassmann流体替换得到饱和岩石模量 Kfl 2.25 Ks Kdry (1.0 - Kdry / Km) ** 2 / ( phi_eval / Kfl (1.0 - phi_eval) / Km - Kdry / Km ** 2 ) Gs Gdry # Gassmann方程中剪切模量不变 # 计算纵横波速度 rho_m, rho_f 2.65, 1.00 rho (1.0 - phi_eval) * rho_m phi_eval * rho_f Vp_dry np.sqrt((Kdry 4.0 / 3.0 * Gdry) * 1e9 / (rho * 1e3)) Vs_dry np.sqrt(Gdry * 1e9 / (rho * 1e3)) Vp_sat np.sqrt((Ks 4.0 / 3.0 * Gs) * 1e9 / (rho * 1e3)) Vs_sat np.sqrt(Gs * 1e9 / (rho * 1e3)) plt.figure(figsize(10, 4)) plt.subplot(1, 2, 1) plt.plot(phi_eval, Kdry, -, labelK_dry) plt.plot(phi_eval, Ks, --, labelK_sat) plt.xlabel(porosity) plt.ylabel(bulk modulus (GPa)) plt.legend() plt.subplot(1, 2, 2) plt.plot(phi_eval, Vp_dry, -, labelVp_dry) plt.plot(phi_eval, Vp_sat, --, labelVp_sat) plt.plot(phi_eval, Vs_dry, :, labelVs) plt.xlabel(porosity) plt.ylabel(velocity (m/s)) plt.legend() plt.show()运行这段代码你会看到干岩石体积模量、剪切模量随着孔隙度增加单调下降饱和岩石体积模量高于干岩石并且Vp对孔隙度的敏感程度明显高于Vs。这正是砂岩储层岩石物理模板里最常见的趋势。3.3 读结果时需要注意的一点上面代码的积分起点是纯石英终点是孔隙度30%。在这个范围内球形孔隙DEM曲线通常平滑K和G不会出现非物理突变。但如果把孔隙度推到45%以上细孔的连通效应会变得显著DEM模型的假设就开始失效曲线可能出现过度刚硬的现象。实际项目中我通常把DEM正演的适用孔隙度控制在25%以内再高需要换其他模型或做实测数据标杆。另外球孔的DEM结果会明显高于相同孔隙度、更扁孔隙的计算结果。很多初学者第一次看到球孔曲线时会觉得“砂岩孔隙度20%体积模量怎么还剩下这么多”这其实是形状假设不同导致的。孔隙纵横比才是决定软硬的灵魂参数我建议你拿到岩心或铸体薄片资料后先用图像统计孔隙形状再去设定DEM中的α值。4. 常见问题与排错速查DEM模型在工程实现中容易踩的坑远不止代码本身。这里整理一些我在实际项目中反复遇到的问题。4.1 搜索结果端到端对不上先看孔隙纵横比最典型的场景是用DEM预测Vp发现结果和实测井曲线偏差很大。我第一步不是怀疑DEM方程而是检查输入的孔隙纵横比α。同样20%孔隙度α1的球孔和α0.05的扁孔预测的纵波速度可以差出接近20%。如果手头没有实测岩心数据至少可以做一个孔隙纵横比敏感性分析把α从0.01扫描到1落在实测范围内的α区间就是你的合理选择。4.2 干岩石Gassmann和“直接加流体”的DEM结果不一样这个问题最容易让人困惑。DEM理论上可以直接把水当作夹杂物加入矿物基质从而得到饱和岩石模量也可以先加入刚度为零的孔隙得到干岩石再用Gassmann方程做流体替换。两种做法结果并不完全一致因为DEM在加流体时Eshelby夹杂解和Gassmann的液体压强假设不同。工程实践上我更推荐后者先干岩石DEM再Gassmann流体替换。这样计算流程清晰也更方便对干岩石实测数据做校准。4.3 多矿物基质怎么设置初始宿主实际岩石很少是纯石英泥质砂岩、灰质砂岩都很常见。DEM起点可以不是单矿物而是一个多矿物混合基质的等效模量。通常做法是先用Voigt-Reuss-Hill平均算矿物混合物的等效模量把它作为DEM的初始宿主。要注意的是矿物混合体本身仍假设为各向同性、均匀的如果存在定向排列的黏土矿物还需要单独考虑各向异性问题。4.4 积分失败或负模量如果孔隙度过高或孔隙纵横比太小比如极薄的裂缝DEM微分方程可能变得刚硬显式欧拉步进容易产生负模量。我在代码里用solve_ivp的LSODA方法就是为了避免这个问题。如果你自己写循环做显式积分要特别注意步长最好设置动态步长或者把最大孔隙度限制在0.35以内。最后再分享一个多年实践的体会DEM不是“放之四海而皆准”的万能模板它最大的价值是提供了一个受控的物理框架让你能把“孔隙形状”“矿物模量”“流体类型”这几个地质参数定量地投射到弹性波速度上。我用它做储层正演模拟从来不追求一次算准而是拿着与测井数据做匹配反推出该层段可能的主导孔隙形状这比单纯算一条曲线有用得多。下次看到“DEM”搜索热词记得先停下来确认一下你要加载的是地形栅格还是准备向矿物里加孔隙的微分等效介质。本文还有配套的精品资源点击获取
