弹塑性变形油藏中多相渗流的数值模拟

弹塑性变形油藏中多相渗流的数值模拟 第 16 卷第 1 期 V o l. 16 N o. 1 计 算 力 学 学 报 1999 年 2 月 CH IN ESE JOU RNAL O F COM PU TA T IONAL M ECHAN ICS February 1999 X 弹塑性变形油藏中多相渗流的数值模拟 冉启全 顾小芸 中国石油天然气集团公司石油勘探开发科学研究院, 北京, 100083 中国科学院力学研究所, 北京, 100080 摘 要 基于流固耦合力学理论, 建立了弹塑性变形油藏中多相渗流的数学模型, 该模 型考虑了渗流与变形的耦合作用, 以及注采交变载荷作用下油藏多孔介质的弹塑性变 形特征; 给出了耦合数值模拟方法和算例。 关键词 弹塑性变形; 多相渗流; 耦合作用; 数值模拟; 油藏 分类号 319; 35713 T E O 1 前 言 计算渗流力学在石油工程中得到了相当重视和广泛应用。油藏渗流的计算机数值计算 即 油藏数值模拟已成为油藏工程师进行油藏工程 分析 定性数据统计分析pdf销售业绩分析模板建筑结构震害分析销售进度分析表京东商城竞争战略分析 必不可少的有力工具。但 长期以来, 在研 究油藏开采过程中的渗流问题时, 均采用纯渗流力学理论进行研究, 不考虑油藏渗流与岩土变 形的耦合作用及其影响。近年来, 随着流固耦合力学理论和数值计算技术的发展, 人们对油藏 开采过程中的力学机制有了更深入的了解和认识: 岩土在多相流体渗流作用下会发生变形, 而 岩土变形反过来又影响多相流体的渗流。因此, 油藏开采过程中的多相渗流问题是流固耦合问 题, 而应用计算机对油藏开采过程中的流固耦合渗流问题进行数值计算和模拟, 称为流固耦合 油藏数值模拟。 [1~ 3 ] 国外在流固耦合油藏数值模拟方面已作过研究工作, 不少学者作了有益的探索 。但这 [4, 5 ] 些研究均将岩土介质视为弹性介质, 对岩土受力变形过程作简化处理。大量的实验研究 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 明, 高温高压下的油藏岩土, 特别是在注采交变载荷作用下其受力变形特性极为复杂, 不仅有 可恢复的弹性变形, 更多的表现为不可恢复的塑性变形。为此, 本文将渗流力学与岩土力学相 结合, 充分考虑注采交变载荷作用下油藏岩土的弹塑性变形特性, 建立了体现多相渗流与弹塑 性变形的流固耦合模型, 用有限差分与有限元相结合的耦合数值计算方法对模型进行求解。 2 基本理论 流固耦合问题由流固耦合方程来描述。这组方程包括流固耦合渗流方程和流固耦合变形 方程, 而未知量含有描述渗流现象的变量和描述变形现象的变量, 流固耦合效应在两类方程中 均要得到体现。 211 渗流方程 油藏可以看成是势量所形成的场。当油藏内各处势量存在差异时就会发生渗流。首先引 入流动势的概念, 定义下式: 5 P - CD 1 A A A X 收稿日期: 1997206203 冉启全: 男, 1965 年生, 博士后, 高级工程师1 期 冉启全等: 弹塑性变形油藏中多相渗流的数值模拟 25 式中 A o和 X分别代表油相和水相, 5 A为流体相 A的流动势, P A为流体相 A的压力, C A为流体 相 A的重度, D 为由某一基准面算起的深度。 由于渗流发生在弹塑性变形多孔介质中, 因而不但流体质点具有一定的渗 流速度, 岩土质 点也有一定的运动速度, 所以流体速度可以描述为 ? ? ? v v + v 2 A rA s ? ? ? 式中 v A为流体相 A的绝对速度; v s 为岩土质点的绝对速度; v rA为流 体相 A的相对速度, 其表达 式为 ? 1 v V 3 rA A U S A 式中 U为孔隙度; S A为流体相 A的饱和度; V A为流体相 A的达西速度, 其表达式为 K K rA V A - grad5 A 4 LA 式中 K 为绝对渗透率; K 为流体相 A的相对渗透率; L 为流体相 A的粘 度。 rA A 根据质量守恒定律, 并考虑源汇项, 流固耦合渗流的微分方程为 ? K K roQ o 5 UQ oS o div grad5 - div U S Qv + q 5 o o o s o L o 5t ? K K rw Q w 5 U Q w S w div grad5 w - div U S w Q w v s + qw 6 L w 5t 式中 q 为单位时间单位体积内的产出或注入流体质量。 A 由于油藏为油、水两相渗流, 因此还应有下列两方程成立, 即: 毛管压力方程 P cow P o - P w 7 饱和度方程 S o + S w 1 8 式中 P cow 为油水间的毛管压力。 由于流固耦合效应, 渗流微分方程的每一项均含有体现流固耦合效应的未 知参数 U、 K 和 ? v s , 它们受岩土变形或运动的影响而不断变化, 因此流固耦合渗流方程具 有高度的非线性特 性, 只有联立流固耦合变形方程才能进行求解。根据流固耦合变形方程的求解结果, 就可以用 ? 下面各式计算 U、 K 和 v s 值: ? 3 ? U o + E v 1 + E v U 0 5 U U , K K , v 9 0 s 1 + E v 1 + E v 5t ? 式中 E v 为体积应变, U 为岩土质点的位移矢量, t 为时间。 212 变形方程 油藏岩土在注采交变载荷及流固耦合作用下, 发生显著的弹塑性变形, 其本构方程需要根 据塑性增量理论确定。塑性增量理论认为: 在弹塑性变形过程中, 任何阶段应变增量 dE ij 可以 e p 线性地分为弹性应变增量 dE ij 与塑性应变增量 dE ij, 即: e p dE dE + dE 10 ij ij ij 弹性应变增量与应力增量 dRij 之间的关系可通过广义虎克定律来表示, 即: e e dRij D ij lddE k l 11 e 式中D ijk l 为弹性刚度矩阵。 塑性应变增量要根据塑性增量理论计算。屈服准则、流动法则及硬化规律是塑性增量理论计 算 力 学 学 报 16 卷 26 的三个组成部分。 首先要确定油藏岩土某点应力达到弹性极限后出现塑性变形的屈服条件, 这里选用 D rucker2P rager 屈服条件, 即: F AI + J - k 0 12 1 2 式中 F 为屈服函数; I 1、J 2 分别为应力张量第一不变量和应力偏量第二不变量; A、 k 为屈服函 数参数, 其表达式为 2sinU 6cco sU A , k 13 3 3 - sinU 3 3 - sinU 式中 c 和 U分别为内聚力和摩擦角。 其次选取岩土所服从的流动法则, 塑性应变增量的方向就由流动法则确定, 其表达式为: 5 Q P dE ij dK 14 5Rij 式中 dK为塑性因子, Q 为塑性势函数。这样就确定了塑性应变增量的方向, 也就确定了塑性应 变增量各分量的比值。 最后还要确定岩土的硬化规律, 这样塑性应变增量的大小就可根据硬化规 律来计算。 根据上述塑性增量理论可以得出塑性应变增量的表达式: 5F e D k lm ndE m n 5R k l 5 Q p dE ij 15 5F 5 Q 5Rcd e H + D A abcd 5R 5R ab cd 并可得到应力增量与总应变增量之间关系的表达式: 5 Q 5F e e D ij rsD m nk l 5Rrs 5R m n e dR D - dE 16 ij ijk l k l 5F 5 Q e H A + D abcd 5Rab 5Rcd 该式给出了一般形式的弹塑性本构方程。其弹塑性刚度矩阵为 5 Q 5F e e D ij rsD m nk l 5R 5R rs m n ep e D ijk l D ijk l - 17 5F 5 Q e H + D A abcd 5Rab 5Rcd 确定了岩土的弹塑性刚度矩阵, 也就确定了岩土的弹塑性本构关系。 设岩土骨架所发生的变形是小变形, 则几何方程为 1 E ij u i, j + u j, i 18 2 式中 E 为应变分量, u 为位移分量。 ij i, j T 由于单元体处于力的平衡状态, 因此按总应力 R 表示的平衡微分方程为 ij T R + [ 1 - U Q + U S Q + U S Q ]f 0 19 ij, j s o o w w i T 式中 Q s 为岩土骨架密度, f i [ 0 0 g ] , g 为重力加速度。 有效应力 Rij 可由太沙基公式给出: T Rij Rij - P Dij 20 式中 P 为等效孔隙压力, Dij 为 K ronecker 常数。因此以有效应力表示的平衡微分方程为1 期 冉启全等: 弹塑性变形油藏中多相渗流的数值模拟 27 Rij, j + P Dij , j + [ 1 - U Q s + U S oQ o + U S w Q w ]f i 0 21 由于流固耦合效应, 平衡微分方程中含有孔隙流体压力和饱和度项, 它们体现了流固耦合 效应对岩土受力变形的影响, 只有联立渗流方程才能求解。 3 数值模拟 渗流方程和变形方程均含有体现流固耦合效应的参数或变量, 因而不能单独进行求解, 必 须将两者联合起来才能求解。因此, 本文采用有限差分与有限元耦合方法进行交替求解, 即渗 流方程采用有限差分法求解, 变形方程采用有限元方法求解, 并在两者间交替进行。 311 渗流方程的求解 首先对渗流方程进行差分离散, 形成差分方程。对渗流方程进行空间和时间项的差分离 散, 并在等式两边同乘以 $x i$y j$z k, 同时令: $y $z K K Q $y $z K K Q j k rA A j k rA A 1 1 T x A i+ , T x A i- 2 1 1 2 1 1 $x i+ u A i+ $x i- u A i- 2 2 2 2 $x $z K K Q $x $z K K Q i k rA A i k rA A 1 1 T , T y A j+ y Aj- 22 2 1 1 2 1 1 $y u $y u j+ A j+ j- A j- 2 2 2 2 $x i$y j K K rAQ A $x i$y j K K rAQ A 1 1 T , T z A k+ z A k- 2 1 2 1 1 1 $z k+ u A $z k- u A k+ k- 2 2 2 2 V $x $y $z , Q q V 23 ijk i j k A A ijk 5Q V C V UQ 5U A ijk tA ijk A C tA Q AS A + U S A , C 1A , C 2A 24 5P 5P $ t $ t A ? 流固耦合项 div U S AQ Av s 与时间无关, 所以对其进行空间离散后, 可 以表示为: U S AQ A i+ 1 - U S AQ A i- 1 U S AQ A j+ 1 - U S AQ A j- 1 C sA v sx + v sy + 1 1 1 1 $x i- + $x i+ $y j- + $y j+ 2 2 2 2 U S AQ A k+ 1 - U S AQ A k- 1 $E v v sz + U S AQ A 25 1 1 $z k- + $z k+ $ t 2 2 为简化方程, 引入线性差分算子: 1 1 $ x T x A$ x P A 3 x Ai+ P A i+ 1 - P A i + 3 x A i- P A i- 1 - P A i 2 2 1 1 $ y T y A$ y P A 3 y A j+ P A j+ 1 - P A j + 3 y A j- P A j- 1 - P A j 2 2 1 1 $ zT z A$ zP A 3 z A k+ P A k+ 1 - P A k + 3 z A k- P A k- 1 - P A k 2 2 26 1 1 $ x T x AC A$ xD 3 x AC A i+ D i+ 1 - D i + T x AC A i+ D i- 1 - D i 2 2 1 1 $ y T y AC A$ yD 3 y AC A j+ D j+ 1 - D j + T y AC A j+ D j- 1 - D j 2 2 1 1 $ T C$ D 3 C D - D + T C D - D z z A A z z A A k+ k+ 1 k z A A k+ k- 1 k 2 2 通过化简、整理, 最后得到的差分方程为: n n+ 1 n n+ 1 n n+ 1 $ x T x o$ x P o + $ y T y o$ y P o + $ zT z o$ zP o - n n n n n n n n $ T C$ D - $ T C$ D - $ T C$ D - C + Q x x o o x y y o o y z z o o z so o n n+ 1 n n n+ 1 n C 1o P oM - P oM + C 2o S o M - S oM 27 n n+ 1 n n+ 1 n n+ 1 $ x T xw $ x P w + $ y T yw $ y P w + $ zT zw $ zP w -计 算 力 学 学 报 16 卷 28 n n n n n n n n $ T C $ D - $ T C $ D - $ T C $ D - C + Q x xw w x y yw w y z zw w z sw w n n+ 1 n n n+ 1 n C 1w P wM - P wM + C 2w S wM - S wM 28 在形成差分方程后, 就可以对差分方程进行求解。本文采用的求解方法为隐压显饱法, 即 隐式求解压力方程。显示求解饱和度方程。其基本思想为: 1 根据饱和度方程 8 , 通过乘以 适当的系数, 合并油方程和水方程, 以消去差分方程中的变量 S o和S w , 从而得到一个只含变量 P 和 P 的压力方程。2 由毛管压力方程 7 将 P 表示为 P 和 P 的形式, 得到一个只含变 o w w o cow n+ 1 n+ 1 量 P 的压力差分方程。3 隐式求解压力差分方程得到 P , 并将 P 代入毛管压力方程 o o o n+ 1 n+ 1 n+ 1 。4 将 P 代入水相差分方程 28 , 显式求解出水相饱和度 S , 再代入 7 , 即可求出 P w w w n+ 1 饱和度方程 8 可求出油相饱和度 S 。 o 312 变形方程的求解 利用有效应力原理, 根据虚位移原理及虚功等效原则, 可以建立以有效应 力为基础的单元 e e 等效结点力F 和单元结点位移D 之间的关系式, 即单元平衡方程: e e e e e K D F g + F p - F R 29 0 e 式中 [K ] 是单元刚度矩阵, 其数学表达式为: e T T [K ] [B ] [D ] [B ]dV [B ] [D ] [B ]dx dy dz 30 m m V V 式中 [D ] 为弹性本构矩阵, [B ] 为应变矩阵。 29 式等号右边第一项是单元的自重载荷等效移置到结点上的等效结点力, 表达式为 e T F g [N ] p vdV 31 m V 式中 [N ] 为形函数矩阵, p v 为单位体积的体积力。 设油藏岩土的密度为 Q , 由于油藏岩土为饱含油、水两相流体的多孔介质, 因而其密度为 Q 1 - U Q s + U S oQ o + U S w Q w 32 T 所以有 p [ 0 0 - Q g ] v 33 则自重载荷的等效结点力为 e 0 F gx i 0 e F gy i 34 e - N iQ g dx dy dz F gz i m v 29式等号右边第二项是单元孔隙压力载荷等效移置到结点上的等效结 点力, 表达式为 e T F [B ] P m dV 35 p m V T 因为 m [ 1 1 1 0 0 0 ] 36 N i, x 0 0 N i, y 0 N i, z T [B ] 0 N i, y 0 N i, x N i, z 0 37 i 0 0 N 0 N N i, z i, y i, x 所以1 期 冉启全等: 弹塑性变形油藏中多相渗流的数值模拟 29 5 N i p dx dy dz e m v 5x F p x i 5 N i e F p y i p dx dy dz 38 m v 5y e F p z i 5 N i p dx dy dz m v 5z 0 29 式等号右边第三项是初始应力 R 等效移置到结点上的等效结点力, 表达式为 e T 0 F R [B ] R dV 39 0 m V 由于 0 0 0 0 0 0 0 T R [R R R S S S ] 40 x y z y z z x x y 所以有 5 N i 5 N i 5 N i 0 0 0 Rx + S z x + S x y e 5x 5z 5y F R x i 0 e 5 N i 5 N i 5 N i 0 0 0 F R y i Ry + S y z + S x y dx dy dz 41 0 m V 5y 5z 5x e F R z i 0 5 N 5 N 5 N i i i 0 0 0 Rz + S y z + S z x 5z 5y 5x 在建立了单元刚度矩阵和单元等效结点力, 形成单元平衡方程后, 将这些 方程集合起来, 就可形成总体平衡方程: [K ]D F 42 式中[K ]、 D、 F 分别是总体刚度矩阵、结点位移和结点载荷列阵。对总 体平衡方程 42 式 采用增量初应力法进行求解, 即可得到结点位移。由结点位移, 根据几何方 程即可求出应变。再 由应变, 根据本构方程就可解出应力。 313 耦合求解过程 由于渗流方程和变形方程均含有体现流固耦合效应的参数或变量, 不能对它们进行单独 求解, 故只能联立起来进行耦合求解。耦合求解过程如图 1 所示。 4 数值算例 算例基本参数: 油藏采用衰竭式开采, 油藏中心有一口产油井, 产油井的工作 制度 关于办公室下班关闭电源制度矿山事故隐患举报和奖励制度制度下载人事管理制度doc盘点制度下载 为定压 生产, 井底流压为 5 M Pa。油藏的油层顶深为 2537152m , 油层厚度为 15125m , 开采面积为 - 3 2 500m × 500m , 孔隙度为 25% , 渗透率为 100 × 10 L m , 含油饱和度 88% , 含水饱和度为 12% ,油层压力 30134M Pa。岩石力学特性参数为: 弹性模量 130010M Pa, 泊松比 0115, 内聚力 010 M Pa, 摩擦角 30? , 硬化参数 010。渗流边界条件为封闭边界; 变形边界条件为: 油藏底部垂 直位移为零, 外边界的水平位移为零, 顶部边界为恒定的应力边界 垂直应力为 6210M Pa, 水 平应力为 4314 M Pa 。 利用本文方法对弹塑性变形油藏中油水渗流及开采动态进行了模拟, 同时 将油藏视为线 弹性变形, 进行了对比计算。图 2、图 3和图 4分别是体积应变、孔隙度、渗透率随生产时间变化 的曲线图, 从图中可以看出: 油藏多孔介质在流固耦合作用下要发生变形, 并导致孔隙度、渗透 率的改变, 从而影响油藏的渗流和开采动态。而油藏多孔介质在线弹性变形和弹塑性变形情况 下的体积应变、孔隙度、渗透率改变有较大的差别, 因而直接影响它们的开采动态 图 5 : 模拟计 算 力 学 学 报 16 卷 30 结束时线弹性变形的采出程度为 31196% , 而弹塑性变形的采出程度为 26165%。 图 1 耦合求解过程示意图 图 2 井点网格体积应变与生产时间关系曲线 图 3 井点网格孔隙度与生产时间关系曲线 图 4 井点网格渗透率与生产时间关系曲线 图 5 采出程度与生产时间关系曲线 5 结 论 1 在油藏开采过程中, 流固耦合作用极强, 在用计算机对油藏渗流问题进行数值计算和 模拟时, 不能忽略渗流与变形的耦合效应。 2油藏岩土在注采交变载荷作用下表现为极强的弹塑性特性, 不能将其视 为简单的线弹 性来处理, 否则计算结果将会与实际情况有很大的误差。 3本文给出的弹塑性变形油藏数值模拟模型和方法, 经算例表明是有效的, 对油田开发 有实用价值。