一种瓦斯抽采钻机系统的多体动力学分析方法与流程

本发明涉及刻钻孔失效的风险预防技术领域，尤其涉及一种瓦斯抽采钻机系统的多体动力学分析方法。

背景技术：

煤炭资源是世界主要能源之一。近年来，数次煤炭开采工程所引发的瓦斯爆炸灾难引发了广泛关注。有效的抽采煤层瓦斯可以将瓦斯含量和压力降低到安全水平，从根本上消除瓦斯灾害。

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

技术实现要素：

针对上述缺陷或不足，本发明的目的在于提供一种瓦斯抽采钻机系统的多体动力学分析方法，使瓦斯抽采技术的优化或控制更容易实现，为瓦斯抽采向智能化发展提供数据支持。

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

一种瓦斯抽采钻机系统的多体动力学分析方法，包括以下步骤：

s1、静力学分析：构建瓦斯抽采钻机系统的静力学方程，然后利用牛顿-拉夫逊newton-raphson法和差分法求出静力学方程的近似解，并通过静力学分析得到广义坐标；

s2、运动学分析：在静力学分析的前提下，构建瓦斯抽采钻机系统的运动学方程，然后利用代入法求解运动学方程，并通过运动学分析得到广义速度；

s3、动力学分析：构建瓦斯抽采钻机系统的动力学方程，利用直接积分和构造状态方程法求解动力学方程，并通过动力学分析得到每一时刻的广义坐标、广义速度、广义加速度，从而得到钻机钻进过程中的一系列动态响应。

具体地，所述步骤s1具体包括：

当系统静止时，构建的瓦斯抽采钻机系统的静力学方程：

设

re＝qe+fmz-kq(2)

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

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

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

式中，为nc×n-nc阶独立坐标约束的雅可比矩阵，为nc×nc阶相关坐标约束的雅可比矩阵，nc表示相互独立的约束方程数，n表示广义坐标数，总体广义坐标的虚位移表示为：

bdi为n×(n-nc)阶矩阵，将式(6)代入式(5)，两边乘以虚位移得：

求解式(7)即得到qi；

设

则有：

求解式(9)得到δqi，其中利用差分法求出，以此来确定全部独立坐标qi，同时，利用牛顿-拉夫逊newton-raphson法求解约束方程得到相关坐标qd，则总体广义坐标矢量为：

具体地，所述步骤s2具体包括：

将广义速度分块为：

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

参考静力学分析的计算方法，参考(1)和(3)，由式(12)推导出相关速度与独立速度的关系为：

独立速度由静力学分析得出的独立坐标直接求导而得，将其代入式(6)求得相关速度，按照式(11)的形式整合独立速度和相关速度，即得到整个系统的初始广义速度。

具体地，所述步骤s3具体包括：

构建瓦斯抽采钻机系统的运动学方程：

式中，m和k分别为系统的质量矩阵和刚度矩阵；和λ分别为约束的雅可比矩阵和拉格朗日乘子；qe和qv分别表示系统的广义外力矢量和速度的二次矢量；为相互独立的约束方程组；为与广义坐标q和时间t相关的函数，nc表示约束方程的个数；

将式(13)对时间进行求导，即对约束方程组求二阶导数得：

式中，下标q表示雅可比矩阵，如的雅可比矩阵；下标t表示对时间的导数，如为对时间的二阶偏导数；和分别为系统的广义速度和广义加速度；

式(15)整理为：

设：

则式(16)写为：

将式(14)中的动力学方程和式(18)整合成矩阵形式为：

式中：ctt表示约束方程对时间的二次导数，cq表示约束方程对广义坐标的一次导数；

设：

式(19)写成紧凑形式为：

求解式(21)得：

将式(22)写为分块形式为：

式中：和均为与广义坐标相关的函数，表示的分量，其中，为mλ的逆矩阵；和均为与广义坐标、广义速度和时间相关的函数；

求解式(23)得加速度矢量和拉格朗日乘子为：

由式(24)知广义加速度与广义坐标、广义速度和时间相关，因此将式(24)写为更加简洁的形式为：

将广义加速度写成相应的分块形式为：

构造出状态方程：

对式(29)直接积分求得下一时刻的独立坐标和独立速度，其中初始时刻的广义坐标、速度和加速度分别由静力学分析、运动学分析和动力学分析得到；

随后，将所求得的独立坐标和独立速度分别代入式(6)和式(13)，得该时刻的相关坐标和相关速度；

整合所求得的独立坐标和相关坐标即得到该时刻的整体广义坐标，整合所求得的独立速度和相关速度即得到该时刻的整体广义速度；此时，若运算时间大于所需要的终止时间，则可整合输出全部时刻的整体广义坐标、速度和加速度；若运算时间小于所需要的终止时间，则将该时刻的整体广义坐标代入步骤s2中的运动学分析中，再次对系统进行运动学和动力学分析，直至运算时间大于所需要的终止时间。

与现有技术比较，本发明提供了一种瓦斯抽采钻机系统的多体动力学分析方法，通过对瓦斯抽采钻机系统静力学、运动学和动力学分析，提高瓦斯抽采钻机系统的动力学模型求解效率，精确求解结果，使瓦斯抽采技术的优化或控制更容易实现，为瓦斯抽采向智能化发展提供数据支持；另外，通过对瓦斯抽采钻机系统刚柔耦合多体动力学模型进行求解，求解结果体现了瓦斯抽采工况的复杂性，为瓦斯抽采工艺的优化提供了方法支持，进而根据本发明的分析方法可以预先调整钻进参数来降低该时刻钻孔失效的风险，从而保证钻进过程中钻孔的稳定性。

附图说明

图1是本发明一种实施例的方法流程图。

图2是本发明一种实施例的初始位置的离散图。

图3是本发明一种实施例的第1根钻杆的z向位移的结果图。

图4是本发明一种实施例的第1根钻杆的z向速度的结果图。

图5是本发明一种实施例的第1根钻杆的z向加速度的结果图。

具体实施方式

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

为了更加确切验证本发明，本实施例以图2所示的一种瓦斯抽采钻机系统为例，柔性体以第1根钻杆为例，刚性体以钻头为例，其余物体的预估值同理可得，钻杆和钻头的初始位置图如2所示，其中，将每个离散柔性物体均用4个节点来标识，将每个刚性物体力学性质均简化到质心上，以第1根钻杆为起始物体，从节点1开始标识，节点1为钻机和第1根钻杆的连接点，节点2、3为第一根钻杆离散后的标识点，节点4为第1根钻杆和第2根钻杆的连接点，以此类推，钻头利用节点0表示。为第i根钻杆以及钻头的物体坐标系，也可以用一组表示位置和方向的数组代替，其中，表示第i根钻杆和钻头的物体坐标系相对于整体坐标系的位置坐标，θⁱ表示第i根钻杆和钻头在整体坐标系中的方位。为第i根钻杆第j个单元的单元坐标系，同样可以用一个数组(ql，qm，qn)来代替，例如(q4，q5，q6)表示节点2的三个弹性位移方向。

在上述如图2所示的瓦斯抽采钻机系统的基础上，参照图1所示，本实施例的一种瓦斯抽采钻机系统的多体动力学分析方法，其具体步骤如下：

1、预估每个物体初始位置的广义坐标和广义速度。以图2中的第1根钻杆为例进行数值计算，其初始广义坐标和广义速度的预估值可为表示：

式中，q的单位为m，的单位为rad/s

系统的静力学模型为：

设

re＝qe+fmz-kq(4)

则式(3)的静力学方程可改写为：

广义坐标的分块形式可写为：

式中，qi和qd分别表示独立坐标和相关坐标。

将式(5)的两边分别乘以虚位移δq可得：

式(3)中约束方程对广义坐标求导的紧凑形式为：

将其按独立坐标和相关坐标进行分块可得：

式中，和分别为nc×(n-nc)阶独立坐标约束的雅可比矩阵和nc×nc阶相关坐标约束的雅可比矩阵，nc表示相互独立的约束方程的个数，n表示广义坐标的总数。

由于式(9)中约的束方程均具有独立性，因此，其雅可比矩阵均为满秩矩阵，即存在逆矩阵，故式(9)可改写为：

令

则式(10)可写为：

考虑到式(12)，总体虚位移可表示为：

其中，

式中，bdi为n×(n-nc)阶矩阵。

将式(13)代入式(7)可得：

考虑到qi的独立性，式(15)可改写为：

考虑到式(11)和式(14)，式(16)左边展开后的第一项可写为：

因此，将式(17)代入式(16)可得

-rebdi＝0(18)

将式(4)和式(14)代入式(18)，可得

考虑到式(6)将广义坐标写成分块形式，并代入式(12)消去方程组中的相关坐标项，即可确定初始广义独立坐标qi，从而通过式(12)求出初始广义相关坐标qd，整合独立坐标和相关坐标得出初始广义坐标，记为q⁽⁰⁾。

因此，设

则有

由差分法可知，当选用hx作为等步距时，式(21)中的系数矩阵为：

将式(22)的系数矩阵代入式(21)后，即可通过求解式(21)得到初始差值δq⁽⁰⁾。

考虑到式(3)中的约束方程，由牛顿-拉夫逊法可知，当非线性约束方程中的q为单变量时，设q^*为该约束方程的准确解，q⁽ⁱ⁾为q^*的近似值，即第i次迭代的近似解，此时有：

式中，为系统的非线性约束方程，其中，q⁽ⁱ⁾为约束方程的解；表示约束方程的雅可比矩阵。

分别设εe和εs为方程和解的误差范围，当且|q⁽ⁱ⁾-q^(i-1)|<εs时，迭代终止。当而时，则需要重新确定系统的初值，在修改好初值后，重新进行迭代，直至求得满足条件的的解。若该算法是收敛点，则该迭代的收敛阶数为2，即

|q⁽ⁱ⁺¹⁾-q^*|<β|q⁽ⁱ⁾-q^*|²(24)

式中，β为与i无关的常数；i为迭代的阶数。

当非线性约束方程组中的q为多变量时,则有：

设设为该约束方程组的准确解，q^(k)∈d为q^*的近似值，即第k次迭代的近似解。因此，可将q^(k)的仿射映像lk：rⁿ→rⁿ定义为：

式中，ak为n×n阶的非奇异矩阵。

考虑到式(26)，将线性方程组的解q＝q^(k+1)设为非线性约束方程组的解q^*的新近似值，即非线性约束方程组的线性迭代法可表示为：

考虑到n维平行弦法，若选取全部ak≡a，则式(27)可改写为：

此时，在欧式空间rⁿ⁺¹中的n个超平面i＝1，2，…，n与超平面z＝0的交点就是q^(k+1)，其中，a＝a^ij为n×n阶矩阵。上标i表示第i个超平面，上标j表示第j个单元。

考虑到简化newton法，当选取时，约束方程组的解q^*的近似值q^(k+1)可由式(28)可改写为：

为的雅可比矩阵，可写为：

由于上述newton-raphson方法可知，每次循环都要计算式(29)中的逆矩阵，n的值过大将会导致求解困难，因此，在实际计算时，利用n阶线性方程组来代替逆矩阵，式(29)可改写为：

式中，δq^(k)表示第k次迭代的矫正量，通常被称为newton差。

将已求得的初值q⁽⁰⁾和δq⁽⁰⁾代入式(31)，通过k＝0，1，2，…次迭代求解式(3.31)中的近似值q^(k+1)，即可确定静态系统的整体广义坐标q。期间，与单变量系统相似，若两次相邻迭代解的误差满足且j＝1，2，…，n，则迭代终止。否则，则需要在是奇异矩阵且的前提下重新估算初值并进行迭代。

2、在静力学分析的前提下，通过对瓦斯抽采钻机系统进行运动学分析，来求解系统运动的广义速度。此时约束方程组可写为：

将其对时间进行求导可得：

式中，为系统的广义速度；为约束方程(3.32)对时间的偏导数。

考虑到式(6)，将广义速度写成相应的分块形式为：

将式(34)代入式(33)可得：

参考静力学分析的计算方法，按式(9)到(12)的步骤，可由式(35)推导出相关速度与独立速度的关系为：

独立速度可由静力学分析得出的独立坐标直接求导而得，将其代入式(12)可求得相关速度，按照式(34)的形式整合独立速度和相关速度，即可得到整个系统的初始广义速度。

3、在静力学分析和运动学分析的前提下，通过对瓦斯抽采钻机系统进行动力学分析，来求解系统在运动时每一时刻的广义坐标、广义速度和广义加速度。系统的动力学模型为：

式中，m和k分别为系统的质量矩阵和刚度矩阵；和λ分别为约束的雅可比矩阵和拉格朗日乘子；qe和qv分别表示系统的广义外力矢量和速度的二次矢量；为相互独立的约束方程组；为与广义坐标q和时间t相关的函数，nc表示约束方程的个数。

将式(36)对时间进行求导，即对约束方程组求二阶导数可得：

式中，下标q表示雅可比矩阵，如为的雅可比矩阵；下标t表示对时间的导数，如为对时间的二阶偏导数；和分别为系统的广义速度和广义加速度。

式(38)可整理为：

设

则式(39)可写为：

将式(37)中的动力学方程和式(41)整合成矩阵形式为：

式中：ctt表示约束方程对时间的二次导数，cq表示约束方程对广义坐标的一次导数。

设：

式(42)可写成紧凑形式为：

求解式(44)可得：

将其写为分块形式为：

式中，和均为与广义坐标相关的函数，表示的分量，其中，为mλ的逆矩阵；和均为与广义坐标、广义速度和时间相关的函数。

求解式(46)可得加速度矢量和拉格朗日乘子为：

由式(47)可知广义加速度与广义坐标、广义速度和时间相关，因此可将式(47)写为更加简洁的形式为：

考虑到式(6)和式(34)，将广义加速度写成相应的分块形式为：

因此，可将状态方程构造为：

对式(52)直接积分求得下一时刻的独立坐标和独立速度，其中初始时刻的广义坐标、速度和加速度分别由静力学分析、运动学分析和动力学分析得到。随后，将所求得的独立坐标和独立速度分别代入式(12)和式(36)，可得该时刻的相关坐标和相关速度。整合所求得的独立坐标和相关坐标即可得到该时刻的整体广义坐标，整合所求得的独立速度和相关速度即可得到该时刻的整体广义速度。此时，若运算时间大于所需要的终止时间，则可整合输出全部时刻的整体广义坐标、速度和加速度。若运算时间小于所需要的终止时间，则将该时刻的整体广义坐标代入运动学分析中，再次对系统进行运动学和动力学分析，直至运算时间大于所需要的终止时间。以第1根钻杆为例，利用matlab软件计算了其在下钻过程中的z向位移、速度和加速度分别如图3、4和5所示。

由图可知，刚柔耦合模型中钻杆的z向位移受钻杆重力的影响较为严重，而z向速度和加速度受钻杆重力的影响相对较小。此外，三条曲线对于振动和摩擦的耦合影响均较为敏感，曲线的波动无明显规律，幅值的突变现象时有发生。由牛顿第二定律可知，力和加速度成正比，因此，由加速度曲线可以更直观的看出力的突变发生在哪一时刻，从而通过预先调整钻进参数来降低该时刻钻孔失效的风险，从而保证钻进过程中钻孔的稳定性。

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