基于ANCF的薄膜太阳帆自旋展开动力学模拟

薛鹏聪，刘 铖,2,3，水小平

（1. 北京理工大学 宇航学院，北京 100081；2. 中物院高性能数值模拟软件中心，北京 100088；3. 北京应用物理与计算数学研究所，北京 100088）

引 言

膜结构能够在保证结构可靠性的同时耗费更少的材料包围更大的空间，并具有自重更轻、造型优美、便于拆卸等优点，因此在航天工程领域得到了广泛的关注，被期望作为执行深空探测任务的新一代创新空间结构，如薄膜太阳帆推进系统、通讯卫星中大型薄膜天线反射器及太阳能电站中电池阵列结构等[1]。这类空间结构在发射前通常处于折叠状态，入轨后在驱动装置的作用下展开，并进入工作状态。然而对于在太空中展开太阳帆这类具有高收纳比的大型薄膜结构，无疑是一项复杂的操作，且展开后的有效面积及膜面精度均会对其所承担的功能产生一定影响。因此，进入太空中的薄膜航天器能否顺利展开至理想工作状态是这类航天器折展系统设计中的关键技术。考虑到地面实验的成本过大，且难以针对大型薄膜结构提供在轨微重力真空环境[2]。于是，采用高效的数值计算方法对空间帆膜结构的动力学展开过程进行仿真模拟是此类薄膜航天器结构设计的重要手段。Shirasawa等[3]采用多质点模型对“伊卡洛斯”（Icarus）太阳帆展开动力学进行数值模拟，采用由弹簧和阻尼器连接的质点代替传统薄膜单元，近似简化了膜结构建模，降低了计算成本。Miyazaki[4]提出了一种刚度退化的褶皱模型，通过引入刚度折减系数描述薄膜压缩时的本构关系，建立了薄膜太阳帆的柔性多体系统有限元模型，并将其应用于薄膜太阳帆的展开动力学分析。

针对类似薄膜结构的柔性多体系统的大转动、大变形问题，柔性多体系统动力学领域提出了绝对节点坐标方法（Absolute Nodal Coordinate Formulation，ANCF）[5]。基于ANCF方法，国内学者针对薄膜太阳帆等多柔体航天器的展开动力学进行了较为深入的研究[6-9]。赵将等[7]采用ANCF薄板单元，对考虑粘弹性阻尼的“伊卡洛斯”薄膜太阳帆进行了展开动力学分析。Liu等[8]系统地对ANCF板/壳单元进行了研究，结合SRM模型提出了缩减的ANCF薄膜单元，并通过薄膜太阳帆的自旋展开动力学模拟验证了所提出方法的有效性。然而，由于薄膜结构在宏观响应中拉压模量呈现的巨大差异，系统非线性特性显著，导致薄膜结构在静力学/动力学模拟中极易发生计算不稳定现象。

针对拉压不同模量弹性结构，前苏联学者阿姆巴book=218,ebook=122尔楚米扬（Ambartsumyan）[10]系统地提出了不同模量弹性理论，并给出了其本构关系的基本假定，推导了二维平面问题和三维空间问题的广义弹性定律。之后许多研究者基于不同角度对Ambartsumyan的不同模量弹性理论进行了发展和修正[11]。近年来，针对双模量结构的数值算法得到了发展，张亮等[12]基于参数变分原理和有限元方法，在主应力方向上建立了统一的含参数变量的本构方程，将不同模量问题演化为一个易于求解的标准互补二次规划问题，避免了在有限元迭代过程中不断根据应力状态更新刚度矩阵，并将其应用于二维薄膜的起褶分析，呈现出不错的收敛效率。杜宗亮等[13]基于不同模量弹性理论推导了一种具有渐进二次收敛效率的切线刚度矩阵，通过将其应用于二维薄膜的应力水平和起褶区域的预测，验证了该方法在处理不同模量薄膜结构时的高效性、稳定性。

本文旨在基于ANCF建模方法，整合张力场与不同模量弹性理论[10]，提出一种考虑不同模量的ANCF薄膜单元，以实现空间薄膜结构的高效动力学模拟。为提高单元在处理拉压不同模量薄膜问题时的算法稳定性，推导了单元精确的切线刚度矩阵，并通过若干静力学算例验证了单元的高效性。此外，针对一正六边形薄膜太阳帆，建立其展开过程多柔体系统动力学模型，分析了薄膜太阳帆的结构参数对其展开过程的影响，数值仿真结果为这类空间帆膜结构的设计和优化提供了理论指导。

1 基于ANCF描述的薄膜单元

如图1所示，点P(ξ,η)为初始构形中薄膜单元上一点，该点的位置矢量r0(ξ,η)表示为

book=219,ebook=123

2 基于不同模量弹性理论的弹性力及切线刚度

2.1 薄膜张力场理论与不同模量弹性理论

张力场理论[16]认为薄膜的弯曲刚度为零，当薄膜内出现压应力时，一般通过面外变形即薄膜起褶的形式释放压应力，张力场理论将其处理为薄膜的面内收缩。起褶后的薄膜处于单轴拉伸应力状态，褶皱的方向平行于大主应力的方向，且垂直于褶皱方向的小主应力为零。

book=220,ebook=124

2.2 考虑不同模量的ANCF薄膜单元的切线刚度

3 基于广义-α算法的多体系统动力学方程

图2给出了薄膜结构的动力学分析流程图，在预处理模块中，根据薄膜结构的初始构形以及材料参数和单元信息等建立了精确的ANCF薄膜单元，给出了系统的常数质量矩阵。然后，采用高斯消元法消除冗余book=221,ebook=125约束，得到系统的稀疏矩阵。通过假定初始时刻系统的应力场确定初始弹性矩阵，并组装系统的弹性力和雅克比矩阵，然后代入到式（36）～（39）的广义-α算法的求解过程中，根据求解结果更新系统位移场和弹性矩阵，如此循环迭代直至达到给定的收敛准则，求解完成。

图2 膜结构动力学分析计算方法Fig. 2 Dynamic analysis and calculation method of membrane structure

4 数值仿真验证

为便于区分，将式（29）对应的线性近似的切线刚度模型记为ANCF LM，将式（30）对应的切线刚度模型记为ANCF TM。基于本文提出的计算方法，本节主要针对薄膜结构的一些经典案例进行了考察。4.1节、4.2节和4.3节分别验证文所提出的计算框架在求解不同模量小变形、大变形，以及面外变形问题时的高效性和精确性。

4.1 基于ANCF的薄膜单元小变形分析

图3所示为一个验证薄膜起褶模型各种数值算法的经典算例，矩形薄膜宽度H，厚度t，其上下边界施加有大小σ0的均匀预应力，左右两端施加有轴向集中力P=σ0tH，同时左右边界均施加有集中力矩M。随着力矩的增加，膜的下边缘开始出现一定宽度的褶皱区域。考虑采用ANCF DM薄膜单元计算框架，对上述薄膜结构进行了数值模拟，仿真过程中，不考虑薄膜的压缩模量。

图3 纯弯曲矩形薄膜受力示意图Fig. 3 Stress diagram of Pure Bending Rectangular membrane

文献[19]针对上述模型，采用解析方法对其起褶区域以及面内变形进行了研究，并得出一系列关系式。其中，式（40）给出了矩形膜起褶区域的宽度b与外力矩M的解析关系式

式（41）给出了不同工况下，薄膜右边界上一点的柯西应力与外力矩的解析关系式

式（42）给出了薄膜变形后曲率κ与外力矩M之间的解析关系式

图4给出了起褶区域宽度的解析解和数值试验结果，不难看出即使在起褶区域超过90%时，数值模拟结果仍然具有很高的精度。

图4 褶皱带宽与外力矩的关系Fig. 4 Relationship between fold bandwidth and external torque

book=222,ebook=126图5则给出了不同外力矩作用下，右边界上一点柯西应力的解析解与数值试验结果，两者的吻合性比较好。

图5 右边界水平应力与外力矩的关系Fig. 5 Relationship between horizontal stress at right boundary and external moment

图6给出了薄膜结构曲率的解析解和数值解随外力矩的变化趋势，分析得知，两者的趋势是非常一致的。

图6 纯弯曲矩形薄膜的弯矩–曲率关系Fig. 6 Moment-curvature relationship of rectangular membrane

通过上述各项数值试验结果与解析结果的对比可以得出，本文所提出的方法在研究薄膜起褶的问题时具有很高的精度。

4.2 基于ANCF的不同模量结构大变形分析

如图7所示为一L形梁结构，其左端为固定约束，自由端受水平方向集中力F= 40 000 N。固定其受拉方向的弹性模量为Et= 3×107Pa，泊松比为vt= 0.3，依次分别取压缩模量为Ec=3×107Pa、Ec=1.5×107Pa、Ec= 6 ×106Pa、Ec=3×106Pa 以及Ec= 1.5 × 106Pa 时的工况，对应的拉压比分别为Et/Ec=1、2、5、10、20。

为验证本文所提出计算方法的有效性，本节基于ABAQUS平台的扩展程序接口进行了二次开发，并将UMAT的结果作为ANCF描述的不同模量问题的参考解。需要指出的是，UMAT这里采用CPS4I平面应力四边形单元和自适应时间增量步进行计算，在拉压比较大的情况下需要较多的加载步。而本文所提出的ANCF TM模型仅需一个增量步即可得到收敛解。

图7 L形梁结构Fig. 7 L-shaped beam structure

图8给出了不同拉压模量比下，L形梁的变形后构形图，表1给出了采用不同技术途径在不同拉压模量比条件下得到的节点P的水平位移以及计算收敛耗时。以上仿真结果表明，即使在求解拉压模量相差20倍的大变形问题时，本文所提出的计算框架仍能高效收敛，且能够保证较高的精度。

图8 变形后构形图Fig. 8 Configuration diagram after deformation

表1 节点P的水平位移与计算耗时Table 1 Horizontal displacement of nodePand calculation convergence time

book=223,ebook=127此外，如表2所示，在采用ANCF-LM模型求解不同模量大变形问题时，无法得到收敛的结果。值得注意的是，即使拉压模量相同的条件下，ANCF-LM模型仍然无法得到收敛的结果，这是由于在相同模量时，该模型无法退化为经典的线性弹性本构模型。而本文提出的ANCF-TM模型求解过程中牛顿迭代收敛性良好，迭代步数和迭代耗时T随拉压模量比的增大略有增加。

表2 不同模型下节点P的求解效率Table 2 Solution efficiency of nodeP under different models

4.3 安全气囊的充气仿真过程

图9所示为一方形安全气囊的初始构形，薄膜材料的杨氏模量为588.0 Mpa，泊松比为0.4。方膜对角线长度为L= 1 200 mm，厚度t= 1mm。考虑将气囊的充气过程近似为一个准静态过程，气囊在内部递增的均匀压力作用下发生变形，直至内部压强P最终增至5 000 Pa，充气过程结束。仿真过程中，不考虑薄膜的压缩模量。

图9 安全气囊的初始构形示意图Fig. 9 Schematic diagram of initial configuration of airbag

表3给出了基于本文所提出的计算方法在不同网格划分密度下得到的安全气囊上节点O、C、B的位移以及计算耗时T。结果表明，基于本文计算框架得到的位移收敛解与文献[20]的研究结果基本吻合。此外，采用10×10的网格即可得到位移收敛解，相应的计算耗时随着网格的加密有大幅增加。

表3 不同网格密度下的薄膜节点的位移Table 3 Displacement of membrane nodes under different mesh densitied

图10出了安全气囊充气完成后的应力分布示意图，根据应力准则[13]，对单元内每个积分点的应力状态进行标记，用黑点标记的薄膜区域处于褶皱状态，用圆标记的薄膜区域处于松弛状态，其余未标记的区域则处于张紧状态，这一结果与文献[21]的结论基本一致。综上结论再次验证了本文所提出的计算框架在求解薄膜面外变形问题时的可行性。

图10 安全气囊变形后应力分布示意图Fig. 10 Schematic diagram of stress distribution of airbag after inflation

5 薄膜太阳帆的自旋展开动力学模拟

作为一种新型空间可展开结构，能够通过自旋展开并最终稳定的大型帆膜结构在未来的空间探索技术应用中具有巨大潜力。能够准确地对薄膜展开过程的动态响应进行预测和模拟是薄膜太阳帆展开系统设计中的关键技术之一。这里将采用本文所提出的计算框架对薄膜太阳帆的自旋展开过程进行分析，通过确定影响展开的关键因素来指导太阳帆的结构设计。

初始时刻，薄膜太阳帆按照图11所示的折痕折叠后缠绕在与其同心的正六边形彀轮之上，彀轮驱动后，太阳帆随之旋转，并在均布于其六个角上的集中质量块的离心力作用下自旋展开，完全展开后的太阳帆是一个正六边形。此时太阳帆的边长为L= 0.635 m，厚度为t= 0.1 mm。薄膜材料的密度为ρ= 1 420 kg/m3，book=224,ebook=128杨氏模量为E= 2.5×109Pa，泊松比为v= 0.3。彀轮的边长为d= 0.1 m，质量为m0=1.0 kg，相对于过质心的Z轴的惯性矩为Iz= 5.0×10–3kg m2。薄膜边界角点处的集中质量块的质量为m= 0.05 kg，初始时刻系统的角速度为w= 0.628 rad/s。

图11 薄膜太阳帆的有限元模型Fig. 11 Finite element model of membrane solar sail

考虑到结构的对称性，本文仅建立1/6模型，并施加旋转边界条件，该有限元模型中包括35个薄膜片，薄膜片之间通过旋转铰连接，具体建模方法可以参考文献[8]。本文采用考虑拉压不同模量的ANCF薄膜单元和第一类拉格朗日方程对太阳帆进行动力学精确建模，然后利用可控数值耗散的广义-α算法对太阳帆的动力学展开进行求解分析。

5.1 压缩刚度折减系数对展开过程的影响

在对太阳帆进行仿真时考察了3种不同的刚度折减的薄膜模型（Ec=s×Et），其中s为刚度折减系数。模型一为不考虑压缩刚度的不可压ANCF薄膜单元（s=0）；模型二为考虑微弱抗压刚度的ANCF薄膜单元（s=10–5）；模型三同样为引入微弱抗压刚度的ANCF薄膜单元（s= 10–4）。

图12给出了采用模型一进行仿真时薄膜太阳帆的展开过程图，结果表明，此时的薄膜太阳帆无法完全展开并维持展开状态，太阳帆在12 s左右展开至最大后又重新收缩至中心彀轮附近。

图12 太阳帆的展开过程图（s=0）Fig. 12 Deployment process diagram of solar sail （s=0）

图13和图14分别给出了采用模型一进行仿真时太阳帆展开过程中的能量与角动量演化过程。图13表明，即使展开过程中薄膜的应变能出现高频振荡，但是系统的总能量仍然是守恒的。在12 s左右时太阳帆系统的应变能达到了峰值，此时太阳帆处于最大张紧状态。图14表明，相比于中心彀轮和集中质量块，薄膜在整个仿真过程中的角动量仅发生小范围的波动，在仿真时间达到12 s之后，太阳帆系统各部件的角动量均趋于稳定。

图13 太阳帆展开过程中系统的应变能及总能量（s=0）Fig. 13 Strain energy and total energy of the system during solar sail deployment （s= 0）

图14 太阳帆展开过程中系统各部件的角动量（s=0）Fig. 14 Angular momentum of each component of the system during solar sail deployment （s=0）

图15给出了采用模型二进行仿真时薄膜太阳帆的展开过程图，此时的薄膜太阳帆同样无法维持最大展开状态，太阳帆在11 s左右展开到最大半径，之后在边缘集中质量区域出现较大幅度的面外变形。根据图16book=225,ebook=129和图17给出的太阳帆系统展开过程中的能量与和角动量的演化过程，可以得知，系统的应变能在12 s之后稳定在一较小的范围内。此外，在仿真时间达到11 s左右时，太阳帆系统各部件的角动量出现突变，之后逐渐稳定，展开过程中系统的总能量和总角动量均是守恒的。

图15 太阳帆的展开过程图（s=10-5）Fig. 15 Deployment process diagram of solar sail（s=10-5）

图16 太阳帆展开过程中系统的应变能及总能量（s=10-5）Fig. 16 Strain energy and total energy of the system during solar sail deployment（s=10-5）

图17 太阳帆展开过程中系统各部件的角动量（s=10-5）Fig. 17 Angular momentum of each component of the system during solar sail deployment （s=10-5）

图18给出了采用模型三进行仿真时薄膜太阳帆的展开过程图，这时的薄膜太阳帆可以完全展开并维持展开状态。太阳帆在13 s左右展开至最大半径，之后跟随彀轮一起稳定旋转，同时在边缘集中质量区域发生微弱的面外变形。同时根据图19和图20分别给出其展开过程中的能量和角动量的演化过程可以得知，此时薄膜太阳帆系统的总能量和总角动量均是守恒的。

图18 太阳帆的展开过程图(s=10-4)Fig. 18 Deployment process diagram of solar sail(s=10-4)

图19 太阳帆展开过程中系统的应变能及总能量（s=10-4）Fig. 19 Strain energy and total energy of the system during solar sail deployment（s=10-4）

图20 太阳帆展开过程中系统各部件的角动量（s=10-4）Fig. 20 Angular momentum of each component of the system during solar sail deployment（s=10-4）

图21和图22给出了基于不同压缩刚度折减系数的薄膜模型仿真得到的太阳帆系统各部件的动力学响应。图21表明，当刚度折减系数为0和10–5时，太阳帆系统的彀轮和边缘质量块无法进入同步旋转。而在刚度折减系数为10–4的情况下，在仿真时间达到13 s左右时彀轮和边缘质量块的旋转角度开始相等，并进入同步旋转。图22表明，在考虑压缩刚度的情况下，当刚book=226,ebook=130度折减系数较大时，薄膜太阳帆展开至一定半径大小之后趋于稳定。而当刚度折减系数较小或者不考虑压缩刚度时，薄膜太阳帆无法维持最大展开状态。综上分析可知，在对薄膜太阳帆系统进行展开动力学分析时，考虑一定的压缩刚度是必要的，这与文献[4]得出的结果基本是一致的。

图21 不同刚度折减系数时边界质量块和中心彀轮的转角Fig. 21 Rotation angle of edge mass and central hub under different stiffness reduction factors

图22 不同刚度折减系数时的太阳帆展开半径Fig. 22 Deployment radius of solar sail under different stiffness reduction factors

book=227,ebook=131

5.2 集中质量和初始转速对展开过程影响分析

为进一步确定薄膜太阳帆的一些结构参数对其展开动力学特性的影响，本文又分别仿真了当边界集中质量和系统初始转速变化时，太阳帆系统各部件的动力学响应。

图23和图24分别给出了太阳帆系统边界集中质量为0.03、0.05、0.07 kg时其各部件的转角演化过程和太阳帆展开半径演化过程。分析得知，太阳帆分别在16、13、11 s左右展开至最大半径，同时中心彀轮和边界质量块的旋转角度开始趋于同步，展开过程的构形演化图类似图18。此外，综合分析得知，边界集中质量对太阳帆系统的最终展开半径大小以及完全展开后发生的面外变形均会产生一定影响。

图23 不同集中质量时边界质量块和中心彀轮的转角Fig. 23 Rotation angle of edge mass and central hub under different concentrated masses

图24 不同集中质量时的太阳帆展开半径Fig. 24 Deployment radius of solar sail with different concentrated masses

图25和图26给出了系统初始转速分别为0.471、0.628、0.785 rad/s时太阳帆展开过程中各部件的动力学响应。结果表明，太阳帆分别在16.5、13、10 s左右展开至最大半径，同时中心彀轮和边界质量块开始进入同步旋转，其展开过程的构形演化类似图18。综合分析得知，系统初始转速对太阳帆展开至最大半径所耗费的时间以及完全展开后发生的面外变形均会产生显著影响。

图25 不同初始转速时边界质量块和中心彀轮的转角Fig. 25 Rotation angle of edge mass and central hub at different initial speeds

图26 不同初始转速时的太阳帆展开半径Fig. 26 Deployment radius of solar sail at different initial speeds

6 结 论

本文整合ANCF建模方法与张力场理论，针对薄膜结构拉压不同模量的宏观特性，提出了一种高效的建模与计算方法，为空间帆膜结构的展开动力学模拟提供了有效途径。

1）针对薄膜太阳帆系统展开过程中显著的非线性动力学特性，基于张力场理论与不同模量弹性理论，推导了单元精确的切向刚度矩阵，提出了一种考虑不同模量的ANCF薄膜单元计算框架。其中考虑了薄膜单元应力场对其本构关系的影响，准确地描述了薄膜结构在拉压过渡区的非线性特性，有效提高了单元在牛顿迭代过程中的收敛性。

2）基于上述框架，采用第一类拉格朗日方程建立了六边形薄膜太阳帆的自旋展开动力学模型，并通过可控数值耗散的广义-α算法求解运动方程。根据数值模拟结果，分析了太阳帆不同结构设计参数对其展开book=228,ebook=132动力学特性的影响，为太阳帆类空间帆膜结构的系统设计和折展动力学仿真提供了一定的理论指导。

3）未来值得进一步研究的问题和方向包括：薄膜太阳帆完全展开之后，帆面产生的微弱面外变形以及振荡过程仍未得到有效解决；现有的针对薄膜太阳帆自旋展开的动力学建模与仿真方法在计算效率和分析精度方面仍然存在不足，随着薄膜单元网格的加密，难以针对大型帆膜结构进行高效模拟。