BIFURCATION AND BASINS OF ATTRACTION IN A FOUR-DEGREE-OF-FREEDOM DRIFTER NONLINEAR ROCK MODEL
-
摘要:
将凿岩机钻进岩石的过程,建立成四自由度凿岩机-非线性岩石模型。利用直接数值积分方法,将液压作用力的角频率和振幅作为控制参数,进行分岔和吸引域分析。结果表明:当ω<2.08时,模型为粘滞模式。当2.08≤ω≤6.5时,模型为非粘滞模式。当4.94≤ω≤5.36时,周期-3解表现出了比周期-1解和准周期解更大的全局稳定性。当0<a ≤1时,随着振幅的增加,周期-1-2解的吸引域逐渐增加,周期-3解的吸引域逐渐减小,周期-1-1解的吸引域起初增加,然后变为0。比较了从模型和实验结果获得的活塞位移和速度数据,实验和模型之间存在很强的相关性。
Abstract:The drilling process into rock has been modeled as a four-degree-of-freedom rock drill with a nonlinear rock model. Using direct numerical integration, the angular frequency and amplitude of the hydraulic force were employed as control parameters for bifurcation and basins of attraction analysis. The results showed that, when ω<2.08, the model operated in the stick mode; when 2.08≤ω≤6.5, the model operated in the non-stick mode. The period-3 solution exhibited greater global stability than periodic-1 and quasi-period solutions when 4.94≤ω≤5.36. When 0< a≤1, the stability of the periodic-1-2 solution increased with the amplitude, while the stability of the periodic-3 solution decreased. The basins of attraction of the period-1-1 solution increased at first and then became zero. A comparison of piston displacement and velocity data between the model and the experiment results revealed a strong correlation.
-
Keywords:
- hydraulic drifter /
- nonlinear rock model /
- stick /
- bifurcation /
- basins of attraction
-
液压凿岩机是凿岩台车、锚杆台车、露天钻机、硬岩掘进机(TBM)等特种工程机械装备上实现凿岩钻孔功能的核心部件,被广泛应用于矿山开采、交通建设、水利水电等领域[1 − 2]。学者们已进行了多项研究以分析和改进液压凿岩机的性能。
Seo[3]对凿岩机的冲击性能进行了分析,采用仿真建模和实验验证来增强其冲击效果。郭勇等[4]提出了一个高频液压凿岩机带套分配器的冲击系统液压机械耦合模型。杨书仪等[5]研究了凿岩机活塞腔的压力变化特性曲线,并探讨了流量、充压压力和溢流阀设置对冲击性能和效率的影响。耿晓光等[6, 7]研究了液压凿岩机双缓冲系统的活塞的汽蚀现象及其影响因素,提出了减轻汽蚀现象的解决方案。杨桢毅等[8]研究了液压凿岩机冲击过程中活塞的振动行为,建立了动态模型来分析一个冲击周期内仅存在一次初级振动的范围。李叶林等[9]分析了缓冲机构对冲击能量、频率和功率的影响,并提出了优化冲击性能的方法。Shah[10]提出了在凿岩机中使用电动阀改善可控性、性能和操作范围的可能性。
这些参考文献涵盖了一系列集中在液压凿岩机分析、建模和优化方面的研究。然而,他们在研究中忽略了将液压凿岩机与岩石整合在一起。
在冲击钻进领域中,Hashiba等[11]专注于通过实验测试和数值模拟对钻头在冲击进入岩石过程中的力-穿透曲线进行建模。Lundberg和Collet[12]研究了使用整体钻杆和可拆卸钻头的冲击钻进波形优化。他们发现通过优化波形可以实现更高的效率,特别是对于更长的波形。廖茂林等[13]对钻头-岩石振动冲击系统进行了研究,以提高钻进效率。
岩石钻进是一个复杂的过程,需要液压凿岩机的高效稳定性能。先前的研究侧重于凿岩机本身,但在凿岩机钻进岩石过程中的稳定性方面缺乏研究。
1 力学模型
1.1 岩石接触模型
本文采用Ajibose提出的一种描述钻头与岩石相互作用的方法[14]。构建的岩石接触模型如下:
kl(X∗n−Xpn)nl=ku(X∗n−Xfn)nu, (1) 式中:kl为加载刚度;ku为卸载刚度;Xpn、Xfn和X∗n分别为n次冲击周期中岩石表面位置的初始加载位置、最终卸载位置和最大钻进点。
图1为非线性岩石接触模型的力-钻进位移关系。图1中标注了加载和卸载阶段。
对于图1(a)所示的周期-1单次冲击,通过式(1)可以预测卸载阶段的结束位置Xfn为:
Xfn=X∗n−[kl(X∗n−Xpn)nlku]1/nu (2) 式中,n=1,2。此时,前一个冲击周期卸载阶段的结束位置,实际上也是下一个冲击周期加载阶段的起始位置。
另一种情况是多次冲击,相应的力-钻进位移关系如图1(b)所示。在卸载过程中,达到预测值Xf1之前,速度的符号已经发生了改变。此时需要根据速度方向改变的位置Xc来计算虚拟起始位置Xp2:
Xp2=Xc−[ku(Xc−Xf1)nukl]1/nl (3) 在这种情况下,第二个冲击周期的虚拟起始位置Xp2与第一个冲击周期的预测结束位置Xf1不重合[15]。
因此,在接下来的动态研究中,当发生单次冲击时,使用式(2);而发生多次冲击时,则使用式(3)。
根据Ajibose的研究,加载指数nl和卸载指数nu与镶嵌在钻头上的齿形相关[14]。如果式(1)的岩石接触模型是线性的,则加载指数和卸载指数都为1,即nl=nu=1,图1中的加载和卸载阶段将变为直线;本文采用非线性拟合,钻头上为圆锥齿,采用nl=nu=2。
图2为四自由度凿岩机-岩石的力学模型,描述了将凿岩机冲击系统的作业过程简化为活塞撞击钻具,能量通过钻具传递到岩石的过程。机械振动的自由度定义为系统中可以独立运动的部件的数量[16, 17]。图2中可以独立运动的部件是冲击活塞、钻具、钻头和岩石表面,其绝对位移用Xi(i=1,2,3,4)表示,分别标注在每个部件的左侧。
活塞m1由简谐液压作用力F(t)=Fasin(Ωt)+Fb驱动,包含振幅Fa、角频率Ω和垂直偏移量Fb。液压油具有由刚度k1和阻尼系数c1表征的线性粘弹性。钻具由钎尾、连接套和钻杆组成,总质量为m2。钻头的质量为m3。岩石通过刚度k2和阻尼系数c2的线性弹簧连接到钻杆上。在t=0时,活塞和钻具之间存在一个间隙G。物理量(X1−X2−G)用于监测该距离。假设模型在水平面上运行,忽略重力的作用。表1为凿岩机-岩石模型的量纲参数,与自主研发的液压凿岩机和凿岩台车中提供的液压参数相对应。
表 1 凿岩机-岩石模型的量纲参数Table 1. Dimension parameters in the drifter-rock model含义 数值 液压力的角频率Ω/(rad/s) 240 液压油阻尼c1/(N/(m/s)) 1×103 结构阻尼c2/(N/(m/s)) 2.45×103 液压力的振幅Fa/N 1.1417×104 液压力的垂直偏移量Fb/N 6.7154×104 活塞和钻具的距离G/m 0.025 液压油刚度k1/(N/m) 2×105 结构刚度k2/(N/m) 6.885×104 加载刚度kl/(N/m) 11.46×106 卸载刚度ku/(N/m) 13.06×106 活塞质量m1/kg 11.5 钻具质量m2/kg 32.6 钻头质量m3/kg 1 模型在六种模式中运行:“非粘滞-非钻进”、“非粘滞-加载钻进”、“非粘滞-卸载钻进”、“粘滞-非钻进”、“粘滞-加载钻进”和“粘滞-卸载钻进”。以下为每个模式的简要概述。
1.2 非粘滞-非钻进
在t=0时,假设活塞和钻具之间的初始距离为间隙G。当X1−X2<G时,活塞和钻具没有碰撞。活塞、钻具和钻头的动力学方程为:
{m1¨X1=Fasin(Ωt)+Fb−c1˙X1−k1X1,m2¨X2=−c2(˙X2−˙X3)−k2(X2−X3),m3¨X3=c2(˙X2−˙X3)+k2(X2−X3). (4) 当X1−X2=G时,活塞撞击钻具。如果碰撞前的瞬间,活塞和钻具的速度差值比较大,那么碰撞后两者就会分离,称为“非粘滞”模式。在这种情况下,动力学方程可以通过动量守恒推导出来。活塞和钻具之间的碰撞关系为:
{X1+=X2−+G,˙X1+=m1−m2em1+m2˙X1−+(1+e)m2m1+m2˙X2−,X2+=X2−,˙X2+=(1+e)m1m1+m2˙X1−+m2−m1em1+m2˙X2−. (5) 式中:下角标(-)和(+)分别为碰撞前和碰撞后;e为活塞和钻具之间的碰撞恢复系数[18 − 19]。
“非钻进”模式指的是钻头不与岩石表面接触。因此,当钻头的位移小于岩石表面的初始位置(X3−Xp<0)时,模型处于这种模式。
当没有钻进时,岩石表面保持静止,其速度为零。
˙X4=0 (6) 在非钻进模式中,Xp和Xf都保持不变:
{˙Xp=0,˙Xf=0. (7) 结合式(4)、式(6)和式(7),得到“非粘滞-非钻进”模式的方程为:
{m1¨X1=Fasin(Ωt)+Fb−c1˙X1−k1X1,m2¨X2=−c2(˙X2−˙X3)−k2(X2−X3),m3¨X3=c2(˙X2−˙X3)+k2(X2−X3),˙X4=0,˙Xp=0,˙Xf=0. (8) 1.3 非粘滞-加载钻进
当钻头与岩石表面接触并向前压缩时,为“加载钻进”模式。这种模式需要同时满足两个条件:m3的位移应大于或等于岩石表面的初始位置(X3−Xp⩾);{m_3}的速度应为正值({\dot X_3} {\geqslant} 0)。
活塞、钻具和钻头的动力学方程为:
\left\{ \begin{aligned} & {{m_1}{{\ddot X}_1} = {F_a}\sin (\varOmega t) + {F_b} - {c_1}{{\dot X}_1} - {k_1}{X_1},} \\ & {{m_2}{{\ddot X}_2} = - {c_2}( {{{\dot X}_2} - {{\dot X}_3}} ) - {k_2}( {{X_2} - {X_3}} ),} \\ & {m_3}{{\ddot X}_3} = {c_2}( {{{\dot X}_2} - {{\dot X}_3}} ) + {k_2}( {{X_2} - {X_3}} ) - \\ &\qquad\;\quad{k_l}{( {{X_3} - {X_p}} )^{{n_l}}}. \\ \end{aligned} \right. (9) 在钻进过程中,钻头和岩石表面的速度相等,即:
{\dot X_4} = {\dot X_3} (10) 在加载钻进模式中,{X_p}的值已确认且保持恒定,而{X_f}的值未知。在准备加载钻进模式时,需要不断计算{X_f}的值,直到{m_3}达到最大钻进位置X_3^*。参考式(1)提出的概念,控制加载钻进模式的微分方程式可表示为:
\left\{ \begin{aligned} &{{{\dot X}_p} = 0,} \\ &{{{\dot X}_f} = \left[ {1 - \dfrac{{{k_l}{n_l}{{( {{X_3} - {X_p}} )}^{{n_l} - 1}}}}{{{k_u}{n_u}{{( {{X_3} - {X_f}} )}^{{n_u} - 1}}}}} \right]{{\dot X}_3}.} \end{aligned} \right. (11) 结合式(9)、式(10)和式(11),得到“非粘滞-加载钻进”模式的方程为:
\left\{ {\begin{aligned} & {{m_1}{{\ddot X}_1} = {F_a}\sin (\varOmega t) + {F_b} - {c_1}{{\dot X}_1} - {k_1}{X_1},} \\[-1.5pt] & {{m_2}{{\ddot X}_2} = - {c_2}( {{{\dot X}_2} - {{\dot X}_3}} ) - {k_2}( {{X_2} - {X_3}} ),} \\[-1.5pt] & {m_3}{{\ddot X}_3} = {c_2}( {{{\dot X}_2} - {{\dot X}_3}} ) + {k_2}( {{X_2} - {X_3}} ) - \\[-1.5pt] &\qquad\quad\;{k_l}{( {{X_3} - {X_p}} )^{{n_l}}}, \\[-1.5pt] & {{{\dot X}_4} = {{\dot X}_3},} \\[-1.5pt] & {{{\dot X}_p} = 0,} \\[-1.5pt] &{{{\dot X}_f} = \left[ {1 - \dfrac{{{k_l}{n_l}{{( {{X_3} - {X_p}} )}^{{n_l} - 1}}}}{{{k_u}{n_u}{{({{X_3} - {X_f}} )}^{{n_u} - 1}}}}} \right]{{\dot X}_3}.} \end{aligned}} \right. (12) 参考式(2),可以在加载钻进模式结束时预测卸载钻进模式的终点{X_f}:
X_f=X_3^*-\left[\frac{k_l(X_3^*-X_p)^{n_l}}{k_u}\right]^{1 / n_u}, (13) 式中,X_3^*为加载钻进模式结束时的位移,表示达到最大钻进的位置。
1.4 非粘滞-卸载钻进
在“非粘滞-卸载钻进”模式中,钻头仍然与岩石表面保持接触,但岩石表面向后移动直到达到终点{X_f}。这种情况必须同时满足三个条件:{m_3}的位移应大于或等于岩石表面的初始位置({X_3} - {X_p} {\geqslant} 0);{m_3}的速度应为负值({\dot X_3} < 0);{m_3}的位移应大于或等于卸载钻进模式的终点位置({X_3} - {X_f} {\geqslant} 0)。一旦满足这三个条件,控制活塞、钻具、钻头和岩石表面的微分方程可以描述如下:
\left\{ { \begin{aligned} &{{m_1}{{\ddot X}_1} = {F_a}\sin (\varOmega t) + {F_b} - {c_1}{{\dot X}_1} - {k_1}{X_1},} \\ &{{m_2}{{\ddot X}_2} = - {c_2}( {{{\dot X}_2} - {{\dot X}_3}} ) - {k_2}( {{X_2} - {X_3}} ),} \\ & \begin{gathered} {m_3}{{\ddot X}_3} = {c_2}( {{{\dot X}_2} - {{\dot X}_3}} ) + {k_2}( {{X_2} - {X_3}} ) - \\ \quad\quad\quad\;{k_u}{( {{X_3} - {X_f}} )^{{n_u}}}, \\ \end{gathered} \\ & {{{\dot X}_4} = {{\dot X}_3}.} \end{aligned}} \right. (14) 在“卸载钻进”模式中,{X_f}的值在前一加载钻进模式中确定并保持恒定。然而,在即将到来的冲击期间,{X_p}的新值是未知的,需要持续重新计算,直到{m_3}达到预测的最终位置{X_f}或速度方向改变的位置{X_{3c}}。根据式(1)中描述的基本概念,控制卸载钻进模式中{X_p}和{X_f}的微分方程可以表示如下:
\left\{ \begin{aligned} &{{{\dot X}_p} = \left( {1 - \dfrac{{{k_u}{n_u}{{( {{X_3} - {X_f}} )}^{{n_u} - 1}}}}{{{k_l}{n_l}{{( {{X_3} - {X_p}} )}^{{n_l} - 1}}}}} \right){{\dot X}_3},} \\ & {{{\dot X}_f} = 0.} \end{aligned} \right. (15) 结合式(14)和式(15),得到“非粘滞-卸载钻进”模式的方程为:
\left\{ { \begin{aligned} &{{m_1}{{\ddot X}_1} = {F_a}\sin (\varOmega t) + {F_b} - {c_1}{{\dot X}_1} - {k_1}{X_1},} \\ &{{m_2}{{\ddot X}_2} = - {c_2}( {{{\dot X}_2} - {{\dot X}_3}} ) - {k_2}( {{X_2} - {X_3}} ),} \\ &\begin{gathered} {m_3}{{\ddot X}_3} = {c_2}( {{{\dot X}_2} - {{\dot X}_3}} ) + {k_2}( {{X_2} - {X_3}} ) - \\ \quad\quad\quad\;{k_u}{( {{X_3} - {X_f}} )^{{n_u}}}, \\ \end{gathered} \\ &{{{\dot X}_4} = {{\dot X}_3},} \\ & {{{\dot X}_p} = \left[ {1 - \dfrac{{{k_u}{n_u}{{( {{X_3} - {X_f}} )}^{{n_u} - 1}}}}{{{k_l}{n_l}{{( {{X_3} - {X_p}} )}^{{n_l} - 1}}}}} \right]{{\dot X}_3},} \\ & {{{\dot X}_f} = 0.} \end{aligned} } \right. (16) 此外,卸载钻进模式的终点位置可以确定后续加载钻进模式的新初始位置。这个新的初始位置等同于当前卸载钻进模式的终点位置,表示为{X_p} = {X_f}。
在这种情况下,模型将直接返回加载钻进模式,而不是恢复到非钻进模式。这是终止卸载钻进模式的另一种方法,特别是当{m_3}的速度在到达位置{X_f}之前改变方向时。在这种情况下,模型的行为再次由式(12)控制,但增加了两个额外的条件:{X_3} - {X_p} {\geqslant} 0和{\dot X_3} {\geqslant} 0。参考式(3),直接连接的加载模式的虚拟初始位置可以计算如下:
X_p=X_{3 c}-\left[\frac{k_u\left(X_{3 c}-X_f\right)^{n_u}}{k_l}\right]^{1 / n_l}, (17) 式中,{X_{3c}}为{m_3}速度变为零的位置。
1.5 粘滞-非钻进
图3为粘滞情况下的四自由度凿岩机-岩石的力学模型。当活塞和钻具碰撞前的瞬间,两者速度接近时,碰撞后会一起移动,称为“粘滞”现象。值得注意的是,凿岩机的粘滞仅存在于理论模型,因为实际工况中,碰撞前的瞬间,活塞和钻具的速度差值很难达到粘滞条件。
当活塞保持与钻具接触并同时移动时,就会出现“粘滞”现象。活塞、钻具和钻头的运动控制方程为:
\left\{ \begin{gathered} ( {{m_1} + {m_2}} ){{\ddot X}_1} = {F_a}\sin (\varOmega t) + {F_b} - {c_1}{{\dot X}_1} - \\ \quad{c_2}( {{{\dot X}_2} - {{\dot X}_3}} ) - {k_1}{X_1} - {k_2}( {{X_2} - {X_3}} ), \\ ( {{m_1} + {m_2}} ){{\ddot X}_2} = {F_a}\sin (\varOmega t) + {F_b} - {c_1}{{\dot X}_1} - \\ \quad{c_2}( {{{\dot X}_2} - {{\dot X}_3}} ) - {k_1}{X_1} - {k_2}( {{X_2} - {X_3}} ), \\ {m_3}{{\ddot X}_3} = {c_2}( {{{\dot X}_2} - {{\dot X}_3}} ) + {k_2}( {{X_2} - {X_3}} ). \\ \end{gathered} \right. (18) 粘滞条件为AC{C_1} {\geqslant} AC{C_2}、{X_1} = {X_2} + G和{\dot X_1} = {\dot X_2}。其中,AC{C_1}定义为非粘滞模式中{m_1}的加速度:
AC{C_1} = [ {{F_a}\sin (\varOmega t) + {F_b} - {c_1}{{\dot X}_1} - {k_1}{X_1}} ]/{m_1} (19) AC{C_2}定义为粘滞模式中({m_1} + {m_2})的加速度:
\begin{split} AC{C_2} =& [{F_a}\sin (\varOmega t) + {F_b} - {c_1}{{\dot X}_1} - {c_2}( {{{\dot X}_2} - {{\dot X}_3}} ) - \\& {k_1}{X_1} - {k_2}( {{X_2} - {X_3}} )]/( {{m_1} + {m_2}} ) \end{split} (20) 结合式(6)、式(7)和式(18),得到“粘滞-非钻进”模式的方程为:
\left\{ {\begin{aligned} &\begin{gathered} ( {{m_1} + {m_2}} ){{\ddot X}_1} = {F_a}\sin (\varOmega t) + {F_b} - {c_1}{{\dot X}_1} - \\ \quad{c_2}( {{{\dot X}_2} - {{\dot X}_3}} ) - {k_1}{X_1} - {k_2}( {{X_2} - {X_3}} ), \\ \end{gathered} \\ & \begin{gathered} ( {{m_1} + {m_2}} ){{\ddot X}_2} = {F_a}\sin (\varOmega t) + {F_b} - {c_1}{{\dot X}_1} - \\ \quad{c_2}( {{{\dot X}_2} - {{\dot X}_3}} ) - {k_1}{X_1} - {k_2}( {{X_2} - {X_3}} ), \\ \end{gathered} \\ & {{m_3}{{\ddot X}_3} = {c_2}( {{{\dot X}_2} - {{\dot X}_3}} ) + {k_2}( {{X_2} - {X_3}} ),} \\ & {{{\dot X}_4} = 0,} \\ &{{{\dot X}_p} = 0,} \\ & {{{\dot X}_f} = 0.} \end{aligned}} \right. (21) 1.6 粘滞-加载钻进
当钻头与岩石表面接触向前压缩,并且活塞和钻具一起运动时,为“粘滞-加载钻进”模式。结合式(10)、式(11)和式(18),得到“粘滞-加载钻进”模式的方程为:
\left\{ {\begin{aligned} &\begin{gathered} ( {{m_1} + {m_2}} ){{\ddot X}_1} = {F_a}\sin (\varOmega t) + {F_b} - {c_1}{{\dot X}_1} - \\ \quad{c_2}( {{{\dot X}_2} - {{\dot X}_3}} ) - {k_1}{X_1} - {k_2}( {{X_2} - {X_3}} ), \\ \end{gathered} \\ &\begin{gathered} ( {{m_1} + {m_2}} ){{\ddot X}_2} = {F_a}\sin (\varOmega t) + {F_b} - {c_1}{{\dot X}_1} - \\ \quad{c_2}( {{{\dot X}_2} - {{\dot X}_3}} ) - {k_1}{X_1} - {k_2}( {{X_2} - {X_3}} ), \\ \end{gathered} \\ &{{m_3}{{\ddot X}_3} = {c_2}( {{{\dot X}_2} - {{\dot X}_3}} ) + {k_2}( {{X_2} - {X_3}} ),} \\ & {{{\dot X}_4} = {{\dot X}_3},} \\ & {{{\dot X}_p} = 0,} \\ &{{{\dot X}_f} = \left[ {1 - \dfrac{{{k_l}{n_l}{{( {{X_3} - {X_p}} )}^{{n_l} - 1}}}}{{{k_u}{n_u}{{({{X_3} - {X_f}} )}^{{n_u} - 1}}}}} \right]{{\dot X}_3}.} \end{aligned}} \right. (22) 1.7 粘滞-卸载钻进
在“粘滞-卸载钻进”模式中,钻头仍与岩石表面保持接触,岩石表面向后移动直到达到终点{X_f}。在这种模式下,活塞与钻具一起运动。结合式(10)、式(15)和式(18),得到“粘滞-卸载钻进”模式的方程为:
\left\{ {\begin{aligned} &\begin{gathered} ( {{m_1} + {m_2}} ){{\ddot X}_1} = {F_a}\sin (\varOmega t) + {F_b} - {c_1}{{\dot X}_1} - \\ \quad{c_2}( {{{\dot X}_2} - {{\dot X}_3}} ) - {k_1}{X_1} - {k_2}( {{X_2} - {X_3}} ), \\ \end{gathered} \\ & \begin{gathered} ( {{m_1} + {m_2}} ){{\ddot X}_2} = {F_a}\sin (\varOmega t) + {F_b} - {c_1}{{\dot X}_1} - \\ \quad{c_2}( {{{\dot X}_2} - {{\dot X}_3}} ) - {k_1}{X_1} - {k_2}( {{X_2} - {X_3}} ), \\ \end{gathered} \\ & \begin{gathered} {m_3}{{\ddot X}_3} = {c_2}( {{{\dot X}_2} - {{\dot X}_3}} ) + {k_2}( {{X_2} - {X_3}} ) - \\ \quad{k_u}{( {{X_3} - {X_f}} )^{{n_u}}}, \\ \end{gathered} \\ & {{{\dot X}_4} = {{\dot X}_3},} \\ & {{{\dot X}_p} = \left[ {1 - \dfrac{{{k_u}{n_u}{{( {{X_3} - {X_f}} )}^{{n_u} - 1}}}}{{{k_l}{n_l}{{( {{X_3} - {X_p}} )}^{{n_l} - 1}}}}} \right]{{\dot X}_3},} \\ & {{{\dot X}_f} = 0.} \end{aligned}} \right. (23) 2 无量纲化
凿岩机-岩石模型的无量纲参数用\eta : = (\omega ,a,b, \alpha , \beta ,\gamma ,\lambda ,\kappa ,\mu ,\zeta ,g,{n_l},{n_u},e)表示。具体含义和数值见表2。
凿岩机-岩石模型的无量纲状态变量用u: = \left( {y_1}, {y_2}, {y_3},{y_4},{y_5},{y_6},{y_7},{y_8},{y_9},{y_{10}} \right)^{\text{T}}表示。{y_1}为活塞位移,{y_2}为活塞速度,{y_3}为钻具位移,{y_4}为钻具速度,{y_5}为钻头位移,{y_6}为钻头速度,{y_7}为钻进位移,{y_1}\sim{y_7} \in ( - \infty , + \infty );{y_8}为岩石表面初始位置,{y_9}为岩石表面终点位置,{y_8},{y_9} \in [0, + \infty );{y_{10}}为相位,{y_{10}} \in [0,2\pi )。{y'_1}\sim{y'_{10}}分别为变量{y_1}\sim{y_{10}}关于无量纲时间\tau 的微分。
结合式(8)、式(12)和式(16),经过无量纲化,得到非粘滞模式的紧凑形式:
\begin{gathered} {\boldsymbol{u}}' = {{\boldsymbol{f}}_{{\text{NS}}}}(\tau ,{\boldsymbol{u}},{\boldsymbol{\eta}} ) = \\ \left\{ \begin{gathered} {y_2}, \\ \frac{1}{\alpha }\left[ {a\sin \left( {{y_{10}}} \right) + b - \gamma 2\zeta {y_2} - \beta {y_1}} \right], \\ {y_4}, \\ - 2\zeta \left( {{y_4} - {y_6}} \right) - \left( {{y_3} - {y_5}} \right), \\ {y_6}, \\ \mu [2\zeta \left( {{y_4} - {y_6}} \right) + \left( {{y_3} - {y_5}} \right) - {h_1}{h_2}\kappa {({y_5} - {y_8})^{{n_l}}} - \\ \quad{h_1}\left( {1 - {h_2}} \right){h_3}\lambda {({y_5} - {y_9})^{{n_u}}}], \\ \left[ {{h_1}{h_2} + {h_1}\left( {1 - {h_2}} \right){h_3}} \right]{y_6}, \\ {h_1}\left( {1 - {h_2}} \right){h_3}\left[ {1 - \frac{{\lambda {n_u}{{({y_5} - {y_9})}^{{n_u} - 1}}}}{{\kappa {n_l}{{({y_5} - {y_8})}^{{n_l} - 1}}}}} \right]{y_6}, \\ {h_1}{h_2}\left[ {1 - \frac{{\kappa {n_l}{{({y_5} - {y_8})}^{{n_l} - 1}}}}{{\lambda {n_u}{{({y_5} - {y_9})}^{{n_u} - 1}}}}} \right]{y_6}, \\ \omega . \\ \end{gathered} \right. \\ \end{gathered} (24) 结合式(21)、式(22)和式(23),经过无量纲化,得到粘滞模式的紧凑形式:
\begin{gathered} {\boldsymbol{u}}' = {{\boldsymbol{f}}_{\text{S}}}(\tau ,{\boldsymbol{u}},{\boldsymbol{\eta}} ) = \\ \left\{ \begin{gathered} {y_2}, \\ \frac{1}{{(\alpha + 1)}}[a\sin \left( {{y_{10}}} \right) + b + c - (\gamma + 1)2\zeta {y_2} - (\beta + 1){y_1} + {y_5}], \\ {y_2}, \\ \frac{1}{{(\alpha + 1)}}[a\sin \left( {{y_{10}}} \right) + b + c - (\gamma + 1)2\zeta {y_2} - (\beta + 1){y_1} + {y_5}], \\ {y_6}, \\ \mu [2\zeta \left( {{y_4} - {y_6}} \right) + \left( {{y_3} - {y_5}} \right) - {h_1}{h_2}\kappa {({y_5} - {y_8})^{{n_l}}} - \\ \quad{h_1}\left( {1 - {h_2}} \right){h_3}\lambda {({y_5} - {y_9})^{{n_u}}}], \\ \left[ {{h_1}{h_2} + {h_1}\left( {1 - {h_2}} \right){h_3}} \right]{y_6}, \\ {h_1}\left( {1 - {h_2}} \right){h_3}\left[ {1 - \frac{{\lambda {n_u}{{({y_5} - {y_9})}^{{n_u} - 1}}}}{{\kappa {n_l}{{({y_5} - {y_8})}^{{n_l} - 1}}}}} \right]{y_6}, \\ {h_1}{h_2}\left[ {1 - \frac{{\kappa {n_l}{{({y_5} - {y_8})}^{{n_l} - 1}}}}{{\lambda {n_u}{{({y_5} - {y_9})}^{{n_u} - 1}}}}} \right]{y_6}, \\ \omega . \\ \end{gathered} \right. \\ \end{gathered} (25) 在式(24)和式(25)中, {h}_{1}、{h}_{2} 和{h_3}为单位阶跃函数,定义如下:
{h_1} = \left\{ \begin{gathered} 1,{\text{ }}{y_5} - {y_8} {\geqslant} 0 \\ 0,{\text{ }}{y_5} - {y_8} < 0 \\ \end{gathered} \right. (26) {h_2} = \left\{ \begin{gathered} 1,{\text{ }}{y_6} {\geqslant} 0 \\ 0,{\text{ }}{y_6} < 0 \end{gathered} \right.\quad \quad (27) {h_3} = \left\{ \begin{gathered} 1,{\text{ }}{y_5} - {y_9} {\geqslant} 0 \\ 0,{\text{ }}{y_5} - {y_9} < 0 \\ \end{gathered} \right. (28) 在式(24)和式(25)中,变量和参数如下所示:
\begin{split} {y_1} &= \frac{{{X_1}}}{{{x_s}}},{y_2} = {{y'}_1},{y_3} = \frac{{{X_2}}}{{{x_s}}},{y_4} = {{y'}_3},{y_5} = \frac{{{X_3}}}{{{x_s}}},{y_6} = {{y'}_5}, \\ {y_7} &= \frac{{{X_4}}}{{{x_s}}},{y_8} = \frac{{{X_p}}}{{{x_s}}},{y_9} = \frac{{{X_f}}}{{{x_s}}},{y_{10}} = \omega \tau ,g = \frac{G}{{{x_s}}},\tau = \frac{t}{{{t_s}}}, \\ \zeta &= \frac{{{c_2}}}{{2\sqrt {{m_2}{k_2}} }},{\varOmega _n} = \sqrt {\frac{{{k_2}}}{{{m_2}}}} ,\omega = \frac{\varOmega }{{{\varOmega _n}}},\alpha = \frac{{{m_1}}}{{{m_2}}},\beta = \frac{{{k_1}}}{{{k_2}}}, \\ \gamma &= \frac{{{c_1}}}{{{c_2}}},\mu = \frac{{{m_2}}}{{{m_3}}},a = \frac{{{F_a}}}{{{k_2}{x_s}}},b = \frac{{{F_b}}}{{{k_2}{x_s}}},\kappa = \frac{{{k_l}}}{{{k_2}}}x_s^{{n_l} - 1}, \\ \lambda &= \frac{{{k_u}}}{{{k_2}}}x_s^{{n_u} - 1} \\ \end{split} 式中:{x_s}为参考位移({x_s} = 1 m);{t_s}为参考时间({t_s} = \sqrt {{m_2}/{k_2}} )。
结合式(24)和式(25),凿岩机-岩石模型的向量场方程以紧凑形式表达如下:
{\boldsymbol{u}}' = \left\{ {\begin{array}{*{20}{l}} {{{\boldsymbol{f}}_{{\text{NS}}}}(\tau ,{\boldsymbol{u}},{\boldsymbol{\eta}} ),{\text{ }}ac{c_1} < ac{c_2}} \\ {{{\boldsymbol{f}}_{\text{S}}}(\tau ,{\boldsymbol{u}},{\boldsymbol{\eta}} ),{\text{ }}ac{c_1} {\geqslant} ac{c_2}} \end{array}} \right. (29) 式中,ac{c_1}为非粘滞模式下{m_1}的加速度,由式(19)经过无量纲化后得到:
ac{c_1} = \left[ {a\sin \left( {{y_{10}}} \right) + b - 2\zeta \gamma {y_2} - \beta {y_1}} \right]/\alpha (30) ac{c_2}为非粘滞模式下\left( {{m_1} + {m_2}} \right)的加速度,由式(20)经过无量纲化后得到:
\begin{split} ac{c_2} = &[a\sin \left( {{y_{10}}} \right) + b - (\gamma + 1)2\zeta {y_2} - \\& (\beta + 1){y_1} + {y_5}]/(\alpha + 1) \end{split} (31) 由式(5)经过无量纲化,得到活塞和钻具之间的碰撞关系为:
\left\{ \begin{gathered} {y_{1 + }} = {y_{3 - }} + g, \\ {y_{2 + }} = \frac{{(\alpha - e){y_{2 - }} + (1 + e){y_{4 - }}}}{{\alpha + 1}}, \\ {y_{3 + }} = {y_{3 - }}, \\ {y_{4 + }} = \frac{{(1 + e)\alpha {y_{2 - }} + (1 - \alpha e){y_{4 - }}}}{{\alpha + 1}}, \\ \end{gathered} \right. (32) 表 2 凿岩机-岩石模型的无量纲参数Table 2. Dimensionless parameters in the drifter-rock model序号 符号 含义 数值 1 \omega 液压作用力的角频率 (0,6.5] 2 a 液压作用力的振幅 (0,1] 3 b 液压作用力的垂直偏移量 0.0975 4 \alpha 活塞和钻具的质量比 0.35 5 \beta 刚度比 2.9 6 \gamma 阻尼系数比 0.4 7 \lambda 卸载刚度和结构刚度的比值 189.7 8 \kappa 加载刚度和结构刚度的比值 166.4 9 \mu 钻具和钻头的质量比 32.6 10 \zeta 阻尼比 0.8 11 g 活塞和钻具的间隙 0.025 12 {n_l} 加载指数 2 13 {n_u} 卸载指数 2 14 e 碰撞恢复系数 0.3 3 分岔和吸引域分析
3.1 角频率的影响
角频率在范围\omega \in \left( {0,6.5} \right]内取值,其余参数为a = 0.1658,以及在表2中的序号3~14的参数。使用直接数值积分法构建分岔图。选择适当的初始值,利用MATLAB中的求解器ode45对式(29)进行足够多的周期的数值积分,以衰减瞬态的影响。使用odeset函数的“Events”选项指定事件函数,在活塞碰撞钻具(即{y_1} - {y_3} - g = 0)时停止积分,然后使用新的初始条件式(32)重新开始积分。然后进行额外的100个周期,并在时间点\tau = iT = i\dfrac{{2\pi }}{\omega }处绘制解,其中i = 1,2,\cdots,100。随后,缓慢增大或减小\omega 并重复该过程,将前一步骤的最终解作为新的初始值。
图4为关于角频率的单参数分岔图。(△)点为粘滞生成的解,(○)点为非粘滞正向扫描生成的解,(●)点为非粘滞反向扫描生成的解。粘滞解的范围是\omega \in \left( {0,2.08} \right],非粘滞解的范围是\omega \in \left( 2.08, 6.5 \right]。大部分正向和反向生成的解是重合的。从放大视图可以看出,当\omega \in [4.94,5.36]时,正向和反向生成的解没有重合,也就是说在此范围内,模型具有多稳态性,取两个典型的角频率{\omega _1} = 5.22和{\omega _2} = 5.33,用于后续分析多稳态的吸引域。
吸引域是指系统在一定条件下能够吸引相邻轨道向其靠近的区域[20, 21]。图5为不同角频率的吸引域。计算参数为a = 0.1658,以及在表2中的序号3~14的参数。如图4中多稳态的放大视图所示,分别取{\omega _1} = 5.22和{\omega _2} = 5.33。图5(a)为{\omega _1} = 5.22时的吸引域。图5(b)为{\omega _1} = 5.22时周期-1和周期-3的相位。将活塞的初始位移和初始速度{y_1} \times {y_2} = [ - 0.06,0.06] \times [ - 0.2,0.2]作为初始值平面,其余初始值{y_3}\sim {y_{10}}均为0。取平面内均匀随机分布的10000个点,通过式(29)和式(32)的计算,经过足够长的时间来消除瞬态的影响,得到两种不同的稳定周期解,如图5(b)所示,分别标记为“周期-1”和“周期-3”。如图5(a)所示,初始值平面被划分为两个区域。周期-1和周期-3吸引域的面积分别占总面积的9.7%和90.3%。图5(c)为{\omega _2} = 5.33时的吸引域。图5(d)为{\omega _2} = 5.33时周期-3和准周期的相位。取{y_1} \times {y_2}平面内均匀随机分布的10000个点,通过式(29)和式(32)的计算,经过足够长的时间来消除瞬态的影响,得到两种不同的稳定周期解,如图5(d)所示,分别标记为“周期-3”和“准周期”。如图5(c)所示,初始值平面被划分为两个区域。周期-3和准周期吸引域的面积分别占总面积的95.3%和4.7%。
从图5可以看出,当\omega \in [4.94,5.36]时,角频率具有多稳态,周期-3占据初始值平面的面积更大,所以周期-3比周期-1和准周期具有更大的全局稳定性。
3.2 振幅的影响
另一个影响系统行为的参数是液压作用力的振幅a。振幅在范围a \in (0,1]内取值,其余参数为\omega = 5.22,以及在表2中的序号3~14的参数。图6为关于振幅的单参数分岔图。可以看到相同数值a对应多个稳定周期解。为了研究多稳态性,取数值{a_1} = 0.5,{a_2} = 0.66和{a_3} = 0.8。
图7为不同振幅的吸引域。计算参数为\omega = 5.22,以及在表2中的序号3~序号14的参数。如图6所示,取{a_1} = 0.5、{a_2} = 0.66和{a_3} = 0.8。图7(a)为{a_1} = 0.5时的吸引域。图7(b)为{a_1} = 0.5时周期-1-1和周期-3的相位。将活塞的初始位移和初始速度{y_1} \times {y_2} = [ - 0.15,0.15] \times [ - 0.6,0.6]作为初始值平面,其余初始值{y_3}\sim {y_{10}}均为0。取平面内均匀随机分布的10000个点,通过式(29)和式(32)的计算,经过足够长的时间来消除瞬态的影响,得到两种不同的稳定周期解,如图7(b)所示,分别标记为“周期-1-1”和“周期-3”。如图7(a)所示,初始值平面被划分为两个区域。周期-1-1和周期-3吸引域的面积分别占总面积的13.9%和86.1%。图7(c)为{a_2} = 0.66时的吸引域。图7(d)为{a_2} = 0.66时周期-1-1、周期-1-2和周期-3的相位。取{y_1} \times {y_2}平面内均匀随机分布的10000个点,通过式(29)和式(32)的计算,经过足够长的时间来消除瞬态的影响,得到三种不同的稳定周期解,如图7(d)所示,分别标记为“周期-1-1”、“周期-1-2”和“周期-3”。如图7(c)所示,初始值平面被划分为三个区域。周期-1-1、周期-1-2和周期-3吸引域的面积分别占总面积的33.4%、5.4%和61.2%。图7(e)为{a_3} = 0.8时的吸引域。图7(f)为{a_3} = 0.8时周期-1-2和周期-3的相位。取{y_1} \times {y_2}平面内均匀随机分布的10000个点,通过式(29)和式(32)的计算,经过足够长的时间来消除瞬态的影响,得到两种不同的稳定周期解,如图7(f)所示,分别标记为“周期-1-2”和“周期-3”。如图7(e)所示,初始值平面被划分为两个区域。周期-1-2和周期-3吸引域的面积分别占总面积的64%和36%。
从图7可以看出,当0 < a {\leqslant} 1时,随着a的增加,周期-1-2解的吸引域稳定性逐渐增加,而周期-3解的吸引域稳定性逐渐减小。周期-1-1解的吸引域稳定性起初增加,然后变为0。
4 实验和模型的比较
图8为液压凿岩机钻凿岩石的过程和数据采集系统。如图8(a)所示,从左往右依次是活塞、钎尾、连接套、钻杆、钻头和岩石。实验现场如图8(b)所示,激光位移传感器安装在活塞后面,测量其位移和速度。测量结果进入电脑和数据采集系统。采样频率为200 kHz,并使用了适当的抗混叠滤波器处理信号。每个运行条件下的测量时间为2 s~4 s,仅在凿岩机运行平稳状态下收集数据。如图8(c)所示,测量结果进入电脑和数据采集系统。
图9为实验和模型的活塞运动的比较,展示了2个工作循环,包括4个行程和5个特殊点。图9(a)为模型和实验的活塞位移曲线,描述了两个完整的工作循环。循环1的时间为0 s~0.0262 s,循环2的时间为0.0262~0.052 s。点SP1为模型的活塞起始位置,此时位移和速度均为0。图9(b)为模型和实验的活塞速度曲线,展示了四个行程:模型冲程加速(Model Impact Acceleration, MIA),模型冲程减速(Model Impact Deceleration, MID),模型回程加速(Model Return Acceleration, MRA)和模型回程减速(Model Return Deceleration, MRD)。MIA行程在t = 0.004 s时开始,此时活塞静止,活塞加速朝着钻具运动。在t = 0.0109 s时,活塞达到最大速度(点SP2)6.29 m/s。之后,活塞开始减速,进入MID行程。在MID行程中,液压作用力与运动方向相反。在t = 0.016 s时,活塞与钻具碰撞。碰撞前活塞速度(点SP3)为2.85 m/s,碰撞后活塞速度(点SP4)为-0.45 m/s。活塞与钻具之间的碰撞关系由式(5)描述。在t = 0.0233 s时,活塞在MRA行程中达到绝对值最大速度(点SP5)-5.8 m/s。随后,活塞减速进入MRD行程,MRD行程在t = 0.03 s时结束。
从图9可以看出,模型和实验的动态特性基本保持一致,实验初步验证了模型的正确性。用实验验证周期-1、周期-3和准周期的状态,是在此研究基础上将来重点研究的方向。
5 结论
本文将凿岩机钻进岩石的过程,建立成四自由度力学模型;利用直接数值积分方法,将液压作用力的角频率和振幅作为控制参数,进行分岔和吸引域分析;介绍液压凿岩机钻凿岩石的数据采集系统,比较模型和实验测得的活塞的位移和速度。得到主要结论如下:
(1)将液压力的角频率作为控制参数时,发现当ω < 2.08时,模型为粘滞模式。当2.08≤ω≤6.5时,模型为非粘滞模式。当4.94≤ω≤5.36时,模型为多稳态。将液压力的振幅作为控制参数时,发现当0<a≤1时,模型为多稳态。
(2)分析角频率的吸引域,发现周期-3解表现出了比周期-1解和准周期解更大的全局稳定性。分析振幅的吸引域,发现随着振幅的增加,周期-1-2解的吸引域逐渐增加,而周期-3解的吸引域逐渐减小。周期-1-1解的吸引域起初增加,然后变为0。
(3)比较了模型和实验测得的活塞的位移和速度。实验和模型数据基本吻合,验证了模型和实验动态特性的一致性。
-
表 1 凿岩机-岩石模型的量纲参数
Table 1 Dimension parameters in the drifter-rock model
含义 数值 液压力的角频率\varOmega /(rad/s) 240 液压油阻尼c1/(N/(m/s)) 1×103 结构阻尼c2/(N/(m/s)) 2.45×103 液压力的振幅Fa/N 1.1417×104 液压力的垂直偏移量Fb/N 6.7154×104 活塞和钻具的距离G/m 0.025 液压油刚度k1/(N/m) 2×105 结构刚度k2/(N/m) 6.885×104 加载刚度kl/(N/m) 11.46×106 卸载刚度ku/(N/m) 13.06×106 活塞质量m1/kg 11.5 钻具质量m2/kg 32.6 钻头质量m3/kg 1 表 2 凿岩机-岩石模型的无量纲参数
Table 2 Dimensionless parameters in the drifter-rock model
序号 符号 含义 数值 1 \omega 液压作用力的角频率 (0,6.5] 2 a 液压作用力的振幅 (0,1] 3 b 液压作用力的垂直偏移量 0.0975 4 \alpha 活塞和钻具的质量比 0.35 5 \beta 刚度比 2.9 6 \gamma 阻尼系数比 0.4 7 \lambda 卸载刚度和结构刚度的比值 189.7 8 \kappa 加载刚度和结构刚度的比值 166.4 9 \mu 钻具和钻头的质量比 32.6 10 \zeta 阻尼比 0.8 11 g 活塞和钻具的间隙 0.025 12 {n_l} 加载指数 2 13 {n_u} 卸载指数 2 14 e 碰撞恢复系数 0.3 -
[1] 耿晓光. 双缓冲液压凿岩机冲击特性及影响因素研究[D]. 北京: 北京科技大学, 2019. GENG Xiaoguang. Research on the influence factors of percussion characteristics for hydraulic rock drill with double damper system [D]. Beijing: University of Science and Technology Beijing, 2019. (in Chinese)
[2] MA W, GENG X G, JIA C Z, et al. Percussion characteristic analysis for hydraulic rock drill with no constant-pressurized chamber through numerical simulation and experiment [J]. Advances in Mechanical Engineering, 2019, 11(4): 1 − 11. (查阅网上资料, 不确定页码信息是否正确, 请确认) .
[3] SEO J, NOH D K, LEE G H, et al. A percussion performance analysis for rock-drill drifter through simulation modeling and experimental validation [J]. International Journal of Precision Engineering and Manufacturing, 2016, 17(2): 163 − 170. doi: 10.1007/s12541-016-0021-0
[4] YONG G, DESHUN L, SHUYI Y, et al. Hydraulic–mechanical coupling modeling by bond graph for impact system of a high frequency rock drill drifter with sleeve distributor [J]. Automation in Construction, 2016, 63: 88 − 99. doi: 10.1016/j.autcon.2015.12.012
[5] YANG S Y, OU Y B, GUO Y, et al. Analysis and optimization of the working parameters of the impact mechanism of hydraulic rock drill based on a numerical simulation [J]. International Journal of Precision Engineering and Manufacturing, 2017, 18(7): 971 − 977. doi: 10.1007/s12541-017-0114-4
[6] GENG X G, MA W, MA F, et al. Cavitation erosion of the damping piston in double damping system of hydraulic rock drill and its influencing factors [J]. AIP Advances, 2018, 8(11): 115025. doi: 10.1063/1.5058162
[7] GENG X, MA W, MA F, et al. Design on key parameters of double damper system for hydraulic rock drill [J]. International Applied Mechanics, 2019, 55(6): 700 − 705. doi: 10.1007/s10778-019-00989-5
[8] YANG Z Y, JI H, LI Y G. Analysis on the main design parameters influencing the impact efficiency of dual-chamber-controlled hydraulic drifter [J]. International Journal of Precision Engineering and Manufacturing, 2018, 19(12): 1781 − 1791. doi: 10.1007/s12541-018-0207-8
[9] LI Y L, BIN Z, TU Y. Research on impact performance of hydraulic rock drill with floating characteristics of double damping system [J]. Shock and Vibration, 2019, 2019(1): 2942890. doi: 10.1155/2019/2942890
[10] SHAH K H. Electrically controlled hydraulic rock drill [D]. Tampere: Tampere University, 2019.
[11] HASHIBA K, FUKUI K, LIANG Y Z, et al. Force-penetration curves of a button bit generated during impact penetration into rock [J]. International Journal of Impact Engineering, 2015, 85: 45 − 56. doi: 10.1016/j.ijimpeng.2015.06.007
[12] LUNDBERG B, COLLET P. Optimal wave shape with respect to efficiency in percussive drilling with detachable drill bit [J]. International Journal of Impact Engineering, 2015, 86: 179 − 187. doi: 10.1016/j.ijimpeng.2015.06.021
[13] LIAO M L, WIERCIGROCH M, SAYAH M, et al. Experimental verification of the percussive drilling model [J]. Mechanical Systems and Signal Processing, 2021, 146: 107067. doi: 10.1016/j.ymssp.2020.107067
[14] AJIBOSE O K, WIERCIGROCH M, PAVLOVSKAIA E, et al. Drifting impact oscillator with a new model of the progression phase [J]. Journal of Applied Mechanics, 2012, 79(6): 061007. doi: 10.1115/1.4006379
[15] LIAO M L. Dynamic methods of stiffness identification in impacting systems for rotary-percussive drilling applications [D]. Aberdeen: University of Aberdeen, 2016.
[16] 丁兰, 丁彪, 吴巧云, 等. 含双自由度周期振子的平行并联梁弯曲振动带隙特性[J]. 工程力学, 2023, 40(10): 1 − 10,57. doi: 10.6052/j.issn.1000-4750.2022.01.0086 DING Lan, DING Biao, WU Qiaoyun, et al. Flexural vibration bandgap characteristics of double beams in parallel with oscillators with two degrees of freedom [J]. Engineering Mechanics, 2023, 40(10): 1 − 10,57. (in Chinese) doi: 10.6052/j.issn.1000-4750.2022.01.0086
[17] 李鑫, 唐贞云, 杜修力. 半无限介质多自由度频响函数离散时间有理近似的稳定参数识别[J]. 工程力学, 2024, 41(8): 1 − 10. doi: 10.6052/j.issn.1000-4750.2022.05.0446 LI Xin, TANG Zhenyun, DU Xiuli. Identification of stable parameters for discrete-time rational approximation of mdof frequency response functions in semi-infinite media [J]. Engineering Mechanics, 2024, 41(8): 1 − 10. (in Chinese) doi: 10.6052/j.issn.1000-4750.2022.05.0446
[18] 武双双, 金映丽, 闫明, 等. 反弹对撞式冲击放大器机理研究[J]. 工程力学, 2023, 40(8): 235 − 242. doi: 10.6052/j.issn.1000-4750.2021.12.0957 WU Shuangshuang, JIN Yingli, YAN Ming, et al. Research on the mechanism of rebound collision shock amplifier [J]. Engineering Mechanics, 2023, 40(8): 235 − 242. (in Chinese) doi: 10.6052/j.issn.1000-4750.2021.12.0957
[19] MA W, MAPURANGA T. Bifurcation and path-following continuation analysis of periodic orbits of an extended fermi oscillator model [J]. Nonlinear Dynamics, 2023, 111(9): 7993 − 8020. doi: 10.1007/s11071-022-08080-4
[20] 吕小红, 罗冠炜. 小型振动冲击式打桩机的非线性动力学分析[J]. 工程力学, 2013, 30(11): 227 − 232. doi: 10.6052/j.issn.1000-4750.2012.07.0523 LÜ Xiaohong, LUO Guanwei. Analysis of nonlinear dynamics of a small vibro-impact driver [J]. Engineering Mechanics, 2013, 30(11): 227 − 232. (in Chinese) doi: 10.6052/j.issn.1000-4750.2012.07.0523
[21] STENDER M, HOFFMANN N. BSTAB: An open-source software for computing the basin stability of multi-stable dynamical systems [J]. Nonlinear Dynamics, 2022, 107(2): 1451 − 1468. doi: 10.1007/s11071-021-06786-5