CN113324840A - 一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法 - Google Patents

一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法 Download PDF

Info

Publication number
CN113324840A
CN113324840A CN202110604358.5A CN202110604358A CN113324840A CN 113324840 A CN113324840 A CN 113324840A CN 202110604358 A CN202110604358 A CN 202110604358A CN 113324840 A CN113324840 A CN 113324840A
Authority
CN
China
Prior art keywords
rock
well
solid
finite element
heterogeneous
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
CN202110604358.5A
Other languages
English (en)
Other versions
CN113324840B (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.)
Southwest Petroleum University
Original Assignee
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 Southwest Petroleum University filed Critical Southwest Petroleum University
Priority to CN202110604358.5A priority Critical patent/CN113324840B/zh
Publication of CN113324840A publication Critical patent/CN113324840A/zh
Application granted granted Critical
Publication of CN113324840B publication Critical patent/CN113324840B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N3/00—Investigating strength properties of solid materials by application of mechanical stress
    • G01N3/08—Investigating strength properties of solid materials by application of mechanical stress by applying steady tensile or compressive forces
    • G01N3/10—Investigating strength properties of solid materials by application of mechanical stress by applying steady tensile or compressive forces generated by pneumatic or hydraulic pressure
    • G01N3/12—Pressure testing
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2203/00—Investigating strength properties of solid materials by application of mechanical stress
    • G01N2203/0014—Type of force applied
    • G01N2203/0016—Tensile or compressive
    • G01N2203/0019—Compressive
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2203/00—Investigating strength properties of solid materials by application of mechanical stress
    • G01N2203/003—Generation of the force
    • G01N2203/0042—Pneumatic or hydraulic means
    • G01N2203/0048—Hydraulic means
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2203/00—Investigating strength properties of solid materials by application of mechanical stress
    • G01N2203/02—Details not specific for a particular testing method
    • G01N2203/025—Geometry of the test
    • G01N2203/0252—Monoaxial, i.e. the forces being applied along a single axis of the specimen
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2203/00—Investigating strength properties of solid materials by application of mechanical stress
    • G01N2203/02—Details not specific for a particular testing method
    • G01N2203/025—Geometry of the test
    • G01N2203/0256—Triaxial, i.e. the forces being applied along three normal axes of the specimen

Landscapes

  • Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • General Health & Medical Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Investigation Of Foundation Soil And Reinforcement Of Foundation Soil By Compacting Or Drainage (AREA)

Abstract

本发明公开一种非均质地层井壁渐进坍塌过程的流‑固‑热耦合模拟方法，包括确定地层地应力、孔隙压力、井筒压力；确定岩石力学参数样本数量、岩石基础物性参数和岩石力学参数；确定每个参数的最佳分布函数及对应的特征参数；建立井周应力分布的流‑固‑热耦合有限元数学模型；计算井周应力分布；计算损伤变量F；损伤变量F确定破坏区域；重复步骤进行迭代，直到损伤变量F≤0，并绘制井壁失稳区域图，确定每次迭代步下的NYZA值。本发明使用输入参数的不确定性，表示了岩石力学参数的非均质性；采用连续损伤理论，更改已破坏区域的力学参数，视为该区域等效挖去，模拟了井壁渐进破坏发展的四个阶段：裂纹萌生、扩展、崩落形成和稳定。

Description

一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法
技术领域
本发明涉及一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法，属于石油钻井工程、岩石力学与工程领域。
背景技术
钻开地层形成井眼会破坏原应力平衡，导致井壁应力集中，钻井液压力不足以支持井壁围压发生井壁失稳。井壁失稳是一个渐进的过程，最初圆形井眼向最小主应力方向的逐级扩展破坏进而释放应力集中。而破坏程度较低并能保持正常的钻井作业，则称为井壁崩落；破坏程度较高，则称为井壁坍塌。Zoback 指出直井的井壁崩落的方向与最小水平地应力方向一致，室内实验也证明了井壁崩落的方向与最小应力方向一致。
最早为了模拟井壁崩落这现象，利用Kirsch方程计算应力并用摩尔库伦准则判断井周围岩是否破坏，这一方法可以精确的模拟井壁开始破坏阶段。而实际井壁崩落过程，井壁围岩破坏后，井周应力会发生改变，新的岩石将会受到新的应力集中，形成新的破坏区，重复这一过程，直到井壁不再失稳。而这一渐进破坏的过程很难使用解析解计算。现有的数值模拟(有限元法、离散元法和边界元法)研究了井壁渐进破坏发展的四个阶段：裂纹萌生、扩展、崩落形成和稳定。结果表明井壁崩落深度在逐渐增大直至稳定，而崩落宽度在第一次破坏的基础上保持不变；井壁崩落是一系列平行于最小主应力方向的剪切破坏引起的。但岩石都被假设为线弹性、均质和各向同性的材料，实际上，岩石作为天然非均质材料，在其力学行为具有非均质性。此外，现有的井壁渐进破坏数值模拟(有限元法、离散元法和边界元法)仅仅考虑了岩石材料变形破坏特征，忽略了流-固-热多场耦合的影响。
为此，本发明提出一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法。使用输入岩石力学参数的不确定性，表示了岩石力学参数的非均质性；采用连续损伤理论，更改已破坏区域的力学参数，视为该区域等效挖去，模拟了井壁渐进破坏发展的四个阶段：裂纹萌生、扩展、崩落形成和稳定。
发明内容
本发明主要是克服现有技术中的不足之处，提出一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法。
本发明解决上述技术问题所提供的技术方案是：一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法，包括以下步骤：
步骤S10、根据室内实验和测井资料确定目标井的地层地应力、孔隙压力、井筒压力；
步骤S20、根据室内岩石力学实验确定目标井地层岩石的岩石力学参数样本数量、岩石基础物性参数和岩石力学参数；
步骤S30、根据岩石基础物性参数和岩石力学参数确定每个参数的最佳分布函数及对应的特征参数；
步骤S40、建立井周应力分布的流-固-热耦合有限元数学模型；
步骤S50、根据井周应力分布流-固耦合有限元数学模型计算井周应力分布；
步骤S60、根据井周应力分布计算损伤变量F；
Figure RE-GDA0003131971070000021
F1＝-σ3-St
Figure RE-GDA0003131971070000022
式中：C为岩石内聚力；
Figure RE-GDA0003131971070000023
为岩石内摩擦角；St是岩石抗拉强度；F1为最大拉应力状态函数；F2为Mohr-Coulomb状态函数；F为损伤变量；σ1为最大主应力；σ3为最小主应力；
步骤S70、损伤变量F确定破坏区域，更新破坏区域的杨氏模量、孔隙度、渗透率；
步骤S80、重复步骤S50至步骤S70进行迭代，直到损伤变量F≤0，此时，无新的损伤区出现，井眼最终趋于稳定，实现非均质地层井壁渐进坍塌过程的数值模拟，并绘制井壁失稳区域图，确定每次迭代步下的NYZA值。
进一步的技术方案是，所述地层地应力包括最大水平地应力、最小水平地应力。
进一步的技术方案是，所述岩石力学实验包括单轴压缩和三轴压缩岩石力学实验；所述岩石基础物性参数包括孔隙度、渗透率和Biot系数；所述岩石力学参数包括杨氏模量、泊松比、内摩擦角、内聚力。
进一步的技术方案是，所述步骤S40的具体过程为：
步骤S41、建立控制方程：
Figure RE-GDA0003131971070000031
Figure RE-GDA0003131971070000032
Figure RE-GDA0003131971070000033
式中，G和λ是拉梅常数；k是多孔介质渗透率；u是流体的粘度；u、p和T 分别是多孔介质的位移、孔隙压力和温度；下标t表示时间的导数；φ是多孔介质的孔隙度；Kf、Km分别为流体和岩石的体积模量；IT＝[1,1,1,0,0,0]；D是弹性刚度矩阵；λT是多孔介质基质导热系数；ρscs是多孔介质基质热熔；ρfcf是流体热熔；βs是多孔介质基质热膨胀系数；βf是流体热膨胀系数；
步骤S41、使用伽辽金法有限元法逼近控制方程，可以得到控制方程的有限元求解格式：
Figure RE-GDA0003131971070000041
其中，
M＝∫VBTDBdV
Figure RE-GDA0003131971070000042
Figure RE-GDA0003131971070000043
Figure RE-GDA0003131971070000044
Figure RE-GDA0003131971070000045
Figure RE-GDA0003131971070000046
Figure RE-GDA0003131971070000047
Figure RE-GDA0003131971070000048
Figure RE-GDA0003131971070000049
Figure RE-GDA00031319710700000410
B＝LNu
式中，B为应变和位移相关的应变矩阵；上标T为矩阵转置；N有限元形函数；u、p和T分别为未知变量u、p和T的向量；ut、pt、Tt分别为未知变量u、p和T的时间导数；fu、fp、fT分别是节点载荷向量、流体源汇向量和热源向量；L是微分算子；Nu是位移形函数。
进一步的技术方案是，所述步骤S50的具体过程为：根据井周应力分布流- 固耦合有限元数学模型建立非均质地层井壁渐进坍塌过程数值模拟的有限元模型，分别对非均质地层有限元模型赋值非均质材料参数、施加边界条件和划分有限元网格，最后计算井周应力分布，得到最大主应力σ1和最小主应力σ3。
进一步的技术方案是，根据井周应力分布流-固耦合有限元数学模型建立非均质地层井壁渐进坍塌过程数值模拟的有限元模型包括：根据平面应变和轴对称条件、井周应力分布流-固耦合有限元数学模型，建立非均质地层井壁渐进坍塌过程数值模拟的有限元模型，其井筒半径为R，整个模型为50R×50R的几何模型，以排除边界效应的影响，并为了减少计算量。
进一步的技术方案是，所述非均质地层井壁渐进坍塌过程数值模拟有限元模型的具体非均质材料参数赋值过程为：
首先选取合适每个参数的分布方式并编译MATLAB外部随机分布函数，其中随机分布函数以空间位置作为因变量；
然后把MATLAB外部随机分布函数放置COMSOL非均质地层有限元程序同一文件夹；
再点击COMSOL Link MATLAB，把MATLAB与COMSOL链接起，即可以调用已编译好的MATLAB函数；
最后在COMSOL中固体力学接口中，给内聚力、内摩擦角、杨氏模量、泊松比和Biot系数输入相对应的MATLAB函数名称；在COMSOL中流体力学接口中给孔隙度和渗透率输入相对应的MATLAB函数名称。
进一步的技术方案是，所述边界条件包括：内边界施加钻井液压力Pm，外边界分别施加最大水平地应力σH和最小水平地应力σh，整个域上施加孔隙压力 Pp。
进一步的技术方案是，在划分有限元网格中由于井壁渐进破坏的方向总与最小水平地应力方向平行，为了提高计算精度，特别是最小水平地应力方向进行网格加密，远离井筒，网格逐渐变粗，以减少要计算的单元数量。
进一步的技术方案是，所述步骤S70中的计算公式为：
E＝(1-F)E0
Figure RE-GDA0003131971070000061
Figure RE-GDA0003131971070000062
式中：E、E0分别为岩石损伤前后的杨氏模量；φ、φ0分别为岩石损伤前后的孔隙度；K、K0分别为岩石损伤前后的渗透率；F为损伤变量。
进一步的技术方案是，所述NYZA的计算公式为：
Figure RE-GDA0003131971070000063
式中：NYZA为归一化屈服区面积；A1为井壁周围屈服区面积；A2为井眼初始面积。
本发明具有以下有益效果：本发明使用输入参数的不确定性，表示了岩石力学参数的非均质性；采用连续损伤理论，更改已破坏区域的力学参数，视为该区域等效挖去，模拟了井壁渐进破坏发展的四个阶段：裂纹萌生、扩展、崩落形成和稳定。
附图说明
图1为实施流程图；
图2为岩石力学实验结果图；
图3为物性参数和岩石力学参数最适合正态分布图；
图4为四分之一几何模型图；
图5为非均质地层参数赋值结果图；
图6为边界条件施加结果图；
图7为非均质地层有限元模型划分有限元网格图；
图8为井壁崩落区域图；
图9为NYZA数值图。
具体实施方式
下面结合实施例和附图对本发明做更进一步的说明。
如图1所示，本发明的一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法，包括以下步骤：
步骤S10、根据室内实验和测井资料，确定井壁坍塌渐进破坏过程模拟目标井地层的地应力、孔隙压力、井筒压力等基础参数；所述地层地应力包括最大水平地应力、最小水平地应力；
步骤S20、根据室内岩石力学实验，确定井壁坍塌渐进破坏过程模拟目标井地层岩石的岩石力学参数样本数量、岩石基础物性参数和岩石力学参数；所述岩石力学实验包括单轴压缩和三轴压缩岩石力学实验；所述岩石基础物性参数包括孔隙度、渗透率和Biot系数；所述岩石力学参数包括杨氏模量、泊松比、内摩擦角、内聚力；
步骤S30、根据步骤S20基本物性参数和岩石力学参数确定每个参数的最佳分布函数及对应的特征参数；
所述的参数分布函数包括均匀分布、三角分布、正态分布、对数正态分布、Beta分布、Gumbel分布、Weibull分布、Gamma分布等多种分布函数。具体分布见表1。
表1参数分布
Figure RE-GDA0003131971070000071
Figure RE-GDA0003131971070000081
其中，x是随机变量；a、b和c又分别称为最小值、最大值和众数；λλ是比例参数(scale)；kk是形状参数；σ是标准差；U是平均值；
步骤S40、建立井周应力分布的流-固-热耦合有限元数学模型；
对于流经饱和非等温变形多孔介质的可压缩流体，以位移、压力和温度为未知量时，控制方程可描述为(忽略体力)，具体过程如下：
步骤S41、建立控制方程：
Figure RE-GDA0003131971070000082
Figure RE-GDA0003131971070000083
Figure RE-GDA0003131971070000084
式中，G和λ是拉梅常数；k是多孔介质渗透率；u是流体的粘度；u、p和T 分别是多孔介质的位移、孔隙压力和温度；下标t表示时间的导数；φ是多孔介质的孔隙度；Kf、Km分别为流体和岩石的体积模量；IT＝[1,1,1,0,0,0]；D是弹性刚度矩阵；λT是多孔介质基质导热系数；ρscs是多孔介质基质热熔；ρfcf是流体热熔；βs是多孔介质基质热膨胀系数；βf是流体热膨胀系数；
步骤S41、使用伽辽金法有限元法逼近控制方程，可以得到控制方程(1)- (3)的有限元求解格式：
Figure RE-GDA0003131971070000091
其中，
M＝∫VBTDBdV------------------------------------(4)
Figure RE-GDA0003131971070000092
Figure RE-GDA0003131971070000093
Figure RE-GDA0003131971070000094
Figure RE-GDA0003131971070000095
Figure RE-GDA0003131971070000096
Figure RE-GDA0003131971070000097
Figure RE-GDA0003131971070000098
Figure RE-GDA0003131971070000099
式中，B为应变和位移相关的应变矩阵；上标T为矩阵转置；N有限元形函数；u、p和T分别为未知变量u、p和T的向量；ut、pt、Tt分别为未知变量u、 p和T的时间导数；fu、fp、fT分别是节点载荷向量、流体源汇向量和热源向量；
Figure RE-GDA00031319710700000910
B＝LNu----------------------------------------(14)
式中，L是微分算子；Nu是位移形函数。
步骤S50、根据井周应力分布流-固耦合有限元数学模型，建立非均质地层井壁渐进坍塌过程数值模拟的有限元模型，对非均质地层有限元模型赋值材料参数、施加边界条件和划分有限元网格，计算井周应力分布；
步骤S51、根据平面应变和轴对称条件，建立非均质地层井壁渐进坍塌过程有限元数值模拟的1/4几何模型，其井筒半径为R，整个模型尺寸为50R×50R 的几何模型，以排除边界效应的影响，并为了减少计算量。
步骤S52、对步骤S51建立的非均质地层有限元数值模拟几何模型求解域赋值非均质材料参数，设定求解域的岩石基础物性参数和岩石力学参数；
具体非均质赋值过程：
(1)选取合适每个参数的分布方式并编译MATLAB外部随机分布函数，其中随机分布函数以空间位置(X，Y)作为因变量；
(2)把MATLAB外部随机分布函数放置COMSOL非均质地层有限元程序同一文件夹；
(3)点击COMSOL Link MATLAB，把MATLAB与COMSOL链接起，即可以调用已编译好的MATLAB函数；
(4)在COMSOL中固体力学接口中，给内聚力、内摩擦角、杨氏模量、泊松比和Biot系数输入相对应的MATLAB函数名称；在COMSOL中流体力学接口中给孔隙度和渗透率输入相对应的MATLAB函数名称；
S53、对非均质地层有限元模型施加边界条件，其中，内边界(井筒)施加钻井液压力Pm，外边界分别施加最大水平地应力σH和最小水平地应力σh，整个域上施加孔隙压力Pp；
S54、对非均质地层有限元模型划分有限元网格，由于井壁渐进破坏的方向总与最小水平地应力方向平行，为了提高计算精度，特别是最小水平地应力方向进行网格加密，远离井筒，网格逐渐变粗，以减少要计算的单元数量；
S35、对非均质地层有限元模型网格模型进行求解，计算井周应力分布，从而得到最大主应力σ1和最小主应力σ3；
步骤S60、根据步骤S50得到的应力计算损伤变量F；
步骤S61、确定地层岩石破坏时的应力状态函数：最大拉应力状态函数和 Mohr-Coulomb应力状态函数；
F1＝-σ3-St--------------------------------------(15)
Figure RE-GDA0003131971070000111
式中，C为岩石内聚力；
Figure RE-GDA0003131971070000112
为岩石内摩擦角；St是岩石抗拉强度；F1为最大拉应力状态函数；F2为Mohr-Coulomb状态函数；
步骤S62、根据最大拉应力状态函数F1和Mohr-Coulomb应力状态函数F2，计算地层岩石的损伤变量F；
Figure RE-GDA0003131971070000113
式中，F为损伤变量；当F＝0对应无损伤单元状态，F＝1对应完全损伤状态，即岩石发生断裂或破坏损伤；F1≥0表示地层岩石发生拉伸破坏损伤；F2≥0表示地层岩石发生剪切破坏损伤；F1'>0和F2'>0分别为两种损伤后的继续加载状态，可引起损伤变量的增加；当F1'<0和F2'<0表示卸载状态，不产生新的损伤，损伤变量保持上一个加载步(时间步)的数值；
步骤S70、根据步骤S60的损伤变量F确定破坏区域，更新破坏区域地层参数杨氏模量、孔隙度、渗透率，使得破坏区域内的单元丧失承载应力能力；
步骤S71、根据各向同性弹性损伤理论，单元的弹性模量会随着损伤过程逐渐退化，损伤后的杨氏模量可表示为：
E＝(1-F)E0------------------------------------(18)
式中：E、E0分别为岩石损伤前后的杨氏模量；由于此处采用的是弹性损伤理论，当F＝1时岩石发生破坏性损伤，为了消除数值计算带来的误差问题，假定F＝1的单元体杨氏模量由万分之一的初始弹性模量代替，即更新F＝1的单元体杨氏模量为0.0001E0；
步骤S72、受应力集中破坏的岩石受钻井液循环而带离原有位置，破坏区域应等效为完全由钻井液支撑，即破坏区域被视为新的内边界。
因此，有必要加入孔隙度和渗透率与岩石单元损伤的函数，其物理场参数如式(19)和(20)所示：
Figure RE-GDA0003131971070000121
Figure RE-GDA0003131971070000122
对于破坏区域，渗透率应无限大，对于数值模拟无法输入无限大这个变量，对于渗透率大于4达西对渐进破坏结果无明显影响；
因此后续数值模拟时，发生井壁破坏时，渗透率至少大于4达西；
步骤S80、重复步骤S50至步骤S70进行迭代，直到损伤变量F≤0，此时，无新的损伤区出现，井眼最终趋于稳定，实现非均质地层井壁渐进坍塌过程的数值模拟，并绘制井壁失稳区域图，确定每次迭代步下的NYZA值；
NYZA的计算：在本文中使用归一化屈服区面积(NYZA)法判断合适的钻井液压力。NYZA的计算方法是将井壁周围屈服区面积除以初始井面积。
Figure RE-GDA0003131971070000123
其中NYZA为归一化屈服区面积，A1为井壁周围屈服区面积，A2为井眼初始面积。
Hawkes和McLellan指出，当NYZA小于1时，井筒在弹性条件下是稳定的。因此，为了确定合理的井筒压力，本文将NYZA＝1作为钻井液压力的阈值。
实施例
步骤S10、根据室内实验和测井资料，确定井壁坍塌渐进破坏过程模拟目标井地层的地应力、孔隙压力、井筒压力等基础参数；
步骤S20、根据室内岩石力学实验，确定井壁坍塌渐进破坏过程模拟目标井地层岩石的岩石力学参数样本数量27个、基本物性参数和岩石力学参数；
所述的岩石力学实验包括单轴压缩和三轴压缩岩石力学实验，实验结果如图2。不难看出，弹性模量和抗压强度均随围压的增加而增加，但泊松比并为随围压增加而不发生明显的变化；误差棒相对较长，说明实验结果波动性较大，也表明火山岩的非均质性强；
步骤S30、根据步骤S20基本物性参数和岩石力学参数确定其分布规律和合适的分布方式。物性参数和岩石力学参数最适合正态分布，正态分布如图3 所示，表2为正态分布特征参数。
表2正态分布特征参数
Figure RE-GDA0003131971070000131
步骤S40、建立井周应力分布的流-固-热耦合有限元数学模型。对于流经饱和非等温变形多孔介质的可压缩流体，以位移、压力和温度为未知量时，控制方程可描述为(忽略体力)；
步骤S50、根据步骤S40的井周应力分布流-固耦合有限元数学模型，建立非均质地层井壁渐进坍塌过程数值模拟的有限元几何模型，对非均质地层有限元模型赋值材料参数、施加边界条件和划分有限元网格，计算井周应力分布；
步骤S51、根据平面应变和轴对称条件，如图4建立井壁坍塌渐进破坏过程有限元数值模拟的1/4几何模型，其井筒半径为R，整个模型尺寸为50R×50R 的几何模型，以排除边界效应的影响，并为了减少计算量。
步骤S52、对建立的非均质地层有限元数值模拟几何模型求解域赋值非均质材料参数，设定求解域的基础物性参数和岩石力学参数，包括孔隙度、渗透率和Biot系数、杨氏模量、泊松比、内摩擦角、内聚力等基础参数。具体非均质赋值过程：(1)选取正态分布方式并编译MATLAB外部随机分布函数，其中随机分布函数以空间位置(X，Y)作为因变量。(2)把MATLAB外部随机分布函数放置COMSOL非均质地层有限元程序同一文件夹；(3)点击COMSOLLink MATLAB，把MATLAB与COMSOL链接起，即可以调用已编译好的 MATLAB函数；(4)在COMSOL中固体力学接口中，给内聚力、内摩擦角、杨氏模量、泊松比和Biot系数输入相对应的MATLAB函数名称；在COMSOL 中流体力学接口中给孔隙度和渗透率输入相对应的MATLAB函数名称。非均质地层参数赋值结果如图5所示。
步骤S53、对非均质地层有限元模型施加边界条件，其中，内边界(井筒) 施加钻井液压力Pm，外边界分别施加最大水平地应力σH和最小水平地应力σh，整个域上施加孔隙压力Pp。外边界和整个域上施加地层温度T0，内边界(井筒) 施加钻井液温度Tm。
边界条件施加结果如图6所示。所施加的边界条件及基础参数如表3所示。
表3边界条件参数以及基础参数
No. 参数 符号 值 单位
1 最大水平地应力 σ<sub>H</sub> 164.4 MPa
2 最小水平地应力 σ<sub>h</sub> 137.1 MPa
3 钻井液压力 P<sub>m</sub> 98 MPa
4 孔隙压力 P<sub>p</sub> 102 MPa
5 地层温度 T<sub>0</sub> 350 K
6 钻井液温度 T<sub>m</sub> 300 K
步骤S54、对非均质地层有限元模型划分有限元网格如图7，由于井壁渐进破坏的方向总与最小水平地应力方向平行，为了提高计算精度，特别是最小水平地应力方向进行网格加密，远离井筒，网格逐渐变粗，以减少要计算的单元数量。
步骤S55、对非均质地层有限元模型网格模型进行求解，计算井周应力分布，从而得到最大主应力σ1、中间主应力σ2和最小主应力σ3。
步骤S60、根据步骤S50得到的应力计算损伤变量F；
步骤S70、根据步骤S60的损伤变量F确定破坏区域，更新破坏区域地层参数杨氏模量、孔隙度、渗透率，使得破坏区域内的单元丧失承载应力能力；
步骤S80、重复步骤S50至步骤S70进行迭代，直到损伤变量F≤0，此时，无新的损伤区出现，井眼最终趋于稳定，实现非均质地层井壁渐进坍塌过程的数值模拟，并绘制井壁失稳区域图，确定每次迭代步下的NYZA值。
图8和图9分别为每个迭代步下得到的井壁崩落区域和对应的NYZA数值。不难看出：破坏区域主要发生在最小水平应力(σh)方向；因为在这个方向上有较大的压应力集中，井筒压力不足以支撑围岩发生剪切破坏；破坏区域不连续，因为地层的岩石力学参数随机分布；图8在第2步进行破坏区域预测时，只有少量岩石破坏，井壁保持较好的圆形；在第2步至第4步之间，由于地层渐进破坏，井眼快速扩大，在第4步之后，出砂区域沿着之前的破坏区域向远处缓慢扩展。根据NYZA公式计算井周屈服区域来判断井壁是否稳定，其NYZA是0.5，其生产压差满足井壁稳定；统计破坏区域(NYZA)，随每次迭代步的增长，每一步的NYZA区域如图9所示。第1步的NYZA是0.20，即破坏区域十分小，第4步的NYZA是0.31，可以从图9中看出前4步NYZA的增长率远高于4步到27步，后23步的NYZA增长率并不满足线性增长，这也表示井周渐进出砂的迭代导致井周应力集中的释放，破坏区域在逐步减小，即井眼趋于稳定。
以上所述，并非对本发明作任何形式上的限制，虽然本发明已通过上述实施例揭示，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容作出些变动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims (10)

1.一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法，其特征在于，包括以下步骤：
步骤S10、根据室内实验和测井资料确定目标井的地层地应力、孔隙压力、井筒压力；
步骤S20、根据室内岩石力学实验确定目标井地层岩石的岩石力学参数样本数量、岩石基础物性参数和岩石力学参数；
步骤S30、根据岩石基础物性参数和岩石力学参数确定每个参数的最佳分布函数及对应的特征参数；
步骤S40、建立井周应力分布的流-固-热耦合有限元数学模型；
步骤S50、根据井周应力分布流-固耦合有限元数学模型计算井周应力分布；
步骤S60、根据井周应力分布计算损伤变量F；
Figure FDA0003093848490000011
F1＝-σ3-St
Figure FDA0003093848490000012
式中：C为岩石内聚力；
Figure FDA0003093848490000013
为岩石内摩擦角；St是岩石抗拉强度；F1为最大拉应力状态函数；F2为Mohr-Coulomb状态函数；F为损伤变量；σ1为最大主应力；σ3为最小主应力；
步骤S70、损伤变量F确定破坏区域，更新破坏区域的杨氏模量、孔隙度、渗透率；
步骤S80、重复步骤S50至步骤S70进行迭代，直到损伤变量F≤0，此时，无新的损伤区出现，井眼最终趋于稳定，实现非均质地层井壁渐进坍塌过程的数值模拟，并绘制井壁失稳区域图，确定每次迭代步下的NYZA值。
2.根据权利要求1所述的一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法，其特征在于，所述地层地应力包括最大水平地应力、最小水平地应力。
3.根据权利要求1所述的一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法，其特征在于，所述岩石力学实验包括单轴压缩和三轴压缩岩石力学实验；所述岩石基础物性参数包括孔隙度、渗透率和Biot系数；所述岩石力学参数包括杨氏模量、泊松比、内摩擦角、内聚力。
4.根据权利要求1所述的一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法，其特征在于，所述步骤S40的具体过程为：
步骤S41、建立控制方程：
Figure FDA0003093848490000021
Figure FDA0003093848490000022
Figure FDA0003093848490000023
式中，G和λ是拉梅常数；k是多孔介质渗透率；u是流体的粘度；u、p和T分别是多孔介质的位移、孔隙压力和温度；下标t表示时间的导数；φ是多孔介质的孔隙度；Kf、Km分别为流体和岩石的体积模量；IT＝[1,1,1,0,0,0]；D是弹性刚度矩阵；λT是多孔介质基质导热系数；ρscs是多孔介质基质热熔；ρfcf是流体热熔；βs是多孔介质基质热膨胀系数；βf是流体热膨胀系数；
步骤S41、使用伽辽金法有限元法逼近控制方程，可以得到控制方程的有限元求解格式：
Figure FDA0003093848490000024
其中，
M＝∫VBTDBdV
Figure FDA0003093848490000031
Figure FDA0003093848490000032
Figure FDA0003093848490000033
Figure FDA0003093848490000034
Figure FDA0003093848490000035
Figure FDA0003093848490000036
HTT＝λT∫V(▽N)(▽N)TdV
Figure FDA0003093848490000037
Figure FDA0003093848490000038
B＝LNu
式中，B为应变和位移相关的应变矩阵；上标T为矩阵转置；N有限元形函数；u、p和T分别为未知变量u、p和T的向量；ut、pt、Tt分别为未知变量u、p和T的时间导数；fu、fp、fT分别是节点载荷向量、流体源汇向量和热源向量；L是微分算子；Nu是位移形函数。
5.根据权利要求4所述的一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法，其特征在于，所述步骤S50的具体过程为：根据井周应力分布流-固耦合有限元数学模型建立非均质地层井壁渐进坍塌过程数值模拟的有限元模型，分别对非均质地层有限元模型赋值非均质材料参数、施加边界条件和划分有限元网格，最后计算井周应力分布，得到最大主应力σ1和最小主应力σ3。
6.根据权利要求5所述的一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法，其特征在于，根据井周应力分布流-固耦合有限元数学模型建立非均质地层井壁渐进坍塌过程数值模拟的有限元模型包括：根据平面应变和轴对称条件、井周应力分布流-固耦合有限元数学模型，建立非均质地层井壁渐进坍塌过程数值模拟的有限元模型，其井筒半径为R，整个模型为50R×50R的几何模型，以排除边界效应的影响，并为了减少计算量。
7.根据权利要求5所述的一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法，其特征在于，所述非均质地层井壁渐进坍塌过程数值模拟有限元模型的具体非均质赋值过程为：
首先选取合适每个参数的分布方式并编译MATLAB外部随机分布函数，其中随机分布函数以空间位置作为因变量；
然后把MATLAB外部随机分布函数放置COMSOL非均质地层有限元程序同一文件夹；
再点击COMSOL Link MATLAB，把MATLAB与COMSOL链接起，即可以调用已编译好的MATLAB函数；
最后在COMSOL中固体力学接口中，给内聚力、内摩擦角、杨氏模量、泊松比和Biot系数输入相对应的MATLAB函数名称；在COMSOL中流体力学接口中给孔隙度和渗透率输入相对应的MATLAB函数名称。
8.根据权利要求5所述的一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法，其特征在于，所述边界条件包括：内边界施加钻井液压力Pm，外边界分别施加最大水平地应力σH和最小水平地应力σh，整个域上施加孔隙压力Pp。
9.根据权利要求5所述的一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法，其特征在于，所述步骤S70中的计算公式为：
E＝(1-F)E0
Figure FDA0003093848490000051
Figure FDA0003093848490000052
式中：E、E0分别为岩石损伤前后的杨氏模量；φ、φ0分别为岩石损伤前后的孔隙度；K、K0分别为岩石损伤前后的渗透率；F为损伤变量。
10.根据权利要求9所述的一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法，其特征在于，所述NYZA的计算公式为：
Figure FDA0003093848490000053
式中：NYZA为归一化屈服区面积；A1为井壁周围屈服区面积；A2为井眼初始面积。
CN202110604358.5A 2021-05-31 2021-05-31 一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法 Active CN113324840B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110604358.5A CN113324840B (zh) 2021-05-31 2021-05-31 一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110604358.5A CN113324840B (zh) 2021-05-31 2021-05-31 一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法

Publications (2)

Publication Number Publication Date
CN113324840A true CN113324840A (zh) 2021-08-31
CN113324840B CN113324840B (zh) 2022-04-22

Family

ID=77422771

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110604358.5A Active CN113324840B (zh) 2021-05-31 2021-05-31 一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法

Country Status (1)

Country Link
CN (1) CN113324840B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114707367A (zh) * 2022-06-06 2022-07-05 浙江陆特能源科技股份有限公司 中深层地埋管换热器传热分析模型
CN117268617A (zh) * 2023-09-13 2023-12-22 中国科学院武汉岩土力学研究所 一种包含三层介质模型的应力张量确定方法

Citations (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105401939A (zh) * 2015-11-30 2016-03-16 中国石油大学(北京) 一种多因素耦合作用下的煤层井壁稳定性分析方法
CN106682757A (zh) * 2015-12-31 2017-05-17 中国石油天然气股份有限公司 一种井壁坍塌压力确定方法和装置
CN107220467A (zh) * 2017-07-07 2017-09-29 中国水利水电科学研究院 蓄水期库岸岩质边坡变形的预测方法
CN107526892A (zh) * 2017-08-30 2017-12-29 广州海洋地质调查局 一种海洋天然气水合物试采储层的稳定性评估方法
CN108830020A (zh) * 2018-07-12 2018-11-16 西南石油大学 一种基于热流固耦合理论的模拟海上油田微压裂增注裂缝扩展的方法
CN109356567A (zh) * 2018-05-04 2019-02-19 中国石油集团海洋工程有限公司 深水浅部地层井壁稳定性预测方法
US20190163855A1 (en) * 2015-05-20 2019-05-30 Saudi Arabian Oil Company Parallel Solution for Fully-Coupled Fully-Implicit Wellbore Modeling in Reservoir Simulation
CN109858147A (zh) * 2019-01-30 2019-06-07 西南石油大学 一种基于可靠度理论的井壁失稳风险定量评价方法
US10378941B2 (en) * 2013-04-30 2019-08-13 Iphase Limited Method and apparatus for monitoring the flow of mixtures of fluid in a pipe
US20200149393A1 (en) * 2017-04-03 2020-05-14 Repsol, S.A. Method of estimating the region of damage due to collapse in the wall of a borehole during the drilling operation
CN111366464A (zh) * 2020-04-19 2020-07-03 长江大学 一种破碎性地层岩石力学参数确定方法
CN112036096A (zh) * 2020-09-07 2020-12-04 西南石油大学 裂缝性地层井壁强化效果评价的流固耦合数值模拟方法

Patent Citations (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10378941B2 (en) * 2013-04-30 2019-08-13 Iphase Limited Method and apparatus for monitoring the flow of mixtures of fluid in a pipe
US20190163855A1 (en) * 2015-05-20 2019-05-30 Saudi Arabian Oil Company Parallel Solution for Fully-Coupled Fully-Implicit Wellbore Modeling in Reservoir Simulation
CN105401939A (zh) * 2015-11-30 2016-03-16 中国石油大学(北京) 一种多因素耦合作用下的煤层井壁稳定性分析方法
CN106682757A (zh) * 2015-12-31 2017-05-17 中国石油天然气股份有限公司 一种井壁坍塌压力确定方法和装置
US20200149393A1 (en) * 2017-04-03 2020-05-14 Repsol, S.A. Method of estimating the region of damage due to collapse in the wall of a borehole during the drilling operation
CN107220467A (zh) * 2017-07-07 2017-09-29 中国水利水电科学研究院 蓄水期库岸岩质边坡变形的预测方法
CN107526892A (zh) * 2017-08-30 2017-12-29 广州海洋地质调查局 一种海洋天然气水合物试采储层的稳定性评估方法
CN109356567A (zh) * 2018-05-04 2019-02-19 中国石油集团海洋工程有限公司 深水浅部地层井壁稳定性预测方法
CN108830020A (zh) * 2018-07-12 2018-11-16 西南石油大学 一种基于热流固耦合理论的模拟海上油田微压裂增注裂缝扩展的方法
CN109858147A (zh) * 2019-01-30 2019-06-07 西南石油大学 一种基于可靠度理论的井壁失稳风险定量评价方法
CN111366464A (zh) * 2020-04-19 2020-07-03 长江大学 一种破碎性地层岩石力学参数确定方法
CN112036096A (zh) * 2020-09-07 2020-12-04 西南石油大学 裂缝性地层井壁强化效果评价的流固耦合数值模拟方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
MA TIANSHOU: "Wellbore stability analysis for arbitrary inclined well in anisotropic formations", 《IOP CONFERENCE SERIES EARTH AND ENVIRONMENTAL SCIENCE》 *
冯雨实: "煤层气水平井井周围岩热流固耦合数值分析", 《煤矿安全》 *
李玉梅: "欠平衡钻井煤层井壁稳定有限元数值计算研究", 《系统仿真学报》 *
贾善坡: "基于热-流-固耦合模型的石油钻井施工过程数值分析", 《岩土力学》 *
陈颖杰: "井壁失稳风险的可靠度理论评价方法", 《钻井工程》 *
雷家蔚: "温压耦合影响下科学钻探深井结晶岩层井壁稳定性分析", 《地质与勘探》 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114707367A (zh) * 2022-06-06 2022-07-05 浙江陆特能源科技股份有限公司 中深层地埋管换热器传热分析模型
CN117268617A (zh) * 2023-09-13 2023-12-22 中国科学院武汉岩土力学研究所 一种包含三层介质模型的应力张量确定方法
CN117268617B (zh) * 2023-09-13 2024-05-28 中国科学院武汉岩土力学研究所 一种包含三层介质模型的应力张量确定方法

Also Published As

Publication number Publication date
CN113324840B (zh) 2022-04-22

Similar Documents

Publication Publication Date Title
CN113324840B (zh) 一种非均质地层井壁渐进坍塌过程的流-固-热耦合模拟方法
Yan et al. A two-dimensional coupled hydro-mechanical finite-discrete model considering porous media flow for simulating hydraulic fracturing
Papamichos et al. Volumetric sand production model and experiment
Tolooiyan et al. Modelling the cone penetration test in sand using cavity expansion and arbitrary Lagrangian Eulerian finite element methods
Pucker et al. Numerical simulation of the installation process of full displacement piles
Chen et al. Wellbore stability analysis using strain hardening and/or softening plasticity models
Xu et al. A numerical meso-scale elasto-plastic damage model for modeling the deformation and fracturing of sandstone under cyclic loading
Pattillo et al. Analysis of horizontal casing integrity in the valhall field
Wang et al. Numerical and experimental studies of pressure-controlled cavity expansion in completely decomposed granite soils of Hong Kong
Suzuki et al. Analysis of CPT end resistance at variable penetration rates using the spherical cavity expansion method in normally consolidated soils
Nishanthan et al. Shielding effect in pile groups adjacent to deep unbraced and braced excavations
Yuan et al. HPHT gas well cementing complications and its effect on casing collapse resistance
Vakili et al. A practical approach to constitutive models for the analysis of geotechnical problems
Arson et al. Numerical study of a thermo-hydro-mechanical damage model for unsaturated porous media
Cardoso Bernardes et al. Coupling hardening soil model and Ménard pressuremeter tests to predict pile behavior
Das et al. Analysis of deformational behavior of circular underground opening in soft ground using three-dimensional physical model
Shahin et al. Numerical simulation of localized compaction creep in heterogeneous porous rock
Pyrah et al. Effects of test procedure on constant rate of strain pressuremeter tests in clay
CN113343336B (zh) 一种井壁坍塌渐进破坏过程的数值模拟方法
Hernández-Hernández et al. Experimental and numerical analysis of triaxial compression test for a clay soil
Sarvaramini Quantifying variability in hydraulic fracture geometry using stage-by-stage ISIP analysis in geomechanically heterogeneous rocks
Khatibi et al. A method to find optimum mud weight in zones with no safe mud weight windows
Choi et al. Numerical analysis of laterally loaded piles affected by bedrock depth
Carvajal-Cardenas et al. Physical modeling of desiccated slopes in fine soil using a geotechnical centrifuge
Abdlrahem et al. Numerical Modelling of Axial Behaviour of Single and Grouped Hollow Bar Micropiles in Cohesionless Soils

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