CN113051678B - 一种基于线性插值形函数的柔性绳索建模方法 - Google Patents

Abstract

本发明提供了一种基于线性插值形函数的柔性绳索建模方法，该方法从计算时间和计算精度两方面考虑，克服了现有技术中如集中质量法的计算速度最快但精度较低、基于绝对节点坐标的非线性有限元方法的计算耗时长但精度较高等的缺陷，方法中采用两节点的非线性绳索单元具有适中的计算时间及较高的计算精度，兼具了计算时间经济性和高精度特性，能够满足于工程化应用需求，是柔性拦截网绳索建模较为合理的选择。

Description

一种基于线性插值形函数的柔性绳索建模方法

技术领域

本发明属于绳索建模技术领域，具体涉及一种基于线性插值形函数的柔性绳索建模方法。

背景技术

柔性绳索具有质量轻、强度大和易折叠储存等优点。由柔性绳索编织而成的柔性拦截网展开面积大且极易与目标发生钩挂和缠绕，对目标外形及运动状态的适应性远远大于常规捕获方式，柔性拦截网在低慢小类无人机目标抓捕方面获得了广泛应用。通过建模的方式有助于对绳索的横向、纵向运动以及局部受力等进行动力学求解分析，但由于柔性绳索通常由有机纤维编织而成，其内部结构十分复杂，如采用过于精细化的建模方法会耗费大量的计算时间并且在精度提高方面成效不够显著。目前，现有技术中常用的绳索建模方法有弹簧阻尼模型和基于绝对节点坐标的非线性有限元单元，前者耗时短但精度相对较差，后者精度虽然很高但是耗时久，对计算资源要求高，难以满足工程化应用需求。

发明内容

为克服现有技术中的问题，本发明提供了一种基于线性插值形函数的柔性绳索建模方法，其具体步骤为：

第一步，以某个预定时刻的绳索构型作为参考构型空间，标识当前时刻的绳索构型上任意绳索点在空间中的位置；

第二步，定义绳索的格林应变张量；

第三步，计算当前时刻绳索构型上任意绳索点的应力；

第四步，基于所述参考构型、格林应变张量以及应力，并考虑绳索的运动情况、材料以及受力，建立绳索构型的控制方程；

第五步，确定所述格林应变张量与所述应力之间的关系；

第六步，选取绳索上的某段作为一个单元，基于线性插值函数确定所述单元两端节点间的绳索上任意一点的位移；

第七步，利用所述第二步得到的格林应变张量以及第五步得到的所述关系，确定绳索单元的应变能；

第八步，针对绳索单元建立两节点有限元动力学方程；

第九步，利用数值仿真算法求所述两节点有限元动力学方程，实现对柔性绳索的动力学分析。

进一步地，所述第一步中，以某个预定时刻作为初始时刻t＝0，令绳索在此时的外形作为参考构型空间Ω₀，在当前时刻t的变形构型域为Ω_t，绳索变形域的边界为Γ，采用三维空间单位正交矢量e₁、e₂和e₃作为参考系。在参考构型中，以任意绳索点的位置矢量X作为材料坐标，该坐标提供了材料点的标识。以材料坐标X和时间t作为函数变量，可将任意材料点在空间中的位置标识为：

式中，X为参考构形中材料点的位置矢量；x为当前时刻构型中材料点的位置矢量；φ(X,t)是从参考构形到当前时刻构型的映射函数；e₁、e₂和e₃为三维空间参考系的单位正交基矢量。

进一步地，所述第二步中，考虑在初始构型中的一个无限小段dX，采用了爱因斯坦求和约定(Einstein summation convention)，即若同一指标成对出现，则遍历其取值范围求和，因此，无限小段dX在当前时刻构型下对应微段dx可表示为：

或/>

考虑度量物体的变形特征，基于连续介质力学理论，运动的变形梯度F为

格林(Green)应变张量E定义为：

dx²-dX²＝2dX·E·dX (4)

结合式(3)与式(4)可得：

在几何非线性力学问题中，具有多种应力度量，由于第二类柯西(Poila-Kirchhoff)应力S与格林应变张量E在能量上是耦合的，故在所述第三步中，具体选取第二类柯西应力S作为应力度量。在当前时刻构形下变形体某一微元截面上内力的合力为df，截面积为dΓ，当前状态下的变形梯度为F，微元截面在初始构形下的截面积为dΓ₀，微元截面在初始时刻的外法线单位矢量为n₀，则定义第二类柯西应力S为：

n₀SdΓ₀＝F^-1df (6)

进一步地，所述第四步建立绳索构型的控制方程具体过程为：

在当前构形上建立平衡方程，转换至参考构形后可得变形体控制方程为：

其中，S_ij是第二类柯西应力S的分量，δE_ij是格林应变张量E的变分分量；Ω₀是初始时刻参考构型空间，是加速度矢量分量，ρ为材料密度，p_idS₀为作用在变形体表面的外力。

进一步地，所述第五步中确定格林应变张量与应力之间关系具体为：

从一般性考虑，应力与应变的关系为：

S_ij＝C_ijklE_kl (8)

其中，C_ijkl为弹性模量的四阶张量；E_kl是格林应变张量E的分量。

在绳索单元中，当仅考虑单元的轴向力时，应力应变关系可简化为：

S₁₁＝E_cE₁₁ (9)

其中，E_c为绳索轴向拉伸杨氏模量；E₁₁是格林应变张量E的x方向的正应变。

进一步地，所述第六步具体包括：

设单元上两节点间的插值形函数是线性的，选取单元节点i的位置坐标r_i和节点j的位置坐标r_j，令w＝(r_i r_j)^T，在参考构形上距离节点(Node)i为s的点，其在当前时刻的位移为：

其中，r(s)是s点的空间位置；s是物质点的参考构型坐标；L为单元的长度；I₃是3阶单位阵，w是单元节点坐标向量；N(s)为单元的插值形函数。

进一步地，所述第七步中确定应变能具体包括：

记r(s)＝(x y z)^T，则单元的变形梯度可以表示为：

结合式(5)与式(8)得到，单元在r(s)点处的格林应变张量为：

令N_i(s)为插值形函数矩阵N(s)的第i行(i＝1,2,3)，则结合式(9)可得到格林应变张量E的x方向的正应变为：

其中，转换矩阵

进一步地，所述第八步中两节点有限元动力学方程的建立过程具体包括：

仅考虑绳索轴向变形产生的应变能，结合应力应变关系式(8)，并代入控制方程式(7)，可得到：

其中，E_c是杨氏模量；w是单元节点坐标向量；是单元节点的加速度向量；δw表示单元节点坐标的变化量；B是一个6阶转换矩阵；s是物质点的参考构型坐标；ρ为材料密度；N为单元的形函数；p是外部作用力向量；

消去变分项，可得线性插值形函数下的两节点有限元动力学方程为：

进一步地，所述第九步中具体采用显式纽马克(Newmark)方法求解动力学方程：

其中，α和β为纽马克系数；Δt为时间步长；u、分别表示单元在当前时刻的位移和速度；u_p、/>和/>分别表示单元在前一时刻结束时的位移、速度和加速度。

至此，完成了基于线性插值形函数的柔性绳索建模方法，该方法可用于柔性拦截网绳索的动力学求解。

有益效果

上述本发明所提供的基于线性插值形函数的柔性绳索建模方法，从计算时间和计算精度两方面考虑，相较于集中质量法的计算速度最快但精度较低、基于绝对节点坐标的非线性有限元方法的计算耗时长但精度较高，两节点的非线性绳索单元具有适中的计算时间及较高的计算精度，兼具了计算时间经济性和高精度特性，能够满足于工程化应用需求，是柔性拦截网绳索建模较为合理的选择。

附图说明

图1示出了本发明的方法中所定义的初始时刻参考构型与当前时刻构型；

图2示出了变形体变形前后作用力的相对变化；

图3示出了本发明的方法中所选取的两节点线性绳索单元；

图4示出了基于本发明的实例与现有技术之间的技术效果对比。

具体实施方式

下文将结合具体实施例对本发明做更进一步的详细说明。应当理解，下列实施例仅为示例性地说明和解释本发明，而不应被解释为对本发明保护范围的限制。凡基于本发明上述内容所实现的技术均涵盖在本发明旨在保护的范围内。

本发明提供的一种基于线性插值形函数的柔性绳索建模方法，其具体步骤为：

第一步，以某个预定时刻的绳索构型作为参考构型空间，标识当前时刻的绳索构型上任意绳索点在空间中的位置；

第二步，定义绳索的格林应变张量；

第三步，计算当前时刻绳索构型上任意绳索点的应力；

第四步，基于所述参考构型、格林应变张量以及应力，并考虑绳索的运动情况、材料以及受力，建立绳索构型的控制方程；

第五步，确定所述格林应变张量与所述应力之间的关系；

第六步，选取绳索上的某段作为一个单元，基于线性插值函数确定所述单元两端节点间的绳索上任意一点的位移；

第七步，利用所述第二步得到的格林应变张量以及第五步得到的所述关系，确定绳索单元的应变能；

第八步，针对绳索单元建立两节点有限元动力学方程；

第九步，利用数值仿真算法求所述两节点有限元动力学方程，实现对柔性绳索的动力学分析。

在本发明的一个优选实施方式中，所述第一步中，以某个预定时刻作为初始时刻t＝0，令绳索在此时的外形作为参考构型空间Ω₀，在当前时刻t的变形构型域为Ω_t，绳索变形域的边界为Γ，采用三维空间单位正交矢量e₁、e₂和e₃作为参考系，如图1所示。在参考构型中，以任意绳索点的位置矢量X作为材料坐标，该坐标提供了材料点的标识。以材料坐标X和时间t作为函数变量，可将任意材料点在空间中的位置标识为：

式中，X为参考构形中材料点的位置矢量；x为当前时刻构型中材料点的位置矢量；φ(X,t)是从参考构形到当前时刻构型的映射函数；e₁、e₂和e₃为三维空间参考系的单位正交基矢量。

所述第二步中，考虑在初始构型中的一个无限小段dX，采用了爱因斯坦求和约定(Einstein summation convention)，即若同一指标成对出现，则遍历其取值范围求和，因此，无限小段dX在当前时刻构型下对应微段dx可表示为：

或/>

考虑度量物体的变形特征，基于连续介质力学理论，运动的变形梯度F为

格林(Green)应变张量E定义为：

dx²-dX²＝2dX·E·dX (4)

结合式(3)与式(4)可得：

在几何非线性力学问题中，具有多种应力度量，由于第二类柯西(Poila-Kirchhoff)应力S与格林应变张量E在能量上是耦合的，故在所述第三步中，具体选取第二类柯西应力S作为应力度量，如图2所示。在当前时刻构形下变形体某一微元截面上内力的合力为df，截面积为dΓ，当前状态下的变形梯度为F，微元截面在初始构形下的截面积为dΓ₀，微元截面在初始时刻的外法线单位矢量为n₀，则定义第二类柯西应力S为：

n₀SdΓ₀＝F^-1df (6)

所述第四步建立绳索构型的控制方程具体过程为：

在当前构形上建立平衡方程，转换至参考构形后可得变形体控制方程为：

其中，S_ij是第二类柯西应力S的分量，δE_ij是格林应变张量E的变分分量；Ω₀是初始时刻参考构型空间，是加速度矢量分量，ρ为材料密度，p_idS₀为作用在变形体表面的外力。

所述第五步中确定格林应变张量与应力之间关系具体为：

从一般性考虑，应力与应变的关系为：

S_ij＝C_ijklE_kl (8)

其中，C_ijkl为弹性模量的四阶张量；E_kl是格林应变张量E的分量。

在绳索单元中，当仅考虑单元的轴向力时，应力应变关系可简化为：

S₁₁＝E_cE₁₁ (9)

其中，E_c为绳索轴向拉伸杨氏模量；E₁₁是格林应变张量E的x方向的正应变。

如图3所示，所述第六步具体包括：

设单元上两节点间的插值形函数是线性的，选取单元节点i的位置坐标r_i和节点j的位置坐标r_j，令w＝(r_i r_j)^T，在参考构形上距离节点(Node)i为s的点，其在当前时刻的位移为：

其中，r(s)是s点的空间位置；s是物质点的参考构型坐标；L为单元的长度；I₃是3阶单位阵，w是单元节点坐标向量；N(s)为单元的插值形函数。

所述第七步中确定应变能具体包括：

记r(s)＝(x y z)^T，则单元的变形梯度可以表示为：

结合式(5)与式(8)得到，单元在r(s)点处的格林应变张量为：

令N_i(s)为插值形函数矩阵N(s)的第i行(i＝1,2,3)，则结合式(9)可得到格林应变张量E的x方向的正应变为：

其中，转换矩阵

所述第八步中两节点有限元动力学方程的建立过程具体包括：

仅考虑绳索轴向变形产生的应变能，结合应力应变关系式(8)，并代入控制方程式(7)，可得到：

其中，E_c是杨氏模量；w是单元节点坐标向量；是单元节点的加速度向量；δw表示单元节点坐标的变化量；B是一个6阶转换矩阵；s是物质点的参考构型坐标；ρ为材料密度；N为单元的形函数；p是外部作用力向量；

消去变分项，可得线性插值形函数下的两节点有限元动力学方程为：

所述第九步中具体采用显式纽马克(Newmark)方法求解动力学方程：

其中，α和β为纽马克系数；Δt为时间步长；u、分别表示单元在当前时刻的位移和速度；u_p、/>和/>分别表示单元在前一时刻结束时的位移、速度和加速度。

表1和图1示出了基于本发明的柔性绳索建模方法，与对基于绝对节点坐标的非线性有限元建模方法(Absolute coordinate nodal method，ANCF)、基于线性插值形函数非线性绳索单元建模方法(Nonliner cable element，NCE)、集中质量建模方法(Consentrated mass model，CMM)等现有技术在计算时间上的对比结果，可以看出本发明所提供的建模方式在耗时与动力学响应方面均优于上述现有方法。

表1本发明所提供方法与现有技术在自由度相同情况下的计算时间对比

以上对本发明的实施方式进行了说明。但是，本发明不限定于上述实施方式。凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种基于线性插值形函数的柔性绳索建模方法，其特征在于：其具体步骤为：

第一步，以某个预定时刻的绳索构型作为参考构型空间，标识当前时刻的绳索构型上任意绳索点在空间中的位置；

第二步，定义绳索的格林应变张量；

第三步，计算当前时刻绳索构型上任意绳索点的应力；

第四步，基于所述参考构型、格林应变张量以及应力，并考虑绳索的运动情况、材料以及受力，建立绳索构型的控制方程；

第五步，确定所述格林应变张量与所述应力之间的关系；

第六步，选取绳索上的某段作为一个单元，基于线性插值函数确定所述单元两端节点间的绳索上任意一点的位移；

第七步，利用所述第二步得到的格林应变张量以及第五步得到的所述关系，确定绳索单元的应变能；

第八步，针对绳索单元建立两节点有限元动力学方程；

第九步，利用数值仿真算法求所述两节点有限元动力学方程，实现对柔性绳索的动力学分析；

所述第一步中，以某个预定时刻作为初始时刻t＝0，令绳索在此时的外形作为参考构型空间Ω₀，在当前时刻t的变形构型域为Ω_t，绳索变形域的边界为Γ，采用三维空间单位正交矢量e₁、e₂和e₃作为参考系；在参考构型中，以任意绳索点的位置矢量X作为材料坐标，以材料坐标X和时间t作为函数变量，则任意材料点在空间中的位置标识为：

式中，X为参考构形中材料点的位置矢量；x为当前时刻构型中材料点的位置矢量；φ(X,t)是从参考构形到当前时刻构型的映射函数；e₁、e₂和e₃为三维空间参考系的单位正交基矢量；

所述第二步中，考虑在初始构型中的一个无限小段dX，采用爱因斯坦求和约定将无限小段dX在当前时刻构型下对应微段dx表示为：

考虑度量物体的变形特征，基于连续介质力学理论，运动的变形梯度F为：

格林应变张量E定义为：

dx²-dX²＝2dX·E·dX (4)

结合式(3)与式(4)可得：

所述第四步建立绳索构型的控制方程具体过程为：

在当前构形上建立平衡方程，转换至参考构形后得到变形体控制方程为：

其中，S_ij是第二类柯西应力S的分量，δE_ij是格林应变张量E的变分分量；Ω₀是初始时刻参考构型空间，是加速度矢量分量，ρ为材料密度，p_idS₀为作用在变形体表面的外力；

所述第八步中两节点有限元动力学方程的建立过程具体包括：

仅考虑绳索轴向变形产生的应变能，结合应力应变关系式(8)，并代入控制方程式(7)，可得到：

其中，E_c是杨氏模量；w是单元节点坐标向量；是单元节点的加速度向量；δw表示单元节点坐标的变化量；B是一个6阶转换矩阵；s是物质点的参考构型坐标；ρ为材料密度；N为单元的形函数；p是外部作用力向量；

消去变分项，可得线性插值形函数下的两节点有限元动力学方程为：

2.如权利要求1所述的方法，其特征在于：所述第三步中，具体选取第二类柯西应力S作为应力度量；在当前时刻构形下变形体某一微元截面上内力的合力为df，截面积为dΓ，当前状态下的变形梯度为F，微元截面在初始构形下的截面积为dΓ₀，微元截面在初始时刻的外法线单位矢量为n₀，则定义第二类柯西应力S为：

n₀SdΓ₀＝F^-1df。 (6)

3.如权利要求1所述的方法，其特征在于：所述第五步中确定格林应变张量与应力之间关系具体为：

从一般性考虑，应力与应变的关系为：

S_ij＝C_ijklE_kl (8)

其中，C_ijkl为弹性模量的四阶张量；E_kl是格林应变张量E的分量；

在绳索单元中，当仅考虑单元的轴向力时，应力应变关系可简化为：

S₁₁＝E_cE₁₁ (9)

其中，E_c为绳索轴向拉伸杨氏模量；E₁₁是格林应变张量E的x方向的正应变。

4.如权利要求3所述的方法，其特征在于：所述第六步具体包括：

设单元上两节点间的插值形函数是线性的，选取单元节点i的位置坐标r_i和节点j的位置坐标r_j，令w＝(r_i r_j)^T，在参考构形上距离节点i为s的点，其在当前时刻的位移为：

其中，r(s)是s点的空间位置；s是物质点的参考构型坐标；L为单元的长度；I₃是3阶单位阵，w是单元节点坐标向量；N(s)为单元的插值形函数。

5.如权利要求4所述的方法，其特征在于：所述第七步中确定应变能具体包括：

记r(s)＝(x y z)^T，则单元的变形梯度可以表示为：

结合式(5)与式(8)得到，单元在r(s)点处的格林应变张量为：

令N_i(s)为插值形函数矩阵N(s)的第i行(i＝1,2,3)，则结合式(9)可得到格林应变张量E的x方向的正应变为：

其中，转换矩阵

6.如权利要求1所述的方法，其特征在于：所述第九步中具体采用显式纽马克方法求解动力学方程：

其中，α和β为纽马克系数，Δt为时间步长，u_p、和/>分别表示单元在前一时刻结束时的位移、速度和加速度；u、/>分别表示单元在当前时刻的位移和速度。