CN104484498B - 确定有效推进加速度折损的太阳帆瞬态动力学分析方法 - Google Patents

Abstract

确定有效推进加速度折损的太阳帆瞬态动力学分析方法，(1)建立太阳帆有限元模型；(2)进入第t个求解时刻，t从t₀+Δt开始；(3)根据每个帆面单元与光压方向时刻的夹角，计算第t时刻、含预应力下的每个帆面单元的刚度矩阵与质量矩阵，进而得到每个帆面单元的受光压面载荷；(4)得到含所有p个单元的总体刚度矩阵K_t与载荷列向量F_t，以及整体质量m和总体质量矩阵M_t；(5)根据第t‑Δt时刻的末位移、末速度、末加速度，即第t个时刻的初始条件，以及第t个时刻的总体刚度矩阵与载荷列向量以及总体质量矩阵，计算第t时刻的末位移、末速度、末加速；(6)开始进行第t+Δt时刻，重复(3)～(5)，直至最终求解时间t_n；(7)根据上述结果得到加速度的折损。

Description

确定有效推进加速度折损的太阳帆瞬态动力学分析方法

技术领域

本发明涉及一种确定有效推进加速度折损的基于精确时变光压力模型的太阳帆瞬态动力学分析方法，尤其适用于含褶皱具有无法度量面型精度以及所处不同轨道与姿态的各种的大型超柔度太阳帆结构的精确有效加速度分析确定。

背景技术

太阳帆是一种先进的无工质长期在轨运行航天器，由于它仅利用太阳在大面积薄膜上的反射光压提供飞行动力，凡是太阳光存在的地方，太阳帆航天器便可以源源不断的获得动力，最终实现普通航天器无法达到的速度。因为太阳帆不受燃料的约束，可以在宇宙中长时间运行，因此特别适用于深空探测任务。

太阳帆通过控制其与太阳光线的夹角改变光压的大小，可以使光压成为动力或者阻力，从而控制太阳帆运行的飞行速度和角度，以及轨道位置和姿态位置。尽管太阳帆产生的有效推进加速度仅仅约为0.001m/s^-2，但因为涉及在轨姿态调节与轨道保持，对其加速度的要求较为苛刻，需要各方向较为精确的平动与转动加速度值。对与太阳帆而言，毫无疑问，其法向推进加速度至关重要，该数值即为太阳帆最重要的评价指标。

但是，国内外太阳帆研究人员通常在计算瞬态有效推进加速度时，都未考虑光压的时变效应。如图2和图3所示，传统太阳帆帆面初始状态的非理想平面以及随着帆面在运行过程中的变形两大因素，非时变光压模型的入射载荷总是与帆面微元垂直(细实线)，而实际的光压模型的入射光压方向始终不变(粗实线)，这就会在最终计算有效推进光压力时漏乘了一次微元与光压方向夹角的方向余弦，导致最终计算得到的有效推进载荷(细虚线)比真实时变情况(粗虚线)略有增加，并且可以注意到，帆面变形越大，对非时变光压模型的计算误差越大。即最终计算得到的有效推进加速度比实际时变状态略有增加。

在真实受力情况下，太阳帆每一帆面微元所受的光压力初始载荷经过两次投影折损相乘计算即为太阳帆真实的有效推进力方向，这两次投影折损计算分别为：帆面入射载荷与帆面法向之间的投影计算，即为第一次折损；每一帆面微元载荷与太阳帆帆面整体法向之间的投影计算，即为第二次折损。正是考虑到这种不真实的非时变光压模型的不准确性带来的加速度非真实增加，研究建立更符合真实状态的时变光压模型，计算时变光压加速度的折损。

因此，利用一种确定有效推进加速度折损的基于精确时变光压力模型的太阳帆瞬态动力学分析方法成为该领域目前亟待解决的技术问题。

发明内容

本发明的技术解决问题是：克服现有分析方法的不足，提供一种确定有效推进加速度折损的基于精确时变光压力模型的太阳帆瞬态动力学分析方法。

本发明的技术解决方案是：确定有效推进加速度折损的太阳帆瞬态动力学分析方法，步骤如下：

(1)建立太阳帆有限元模型，即建立支撑臂一维梁单元与帆面二维三角形壳单元模型，设置材料参数与截面属性，施加预紧力，建立支撑臂与帆面连接节点约束；设置光压力载荷P，设置动力学分析的初始时间t₀，最终时间t_n与时间步长Δt以及所有节点的初始速度及加速度列向量v₀与a₀；对帆面二维三角形壳单元模型中的帆面单元进行统一单元编号从1,2,…,p共p个单元；对太阳帆有限元模型中节点统一编号从1,2,…,q共q个节点；

(2)进入第t个求解时刻，t从t₀+Δt开始；

(3)根据每个帆面单元与光压方向时刻的夹角，计算第t时刻、含预应力下的的每个帆面单元的刚度矩阵与质量矩阵，进而得到每个帆面单元的受光压面载荷；

(4)根据步骤(3)的结果，得到含所有p个单元的总体刚度矩阵K_t与载荷列向量F_t，以及整体质量m和总体质量矩阵M_t；

(5)根据第t-Δt时刻的末位移、末速度、末加速度，即第t个时刻的初始条件，以及步骤(4)第t个时刻的总体刚度矩阵与载荷列向量以及总体质量矩阵，计算第t时刻的末位移、末速度、末加速；

(6)开始进行第t+Δt时刻，重复步骤(3)～(5)，直至最终求解时间tn；

(7)将(6)计算得到的时间历程下，太阳帆中心点所有时刻的末加速度，与现有的非时变光压力模型得到的瞬时有效推进加速度数值进行比较，并根据最后时刻末位移求得平均有效推进加速度，进而得到加速度的折损。

步骤(7)中可以提取求解历程中太阳帆特征点位移与应力数值，分析其瞬态动力学特性，确定最大变形。

所述步骤(3)具体实现步骤如下：

(3.1)从步骤(1)中建立的太阳帆有限元模型中读取第i个帆面单元的节点初始坐标Nⁱ＝(xⁱ,yⁱ,zⁱ,αⁱ,βⁱ,γⁱ,)；

(3.2)根据第i个帆面单元的在第t-Δt时刻的位移以及(3.1)的节点初始坐标Nⁱ＝(xⁱ,yⁱ,zⁱ,αⁱ,βⁱ,γⁱ,)，计算第t时刻的初始位置第i个帆面单元的面积总体坐标系下的方向向量含预应力下的第t时刻的帆面单元刚度矩阵与质量矩阵求解初始时间t₀时刻的位移为0；

(3.3)根据步骤(3.2)的结果计算第i个帆面单元的受光压面载荷并改写载荷节点力列向量

(3.4)开始读取第i+1个帆面单元，重复步骤(3.1)～(3.3)，直至全部帆面单元；

(3.5)组成含所有p个单元的第t时刻的总体刚度矩阵K_t和总体质量矩阵M_t，并组成载荷列向量

本发明与现有技术相比有益效果为：

(1)时变光压力模型更加准确的反映了太阳帆在轨的真实状态，并且完全覆盖了太阳帆位于不同轨道和所处不同姿态及不同光压载荷的各种在轨运行状态。应用该方法只需知道太阳帆某一帆面单元与光压方向时刻的夹角即可计算任意状态的太阳帆实际所受光压载荷，并且这种时变光压载荷是与真实情况完全符合的，即把每一微元当做理想帆面，而不是将太阳帆整体结构当做理想帆面。

(2)考虑太阳帆受光压载荷下薄膜的法向变形，使用传统时不变光压力模型无法体现该状态的两次光压力折减，时变光压模型可以有效的计算每时刻的有效光压力的变化，严格计算两步方向余弦，从而得到真实的平均加速度的折减量。

(3)考虑太阳帆薄膜的褶皱情况以及所受光压载荷下薄膜的法向变形两者同时叠加的时变效应，将会导致极其不平整且凹陷的理想帆面，时变光压模型可以有效的计算各种初始条件和过程中帆面单元形函数，构建每一时刻的总体刚度矩阵与总体质量矩阵。

(4)对光压力载荷的施加不需如传统的计算方法一样，必须已知确定的光压力，而是只需知道当前位置太阳光压系数，即可通过每个帆面单元节点的具体坐标及其方向向量和单元面积求得真实施加在该单元的太阳光压力。

(5)相较于传统非时变光压力模型，载荷因为漏乘一次方向余弦导致施加的载荷过大，造成各分析结果的不准确，以致影响设计人员对太阳帆的前期设计工作过于保守。而通过时变光压模型，真实确定载荷，通过瞬态动力学计算便可以确定求解历程中太阳帆特征点位移与应力数值，分析其瞬态动力学特性，确定真实的最大变形。

(6)鉴于仅从常规的平均加速度数值无法有效描述该方法的太阳帆推进加速度真实折损的力学特性，本发明特别针对用户关注的该加速度数值，绘制太阳帆中心点的加速度历程时间曲线，并与传统非时变光压和理论加速度值进行对比，方便的给出设计人员关心的真实太阳帆推进加速度折损量，便于方案对比。

附图说明

图1为本发明流程图；

图2为本发明与非时变光压模型的比较；

图3为太阳帆微元载荷示意图与等效节点力；

图4为太阳帆有限元图；

图5为本发明(精确时变光压模型)与传统非时变光压模型太阳帆中心位置位移-时间历程曲线；

图6为精确时变光压模型与传统非时变光压模型太阳帆位移折损-时间历程曲线；

图7为时变、非时变模型平均加速度与理论加速度的折损-时间历程曲线；

图8为时变相对于非时变模型平均加速度折损-时间历程曲线；

具体实施方式

下面结合附图及实施例对本发明做详细说明，本发明一种确定有效推进加速度折损的基于精确时变光压力模型的太阳帆瞬态动力学分析方法，如图1所示，步骤如下：

(1)建立太阳帆有限元模型，主要包括：分别建立支撑臂一维梁单元与帆面二维三角形壳单元模型，设置材料参数与截面属性，施加预紧力，建立支撑臂与帆面连接节点约束；设置光压力载荷P，设置求解初始时间t₀，最终时间t_n与时间步长Δt，以及初始条件每个节点的速度、加速度列向量v₀与a₀；假设在整个分析过程中单元质量建立统计单元编号，从1,2,…,p共p个单元，统计节点编号，从1,2,…,q共q个节点。

(2)进入第t个求解时刻，t从t₀+Δt开始；

(3)读取第i个帆面单元的节点初始坐标Nⁱ＝(xⁱ,yⁱ,zⁱ,αⁱ,βⁱ,γⁱ,)；

单元坐标系如图3所示的o_ex_ey_ez_e坐标系，单元坐标原点位于该单元最小节点号位置，单元坐标系x轴为从该节点指向该单元最大节点位置，单元坐标系z轴为该单元面外法向并与总体坐标系z轴夹角成小于90度一侧，单元坐标系y轴成右手法则。Nⁱ＝(xⁱ,yⁱ,zⁱ,αⁱ,βⁱ,γⁱ,)中的六个量分别为三个平动自由度与三个转动自由度。

(4)将上一求解时刻的第i个帆面单元的在第t-Δt时刻的位移作为下一步的初始条件，求解初始时间t₀时刻的位移为0；

每一求解时刻的计算得到的瞬时加速度用来绘制加速度-时间历程曲线，最后时刻计算得到的末位移用来求解平均有效推进加速度。

(5)根据步骤(3)的节点初始坐标Nⁱ＝(xⁱ,yⁱ,zⁱ,αⁱ,βⁱ,γⁱ,)和步骤(4)第t-Δt时刻的位移计算第t时刻的初始位置第i个单元的面积总体坐标系下的方向向量含预应力下的第t时刻的单元刚度矩阵与质量矩阵

面积和方向向量都属于公知常识，单元刚度矩阵与质量矩阵参见，王勖成.有限单元法.清华大学出版社,2003。以上四个量都是根据上一步的计算结果再次更新计算得到的。

光压力与单元坐标系下的三个轴都不平行，但与总体坐标系下的z轴平行。在单元坐标系中构建单元刚度(质量)矩阵，将单元阵扩阶、组装与坐标转换后变为总体坐标系下的总体刚度(质量)阵。单元坐标系为图3中的o_ex_ey_ez_e，总体坐标系为oxyz，总体坐标系原点为太阳帆初始位置中心点，延两支撑臂方向分别为x、y轴，光压方向为总体坐标系中的z轴正向，三方向成右手法则。

(6)根据步骤(5)计算第i个单元的受光压面载荷，并按单元坐标系和总体坐标系的坐标转换公式(坐标转换公式为公知)计算总体坐标系下的单元节点力并改写载荷节点力列向量

其中，为组成第i个单元的第1个节点在t时刻的总体坐标系x方向的单元节点力，为组成第i个单元的第3个节点在t时刻的总体坐标系z方向的单元节点力矩，其他参数含义依次类推。

(7)开始读取第i+1个单元，重复步骤(3)～(6)，直至全部单元；

(8)用所有p个帆单元的上述计算结果组成第t时刻的总体刚度矩阵K_t和总体质量矩阵M_t，并组成载荷列向量

第t时刻的等效总体刚度矩阵其中为等效总体刚度矩阵，α₀、α₁为积分常数，C_t为第t时刻比例阻尼矩阵，即C_t＝αK_t+βM_t，α与β分别为Rayleigh阻尼常系数。

总体质量矩阵M_t与总体刚度矩阵K_t为按单元质量(刚度)阵的自由度组装和扩阶的，属于有限元的公知。

(9)根据第t-Δt个时刻的末位移、末速度、末加速度，即第t个时刻的初始条件，以及步骤(8)第t个时刻的总体刚度矩阵K_t与载荷列向量F_t，以及总体质量矩阵M_t，利用Newmark动力学求解方法计算第t个时刻的末位移、末速度、末加速；

具体步骤为

根据t-Δt时刻的位移δ_t-Δt、速度和加速度计算第t时刻有效载荷

求解第t时刻的位移，

计算t时刻的加速度与速度

上述，α₀-α₇是Newmark方法的积分常数，属于公知。参见，王勖成.有限单元法.清华大学出版社,2003。

(10)开始进行第t+Δt时刻，重复步骤(3)～(9)，直至最终求解时间tn。

(11)将(10)计算得到的时间历程下太阳帆中心点所有时刻瞬时有效推进加速度数值(即上述计算的所有时刻的末加速度)，与现有的非时变光压力模型得到的瞬时有效推进加速度数值进行比较，分别绘制两种光压力模型的有效推进加速度时间历程曲线，利用最后时刻的末位移求得平均有效推进加速度进而得到加速度折损。

非时变光压力模型没有考虑帆面每一时刻的变形，导致光压力输入载荷的人为增加，与真实情况不符。参见，McInnes C R.Solar Sailing Technology，Dynamics andMission Applications.Springer-Verlag,Berlin,1999.洪延姬等.先进航天推进技术.国防工业出版社,2012.

(12)提取求解历程中太阳帆特征点位移与应力数值，分析其瞬态动力学特性，确定最大变形。

实施例：160m太阳帆结构静力学有限元仿真分析

1)建立太阳帆有限元模型，如4所示，主要包括：分别建立支撑臂一维梁单元与帆面二维三角形壳单元模型，设置方形帆边长为160m，设置中心体质量m_center＝100kg，端点配重m_vertex＝12.4kg、中心体转动惯量J_center＝0kg·m²，支撑臂外径D＝300mm、壁厚b＝0.1mm，帆面厚度h＝1.5μm；设置支撑臂材料属性E_beam＝210GPa、γ_beam＝0.33、ρ_beam＝848kg/m³，帆面薄膜材料属性E_sail＝2.5GPa、γ_sail＝0.34、ρ_sail＝1420kg/m³；设置预应力方向为帆面内的两个正应力，大小为50000Pa；设置求解初始时间t₀＝0s，最终时间t_n＝400s与时间步长Δt＝0.5s，以及初始条件各节点的速度、加速度列向量v₀＝0m/s与a₀＝0m/s²，节点总数q＝3069、单元总数p＝5633；设置地球附近全反射太阳帆光压载荷数值P＝9.12×10^-6Pa；最终确定l＝112m、面积S＝25088m²、总质量为m＝205.6045kg、总转动惯量为I_x＝241961kg·m²、I_y＝241961kg·m²、I_z＝488602kg·m²。

2)进入第t个求解时刻；

3)读取第i个帆面单元的节点初始坐标Nⁱ＝(xⁱ,yⁱ,zⁱ,αⁱ,βⁱ,γⁱ,)；

4)读取第i个帆面单元的在第t-Δt时刻的位移

5)根据步骤(3)的节点初始坐标Nⁱ＝(xⁱ,yⁱ,zⁱ,αⁱ,βⁱ,γⁱ,)和步骤(4)第t-Δt时刻的位移计算第t时刻的初始位置第i个单元的面积总体坐标系下的方向向量含预应力下的第t时刻的单元刚度矩阵与质量矩阵

6)根据步骤(5)计算第i个单元的受光压面载荷并改写成载荷节点力列向量，如图3所示；

7)开始读取第i+1个单元，重复步骤(3)～(6)，直至全部单元；

8)用所有p＝5633个单元刚度矩阵组成第t时刻的总体刚度矩阵K_t和总体质量矩阵M_t，并组成载荷列向量

9)根据第t-Δt个时刻的末位移、末速度、末加速度，即第t个时刻的初始条件，以及步骤(8)第t个时刻的总体刚度矩阵K_t与载荷列向量F_t，以及总体质量矩阵M_t，利用Newmark动力学求解方法计算第t个时刻的末位移、末速度、末加速；

具体步骤为

形成第t时刻的等效总体刚度矩阵Ct为第t时刻比例阻尼矩阵，即C_t＝αK_t+βM_t，α与β分别为Rayleigh阻尼常系数。

根据t-Δt时刻的位移δ_t-Δt、速度和加速度计算第t时刻有效载荷

求解第t时刻的位移，

计算t时刻的加速度与速度

其中α₀～α₆为与时间步长Δt相关的积分常数为通用的Newmark方法默认值。

10)开始进行第t+Δt时刻，重复步骤(3)～(9)，直至最终求解时间t_n＝400s。

11)将计算求解得到的时间历程前400s太阳帆中心点所有时刻瞬时有效推进加速度数值与非时变光压力模型的瞬时有效推进加速度数值绘制图5～图8时间历程曲线图。将计算求解得到的400s末太阳帆中心点位移计算平均有效推进加速度数值a_TV＝0.00113144875m/s²，并与非时变光压力模型的平均有效推进加速度a_TIV＝0.0011315075m/s²以及理论加速度计算折损，分别得到时变光压模型(本发明)相对非时变光压模型的折损以及时变光压模型相对理论加速度的折损

12)提取求解历程中太阳帆特征点位移与应力数值，分析其瞬态动力学特性，确定最大变形，结果为帆面边缘最大变形1.03m，帆面中心最大0.65m，帆面直边中点变形0.44m，支撑臂端点最大变形0.15m，太阳帆中心薄膜处最大Von Mises应力为0.6MPa，支撑臂端点处薄膜最大Von Mises应力为1.5MPa，帆面中心处薄膜处薄膜最大Von Mises应力为0.02MPa。

本发明未详细说明部分属于本领域技术人员公知常识。

Claims (2)

1.确定有效推进加速度折损的太阳帆瞬态动力学分析方法，其特征在于步骤如下：

(1)建立太阳帆有限元模型，即建立支撑臂一维梁单元与帆面二维三角形壳单元模型，设置材料参数与截面属性，施加预紧力，建立支撑臂与帆面连接节点约束；设置光压力载荷P，设置动力学分析的初始时间t₀，最终时间t_n与时间步长Δt以及所有节点的初始速度及加速度列向量v₀与a₀；对帆面二维三角形壳单元模型中的帆面单元进行统一单元编号，从1,2,…,p共p个单元；对太阳帆有限元模型中节点统一编号从1,2,…,q共q个节点；

(2)进入第t个求解时刻，t从t₀+Δt开始；

(3)根据每个帆面单元与光压方向时刻的夹角，计算第t时刻、含预应力下的的每个帆面单元的刚度矩阵与质量矩阵，进而得到每个帆面单元的受光压面载荷；具体实现步骤如下：

(3.1)从步骤(1)中建立的太阳帆有限元模型中读取第i个帆面单元的节点初始坐标Nⁱ＝(xⁱ,yⁱ,zⁱ,αⁱ,βⁱ,γⁱ,)；

(3.2)根据第i个帆面单元的在第t-Δt时刻的位移以及(3.1)的节点初始坐标Nⁱ＝(xⁱ,yⁱ,zⁱ,αⁱ,βⁱ,γⁱ,)，计算第t时刻的初始位置第i个帆面单元的面积总体坐标系下的方向向量含预应力下的第t时刻的帆面单元刚度矩阵与质量矩阵求解初始时间t₀时刻的位移为0；

(3.3)根据步骤(3.2)的结果计算第i个帆面单元的受光压面载荷并改写载荷节点力列向量

其中，为组成第i个单元的第1个节点在t时刻的总体坐标系x方向的单元节点力，为组成第i个单元的第3个节点在t时刻的总体坐标系z方向的单元节点力矩，其他参数含义依次类推；P是太阳帆光压载荷，单位Pa；

(3.4)开始读取第i+1个帆面单元，重复步骤(3.1)～(3.3)，直至全部帆面单元；

(3.5)组成含所有p个单元的第t时刻的总体刚度矩阵K_t和总体质量矩阵M_t，并组成载荷列向量

(4)根据步骤(3)的结果，得到含所有p个单元的总体刚度矩阵K_t与载荷列向量F_t，以及整体质量m和总体质量矩阵M_t；

(5)根据第t-Δt时刻的末位移、末速度、末加速度，即第t个时刻的初始条件，以及步骤(4)第t个时刻的总体刚度矩阵与载荷列向量以及总体质量矩阵，计算第t时刻的末位移、末速度、末加速；

(6)开始进行第t+Δt时刻，重复步骤(3)～(5)，直至最终求解时间t_n；

(7)将(6)计算得到的时间历程下，太阳帆中心点所有时刻的末加速度，与现有的非时变光压力模型得到的瞬时有效推进加速度数值进行比较，并根据最后时刻末位移求得平均有效推进加速度，进而得到加速度的折损。

2.根据权利要求1所述的确定有效推进加速度折损的太阳帆瞬态动力学分析方法，其特征在于：步骤(7)中可以提取求解历程中太阳帆特征点位移与应力数值，分析其瞬态动力学特性，确定最大变形。