-
膜结构能够在保证结构可靠性的同时耗费更少的材料包围更大的空间，并具有自重更轻、造型优美、便于拆卸等优点，因此在航天工程领域得到了广泛的关注，被期望作为执行深空探测任务的新一代创新空间结构，如薄膜太阳帆推进系统、通讯卫星中大型薄膜天线反射器及太阳能电站中电池阵列结构等[1]。这类空间结构在发射前通常处于折叠状态，入轨后在驱动装置的作用下展开，并进入工作状态。然而对于在太空中展开太阳帆这类具有高收纳比的大型薄膜结构，无疑是一项复杂的操作，且展开后的有效面积及膜面精度均会对其所承担的功能产生一定影响。因此，进入太空中的薄膜航天器能否顺利展开至理想工作状态是这类航天器折展系统设计中的关键技术。考虑到地面实验的成本过大，且难以针对大型薄膜结构提供在轨微重力真空环境[2]。于是，采用高效的数值计算方法对空间帆膜结构的动力学展开过程进行仿真模拟是此类薄膜航天器结构设计的重要手段。Shirasawa等[3]采用多质点模型对“伊卡洛斯”（Icarus）太阳帆展开动力学进行数值模拟，采用由弹簧和阻尼器连接的质点代替传统薄膜单元，近似简化了膜结构建模，降低了计算成本。Miyazaki[4]提出了一种刚度退化的褶皱模型，通过引入刚度折减系数描述薄膜压缩时的本构关系，建立了薄膜太阳帆的柔性多体系统有限元模型，并将其应用于薄膜太阳帆的展开动力学分析。
针对类似薄膜结构的柔性多体系统的大转动、大变形问题，柔性多体系统动力学领域提出了绝对节点坐标方法（Absolute Nodal Coordinate Formulation，ANCF）[5]。基于ANCF方法，国内学者针对薄膜太阳帆等多柔体航天器的展开动力学进行了较为深入的研究[6-9]。赵将等[7]采用ANCF薄板单元，对考虑粘弹性阻尼的“伊卡洛斯”薄膜太阳帆进行了展开动力学分析。Liu等[8]系统地对ANCF板/壳单元进行了研究，结合SRM模型提出了缩减的ANCF薄膜单元，并通过薄膜太阳帆的自旋展开动力学模拟验证了所提出方法的有效性。然而，由于薄膜结构在宏观响应中拉压模量呈现的巨大差异，系统非线性特性显著，导致薄膜结构在静力学/动力学模拟中极易发生计算不稳定现象。
针对拉压不同模量弹性结构，前苏联学者阿姆巴尔楚米扬（Ambartsumyan）[10]系统地提出了不同模量弹性理论，并给出了其本构关系的基本假定，推导了二维平面问题和三维空间问题的广义弹性定律。之后许多研究者基于不同角度对Ambartsumyan的不同模量弹性理论进行了发展和修正[11]。近年来，针对双模量结构的数值算法得到了发展，张亮等[12]基于参数变分原理和有限元方法，在主应力方向上建立了统一的含参数变量的本构方程，将不同模量问题演化为一个易于求解的标准互补二次规划问题，避免了在有限元迭代过程中不断根据应力状态更新刚度矩阵，并将其应用于二维薄膜的起褶分析，呈现出不错的收敛效率。杜宗亮等[13]基于不同模量弹性理论推导了一种具有渐进二次收敛效率的切线刚度矩阵，通过将其应用于二维薄膜的应力水平和起褶区域的预测，验证了该方法在处理不同模量薄膜结构时的高效性、稳定性。
本文旨在基于ANCF建模方法，整合张力场与不同模量弹性理论[10]，提出一种考虑不同模量的ANCF薄膜单元，以实现空间薄膜结构的高效动力学模拟。为提高单元在处理拉压不同模量薄膜问题时的算法稳定性，推导了单元精确的切线刚度矩阵，并通过若干静力学算例验证了单元的高效性。此外，针对一正六边形薄膜太阳帆，建立其展开过程多柔体系统动力学模型，分析了薄膜太阳帆的结构参数对其展开过程的影响，数值仿真结果为这类空间帆膜结构的设计和优化提供了理论指导。
-
如图1所示，点P(ξ,η)为初始构形中薄膜单元上一点，该点的位置矢量r0(ξ,η)表示为
图 1 基于ANCF描述的薄膜单元
Figure 1. Membrane element based on ANCF description
$$ {{\boldsymbol{r}}_0}\left( {\xi ,\eta } \right) = {\boldsymbol{S}}\left( {\xi ,\eta } \right){{\boldsymbol{e}}_0} $$ (1) 其中：0 ≤ ξ ≤ 1和0 ≤ η ≤ 1为点P所在曲面的正则参数；e0为膜单元初始广义坐标；S为膜单元的形函数矩阵，具体形式见文献[14]。
对于薄膜单元，一般认为其截面厚度保持不变。与全参数板单元[15]有所不同，薄膜单元取节点处的位置向量及其在长度和宽度两个方向上的梯度向量作为广义坐标
$$ \left\{ \begin{gathered} {{\boldsymbol{e}}_A} = {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{r}}_A^{\text{T}}}&{{\boldsymbol{r}}_{A,x}^{\text{T}}}&{{\boldsymbol{r}}_{A,y}^{\text{T}}} \end{array}} \right]^{\text{T}}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\boldsymbol{e}}_B} = {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{r}}_B^{\text{T}}}&{{\boldsymbol{r}}_{B,x}^{\text{T}}}&{{\boldsymbol{r}}_{B,y}^{\text{T}}} \end{array}} \right]^{\text{T}}} \hfill \\ {{\boldsymbol{e}}_C} = {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{r}}_C^{\text{T}}}&{{\boldsymbol{r}}_{C,x}^{\text{T}}}&{{\boldsymbol{r}}_{C,y}^{\text{T}}} \end{array}} \right]^{\text{T}}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {{\boldsymbol{e}}_D} = {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{r}}_D^{\text{T}}}&{{\boldsymbol{r}}_{D,x}^{\text{T}}}&{{\boldsymbol{r}}_{D,y}^{\text{T}}} \end{array}} \right]^{\text{T}}} \hfill \\ {{\boldsymbol{e}}_0} = {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{e}}_A^{\text{T}}}&{{\boldsymbol{e}}_B^{\text{T}}}&{{\boldsymbol{e}}_C^{\text{T}}}&{{\boldsymbol{e}}_D^{\text{T}}} \end{array}} \right]^{\text{T}}} \hfill \\ \end{gathered} \right. $$ (2) 为了对薄膜单元的变形进行精确描述，考虑在初始构形及当前构形中建立局部的曲线坐标系和对应的笛卡尔坐标系。其中，初始构形下的两个坐标系定义为
$$ \left\{ \begin{gathered} {\left( {{{\boldsymbol{g}}_0}} \right)_1} = \frac{{\partial {{\boldsymbol{r}}_0}}}{{\partial \xi }},{\kern 1pt} {\left( {{{\boldsymbol{g}}_0}} \right)_2} = \frac{{\partial {{\boldsymbol{r}}_0}}}{{\partial \eta }},{\kern 1pt} {\left( {{{\boldsymbol{g}}_0}} \right)_3} = \frac{{{{\left( {{{\boldsymbol{g}}_0}} \right)}_1} \times {{\left( {{{\boldsymbol{g}}_0}} \right)}_2}}}{{\left| {{{\left( {{{\boldsymbol{g}}_0}} \right)}_1} \times {{\left( {{{\boldsymbol{g}}_0}} \right)}_2}} \right|}} \hfill \\ {\left( {{{\boldsymbol{i}}_0}} \right)_1} = \frac{{{{\left( {{{\boldsymbol{g}}_0}} \right)}_1}}}{{\left| {{{\left( {{{\boldsymbol{g}}_0}} \right)}_1}} \right|}},{\left( {{{\boldsymbol{i}}_0}} \right)_2} = {\left( {{{\boldsymbol{i}}_0}} \right)_3} \times {\left( {{{\boldsymbol{i}}_0}} \right)_1},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\left( {{{\boldsymbol{i}}_0}} \right)_3} = {\left( {{{\boldsymbol{g}}_0}} \right)_3} \hfill \\ \end{gathered} \right.{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} $$ (3) 在膜单元的初始构形和当前构形中分别取微元dX和微元dx，根据连续介质力学知识，两段微元通过位移梯度张量G建立联系
$$ {\text{d}}{\boldsymbol{x}} = {\left[ {\begin{array}{*{20}{c}} {{\text{d}}x}&{{\text{d}}y} \end{array}} \right]^{\text{T}}} = {\boldsymbol{G}} \cdot {\text{d}}{\boldsymbol{X}} = {\boldsymbol{G}} \cdot {\left[ {\begin{array}{*{20}{c}} {{\text{d}}X}&{{\text{d}}Y} \end{array}} \right]^{\text{T}}} $$ (4) 初始构形中，微元dX在局部的曲线坐标系和笛卡尔坐标系中分量分别为[dξ, dη]T和[dX, dY]T，两者之间的转换关系为
$$ \begin{aligned}[b] {\text{d}}{\boldsymbol{\xi }} = & \left[ {\begin{array}{*{20}{l}} {{\text{d}}\xi } \\ {{\text{d}}\eta } \end{array}} \right] = {\left[ {\begin{array}{*{20}{c}} {{{\left( {{{\boldsymbol{g}}_0}} \right)}_1} \cdot {{\left( {{{\boldsymbol{i}}_0}} \right)}_1}}&{{{\left( {{{\boldsymbol{g}}_0}} \right)}_1} \cdot {{\left( {{{\boldsymbol{i}}_0}} \right)}_2}} \\ {{{\left( {{{\boldsymbol{g}}_0}} \right)}_2} \cdot {{\left( {{{\boldsymbol{i}}_0}} \right)}_1}}&{{{\left( {{{\boldsymbol{g}}_0}} \right)}_2} \cdot {{\left( {{{\boldsymbol{i}}_0}} \right)}_2}} \end{array}} \right]^{ - {\text{T}}}}\\ & {\left[ {\begin{array}{*{20}{c}} {{\text{d}}X} \\ {{\text{d}}Y} \end{array}} \right] = {{\boldsymbol{T}}_0}{\text{d}}{\boldsymbol{X}}} \end{aligned}$$ (5) 同理，在当前构形中有
$$ {\text{d}}{\boldsymbol{\xi }} = \left[ {\begin{array}{*{20}{c}} {{\text{d}}\xi } \\ {{\text{d}}\eta } \end{array}} \right] = {\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{g}}_1} \cdot {{\boldsymbol{i}}_1}}&{{{\boldsymbol{g}}_1} \cdot {{\boldsymbol{i}}_2}} \\ {{{\boldsymbol{g}}_2} \cdot {{\boldsymbol{i}}_1}}&{{{\boldsymbol{g}}_2} \cdot {{\boldsymbol{i}}_2}} \end{array}} \right]^{ - {\text{T}}}}\left[ {\begin{array}{*{20}{c}} {{\text{d}}x} \\ {{\text{d}}y} \end{array}} \right] = {\boldsymbol{T}}{\text{d}}{\boldsymbol{x}} $$ (6) 联立式（5）和式（6），得到局部坐标系下的位移梯度张量G的表达式G=T–1T0，由介质力学知识可知，Green-Langrange应变张量的表达式为
$$ {\boldsymbol{E}} = \frac{1}{2}\left( {{{\boldsymbol{G}}^{\text{T}}}{\boldsymbol{G}} - {\boldsymbol{G}}_0^{\text{T}}{{\boldsymbol{G}}_0}} \right) = \frac{1}{2}{\boldsymbol{T}}_0^{\text{T}}\left[ {\begin{array}{*{20}{c}} {{g_{11}} - {{\left( {{g_0}} \right)}_{11}}}&{{g_{12}} - {{\left( {{g_0}} \right)}_{12}}} \\ {{g_{12}} - {{\left( {{g_0}} \right)}_{12}}}&{{g_{22}} - {{\left( {{g_0}} \right)}_{22}}} \end{array}} \right]{{\boldsymbol{T}}_0} $$ (7) 其中：
${g_{\alpha \beta }} = {{\boldsymbol{g}}_\alpha } \cdot {{\boldsymbol{g}}_\beta }{\kern 1pt} {\kern 1pt} {\kern 1pt} \left( {\alpha = 1,2;{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \beta = 1,2} \right)$ 。根据单元的应变张量，给出薄膜单元内力所做的虚功为$$ {\delta}U = \int_{{V^e}} {{\delta}{{\boldsymbol{\varepsilon }}^{\text{T}}}{\boldsymbol{\sigma }}{\text{d}}{V^e}} = h\int_A {{\delta}{{\boldsymbol{\varepsilon }}^{\text{T}}}{\boldsymbol{\sigma }}J{\text{d}}A} $$ (8) 其中：A为母单元面积；h为膜单元厚度；σ为应力张量对应的Vogit向量；ε为应变张量E对应的Vogit向量；
$J = \left| {{{\left( {{{\boldsymbol{g}}_0}} \right)}_1} \times {{\left( {{{\boldsymbol{g}}_0}} \right)}_2}} \right|$ 。$$ {{\delta}}U = h\int_A {{{\delta }}{{\boldsymbol{\varepsilon }}^{\text{T}}}{\boldsymbol{\sigma }}J{\text{d}}A} = {{\delta }}{{\boldsymbol{e}}^{\text{T}}}\left( {h\int_A {{{\left( {\frac{{\partial {\boldsymbol{\varepsilon }}}}{{\partial {\boldsymbol{e}}}}} \right)}^{\text{T}}}{\boldsymbol{\sigma }}J{\text{d}}A} } \right) $$ (9) 其中：e为单元广义坐标。由式（9）可以给出薄膜单元的弹性力为
$$ {{\boldsymbol{F}}_e} = h\int_A {\frac{{\partial {{\boldsymbol{\varepsilon }}^{\text{T}}}}}{{\partial {\boldsymbol{e}}}}{\boldsymbol{\sigma }}J{\text{d}}A} $$ (10) 对式（9）求二次变分可得到
$$ \Delta \left( {\delta U} \right) = \delta {{\boldsymbol{e}}^{\text{T}}}\left( {h\int_A {\left( {{{\left( {\frac{{{\partial ^2}{\boldsymbol{\varepsilon }}}}{{\partial {{\boldsymbol{e}}^2}}}} \right)}^{\text{T}}}{\boldsymbol{\sigma }} + {{\left( {\frac{{\partial {\boldsymbol{\varepsilon }}}}{{\partial {\boldsymbol{e}}}}} \right)}^{\text{T}}}\frac{{\partial \sigma }}{{\partial \varepsilon }}\frac{{\partial {\boldsymbol{\varepsilon }}}}{{\partial {\boldsymbol{e}}}}} \right)J{\text{d}}A} } \right)\delta {\boldsymbol{e}} $$ (11) 进而得到薄膜单元的切线刚度矩阵的表达式为
$$ {\boldsymbol{K}} = h\int_A {\left( {\frac{{{\partial ^2}{{\boldsymbol{\varepsilon }}^{\text{T}}}}}{{\partial {{\boldsymbol{e}}^2}}}{\boldsymbol{\sigma }} + \frac{{\partial {{\boldsymbol{\varepsilon }}^{\text{T}}}}}{{\partial {\boldsymbol{e}}}}\frac{{\partial {\boldsymbol{\sigma }}}}{{\partial {\boldsymbol{\varepsilon }}}}\frac{{\partial {\boldsymbol{\varepsilon }}}}{{\partial {\boldsymbol{e}}}}} \right)J{\text{d}}A} $$ (12) 薄膜单元运动过程中的动能表达式为
$$ {E_k} = \frac{1}{2}\int_{{V^e}} \rho {{\boldsymbol{\dot r}}^{\text{T}}}{\boldsymbol{\dot r}}{\text{d}}{V^e} = \frac{h}{2}\int_A \rho {{\boldsymbol{\dot r}}^{\text{T}}}{\boldsymbol{\dot r}}J{\text{d}}A $$ (13) 其中；ρ为材料密度；
${\boldsymbol{\dot r}}$ 为绝对速度。从中可以提取出薄膜单元的质量矩阵Me的表达式为$$ {{\boldsymbol{M}}_{\text{e}}} = h\int_A \rho {{\boldsymbol{S}}^{\text{T}}}{\boldsymbol{S}}J{\text{d}}A $$ (14) -
张力场理论[16]认为薄膜的弯曲刚度为零，当薄膜内出现压应力时，一般通过面外变形即薄膜起褶的形式释放压应力，张力场理论将其处理为薄膜的面内收缩。起褶后的薄膜处于单轴拉伸应力状态，褶皱的方向平行于大主应力的方向，且垂直于褶皱方向的小主应力为零。
考虑到基于张力场理论的大部分方法在处理薄膜问题时，主要根据材料点的主应力状态对薄膜本构进行修正并确定褶皱区域。其本质上可以理解为薄膜的拉压不同模量问题。根据Ambartsumyan[10]建立的不同模量弹性理论，薄膜单元的广义弹性定律可表示为
$$ {{\boldsymbol{\varepsilon }}_{\text{I}}}{\text{ = }}{{\boldsymbol{a}}_{\text{I}}}{{\boldsymbol{\sigma }}_{\text{I}}} $$ (15) 其中：aI为柔度矩阵，主应力与主应变列阵分别记为
$$ {\boldsymbol \sigma }_{\rm I}={\left[\begin{array}{ccc}{\sigma }_{1}& {\sigma }_{2}& {\sigma }_{3}\end{array}\right]}^{\text{T}}\text{，}{\boldsymbol \varepsilon }_{I}={\left[\begin{array}{ccc}{\varepsilon }_{1}& { \varepsilon }_{2}& {\varepsilon }_{3}\end{array}\right]}^{\text{T}} $$ (16) 广义弹性定律的展开形式为
$$ \left\{ \begin{gathered} {\varepsilon _1}{\text{ = }}{a_{11}}{\sigma _1} + {a_{12}}{\sigma _2} + {a_{13}}{\sigma _3},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\varepsilon _{13}} = 0 \hfill \\ {\varepsilon _2}{\text{ = }}{a_{21}}{\sigma _1} + {a_{22}}{\sigma _2} + {a_{23}}{\sigma _3},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\varepsilon _{23}} = 0 \hfill \\ {\varepsilon _3}{\text{ = }}{a_{31}}{\sigma _1} + {a_{32}}{\sigma _2} + {a_{33}}{\sigma _3},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\varepsilon _{12}} = 0 \hfill \\ \end{gathered} \right. $$ (17) 其中：aij为柔度系数，每一项的具体形式取决于与之对应相乘的主应力的符号，即
$$ \left\{ \begin{array}{l}{\sigma }_{j} > 0\Rightarrow {a}_{jj}=1/{E}_{\rm t},{a}_{ij}=-{\upsilon }_{\rm t}/{E}_{\rm t}\\ {\sigma }_{j}\leqslant 0\Rightarrow {a}_{jj}=1/{E}_{\rm c},{a}_{ij}=-{\upsilon }_{\rm c}/{E}_{\rm c}\end{array}\left(i,j=1,2,3;且i\ne j\right) \right.$$ (18) 本文中采用Et、υt和Ec、υc分别表示拉伸和压缩状态下的杨氏模量和泊松比。二维工况下的广义弹性定律为
1）当σ1>0，σ2>0时
$$ \left[ {\begin{array}{*{20}{c}} {{\varepsilon _1}} \\ {{\varepsilon _2}} \end{array}} \right]=\left[ {\begin{array}{*{20}{c}} {{1 \mathord{\left/ {\vphantom {1 {{E_{{\rm t}}}}}} \right. } {{E_{\rm{t}}}}}}&{{{ - {\upsilon _{\rm t}}} \mathord{\left/ {\vphantom {{ - {\upsilon _{\rm {t}}}} {{E_{\rm tz]}}}}} \right. } {{E_{\rm{t}}}}}} \\ {{{ - {\upsilon _{\rm{t}}}} \mathord{\left/ {\vphantom {{ - {\upsilon _{\rm t}}} {{E_{\rm t}}}}} \right. } {{E_{\rm{t}}}}}}&{{1 \mathord{\left/ {\vphantom {1 {{E_{\rm{t}}}}}} \right. } {{E_{\rm{t}}}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\sigma _1}} \\ {{\sigma _2}} \end{array}} \right] $$ (19) 2）当σ1≤0，σ2≤0时
$$ \left[ {\begin{array}{*{20}{c}} {{\varepsilon _1}} \\ {{\varepsilon _2}} \end{array}} \right]=\left[ {\begin{array}{*{20}{c}} {{1 \mathord{\left/ {\vphantom {1 {{E_{\rm c}}}}} \right. } {{E_{\rm{c}}}}}}&{{{ - {\upsilon _{\rm{c}}}} \mathord{\left/ {\vphantom {{ - {\upsilon _{\rm{c}}}} {{E_{\rm{c}}}}}} \right. } {{E_{\rm{c}}}}}} \\ {{{ - {\upsilon _{\rm{c}}}} \mathord{\left/ {\vphantom {{ - {\upsilon _{\rm{c}}}} {{E_{\rm{c}}}}}} \right. } {{E_{\rm{c}}}}}}&{{1 \mathord{\left/ {\vphantom {1 {{E_{\rm{c}}}}}} \right. } {{E_{\rm{c}}}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\sigma _1}} \\ {{\sigma _2}} \end{array}} \right] $$ (20) 3）当σ1>0，σ2≤0时
$$ \left[ {\begin{array}{*{20}{c}} {{\varepsilon _1}} \\ {{\varepsilon _2}} \end{array}} \right]=\left[ {\begin{array}{*{20}{c}} {{1 \mathord{\left/ {\vphantom {1 {{E_{\rm{t}}}}}} \right. } {{E_{\rm{t}}}}}}&{{{ - {\upsilon _{\rm{c}}}} \mathord{\left/ {\vphantom {{ - {\upsilon _{\rm{c}}}} {{E_{\rm{c}}}}}} \right. } {{E_{\rm{c}}}}}} \\ {{{ - {\upsilon _{\rm{t}}}} \mathord{\left/ {\vphantom {{ - {\upsilon _{\rm{t}}}} {{E_{\rm{t}}}}}} \right. } {{E_{\rm{t}}}}}}&{{1 \mathord{\left/ {\vphantom {1 {{E_{\rm{c}}}}}} \right. } {{E_{\rm{c}}}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\sigma _1}} \\ {{\sigma _2}} \end{array}} \right] $$ (21) 4）当σ1≤0，σ2>0时
$$ \left[ {\begin{array}{*{20}{c}} {{\varepsilon _1}} \\ {{\varepsilon _2}} \end{array}} \right]=\left[ {\begin{array}{*{20}{c}} {{1 \mathord{\left/ {\vphantom {1 {{E_{\rm{c}}}}}} \right. } {{E_{\rm{c}}}}}}&{{{ - {\upsilon _{\rm{t}}}} \mathord{\left/ {\vphantom {{ - {\upsilon _{\rm{t}}}} {{E_{\rm{t}}}}}} \right. } {{E_{\rm{t}}}}}} \\ {{{ - {\upsilon _{\rm{c}}}} \mathord{\left/ {\vphantom {{ - {\upsilon _{\rm{c}}}} {{E_{\rm{c}}}}}} \right. } {{E_{\rm{c}}}}}}&{{1 \mathord{\left/ {\vphantom {1 {{E_t}}}} \right. } {{E_t}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\sigma _1}} \\ {{\sigma _2}} \end{array}} \right] $$ (22) 根据柔度矩阵与弹性矩阵的互逆关系，可得
$$ {{\boldsymbol{D}}_{\text{I}}}={\boldsymbol{a}}_{\text{I}}^{- 1}{\kern 1pt} {\kern 1pt} \Rightarrow {{\boldsymbol{\sigma }}_{\text{I}}}{\kern 1pt} ={{\boldsymbol{D}}_{\text{I}}}{{\boldsymbol{\varepsilon }}_{\text{I}}}{\kern 1pt} {\kern 1pt} $$ (23) 其中：DI为主应力空间的弹性矩阵。在直角坐标系中，广义胡克定律表示为
$$ {\boldsymbol{\sigma }}{\kern 1pt} {\text{ = }}{\boldsymbol{D\varepsilon }}{\kern 1pt} $$ (24) 其中：D为直角坐标系中的弹性矩阵。根据弹性力学中的转轴公式，给出如下转换关系
$$ {\boldsymbol{\sigma }}={\rm{R}}{{\boldsymbol{\sigma}} }_{\text{I}}\text{，}{{\boldsymbol \varepsilon }}_{\text{I}}={\boldsymbol R}^{\text{T}}{\boldsymbol{\varepsilon}} $$ (25) 其中：R为转换矩阵，其表达式为
$$ {\boldsymbol{R}} = \left[ {\begin{array}{*{20}{c}} {l_1^2}&{l_2^2} \\ {m_1^2}&{m_2^2} \\ {{l_1}{m_1}}&{{l_2}{m_2}} \end{array}} \right] $$ (26) 其中：li、mi为主应力方向余弦，i=1，2。联立式（23）~（25），可得
$$ {\boldsymbol{\sigma }}{\text{ = }}{\boldsymbol{R}}{{\boldsymbol{D}}_I}{{\boldsymbol{R}}^{\text{T}}}{\boldsymbol{\varepsilon }}{\kern 1pt} {\kern 1pt} {\kern 1pt} \Rightarrow {\boldsymbol{D}}{\text{ = }}{\boldsymbol{R}}{{\boldsymbol{D}}_I}{{\boldsymbol{R}}^{\text{T}}} $$ (27) 在传统的拉压不同模量算法中，确定应力状态下的不同模量本构关系通常被近似当作线性本构进行处理，即直接采用式（27）中所推导的弹性矩阵D作为应力对应变的偏导数， ANCF薄膜单元的弹性力和切线刚度矩阵可以表示为
$$ {{\boldsymbol{F}}_e} = h\int_A {\frac{{\partial {{\boldsymbol{\varepsilon }}^{\text{T}}}}}{{\partial {\boldsymbol{e}}}}{\boldsymbol{\sigma }}J{\text{d}}A} {\text{ = }}h\int_A {\frac{{\partial {{\boldsymbol{\varepsilon }}^{\text{T}}}}}{{\partial {\boldsymbol{e}}}}{\boldsymbol{D\varepsilon }}J{\text{d}}A} $$ (28) $$ {\boldsymbol{K}} = h\int_A {\left( {\frac{{{\partial ^2}{{\boldsymbol{\varepsilon }}^{\text{T}}}}}{{\partial {{\boldsymbol{e}}^{\text{2}}}}}{\boldsymbol{D\varepsilon }} + \frac{{\partial {{\boldsymbol{\varepsilon }}^{\text{T}}}}}{{\partial {\boldsymbol{e}}}}{\boldsymbol{D}}\frac{{\partial {\boldsymbol{\varepsilon }}}}{{\partial {\boldsymbol{e}}}}} \right)J{\text{d}}A} $$ (29) 实际工况中的薄膜结构在拉压过渡区往往会表现出强烈的非线性本构行为，应力会在张紧、松弛和褶皱几种状态下频繁切换，而式（29）不能连续地描述薄膜在这些应力状态下的变化过程，这也是影响薄膜问题计算效率的重要原因之一。张允真等[17]通过经典的双模量问题验证了这种弹性矩阵无法保证计算求解过程的收敛性和稳定性。
-
转换矩阵R依赖于薄膜单元的应变状态，式（27）推导得到的弹性矩阵D也是随着单元应力或者应变状态的变化而不断变化的，而式（29）中第二项中的推导结果显然忽略了弹性矩阵对应变的偏导数项，这种处理是不尽合理的。对式（29）做出进一步整理可得
$$ {\boldsymbol{K}} = h\int_A {\left( {\frac{{{\partial ^2}{{\boldsymbol{\varepsilon }}^{\text{T}}}}}{{\partial {{\boldsymbol{e}}^{\text{2}}}}}{\boldsymbol{D\varepsilon }} + \frac{{\partial {{\boldsymbol{\varepsilon }}^{\text{T}}}}}{{\partial {\boldsymbol{e}}}}{{\boldsymbol{D}}^{\rm{Tan} }}\frac{{\partial {\boldsymbol{\varepsilon }}}}{{\partial {\boldsymbol{e}}}}} \right)J{\text{d}}A} $$ (30) 其中：DTan为切线本构关系矩阵
$$ {{\boldsymbol{D}}^{{\rm Tan} }} = \frac{{\partial {\boldsymbol{\sigma }}}}{{\partial {\boldsymbol{\varepsilon }}}} = \frac{{\partial \left( {{\boldsymbol{D\varepsilon }}} \right)}}{{\partial {\boldsymbol{\varepsilon }}}} = {\boldsymbol{D}} + \frac{{\partial {\boldsymbol{D}}}}{{\partial {\boldsymbol{\varepsilon }}}}{\boldsymbol{\varepsilon }} $$ (31) 当主应力符号组合属于第一种类型时，对应的是纯拉伸或者纯压缩的工况，此时的研究对象可以看作均匀各向同性材料，式（30）中的切线本构关系矩阵DTan采用经典弹性理论中平面应力问题的弹性矩阵，即
$$ {{\boldsymbol{D}}^{{\rm Tan} }} = \frac{E}{{1 - {\upsilon ^2}}}\left[ {\begin{array}{*{20}{c}} 1&\upsilon &0 \\ \upsilon &1&0 \\ 0&0&{{{\left( {1 - \upsilon } \right)} \mathord{\left/ {\vphantom {{\left( {1 - \upsilon } \right)} 2}} \right. } 2}} \end{array}} \right]{\kern 1pt} {\kern 1pt} {\kern 1pt} $$ (32) 当主应力符号组合属于第二种类型时，对应的是拉压不同模量的工况，杜宗亮等[13]基于不同模量广义弹性定律给出了切线本构关系矩阵的具体形式
$$ {{\boldsymbol{D}}^{\rm{Tan} }} = {\boldsymbol{D}} + \left( {\frac{{\partial {\boldsymbol{R}}}}{{\partial {\boldsymbol{\varepsilon }}}}{{\boldsymbol{D}}_{\rm I}}{{\boldsymbol{R}}^{\text{T}}} + {\boldsymbol{R}}{{\boldsymbol{D}}_{\rm I}}\frac{{\partial {{\boldsymbol{R}}^{\text{T}}}}}{{\partial {\boldsymbol{\varepsilon }}}}} \right){\boldsymbol{\varepsilon }} $$ (33) 与式（27）相比，式（33）中多出的项包含了弹性矩阵在有限元迭代过程中的切线信息[13]。转换矩阵R对应变向量的偏导数可以通过对特征值和特征向量的灵敏度分析得到。
$$ \left\{ \begin{gathered} \frac{{\partial {l_1}}}{{\partial {\varepsilon _x}}} = \frac{{l_2^2{l_1}}}{{{\varepsilon _1} - {\varepsilon _2}}}{\kern 1pt} ,{\kern 1pt} \frac{{\partial {m_1}}}{{\partial {\varepsilon _x}}} = \frac{{{l_1}{l_2}{m_2}}}{{{\varepsilon _1} - {\varepsilon _2}}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \hfill \\ \frac{{\partial {l_1}}}{{\partial {\varepsilon _y}}} = \frac{{{l_2}{m_1}{m_2}}}{{{\varepsilon _1} - {\varepsilon _2}}},\frac{{\partial {m_1}}}{{\partial {\varepsilon _y}}} = \frac{{m_2^2{m_1}}}{{{\varepsilon _1} - {\varepsilon _2}}} \hfill \\ \frac{{\partial {l_1}}}{{\partial {\gamma _{xy}}}} = \frac{{{l_1}{l_2}{m_2} + l_2^2{m_1}}}{{2\left( {{\varepsilon _1} - {\varepsilon _2}} \right)}},\frac{{\partial {m_1}}}{{\partial {\gamma _{xy}}}} = \frac{{{m_1}{m_2}{l_2} + {l_1}m_2^2}}{{2\left( {{\varepsilon _1} - {\varepsilon _2}} \right)}} \hfill \\ \end{gathered} \right. $$ (34) $$ \left\{ \begin{gathered} \frac{{\partial {l_2}}}{{\partial {\varepsilon _x}}} = \frac{{l_1^2{l_2}}}{{{\varepsilon _2} - {\varepsilon _1}}}{\kern 1pt} ,\frac{{\partial {m_2}}}{{\partial {\varepsilon _x}}} = \frac{{{l_1}{l_2}{m_1}}}{{{\varepsilon _2} - {\varepsilon _1}}}{\kern 1pt} {\kern 1pt} {\kern 1pt} \hfill \\ {\kern 1pt} \frac{{\partial {l_2}}}{{\partial {\varepsilon _y}}} = \frac{{{l_1}{m_1}{m_2}}}{{{\varepsilon _2} - {\varepsilon _1}}},{\kern 1pt} {\kern 1pt} \frac{{\partial {m_2}}}{{\partial {\varepsilon _y}}} = \frac{{m_1^2{m_2}}}{{{\varepsilon _2} - {\varepsilon _1}}}{\kern 1pt} \hfill \\ \frac{{\partial {l_2}}}{{\partial {\gamma _{xy}}}} = \frac{{{l_1}{l_2}{m_1} + l_1^2{m_2}}}{{2\left( {{\varepsilon _2} - {\varepsilon _1}} \right)}},\frac{{\partial {m_2}}}{{\partial {\gamma _{xy}}}} = \frac{{{m_1}{m_2}{l_1} + {l_2}m_1^2}}{{2\left( {{\varepsilon _2} - {\varepsilon _1}} \right)}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \hfill \\ \end{gathered} \right. $$ (35) 通过经典的双模量算例[13]，上述模型在求解拉压不同模量小变形和大变形问题时均呈现出不错的收敛效率和计算精度。
-
根据刚柔多体系统的的约束条件，利用第一类拉格朗日方程可建立一组具有常数质量矩阵的微分代数方程
$$ \left\{ \begin{aligned} & {\boldsymbol{M}}{{{\boldsymbol{\ddot q}}}_{n + 1}} + {\boldsymbol{\varPhi }}_{\boldsymbol{q}}^{\rm T}{\boldsymbol{\lambda }} + {\boldsymbol{F}}\left( {{{\boldsymbol{q}}_{n + 1}}} \right) - {\boldsymbol{Q}}\left( {{{\boldsymbol{q}}_{n + 1}},{{{\boldsymbol{\dot q}}}_{n + 1}}} \right) = 0\\ & {\boldsymbol{\varPhi }}\left( {{{\boldsymbol{q}}_{n + 1}},{t_{n + 1}}} \right) = 0 \end{aligned} \right. $$ (36) 其中：M为系统质量矩阵；q为系统广义坐标向量；F为系统弹性力向量；Q为广义外力向量；λ为拉格朗日乘子向量；Φ和Φq分别为系统约束方程及其对广义坐标的Jacobi矩阵。
采用能够精确控制数值耗散的广义-α算法[18]对动力学方程进行离散时，需要在Newton-Raphon （N-R）迭代过程中多次求解如下高维线性代数方程组，式中Δq和Δλ为当前时刻系统的广义坐标向量和拉格朗日乘子向量的增量。
$$ \left[ {\begin{array}{*{20}{c}} {\boldsymbol{\varPsi }}&{{\boldsymbol{\varPhi }}_{\boldsymbol{q}}^{\rm T}} \\ {{\boldsymbol{\varPhi }}{}_{\boldsymbol{q}}}&{\boldsymbol{0}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\Delta {\boldsymbol{q}}} \\ {\Delta {\boldsymbol{\lambda }}} \end{array}} \right] = - \left[ {\begin{array}{*{20}{c}} {\boldsymbol{\varGamma }} \\ {\boldsymbol{\varPhi }} \end{array}} \right] $$ (37) 其中：
$$ {\boldsymbol{\varPsi }} = {\boldsymbol{M}}{\hat \beta} + {\left( {{\boldsymbol{\varPhi }}_{\boldsymbol{q}}^{\rm T}{\boldsymbol{\lambda }}} \right)_{\boldsymbol{q}}} + {\left[{{\boldsymbol{F}}\left( {\boldsymbol{q}} \right) - {\boldsymbol{Q}}\left( {{\boldsymbol{q}},{\boldsymbol{\dot q}}} \right)} \right]_{\boldsymbol{q}}} $$ (38) $$ {\boldsymbol{\varGamma }} = {\boldsymbol{M\ddot q}} + {\boldsymbol{\varPhi }}_{\boldsymbol{q}}^{\rm T}{\boldsymbol{\lambda }} + {\boldsymbol{F}}\left( {\boldsymbol{q}} \right) - {\boldsymbol{Q}}\left( {{\boldsymbol{q}},{\boldsymbol{\dot q}}} \right) $$ (39) 图2给出了薄膜结构的动力学分析流程图，在预处理模块中，根据薄膜结构的初始构形以及材料参数和单元信息等建立了精确的ANCF薄膜单元，给出了系统的常数质量矩阵。然后，采用高斯消元法消除冗余约束，得到系统的稀疏矩阵。通过假定初始时刻系统的应力场确定初始弹性矩阵，并组装系统的弹性力和雅克比矩阵，然后代入到式（36）~（39）的广义-α算法的求解过程中，根据求解结果更新系统位移场和弹性矩阵，如此循环迭代直至达到给定的收敛准则，求解完成。
图 2 膜结构动力学分析计算方法
Figure 2. Dynamic analysis and calculation method of membrane structure
-
为便于区分，将式（29）对应的线性近似的切线刚度模型记为ANCF LM，将式（30）对应的切线刚度模型记为ANCF TM。基于本文提出的计算方法，本节主要针对薄膜结构的一些经典案例进行了考察。4.1节、4.2节和4.3节分别验证文所提出的计算框架在求解不同模量小变形、大变形，以及面外变形问题时的高效性和精确性。
-
图3所示为一个验证薄膜起褶模型各种数值算法的经典算例，矩形薄膜宽度H，厚度t，其上下边界施加有大小σ0的均匀预应力，左右两端施加有轴向集中力P = σ0tH，同时左右边界均施加有集中力矩M。随着力矩的增加，膜的下边缘开始出现一定宽度的褶皱区域。考虑采用ANCF DM薄膜单元计算框架，对上述薄膜结构进行了数值模拟，仿真过程中，不考虑薄膜的压缩模量。
图 3 纯弯曲矩形薄膜受力示意图
Figure 3. Stress diagram of Pure Bending Rectangular membrane
文献[19]针对上述模型，采用解析方法对其起褶区域以及面内变形进行了研究，并得出一系列关系式。其中，式（40）给出了矩形膜起褶区域的宽度b与外力矩M的解析关系式
$$ \frac{b}{H} = \left\{ \begin{gathered} 0, \quad\quad\quad\quad \frac{M}{{PH}} < \frac{1}{6} \hfill \\ \frac{{3M}}{{PH}} - \frac{1}{2}{\kern 1pt} {\kern 1pt} ,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{6} \leqslant \frac{M}{{PH}} < \frac{1}{2} \hfill \\ \end{gathered} \right. $$ (40) 式（41）给出了不同工况下，薄膜右边界上一点的柯西应力与外力矩的解析关系式
$$ \frac{{2M}}{{PH}} = \left\{ \begin{gathered} \frac{1}{3}\frac{{E{H^2}t\kappa }}{{2P}}, \quad\quad\quad\quad \frac{{E{H^2}t\kappa }}{{2P}} \leqslant 1 \hfill \\ 1 - \frac{2}{3}\sqrt {\frac{{2P}}{{E{H^2}t\kappa }}} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{E{H^2}t\kappa }}{{2P}} > 1 \hfill \\ \end{gathered} \right.{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} $$ (41) 式（42）给出了薄膜变形后曲率κ与外力矩M之间的解析关系式
$$ \frac{{2M}}{{PH}} = \left\{ \begin{gathered} \frac{1}{3}\frac{{E{H^2}t\kappa }}{{2P}}, \quad\quad\quad\quad \frac{{E{H^2}t\kappa }}{{2P}} \leqslant 1 \hfill \\ 1 - \frac{2}{3}\sqrt {\frac{{2P}}{{E{H^2}t\kappa }}} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{E{H^2}t\kappa }}{{2P}} > 1 \hfill \\ \end{gathered} \right.{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} $$ (42) 图4给出了起褶区域宽度的解析解和数值试验结果，不难看出即使在起褶区域超过90%时，数值模拟结果仍然具有很高的精度。
图 4 褶皱带宽与外力矩的关系
Figure 4. Relationship between fold bandwidth and external torque
图5则给出了不同外力矩作用下，右边界上一点柯西应力的解析解与数值试验结果，两者的吻合性比较好。
图 5 右边界水平应力与外力矩的关系
Figure 5. Relationship between horizontal stress at right boundary and external moment
图6给出了薄膜结构曲率的解析解和数值解随外力矩的变化趋势，分析得知，两者的趋势是非常一致的。
图 6 纯弯曲矩形薄膜的弯矩–曲率关系
Figure 6. Moment-curvature relationship of rectangular membrane
通过上述各项数值试验结果与解析结果的对比可以得出，本文所提出的方法在研究薄膜起褶的问题时具有很高的精度。
-
如图7所示为一L形梁结构，其左端为固定约束，自由端受水平方向集中力F = 40 000 N。固定其受拉方向的弹性模量为Et = 3×107 Pa，泊松比为 vt = 0.3，依次分别取压缩模量为Ec = 3 × 107 Pa、Ec = 1.5 × 107 Pa、Ec = 6 ×106 Pa、Ec = 3× 106 Pa 以及Ec = 1.5 × 106 Pa 时的工况，对应的拉压比分别为Et/Ec = 1、2、5、10、20。
为验证本文所提出计算方法的有效性，本节基于ABAQUS平台的扩展程序接口进行了二次开发，并将UMAT的结果作为ANCF描述的不同模量问题的参考解。需要指出的是，UMAT这里采用CPS4I平面应力四边形单元和自适应时间增量步进行计算，在拉压比较大的情况下需要较多的加载步。而本文所提出的ANCF TM模型仅需一个增量步即可得到收敛解。
图 7 L形梁结构
Figure 7. L-shaped beam structure
图8给出了不同拉压模量比下，L形梁的变形后构形图，表1给出了采用不同技术途径在不同拉压模量比条件下得到的节点P的水平位移以及计算收敛耗时。以上仿真结果表明，即使在求解拉压模量相差20倍的大变形问题时，本文所提出的计算框架仍能高效收敛，且能够保证较高的精度。
图 8 变形后构形图
Figure 8. Configuration diagram after deformation
表 1 节点P的水平位移与计算耗时
Table 1. Horizontal displacement of node P and calculation convergence time
Et/Ec ANCF TM ABAQUS UMAT Ux/m T /s Ux/m T /s 1 6.774 33 13.8 6.789 95 3.1 2 7.229 50 14.7 7.213 68 7.0 5 7.696 99 17.2 7.645 53 11.3 10 7.970 26 17.8 7.898 76 14.7 20 8.193 32 20.3 8.107 53 17.1 此外，如表2所示，在采用ANCF-LM模型求解不同模量大变形问题时，无法得到收敛的结果。值得注意的是，即使拉压模量相同的条件下，ANCF-LM模型仍然无法得到收敛的结果，这是由于在相同模量时，该模型无法退化为经典的线性弹性本构模型。而本文提出的ANCF-TM模型求解过程中牛顿迭代收敛性良好，迭代步数和迭代耗时T随拉压模量比的增大略有增加。
表 2 不同模型下节点P的求解效率
Table 2. Solution efficiency of node P under different models
Et/Ec ANCF LM ANCF TM Uy/m 迭代步数 Uy/m 迭代步数 1 — — −6.139 83 16 2 — — −7.139 10 17 5 — — −8.237 28 20 10 — — −8.761 38 21 20 — — −9.028 77 24 -
图9所示为一方形安全气囊的初始构形，薄膜材料的杨氏模量为588.0 Mpa，泊松比为0.4。方膜对角线长度为L = 1 200 mm，厚度t = 1mm。考虑将气囊的充气过程近似为一个准静态过程，气囊在内部递增的均匀压力作用下发生变形，直至内部压强P最终增至5 000 Pa，充气过程结束。仿真过程中，不考虑薄膜的压缩模量。
图 9 安全气囊的初始构形示意图
Figure 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
网格密度 ANCF TM Oz/cm Bx/cm Cx/cm T/s 10×10 21.08 −11.63 −3.80 56.7 20×20 21.14 −11.93 −3.57 422.7 30×30 21.18 −11.92 −3.51 2248.9 图10出了安全气囊充气完成后的应力分布示意图，根据应力准则[13]，对单元内每个积分点的应力状态进行标记，用黑点标记的薄膜区域处于褶皱状态，用圆标记的薄膜区域处于松弛状态，其余未标记的区域则处于张紧状态，这一结果与文献[21]的结论基本一致。综上结论再次验证了本文所提出的计算框架在求解薄膜面外变形问题时的可行性。
图 10 安全气囊变形后应力分布示意图
Figure 10. Schematic diagram of stress distribution of airbag after inflation
-
作为一种新型空间可展开结构，能够通过自旋展开并最终稳定的大型帆膜结构在未来的空间探索技术应用中具有巨大潜力。能够准确地对薄膜展开过程的动态响应进行预测和模拟是薄膜太阳帆展开系统设计中的关键技术之一。这里将采用本文所提出的计算框架对薄膜太阳帆的自旋展开过程进行分析，通过确定影响展开的关键因素来指导太阳帆的结构设计。
初始时刻，薄膜太阳帆按照图11所示的折痕折叠后缠绕在与其同心的正六边形彀轮之上，彀轮驱动后，太阳帆随之旋转，并在均布于其六个角上的集中质量块的离心力作用下自旋展开，完全展开后的太阳帆是一个正六边形。此时太阳帆的边长为L = 0.635 m，厚度为t = 0.1 mm。薄膜材料的密度为ρ = 1 420 kg/m3，杨氏模量为E = 2.5×109 Pa，泊松比为v = 0.3。彀轮的边长为d = 0.1 m，质量为m0 = 1.0 kg，相对于过质心的Z轴的惯性矩为Iz = 5.0×10–3 kg m2。薄膜边界角点处的集中质量块的质量为m = 0.05 kg，初始时刻系统的角速度为w = 0.628 rad/s。
图 11 薄膜太阳帆的有限元模型
Figure 11. Finite element model of membrane solar sail
考虑到结构的对称性，本文仅建立1/6模型，并施加旋转边界条件，该有限元模型中包括35个薄膜片，薄膜片之间通过旋转铰连接，具体建模方法可以参考文献[8]。本文采用考虑拉压不同模量的ANCF薄膜单元和第一类拉格朗日方程对太阳帆进行动力学精确建模，然后利用可控数值耗散的广义-α算法对太阳帆的动力学展开进行求解分析。
-
在对太阳帆进行仿真时考察了3种不同的刚度折减的薄膜模型（Ec = s × Et），其中s为刚度折减系数。模型一为不考虑压缩刚度的不可压ANCF薄膜单元（s = 0）；模型二为考虑微弱抗压刚度的ANCF薄膜单元（s = 10–5）；模型三同样为引入微弱抗压刚度的ANCF薄膜单元（s = 10–4）。
图12给出了采用模型一进行仿真时薄膜太阳帆的展开过程图，结果表明，此时的薄膜太阳帆无法完全展开并维持展开状态，太阳帆在12 s左右展开至最大后又重新收缩至中心彀轮附近。
图 12 太阳帆的展开过程图（s=0）
Figure 12. Deployment process diagram of solar sail （s =0）
图13和图14分别给出了采用模型一进行仿真时太阳帆展开过程中的能量与角动量演化过程。图13表明，即使展开过程中薄膜的应变能出现高频振荡，但是系统的总能量仍然是守恒的。在12 s左右时太阳帆系统的应变能达到了峰值，此时太阳帆处于最大张紧状态。图14表明，相比于中心彀轮和集中质量块，薄膜在整个仿真过程中的角动量仅发生小范围的波动，在仿真时间达到12 s之后，太阳帆系统各部件的角动量均趋于稳定。
图 13 太阳帆展开过程中系统的应变能及总能量（s =0）
Figure 13. Strain energy and total energy of the system during solar sail deployment （s = 0）
图 14 太阳帆展开过程中系统各部件的角动量（s =0）
Figure 14. Angular momentum of each component of the system during solar sail deployment （s =0）
图15给出了采用模型二进行仿真时薄膜太阳帆的展开过程图，此时的薄膜太阳帆同样无法维持最大展开状态，太阳帆在11 s左右展开到最大半径，之后在边缘集中质量区域出现较大幅度的面外变形。根据图16和图17给出的太阳帆系统展开过程中的能量与和角动量的演化过程，可以得知，系统的应变能在12 s之后稳定在一较小的范围内。此外，在仿真时间达到11 s左右时，太阳帆系统各部件的角动量出现突变，之后逐渐稳定，展开过程中系统的总能量和总角动量均是守恒的。
图 15 太阳帆的展开过程图（s =10−5）
Figure 15. Deployment process diagram of solar sail（s=10−5）
图 16 太阳帆展开过程中系统的应变能及总能量（s =10−5）
Figure 16. Strain energy and total energy of the system during solar sail deployment（s =10−5）
图 17 太阳帆展开过程中系统各部件的角动量（s =10−5）
Figure 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)
Figure 18. Deployment process diagram of solar sail(s=10−4)
图 19 太阳帆展开过程中系统的应变能及总能量（s =10−4）
Figure 19. Strain energy and total energy of the system during solar sail deployment（s =10−4）
图21和图22给出了基于不同压缩刚度折减系数的薄膜模型仿真得到的太阳帆系统各部件的动力学响应。图21表明，当刚度折减系数为0和10–5时，太阳帆系统的彀轮和边缘质量块无法进入同步旋转。而在刚度折减系数为10–4的情况下，在仿真时间达到13 s左右时彀轮和边缘质量块的旋转角度开始相等，并进入同步旋转。图22表明，在考虑压缩刚度的情况下，当刚度折减系数较大时，薄膜太阳帆展开至一定半径大小之后趋于稳定。而当刚度折减系数较小或者不考虑压缩刚度时，薄膜太阳帆无法维持最大展开状态。综上分析可知，在对薄膜太阳帆系统进行展开动力学分析时，考虑一定的压缩刚度是必要的，这与文献[4]得出的结果基本是一致的。
图 20 太阳帆展开过程中系统各部件的角动量（s =10−4）
Figure 20. Angular momentum of each component of the system during solar sail deployment（s =10−4）
图 21 不同刚度折减系数时边界质量块和中心彀轮的转角
Figure 21. Rotation angle of edge mass and central hub under different stiffness reduction factors
图 22 不同刚度折减系数时的太阳帆展开半径
Figure 22. Deployment radius of solar sail under different stiffness reduction factors
-
为进一步确定薄膜太阳帆的一些结构参数对其展开动力学特性的影响，本文又分别仿真了当边界集中质量和系统初始转速变化时，太阳帆系统各部件的动力学响应。
图23和图24分别给出了太阳帆系统边界集中质量为0.03、0.05、0.07 kg时其各部件的转角演化过程和太阳帆展开半径演化过程。分析得知，太阳帆分别在16、13、11 s左右展开至最大半径，同时中心彀轮和边界质量块的旋转角度开始趋于同步，展开过程的构形演化图类似图18。此外，综合分析得知，边界集中质量对太阳帆系统的最终展开半径大小以及完全展开后发生的面外变形均会产生一定影响。
图 23 不同集中质量时边界质量块和中心彀轮的转角
Figure 23. Rotation angle of edge mass and central hub under different concentrated masses
图 24 不同集中质量时的太阳帆展开半径
Figure 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 不同初始转速时边界质量块和中心彀轮的转角
Figure 25. Rotation angle of edge mass and central hub at different initial speeds
图 26 不同初始转速时的太阳帆展开半径
Figure 26. Deployment radius of solar sail at different initial speeds
-
本文整合ANCF建模方法与张力场理论，针对薄膜结构拉压不同模量的宏观特性，提出了一种高效的建模与计算方法，为空间帆膜结构的展开动力学模拟提供了有效途径。
1）针对薄膜太阳帆系统展开过程中显著的非线性动力学特性，基于张力场理论与不同模量弹性理论，推导了单元精确的切向刚度矩阵，提出了一种考虑不同模量的ANCF薄膜单元计算框架。其中考虑了薄膜单元应力场对其本构关系的影响，准确地描述了薄膜结构在拉压过渡区的非线性特性，有效提高了单元在牛顿迭代过程中的收敛性。
2）基于上述框架，采用第一类拉格朗日方程建立了六边形薄膜太阳帆的自旋展开动力学模型，并通过可控数值耗散的广义-α算法求解运动方程。根据数值模拟结果，分析了太阳帆不同结构设计参数对其展开动力学特性的影响，为太阳帆类空间帆膜结构的系统设计和折展动力学仿真提供了一定的理论指导。
3）未来值得进一步研究的问题和方向包括：薄膜太阳帆完全展开之后，帆面产生的微弱面外变形以及振荡过程仍未得到有效解决；现有的针对薄膜太阳帆自旋展开的动力学建模与仿真方法在计算效率和分析精度方面仍然存在不足，随着薄膜单元网格的加密，难以针对大型帆膜结构进行高效模拟。
Spin Deployment Dynamics Simulation of Membrane Solar Sail Based on ANCF
-
摘要: 执行深空探测任务的薄膜太阳帆在入轨展开过程中，其薄膜结构在宏观力学响应上表现出显著的拉压不对称特性。系统的动力学行为呈现出强烈的非线性特征，使其动力学建模与仿真计算面对巨大挑战。基于绝对节点坐标方法（Absolute Nodal Coordinate Formulation，ANCF），整合张力场和不同模量弹性理论，推导了单元精确的切线刚度矩阵，提出了一种考虑不同模量的ANCF薄膜单元计算方法。在此基础上，采用第一类拉格朗日方程对薄膜太阳帆进行动力学建模，并通过可控数值耗散的广义-α算法求解动力学方程，分析了太阳帆不同结构设计参数对其展开动力学特性的影响。仿真结果验证了该方法在处理空间薄膜问题时的算法稳定性和高效性，对大型帆膜结构的系统设计有理论指导意义。Abstract: During the launch process of the membrane solar sail for deep space exploration mission, the membrane structure presents significant asymmetric characteristics under tension and compression. The dynamic behavior of the system shows strong nonlinear characteristics, which brings great challenges to the dynamic modeling and simulation calculation. Based on the Absolute Nodal Coordinate Formulation （ANCF）, integrating the tension field theory and the elasticity theory of different moduluses, the accurate tangent stiffness matrix of the element was derived, and a calculation method of ANCF membrane element considering different moduluses was proposed. On this basis, the first kind of Lagrange equation was used to model the dynamics of membrane solar sail, the dynamic equations were solved by the generalized-α algorithm with controllable numerical dissipation, and the influence of different structural design parameters of solar sail on its deployment dynamic characteristics was analyzed. The simulation results verify the stability and efficiency of this method in dealing with space membrane problems, and provide theoretical guidance for the system design of large sail membrane structures.
-
Key words:
- ANCF /
- membrane structure /
- elastic theory of different moduluses
Highlights● Based on ANCF， the wrinkling behavior of membrane structure is described from the viewpoint of different moduluses under tension and compression. ● Combined with elastic theory of different moduluses， the tangent stiffness matrix of membrane element based on ANCF is derived. ● A calculation framework of ANCF membrane element considering different moduluses is proposed and applied to the dynamic simulation of membrane solar sail. -
图 5 右边界水平应力与外力矩的关系
Fig. 5 Relationship between horizontal stress at right boundary and external moment
图 10 安全气囊变形后应力分布示意图
Fig. 10 Schematic diagram of stress distribution of airbag after inflation
图 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）
图 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）
图 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 不同刚度折减系数时边界质量块和中心彀轮的转角
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
图 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 不同初始转速时边界质量块和中心彀轮的转角
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
表 1 节点P的水平位移与计算耗时
Table 1 Horizontal displacement of node P and calculation convergence time
Et/Ec ANCF TM ABAQUS UMAT Ux/m T /s Ux/m T /s 1 6.774 33 13.8 6.789 95 3.1 2 7.229 50 14.7 7.213 68 7.0 5 7.696 99 17.2 7.645 53 11.3 10 7.970 26 17.8 7.898 76 14.7 20 8.193 32 20.3 8.107 53 17.1 
下载: 导出CSV
表 2 不同模型下节点P的求解效率
Table 2 Solution efficiency of node P under different models
Et/Ec ANCF LM ANCF TM Uy/m 迭代步数 Uy/m 迭代步数 1 — — −6.139 83 16 2 — — −7.139 10 17 5 — — −8.237 28 20 10 — — −8.761 38 21 20 — — −9.028 77 24
下载: 导出CSV
表 3 不同网格密度下的薄膜节点的位移
Table 3 Displacement of membrane nodes under different mesh densitied
网格密度 ANCF TM Oz/cm Bx/cm Cx/cm T/s 10×10 21.08 −11.63 −3.80 56.7 20×20 21.14 −11.93 −3.57 422.7 30×30 21.18 −11.92 −3.51 2248.9
下载: 导出CSV
-
[1] 彭福军,谢超,张良俊. 面向空间应用的薄膜可展开结构研究进展及技术挑战[J]. 载人航天,2017,23(4):427-437. doi: 10.3969/j.issn.1674-5825.2017.04.001PENG F J,XIE C,ZHANG L J. Research progress and technical challenges of membrane deployable structures for space applications[J]. Manned Spaceflight,2017,23(4):427-437. doi: 10.3969/j.issn.1674-5825.2017.04.001 [2] 周晓俊,周春燕,张新兴,等. 太阳帆自旋展开动力学地面模拟试验研究[J]. 振动工程学报,2015,28(2):175-182.ZHOU X J,ZHOU C Y,ZHANG X X,et al. Ground simulation experiment on spin deployment dynamics of solar sail[J]. Journal of Vibration Engineering,2015,28(2):175-182. [3] SHIRASAWA Y,MORI O,MIYAZAKI Y,et al. Evaluation of membrane dynamics of IKAROS based on flight result and simulation using multi-particle model[J]. Transactions of the Japan Society for Aeronautical and Spaceences Space Technology Japan,2012,10(28):421-426. [4] MIYAZAKI Y. Wrinkle/slack model and finite element dynamics of membrane[J]. International Journal for Numerical Methods in Engineering,2010,66(7):1179-1209. [5] 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 [6] 胡海岩. 太阳帆航天器的关键技术[J]. 深空探测学报(中英文),2016,3(4):334-342.HU H Y. Key technologies of solar sail spacecraft[J]. Journal of Deep Space Exploration,2016,3(4):334-342. [7] 赵将,刘铖,田强,等. 黏弹性薄膜太阳帆自旋展开动力学分析[J]. 力学学报,2013,45(5):746-752. doi: 10.6052/0459-1879-13-002ZHAO J,LIU C,TIAN Q,et al. Spin deployment dynamics analysis of viscoelastic membrane solar sail[J]. Chinese Journal of Theoretical and Applied Mechanics,2013,45(5):746-752. doi: 10.6052/0459-1879-13-002 [8] LIU C,TIAN Q,HU H Y,et al. Dynamic analysis of membrane systems undergoing overall motions,large deformations and wrinkles via thin shell elements of ANCF[J]. Computer Methods in Applied Mechanics and Engineering,2013,258:81-95. doi: 10.1016/j.cma.2013.02.006 [9] 田强,刘铖,李培,等. 多柔体系统动力学研究进展与挑战[J]. 动力学与控制学报,2017,15(5):385-405. doi: 10.6052/1672-6553-2017-039TIAN Q,LIU C,LI P,et al. Research progress and challenges of multi flexible body system dynamics[J]. Journal of Dynamics and Control,2017,15(5):385-405. doi: 10.6052/1672-6553-2017-039 [10] 阿姆巴尔楚米扬 S A. 不同模量弹性理论[M]. 邬瑞锋, 张允真, 译. 北京: 中国铁道出版社, 1986.AMBARTSUMYAN S A. Elasticity theory of different moduli[M]. Beijing: China Railway Press, 1986. [11] 叶志明,陈彤,姚文娟. 不同模量弹性问题理论及有限元法研究进展[J]. 力学与实践,2004,26(2):9-13. doi: 10.3969/j.issn.1000-0879.2004.02.002YE Z M,CHEN T,YAO W J. Research progress on theory and finite element method of elastic problems with different modulus[J]. Mechanics and Practice,2004,26(2):9-13. doi: 10.3969/j.issn.1000-0879.2004.02.002 [12] ZHANG H G,ZHANG L,GAO Q. An efficient computational method for mechanical analysis of bimodular structures based on parametric variational principle[J]. Computers and Structures,2011,89(23-24):2352-2360. doi: 10.1016/j.compstruc.2011.07.008 [13] DU Z L,ZHANG Y P,GUO X,et al. A new computational framework for materials with different mechanical responses in tension and compression and its applications[J]. International Journal of Solids and Structures,2016,100:54-73. [14] DUFVA K,SHABANA A A. Analysis of thin plate structures using the absolute nodal coordinate formulation[J]. Journal of Multi-body Dynamics,2005,219(4):345-355. [15] MIKKOLA A M,SHABANA A A. A non-incremental finite element procedure for the analysis of large deformation of plates and shells in mechanical system applications[J]. Multibody System Dynamics,2003,9(3):283-309. doi: 10.1023/A:1022950912782 [16] 李云良, 田振辉, 谭惠丰. 基于张力场理论的薄膜褶皱研究评述[J]. 力学与实践, 2008, 30(4):8−14.LI Y L, TIAN Z H, TAN H F. Review of methods of wrinkle studies based on tension field theory[J] Mechanics in Engineering, 2008, 30(4):8−14. [17] 张允真,王志锋. 不同拉、压模量弹性力学问题的有限元法[J]. 计算结构力学及其应用,1989,6(1):236-245.ZHANG Y Z,WANG Z F. Finite element method for elastic problems with different tensile and compressive modulus[J]. Computational Structural Mechanics and Application,1989,6(1):236-245. [18] YOU B D,LIANG D,YU X J,et al. Deployment dynamics for flexible deployable primary mirror of space telescope with paraboloidal and laminated structure by using absolute node coordinate method[J]. Chinese Journal of Aeronautics,2020,34(4):1-13. [19] MILLER R K,HEDGEPETH J M,WEINGARTEN V I,et al. Finite element analysis of partly wrinkled membranes[J]. Computers and Structures,1985,20(1):631-639. [20] KANG S,IM S. Finite element analysis of dynamic response of wrinkling membranes[J]. Computer Methods in Applied Mechanics and Engineering,1999,173(1-2):227-240. doi: 10.1016/S0045-7825(98)00271-0 [21] LEE E S,YOUN S K. Finite element analysis of wrinkling membrane structures with large deformations[J]. Finite Elements in Analysis and Design,2006,42(8-9):780-791. doi: 10.1016/j.finel.2006.01.004 -
点击查看大图
计量
- 文章访问数: 317
- HTML全文浏览量: 186
- PDF下载量: 53
- 被引次数: 0
