本发明属于石油与天然气钻井工程领域，具体涉及一种考虑岩屑床厚度的钻柱摩阻计算方法。
背景技术：
1、大位移及大斜度井钻井技术因其增大了储层与井筒的接触面积，提高了钻井效率,广泛应用于油气资源勘探和开发工程。然而,在定向井或水平井钻井过程中，由于井眼轨迹的复杂性，钻井液携带能力不足，井眼净化问题突出。岩屑极易在环空底部积聚成岩屑床,增加了钻柱与井壁之间的摩擦系数和阻尼，使钻柱与井壁之间的摩阻急剧增加。此外，当井斜角过大时，岩屑由于重力的影响将沉积在井壁下放，并且钻柱在钻压下将产生较大的形变，增大了钻柱与井壁的接触面积，井壁对钻柱的磨损也随之增大。从而导致憋泵、卡钻等事故,严重影响了钻井施工安全和钻井时效。
2、因此,针对存在岩屑床的钻柱摩阻问题，有必要对钻柱的摩阻进行研究。
技术实现思路
1、本发明所为了解决背景技术中存在的技术问题，目的在于提供了一种考虑岩屑床厚度的钻柱摩阻计算方法，在传统钻柱摩阻计算方法的基础上，引入了钻柱系统非线性有限元模型及岩屑运移模型，考虑了岩屑床对钻柱的影响以及多种外力边界条件，在迭代过程中不断修正计算结果，可以有效计算钻柱在含岩屑床的井眼中的摩阻特性，计算结果更接近工程实际。避免了卡钻风险，为现场钻井安全提供了理论支持。
2、为了解决技术问题，本发明的技术方案是：
3、一种考虑岩屑床厚度的钻柱摩阻计算方法，所述方法包括：
4、步骤s1：获取钻具组合参数及材料特性参数，建立全局、井眼及单元坐标系，计算全局坐标系下钻柱的井眼轨迹，根据钻具组合和材料特性参数，将钻柱系统离散为n个梁单元，从井口到井底依次对梁单元进行标号；
5、步骤s2：根据水力、泥浆、泵压及钻井参数，构建岩屑运移模型，计算岩屑床的厚度及环空间隙，根据岩屑颗粒碰撞方程计算岩屑间接触力，并拟合钻柱与岩屑之间的摩擦系数；
6、步骤s3：根据步骤s2得到的环空间隙，计算井壁与钻柱的间隙，从而得到含岩屑床的井壁边界条件，输入井口和井底边界条件；
7、步骤s4：基于梁单元的动能、势能和外力功，利用哈密顿原理建立单元坐标系下的动力学方程；
8、步骤s5：结合步骤s3得到的边界条件，迭代求解所述动力学方程，计算节点的位移；
9、步骤s6：根据步骤s5得到的节点位移，判断钻柱是否与井壁接触；
10、步骤s7：若接触，计算接触时的节点接触力，根据接触力计算节点摩阻并保存数据；
11、步骤s8：若没有接触，重复步骤s2-s7，计算下一个钻柱节点(i+1)的摩阻值，依次迭代循环，直至计算出第n个钻柱节点摩阻；最终，全井段钻柱系统所有节点摩阻值计算完毕，结束计算；汇总得到含岩屑床的全井段钻柱摩阻计算结果和分布规律。
12、进一步，所述步骤s1的全局坐标系o-ars描述井眼轨迹的变化，原点o位于井口处，a轴指向正北方向，r轴指向正东方向，s轴垂直向下指向井眼方向；井眼坐标系o1-xyz描述钻柱的运动，原点o1位于井眼中心，x轴指向井眼切线方向，y轴指向井眼轴线法线方向，z由右手定则确定；单元坐标系o-xyz描述钻柱的变形，原点o位于钻柱中心，x轴指向钻柱轴线的切线方向，y轴和z轴都位于钻柱截面内。
13、进一步，所述步骤s1的井眼轨迹数据包括：井斜角、方位角以及测深；所述步骤s1的轨迹计算采用曲率半径法，具体包括：
14、计算相邻两测点间的坐标增量δx、δy、δz，其中：
15、(1)若
16、
17、(2)若
18、
19、(3)若
20、
21、(4)若
22、
23、上式中，δs＝s2-s1表示相邻两测点间的曲线长度，m；α2、α1表示相邻两测点间的井斜角，°；表示相邻两测点间的方位角，°；
24、计算井眼轨迹中任意相邻两点的三维坐标：
25、
26、上式中，x1、x2表示节点1和节点2沿x轴的坐标，m；y1、y2表示节点1和节点2沿y轴的坐标，m；z1、z2表示节点1和节点2沿z轴的坐标，m；
27、按照上述方法依次循环计算i＝1,2,…,n时节点的三维坐标；进而得到全井段井眼轨迹。
28、进一步，所述步骤s2构建的岩屑运移模型由岩屑颗粒重力、浮力、接触力(岩屑-岩屑、岩屑-钻杆、岩屑-岩壁)和相互作用力(拖曳力、剪切升力、旋转升力)决定，表示为：
29、
30、上式中，mp为岩屑p的质量，kg；为岩屑颗粒p和q之间的相互接触力，n；fd为流体拖曳力，n；fs为剪切升力，n；fm为旋转升力，n；fp为流体压力的梯度力，n；
31、岩屑床的厚度h可以表示为：
32、
33、上式中，h为岩屑床厚度，m；d为井眼直径，m；ν为环空平均流速，m/s；n为钻柱转速，r/min；u为偏心度；ρf为钻井液密度，kg/m3；μe为钻井液有效粘度，mpa·s；a1-a15为经验系数，可通过多元非线性数学回归的方法获得；φ为岩屑的形状系数；ds为岩屑的当量直径，m；cangθ与cang90分别表示井斜角为θ和90度时的修正系数；
34、环空间隙表示为：
35、
36、上式中，δ为环空间隙，m。
37、进一步，所述步骤s2中，钻柱与岩屑之间的摩擦系数拟合公式表示为：
38、μ＝tan(ψ)/(1+σnc)
39、上式中，μ为钻柱与岩屑之间的摩擦系数；ψ为摩擦角，rad；σn为正应力或侧压力，n；c为内聚力，n。
40、进一步，所述步骤s3的井口和井底边界条件表示为：
41、
42、上式中，u为第一个梁单元的轴向拉伸长度，m；l0为第一个梁单元的初始长度，m；wn和tn分别为钻头在轴向和扭转方向分别受到钻压和扭矩；e为材料的弹性模量，pa；a为钻柱截面面积，m2；v和w分别表示钻柱节点延y和z方向上的位移，m；θy和θz分别表示梁单元绕y轴和z轴的扭转角度，rad；下标1和b分别表示第一个节点和最后一个节点。
43、进一步，所述步骤s4的动力学方程表示为：
44、
45、上式中，m、c和k分别表示整体坐标系下钻柱系统的质量矩阵、阻尼矩阵和刚度矩阵；q和f分别表示整体坐标系下钻柱系统的加速度向量、速度向量、位移向量和外力向量。
46、进一步，所述外力向量包括梁单元重力、不平衡力、钻柱与井壁的接触力以及钻柱受到的钻井液浮力，具体表示为：
47、(1)梁单元自身重力和钻井液浮力的等效节点力表示为：
48、
49、上式中，le表示梁单元长度，mm；单位长度的梁单元在x轴和y轴的重力分量分别表示为qx和qy；ρ为梁单元材料密度，kg/m3；
50、(2)梁单元由于质量偏心引起的不平衡力表示为：
51、
52、上式中，fy和fz分别表示沿着梁单元y轴和z轴的离心力分量；
53、(3)当钻柱与井壁接触时，钻柱与井壁的接触力向量表示为：
54、fc＝[0 fy1 fz1 tf1 0 0 0 fz1 fz2 tf1 0 0]t
55、上式中，fy和fz分别表示沿着梁单元y轴和z轴的摩擦力分量；tf表示钻柱与井壁的摩擦力矩。
56、进一步，所述步骤s5的节点的横向位移表示为：
57、
58、上式中，和分别表示第j+1次和第j次迭代计算的节点位移；αr、δr、ηr、βr、εr和γr为独立的常系数；ρr为谱半径；δt为计算时间步长，s；为残余力，n。
59、进一步，计算步骤s7钻柱节点与岩屑的接触力，表示为：
60、fn＝kc(r-δ)
61、上式中，kc表示岩石的接触刚度，n/m，迭代计算节点位移和接触力，直至满足容许误差。
62、若节点位移q小于环空间隙δ，则钻柱与岩屑不接触，摩阻值为零；
63、若节点位移q大于等于环空间隙δ，则钻柱与岩屑接触，根据滑动摩擦力公式f＝μfn计算摩阻。
64、一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现上述中任一项所述的一种考虑岩屑床厚度的钻柱摩阻计算方法。
65、一种计算机可读存储介质，所述计算机可读存储介质上存储有计算机程序，该程序被处理器执行时实现上述中任一项所述的一种考虑岩屑床厚度的钻柱摩阻计算方法。
66、与现有技术相比，本发明的优点在于：
67、该方法首先根据实际井眼轨迹和钻具组合参数，将井眼中的钻柱空间离散为n个梁单元。考虑水力、泥浆、钻井等参数建立岩屑运移模型，计算岩屑床厚度及岩屑间接触力，拟合出钻柱与岩屑之间的摩擦系数。利用哈密顿原理建立钻柱系统动力学方程，结合井口、井底及考虑岩屑床的井壁边界条件，迭代计算得到满足精度的解，接着判断钻柱与井壁是否接触，计算接触时的节点摩阻值。依次计算钻柱系统各个节点直至结束。本发明考虑了岩屑床厚度对钻柱钻进时摩阻的影响，引入了多种接触模型，计算结果更接近工程实际。避免了卡钻风险，为现场钻井安全提供了理论支持。
1.一种考虑岩屑床厚度的钻柱摩阻计算方法，其特征在于，所述方法包括：
2.根据权利要求1所述的一种考虑岩屑床厚度的钻柱摩阻计算方法，其特征在于，所述步骤s1的全局坐标系o-ars描述井眼轨迹的变化，原点o位于井口处，a轴指向正北方向，r轴指向正东方向，s轴垂直向下指向井眼方向；井眼坐标系o1-xyz描述钻柱的运动，原点o1位于井眼中心，x轴指向井眼切线方向，y轴指向井眼轴线法线方向，z由右手定则确定；单元坐标系o-xyz描述钻柱的变形，原点o位于钻柱中心，x轴指向钻柱轴线的切线方向，y轴和z轴都位于钻柱截面内。
3.根据权利要求1所述的一种考虑岩屑床厚度的钻柱摩阻计算方法，其特征在于，所述步骤s1的井眼轨迹数据包括：井斜角、方位角以及测深；所述步骤s1的轨迹计算采用曲率半径法，具体包括：
4.根据权利要求1所述的一种考虑岩屑床厚度的钻柱摩阻计算方法，其特征在于，所述步骤s2构建的岩屑运移模型由岩屑颗粒重力、浮力、接触力(岩屑-岩屑、岩屑-钻杆、岩屑-岩壁)和相互作用力(拖曳力、剪切升力、旋转升力)决定，表示为：
5.根据权利要求1所述的一种考虑岩屑床厚度的钻柱摩阻计算方法，其特征在于，所述步骤s2中，钻柱与岩屑之间的摩擦系数拟合公式表示为：
6.根据权利要求1所述的一种考虑岩屑床厚度的钻柱摩阻计算方法，其特征在于，所述步骤s3的井口和井底边界条件表示为：
7.根据权利要求1所述的一种考虑岩屑床厚度的钻柱摩阻计算方法，其特征在于，所述步骤s4的动力学方程表示为：
8.根据权利要求7所述的一种考虑岩屑床厚度的钻柱摩阻计算方法，其特征在于，所述外力向量包括梁单元重力、不平衡力、钻柱与井壁的接触力以及钻柱受到的钻井液浮力，具体表示为：
9.根据权利要求1所述的一种考虑岩屑床厚度的钻柱摩阻计算方法，其特征在于，所述步骤s5的节点的横向位移表示为：
10.根据权利要求1所述的一种考虑岩屑床厚度的钻柱摩阻计算方法，其特征在于，计算步骤s7钻柱节点与岩屑的接触力，表示为：
