CN111210522A - 利用fem在三维非结构网格流场内追踪流线分布的方法 - Google Patents

利用fem在三维非结构网格流场内追踪流线分布的方法 Download PDF

Info

Publication number
CN111210522A
CN111210522A CN202010035803.6A CN202010035803A CN111210522A CN 111210522 A CN111210522 A CN 111210522A CN 202010035803 A CN202010035803 A CN 202010035803A CN 111210522 A CN111210522 A CN 111210522A
Authority
CN
China
Prior art keywords
coordinate
coordinates
real space
masterelement
streamline
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010035803.6A
Other languages
English (en)
Other versions
CN111210522B (zh
Inventor
罗毅
张烈辉
赵玉龙
张博宁
冯胤
刘启国
张东旭
梁海鹏
刘沙
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Chengdu North Oil Exploration Development Technology Co ltd
Southwest Petroleum University
Original Assignee
Chengdu North Oil Exploration Development Technology Co ltd
Southwest Petroleum University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Chengdu North Oil Exploration Development Technology Co ltd, Southwest Petroleum University filed Critical Chengdu North Oil Exploration Development Technology Co ltd
Priority to CN202010035803.6A priority Critical patent/CN111210522B/zh
Publication of CN111210522A publication Critical patent/CN111210522A/zh
Application granted granted Critical
Publication of CN111210522B publication Critical patent/CN111210522B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T19/00—Manipulating 3D models or images for computer graphics
    • G06T19/20—Editing of 3D images, e.g. changing shapes or colours, aligning objects or positioning parts

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Computer Graphics (AREA)
  • Software Systems (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Architecture (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开一种利用FEM在三维非结构网格流场内追踪流线分布的方法，建立三维非结构化网格的致密砂岩储层数值模拟模型，对任意一条流线的起始点，定位它在四面体网格***内的位置，计算起始点Pi0在Master Element空间的坐标和速度；获得质点离开Master Element空间内四面体网格的时间和出口坐标；将出口坐标转换为真实空间坐标，连接真实空间内的起点和终点即可获得该网格内的流线线段；判断Pe是否到达流场内的汇，如果到达则说明该条流线追踪完毕，则可以追踪下一条流线，否则用Pe替换Pi0，重复上述步骤继续追踪该条流线直到该条流线到达汇，然后接着追踪下一条流线。本发明能够准确追踪四面体非结构化网格***内的流线分布情况。

Description

利用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)所在四面体网格的网格编号Cell ID，并获得该网格的四面体网格顶点坐标(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)转换为Master Element中的坐标(ξ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]转化为Master Element中的速度([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、在Master Element中利用四面体网格的Shape Function和四面体网格顶点的速度计算起始点Pi0(x0,y0,z0)的速度；
步骤S6、在Master Element中，利用Shape Function和建立任意质点速度和位置关系方程，求解该方程获得质点在MasterElement中的轨迹方程；
步骤S7、在Master Element中通过Newton-Raphson迭代算法求解Master Element中的轨迹方程获得质点到达四面体单元四个面的时间，取最小正值为该单元内的飞行时间Δτ，并计算出口坐标Pe *(ξe,ηe,ζe)；
步骤S8、利用Shape Function将出口坐标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空间的坐标转换公式如下：
Figure BDA0002365951470000041
其中
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空间中，其转换公式如下：
Figure BDA0002365951470000042
其中Jacobian矩阵为：
Figure BDA0002365951470000051
Figure BDA0002365951470000052
式中：υ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公式如下：
Figure BDA0002365951470000053
本发明具有以下有益效果：
(1)在流线追踪过程中，只需要获得四面体网格顶点的属性(坐标和流量)，而不需要四面体各个面上的流量；
(2)根据流线的入口坐标能够确定流线在四面体内的出口坐标，这个特点与Pollock算法相似；
(3)本发明提出了四面体网格中描述流线轨迹的方程的解析解，能够准确高效的根据入口坐标确定流线的在面体网格内出口坐标和离开网格的时间，累加Δτ从而建立TOF坐标***；
(4)能够实现三维非结构网格(四面体网格)流场内的流线追踪。
附图说明
图1为FEniCS***中真实空间坐标转换到MasterElement空间坐标代码；
图2为FEniCS***中Master Element空间坐标转换到真实空间坐标代码；
图3为FEniCS***中真实空间速度转换到MasterElement空间速度代码；
图4为真实空间和Master Element空间中四面体实体图和坐标分布；
图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所示)：
Figure BDA0002365951470000071
其中
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空间中，其转换公式为
Figure BDA0002365951470000081
其中Jacobian矩阵为：
Figure BDA0002365951470000082
Figure BDA0002365951470000083
式中：(x1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4)分别为真实空间内四面体四个顶点的坐标；υx、υy和υy分别为真实空间中x、y和z方向的速度，υξ、υη和υζ为MasterElement中ξ、η和ζ方向的速度；
步骤S5、在MasterElement中利用四面体网格的Shape Function和四面体网格顶点的速度计算起始点Pi0(x0,y0,z0)的速度；
步骤S6、在Master Element中，利用Shape Function、建立任意质点速度和位置关系方程，求解该方程获得质点在MasterElement中的轨迹方程；
轨迹方程求解过程如下：
Master element空间是一个正四面体，给定任意一个起点以后Pi0(x0,y0,z0)后，Pi0(x0,y0,z0)将会从正四面体的一个面离开，所用时间为离开该网格的飞行时间TOF，用Δτ表示；
所以一阶拉格朗日shape function：
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)
Figure BDA0002365951470000091
令
Figure BDA0002365951470000092
所以方程可以写为
Figure BDA0002365951470000101
换为矩阵形式
Figure BDA0002365951470000102
方程8的通解形式为
X(t)＝E1Φ1(t)+E2Φ2(t)+E3Φ3(t)+C (9)
其中特征函数Φ1(t),Φ2(t),Φ3(t)和系数E1,E2,E3 and C根据矩阵A的特征值确定，根据A的特征值特征值分布情况，X(t)存在9种不同的情况；
Case1三个不相等的非零实数特征根：(0≠r1≠r2≠r3≠0)
Figure BDA0002365951470000103
Case2三个非零实数特征根，其中两个根相等：(0≠r1≠r2＝r3≠0)
Figure BDA0002365951470000111
Figure BDA0002365951470000112
Figure BDA0002365951470000113
Figure BDA0002365951470000114
APc+B＝0 (11)
Case3三个相等非零实数特征根：(0≠r1≠r2＝r3≠0)
Figure BDA0002365951470000115
E1＝P0-Pc
E2＝(A-r1I)(P0-Pc)
Figure BDA0002365951470000116
APc+B＝0 (12)
Case4一个非零实数特征根和两个复数根：(r1≠0,μ+λi,λ≠0)
Figure BDA0002365951470000117
Figure BDA0002365951470000118
Figure BDA0002365951470000119
Figure BDA00023659514700001110
APc+B＝0 (13)
Case5一个为零的实数特征根和两个复数根：(r1＝0,μ+λi,λ≠0)
Figure BDA0002365951470000121
Figure BDA0002365951470000122
Figure BDA0002365951470000123
Figure BDA0002365951470000124
Figure BDA0002365951470000125
G2＝(μA-(μ2-λ2)I)
Figure BDA0002365951470000126
Case6一个为零的实数特征根和两个非零的实数根：(r1＝0,0≠r2≠r3≠0)
Figure BDA0002365951470000127
Figure BDA0002365951470000128
Figure BDA0002365951470000129
Figure BDA00023659514700001210
Case7两个为零的实数特征根和一个个非零的实数根：(r1＝0,r2＝0,r3≠0)
Figure BDA00023659514700001211
E1＝AP0+B
Figure BDA00023659514700001212
Figure BDA00023659514700001213
Case8只有零解：(r1＝0,r2＝0,r3＝0)
X(t)＝E1t+E2t2+E3t3+P0
E1＝AP0+B
Figure BDA0002365951470000131
Figure BDA0002365951470000132
Case9一个为零的实数特征根和两个相等且不为零的实数特征根：(r1＝0,r2＝0,r3＝0)
Figure BDA0002365951470000133
Figure BDA0002365951470000134
Figure BDA0002365951470000135
Figure BDA0002365951470000136
步骤S7、在Master Element中通过Newton-Raphson迭代算法求解Master Element中的轨迹方程获得质点到达四面体单元四个面的时间，取最小正值为该单元内的飞行时间Δτ，并计算出口坐标Pe*(ξe,ηe,ζe)；
Master element为一个正四面体(如图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，用Δτ表示；
根据该质点在Master Element轨迹方程，出口坐标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；
Figure BDA0002365951470000141
Figure BDA0002365951470000142
Figure BDA0002365951470000143
Figure BDA0002365951470000144
第二步：根据计算结果(离开时间Δτ1、Δτ2、Δτ3、Δτ4)找出真实的Δτ；在给定四面体网格顶点坐标和速度的情况下，根据轨迹方程计算出来的TOF只能时最小正值，其余的Δτ是无效的；
所以单元内的飞行时间Δτ＝min(Δτ1，Δτ2，Δτ3，Δτ4)；
第三步：根据单元内的飞行时间Δτ计算Master Element中的出口坐标，将单元内的飞行时间Δτ带入case1公式中可以计算获得Pe*(ξe,ηe,ζe)；
Figure BDA0002365951470000145
Figure BDA0002365951470000146
Figure BDA0002365951470000147
S8、将Pe*(ξe,ηe,ζe)通过Shape Function转换到真实空间，得到真实空间下的出口坐标Pe(xe,ye,ze)，并累加单元内的飞行时间Δτ得到TOF(Time ofFlight)；
将Master Element空间中出口坐标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所示。
以上所述，并非对本发明作任何形式上的限制，虽然本发明已通过上述实施例揭示，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容作出些变动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims (4)

1.利用FEM在三维非结构网格流场内追踪流线分布的方法，其特征在于，包括以下步骤：
步骤S1、利用有限元方法建立三维非结构化网格的致密砂岩储层数值模拟模型，并根据致密砂岩储层的参数获取致密砂岩储层压力场和速度场；
步骤S2、给定流线数目，根据流场内源的位置确定起始点的分布；
步骤S3、对任意一条流线i的起始点Pi0(x0,y0,z0)，定位它在四面体网格***内的位置，获得Pi0(x0,y0,z0)所在四面体网格的网格编号Cell ID，并获得该网格的四面体网格顶点坐标(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)转换为Master Element中的坐标(ξ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]转化为Master Element中的速度([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、在Master Element中利用四面体网格的Shape Function和四面体网格顶点的速度计算起始点Pi0(x0,y0,z0)的速度；
步骤S6、在Master Element中，利用Shape Function和建立任意质点速度和位置关系方程，求解该方程获得质点在MasterElement中的轨迹方程；
步骤S7、在Master Element中通过Newton-Raphson迭代算法求解Master Element中的轨迹方程获得质点到达四面体单元四个面的时间，取最小正值为该单元内的飞行时间Δτ，并计算出口坐标Pe *(ξe,ηe,ζe)；
步骤S8、利用Shape Function将出口坐标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继续追踪该条流线直到该条流线到达汇，然后接着追踪下一条流线。
2.根据权利要求1所述的利用FEM在三维非结构网格流场内追踪流线分布的方法，其特征在于，所述步骤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空间的坐标转换公式如下：
Figure FDA0002365951460000021
其中
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)分别为真实空间内四面体内四个顶点的坐标。
3.根据权利要求1所述的利用FEM在三维非结构网格流场内追踪流线分布的方法，其特征在于，所述步骤S4中四面体网格顶点对应的速度[vx1,vy1,vz1],[vx2,vy2,vz2],[vx3,vy3,vz3],[vx4,vy4,vz4]从真实空间转换到MasterElement空间中，其转换公式如下：
Figure FDA0002365951460000031
其中Jacobian矩阵为：
Figure FDA0002365951460000032
Figure FDA0002365951460000033
式中：υx、υy和υy分别为真实空间中x、y和z方向的速度，υξ、υη和υζ为MasterElement中ξ、η和ζ方向的速度。
4.根据权利要求1所述的利用FEM在三维非结构网格流场内追踪流线分布的方法，其特征在于，所述步骤S7的具体过程为：
步骤S71、质点自起点出发，设定该质点能够从四面体的任何一个面离开四面体，并分别计算Pi0(x0,y0,z0)从面1、面2、面3、面4离开时间Δτ1、Δτ2、Δτ3、Δτ4；
步骤S72、取离开时间Δτ1、Δτ2、Δτ3、Δτ4的最小正值为单元内的飞行时间Δτ；
步骤S73、根据单元内的飞行时间Δτ计算MasterElement中的出口坐标，将单元内的飞行时间Δτ带入case1公式中可以计算获得Pe *(ξe,ηe,ζe)。
CN202010035803.6A 2020-01-14 2020-01-14 利用fem在三维非结构网格流场内追踪流线分布的方法 Active CN111210522B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010035803.6A CN111210522B (zh) 2020-01-14 2020-01-14 利用fem在三维非结构网格流场内追踪流线分布的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010035803.6A CN111210522B (zh) 2020-01-14 2020-01-14 利用fem在三维非结构网格流场内追踪流线分布的方法

Publications (2)

Publication Number Publication Date
CN111210522A true CN111210522A (zh) 2020-05-29
CN111210522B CN111210522B (zh) 2021-03-19

Family

ID=70786736

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010035803.6A Active CN111210522B (zh) 2020-01-14 2020-01-14 利用fem在三维非结构网格流场内追踪流线分布的方法

Country Status (1)

Country Link
CN (1) CN111210522B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116861753A (zh) * 2023-07-28 2023-10-10 长江大学 一种基于模拟有限差分法的油水两相流线模拟新方法

Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102156779A (zh) * 2011-04-13 2011-08-17 北京石油化工学院 地下水流仿真与预测分析方法
CN102323965A (zh) * 2011-08-25 2012-01-18 浙江大学 面向非定常三维流场涡结构的智能分析方法
CN102495427A (zh) * 2011-12-10 2012-06-13 北京航空航天大学 一种基于隐式模型表达的界面感知射线追踪方法
US20140334190A1 (en) * 2013-05-07 2014-11-13 University Of Central Florida Research Foundation, Inc. Variable frequency iteration mppt for resonant power converters
CN105069826A (zh) * 2015-08-26 2015-11-18 中国科学院深圳先进技术研究院 弹性物体变形运动的建模方法
CN108170895A (zh) * 2016-12-06 2018-06-15 富士通株式会社 纹线可视化设备和方法
CN108241777A (zh) * 2017-12-27 2018-07-03 青岛海洋地质研究所 基于非结构网格有限元法计算水合物沉积物中渗流速度场的方法
CN108846224A (zh) * 2018-06-27 2018-11-20 中国人民解放军国防科技大学 一种超声速流道设计方法及装置
CN108875140A (zh) * 2018-05-24 2018-11-23 西安石油大学 一种基于数字岩心模型的稠油油藏沥青质沉积吸附损害模拟方法
WO2018216931A1 (ko) * 2017-05-24 2018-11-29 필드지 주식회사 와류 생성핀을 포함하는 선박
CN109635353A (zh) * 2018-11-17 2019-04-16 天津大学 一种厨师操作对抽油烟机油烟捕集效率影响的评价方法

Patent Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102156779A (zh) * 2011-04-13 2011-08-17 北京石油化工学院 地下水流仿真与预测分析方法
CN102323965A (zh) * 2011-08-25 2012-01-18 浙江大学 面向非定常三维流场涡结构的智能分析方法
CN102495427A (zh) * 2011-12-10 2012-06-13 北京航空航天大学 一种基于隐式模型表达的界面感知射线追踪方法
US20140334190A1 (en) * 2013-05-07 2014-11-13 University Of Central Florida Research Foundation, Inc. Variable frequency iteration mppt for resonant power converters
CN105069826A (zh) * 2015-08-26 2015-11-18 中国科学院深圳先进技术研究院 弹性物体变形运动的建模方法
CN108170895A (zh) * 2016-12-06 2018-06-15 富士通株式会社 纹线可视化设备和方法
WO2018216931A1 (ko) * 2017-05-24 2018-11-29 필드지 주식회사 와류 생성핀을 포함하는 선박
CN108241777A (zh) * 2017-12-27 2018-07-03 青岛海洋地质研究所 基于非结构网格有限元法计算水合物沉积物中渗流速度场的方法
CN108875140A (zh) * 2018-05-24 2018-11-23 西安石油大学 一种基于数字岩心模型的稠油油藏沥青质沉积吸附损害模拟方法
CN108846224A (zh) * 2018-06-27 2018-11-20 中国人民解放军国防科技大学 一种超声速流道设计方法及装置
CN109635353A (zh) * 2018-11-17 2019-04-16 天津大学 一种厨师操作对抽油烟机油烟捕集效率影响的评价方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
FENG YIN 等: "Implementation of streamline simulation based on finite element method in FEniCS", 《COMPUTATIONAL GEOSCIENCES》 *
武强 等: "地下水流场三维流线可视化模拟与实现", 《中国煤炭地质》 *
王学德 等: "三维非结构网格DSMC并行算法及应用研究", 《宇航学报》 *
王昉: "非定常流场拉格朗日拟序结构可视计算关键技术研究", 《中国博士学位论文全文数据库 信息科技辑》 *
田书玲: "基于非结构网格方法的重叠网格算法研究", 《中国博士学位论文全文数据库 信息科技辑》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116861753A (zh) * 2023-07-28 2023-10-10 长江大学 一种基于模拟有限差分法的油水两相流线模拟新方法
CN116861753B (zh) * 2023-07-28 2024-05-03 长江大学 一种基于模拟有限差分法的油水两相流线模拟新方法

Also Published As

Publication number Publication date
CN111210522B (zh) 2021-03-19

Similar Documents

Publication Publication Date Title
CN108984920B (zh) 一种发动机冷却水套的直接流固耦合传热分析方法
CN111709171B (zh) 一种热流强耦合问题的等几何求解及散热拓扑生成方法
US10296672B2 (en) Generating inviscid and viscous fluid-flow simulations over an aircraft surface using a fluid-flow mesh
Majumdar et al. Three-dimensional finite-volume method for incompressible flows with complex boundaries
CN106503396B (zh) 基于有限差分法与有限体积法耦合的多维水力***瞬变模拟方法
CN114077802B (zh) 一种利用形函数插值替代核函数近似的粒子建模方法
CN107341322A (zh) 一种在线监测核级设备和管道疲劳损伤的方法
CN111210522B (zh) 利用fem在三维非结构网格流场内追踪流线分布的方法
CA2964250A1 (en) Junction models for simulating proppant transport in dynamic fracture networks
CN108595782B (zh) 一种离散裂缝中基质与裂缝间的传质计算方法
CN114117877A (zh) 一种基于等几何粒子描述的拓扑优化方法
KR102266279B1 (ko) 비정상상태를 구현하기 위한 차수 감축 모델 구축 방법
Cummings et al. Supersonic, turbulent flow computation and drag optimization for axisymmetric afterbodies
Sergio et al. Optimization methodology assessment for the inlet velocity profile of a hydraulic turbine draft tube: part II—performance evaluation of draft tube model
US9087165B2 (en) Automatic extremum detection on a surface mesh of a component
CN106294913A (zh) 提高零部件热分析计算结果可靠性的方法
CN110717295B (zh) 一种利用有限元方法追踪致密砂岩储层流线分布的方法
Saleel et al. On simulation of backward facing step flow using immersed boundary method
CN110909511B (zh) 一种无曲面体积分的无粘低速绕流数值模拟方法
CN106066912A (zh) 一种多分区结构化网格的生成方法
Lakshminarayan et al. Fully Automated Surface Mesh Adaptation in Strand Grid Framework
CN114707384B (zh) 一种聚变堆包层统一三维模型非结构网格核热耦合计算方法
White et al. An Overview of the Development of Grid Adaptation Related Technologies and Tools for the VULCAN-CFD Code
CN113420392B (zh) 一种基于流道轨迹优化的共轭传热散热器设计方法
CN114036815B (zh) 面向耦合物理场快速求解的真伪双重粒子模型建模方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant