CN111291499A - 一种基于多体动力学的瓦斯抽采钻机建模方法 - Google Patents

Abstract

本发明公开了一种基于多体动力学的瓦斯抽采钻机建模方法，提出了一种新的瓦斯抽采钻机的建模方法，属于钻机建模技术领域，包括：第一，将钻头、孔壁和钻杆连接处视为刚体，采用牛顿‑欧拉法对其进行建模。第二，将每根钻杆均视为柔性体，用第二类拉格朗日方程对其进行建模。第三，考虑振动摩擦等因素对系统的影响，利用浮动节点坐标法建立松软煤中短水平段钻机‑钻杆‑孔壁系统的刚柔耦合多体动力学模型。第四，采用有限元方法对模型进行离散化，并采用虚拟样机技术对数学模型进行求解。该模型可以模拟钻机系统在钻进过程中的一系列动态响应，为后续的多因素耦合建模提供理论支持。

Description

一种基于多体动力学的瓦斯抽采钻机建模方法

技术领域

本发明涉及钻机装置技术领域，尤其涉及一种基于多体动力学的瓦斯抽采钻机建模方法。

背景技术

煤炭资源是世界主要能源之一。近年来，数次煤炭开采工程所引发的瓦斯爆炸灾难引发了广泛关注。统计发现，在过去的10年中，全国共发生煤矿瓦斯事故389起，死亡3307人。瓦斯事故发生频率虽低于运输事故，但死亡率却是其2倍之多。有效的抽采煤层瓦斯可以将瓦斯含量和压力降低到安全水平，从根本上消除瓦斯灾害。

近年来，随着煤层勘探开发规模的扩大，瓦斯抽采的成孔率问题受到了国内外学者的持续关注。但是，由于钻孔失效是由多个因素相互影响、相互促进的复杂的系统性问题，涉及到摩擦学、振动学、运动学、材料科学、化学、传热学等一系列学科，存在许多不确定性和随机性，因此，关于钻井方面的研究工作还需要进一步完善。

发明内容

针对上述缺陷或不足，本发明的目的在于提供一种基于多体动力学的瓦斯抽采钻机建模方法。

为达到以上目的，本发明的技术方案为：

一种基于多体动力学的瓦斯抽采钻机建模方法，包括以下步骤：

1)、将钻机、钻杆连接处、以及钻具均设为刚性体，并使用牛顿-欧拉法分别对钻机、钻杆连接处、以及钻具进行建模，建立整体坐标系，利用质量、惯性张量、广义力、以及广义力描述刚性体，得到刚性模型；

2)、将钻杆设为柔性体，建立参考坐标系和弹性坐标，并使用第二类拉格朗日方程进行建模，得到柔性模型；

3)、利用浮动节点坐标法对刚性模型、柔性模型进行组集并利用有限单元法对模型进行离散，得到离散后的耦合模型；

4)、将约束条件引入离散后的耦合模型，得到动力学模型；

5)、对离散后的离散后的耦合模型进行求解，通过静力学分析得到广义坐标，通过运动学分析得到广义速度，通过动力学分析得到每一时刻的广义坐标、广义速度、广义加速度，从而得到钻机钻进过程中的一系列动态响应；

6)、运用多体动力学，分别建立各个刚性体模型和柔性体模型，并对顶部组合、钻杆、底部钻具组合进行刚性连接，完成刚柔耦合模型的建立。

与现有技术比较，本发明的有益效果为：

本发明提供了一种基于多体动力学的瓦斯抽采钻机建模方法，通过对抽采钻机进行刚性和柔性进行划分和建立瓦斯抽采钻机系统刚柔耦合多体动力学模型，体现了瓦斯抽采工况的复杂性，为瓦斯抽采工艺的优化提供了方法支持；另外，本发明运用多体动力学软件对模型进行求解，提高求解效率，精确求解结果，使瓦斯抽采技术的优化或控制更容易实现，为瓦斯抽采向智能化发展提供数据支持。

附图说明

图1是本发明方法流程图；

图2是本发明刚性体上P在整体系中的位置图；

图3是本发明柔性杆上P^ij在整体系中的位置图；

图4是本发明数值分析流程图；

图5是本发明本发明的仿真分析流程图；

图6是本发明是三种情况下第100根钻杆的动能图；

图7是本发明是三种不同情况下第100根钻杆的X轴位移幅频特性图；

图8是本发明是钻杆长度对系统固有频率的影响图。

具体实施方式

下面将结合附图对本发明做详细描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。

如图1所示，本发明提供了一种基于多体动力学的瓦斯抽采钻机建模方法，包括以下步骤：

1)、将钻机、钻杆连接处、以及钻具均设为刚性体，并使用牛顿-欧拉法分别对钻机、钻杆连接处、以及钻具进行建模，建立整体坐标系，利用质量、惯性张量、广义力、以及广义力描述刚性体，得到刚性模型；

其中，所述钻机包括主机、操纵装置、电机泵组、油箱、车体、稳固装置、履带总成、电磁起动器及泥浆泵等效为一个集中质量块，将钻头简化为一段钻杆，将等效集中质量块、简化钻杆、钻杆连接处、钻具以及孔壁均视为刚性体，并使用牛顿-欧拉法进行建模；

具体的，所述步骤1)中所述刚性体模型建立具体包括：

如图2所示，选定一组物体坐标系

利用参考坐标

表示其在整体坐标系X₁X₂X₃中的位置

式中Rⁱ表示一组确定物体参考系原点位置的笛卡尔坐标，θⁱ表示一组描述物体参考系在整体坐标系位置的转动坐标，

表示P点的局部位置矢量，Aⁱ表示物体系到整体系的变换矩阵；

2)、将钻杆设为柔性体，建立参考坐标系和弹性坐标，并使用第二类拉格朗日方程进行建模，得到柔性模型；

如图3所示，所述步骤2)中所的柔性模型的建立具体包括：

柔性杆的坐标qⁱ用一组广义参考坐标矢量

和一组广义弹性坐标矢量

表示为：

式中Rⁱ、θⁱ表示参考坐标，

表示弹性坐标，柔性杆上任意一点P在整体系中的位置为：

表示点P^ij在整个系统中的位置矢量，

是点P^ij在单元坐标系中的弹性变形矢量，

是单元的节点位移矢量，d^ij是点P^ij在对象系统中的位置，

是点P^ij在物体系统中不变形时的位置，S^ij是单位函数的修正矩阵。

3)、利用浮动节点坐标法对刚性模型、柔性模型进行组集并利用有限单元法对模型进行离散，得到离散后的耦合模型；

所述组集、离散具体包括：

柔性模型离散：根据所选的有限元单元法将柔性体进行离散，分别建立整体坐标系XYZ、物体坐标系XⁱYⁱZⁱ、单元坐标系X^ijY^ijZ^ij以及中间坐标系

四组笛卡尔坐标系；物体系固连于第i根杆的某个单元上，将物体中所有单元的表达方式统一化；单元系固连于第j单元的某一点上，随单元运动；中间系固连于物体系上，其初始位置平行于单元系；

第i根柔性杆中第j个单元上的任意一点P^ij在整体系中的位置为：

表示P^ij点在整体系中的位置矢量，S^ij单元形函数矩阵，

表示单元的节点位移向量；

将式(5)对时间进行求导，得到点P^ij的速度矢量：

Aⁱ由转角坐标θⁱ确定，因此，在有n_r个转动坐标的情况下，式(7)中的

写成：

其中，

表示物体系中的转角坐标对时间的导数，P^ij在整体系中的位置按分块形式改写为：

I为3*3的单位矩阵，第i根柔性杆动能的分块形式定义为：

Tⁱ表示多体系统中，杆i的动能，

表示单元相对于物体系平动的惯性张量，

表示物体系平动和转动之间的惯性耦合，

表示物体系的平动与单元弹性形变的耦合，

表示单元相对于物体系转动的惯性张量，

表示物体系的转动与单元弹性变形之间的惯性耦合，

表示单元相对于物体系弹性变形的惯性张量；

刚性模型离散：利用有限单元法中的单元刚度组集法，并参考坐标和弹性坐标的分块形式，第i根柔性杆应变能的分块形式可定义为：

式中，

表示与第ij单元弹性坐标相关的对称正定刚度矩阵。作用在第i根柔性杆的全部外力虚功可表示为：

表示除杆间无功约束力外，与广义坐标相关的所有广义力矢量，下标e表示外力，将钻杆变形引起的弹性力所作的虚功和作用在杆上的全部外力所作的虚功相加，得到第i根杆上的所有力所作的虚功：

δWⁱ表示作用于第i根杆上的所有力所作的虚功，δUⁱ表示由变形引起的弹性力虚功，

表示外力虚功，根据式(12)和式(14)，式(15)写为：

式中，Qⁱ表示与杆i广义坐标相关的广义力矢量，由杆i的动能及广义力矢量，导出拉格朗日方程为：

Tⁱ表示杆i的动能。可进一步导出柔性杆i的动力学方程为：

表示与杆i广义坐标相关的速度二次矢量。

4)、将约束条件引入离散后的耦合模型，得到动力学模型；

约束条件具体包括：约束方程、摩擦方程、振动方程、耦合方程、以及边界条件。

所述步骤4)将约束条件引入刚性模型和离散后的柔性模型具体包括：

计算因非连续性产生的刚性，在此引入平滑函数，将库仑摩擦公式改写为：

式中，F_f表示摩擦力，μ表示接触点之间的动摩擦系数，μ_s表示最大静摩擦系数，ε、τ表示平滑系数，是无量纲数值常数，q_h表示接触点之间的相对滑动速度，滚动摩擦力为：

F_r＝δ_rF_N (21)

δ_r为钻杆与接触煤岩的滚动摩擦系数，F_N表示接触点之间的法向力，作用在钻杆上的离心力为：

m表示单根钻杆质量，e表示偏心距，ω表示转速。杆横向振动微分方程可写为：

E_wI_p表示弯曲刚度，F_H(q，t)表示施加在钻杆上的横向力，C_s表示应变内阻尼系数，ρ表示钻杆的密度，钻杆轴向振动方程为：

E_lI_p表示拉压刚度，u(q，t)表示截面在时刻t的位移，钻杆扭转振动方程写为：

GI_p表示旋转刚度，θ(q，t)表示截面在时刻t的角位移，摩擦、振动耦合后耦合项可写为：

边界条件为：

其中，+表示碰撞前，—表示碰撞后，下标x、y、z分别表示x、y、z的三个方向，δ为最小正数，

分别为碰撞和恢复阶段钻杆的加速度，

为钻杆单元与岩层接触产生的摩擦力矩，

为碰撞钻孔半径，Δ为煤岩的挤压变形，

为碰撞后的速度；

综上，动力学模型为：

C(R，θ，q_f，t)＝0 (28)

5)、对离散后的离散后的耦合模型进行求解，通过静力学分析得到广义坐标，通过运动学分析得到广义速度，通过动力学分析得到每一时刻的广义坐标、广义速度、广义加速度，从而得到钻机钻进过程中的一系列动态响应；

如图4所示，所述步骤5)具体包括：

当系统静止时，动力学模型为：

设

R_e＝Q_e-Kq (30)

将广义坐标写成分块形式：

q_i、q_d分别表示广义独立坐标和相关坐标，求约束方程对广义坐标的偏导，并写成紧凑形式为：

C_q·δq＝0 (32)

将其写为独立坐标、相关坐标分块形式为：

C_qiδq_i+C_qdδq_d＝0 (33)

C_qd为n_c×n_c阶矩阵，C_qi为n_c×(n-n_c)阶矩阵，n_c表示相互独立的约束方程数，n表示广义坐标数，总体广义坐标的虚位移可表示为：

B_di为n×(n-n_c)阶矩阵，将式(30)代入式(29)，两边乘以虚位移，并考虑到式(34)可得：

求解式(35)即可得到q_i：

设

则有

求解式(37)可得到Δq_i，其中

可利用差分法求出，以此来确定全部独立坐标q_i，同时，利用Newton-Raphson方法求解约束方程可得到相关坐标q_d，则总体广义坐标矢量为：

将广义速度分块为：

求约束方程对时间的偏导，并写成分块形式为：

独立坐标对时间求导即可得到独立广义速度，通过对式(40)进行求解，得到相关广义速度，将独立速度和相关速度合并，得到总体的广义速度矢量；

6)、运用多体动力学，分别建立各个刚性体模型和柔性体模型，并对顶部组合、钻杆、底部钻具组合进行刚性连接，完成刚柔耦合模型的建立。

根据动力学，求约束方程对时间的二次偏导，并与多体动力学方程合并得：

C_tt表示约束方程对时间的二次导数，C_q表示约束方程对广义坐标的一次导数。因此，可由上式导出

与广义坐标相似，可将广义加速度分块为：

因此构造出状态方程

通过式(41)解出广义加速度，将其返回到式(42)和式(44)，利用直接积分法求得下一点的广义速度和广义坐标。

如图5所示，本发明建立模型后，创建旋转副和滑移副，并添加旋转驱动和滑移驱动。定义接触、约束、载荷，添加边界条件、初始条件等参数，通过查询模型是否有过约束信息来验证模型是否正确，若出现错误则返回修改，验证模型正确无误后，分别添加摩擦力和输入通道、振动激励、输出通道，并对模型进行求解、后处理和优化；

从图6、图7可以看出，摩擦和振动耦合不一定会增加钻孔失效概率，有时会抑制振动的幅度。在1.58HZ之前，系统受摩擦振动(摩擦引起的振动)的影响较大。1.58HZ后，系统受偏心振动(偏心引起的振动)的影响较大。耦合振动(驱动激振、摩擦振动、偏心振动耦合)曲线，在1.58HZ前，与对系统影响较大的摩擦振动曲线相符，在1.58HZ以后，与对系统影响较大的偏心振动曲线相符。因此，分析耦合振动更为合理。

从图8可以看出，钻机工作频率范围为6.9e-3-2.1e+4。与每根钻杆的长度相比，更应考虑钻深对系统的影响。此外，浅孔段的固有频率受阶次的影响较大，而深孔段的固有频率受阶次的影响较小。在实际工程中，可以通过控制操纵装置的振动激励，来避免模态频率，从而降低振动和钻井失败的概率。

对于本领域技术人员而言，显然能了解到上述具体事实例只是本发明的优选方案，因此本领域的技术人员对本发明中的某些部分所可能作出的改进、变动，体现的仍是本发明的原理，实现的仍是本发明的目的，均属于本发明所保护的范围。

Claims (8)

1.一种基于多体动力学的瓦斯抽采钻机建模方法，其特征在于，包括以下步骤：

1)、将钻机、钻杆连接处、以及钻具均设为刚性体，并使用牛顿-欧拉法分别对钻机、钻杆连接处、以及钻具进行建模，建立整体坐标系，利用质量、惯性张量、广义力、以及广义力描述刚性体，得到刚性模型；

2)、将钻杆设为柔性体，建立参考坐标系和弹性坐标，并使用第二类拉格朗日方程进行建模，得到柔性模型；

3)、利用浮动节点坐标法对刚性模型、柔性模型进行组集并利用有限单元法对模型进行离散，得到离散后的耦合模型；

4)、将约束条件引入离散后的耦合模型，得到动力学模型；

5)、对离散后的耦合模型进行求解，通过静力学分析得到广义坐标，通过运动学分析得到广义速度，通过动力学分析得到每一时刻的广义坐标、广义速度、广义加速度，从而得到钻机钻进过程中的一系列动态响应；

6)、运用多体动力学，分别建立各个刚性体模型和柔性体模型，并对顶部组合、钻杆、底部钻具组合进行刚性连接，完成刚柔耦合模型的建立。

2.根据权利要求1所述的基于多体动力学的瓦斯抽采钻机建模方法，其特征在于，所述步骤1)中所述刚性体模型建立具体包括：

选定一组物体坐标系

利用参考坐标

表示其在整体坐标系X₁X₂X₃中的位置

式中Rⁱ表示一组确定物体参考系原点位置的笛卡尔坐标，θⁱ表示一组描述物体参考系在整体坐标系位置的转动坐标，

表示P点的局部位置矢量，Aⁱ表示物体系到整体系的变换矩阵；

3.根据权利要求2所述的基于多体动力学的瓦斯抽采钻机建模方法，其特征在于，所述步骤2)中所的柔性模型的建立具体包括：

柔性杆的坐标qⁱ用一组广义参考坐标矢量

和一组广义弹性坐标矢量

表示为：

式中Rⁱ、θⁱ表示参考坐标，

表示弹性坐标，柔性杆上任意一点P在整体系中的位置为：

表示点P^ij在整个系统中的位置矢量，

是点P^ij在单元坐标系中的弹性变形矢量，

是单元的节点位移矢量，d^ij是点P^ij在对象系统中的位置，

是点P^ij在物体系统中不变形时的位置，S^ij是单位函数的修正矩阵。

4.根据权利要求3所述的基于多体动力学的瓦斯抽采钻机建模方法，其特征在于，所述利用浮动节点坐标法对刚性模型、柔性模型进行组集并利用有限单元法对模型进行离散具体包括：

柔性模型离散：根据所选的有限元单元法将柔性体进行离散，分别建立整体坐标系XYZ、物体坐标系XⁱYⁱZⁱ、单元坐标系X^ijY^ijZ^ij以及中间坐标系

四组笛卡尔坐标系；物体系固连于第i根杆的某个单元上，将物体中所有单元的表达方式统一化；单元系固连于第j单元的某一点上，随单元运动；中间系固连于物体系上，其初始位置平行于单元系；

第i根柔性杆中第j个单元上的任意一点P^ij在整体系中的位置为：

表示P^ij点在整体系中的位置矢量，S^ij单元形函数矩阵，

表示单元的节点位移向量；

将式(5)对时间进行求导，得到点P^ij的速度矢量：

Aⁱ由转角坐标θⁱ确定，因此，在有n_r个转动坐标的情况下，式(7)中的

写成：

其中，

表示物体系中的转角坐标对时间的导数，P^ij在整体系中的位置按分块形式改写为：

I为3*3的单位矩阵，第i根柔性杆动能的分块形式定义为：

Tⁱ表示多体系统中，杆i的动能，

表示单元相对于物体系平动的惯性张量，

表示物体系平动和转动之间的惯性耦合，

表示物体系的平动与单元弹性形变的耦合，

表示单元相对于物体系转动的惯性张量，

表示物体系的转动与单元弹性变形之间的惯性耦合，

表示单元相对于物体系弹性变形的惯性张量；

刚性模型离散：利用有限单元法中的单元刚度组集法，并参考坐标和弹性坐标的分块形式，第i根柔性杆应变能的分块形式可定义为：

式中，

表示与第ij单元弹性坐标相关的对称正定刚度矩阵，作用在第i根柔性杆的全部外力虚功可表示为：

表示除杆间无功约束力外，与广义坐标相关的所有广义力矢量，下标e表示外力，将钻杆变形引起的弹性力所作的虚功和作用在杆上的全部外力所作的虚功相加，得到第i根杆上的所有力所作的虚功：

δWⁱ表示作用于第i根杆上的所有力所作的虚功，δUⁱ表示由变形引起的弹性力虚功，

表示外力虚功，根据式(12)和式(14)，式(15)写为：

式中，Qⁱ表示与杆i广义坐标相关的广义力矢量，由杆i的动能及广义力矢量，导出拉格朗日方程为：

Tⁱ表示杆i的动能，可进一步导出柔性杆i的动力学方程为：

表示与杆i广义坐标相关的速度二次矢量。

5.根据权利要求4所述的基于多体动力学的瓦斯抽采钻机建模方法，其特征在于，所述步骤4)约束条件具体包括：约束方程、摩擦方程、振动方程、耦合方程、以及边界条件。

6.根据权利要求1或5所述的基于多体动力学的瓦斯抽采钻机建模方法，其特征在于，所述步骤4)将约束条件引入离散后的耦合模型具体包括：

计算因非连续性产生的刚性，在此引入平滑函数，将库仑摩擦公式改写为：

式中，F_f表示摩擦力，μ表示接触点之间的动摩擦系数，μ_s表示最大静摩擦系数，ε、τ表示平滑系数，是无量纲数值常数，q_h表示接触点之间的相对滑动速度，滚动摩擦力为：

F_r＝δ_rF_N (21)

δ_r为钻杆与接触煤岩的滚动摩擦系数，F_N表示接触点之间的法向力，作用在钻杆上的离心力为：

m表示单根钻杆质量，e表示偏心距，ω表示转速，杆横向振动微分方程可写为：

E_wI_p表示弯曲刚度，F_H(q，t)表示施加在钻杆上的横向力，Cs表示应变内阻尼系数，ρ表示钻杆的密度，钻杆轴向振动方程为：

E_lI_p表示拉压刚度，u(q，t)表示截面在时刻t的位移，钻杆扭转振动方程写为：

GI_p表示旋转刚度，θ(q，t)表示截面在时刻t的角位移，摩擦、振动耦合后耦合项可写为：

边界条件为：

其中，+表示碰撞前，—表示碰撞后，下标x、y、z分别表示x、y、z的三个方向，δ为最小正数，

分别为碰撞和恢复阶段钻杆的加速度，

为钻杆单元与岩层接触产生的摩擦力矩，

为碰撞钻孔半径，Δ为煤岩的挤压变形，

为碰撞后的速度；

综上，动力学模型为：

C(R，θ，q_f，t)＝0 (28)。

7.根据权利要求6所述的基于多体动力学的瓦斯抽采钻机建模方法，其特征在于，所述步骤5)具体包括：

当系统静止时，动力学模型为：

设

R_e＝Q_e-Kq (30)

将广义坐标写成分块形式：

q_i、q_d分别表示广义独立坐标和相关坐标，求约束方程对广义坐标的偏导，并写成紧凑形式为：

C_q·δq＝0 (32)

将其写为独立坐标、相关坐标分块形式为：

C_qiδq_i+C_qdδq_d＝0 (33)

C_qd为n_c×n_c阶矩阵，C_qi为n_c×(n-n_c)阶矩阵，n_c表示相互独立的约束方程数，n表示广义坐标数，总体广义坐标的虚位移可表示为：

B_di为n×(n-n_c)阶矩阵，将式(30)代入式(29)，两边乘以虚位移，并考虑到式(34)可得：

求解式(35)即可得到q_i：

设

则有

求解式(37)可得到Δq_i，其中

可利用差分法求出，以此来确定全部独立坐标q_i，同时，利用Newton-Raphson方法求解约束方程可得到相关坐标q_d，则总体广义坐标矢量为：

将广义速度分块为：

求约束方程对时间的偏导，并写成分块形式为：

独立坐标对时间求导即可得到独立广义速度，通过对式(40)进行求解，得到相关广义速度，将独立速度和相关速度合并，得到总体的广义速度矢量。

8.根据权利要求7所述的基于多体动力学的瓦斯抽采钻机建模方法，其特征在于，所述步骤6)具体包括：

根据动力学，求约束方程对时间的二次偏导，并与多体动力学方程合并得：

C_tt表示约束方程对时间的二次导数，C_q表示约束方程对广义坐标的一次导数，由上式导出

与广义坐标相似，将广义加速度分块为：

构造出状态方程

通过式(41)解出广义加速度，将其返回到式(42)和式(44)，利用直接积分法求得下一点的广义速度和广义坐标。

Images