FAST索牵引并联机器人的动力学建模与仿真
赵雅聪1,2, 王启明1    
1. 中国科学院 国家天文台, 北京 100101;
2. 中国科学院大学, 北京 100049
摘要:500 m口径球面射电望远镜(five-hundred-meter aperture spherical radio telescope,FAST)中馈源支撑系统是一个巨型的索牵引并联机器人,钢索牵引动平台在数百米空间内运动,包含变长度柔性索的运动和动平台的刚体运动。该文以索牵引并联机器人为研究对象,首先,采用绝对节点坐标法建立变长度钢索单元模型,并以此为基础利用虚功原理建立变长度索单元的动力学方程组;其次,对刚性动平台进行动力学分析,利用钢索与动平台间的约束得到六索牵引并联机器人整体的动力学方程组,其中包含柔性索和动平台的自由度;最后,使用MATLAB对建立的动力学模型进行了数值计算及仿真,将所得的仿真数据与对应的FAST实际运行数据进行对比,验证了该模型的可行性。
关键词FAST    索牵引并联机器人    变长度柔索    动力学模型    动力学仿真    
Dynamic modeling of the FAST cable-driven parallel robots
ZHAO Yacong1,2, WANG Qiming1    
1. National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China;
2. University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: The feed support system of the five-hundred-meter aperture spherical radio telescope (FAST) is a kind of giant cable-driven parallel robot that uses flexible cables to move the platform in the working space. Thus, the system combines variable length flexible cables and rigid body motion. This study analyzed the cable-driven parallel robot. The variable length cable elements were analyzed based on the absolute node coordinate method, the dynamic equations of the cable elements were then derived using the virtual work principle. Then, the dynamics of the rigid moving platform were analyzed, the dynamics equations of the cable-driven parallel robot with constraints between the cables and the moving platform were derived. Thus, the model included the degrees of freedom of the flexible cables and the platform. MATLAB was used for the simulations. The results compare well with measured data to verify the model feasibility.
Key words: FAST    cable-driven parallel robots    variable length flexible cable    dynamics model    dynamics simulations    

500 m口径球面射电望远镜(five-hundred-meter aperture spherical radio telescope,FAST)[1]是目前世界上最大、最灵敏的单口径球面射电望远镜。FAST工程于2011年正式开工建设,2016年9月建设完成后进入调试及试观测阶段,并于2020年1月通过国家验收,正式投入运行。FAST与目前国际上的单口径射电望远镜相比,拥有3项自主创新:利用天然的喀斯特地貌作为台址、可实现实时变形的主动反射面系统以及光机电一体化轻型馈源索支撑系统。全新的设计思路使望远镜突破了在地面建造全可动天线的百米工程极限,开创了建造巨型射电望远镜的新模式。

FAST的3项自主创新之一——馈源支撑系统是由多级机构组成的刚柔结合的并联机器人系统,该系统将馈源精准定位于瞬时焦点,实现对天体的高精度指向跟踪观测。FAST馈源支撑系统包含馈源舱(结构框架、AB轴、Stewart平台和馈源接收机)、6根并联钢索及与其对应的支撑塔和卷扬机。钢索一端固定于馈源舱的铰接点上,另一端通过支撑塔顶端和底部的定滑轮转向与卷扬机相连,卷扬机带动6根钢索拖动重达30 t的馈源舱在开口直径约为207 m的球冠面上运动,使馈源接收机对空间中汇集的射电信号进行跟踪,馈源舱的工作区域如图 1所示,O为望远镜反射面的球心,R为反射面半径。

图 1 馈源舱的工作区域示意图

应用于FAST中的六索牵引并联机器人具有大负载和大工作空间等特点,对钢索理论模型的研究决定了整个系统的力学特性。许多学者对FAST馈源支撑系统进行了研究分析,黄亮等[2]基于静力学分析对铰接点和塔顶滑轮间钢索的悬链线模型进行简化,提出了赝曲线模型;Yin等[3]通过重新构造悬链线中的参数关系,提出了一套可以近似解析求解的高速算法;李辉等[4]使用张力均衡分配的优化原则对FAST馈源支撑系统中柔性钢索牵引并联机构进行了静力学分析,对馈源支撑系统进行了1∶1的原型虚拟样机的全程仿真,在索牵引并联机构的动力学分析中使用其静力平衡位置进行了线性简化,选择数个位置点及其附近轨迹进行仿真模拟。以上研究没有充分考虑柔性索的动力学特性对动平台定位精度的影响。Du等[5]使用位移坐标,采用中心差分法对索牵引并联机器人中的变长度柔性钢索进行了动力学分析,通过与静态运动位移相比较,证明了对柔性钢索进行动力学建模的必要性。当馈源舱在三维空间运动时,外界的干扰、舱体位置的误差、索长的变化、速度等都可能引起钢索的振动,而钢索的振动又会导致馈源舱的振动,因此建立索牵引并联机器人的精确动力学模型的前提就是建立精确的钢索动力学模型。

对于空间中运动的具有大变形、大位移的梁单元,Shabana[6]提出了基于有限元和连续介质力学的绝对节点坐标法,使用节点的位移和梯度矢量作为广义坐标来描述空间梁单元的位移和变形;吴懋琦等[7]采用绝对节点坐标法对柔性结构大变形下非线性的平面梁应变-位移关系进行精确描述, 构造了一种逆梯度缩减ANCF平面梁单元,并采用数值优化方法进行位移的求解;Luo等[8]利用绝对节点坐标法对大位移和大变形的变长度系绳进行了描述,采用自然坐标法描述3颗卫星的整体运动,研究了绳系卫星编队的动力学问题;王忠民等[9]基于绝对节点坐标法,建立了一种变长度的Euler-Bernoulli梁单元模型,分析了材料参数和伸展规律对悬臂梁系统的末端非线性挠度响应的影响;Bulín等[10]利用绝对节点坐标法对滑轮绳索系统进行动力学分析,并将模拟数据与实验数据进行对比,验证了该方法的可行性。由于柔性系统的动力学方程刚性问题突出,目前多采用刚性求解器进行结算,如隐式Runge-Kutta法、向后差分公式方法等。Ferravante等[11]通过引入Rayleigh阻尼矩阵,即将质量矩阵和刚度矩阵按比例相加使数值稳定性和计算效率有所提高;齐朝晖等[12]对柔性体的应力进行均匀化,用均匀化后的应力计算柔性体的变形虚功率,即在建模过程中滤除高频振荡分量,提高了求解动力学方程组的计算效率。

本文基于FAST馈源支撑系统中所使用的索牵引并联机器人,提出一种该类型机器人的动力学建模方法:应用绝对节点坐标法描述具有时变性的钢索单元,根据虚功原理对钢索的动态特性进行分析,结合钢索和馈源舱间的约束以及钢索出索点的约束,得到六索牵引并联机构整体的动力学方程组。该动力学模型不仅考虑了馈源舱的动态特性,而且包括了大跨度柔性钢索的振动。使用MATLAB进行计算仿真,将仿真数据与实际观测数据进行对比,验证了所建立模型的有效性。

1 索牵引并联机器人的动力学方程

本文对FAST馈源支撑系统中使用的六索牵引并联机器人进行分析,其中6根钢索通过2个定滑轮由卷扬机带动,长度随时间变化;馈源舱为末端执行机构,可将其视为刚性动平台。分别建立变长度钢索和动平台的运动方程,根据钢索与动平台的约束,建立六索牵引并联机器人系统的动力学方程。以图 1中反射面的球心O点为原点,建立全局坐标系,以此为基础对索牵引并联机器人进行分析。

1.1 变长度钢索单元

假设钢索在运动过程中只受拉力作用,并且不考虑其剪切和扭转变形。基于绝对节点坐标法,在全局坐标系中对钢索进行动力学分析。先对单根钢索进行分析,随着收索或放索的操作,单根钢索任意时刻的长度为

$L(t)=L_0+\Delta L(t) .$ (1)

其中:L0为钢索初始时刻的长度,ΔL(t)为钢索随时间t变化的长度。

将单根钢索进行离散,如图 2所示,钢索被划分为n个单元,有n+1个节点。如果钢索单元长度le发生变化,则第i个钢索单元在任意时刻的长度为

$l_{\mathrm{e}}^i(t)=l_0^i+\Delta l_{\mathrm{e}}^i(t), i=1, 2, \cdots, n .$ (2)
图 2 单根钢索的划分

其中:l0i为钢索单元初始时刻的长度,Δlei(t)为钢索单元随时间t变化的长度。设置钢索单元的长度阈值lminlmax,当钢索单元长度大于lmax时,增加1个单元节点;钢索单元长度小于lmin时,减少1个单元节点。

将钢索单元视为一维2节点单元,如图 3所示。在三维坐标系中每个节点的坐标由该点的位置和切线斜率的6个广义坐标组成,每个单元包含2个节点,即有12个广义坐标。第i个钢索单元的广义坐标矢量可表示为

$\boldsymbol{u}^i=\left[\boldsymbol{e}_{i-1}, \boldsymbol{e}_i\right]^{\mathrm{T}}, i=1, 2, \cdots, n.$ (3)
图 3 钢索单元

其中ei为第i个节点的广义坐标矢量,可表示为

$\boldsymbol{e}_i=\left[r_{i X}, r_{i Y}, r_{i Z}, \frac{\partial r_{i X}}{\partial \chi}, \frac{\partial r_{i Y}}{\partial \chi}, \frac{\partial r_{i Z}}{\partial \chi}\right].$ (4)

其中:riXriYriZ为节点i的位置矢量ri在全局坐标系中的分量;$ \frac{\partial r_{i x}}{\partial \chi}、\frac{\partial r_{i Y}}{\partial \chi}、\frac{\partial r_{i Z}}{\partial \chi}$为节点i处切线斜率,χ为单元内弧长坐标。ri为该钢索单元内任意一点p的全局坐标矢量,可表示为

$\boldsymbol{r}^i=\boldsymbol{S} \boldsymbol{u}^i, $ (5)
$\boldsymbol{S}=\left[\begin{array}{l} \left(1-3 \xi^2+2 \xi^3\right) \boldsymbol{I} \\ l_{\mathrm{e}}^i(t)\left(\xi-2 \xi^2+\xi^3\right) \boldsymbol{I} \\ \left(3 \xi^2-2 \xi^3\right) \boldsymbol{I} \\ l_{\mathrm{e}}^i(t)\left(\xi^3-\xi^2\right) \boldsymbol{I} \end{array}\right]^{\mathrm{T}} .$ (6)

其中:S为钢索单元的形函数矩阵,I为3×3的单元矩阵,ξ=χ/lei(t)。该钢索单元上任意点在全局坐标系中的表达如式(5)所示,将其对时间t求一阶导数和二阶导数,得

$\dot{\boldsymbol{r}}^i=\dot{\boldsymbol{S}} \boldsymbol{u}^i+\boldsymbol{S} \boldsymbol{u}^i.$ (7)
$\ddot{\boldsymbol{r}}^i=\ddot{\boldsymbol{S}} \boldsymbol{u}^i+2 \dot{\boldsymbol{S}} \dot{\boldsymbol{u}}^i+\boldsymbol{S} \ddot{\boldsymbol{u}}^i.$ (8)

如果钢索单元的长度不变,即lei(t)=l0i时,式(7)和(8)中形函数矩阵S对时间t的一阶导数和二阶导数都为0;钢索单元的长度随时间发生变化时,S对时间t的一阶导数和二阶导数可表示为:

$\dot{\boldsymbol{S}}=\frac{\partial \boldsymbol{S}}{\partial \xi} \dot{\xi}+\frac{\partial \boldsymbol{S}}{\partial l_{\mathrm{e}}^i} \dot{l}_{\mathrm{e}}^i, $ (9)
$\ddot{\boldsymbol{S}}=\frac{\partial^2 \boldsymbol{S}}{\partial \dot{\xi}^2} \dot{\xi}^2+\frac{\partial \boldsymbol{S}}{\partial \xi} \ddot{\xi}+2 \frac{\partial^2 \boldsymbol{S}}{\partial t \partial \xi} \dot{\xi}+\frac{\partial^2 \boldsymbol{S}}{\partial l_{\mathrm{e}}^i} \ddot{l}_{\mathrm{e}}^i, $ (10)
$\dot{\xi}=\frac{\dot{l}_{\mathrm{e}}^i}{l_{\mathrm{e}}^i}(1-\xi), $ (11)
$\ddot{\xi}=\frac{\ddot{l}_{\mathrm{e}}^i}{l_{\mathrm{e}}^i}(1-\xi)-2 \frac{\left(\dot{l}_{\mathrm{e}}^i\right)^2}{\left(l_{\mathrm{e}}^i\right)^2}(1-\xi) .$ (12)

已知钢索单元中任意一点的位置、速度和加速度,在虚位移δr钢索单元惯性力的虚功可表示为

$\begin{gathered} \boldsymbol{\delta} W_{\mathrm{eI}}^i=-\int_0^{l_{\mathrm{e}}} \rho A \boldsymbol{\delta}\left(\boldsymbol{r}^i\right)^{\mathrm{T}} \ddot{\boldsymbol{r}}^i \mathrm{~d} \chi= \\ -\int_0^{l_{\mathrm{e}}} \rho A \boldsymbol{\delta}\left(\boldsymbol{r}^i\right)^{\mathrm{T}} \ddot{\boldsymbol{r}}^i \frac{\mathrm{d} \chi}{\mathrm{d} \xi} \mathrm{d} \xi= \\ -\boldsymbol{\delta}\left(\boldsymbol{u}^i\right)^{\mathrm{T}} l_{\mathrm{e}}^i \int_0^1 \rho A \boldsymbol{S}^{\mathrm{T}} \ddot{\boldsymbol{r}}^i \mathrm{~d} \xi= \\ -\boldsymbol{\delta}\left(\boldsymbol{u}^i\right)^{\mathrm{T}}\left(\boldsymbol{M}_{\mathrm{eI}}^i \ddot{\boldsymbol{u}}^i+\boldsymbol{C}_{\mathrm{eI}}^i \dot{\boldsymbol{u}}^i+\boldsymbol{K}_{\mathrm{eI}}^i \boldsymbol{u}^i\right). \end{gathered}$ (13)

其中:A为钢索的横截面积, ρ为钢索的密度。MeIi为广义质量矩阵, CeIi为广义阻尼矩阵, KeIi为广义刚度矩阵,可表示为:

$\boldsymbol{M}_{\mathrm{eI}}^i=l_{\mathrm{e}}^i \int_0^1 \rho A \boldsymbol{S}^{\mathrm{T}} \boldsymbol{S} \mathrm{d} \xi, $ (14)
$\boldsymbol{C}_{\mathrm{eI}}^i=2 l_{\mathrm{e}}^i \int_0^1 \rho A \boldsymbol{S}^{\mathrm{T}} \dot{\boldsymbol{S}} \mathrm{d} \xi,$ (15)
$\boldsymbol{K}_{\mathrm{eI}}^i=l_{\mathrm{e}}^i \int_0^1 \rho A \boldsymbol{S}^{\mathrm{T}} \ddot{\boldsymbol{S}} \mathrm{d} \xi.$ (16)

钢索轴向形变引起的钢索单元弹性力的虚功为

$\begin{gathered} \boldsymbol{\delta} W_{\mathrm{eE}}^i=\int_0^{l_{\mathrm{e}}} E A \delta\left(\varepsilon^i\right)^{\mathrm{T}} \varepsilon^i \mathrm{~d} \chi= \\ \boldsymbol{\delta}\left(\boldsymbol{u}^i\right)^{\mathrm{T}} l_{\mathrm{e}}^i \int_0^1 \frac{\partial \varepsilon^i}{\partial \boldsymbol{u}^i} E A \varepsilon^i \mathrm{~d} \xi= \\ \boldsymbol{\delta}\left(\boldsymbol{u}^i\right)^{\mathrm{T}} l_{\mathrm{e}}^i \int_0^1 E A \varepsilon^i\left(\boldsymbol{S}_\chi^{\mathrm{T}} \boldsymbol{S}_\chi\right) \boldsymbol{u}^i \mathrm{~d} \xi= \\ \boldsymbol{\delta}\left(\boldsymbol{u}^i\right)^{\mathrm{T}} \boldsymbol{Q}_{\mathrm{eE}}^i . \end{gathered}$ (17)

其中:E为钢索的弹性模量,Sχ=S/∂χQeEi=$ l_{\mathrm{e}}^i \int_0^1 E A \varepsilon^i\left(\boldsymbol{S}_\chi^{\mathrm{T}} \boldsymbol{S}_\chi\right) \boldsymbol{u}^i \mathrm{~d} \xi$为钢索单元受到的广义弹性力,εi为钢索单元的轴向拉伸应变,可表示为

$\varepsilon^i=\sqrt{\left(\boldsymbol{r}_\chi^i\right)^{\mathrm{T}} \boldsymbol{r}_\chi^i}-1 .$ (18)

其中rχi=ri/∂χ。对于钢索单元弹性力的计算,Berzeri等[13]建立了纵向弹性力模型L1、L2和L3,其中,L2、L3适用于较大的纵向变形,L1对纵向变形仅计算了相邻节点的直线距离,在单元已具有弯曲变形或单元尺寸较大时可能会发生离散。Wago等[14]改进了弹性力的计算模型,能够更精确地计算钢索单元纵向拉伸应变,同时使计算更容易收敛。钢索单元的轴向拉伸应变也可以表示为

$\varepsilon^i=\frac{l_{\mathrm{es}}^i-l_{\mathrm{e}}^i}{l_{\mathrm{e}}^i},$ (19)
$\begin{gathered} l_{\mathrm{es}}^i=\int_0^{l_{\mathrm{e}}} \sqrt{\left(\boldsymbol{u}^i\right)^{\mathrm{T}} \boldsymbol{S}_\chi^{\mathrm{T}} \boldsymbol{S}_\chi \boldsymbol{u}^i} \mathrm{~d} \chi= \\ \int_0^{l_{\mathrm{e}}} \sqrt{\left(\boldsymbol{u}^i\right)^{\mathrm{T}} \boldsymbol{S}_\chi^{\mathrm{T}} \boldsymbol{S}_\chi \boldsymbol{u}^i-1+1} \mathrm{~d} \chi \approx \\ \int_0^{l_{\mathrm{e}}} 1+\frac{\left(\boldsymbol{u}^i\right)^{\mathrm{T}} \boldsymbol{S}_\chi^{\mathrm{T}} \boldsymbol{S}_\chi \boldsymbol{u}^i-1}{2} \mathrm{~d} \chi= \\ \frac{1}{2} \int_0^{l_{\mathrm{e}}}\left(1+\left(\boldsymbol{u}^i\right)^{\mathrm{T}} \boldsymbol{S}_\chi^{\mathrm{T}} \boldsymbol{S}_\chi \boldsymbol{u}^i\right) \mathrm{d} \chi= \\ \frac{1}{2}\left[\left(\boldsymbol{u}^i\right)^{\mathrm{T}} l_{\mathrm{e}}^i \int_0^1 \boldsymbol{S}_\chi^{\mathrm{T}} \boldsymbol{S}_\chi \mathrm{d} \xi \hat{\xi} \boldsymbol{u}^i+l_{\mathrm{e}}^i\right] . \end{gathered}$ (20)

其中lesi为钢索单元变形后的纵向长度,将式(20)代入式(17)和(19)即可求得钢索单元弹性力的虚功。

在忽略外界干扰的情况下,钢索所受外力为重力,重力虚功δWeFi可表示为

$\begin{gathered} \boldsymbol{\delta} W_{\mathrm{eF}}^i=\int_0^{l_{\mathrm{e}}} \boldsymbol{\delta}\left(\boldsymbol{r}^i\right)^{\mathrm{T}} \boldsymbol{f} \mathrm{d} \chi= \\ \boldsymbol{\delta}\left(\boldsymbol{u}^i\right)^{\mathrm{T}} l_{\mathrm{e}}^i \int_0^1 \boldsymbol{S}^{\mathrm{T}} \boldsymbol{f} \mathrm{d} \xi=\boldsymbol{\delta}\left(\boldsymbol{u}^i\right)^{\mathrm{T}} \boldsymbol{Q}_{\mathrm{eF}}^i. \end{gathered}$ (21)

其中:f=[0,0,-ρg]T; $ \boldsymbol{Q}_{\mathrm{eF}}^i=l_{\mathrm{e}}^i \int_0^1 \boldsymbol{S}^{\mathrm{T}} \boldsymbol{f} \mathrm{d} \xi$为广义重力。

根据虚功原理可知,钢索单元的惯性力、弹性力、广义外力虚功和为0,即

$\boldsymbol{\delta} W_{\mathrm{eI}}^i-\boldsymbol{\delta} W_{\mathrm{eE}}^i+\boldsymbol{\delta} W_{\mathrm{eF}}^i={\bf{0}} .$ (22)

钢索单元的动力学方程可表示为

$\boldsymbol{M}_{\mathrm{eI}}^i \ddot{\boldsymbol{u}}^i+\boldsymbol{C}_{\mathrm{eI}}^i \dot{\boldsymbol{u}}^i+\boldsymbol{K}_{\mathrm{eI}}^i \boldsymbol{u}^i+\boldsymbol{Q}_{\mathrm{eE}}^i=\boldsymbol{Q}_{\mathrm{eF}}^i \text {. }$ (23)

当钢索单元长度不变时,式(23)中质量矩阵MeIi为常数矩阵,CeIiKeIi均消除,即为传统的绝对节点坐标法钢索单元的动力学方程组。

1.2 动平台

以馈源舱质心为原点建立与馈源舱全固连的局部坐标系,对馈源舱进行分析。将馈源舱视为刚性动平台,六索牵引并联机器人是6根钢索和重力组成完全约束并联机构,其动平台的姿态xe可用6个独立坐标来描述

$\boldsymbol{x}_{\mathrm{e}}=[x, y, z, \gamma, \beta, \alpha]^{\mathrm{T}}=\left[\boldsymbol{x}_{\mathrm{e} 1}^{\mathrm{T}}, \boldsymbol{x}_{\mathrm{e} 2}^{\mathrm{T}}\right]^{\mathrm{T}} .$ (24)

其中:xe1为动平台质心在全局坐标系中的位置矢量, xe2为动平台在局部坐标系转动的Euler角。动平台先绕z轴转动γ角,即偏航角;再绕新的y轴转动β角,即俯仰角;最后再绕新的x轴旋转α角,即翻滚角。根据3个Euler角可以得到旋转矩阵R(xe2)。

在忽略外界干扰,刚性动平台在不受外力的情况下动力学方程可写为

$\boldsymbol{M}_{\mathrm{p}} \ddot{\boldsymbol{x}}_{\mathrm{e}}+\boldsymbol{C}_{\mathrm{p}} \dot{\boldsymbol{x}}_{\mathrm{e}}=\boldsymbol{G}_{\mathrm{p}}.$ (25)

其中:Mp为动平台的广义质量矩阵, Cp为广义阻尼矩阵,Gp为广义重力。

每根钢索的第n个节点与动平台铰接,该节点和动平台所受拉力大小相等方向相反,则钢索与动平台存在约束,表示如下:

$\boldsymbol{T}_1=\boldsymbol{e}_n^j-\boldsymbol{x}_{\mathrm{e} 1}-\boldsymbol{R}\left(\boldsymbol{x}_{\mathrm{e} 2}\right) \boldsymbol{b}_j={{\bf{0}}}, j=1, 2, \cdots, 6 .$ (26)

其中:enj为第j根钢索的第n个节点坐标, bj为铰接点在局部坐标系中的位置矢量。

1.3 索牵引并联机器人

1.1和1.2节基于绝对节点坐标法对钢索单元进行了分析,得到式(23)所示的钢索单元的动力学方程。对划分的钢索单元进行组装,即可得到整根钢索的动力学方程。对于六索牵引的并联机器人,六根钢索分别划分为nj(j=1, 2,…, 6)个单元,共有nj+1个节点,每个节点坐标eij包含6个广义位移分量,第j根钢索的节点坐标矢量为

$\boldsymbol{q}_j=\left[\boldsymbol{e}_0^j, \cdots, \boldsymbol{e}_i^j, \cdots, \boldsymbol{e}_{n_j}^j\right]^{\mathrm{T}}.$ (27)

其中eij表示第j根钢索中i(i=0, n1,…, nj)节点广义位移分量。

引入只包含0和1元素的Boole矩阵B,则整根钢索的节点坐标矢量和钢索单元节点坐标矢量转换关系为

$\boldsymbol{u}_j^i=\boldsymbol{B}^i \boldsymbol{q}_j,$ (28)
$\boldsymbol{B}^i=\left[{\bf{0}}_{12 \times 6}, \cdots, \boldsymbol{I}_{12 \times 12}, \cdots, {\bf{0}}_{12 \times 6}\right] .$ (29)

其中I12×12为单位方阵。整根钢索惯性力虚功

$\begin{aligned} \boldsymbol{\delta} W_I^j=& \sum\limits_{i=1}^{n_j} \boldsymbol{\delta} W_{\mathrm{eI}}^i=\sum\limits_{i=1}^{n_j} \boldsymbol{\delta}\left(\boldsymbol{u}^i\right)^{\mathrm{T}} l_{\mathrm{e}}^i \int_0^1 \rho A \boldsymbol{S}^{\mathrm{T}} \ddot{\boldsymbol{r}} \mathrm{d} \xi=\\ & \boldsymbol{\delta}\left(\boldsymbol{u}_j\right)^{\mathrm{T}}\left(\boldsymbol{M}_I^j \ddot{\boldsymbol{q}}_j+\boldsymbol{C}_{\mathrm{I}}^j \dot{\boldsymbol{q}}_j+\boldsymbol{K}_I^j \boldsymbol{q}_j\right), \end{aligned}$ (30)
$\boldsymbol{M}_{\mathrm{1}}^j=\sum\limits_{i=1}^{n_j}\left(\boldsymbol{B}^i\right)^{\mathrm{T}} \boldsymbol{M}_{\mathrm{eI}}^i \boldsymbol{B}^i, $ (31)
$\boldsymbol{C}_{\mathrm{I}}^j=\sum\limits_{{I}=1}^{n_j}\left(\boldsymbol{B}^i\right)^{\mathrm{T}} \boldsymbol{C}_{\mathrm{eI}}^i \boldsymbol{B}^i, $ (32)
$\boldsymbol{K}_{\mathrm{I}}^j=\sum\limits_{i=1}^{n_j}\left(\boldsymbol{B}^i\right)^{\mathrm{T}} \boldsymbol{K}_{\mathrm{eI}}^i \boldsymbol{B}^i.$ (33)

其中:MIj为整根钢索的广义质量矩阵,CIj为广义阻尼矩阵,KIj为广义刚度矩阵。

同理,得到整根钢索的广义弹性力和重力为:

$\boldsymbol{Q}_{\mathrm{E}}^j=\sum\limits_{i=1}^{n_j}\left(\boldsymbol{B}^i\right)^{\mathrm{T}} \boldsymbol{Q}_{\mathrm{eE}}^i,$ (34)
$\boldsymbol{Q}_{\mathrm{F}}^j=\sum\limits_{i=1}^{n_j}\left(\boldsymbol{B}^i\right)^{\mathrm{T}} \boldsymbol{Q}_{\mathrm{eF}}^i .$ (35)

整根钢索的动力学方程为

$\boldsymbol{M}_1^j \ddot{\boldsymbol{q}}_j+\boldsymbol{C}_1^j \dot{\boldsymbol{q}}_j+\boldsymbol{K}_{\mathrm{I}}^j \boldsymbol{q}_j+\boldsymbol{Q}_{\mathrm{E}}^j=\boldsymbol{Q}_{\mathrm{F}}^j .$ (36)

每根钢索的第0个节点与塔顶定滑轮接触,将该节点视为定点,则该点的约束方程为

$\boldsymbol{T}_2=\left[e_0^j(1), e_0^j(2), e_0^j(3)\right]^{\mathrm{T}}-\boldsymbol{P}_j={\bf{0}}.$ (37)

其中Pj为出索点在全局坐标系中的位置矢量。

将钢索的动力学方程与动平台的动力学方程进行组装,六索牵引并联机器人的整体坐标矢量包括6根钢索的节点e和动平台的位姿xe,可表示为

$\begin{gathered} \boldsymbol{q}=\left[\boldsymbol{e}_0^1, \cdots, \boldsymbol{e}_{n 1}^1, \cdots, \boldsymbol{e}_i^j, \cdots, \right. \\ \left.\boldsymbol{e}_0^6, \cdots, \boldsymbol{e}_{n 6}^6, \boldsymbol{x}_e^{\mathrm{T}}\right]^{\mathrm{T}} . \end{gathered}$ (38)

引入Lagrange乘子λ,六索牵引并联机器人的动力学方程组可表示为

$\left\{\begin{array}{l} \boldsymbol{M} \ddot{\boldsymbol{q}}+\boldsymbol{C} \dot{\boldsymbol{q}}+\boldsymbol{K} \boldsymbol{q}+\boldsymbol{T}_{\mathrm{q}}^{\mathrm{T}} \lambda=\boldsymbol{Q}, \\ \boldsymbol{T}=\left[\boldsymbol{T}_1, \boldsymbol{T}_2\right]^{\mathrm{T}}={\bf{0}} . \end{array}\right.$ (39)

其中:M为索牵引并联机构的广义质量矩阵,C为广义阻尼矩阵,K为广义刚度矩阵,Q为广义力,T为全部约束方程,Tq为约束方程对整体坐标矢量的偏导。

以FAST馈源支撑系统中的六索牵引并联机器人为研究对象,分别建立钢索的动力学方程和馈源舱的刚性动力学方程,并基于钢索与馈源舱铰接点的约束和钢索出索点的约束建立了如式(39)所示的六索牵引并联机器人的动力学方程组。本文使用MATLAB中求解刚性微分方程的ode23tb,即向后差分法进行计算。

2 数值算例

为了验证所建立模型的可行性,使用FAST馈源支撑系统的主要参数建立动力学方程组,参数如表 1所示。

表 1 FAST馈源支撑系统主要参数
参数描述 取值
支撑塔分布圆半径/m 300
馈源焦面曲率半径/m 160
支撑塔距原点垂直距离/m 28
馈源舱直径/m 13
馈源舱质心与铰接点垂直距离/m 0.246
两端点间钢索长度/m 140~420
钢索截面直径/mm 50
馈源舱自重/t 30
钢索单位自重/(N·m-1) 68
钢索弹性模量/GPa 150

馈源舱在FAST跟踪模式中缓慢匀速运动,在不加入闭环控制的条件下,馈源舱在初始时刻t0的期望位姿为xe_exp(t0)=[5.077 0,-3.020 2×10-2,-159.878 0,3.932 5×10-7,1.189 2×10-2,-6.613 1×10-5]Ttd时刻终值为xe_exp(td)= [5.182 1,-3.182 0×10-2,-159.873 9,4.244 2×10-7,1.213 9×10-2,-6.992 5×10-5]T,运行时间为10 s,根据钢索直线模型可求得钢索长度的变化ΔL。在初始时刻t0,期望位姿与观测位姿存在误差,观测位姿为xe_obs(t0)=[4.966 6,-3.682 9×10-2,-159.869 0,4.782 9×10-4,1.215 8×10-2,-4.116 4×10-3]T,在动平台位姿为xe_obs(t0)时,根据钢索的抛物线模型[15]计算钢索长度,6根钢索的初始长度分别为318.8,317.7,320.9,325.2,326.3和323.2 m。

将抛物线在全局坐标系OXY平面的投影等分为11个单元,计算相对应的钢索长度,即图 2中的钢索划分为11个单元,nj(j=1, 2,…, 6),得到钢索各节点的广义坐标初值。单根钢索与塔顶定滑轮接触的第1个钢索单元长度随时间变化,第0个节点的坐标可由表 1中参数计算得出。设置长度阈值lmax=45 m,当第1个钢索单元的长度大于lmax时,增加1个节点,即增加1个钢索单元;lmin=15 m,当第1个钢索单元的长度小于lmin时,减少1个节点,其余钢索单元的长度保持不变。以ΔL为输入值,计算索牵引并联机器人动力学方程组,得到动平台的位姿。钢索长度变化缓慢时动平台的期望位姿、观测位姿和模型仿真位姿对比如图 4所示。

图 4 钢索长度变化缓慢动平台位姿对比

由于初始观测位姿的测量误差、钢索节点位移速度的初始值误差,以及钢索弹性振动的存在导致本文中模型仿真结果的振动始终存在。如果将仿真的初始位姿设为期望值的初始位姿,则表现为围绕期望位姿的振动。上述仿真中动平台的运动速度缓慢,约为10 mm/s,钢索长度变化也十分缓慢,即形函数S对时间的一阶导数和二阶导数均可忽略不计,因此选取收放索较快的数据再次进行计算。

在FAST实际运行中,馈源支撑系统加入了闭环控制,此时钢索长度的变化不能直接通过钢索直线模型计算得到。将开环运行中$ \Delta \dot l$和观测钢索长度变化速度进行线性拟合,发现根据直线模型得到的索长变化速度与观测值极为相近。选取闭环控制中观测转速相对较大的数据代入本文所提出的模型再次进行仿真计算,图 5为观测速度的3次插值。在初始时刻t0,动平台的观测位姿为xe_obs(t0)=[13.983 3,2.747 1,-159.319 0,-9.486 2×10-4,-3.693 6×10-2,1.214 7×10-2]T,6根钢索的初始长度分别为311.0,310.3,321.5,323.7,333.3和322.8 m。

图 5 闭环控制中观测转速

钢索长度变化较快时动平台的期望位姿、观测位姿和模型仿真位姿对比如图 6所示。由图 6可知,FAST在实际运行中加入闭环控制后,动平台观测位姿与期望位姿误差明显减小,振动明显减小。而仿真中振动依然存在且振动误差有所增大,模型仿真得到的动平台位姿变化仍与期望位姿趋势和观测位姿趋势相同。

图 6 钢索长度变化较快动平台位姿对比

3 结论

本文针对FAST馈源支撑系统中的六索牵引并联机器人进行动力学分析,首先使用绝对节点坐标法建立了柔性钢索的模型,在此基础上根据牵引索末端与馈源舱的耦合关系,推导出了索牵引并联机器人的动力学模型。将FAST实际运行中钢索的长度、速度和加速度的数据代入所建立动力学方程组中,对建立的动力学模型进行实验验证。对建立的方程组进行数值求解可以得到运动过程中钢索上节点位置的广义向量以及馈源舱的位姿。在建模过程中,为了方便计算,简化了一系列问题,例如使用钢索的抛物线模型,忽略了钢索上不均匀分布的传输光缆,将塔顶定滑轮与钢索的接触点视为定点等,导致所建立的动力学模型与实际模型存在一定误差。但用本文建立的动力学方程组仿真得到的馈源舱位姿与实际运行中馈源舱位姿进行对比,仍然有较好的吻合性,验证了该建模方法的可行性与有效性,对刚-柔耦合的钢索系统建模具有一定的参考价值。

参考文献
[1]
NAN R D. Five-hundred-meter aperture spherical radio telescope (FAST)[J]. Science in China Series G: Physics, Mechanics & Astronomy, 2006, 49(2): 129-148.
[2]
黄亮, 朱文白, 唐晓强. 大射电望远镜巨型柔性并联机构悬索分析及简化[J]. 天文研究与技术, 2013, 10(1): 77-84.
HUANG L, ZHU W B, TANG X Q. Analysis and a simplified model of the cable of a large-scale flexible parallel manipulator for a large radio telescope[J]. Astronomical Research & Technology, 2013, 10(1): 77-84. (in Chinese)
[3]
YIN J N, JIANG P, YAO R. An approximately analytical solution method for the cable-driven parallel robot in FAST[J]. Research in Astronomy and Astrophysics, 2021, 21(2): 46. DOI:10.1088/1674-4527/21/2/46
[4]
李辉, 孙京海, 朱文白, 等. 500m口径球面射电望远镜柔性馈源支撑系统仿真[J]. 计算机辅助工程, 2011, 20(1): 106-112.
LI H, SUN J H, ZHU W B, et al. Simulation on flexible feed support system of five-hundred-meter aperture spherical radio telescope[J]. Computer Aided Engineering, 2011, 20(1): 106-112. (in Chinese)
[5]
DU J L, BAO H, CUI C Z, et al. Dynamic analysis of cable-driven parallel manipulators with time-varying cable lengths[J]. Finite Elements in Analysis and Design, 2012, 48(1): 1392-13. DOI:10.1016/j.finel.2011.08.012
[6]
SHABANA A A. Flexible multibody dynamics: Review of past and recent developments[J]. Multibody System Dynamics, 1997, 1(2): 189-222. DOI:10.1023/A:1009773505418
[7]
吴懋琦, 谭述君, 高飞雄. 基于绝对节点坐标法的平面梁有限变形下变形重构[J]. 力学学报, 2021, 53(10): 2776-2789.
WU M Q, TAN S J, GAO F X. Shape reconstruction of plane beam with finite deformation based on absolute nodal coordinate formulation[J]. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(10): 2776-2789. (in Chinese)
[8]
LUO C Q, SUN J L, WEN H, et al. Dynamics of a tethered satellite formation for space exploration modeled via ANCF[J]. Acta Astronautica, 2020, 177: 882-890. DOI:10.1016/j.actaastro.2019.11.028
[9]
王忠民, 吴力国. 基于变长度单元ANCF的轴向伸展悬臂梁振动分析[J]. 振动与冲击, 2019, 38(3): 186-191.
WANG Z M, WU L G. Vibration analysis of axially deploying cantilever beam based on ANCF with length-varying beam element[J]. Journal of Vibration and Shock, 2019, 38(3): 186-191. (in Chinese)
[10]
BULÍN R, HAJŽMAN M, POLACH P. Nonlinear dynamics of a cable-pulley system using the absolute nodal coordinate formulation[J]. Mechanics Research Communications, 2017, 82: 21-28. DOI:10.1016/j.mechrescom.2017.01.001
[11]
FERRAVANTE V, RIVA E, TAGHAVI M, et al. Dynamic analysis of high precision construction cable-driven parallel robots[J]. Mechanism and Machine Theory, 2019, 135: 54-64. DOI:10.1016/j.mechmachtheory.2019.01.023
[12]
齐朝晖, 曹艳, 王刚. 多柔体系统数值分析的模型降噪方法[J]. 力学学报, 2018, 50(4): 863-870.
QI Z H, CAO Y, WANG G. Model smoothing methods in numerical analysis of flexible multibody systems[J]. Chinese Journal of Theoretical and Applied Mechanics, 2018, 50(4): 863-870. (in Chinese)
[13]
BERZERI M, SHABANA A A. Development of simple models for the elastic forces in the absolute nodal co-ordinate formulation[J]. Journal of Sound and Vibration, 2000, 235(4): 539-565.
[14]
WAGO T, KOBAYASHI N, SUGAWARA Y. Improvement on evaluating axial elastic force in Bernoulli-Euler beam based on the absolute nodal coordinate formulation by accurate mean axial strain measure[C]//Proceedings of ASME 2011 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference. Washington, USA: American Society of Mechanical Engineers, 2011: 963-970.
[15]
姚蕊, 唐晓强, 李铁民, 等. 大型射电望远镜馈源定位3T索牵引并联机构分析与设计[J]. 机械工程学报, 2007, 43(11): 105-109.
YAO R, TANG X Q, LI T M, et al. Analysis and design of 3T cable-driven parallel manipulator for the feedback's orientation of the large radio telescope[J]. Chinese Journal of Mechanical Engineering, 2007, 43(11): 105-109. (in Chinese)