CN111723536A - 一种瓦斯抽采钻机系统的多体动力学分析方法 - Google Patents

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

Info

Publication number
CN111723536A
CN111723536A CN202010550055.5A CN202010550055A CN111723536A CN 111723536 A CN111723536 A CN 111723536A CN 202010550055 A CN202010550055 A CN 202010550055A CN 111723536 A CN111723536 A CN 111723536A
Authority
CN
China
Prior art keywords
equation
generalized
coordinates
speed
independent
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
CN202010550055.5A
Other languages
English (en)
Other versions
CN111723536B (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.)
Lingnan Normal University
Original Assignee
Lingnan Normal 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 Lingnan Normal University filed Critical Lingnan Normal University
Priority to CN202010550055.5A priority Critical patent/CN111723536B/zh
Publication of CN111723536A publication Critical patent/CN111723536A/zh
Application granted granted Critical
Publication of CN111723536B publication Critical patent/CN111723536B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • E—FIXED CONSTRUCTIONS
    • E21—EARTH OR ROCK DRILLING; MINING
    • E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
    • E21B44/00—Automatic control systems specially adapted for drilling operations, i.e. self-operating systems which function to carry out or modify a drilling operation without intervention of a human operator, e.g. computer-controlled drilling systems; Systems specially adapted for monitoring a plurality of drilling variables or conditions
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Theoretical Computer Science (AREA)
  • Fluid Mechanics (AREA)
  • Geology (AREA)
  • Mining & Mineral Resources (AREA)
  • General Physics & Mathematics (AREA)
  • Computing Systems (AREA)
  • Pure & Applied Mathematics (AREA)
  • Geochemistry & Mineralogy (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Algebra (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Environmental & Geological Engineering (AREA)
  • Complex Calculations (AREA)

Abstract

本发明公开了一种瓦斯抽采钻机系统的多体动力学分析方法，包括：利用牛顿‑拉夫逊Newton‑Raphson和差分法求出静力学方程的近似解；利用代入法求解运动学方程；利用直接积分和构造状态方程法求解动力学方程，并通过动力学分析得到每一时刻的广义坐标、广义速度、广义加速度，从而得到钻机钻进过程中的一系列动态响应。与现有技术相比，根据本发明的分析方法可以预先调整钻进参数来降低该时刻钻孔失效的风险，从而保证钻进过程中钻孔的稳定性。

Description

一种瓦斯抽采钻机系统的多体动力学分析方法
技术领域
本发明涉及刻钻孔失效的风险预防技术领域，尤其涉及一种瓦斯抽采钻机系统的多体动力学分析方法。
背景技术
煤炭资源是世界主要能源之一。近年来，数次煤炭开采工程所引发的瓦斯爆炸灾难引发了广泛关注。有效的抽采煤层瓦斯可以将瓦斯含量和压力降低到安全水平，从根本上消除瓦斯灾害。
近年来，随着煤层勘探开发规模的扩大，瓦斯抽采的成孔率问题受到了国内外学者的持续关注。但是，由于钻孔失效是由多个因素相互影响、相互促进的复杂的系统性问题，涉及到摩擦学、振动学、运动学、材料科学、化学、传热学等一系列学科，存在许多不确定性和随机性，因此，关于瓦斯抽采钻进方面的研究工作还需要进一步完善。
发明内容
针对上述缺陷或不足，本发明的目的在于提供一种瓦斯抽采钻机系统的多体动力学分析方法，使瓦斯抽采技术的优化或控制更容易实现，为瓦斯抽采向智能化发展提供数据支持。
为达到以上目的，本发明的技术方案为：
一种瓦斯抽采钻机系统的多体动力学分析方法，包括以下步骤：
S1、静力学分析：构建瓦斯抽采钻机系统的静力学方程，然后利用牛顿-拉夫逊Newton-Raphson法和差分法求出静力学方程的近似解，并通过静力学分析得到广义坐标；
S2、运动学分析：在静力学分析的前提下，构建瓦斯抽采钻机系统的运动学方程，然后利用代入法求解运动学方程，并通过运动学分析得到广义速度；
S3、动力学分析：构建瓦斯抽采钻机系统的动力学方程，利用直接积分和构造状态方程法求解动力学方程，并通过动力学分析得到每一时刻的广义坐标、广义速度、广义加速度，从而得到钻机钻进过程中的一系列动态响应。
具体地，所述步骤S1具体包括：
当系统静止时，构建的瓦斯抽采钻机系统的静力学方程：
Figure BDA0002542201150000021
设
Re＝Qe+Fmz-Kq (2)
将广义坐标写成分块形式：
Figure BDA0002542201150000022
qi、qd分别表示广义独立坐标和相关坐标，求约束方程对广义坐标的偏导，并写成紧凑形式为：
Figure BDA0002542201150000023
将其写为独立坐标、相关坐标分块形式为：
Figure BDA0002542201150000024
式中，
Figure BDA0002542201150000025
为nc×n-nc阶独立坐标约束的雅可比矩阵，
Figure BDA0002542201150000026
为nc×nc阶相关坐标约束的雅可比矩阵，nc表示相互独立的约束方程数，n表示广义坐标数，总体广义坐标的虚位移表示为：
Figure BDA0002542201150000027
Bdi为n×(n-nc)阶矩阵，将式(6)代入式(5)，两边乘以虚位移得：
Figure BDA0002542201150000031
求解式(7)即得到qi；
设
Figure BDA0002542201150000032
则有：
Figure BDA0002542201150000033
求解式(9)得到Δqi，其中
Figure BDA0002542201150000034
利用差分法求出，以此来确定全部独立坐标qi，同时，利用牛顿-拉夫逊Newton-Raphson法求解约束方程得到相关坐标qd，则总体广义坐标矢量为：
Figure BDA0002542201150000035
具体地，所述步骤S2具体包括：
将广义速度分块为：
Figure BDA0002542201150000036
求约束方程对时间的偏导，并写成分块形式为：
Figure BDA0002542201150000037
参考静力学分析的计算方法，参考(1)和(3)，由式(12)推导出相关速度与独立速度的关系为：
Figure BDA0002542201150000038
独立速度由静力学分析得出的独立坐标直接求导而得，将其代入式(6)求得相关速度，按照式(11)的形式整合独立速度和相关速度，即得到整个系统的初始广义速度。
具体地，所述步骤S3具体包括：
构建瓦斯抽采钻机系统的运动学方程：
Figure BDA0002542201150000041
式中，M和K分别为系统的质量矩阵和刚度矩阵；
Figure BDA0002542201150000042
和λ分别为约束的雅可比矩阵和拉格朗日乘子；Qe和Qv分别表示系统的广义外力矢量和速度的二次矢量；
Figure BDA0002542201150000043
为相互独立的约束方程组；
Figure BDA0002542201150000044
为与广义坐标q和时间t相关的函数，nc表示约束方程的个数；
将式(13)对时间进行求导，即对约束方程组求二阶导数得：
Figure BDA0002542201150000045
式中，下标q表示雅可比矩阵，如
Figure BDA0002542201150000046
的雅可比矩阵；下标t表示对时间的导数，如
Figure BDA0002542201150000047
为
Figure BDA0002542201150000048
对时间的二阶偏导数；
Figure BDA0002542201150000049
和
Figure BDA00025422011500000410
分别为系统的广义速度和广义加速度；
式(15)整理为：
Figure BDA00025422011500000411
设：
Figure BDA00025422011500000412
则式(16)写为：
Figure BDA00025422011500000413
将式(14)中的动力学方程和式(18)整合成矩阵形式为：
Figure BDA00025422011500000414
式中：Ctt表示约束方程对时间的二次导数，Cq表示约束方程对广义坐标的一次导数；
设：
Figure BDA00025422011500000415
式(19)写成紧凑形式为：
Figure BDA0002542201150000051
求解式(21)得：
Figure BDA0002542201150000052
将式(22)写为分块形式为：
Figure BDA0002542201150000053
式中：
Figure BDA0002542201150000054
和
Figure BDA0002542201150000055
均为与广义坐标相关的函数，表示
Figure BDA0002542201150000056
的分量，其中，
Figure BDA0002542201150000057
为Mλ的逆矩阵；
Figure BDA0002542201150000058
和
Figure BDA0002542201150000059
均为与广义坐标、广义速度和时间相关的函数；
求解式(23)得加速度矢量和拉格朗日乘子为：
Figure BDA00025422011500000510
Figure BDA00025422011500000511
由式(24)知广义加速度与广义坐标、广义速度和时间相关，因此将式(24)写为更加简洁的形式为：
Figure BDA00025422011500000512
将广义加速度写成相应的分块形式为：
Figure BDA00025422011500000513
构造出状态方程：
Figure BDA00025422011500000514
Figure BDA00025422011500000515
对式(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表示。
Figure BDA0002542201150000071
为第i根钻杆以及钻头的物体坐标系，也可以用一组表示位置和方向的数组
Figure BDA0002542201150000072
代替，其中，
Figure BDA0002542201150000073
表示第i根钻杆和钻头的物体坐标系相对于整体坐标系的位置坐标，θi表示第i根钻杆和钻头在整体坐标系中的方位。
Figure BDA0002542201150000074
为第i根钻杆第j个单元的单元坐标系，同样可以用一个数组(ql，qm，qn)来代替，例如(q4，q5，q6)表示节点2的三个弹性位移方向。
在上述如图2所示的瓦斯抽采钻机系统的基础上，参照图1所示，本实施例的一种瓦斯抽采钻机系统的多体动力学分析方法，其具体步骤如下：
1、预估每个物体初始位置的广义坐标和广义速度。以图2中的第1根钻杆为例进行数值计算，其初始广义坐标和广义速度的预估值可为表示：
Figure BDA0002542201150000075
Figure BDA0002542201150000081
Figure BDA0002542201150000082
式中，q的单位为m，
Figure BDA0002542201150000083
的单位为rad/s
系统的静力学模型为：
Figure BDA0002542201150000084
设
Re＝Qe+Fmz-Kq (4)
则式(3)的静力学方程可改写为：
Figure BDA0002542201150000085
广义坐标的分块形式可写为：
Figure BDA0002542201150000086
式中，qi和qd分别表示独立坐标和相关坐标。
将式(5)的两边分别乘以虚位移δq可得：
Figure BDA0002542201150000087
式(3)中约束方程对广义坐标求导的紧凑形式为：
Figure BDA0002542201150000088
将其按独立坐标和相关坐标进行分块可得：
Figure BDA0002542201150000089
式中，
Figure BDA0002542201150000091
和
Figure BDA0002542201150000092
分别为nc×(n-nc)阶独立坐标约束的雅可比矩阵和nc×nc阶相关坐标约束的雅可比矩阵，nc表示相互独立的约束方程的个数，n表示广义坐标的总数。
由于式(9)中约的束方程均具有独立性，因此，其雅可比矩阵均为满秩矩阵，即存在逆矩阵，故式(9)可改写为：
Figure BDA0002542201150000093
令
Figure BDA0002542201150000094
则式(10)可写为：
Figure BDA0002542201150000095
考虑到式(12)，总体虚位移可表示为：
Figure BDA0002542201150000096
其中，
Figure BDA0002542201150000097
式中，Bdi为n×(n-nc)阶矩阵。
将式(13)代入式(7)可得：
Figure BDA0002542201150000098
考虑到qi的独立性，式(15)可改写为：
Figure BDA0002542201150000099
考虑到式(11)和式(14)，式(16)左边展开后的第一项
Figure BDA00025422011500000910
可写为：
Figure BDA0002542201150000101
因此，将式(17)代入式(16)可得
-ReBdi＝0 (18)
将式(4)和式(14)代入式(18)，可得
Figure BDA0002542201150000102
考虑到式(6)将广义坐标写成分块形式，并代入式(12)消去方程组中的相关坐标项，即可确定初始广义独立坐标qi，从而通过式(12)求出初始广义相关坐标qd，整合独立坐标和相关坐标得出初始广义坐标，记为q(0)。
因此，设
Figure BDA0002542201150000103
则有
Figure BDA0002542201150000104
由差分法可知，当选用hx作为等步距时，式(21)中的系数矩阵为：
Figure BDA0002542201150000105
将式(22)的系数矩阵代入式(21)后，即可通过求解式(21)得到初始差值Δq(0)。
考虑到式(3)中的约束方程，由牛顿-拉夫逊法可知，当非线性约束方程
Figure BDA00025422011500001013
中的q为单变量时，设q*为该约束方程的准确解，q(i)为q*的近似值，即第i次迭代的近似解，此时有：
Figure BDA0002542201150000106
式中，
Figure BDA0002542201150000107
为系统的非线性约束方程，其中，q(i)为约束方程的解；
Figure BDA0002542201150000108
表示约束方程的雅可比矩阵。
分别设εe和εs为方程和解的误差范围，当
Figure BDA0002542201150000109
且|q(i)-q(i-1)|<εs时，迭代终止。当
Figure BDA00025422011500001010
而
Figure BDA00025422011500001011
时，则需要重新确定系统的初值，在修改好初值后，重新进行迭代，直至求得满足条件的
Figure BDA00025422011500001012
的解。若该算法是收敛点，则该迭代的收敛阶数为2，即
|q(i+1)-q*|<β|q(i)-q*|2 (24)
式中，β为与i无关的常数；i为迭代的阶数。
当非线性约束方程组
Figure BDA00025422011500001117
中的q为多变量时,则有：
Figure BDA0002542201150000111
设
Figure BDA0002542201150000112
设
Figure BDA0002542201150000113
为该约束方程组的准确解，q(k)∈D为q*的近似值，即第k次迭代的近似解。因此，可将q(k)的仿射映像Lk：Rn→Rn定义为：
Figure BDA0002542201150000114
式中，Ak为n×n阶的非奇异矩阵。
考虑到式(26)，将线性方程组
Figure BDA0002542201150000115
的解q＝q(k +1)设为非线性约束方程组
Figure BDA0002542201150000116
的解q*的新近似值，即非线性约束方程组
Figure BDA0002542201150000117
的线性迭代法可表示为：
Figure BDA0002542201150000118
考虑到n维平行弦法，若选取全部Ak≡A，则式(27)可改写为：
Figure BDA0002542201150000119
此时，在欧式空间Rn+1中的n个超平面
Figure BDA00025422011500001110
Figure BDA00025422011500001111
i＝1，2，…，n与超平面Z＝0的交点就是q(k+1)，其中，A＝aij为n×n阶矩阵。上标i表示第i个超平面，上标j表示第j个单元。
考虑到简化Newton法，当选取
Figure BDA00025422011500001112
时，约束方程组的解q*的近似值q(k+1)可由式(28)可改写为：
Figure BDA00025422011500001113
Figure BDA00025422011500001114
为
Figure BDA00025422011500001115
的雅可比矩阵，可写为：
Figure BDA00025422011500001116
由于上述Newton-Raphson方法可知，每次循环都要计算式(29)中
Figure BDA0002542201150000121
的逆矩阵，n的值过大将会导致求解困难，因此，在实际计算时，利用n阶线性方程组来代替逆矩阵，式(29)可改写为：
Figure BDA0002542201150000122
式中，Δq(k)表示第k次迭代的矫正量，通常被称为Newton差。
将已求得的初值q(0)和Δq(0)代入式(31)，通过k＝0，1，2，…次迭代求解式(3.31)中的近似值q(k+1)，即可确定静态系统的整体广义坐标q。期间，与单变量系统相似，若两次相邻迭代解的误差满足
Figure BDA0002542201150000123
且
Figure BDA0002542201150000124
Figure BDA0002542201150000125
j＝1，2，…，n，则迭代终止。否则，则需要在
Figure BDA0002542201150000126
是奇异矩阵且
Figure BDA0002542201150000127
的前提下重新估算初值并进行迭代。
2、在静力学分析的前提下，通过对瓦斯抽采钻机系统进行运动学分析，来求解系统运动的广义速度。此时约束方程组可写为：
Figure BDA0002542201150000128
将其对时间进行求导可得：
Figure BDA0002542201150000129
式中，
Figure BDA00025422011500001210
为系统的广义速度；
Figure BDA00025422011500001211
为约束方程(3.32)对时间的偏导数。
考虑到式(6)，将广义速度写成相应的分块形式为：
Figure BDA00025422011500001212
将式(34)代入式(33)可得：
Figure BDA00025422011500001213
参考静力学分析的计算方法，按式(9)到(12)的步骤，可由式(35)推导出相关速度与独立速度的关系为：
Figure BDA00025422011500001214
独立速度可由静力学分析得出的独立坐标直接求导而得，将其代入式(12)可求得相关速度，按照式(34)的形式整合独立速度和相关速度，即可得到整个系统的初始广义速度。
3、在静力学分析和运动学分析的前提下，通过对瓦斯抽采钻机系统进行动力学分析，来求解系统在运动时每一时刻的广义坐标、广义速度和广义加速度。系统的动力学模型为：
Figure BDA0002542201150000131
式中，M和K分别为系统的质量矩阵和刚度矩阵；
Figure BDA0002542201150000132
和λ分别为约束的雅可比矩阵和拉格朗日乘子；Qe和Qv分别表示系统的广义外力矢量和速度的二次矢量；
Figure BDA0002542201150000133
为相互独立的约束方程组；
Figure BDA0002542201150000134
为与广义坐标q和时间t相关的函数，nc表示约束方程的个数。
将式(36)对时间进行求导，即对约束方程组求二阶导数可得：
Figure BDA0002542201150000135
式中，下标q表示雅可比矩阵，如
Figure BDA0002542201150000136
为
Figure BDA0002542201150000137
的雅可比矩阵；下标t表示对时间的导数，如
Figure BDA0002542201150000138
为
Figure BDA0002542201150000139
对时间的二阶偏导数；
Figure BDA00025422011500001310
和
Figure BDA00025422011500001311
分别为系统的广义速度和广义加速度。
式(38)可整理为：
Figure BDA00025422011500001312
设
Figure BDA00025422011500001313
则式(39)可写为：
Figure BDA00025422011500001314
将式(37)中的动力学方程和式(41)整合成矩阵形式为：
Figure BDA00025422011500001315
式中：Ctt表示约束方程对时间的二次导数，Cq表示约束方程对广义坐标的一次导数。
设：
Figure BDA00025422011500001316
Figure BDA0002542201150000141
式(42)可写成紧凑形式为：
Figure BDA0002542201150000142
求解式(44)可得：
Figure BDA0002542201150000143
将其写为分块形式为：
Figure BDA0002542201150000144
式中，
Figure BDA0002542201150000145
和
Figure BDA0002542201150000146
均为与广义坐标相关的函数，表示
Figure BDA0002542201150000147
的分量，其中，
Figure BDA0002542201150000148
为Mλ的逆矩阵；
Figure BDA0002542201150000149
和
Figure BDA00025422011500001410
均为与广义坐标、广义速度和时间相关的函数。
求解式(46)可得加速度矢量和拉格朗日乘子为：
Figure BDA00025422011500001411
Figure BDA00025422011500001412
由式(47)可知广义加速度与广义坐标、广义速度和时间相关，因此可将式(47)写为更加简洁的形式为：
Figure BDA00025422011500001413
考虑到式(6)和式(34)，将广义加速度写成相应的分块形式为：
Figure BDA00025422011500001414
因此，可将状态方程构造为：
Figure BDA00025422011500001415
Figure BDA00025422011500001416
对式(52)直接积分求得下一时刻的独立坐标和独立速度，其中初始时刻的广义坐标、速度和加速度分别由静力学分析、运动学分析和动力学分析得到。随后，将所求得的独立坐标和独立速度分别代入式(12)和式(36)，可得该时刻的相关坐标和相关速度。整合所求得的独立坐标和相关坐标即可得到该时刻的整体广义坐标，整合所求得的独立速度和相关速度即可得到该时刻的整体广义速度。此时，若运算时间大于所需要的终止时间，则可整合输出全部时刻的整体广义坐标、速度和加速度。若运算时间小于所需要的终止时间，则将该时刻的整体广义坐标代入运动学分析中，再次对系统进行运动学和动力学分析，直至运算时间大于所需要的终止时间。以第1根钻杆为例，利用MATLAB软件计算了其在下钻过程中的Z向位移、速度和加速度分别如图3、4和5所示。
由图可知，刚柔耦合模型中钻杆的Z向位移受钻杆重力的影响较为严重，而Z向速度和加速度受钻杆重力的影响相对较小。此外，三条曲线对于振动和摩擦的耦合影响均较为敏感，曲线的波动无明显规律，幅值的突变现象时有发生。由牛顿第二定律可知，力和加速度成正比，因此，由加速度曲线可以更直观的看出力的突变发生在哪一时刻，从而通过预先调整钻进参数来降低该时刻钻孔失效的风险，从而保证钻进过程中钻孔的稳定性。
对于本领域技术人员而言，显然能了解到上述具体事实例只是本发明的优选方案，因此本领域的技术人员对本发明中的某些部分所可能作出的改进、变动，体现的仍是本发明的原理，实现的仍是本发明的目的，均属于本发明所保护的范围。

Claims (4)

1.一种瓦斯抽采钻机系统的多体动力学分析方法，其特征在于，包括以下步骤：
S1、静力学分析：构建瓦斯抽采钻机系统的静力学方程，然后利用牛顿-拉夫逊Newton-Raphson法和差分法求出静力学方程的近似解，并通过静力学分析得到广义坐标；
S2、运动学分析：在静力学分析的前提下，构建瓦斯抽采钻机系统的运动学方程，然后利用代入法求解运动学方程，并通过运动学分析得到广义速度；
S3、动力学分析：构建瓦斯抽采钻机系统的动力学方程，利用直接积分和构造状态方程法求解动力学方程，并通过动力学分析得到每一时刻的广义坐标、广义速度、广义加速度，从而得到钻机钻进过程中的一系列动态响应。
2.根据权利要求1所述的瓦斯抽采钻机系统的多体动力学分析方法，其特征在于，所述步骤S1具体包括：
当系统静止时，构建的瓦斯抽采钻机系统的静力学方程：
Figure FDA0002542201140000011
设
Re＝Qe+Fmz-Kq (2)
将广义坐标写成分块形式：
Figure FDA0002542201140000012
qi、qd分别表示广义独立坐标和相关坐标，求约束方程对广义坐标的偏导，并写成紧凑形式为：
Figure FDA0002542201140000013
将其写为独立坐标、相关坐标分块形式为：
Figure FDA0002542201140000021
式中，
Figure FDA0002542201140000022
为nc×(n-nc)阶独立坐标约束的雅可比矩阵，
Figure FDA0002542201140000023
为nc×nc阶相关坐标约束的雅可比矩阵，nc表示相互独立的约束方程数，n表示广义坐标数，总体广义坐标的虚位移表示为：
Figure FDA0002542201140000024
Bdi为n×(n-nc)阶矩阵，将式(6)代入式(5)，两边乘以虚位移得：
Figure FDA0002542201140000025
求解式(7)即得到qi；
设
Figure FDA0002542201140000026
则有：
Figure FDA0002542201140000027
求解式(9)得到Δqi，其中
Figure FDA0002542201140000028
利用差分法求出，以此来确定全部独立坐标qi，同时，利用牛顿-拉夫逊Newton-Raphson法求解约束方程得到相关坐标qd，则总体广义坐标矢量为：
Figure FDA0002542201140000029
3.根据权利要求2所述的瓦斯抽采钻机系统的多体动力学分析方法，其特征在于，所述步骤S2具体包括：
将广义速度分块为：
Figure FDA00025422011400000210
求约束方程对时间的偏导，并写成分块形式为：
Figure FDA0002542201140000031
参考静力学分析的计算方法，参考(1)和(3)，由式(12)推导出相关速度与独立速度的关系为：
Figure FDA0002542201140000032
独立速度由静力学分析得出的独立坐标直接求导而得，将其代入式(6)求得相关速度，按照式(11)的形式整合独立速度和相关速度，即得到整个系统的初始广义速度。
4.根据权利要求3所述的瓦斯抽采钻机系统的多体动力学分析方法，其特征在于，所述步骤S3具体包括：
构建瓦斯抽采钻机系统的运动学方程：
Figure FDA0002542201140000033
式中，M和K分别为系统的质量矩阵和刚度矩阵；
Figure FDA0002542201140000034
和λ分别为约束的雅可比矩阵和拉格朗日乘子；Qe和Qv分别表示系统的广义外力矢量和速度的二次矢量；
Figure FDA0002542201140000035
为相互独立的约束方程组；为与广义坐标q和时间t相关的函数，nc表示约束方程的个数；
将式(13)对时间进行求导，即对约束方程组求二阶导数得：
Figure FDA0002542201140000037
式中，下标q表示雅可比矩阵，如
Figure FDA0002542201140000038
为
Figure FDA0002542201140000039
的雅可比矩阵；下标t表示对时间的导数，如
Figure FDA00025422011400000310
为
Figure FDA00025422011400000311
对时间的二阶偏导数；
Figure FDA00025422011400000312
和
Figure FDA00025422011400000313
分别为系统的广义速度和广义加速度；
式(15)整理为：
Figure FDA00025422011400000314
设：
Figure FDA0002542201140000041
则式(16)写为：
Figure FDA0002542201140000042
将式(14)中的动力学方程和式(18)整合成矩阵形式为：
Figure FDA0002542201140000043
式中：Ctt表示约束方程对时间的二次导数，Cq表示约束方程对广义坐标的一次导数；
设：
Figure FDA0002542201140000044
式(19)写成紧凑形式为：
Figure FDA0002542201140000045
求解式(21)得：
Figure FDA0002542201140000046
将式(22)写为分块形式为：
Figure FDA0002542201140000047
式中：
Figure FDA0002542201140000048
和
Figure FDA0002542201140000049
均为与广义坐标相关的函数，表示
Figure FDA00025422011400000410
的分量，其中，
Figure FDA00025422011400000411
为Mλ的逆矩阵；
Figure FDA00025422011400000412
和
Figure FDA00025422011400000413
均为与广义坐标、广义速度和时间相关的函数；
求解式(23)得加速度矢量和拉格朗日乘子为：
Figure FDA00025422011400000414
Figure FDA00025422011400000415
由式(24)知广义加速度与广义坐标、广义速度和时间相关，因此将式(24)写为更加简洁的形式为：
Figure FDA0002542201140000051
将广义加速度写成相应的分块形式为：
Figure FDA0002542201140000052
构造出状态方程：
Figure FDA0002542201140000053
Figure FDA0002542201140000054
对式(29)直接积分求得下一时刻的独立坐标和独立速度，其中初始时刻的广义坐标、速度和加速度分别由静力学分析、运动学分析和动力学分析得到；
随后，将所求得的独立坐标和独立速度分别代入式(6)和式(13)，得该时刻的相关坐标和相关速度；
整合所求得的独立坐标和相关坐标即得到该时刻的整体广义坐标，整合所求得的独立速度和相关速度即得到该时刻的整体广义速度；此时，若运算时间大于所需要的终止时间，则可整合输出全部时刻的整体广义坐标、速度和加速度；若运算时间小于所需要的终止时间，则将该时刻的整体广义坐标代入步骤S2中的运动学分析中，再次对系统进行运动学和动力学分析，直至运算时间大于所需要的终止时间。
CN202010550055.5A 2020-06-16 2020-06-16 一种瓦斯抽采钻机系统的多体动力学分析方法 Active CN111723536B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010550055.5A CN111723536B (zh) 2020-06-16 2020-06-16 一种瓦斯抽采钻机系统的多体动力学分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010550055.5A CN111723536B (zh) 2020-06-16 2020-06-16 一种瓦斯抽采钻机系统的多体动力学分析方法

Publications (2)

Publication Number Publication Date
CN111723536A true CN111723536A (zh) 2020-09-29
CN111723536B CN111723536B (zh) 2022-10-04

Family

ID=72567071

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010550055.5A Active CN111723536B (zh) 2020-06-16 2020-06-16 一种瓦斯抽采钻机系统的多体动力学分析方法

Country Status (1)

Country Link
CN (1) CN111723536B (zh)

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101639681A (zh) * 2008-07-29 2010-02-03 深圳市大族激光科技股份有限公司 一种电子装备运动机构性能参数优化方法
NL2007656C2 (en) * 2011-10-25 2013-05-01 Cofely Experts B V A method of and a device and an electronic controller for mitigating stick-slip oscillations in borehole equipment.
CN104239599A (zh) * 2014-07-07 2014-12-24 西安工业大学 一种基于多点定位柔性工装系统的动力学仿真分析方法
CN106695793A (zh) * 2017-01-18 2017-05-24 宁波韦尔德斯凯勒智能科技有限公司 一种xyzr四轴钻孔机器人主动柔顺控制装置及方法
CN107220421A (zh) * 2017-05-18 2017-09-29 北京理工大学 一种空间复杂柔性结构多体系统动力学建模与计算方法
CN109543264A (zh) * 2018-11-12 2019-03-29 天津理工大学 一种基于多维度重构校正的柔性多体机器人建模与求解方法
CN110259433A (zh) * 2019-06-28 2019-09-20 宝鸡石油机械有限责任公司 一种实体钻机数字化监测方法
CN111291499A (zh) * 2020-03-04 2020-06-16 岭南师范学院 一种基于多体动力学的瓦斯抽采钻机建模方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101639681A (zh) * 2008-07-29 2010-02-03 深圳市大族激光科技股份有限公司 一种电子装备运动机构性能参数优化方法
NL2007656C2 (en) * 2011-10-25 2013-05-01 Cofely Experts B V A method of and a device and an electronic controller for mitigating stick-slip oscillations in borehole equipment.
CN104239599A (zh) * 2014-07-07 2014-12-24 西安工业大学 一种基于多点定位柔性工装系统的动力学仿真分析方法
CN106695793A (zh) * 2017-01-18 2017-05-24 宁波韦尔德斯凯勒智能科技有限公司 一种xyzr四轴钻孔机器人主动柔顺控制装置及方法
CN107220421A (zh) * 2017-05-18 2017-09-29 北京理工大学 一种空间复杂柔性结构多体系统动力学建模与计算方法
CN109543264A (zh) * 2018-11-12 2019-03-29 天津理工大学 一种基于多维度重构校正的柔性多体机器人建模与求解方法
CN110259433A (zh) * 2019-06-28 2019-09-20 宝鸡石油机械有限责任公司 一种实体钻机数字化监测方法
CN111291499A (zh) * 2020-03-04 2020-06-16 岭南师范学院 一种基于多体动力学的瓦斯抽采钻机建模方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
ALEXANDER K. BELYAEV ET AL.: "Flexible rod model for the rotation of a drill string in an arbitrary borehole", 《ACTA MECHANICA》 *
WANG ZHE ET AL.: "Dynamics of flexible multibody systems with hybrid uncertain parameters", 《MECHANISM AND MACHINE THEORY》 *
朱吉良: "钻杆自动传送系统结构设计与仿真分析", 《中国优秀博硕士学位论文全文数据库(硕士)工程科技Ⅰ辑》 *
李康健: "旋挖钻机力学特性分析与疲劳寿命研究", 《中国优秀博硕士学位论文全文数据库(硕士)工程科技Ⅱ辑》 *
胡均平等: "螺旋钻机变幅时机液耦合动力学的键合图建模", 《中南大学学报(自然科学版)》 *
陈根良等: "基于广义坐标形式牛顿-欧拉方法的空间并联机构动力学正问题分析", 《机械工程学报》 *

Also Published As

Publication number Publication date
CN111723536B (zh) 2022-10-04

Similar Documents

Publication Publication Date Title
CN111136633B (zh) 针对时变时延下柔性主-从机器人系统的全状态控制方法
Johnson et al. A generalized particle algorithm for high velocity impact computations
CN105404304A (zh) 基于归一化神经网络的航天器容错姿态协同跟踪控制方法
Pappalardo et al. A comparative study of the principal methods for the analytical formulation and the numerical solution of the equations of motion of rigid multibody systems
CN112904728A (zh) 一种基于改进型趋近律的机械臂滑模控制轨迹跟踪方法
CN108628172A (zh) 一种基于扩张状态观测器的机械臂高精度运动控制方法
Piperno et al. Design of efficient partitioned procedures for the transient solution of aeroelastic problems
NO342719B1 (no) Proxymetoder for kostnadskrevende funksjonsoptimalisering med kostnadskrevende ikkelineare begrensninger
Fang et al. A POD reduced‐order 4D‐Var adaptive mesh ocean modelling approach
CN109614580A (zh) 基于在线Xgboost算法的抗震混合试验模型更新方法
Gao et al. Modeling and error compensation of robotic articulated arm coordinate measuring machines using BP neural network
Zahr et al. An adjoint method for a high-order discretization of deforming domain conservation laws for optimization of flow problems
Love et al. An energy‐consistent material‐point method for dynamic finite deformation plasticity
Nejat et al. Adjoint sensitivity analysis of flexible multibody systems in differential-algebraic form
Ding et al. Second order adjoint sensitivity analysis of multibody systems described by differential–algebraic equations
CN108445768A (zh) 空间机器人操作空间轨迹跟踪的增广自适应模糊控制方法
CN111291499B (zh) 一种基于多体动力学的瓦斯抽采钻机建模方法
CN114707385A (zh) 基于Krylov子空间的深部地层导热系数三维预测方法及装置
Farhat et al. Transient aeroelastic computations using multiple moving frames of reference
CN111723536B (zh) 一种瓦斯抽采钻机系统的多体动力学分析方法
CN113276114B (zh) 一种基于终端任务指派的可重构机械臂协同力/运动控制系统与方法
Duruisseaux et al. Lie group forced variational integrator networks for learning and control of robot systems
De Silva et al. The right invariant nonlinear complementary filter for low cost attitude and heading estimation of platforms
CN104763694A (zh) 一种掘进机液压推进系统分区压力设定值优化方法
CN111241728A (zh) 一种欧拉方程的间断伽辽金有限元数值求解方法

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