(19)国家知识产权局 (12)发明 专利 (10)授权公告 号 (45)授权公告日 (21)申请 号 202111303690.4 (22)申请日 2021.11.05 (65)同一申请的已公布的文献号 申请公布号 CN 114021498 A (43)申请公布日 2022.02.08 (73)专利权人 中国矿业大 学 地址 221000 江苏省徐州市铜山区大 学路1 号 (72)发明人 孙超　岳建华　王尚旭　姜志海　 唐跟阳　 (74)专利代理 机构 徐州千秋知识产权代理事务 所(普通合伙) 32556 专利代理师 李翩 (51)Int.Cl. G06F 30/28(2020.01)G06F 30/23(2020.01) G06F 111/10(2020.01) G06F 113/08(2020.01) G06F 119/14(2020.01) (56)对比文件 CN 110263434 A,2019.09.20 CN 109471 168 A,2019.0 3.15 CN 109116420 A,2019.01.01 CN 105607123 A,2016.0 5.25 US 20170 59742 A1,2017.0 3.02 审查员 徐晓 (54)发明名称 一种预测多相孔隙介质弹性模量的高效数 值模拟方法 (57)摘要 本发明公开了一种预测多相孔隙介质弹性 模量的高效数值模拟方法， 包括： 通过流体分布 图像获得流体分布特征； 定义流体密度、 粘度、 体 积模量分布特征； 构建频率域Biot方程的位移 ‑ 压力形式； 约束纵波边界条件； 全局粗化网格剖 分， 局部网格精细剖分求解； 使用有限元求解方 程， 获得应变固体位移场和流体压力场； 计算纵 波模量。 本发 明提供的预测多相孔隙介质弹性模 量的高效数值模拟方法， 相比于以前求解孔弹性 Biot流‑固耦合位移方程， 实现了对部分饱和岩 石样品频散和衰减预测的有限元数值模拟方法， 特别适用于模 型中存在小非均质目标， 可避免常 规Biot固体位移 ‑流体位移方程计算效率低， 四 面体网格剖分内存占用量大， 计算时间长的问 题。 权利要求书2页 说明书7页 附图6页 CN 114021498 B 2022.10.11 CN 114021498 B 1.一种预测多相孔隙介质弹性模量的高效数值模拟方法， 其特 征在于， 包括如下步骤： S1： 通过流体分布图像获得第 相流体分布范围 ； 为混合流体的最大 相数； S2： 定义流体密度 为： ， x,y,z为笛卡尔坐标系 空间坐标； 是第 相流体密度； S3： 定义流体粘度为 ： ； 是第 相流体粘度； S4： 定义流体体积模量 为： ； S5： 构建频率 域Biot方程的位移 ‑压力形式： 其中， ； ； 为饱和岩石密度； ； 为岩石骨架密度， 为孔隙度， 为流体复密度， 为扭曲因子， 为岩 石渗透率； 为总应力， C为干燥岩石骨架柔度， 为应变张量， 为Biot‑Willis系数， 为 流体压力， I 为单位张量矩阵， F为周期震荡的外力； 为角频率， 为固体位移张量； 其中， ， ， 为岩石骨架体积模量； S6： 约束边界条件： 数值模拟的目标是实验室岩石样品， 岩石样品轴向方向是z方向， 垂 直于z方向的平面上存在着正交的x、 y方向 ， 根据纵波模量的定义 ， 横向应变 ， 及流体边界没有流体流出， 因此边界条件为： 下边界为固定边界， 需要同时满足两个条件： 边界法向上流体加速度为0， 固体位移是 0， 因此边界条件是： 权　利　要　求　书 1/2 页 2 CN 114021498 B 2上边界加载周期震荡应力， 对应的边界条件是边界法向应力是 ， 流体流动的加速度 为0， 因此边界条件是： S7： 全局粗 化网格剖分， 局部网格精细剖分求 解； S8： 使用有限元求 解方程， 获得应 变固体位移场和流体压力场； S9： 计算纵波模量M： 其中， 为垂向应力， 为垂向应 变； 则纵波速度 和纵波衰减 为： 其中， 和 是虚步和实步 算子。 2.根据权利要求1所述的预测多相孔隙介质弹性模量的高效数值模拟方法， 其特征在 于， 步骤S1中流体分布图像通过CT扫描获得。 3.根据权利要求1所述的预测多相孔隙介质弹性模量的高效数值模拟方法， 其特征在 于， 步骤S7 具体包括： S71： 利用有限元算法使用粗 化网格求 解方程； S72： 计算所有网格元 素上的偏微分方程的残差； S73： 估计所有网格元素的解中误差， 其中第l 方程的L2范数误差 为 ， h为局 部网格尺寸， 为第l个偏微分方程的残差， 为偏微分方程稳定性估计导数阶数， 默认值 是2， 由于参数C是任意常数， 计算出的误差实际上 是误差指示因子； S74： 如果已进行了预设的优化次数或超过了最大有限元剖分数， 则终止执 行； S75： 建立局部误差指示因子， ， 其中A是网格的面积， 根据局部误差指标 的大小细化有限元； S76： 重复步骤S71 ‑S75。权　利　要　求　书 2/2 页 3 CN 114021498 B 3