利用FEM在三维非结构网格流场内追踪流线分布的方法与流程

本发明涉及利用fem在三维非结构网格流场内追踪流线分布的方法，属于非常规油气勘探开发技术领域。

背景技术：

基于流线的数值模拟是一种描述储层中流体流动动态物理过程的工程方法，这种方法广泛应用于石油工程、地下水问题等方面。从上世纪50年代开始，流线模拟引入到石油领域，主要用于追踪流体流动路径和可视化流场内速度分布。流线追踪的准确性对沿着流线求解传质问题来说非常重要。在过去的几十年中，流线模拟方法在石油工业中的油藏数值模拟、生产优化、历史拟合、预测eor采油产量等方面具有广泛的应用。

在追踪流线方面一般使用较多的方法有两类：runge-kutta流线数值追踪方法和pollock半解析流线追踪方法。但是runge-kutta流线数值追踪方法需要给定一个时间步长后迭代计算流线，从而导致runge-kutta追踪流线的效率非常低，不容易生成tof网格，不利于后续基于流线的数值模拟工作的进行。所以，pollock是一种应用非常广泛的流线追踪方法，pollock追踪流线一般分为两步，通过数值模拟方法求解每个网格中的压力，获得压力场和速度场；然后设定流线的起点，根据速度场获得起点的速度，计算该起点离开网格的时间，然后追踪下一个网格，直到该条流线追踪完毕，连接所有网格内的起点和终点，即可追踪获得该条流线。当假设速度场内各方向的速度呈线性变化，且与其它方向的速度无关，pollock提出了分段解析解，该方法是通过tof推导得来的。对任假设流体粒子，给定起始点坐标和速度，通过该算法就能确定离开该网格的出口坐标和穿过该网格用要的时间(这就是tof，是粒子在一个网格内的滞留时间)。这一类方法依赖于网格的结构、计算精度和效率。但是pollock方法时适用于结构化网格，同时当网格内具有源或汇时该方法不能追踪到流线。虽然非结构化网格中流线模拟被许多学者研究，但是这些算法至少存在一个以下缺点：

(1)这些公式依赖于控制体积的概念，需要用到网格各个面上的流量，以及在每个网格中需要一个确定的流量；

(2)在追踪流线的过程中，由于涉及到微分方程的数值求解，无法有效计算并获得tof网格；

(3)大部分非结构网格都是二维的，三维非结构网格内的流线追踪算法还比较少。

因此，目前方法获得的流线只能用来可视化速度场，不能将一维传质方程映射到流线上，不利于后续基于流线的数值模拟工作的展开。为了能够开展基于流线的数值模拟工作，需要开发一种新的三维非结构网格(四面体)流场内流线追踪方法。

技术实现要素：

针对上述问题，本发明主要是克服现有技术中的不足之处，提出利用fem在三维非结构网格流场内追踪流线分布的方法。

本发明解决上述技术问题所提供的技术方案是：利用fem在三维非结构网格流场内追踪流线分布的方法，包括以下步骤：

步骤s1、利用有限元方法建立三维非结构化网格的致密砂岩储层数值模拟模型，并根据致密砂岩储层的参数获取致密砂岩储层压力场和速度场；

步骤s2、给定流线数目，根据流场内源的位置确定起始点的分布；

步骤s3、对任意一条流线i的起始点pi0(x0,y0,z0)，定位它在四面体网格系统内的位置，获得pi0(x0,y0,z0)所在四面体网格的网格编号cellid，并获得该网格的四面体网格顶点坐标(x1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4)和四面体网格顶点对应的速度[vx1,vy1,vz1],[vx2,vy2,vz2],[vx3,vy3,vz3],[vx4,vy4,vz4]；

步骤s4、通过坐标转换公式把起始点pi0(x0,y0,z0)和四面体网格顶点坐标(x1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4)转换为masterelement中的坐标(ξ1,η1,ζ1)，(ξ2,η2,ζ2)，(ξ3,η3,ζ3)，(ξ4,η4,ζ4)，通过jacobian矩阵将四面体网格顶点对应的速度[vx1,vy1,vz1],[vx2,vy2,vz2],[vx3,vy3,vz3],[vx4,vy4,vz4]转化为masterelement中的速度([vξ1,vη1,vζ1],[vξ2,vη2,vζ2],[vξ3,vη3,vζ3],[vξ4,vη4,vζ4])；

步骤s5、在masterelement中利用四面体网格的shapefunction和四面体网格顶点的速度计算起始点pi0(x0,y0,z0)的速度；

步骤s6、在masterelement中，利用shapefunction和建立任意质点速度和位置关系方程，求解该方程获得质点在masterelement中的轨迹方程；

步骤s7、在masterelement中通过newton-raphson迭代算法求解masterelement中的轨迹方程获得质点到达四面体单元四个面的时间，取最小正值为该单元内的飞行时间δτ，并计算出口坐标pe^*(ξe,ηe,ζe)；

步骤s8、利用shapefunction将出口坐标pe^*(ξe,ηe,ζe)转换为真实空间坐标pe(xe,ye,ze)，累加δτ建立tof坐标，连接真实空间内的起点和终点即可获得该网格内的流线线段；

步骤s9、判断pe(xe,ye,ze)是否到达流场内的汇；

步骤s10、如果到达则说明该条流线追踪完毕，则可以追踪下一条流线，否则用pe(xe,ye,ze)替换pi0(x0,y0,z0)，重复步骤s3～s9继续追踪该条流线直到该条流线到达汇，然后接着追踪下一条流线。

进一步的技术方案是，所述步骤s4中四面体网格顶点坐标(x1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4)在真实空间和masterelement空间相互转换；

其中masterelement空间到真实空间的坐标转换公式如下：

x＝x1(1-ξ-η-ζ)+x2ξ+x3η+x4ζ

y＝y1(1-ξ-η-ζ)+y2ξ+y3η+y4ζ

z＝z1(1-ξ-η-ζ)+z2ξ+z3η+z4ζ

真实空间到masterelement空间的坐标转换公式如下：

其中

a1＝(x1-x2),a2＝(x1-x3),a3＝(x1-x4)

b1＝(y1-y2),b2＝(y1-y2),b3＝(y1-y4)

c1＝(z1-z2),c2＝(z1-z3),c3＝(z1-z4)

式中：x为真实空间的x坐标，y为真实空间的y坐标,z为真实空间的z坐标,ξ为masterelement中的与x坐标相对应坐标；η为masterelement中与y坐标相对应的坐标；ζ为masterelement中与z坐标相对应的坐标；(x1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4)分别为真实空间内四面体内四个顶点的坐标，这四个坐标转换到masterelement中为(0,0,0),(1,0,0),(0,1,0),(0,0,1)。

进一步的技术方案是，所述步骤s4中四面体网格顶点对应的速度[vx1,vy1,vz1],[vx2,vy2,vz2],[vx3,vy3,vz3],[vx4,vy4,vz4]从真实空间转换到masterelement空间中，其转换公式如下：

其中jacobian矩阵为：

式中：υx、υy和υy分别为真实空间中x、y和z方向的速度，υξ、υη和υζ为masterelement中ξ、η和ζ方向的速度。

进一步的技术方案是，所述步骤s7的具体过程为：

步骤s71、质点自起点出发，设定该质点能够从四面体的任何一个面离开四面体，并分别计算pi0(x0,y0,z0)从面1(p2-p1-p4)、面2(p2-p1-p3)、面3(p3-p1-p4)、面4(p2-p3-p4)的离开时间δτ1、δτ2、δτ3、δτ4；

步骤s72、取离开时间δτ1、δτ2、δτ3、δτ4的最小正值为单元内的飞行时间δτ；

步骤s73、根据单元内的飞行时间δτ计算masterelement中的出口坐标，将单元内的飞行时间δτ带入case1公式中可以计算获得pe^*(ξe,ηe,ζe)；

其中case1公式如下：

本发明具有以下有益效果：

(1)在流线追踪过程中，只需要获得四面体网格顶点的属性(坐标和流量)，而不需要四面体各个面上的流量；

(2)根据流线的入口坐标能够确定流线在四面体内的出口坐标，这个特点与pollock算法相似；

(3)本发明提出了四面体网格中描述流线轨迹的方程的解析解，能够准确高效的根据入口坐标确定流线的在面体网格内出口坐标和离开网格的时间，累加δτ从而建立tof坐标系统；

(4)能够实现三维非结构网格(四面体网格)流场内的流线追踪。

附图说明

图1为fenics系统中真实空间坐标转换到masterelement空间坐标代码；

图2为fenics系统中masterelement空间坐标转换到真实空间坐标代码；

图3为fenics系统中真实空间速度转换到masterelement空间速度代码；

图4为真实空间和masterelement空间中四面体实体图和坐标分布；

图5为本发明实施例中致密砂岩储层利用有限元追踪流线方法技术路线图；

图6为本发明实施例中三维流场离散网格图；

图7为本发明实施例中三维流场的压力场图；

图8为本发明实施例中三维非结构网格(四面体)流场内的流线分布图(三维视图)。

具体实施方式

下面结合实施例和附图对本发明做更进一步的说明。

实施例1

本发明利用fem在三维非结构网格流场内追踪流线分布的方法，包括以下步骤：

步骤s1、利用有限元方法建立三维非结构化网格的达西流动数值模拟模型，获取三维流场的压力场和速度场；

步骤s2、给定流线数目为100条，流线的起点均匀的分布在源的周围(径向数量为10，法向数量为10)；

步骤s3、对任意一条流线i的起始点pi0(x0,y0,z0)，定位它在四面体网格系统内的位置，获得pi0(x0,y0,z0)所在四面体网格的网格编号cellid，并获得该网格的四面体网格顶点坐标(x1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4)和对应的速度[vx1,vy1,vz1],[vx2,vy2,vz2],[vx3,vy3,vz3],[vx4,vy4,vz4]；

步骤s4、通过坐标转换公式把起始点pi0(x0,y0,z0)和四面体网格顶点坐标(x1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4)转换为masterelement中的坐标(ξ1,η1,ζ1)，(ξ2,η2,ζ2)，(ξ3,η3,ζ3)，(ξ4,η4,ζ4)，通过jacobian矩阵将四面体网格顶点对应的速度转化为masterelement中的速度([vξ1,vη1,vζ1],[vξ2,vη2,vζ2],[vξ3,vη3,vζ3],[vξ4,vη4,vζ4])，转换代码如图3所示；

四面体网格顶点坐标需要在真实空间和masterelement空间相互转换，

其中masterelement空间到真实空间的坐标转换公式(转换代码如图2所示)：

x＝x1(1-ξ-η-ζ)+x2ξ+x3η+x4ζ

y＝y1(1-ξ-η-ζ)+y2ξ+y3η+y4ζ

z＝z1(1-ξ-η-ζ)+z2ξ+z3η+z4ζ(1)

而真实空间到masterelement空间的坐标转换公式为(转换代码如图1所示)：

其中

a1＝(x1-x2),a2＝(x1-x3),a3＝(x1-x4)

b1＝(y1-y2),b2＝(y1-y2),b3＝(y1-y4)

c1＝(z1-z2),c2＝(z1-z3),c3＝(z1-z4)

式中：x为真实空间的x坐标，y为真实空间的y坐标,z为真实空间的z坐标,ξ为masterelement中的与x坐标相对应坐标；η为masterelement中与y坐标相对应的坐标；ζ为masterelement中与z坐标相对应的坐标；(x1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4)分别为真实空间内四面体内四个顶点的坐标，这四个坐标转换到masterelement中为(0,0,0),(1,0,0),(0,1,0),(0,0,1)，如图4所示；

在流线追踪过程中，四面体网格顶点的速度只需要从真实空间转换到masterelement空间中，其转换公式为

其中jacobian矩阵为：

式中：(x1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4)分别为真实空间内四面体四个顶点的坐标；υx、υy和υy分别为真实空间中x、y和z方向的速度，υξ、υη和υζ为masterelement中ξ、η和ζ方向的速度；

步骤s5、在masterelement中利用四面体网格的shapefunction和四面体网格顶点的速度计算起始点pi0(x0,y0,z0)的速度；

步骤s6、在masterelement中，利用shapefunction、建立任意质点速度和位置关系方程，求解该方程获得质点在masterelement中的轨迹方程；

轨迹方程求解过程如下：

masterelement空间是一个正四面体，给定任意一个起点以后pi0(x0,y0,z0)后，pi0(x0,y0,z0)将会从正四面体的一个面离开，所用时间为离开该网格的飞行时间tof，用δτ表示；

所以一阶拉格朗日shapefunction：

s1(ξ,η,ζ)＝1-ξ-η-ζ

s2(ξ,η,ζ)＝ξ

s3(ξ,η,ζ)＝η

s4(ξ,η,ζ)＝ζ(4)

任意点p(ξ,η)速度υξ为

υξ＝s1(ξ,η,ζ)υξ1+s2(ξ,η,ζ)υξ2+s3(ξ,η,ζ)υξ3+s4(ξ,η,ζ)υξ4

υη＝s1(ξ,η,ζ)υη1+s2(ξ,η,ζ)υη2+s3(ξ,η,ζ)υη3+s4(ξ,η,ζ)υη4

υζ＝s1(ξ,η,ζ)υζ1+s2(ξ,η,ζ)υζ2+s3(ξ,η,ζ)υζ3+s4(ξ,η,ζ)υζ4(5)

令

所以方程可以写为

换为矩阵形式

方程8的通解形式为

x(t)＝e1φ1(t)+e2φ2(t)+e3φ3(t)+c(9)

其中特征函数φ1(t),φ2(t),φ3(t)和系数e1,e2,e3andc根据矩阵a的特征值确定，根据a的特征值特征值分布情况，x(t)存在9种不同的情况；

case1三个不相等的非零实数特征根：(0≠r1≠r2≠r3≠0)

case2三个非零实数特征根，其中两个根相等：(0≠r1≠r2＝r3≠0)

apc+b＝0(11)

case3三个相等非零实数特征根：(0≠r1≠r2＝r3≠0)

e1＝p0-pc

e2＝(a-r1i)(p0-pc)

apc+b＝0(12)

case4一个非零实数特征根和两个复数根：(r1≠0,μ+λi,λ≠0)

apc+b＝0(13)

case5一个为零的实数特征根和两个复数根：(r1＝0,μ+λi,λ≠0)

g2＝(μa-(μ²-λ²)i)

case6一个为零的实数特征根和两个非零的实数根：(r1＝0,0≠r2≠r3≠0)

case7两个为零的实数特征根和一个个非零的实数根：(r1＝0,r2＝0,r3≠0)

e1＝ap0+b

case8只有零解：(r1＝0,r2＝0,r3＝0)

x(t)＝e1t+e2t²+e3t³+p0

e1＝ap0+b

case9一个为零的实数特征根和两个相等且不为零的实数特征根：(r1＝0,r2＝0,r3＝0)

步骤s7、在masterelement中通过newton-raphson迭代算法求解masterelement中的轨迹方程获得质点到达四面体单元四个面的时间，取最小正值为该单元内的飞行时间δτ，并计算出口坐标pe*(ξe,ηe,ζe)；

masterelement为一个正四面体(如图4所示)，四个顶点为p1(0,0,0),p2(1,0,0),p3(0,1,0),p4(0,0,1),四个面为：面1(p2-p1-p4),面2(p2-p1-p3),面3(p3-p1-p4),面4(p2-p3-p4)，质点从起点pi0(x0,y0,z0)出发后，经过一定时间后将会从四面体的一个面的一点pe离开四面体，所用的时间即为tof，用δτ表示；

根据该质点在masterelement轨迹方程，出口坐标pe*(ξe,ηe,ζe)和δτ可以通过以下四步计算：

第一步：质点自起点出发，假设该该质点能够从四面体的任何一个面离开四面体，因此分别计算pi0(x0,y0,z0)从面1(p2-p1-p4),面2(p2-p1-p3),面3(p3-p1-p4),面4(p2-p3-p4)的离开时间δτ1、δτ2、δτ3、δτ4；

以case1为例，给定任意一个位置，通过newton-raphson算法就能获得起点pi0(x0,y0,z0)到达该位置的时间，给定ξ＝0，可以确定δτ1；η＝0，可以确定δτ2；ζ＝0，可以确定δτ3；ξ+η+ζ＝1，可以确定δτ4；

第二步：根据计算结果(离开时间δτ1、δτ2、δτ3、δτ4)找出真实的δτ；在给定四面体网格顶点坐标和速度的情况下，根据轨迹方程计算出来的tof只能时最小正值，其余的δτ是无效的；

所以单元内的飞行时间δτ＝min(δτ1，δτ2，δτ3，δτ4)；

第三步：根据单元内的飞行时间δτ计算masterelement中的出口坐标，将单元内的飞行时间δτ带入case1公式中可以计算获得pe*(ξe,ηe,ζe)；

s8、将pe*(ξe,ηe,ζe)通过shapefunction转换到真实空间，得到真实空间下的出口坐标pe(xe,ye,ze)，并累加单元内的飞行时间δτ得到tof(timeofflight)；

将masterelement空间中出口坐标pe*(ξe,ηe,ζe)转换到真实空间出口坐标pe(xe,ye,ze)的公式：

xe＝x1(1-ξe-ηe-ζe)+x2ξe+x3ηe+x3ζe

ye＝y1(1-ξe-ηe-ζe)+y2ξe+y3ηe+y3ζe

ze＝z1(1-ξe-ηe-ζe)+z2ξe+z3ηe+z3ζe；

s9、判断pe(xe,ye,ze)是否到达流场内的汇；

s10、如果到达则说明该条流线追踪完毕，则可以追踪下一条流线，否则用pe(xe,ye,ze)替换pi0(x0,y0,z0)，重复步骤s3到s9继续追踪该条流线直到该条流线到达汇，然后接着追踪下一条流线。

以上步骤可以总结为如图5所示的技术路线图，获得三维流场后(图6为三维流场的离散化网格图，图7为三维流场的压力场图)，以注入井为起点，追踪了100条流线，流线结果如图8所示。

以上所述，并非对本发明作任何形式上的限制，虽然本发明已通过上述实施例揭示，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容作出些变动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。