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

本发明专利技术公开一种利用FEM在三维非结构网格流场内追踪流线分布的方法，建立三维非结构化网格的致密砂岩储层数值模拟模型，对任意一条流线的起始点，定位它在四面体网格系统内的位置，计算起始点P

Using FEM to trace the streamline distribution in three-dimensional unstructured grid flow field

【技术实现步骤摘要】
利用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)。进一步的技术本文档来自技高网...

【技术保护点】
1.利用FEM在三维非结构网格流场内追踪流线分布的方法，其特征在于，包括以下步骤：/n步骤S1、利用有限元方法建立三维非结构化网格的致密砂岩储层数值模拟模型，并根据致密砂岩储层的参数获取致密砂岩储层压力场和速度场；/n步骤S2、给定流线数目，根据流场内源的位置确定起始点的分布；/n步骤S3、对任意一条流线i的起始点P

【技术特征摘要】
1.利用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继续追踪该条流线直到该条流线到达汇，然后接着追踪下...

【专利技术属性】
技术研发人员：罗毅，张烈辉，赵玉龙，张博宁，冯胤，刘启国，张东旭，梁海鹏，刘沙，
申请(专利权)人：西南石油大学，成都北方石油勘探开发技术有限公司，
类型：发明
国别省市：四川;51