一种拖曳系统非线性动力学响应求解方法

1.本发明涉及柔性缆索拖曳系统，具体涉及一种拖曳系统非线性动力学响应求解方法。


背景技术：

2.在工程结构中，空间柔性缆索结构是一种典型构件，从物理角度看，其纵向尺寸远远大于横向尺寸，表现出大柔性。从力学角度看，柔性索承受拉伸荷载,抗压缩、弯曲、剪切及扭转的能力极低。柔性缆索结构耗材少，强度高、低柔度等特性，极易连接远距离的两个系统。作为信号传输和结构牵引的媒介，已经广泛应用于空间绳系卫星系统，空中加油系统，无人机、微型飞行器系统，无人遥控潜水器系统，海洋工程，以及大跨度悬索桥结构，其应用范围持续扩大，与现代科技联系紧密。柔性缆索系统动力学受工作环境的影响显著，呈现强非线性特征，动力学响应行为复杂。理论解析法的应用，只限于非常简单的工况和准静态情形。在大多数情况下，尚没有一套完整的理论解体系，可以综合介质阻力效应和缆索系统运动，主要采用数值方法进行问题的求解。缆索动力学行为的研究，一直是一个热点问题。研究缆索动力学的方法中，有限单元法是基于经典的虚功原理，对缆索微元从虚功平衡方程进行推导，构造出缆索系统的动力学控制方程，控制方程为二阶常微分方程组。对于不同类型的结构单元，通过选取适当阶次的插值形函数，构造相应形式的虚功原理平衡方程，并利用数值算法进行积分求解。有限元方法具有处理复杂边界条件和不同索元材料特性的优势，被广泛应用于柔性缆索拖曳系统的分析，但是该方法求解结果无法良好解决长时程积分当中的能量耗散问题。传统的数值积分方法(如龙格-库塔法，newmark-法等)能够有效预测缆索结构的动力学响应，但是，对于经历长时程、大位移运动的柔性缆索结构，因为低频刚体运动和高频的弹性变形的长期相耦合，使得传统积分算法逐渐出现数值不稳定性。
3.2019年，阚子云等人公开了《一种基于多体系统滑移绳索单元的聚合式张拉整体结构动力响应分析方法》，在多体系统的传统绳索单元基础上，建立多体系统滑移绳索单元，求解多体动力学微分代数方程组，以获得聚合式张拉整体结构的动力响应，此方法无法克服求解大位移动力学控制方程中的非线性问题；g li等人于《characteristics of coupled orbital-attitude dynamics offlexible electric solar wind sail》一文中针对太阳风帆的绳索拖曳问题，采用有限元法进行了缆索系统的动力学响应，但是求解结果无法良好解决长时程积分当中的能量耗散问题。


技术实现要素：

4.本发明的目的在于提供一种拖曳系统非线性动力学响应求解方法，能够对柔性缆索系统在流体介质中的动力学性能进行详细准确的求解，有效降低长时程计算分析中带来的累积误差，提高求解精度。
5.实现本发明目的的技术方案为：一种拖曳系统非线性动力学响应求解方法，包括步骤：
6.时间和计算步数初始化，初始化值为0；速度与位移数值初始化，初始化值为给定的初始值；
7.将时间步i的速度位移值作为时间步i+1的速度位移值进行下一时间步的初始化；
8.基于节点坐标有限元法构建拖曳系统单元阵；
9.基于欧拉角变换构建从局部坐标转换到整体坐标的转换矩阵，确定整体坐标系下的拖曳系统单元阵；
10.考虑外力确定外力阵，基于整体坐标系下拖曳系统单元阵，构建拖曳系统动力学控制方程；
11.通过辛差分算法进行动力学控制方程响应求解；
12.若求解的位移速度收敛性未达到要求，则更新位移速度阵，继续通过辛差分算法进行动力学控制方程响应求解；否则判别计算时间是否达到设定的总计算时长，若达到总计算时长，则输出结果，求解结束；若未达到总计算时长，则进入下一时间步进行新一轮计算。
13.进一步的，所述基于节点坐标有限元法构建拖曳系统单元阵具体包括：
14.基于节点坐标有限元法构建单元插值形函数n；
15.以节点坐标为变量表示单元的能量，构造单元矩阵。
16.进一步的，所述单元矩阵包括：
17.单元质量阵me：
[0018][0019]
其中，l是单元长度，ρ是单元的材料密度，a是单元的截面积；
[0020]
单元刚度阵ke：
[0021][0022]
其中，l0是单元初始长度，θ
x
、θy、θz分别是局部坐标系和全局坐标系x、y、z轴的夹角，即欧拉角，e是单元杨氏模量；
[0023]
单元位移阵r：
[0024]
r＝nxe[0025]
其中，r是单元上任意点的位置向量，xe是单元两个节点的位置向量；
[0026]
单元阻尼阵ce：
[0027]ce
＝αme+βke[0028]
其中，α,β为阻尼系数。
[0029]
进一步的，基于欧拉角变换构建从局部坐标转换到整体坐标的转换矩阵为：
[0030]
通过三个欧拉角θ
x-滚动，θ
y-俯仰，θ
z-偏航描述局部坐标系的基矢量(e
x
,ey,ez)；从全局坐标系到局部坐标系的转换顺序分别定义为偏航、俯仰和滚动；转换分为三个步骤
进行：
[0031]
(1)逆时针旋转oz轴，角度为θz，0≤θz《π，使oy轴和局部坐标系on轴重合；
[0032]
(2)顺时针旋转新的oy轴，角度为θy，0≤θy《π，使整体坐标系中的ox轴和局部坐标系的ox轴重合；
[0033]
(3)逆时针旋转新的ox轴，角度为θz，0≤θz《π，使整体坐标系中的oy轴和局部坐标系的oy轴重合，同时，oz轴和局部坐标系oz轴重合；
[0034]
通过以上旋转，得到整体坐标系到局部坐标系的转换矩阵t。
[0035]
进一步的，所述外力阵f为：
[0036]
f＝f
ge
+f
de
[0037][0038][0039][0040]
其中，f
ge
为整体坐标系下单元重力矢量，f
de
为整体坐标系下单元阻力，cd是阻力系数，ρf是介质密度，l是单元长度，ng为线性插值函数，vf是该点的索结构的介质流速。
[0041]
进一步的，所述线性插值函数为：
[0042][0043]
其中，ξ是单元局部坐标。
[0044]
进一步的，所述构建拖曳系统动力学控制方程包括：
[0045]
建立拖曳系统单元的哈密尔顿正则方程；
[0046]
基于节点位移，速度的连续性条件，将整体坐标系中单元阵组装成整体拖曳系统的动力学控制方程。
[0047]
进一步的，所述拖曳系统单元的哈密尔顿正则方程为：
[0048][0049]
其中，单元的共轭动量单元的速度矢量为单元的动能为he为单元的哈密尔顿函数，f
ke
和f
ge
分别是单元两节点内力阵向量和重力向量。
[0050]
进一步的，所述拖曳系统动力学控制方程为：
[0051][0052]
其中，m是系统质量矩阵，k是系统刚度矩阵，c是系统阻尼矩阵，fk是系统内力矩阵，fg是系统重力矩阵。
[0053]
进一步的，所述辛差分算法将哈密尔顿正则方程在k+1时间步分解，得到矩阵为：
[0054][0055][0056]
其中，h是计算时间步。
[0057]
本发明与现有技术相比，其显著效果为：
[0058]
(1)采用节点坐标有限元法，相比于传统有限元方法，采用全量格式，表达简单、求解精度高，有效改进了传统有限元方法冗杂的求解过程；
[0059]
(2)采用辛差分算法求解缆索系统的动力学控制方程，相比于传统积分算法，避免了随时间的逐渐增大，出现的累积计算误差会明显降低动力学响应的预测准确度的问题，提高了求解精度。
附图说明
[0060]
图1为本发明基于节点坐标有限元法和辛差分算法构建拖曳系统矩阵和动力学控制方程的推导流程图。
[0061]
图2为本发明拖曳系统非线性动力学响应求解方法的流程图。
具体实施方式
[0062]
下面将结合附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0063]
为了有效求解绳系拖曳系统在流体介质中的动力学响应，本发明结合了具有表达简单、求解精度高特点的节点坐标有限元法和在长时程仿真中保持高精度的辛差分算法。包含如下步骤：
[0064]
结合图1，本发明所提供的基于节点坐标有限元法和辛差分算法构建拖曳系统矩阵和动力学控制方程的方法，包括如下步骤：
[0065]
a1、基于节点坐标有限元法构建单元插值形函数
[0066]
具体的，考虑一两个节点的直线缆索单元，该单元的位置由其节点坐标来描述。缆索单元中任意点的位置，可以用整体框架中单元的形函数和单元两个节点的坐标来表示
[0067]
r＝nxe[0068]
其中，r是柔性索单元上任意点的位置向量，xe是单元两个节点的位置向量，n是插值形函数。相类似地，单元上任意点的速度、加速度，也将由相同的单元形函数和进行插值。
[0069][0070][0071]
其中，v和a分别是缆索单元上任意点的速度和加速度向量，和分别是单元两个节点的速度和加速度向量。
[0072]
a2、构造单元矩阵
[0073]
具体的，以节点坐标为变量表示单元的能量，分别有动能、弹性势能、重力势能。
[0074][0075][0076][0077]
其中，te和ue、u
ge
分别是单元动能、弹性势能、重力势能，l和ε分别是单元长度和应变，ρ和e是单元的材料密度和杨氏模量，a是单元的截面积，g是重力加速度。在上述表达式中，me和ke分别是单元的质量阵和刚度阵。f
ke
和f
ge
分别是单元两节点内力阵向量和重力向量。由单元能量的力学表达式可得单元的质量阵和刚度阵。
[0078][0079][0080]
其中，l0是单元初始长度，θ
x
、θy、θz分别是局部坐标系和全局坐标系x、y、z轴的夹角。
[0081]
a3、建立哈密尔顿正则方程
[0082]
具体的，定义单元的共轭动量为单元的速度矢量为单元的动能为
[0083]
单元的哈密尔顿函数为
[0084][0085]
可以得到缆索单元的哈密尔顿正则方程
[0086][0087]
其中，ce＝αme+βke，α,β为阻尼系数。
[0088]
a4、基于哈密尔顿系统采用辛差分方程
[0089]
具体的，缆索系统动力学在构型空间中的n个二阶微分方程被哈密顿形式简化为相空间中的2n个一阶正则方程。所得到的正则方程可以很自然地用辛积分格式积分求解。可以将哈密尔顿正则方程在k+1时间步分解，并得到如下的二阶辛方案。
[0090][0091][0092]
其中，h是计算时间步。h
x
和h
p
分别表示系统哈密尔顿函数对正则变量x,p的偏微分。单元矩阵形式为：
[0093][0094][0095]
a5、局部坐标到整体坐标的转换
[0096]
具体的，o-xyz定义了位于原点o的具有单位基向量(i,j,k)的全局笛卡儿惯性坐标系。对于任意单元e，局部坐标系o-xyz的原点位于该单元的起始点。局部坐标系的基矢量(e
x
,ey,ez)由三个欧拉角(θ
x-滚动，θ
y-俯仰，θ
z-偏航)来描述。从全局坐标系到局部坐标系的转换顺序分别定义为偏航、俯仰和滚动。转换分为三个步骤进行：
[0097]
(1)逆时针旋转oz轴，角度为θz(0≤θz《π)，使得oy轴和on重合。
[0098]
(2)顺时针旋转新的oy轴，角度为θy(0≤θy《π)，使得整体坐标系中的ox轴和局部坐标系的ox轴重合。
[0099]
(3)逆时针旋转新的ox轴，角度为θz(0≤θz《π)，使得整体坐标系中的oy轴和局部坐标系的oy轴重合，同时，oz轴和oz轴重合。
[0100]
通过以上旋转，可得整体坐标系到局部坐标系的转换矩阵t
[0101]
一旦系统整体坐标系到局部坐标系的转换矩阵确定，单元的动能te和弹性势能ue可以在整体坐标系中表示，在单元方程中施加位移、速度连续的边界约束，组装成缆索系统的整体质量阵和刚度阵。根据式的转换关系，局部坐标系下的单元两节点速度矢量可以转换为全局坐标系中的表达形式，如下所示
[0102][0103]
其中，ve是单元两节点速度矢量在整体坐标系中的表达形式，te是6
×
6的单元转换
矩阵。可得单元动能在整体坐标系下的表示
[0104][0105]
并通过此式求得整体坐标系下的单元质量阵
[0106][0107]
a6、考虑外力做功
[0108]
具体的，对于缆索结构来说，缆索结构本身的重量是不可忽略的。整体坐标系的oz轴沿竖直方向，并与重力方向相反。显然，在全局坐标系中可以方便地导出单元节点的重力矢量。假定未变形缆索结构自由下垂时，末节点所在的水平平面为零势能平面。单元的重力势能可以在整体坐标系下表示为
[0109][0110][0111]
其中，ng为线性插值函数
[0112][0113]
g,f
ge
和l
total
分别是重力加速度，单元重力矢量和初始索结构的总长度。ξ是单元局部坐标。xe是整体坐标系下的单元节点位置矢量。采用莫里森方程计算介质阻力，在考虑受阻面积的情况下，当直径尺寸远小于缆索单元长度时，可直接在全局坐标系中得到相应的单元阻力矢量。应用相同的单元线性插值形函数n，沿单元弧长插值索结构的速度矢量和介质的流速，可得
[0114]vf
＝nv
fe
[0115]
其中，和vf分别是该点的索结构的运动速度和介质流速。v
fe
是单元两节点处的介质流速矢量。由拖曳力在全局坐标下所做的虚功可将单元阻力式表示为
[0116][0117][0118]
其中，cd是阻力系数，ρf是介质密度，l是单元长度。
[0119]
a7、单元组装
[0120]
具体的，基于节点位移，速度的连续性条件，将系统单元阵组装成整体坐标系中的整体缆索系统的动力学控制方程。
[0121][0122]
m是系统质量矩阵，k是系统刚度矩阵，c是系统阻尼矩阵，fk是系统内力矩阵，fg是系统重力矩阵
[0123]
在图1构建的拖曳系统非线性动力学控制方程基础上，结合图2，对本发明拖曳系统非线性动力学响应求解方法进行描述，具体步骤如下：
[0124]
s1：对时间和计算步数进行初始化，初始化值为0；对速度与位移数值进行初始化，初始化值为赋予的初始值。
[0125]
s2:将时间步i的速度位移值作为时间步i+1的速度位移值进行下一时间步的初始化。
[0126]
s3：基于欧拉角变换构建从局部坐标转换到整体坐标的转换矩阵。
[0127]
s4:以节点坐标为变量来描述单元能量，进而进一步推导出单元的质量阵、刚度阵、位移阵、阻尼阵、外力阵。
[0128]
质量阵
[0129][0130]
刚度阵
[0131][0132]
位移阵
[0133]
r＝nxe[0134]
阻尼阵
[0135]ce
＝αme+βke[0136]
外力阵
[0137]
f＝f
ge
+f
de
[0138]
s5：根据不同单元在连接节点处的位移速度应保持一致的连续条件，将单元组装成缆索拖曳系统，并采用二阶辛差分算法进行求解。
[0139]
s6：对求解的位移速度结果收敛性进行判别，若收敛性达到要求，进入s7步骤；若未达到要求，则更新位移速度阵回到s5继续新一轮求解。
[0140]
s7：判别计算时间是否达到总计算时长，若达到总计算时长，则输出结果，求解结束；若未达到总计算时长，则进入下一时间步进行新一轮计算。