一种钻井液瞬态波动压力确定方法及系统与流程

本发明涉及钻井技术领域，具体涉及一种钻井液瞬态波动压力确定方法及系统。

背景技术：

这里的陈述仅提供与本发明相关的背景技术，而不必然地构成现有技术。

起下钻或下套管作业引起的波动压力是钻进安全钻井液密度窗口窄地层时井下压力控制的重要敏感因素，引起波动压力的主要因素为：钻井液胶凝强度、钻井液粘性和钻井液惯性。波动压力计算方法有稳态法和瞬态法。发明人发现，稳态法将钻井液看做不可压缩流体，不考虑流道弹性和惯性力，依据管柱运动引起的钻井液流量计算得到的摩阻压降作为波动压力，计算方法相对简单稳定，但计算结果比实测数据偏高30％-50％。瞬态法计算结果与实测数据吻合程度较高，但常规瞬态波动压力计算方法对时间步长和距离步长划分要求高，现场工况一般难以满足时空比要求，计算经常溢出、发散。

技术实现要素：

本发明的目的是为克服现有技术的不足，提供一种钻井液瞬态波动压力确定方法，保证了计算过程的稳定收敛。

为实现上述目的，本发明采用如下技术方案：

第一方面，本发明的实施例提供了一种钻井液瞬态波动压力确定方法，包括以下步骤：

对波动压力物理模型的流动区域进行按照设定的时间步长和距离步长进行网格划分，得到多个网格节点；

获取网格节点原始控制模型的特征线模型，根据得到的特征线模型获得原始控制模型的常微分计算模型；

获取特征线模型的在待求解时刻的前一时刻对应的特征节点的位置信息，根据与特征节点相邻的网格节点的速度和压力获取特征节点速度和压力，将获取的特征节点速度和压力带入常微分计算模型中，得到待计算时刻的网格节点的波动压力计算模型；

利用波动压力计算模型得到不同网格节点在不同时刻的瞬态波动压力。

可选的，根据与特征节点相邻的网格节点的压力和速度信息采用插值算法得到特征节点的压力和速度信息。

可选的，计算边界条件，对流动区域施加边界条件，设定初始条件，对网格节点进行求解变量初始化并按照时间步长，利用波动压力计算模型推演得到不同时刻、不同网格节点处的瞬态波动压力。

可选的，所述流动区域划分为管柱内流动区域、环空流动区域及管柱下方流动区域，运动管柱流动区域为管柱模型内的区域，环空流动区域为管柱模型与井壁模型之间的区域，管柱下方流动区域为管柱模型下方的流动区域。

可选的，计算管柱内流动区域井口节点的压力、速度信息、计算环空流动区域井口节点的压力速度信息、计算管柱内流动区域、环空流动区域和管柱下方流动区域交汇处节点的压力和速度信息，将计算得到的压力和速度信息作为施加的边界条件。

可选的，所述原始控制模型包括质量守恒模型和动量守恒模型，相应的，分别根据质量守恒模型和动量守恒模型获取第一特征线模型和第二特征线模型，根据第一特征线模型得到第一常微分计算模型，根据第二特征线模型得到第二常微分计算模型，根据第一常微分计算模型和第二常微分计算模型得到波动压力计算模型。

可选的，时间步长和距离步长的确定方法为：

获取管柱模型的单管压力波速和流道弹性系数；

对单管距离步长进行初步计算；

对单管距离步长进行调节，根据调节后的单管距离步长得到时间步长；

根据调节后的单管距离部长和管柱的单管长度确定节点数；

根据确定的节点数得到距离步长；

可选的，根据管柱模型运动距离、管柱模型最大速度和加速度确定管柱模型运动时间，根据管柱模型运动时间和时间步长确定时间计算次数。

可选的，采用厚壁筒弹性理论得到管柱模型的流道弹性系数。

一种钻井液瞬态波动压力确定系统，包括：

网格划分模块：被配置为能够对波动压力物理模型的流动区域进行按照设定的时间步长和距离步长进行网格划分；

常微分计算模型获取模块：被配置为能够获取网格节点原始控制模型的特征线模型，根据得到的特征线模型获得原始控制模型的常微分计算模型；

波动压力计算模型获取模块：被配置为能够获取特征线模型的在待求解时刻的前一时刻对应的特征节点的位置信息，根据与特征节点相邻的网格节点的速度和压力信息获取特征节点速度和压力信息，将获取的特征节点速度和压力信息带入常微分计算模型中，得到待计算时刻的网格节点的波动压力计算模型。

波动压力计算模块：被配置为能够利用波动压力计算模型得到不同网格节点在不同时刻的瞬态波动压力。

本发明的有益效果：

1.本发明的方法，利用特征节点的速度和压力信息得到网格节点的波动压力计算模型，无需网格节点全部落在特征线上，允许特征节点的位置不必与网格节点的位置相同，保证了数值计算过程的稳定收敛，提高了计算方法的健壮性，降低了网格划分的要求，避免了数值计算崩溃溢出。

2.本发明的方法，计算边界条件时，三流域汇交部分是三个流动区域的耦合，距离步长和时间步长设定时考虑了流道弹性系数，相较于稳态法将钻井液看做不可压缩和将管柱看做完全刚性体，解决了稳态法计算结果比实测数据偏高30％-50％的问题，计算精度可达90％以上。

附图说明

构成本申请的一部分的说明书附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的限定。

图1为本发明实施例1波动压力物理模型示意图；

图2为本发明实施例1管柱内流动区域网格节点压力和速度计算示意图；

图3为本发明实施例1环空流动区域网格节点压力和速度计算示意图；

图4为本发明实施例1管柱下方流动区域网格节点压力和速度计算示意图；

图5为本发明实施例1钻井系统流动示意图；

图6为本发明实施例1时间步长和距离步长确定流程图；

图7为本发明实施例1流道内节点布置示意图；

图8为本发明实施例1管柱井口节点压力和速度计算示意图；

图9为本发明实施例1管柱底部节点压力和速度计算示意图；

图10为本发明实施例1环空流动区域顶部节点压力和速度计算示意图；

图11为本发明实施例1环空流动区域底部节点压力和速度计算示意图；

图12为本发明实施例1管柱下方流动区域顶部节点压力和速度计算示意图；

图13为本发明实施例1管柱下方流动区域井底节点压力和速度计算示意图；

图14为本发明实施例1波动压力计算流程图；

图15为本发明实施例1变截面连接节点压力和速度计算示意图；

具体实施方式

实施例1

本实施例公开了一种钻井液瞬态波动压力确定方法，包括以下步骤：

步骤1：建立波动压力物理模型，如图1所示，所述波动压力物理模型包括钻井模型，位于钻井模型内的管柱模型及填充在管柱模型和钻井模型之间的钻井液模型，所述波动压力物理模型具有三个流动区域，分别为位于管柱模型内的管柱内流动区域，位于管柱模型下方的管柱下方流动区域及位于管柱和钻井模型井壁之间的环空流动区域。

以管柱底部为z轴原点，建立两套坐标系。管柱内和环空流体z轴向上，管柱下部流体z轴向下

本实施例中，所述波动压力物理模型基于以下假设：

(1)井内水力系统各流道中，泥浆的流动均为一元流动；

(2)钻井液可压缩，钻井液物性受温度和压力影响而改变；

(3)流道和泥浆均为线弹性的，即应力与应变成正比；

(4)3个流动区域的横截面积均受固壁材料的弹性力学参数和流域内压力影响而变化；

(5)管柱底部运动速度受管柱轴向弹性和流体作用在运动管柱壁面上的作用力影响，使得运动管柱底部速度不同于井口速度；

(6)井壁弹性受地层弹性、套管弹性和水泥弹性复合影响；

(7)计算流道中稳定流摩阻压降的公式适用于瞬变流。

步骤2：对波动压力物理模型的流动区域案子设定的时间步长和距离步长划分网格，得到多个网格节点。

对获取波动压力物理模型中流动区域原始控制模型的特征线模型，根据得到的特征线模型获得原始控制模型的常微分计算模型，根据得到的常微分计算模型得到波动压力计算模型。

具体的，原始控制模型包括质量守恒模型和动量守恒模型：

针对管柱内流动区域：

质量守恒模型：

动量守恒模型：

针对环空流动区域：

质量守恒模型：

动量守恒模型：

针对管柱下方流动区域：

质量守恒模型：

动量守恒模型：

其中v为流速；p为总压力；a为流道横截面积；z为轴向坐标；ρ为钻井液密度；δp为摩阻压降，是流体流速和管柱运动速度vp的函数；vp为管柱下入速度，是时间的函数；k为钻井液体积模量，下标1、2、3分别表示三个不同流动区域。

公式(1)-(6)构成了三对一阶拟线性双曲型二元偏微分方程组。

针对管柱内流动区域的原始控制模型：

对原始控制模型进行离散：

原始控制模型中，记则控制方程展开为：

对于管柱内流动区域来说，假定p2已知，即和为常数。

应用双曲型偏微分方程特征方向分析方法，得该方程组的第一特征线模型和第二特征线模型为：

记物理意义为压力波在钻井液中的传播速度。

现场起下管柱引起的管内钻井液运动速度一般为10m/s数量级，一般为1000m/s数量级。故，计算特征方向时可忽略v1项。

将两个原始控制模型进行线性组合，沿特征线模型的特征线方向有；

dz-(v1+c1)dt＝0(10)

a1dp1+[a1ρ1(v1+c1)-a1ρ1v1]dv1+[k(v1+c1)-h]dt＝0(11)

dz-(v1-c)dt＝0(12)

a1dp1+[a1ρ1(v1-c1)-a1ρ1v1]dv1+[k(v1-c1)-h]dt＝0(13)

式中，k＝a1[δp(v1，vp)-ρ1gcosθ]

其中公式(11)和公式(13)分别为第一常微分计算模型和第二常微分计算模型。

将公式(10)-(13)进行简化有：

dz-(v1+c)dt＝0(15)

dz-(v1-c)dt＝0(17)

获取特征线模型的在待求解时刻的前一时刻对应的特征节点的位置信息，根据位置信息获取特征节点速度和压力信息，将获取的特征节点速度和压力信息带入常微分计算模型中，得到待计算时刻的网格节点的波动压力计算模型。

如图2所示，p点为网格节点，tⁿ⁺¹时刻为其待求解时刻，a、q、b分别为tⁿ时刻已知压力和速度的三个相邻的网格节点，设自p点的特征线c+和c-分别与上一时刻交于r和s点，即r和s点分别为特征线模型在待求解时刻的前一时刻对应的特征节点

由特征线模型得，得到r点和s点的z轴坐标即位置信息。

采用插值法获得r点和s点的压力p和速度v

式中，θ＝δtδz，为松弛时空比例因子。

则第一常微分计算模型和第二常微分计算模型一阶近似为：

参数a、b、c、ρ、v等参数下标i表示在节点i和节点i+1形成的流道的参数。

记

ui-1＝a1，i-1ρ1，i-1c1，i-1

ui+1＝a1，i+1ρ1i+1c1，i+1

则波动压力计算模型为：

针对环空流动区域：

质量守恒模型：

动量守恒模型：

对原始控制模型进行离散

原始控制模型中，记则控制模型展开为：

对于环空流动区域来说，假定p1已知，即和为常数。

应用双曲型偏微分方程特征方向分析方法，得到该方程组的第一特征线模型和第二特征线模型为：

记物理意义为压力波在钻井液中的传播速度。

现场下管柱引起的管内钻井液运动速度一般为10m/s数量级，一般为1000m/s数量级，因此，计算特征线模型时可忽略v2项。

将原始控制模型进行线性组合，沿特征线模型的特征线方向有：

dz-(v2+c2)dt＝0(37)

a2dp2+[a2ρ2(v2+c2)-a2ρ2v2]dv2+[k(v2+c2)-h]dt＝0(38)

dz-(v2-c2)dt＝0(39)

a2dp2+[a2ρ2(v2-c2)-a2ρ2ν2]dv2+[k(v2-c2)-h]dt＝0(40)

式中，k＝a2[δp(v2，vp)-ρ2gcosθ]

公式(38)和公式(40)为第一常微分计算模型和第二常微分计算模型。

公式(37)-(40)简化后：

dz-(v2+c2)dt＝0(41)

dz-(v2-c2)dt＝0(43)

如图3所示，p点为网格节点，tⁿ⁺¹时刻为其待求解时刻，a、q、b分别为tⁿ时刻已知压力和速度的三个相邻的网格节点，设自p点的特征线c+和c-分别与上一时刻交于r和s点，即r和s点分别为特征线模型在待求解时刻的前一时刻对应的特征节点

由特征线模型得，得到r点和s点的z轴坐标即位置信息。

采用插值法获得r点和s点的压力p和速度v

式中，θ＝δtδz，为环空松弛时空比例因子。

则第一常微分计算模型和第二常微分计算模型一阶近似为：

参数a、b、c、ρ、v等参数下标i表示在节点i和节点i+1形成的流道的参数。

记

ui-1＝a2，i-1ρ2，i-1c2，i-1

ui+1＝a2，i+1ρ2，i+1c2，i+1

则波动压力计算模型为：

针对管柱下方流动区域

原始控制模型的离散：

该区域内节点原始控制模型的离散形式与管柱内流动区域的离散形式相同，唯一不同之处在于该区域无b1相关项，因此在此不进行重复叙述。

如图4所示，p、a、r、q、s、b节点的设置方式与管柱内流动区域的设置方式相同。

由特征线模型得，得到r点和s点的z轴坐标即位置信息。

则采用插值法获取r点和s点的压力速度信息。

式中θ＝δtδz为松弛时空比例因子。

方程组简化为：

p3，0-p3，s-σ3，1c3，1(v3，0-v3，s)+[-c3，1δp(v3，1)+c3，1ρ3，1gcosθ]δt＝0(65)

dp3+ρ3c3dv3+[c3δp(v3，vp)-c3ρ3gcosθ]dt＝0(66)

dp3-ρ3c3dv3+[-c3δp(v3，vp)+c3ρ3gcosθ]dt＝0(67)

记

ui-1＝p3，i-1c3，i-1

ui+1＝ρ3，i+1c3，i+1

则波动压力计算模型为：

pp＝ps+ui+1vp-ui+1vs-n(69)

划分网格的具体方法为：如图5所示，根据井身结构和运动管柱组合，分析三个流动区域内的流道组成情况。根据套管下入深度在井中的位置和井身结构类型，确定每个流道由几个单管组成，记录单管的外径、内径和长度等信息。如图所示，点划线表示运动管柱下入深度，裸眼井较简单，只有一种情况；套管+裸眼井有三种情况；套管衬管+裸眼井身结构最复杂，共5种情况。以courant准则(δzi/δt≤1)为约束条件，基于用户设定的时间步长，利用系统流道组合分析确定的数据，根据下述算法确定单管i的节点数、距离步长和调整后的时间步长。

如图6所示，所述设定的时间步长和距离步长的确定方法包括以下步骤：

获取管柱模型的单管压力波速ci和流道弹性系数；

对单管距离步长进行初步计算δzi＝(ci+50)δt；

若li＜2δzi，则对单管距离步长进行调节δzi＝li/2，根据调节后的单管距离步长得到时间步长：δt＝δzi/(ci+50)；否则不进行调整。

根据调节后的单管距离部长和管柱的单管长度确定节点数ni＝(int)li/δzi；

根据确定的节点数得到距离步长：δzi＝li/ni；

其中，li为管柱模型的单管i长度，ni为单管节点数，δzi为单管i距离步长。

基于管柱运动距离、管柱运动最大速度和加速度确定管柱运动时间，基于时间步长确定时间计算次数即需要计算多少个时间部长。

nt＝(int)(tt+3)/δt(143)

式中，l为管柱运动距离，vmax为管柱运动最大速度，a为管柱运动加速度。

本实施例中，离散方程为泥浆的压缩系数，为管道的弹性系数，泥浆的压缩系数与泥浆性质，温度，压力有关，一般应通过实测得特定条件下的压缩系数。研究表明，压力在0～49mpa范围内变化时，泥浆压缩系数对波动压力的影响很小，因此，为方便应用，实际计算环空波动压力时，以水在50℃，50mpa下的压缩系数代替泥浆的实际压缩系数：

k＝0，39×10^-9(pa)^-1(137)

流道内节点布置如图7所示，对于流道弹性系数，采用厚壁筒弹性理论进行计算，

1.管柱内流道的弹性系数，仅考虑波动压力引起的流道变形，而管内外静压平衡：

式中r＝d2/d1

2.裸眼井流道弹性系数：

3.套管与套管间环空流道的弹性系数，假设套管材料相同：

4.裸眼井壁与套管间环空流道的弹性系数：

式中es，μs，ef，μf分别为套管和地层的弹性模量和泊松比，计算中取es＝0.2068×10¹²pa，μs＝0.3，ef＝0.17237×10¹¹，μf＝0.28。

本实施例中，考虑了钻井液的可压缩性、管柱和井壁的弹性，相对于稳态法提高了计算精度。

步骤3：计算边界条件，对流动区域施加边界条件。

具体的，计算管柱内流动区域井口节点的压力、速度信息、计算环空流动区域井口节点的压力速度信息、计算管柱内流动区域、环空流动区域和管柱下方流动区域交汇处节点的压力和速度信息，将计算得到的压力和速度信息作为施加的边界条件。

管柱井口节点：

如图8所示，井口节点z方向索引号为is，该节点为求解区域的下游边界，c+为过p点的特征线，可根据该节点原始控制模型得到，根据特征线得到r点的z轴坐标：

进而利用a点和q点已知的压力和速度信息采用插值得到r点的压力和速度信息。

则

井口压力为大气压，表压为0，以表压表示压力，则井口节点的压力和速度即管柱井口的边界条件为：

管柱底部节点

管柱底部节点编号为0，为求解区域的上游边界点，求解节点如图9所示：

c-为过p点的特征线，根据p点的原始控制模型求得，特征节点s的在轴坐标为：

则利用q点和b点的压力和速度信息采用插值法得到s点的压力和速度信息。

则

该节点为运动管柱内流动区域、环空流动区域和运动管柱下部流动区域的汇交点，需要应用汇交点边界条件求解该节点出的压力和流速。

环空流动区域井口节点

如图10所示，c+为过p点的特征线，由p点的原始控制模型得到，特征节点r点的z轴坐标为：

采用插值法获得r点的压力和速度：

则

井口压力为大气压，表压为0，以表压表示压力，则井口节点的边界条件为：

环空流动区域底部节点

环空底部节点z方向索引号为0，节点如图11所示。

过p点的特征线为c-，由p点的原始控制模型获得，特征节点s的z轴坐标为：

采用插值法获得s点的速度和压力：

则

该节点为运动管柱内流动区域、环空流动区域和运动管柱下部流动区域的汇交点，需要应用汇交点边界条件求解该节点出的压力和流速。

管柱下方流动区域顶部节点

运动管柱下方区域的顶部也是运动管柱的底部，设该位置出的节点z方向索引号为0，

如图12所示，过p点的特征线为c-，由p点的原始控制模型求得，特征节点s的z轴坐标为：

采用插值法得到s点的压力和速度：

则

管柱下方流动区域井底节点：

井底节点z方向的索引号为ib，如图13所示，

过p点的特征线为c+，由p点的原始控制模型得到，则特征节点r的z轴坐标为：

采用插值法得到r点的压力和速度：

则

井底处流速vp＝0(vp即为)，得

管柱底部三流动区域汇交区域

管柱流动区域底部节点边界条件满足：

环空流动区域底部节点边界条件满足：

管柱下方流动区域顶部节点边界条件满足：

p3，0-p3，s-ρ3，1c3，1(v3，0-v3，s)+[-c3，1δp(v3，1)+c3，1ρ3，1gcos|δt＝0

井底汇交点处三个流动区域的流量之和为运动管柱单位时间排开的钻井液体积，即

v1，0a1，0+v2，0a2，0+v3，0a3，0＝vp(ap-ae)(101)

环空流动区域底部节点压力等于管柱下方区域顶部节点压力减去钻头和井眼形成的窄环空压降，即

运动管柱内区域底部节点压力等于运动管柱下方区域顶部节点压力减去钻头喷嘴压降，即

式中，aa是钻头与井眼之间的间隙面积，aε为钻头喷嘴截面面积。

公式(99)-公式(103)的六个方程含有六个未知量，p1，0、p2，0、p3，0、v1，0、v2，0、v3，0，且摩阻压降也是未知量的函数，为降低计算难度，用上一时间节点的速度计算摩阻压降梯度，令：

k3＝[-c3，1δp(v3，1)+c3，1ρ3，1gcosθ]δt(106)

则上述六个方程组简化为：

a1，1(p1，0-p1，s)-a1，1ρ1，1c1，1(v1，0-v1，s)+k1+l1v1，0＝0(109)

a2，1(p2，0-p2，s)-a2，1ρ2，1c2，1(v2，0-v2，s)+k2+l2v2，0＝0(110)

p3，0-p3，s-u3，1(v3，0-v3，s)+k3＝0(111)

v1，0a1，0+v2，0a2，0+v3，0a3，0＝vp(ap-ae)(112)

令b1＝-u1，1+l1，b2＝-u2，1+l2，b3＝-u3，1，c1＝u1，1v1，s+k1-a1，1p1，s，c2＝u2，i1v2，s+k2-a2，1p2，s，qt＝vp(ap-ae)则方程(109)-方程(114)简化为：

v1，0a1，0+v2，0a2，0+v3，0a3，0＝vp(ap-ae)(118)

取v2，0最小值vmin为0，最大值vmax为vp(ap-ae)/a2，0，用二分法求解上述6个代数方程组，步骤如下：

①v2，0＝(vmin+vmax)/2；

②将v2，0代入式(116)求得p2，0；

③将p2，0代入式(119)，求得p3，0

④将p3，0代入式(117)，求得v3，0；

⑤将v2，0和v3，0代入式(118)求得v1，0

⑥将v1，0代入式(115)求得

⑦将p3，0和v1，0代入式(120)

⑧比较和若则vmin＝ν2，0，否则，vmax＝v2，0，重复步骤①～⑧，直至

步骤4：如图14所示，输入井身数据(包括各井段套管顶深、套管下深、套管外径、套管内径、井眼尺寸和井眼深度，如果是裸眼段则套管外径输入“-1”。)、输入管柱数据(运动管柱组合外径、内径、长度、运动最大速度vmax和加速度a。)、输入钻井液数据(钻井液密度，钻井液流变性：k、n，钻井液等温压缩系数(缺省取0.39×10^-9pa^-1)。)、输入材料参数(地层弹性模量(缺省取ef＝0.17237×10¹¹pa)，地层泊松比(缺省取0.28)，钢弹性模量(缺省取es＝0.2068×10¹²pa)，钢泊松比(缺省取0.3))。

初始化初始条件，节点压力为依据垂深计算的静液柱压力，节点速度为0，根据划分的网格和施加的边界条件，利用步骤2得到的三个流动区域的波动压力计算模型循环迭代求解各个网格节点在不同时刻的速度和压力，压力即为波动压力。

本实施例中，如图15所示，管柱内流动区域、环空流动区域和管柱下方流动区域存在变截面。

管柱内流动区域和环空流动区域变截面处连接点的离散方程相同，如图11所示，节点i位于流道变截面处，该节点处存在两个不同特征线确定的流速vpu和vpd以及1个压力pp，沿特征线得差分方程：

auvpv＝advpd-(au-ad)vp(123)

解得波动压力计算模型为：

管柱下方流动区域连接点离散方程：

auvpu＝advpd(129)

pp＝ps+ui+1vpd-ui+1vs-n(132)

ui-1＝ρ3，i-1c3，i-1(135)

ui+1＝ρ3，i+1c3，i+1(136)

实施例2：

本实施例公开了一种钻井液瞬态波动压力确定系统，包括：

网格划分模块：被配置为能够对波动压力物理模型的流动区域进行按照设定的时间步长和距离步长进行网格划分；

常微分计算模型获取模块：被配置为能够获取网格节点原始控制模型的特征线模型，根据得到的特征线模型获得原始控制模型的常微分计算模型；

波动压力计算模型获取模块：被配置为能够获取特征线模型的在待求解时刻的前一时刻对应的特征节点的位置信息，根据与特征节点相邻的网格节点的速度和压力信息获取特征节点速度和压力信息，将获取的特征节点速度和压力信息带入常微分计算模型中，得到待计算时刻的网格节点的波动压力计算模型；

波动压力计算模块：被配置为能够利用波动压力计算模型得到不同网格节点在不同时刻的瞬态波动压力。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。