复合材料薄壁带簧刚度特性计算及折叠行为分析

Stiffness Properties and Folding Analysis of Composite Thin-walled Tape Spring

摘要

复合材料薄壁带簧结构被广泛应用于航空航天领域，开展复合材料层合板刚度特性和薄壁带簧结构折叠行为的研究具有重要工程价值。首先，基于经典层合板理论，推导描述复合材料层合板刚度特性的ABD刚度矩阵的理论计算方法。随后，采用商用有限元软件ABAQUS建立复合材料层合板代表性体积单元，通过施加周期性边界条件模拟无限大的复合材料层合板结构。根据推导的ABD刚度矩阵计算公式，向代表性体积单元施加不同的载荷，结合均匀化理论依次获取ABD刚度矩阵内各元素数值。最终，建立复合材料层合板ABD刚度矩阵的有限元快速计算方法，实现任意铺层顺序、铺层角度下ABD刚度矩阵的快速获取，且计算结果与理论计算结果一致。在获取复合材料层合板刚度特性后，开展带簧结构折叠行为数值模拟研究，探究使用ABAQUS中Static General、Dynamic Implicit两种隐式求解器和Dynamic Explicit显式求解器对复合材料薄壁带簧结构大变形弯曲折叠过程力学响应的影响。结果表明，两种隐式求解器和显式求解器均能获取较好的数值模拟结果，复合材料薄壁带簧结构力学响应与实验结果吻合。但其中显式求解方法不需要迭代和收敛准则，避免了大变形等非线性因素下计算收敛困难等问题，对于复合材料薄壁带簧弯曲、折叠等过程的数值模拟更具优势。

Abstract

Composite materials possess the advantages of light weight, high specific strength, and high specific stiffness, which have attracted considerable attention in recent years. Thin-walled flexible deployable structures made of composite materials have become one of the research hotspots in aerospace applications. The composite thin-walled tape spring is one of the typical space deployable structures that can serve as supporting components for solar sails, wrapped-rib antennas, and similar systems. The composite thin-walled tape spring stores strain energy during the folding process and uses this energy to drive itself and other components to achieve deployment in orbit. The stiffness properties of the composite laminate determine the mechanical response of the tape spring and affect the stored strain energy. In addition, the folding process of the tape spring is a complex process involving large deformations. Objective The investigation of the stiffness properties of the composite laminate and the folding behavior of the composite thin-walled tape spring holds significant engineering value. First, the theoretical calculation method of the ABD matrix describing the stiffness properties of composite laminates is derived based on classical laminated plate theory. However, when the layering order and layering angle of the composite lamina are changed, the ABD matrix must be recalculated to obtain the values of each element, which is a time-consuming process. Therefore, to make the process faster and more convenient for engineering applications, a finite element method for obtaining the ABD matrix of composite laminates is proposed. Methods The commercial finite element software ABAQUS was utilized to create a representative volume element of the composite laminate, which was modeled as a square structure with two sets of opposing edges. During the model discretization, it was necessary to ensure that the mesh nodes on the opposite edges of the representative volume element corresponded so that periodic boundary conditions can be correctly applied in subsequent steps. Two reference points were then established, and the nodes on the two sets of opposite edges of the representative volume element were correlated with these reference points by imposing the degree of freedom equation. Thus, the displacement relationship between each edge node and the corresponding opposite edge node was constrained through the reference points. The successful implementation of the periodic boundary conditions was achieved by establishing a linear equation on the reference points, ensuring that the structure deformed with identical shapes on the opposite edges, meaning that the overall structure represented by the representative volume element was ideal with an infinite array in space and no boundaries. After modeling, based on the derived ABD matrix formula, six different loads were applied to the reference points of the representative volume elements, corresponding to the six kinds of unit strain loading conditions. The reaction force and reaction moment on the reference points under six different loading conditions were then solved through static analysis. For each loading condition, the obtained reaction forces and moments were utilized to calculate the value of each element in a column of the ABD matrix. Finally, the values of all elements in the ABD matrix were sequentially obtained using homogenization theory. The results of the finite element method were compared to the theoretical results, and the maximum deviation of each element value in the ABD matrix was approximately 0.73%, verifying the correctness of the numerical calculation method developed for obtaining the ABD matrix. This approach can be extended to calculate the ABD matrix of composite laminates with any material parameters, ply sequences, or ply angles. After obtaining the stiffness properties of composite laminates, a numerical simulation of the folding behavior of the composite thin-walled tape spring was performed. Results and Discussions The tape spring structure has a length of 252 mm and a width of 40 mm, and was discretized using the shell element S8R5 in ABAQUS. The material properties of the tape spring were defined using the ABD matrix obtained previously through the developed method. Two reference points were created for kinematic coupling with the left and right edges of the tape spring to define the boundary conditions of the folding process. The rotational degrees of freedom of these reference points were linked to a virtual reference point, and the folding process was achieved by applying a rotational displacement to the virtual reference point to ensure uniform bending of both edges of the tape spring during the folding process. The influence of two implicit solvers (Static General and Dynamic Implicit) and one explicit solver (Dynamic Explicit) in the finite element software ABAQUS on the mechanical response of large deformation folding processes of the composite thin-walled tape spring was investigated. The implicit algorithm approached the equilibrium solution through continuous iteration, which has a very high computational cost for completing one iteration, and often encountered convergence difficulties for strongly nonlinear problems such as large deformations and contacts. The explicit algorithm employed the central difference method and used the result of one incremental step as the initial condition for calculating the next incremental step. This method did not require iteration or a convergence criterion and handled complex contacts and extreme discontinuities more effectively. However, the time increment step of this method needed to be sufficiently small to ensure the accuracy of the numerical simulation results. The change curve of the folding moment and curvature at the center of the composite thin-walled tape spring during the folding process was obtained through numerical simulation analysis under different solvers. Conclusions The numerical simulation results indicate that both implicit and explicit solvers can achieve accurate outcomes, and the mechanical responses of the composite thin-walled tape spring correspond well with the experimental results. Although the implicit solution method effectively addresses the numerical calculation challenges of the thin-walled tape spring under large deformation, it requires numerous iterations during the solution process. In contrast, the explicit solution method does not require iteration or convergence criteria, avoiding computational convergence issues caused by nonlinear factors such as large deformation. Therefore, it is more advantageous for numerically simulating the bending and folding processes of composite thin-walled tape springs.

Graphical abstract

关键词

复合材料 / 薄壁结构 / ABD刚度矩阵 / 有限元分析

Key words

composites / thin-walled structures / ABD matrix / finite element analysis

本刊网刊

复合材料具有轻量化、高比强度和高比刚度等特点，被广泛应用于航天航空领域^[1]。近年来，为了满足大型空间结构如卫星星载大口径天线等的发展需求,由复合材料制成的薄壁柔性可展开结构受到了广泛关注，其结构特点是能在不发生复合材料失效破坏的前提下折叠收拢，减小结构体积，并在在轨后展开，提供一定的承载能力^[2]。其中，复合材料薄壁带簧结构是一种典型的空间可展开单元^[3‒5]，可用于太阳帆^[6‒7]、缠绕肋天线等结构的支撑部件^[8‒10]，并依靠带簧结构自身折叠过程中存储的应变能驱动太阳帆薄膜以及天线反射面等结构进行展开。此外，带簧结构也可作为伸展臂，实现大尺寸天线等杆件的展开^[11‒12]。2003年6月2日，由欧洲航天局ESA负责研制的“火星快车”号火星探测飞行器成功发射，其中搭载的雷达携带有由两根复合材料薄壁带簧铰链制成的偶极子天线杆^[13]，但其在轨展开过程中出现展开停滞等问题，并未成功展开到位。此次事故引发了对复合材料力学性能以及复合材料薄壁柔性可展开结构折叠收拢/展开过程的深入研究。

复合材料薄壁柔性可展开结构一般由复合材料层合板制作完成。针对复合材料层合板最为关键的性能表征，提出一种准确预测复合材料层合板刚度特性ABD刚度矩阵的方法成为了众多学者的研究方向。现有的预测模型可以大致分为两类：解析模型与数值模型。Ishikawa和Chou^[14‒15]将经典层合板理论结合等应变或等应力假设提出了编织复合材料均匀化ABD刚度矩阵预测模型，后续学者在此基础上又进一步提出了一些解析模型^[16‒17]。尽管解析模型具有预测快、节约计算成本的优点，但是大部分解析模型只在面内力学性能预测上有较好表现，在面外或者剪切方向的预测精度不足。相较于解析模型，数值模型对计算资源要求更高，求解时间更长，但是随着计算机技术的不断发展，基于代表性体积单元（RVE）的多尺度有限元分析的数值模型在复合材料层合板刚度特性的预测中受到了广泛的关注。An等^[18]基于商用有限元软件ABAQUS开发了一种包含纤维和基体的复合材料微观尺度建模方法，通过施加周期性边界条件结合均匀化方法成功获取了复合材料单层板的材料参数。Omairey^[19]、Daniel^[20]和Quek^[21]等将编织复合材料层合板结构划分为3个尺度：包含纤维和基体的微观尺度、包含纱线的中观尺度，以及整体结构的宏观尺度。整个复合材料层合板ABD刚度矩阵预测过程可划分为两次均匀化：第一次均匀化获得单根纱线的正交各向异性本构关系；第二次均匀化将微观尺度均匀化获得的纱线工程常数输入至中观尺度下的代表性体积单元有限元模型。随后，施加周期性边界条件并进行有限元计算来获得ABD刚度矩阵。最后，将ABD刚度矩阵传递至宏观尺度下的壳单元有限元模型后完成对整体结构力学响应的分析。

在准确获取复合材料层合板刚度特性后，针对复合材料薄壁带簧结构的弯曲折叠等大变形问题，大量学者开展了实验和数值模拟研究。Mallikarachchi等^[22‒23]首先通过实验研究了薄壁带簧结构在准静态折叠过程中力矩随转角变化情况以及动力学展开过程中带簧结构的响应，如展开角度随时间变化情况等。闫中曦等^[24]建立了带簧结构弯曲时的应变能与弯矩的解析模型，通过理论模型与仿真分析两种方法计算了带簧结构在纯弯作用下的力矩响应，并研究了几何尺寸等参数对力矩和应力的影响规律。Leclerc等^[25‒26]采用商用有限元软件ABAQUS中用于非线性静力分析的牛顿‒辛普森隐式求解算法，进行了复合材料薄壁结构在纯弯作用下的屈曲分析，数值计算发现在结构后屈曲计算末段出现了数值发散，难以计算复合材料薄壁结构坍塌过程的力学响应。复合材料薄壁结构弯曲折叠过程中涉及复杂的不稳定后屈曲和坍塌等问题，存在复杂的接触和显著的不连续性，难以采用隐式方法迭代解决，而显式求解算法因其不需要迭代和收敛准则，在处理此类问题上更具优势。Seffen等^[27]采用显式动力学求解技术进行了复合材料薄壁带簧结构大变形折叠过程数值模拟，验证了基于显式动力学的数值模拟方法在求解复合材料薄壁结构大变形等非线性问题上的可靠性。此外，在对复合材料薄壁带簧结构折叠及展开过程数值模拟的基础上，Ye^[28‒29]和Jin^[30‒31]等学者针对折叠和动力学展开过程中的关键力学响应开展了带簧结构的尺寸优化和形状优化工作。

本文首先从经典层合板理论出发，推导描述复合材料层合板刚度特性的ABD刚度矩阵，随后在此基础上建立了一种基于有限元数值模拟方法的复合材料层合板ABD刚度矩阵快速计算方法，并与解析结果进行对比验证。利用上述方法计算获取的ABD刚度矩阵建立复合材料薄壁带簧有限元模型，采用商用有限元软件ABAQUS中隐式和显式求解方法进行大变形弯曲折叠过程数值模拟，探究不同求解器对数值模拟计算的影响。

1 复合材料层合板刚度特性计算

1.1 经典层合板理论

复合材料层合板由两层或两层以上的复合材料单层板粘合在一起成为整体的层合式结构元件。复合材料层合板可以由相同或不同材质的复合材料单层板沿着一定的铺设方向构成，其受力特性和各单层板紧密关联，呈现出各向异性的性质。复合材料单层板本构关系，即应力‒应变关系如方程（1）所示：

\left[\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right]=\left[\mathit{Q}\right]\left[\begin{array}{c}{\epsilon }_x\\ {\epsilon }_y\\ {\gamma }_{xy}\end{array}\right]（1）

式中，{\sigma }_x、{\sigma }_y、{\tau }_{xy}为应力分量， Q 为自然坐标系下的刚度矩阵，{\epsilon }_x^{}、{\epsilon }_y^{}、{\gamma }_{xy}^{}为应变分量。

图1所示的复合材料层合板由N层任意铺层角度的复合材料单层板构成。其自然坐标系中，xoy面与层合板中面重合，z轴垂直于层合板，层合板总厚度为t。图1中，u、v、w为位移分量，z_k 为中性面至第k层合板厚度。

由弹性力学可知，沿复合材料层合板厚度变化的应变方程为：

\left[\begin{array}{c}{\epsilon }_x\\ {\epsilon }_y\\ {\gamma }_{xy}\end{array}\right]=\left[\begin{array}{c}{\epsilon }_x^{\mathrm{0}}\\ {\epsilon }_y^{\mathrm{0}}\\ {\gamma }_{xy}^{\mathrm{0}}\end{array}\right]+z\left[\begin{array}{c}{\kappa }_x\\ {\kappa }_y\\ {\kappa }_{xy}\end{array}\right]=\left[\begin{array}{c}\frac{\partial u_{\mathrm{0}}}{\partial x}\\ \frac{\partial v_{\mathrm{0}}}{\partial y}\\ \frac{\partial u_{\mathrm{0}}}{\partial y}+\frac{\partial v_{\mathrm{0}}}{\partial x}\end{array}\right]-z\left[\begin{array}{c}\frac{{\partial }^{\mathrm{2}}w}{\partial x^{\mathrm{2}}}\\ \frac{{\partial }^{\mathrm{2}}w}{\partial y^{\mathrm{2}}}\\ \frac{\mathrm{2}{\partial }^{\mathrm{2}}w}{\partial x\partial y}\end{array}\right]（2）

式中，{\epsilon }_x^{\mathrm{0}}、{\epsilon }_y^{\mathrm{0}}、{\gamma }_{xy}^{\mathrm{0}}为中面应变分量，u_{\mathrm{0}}、v_{\mathrm{0}}、w为中面位移分量，{\kappa }_{xy}、{\kappa }_y、{\kappa }_{xy}为中面曲率分量。

将式（2）代入式（1），可得复合材料层合板第k层应力使用中面应变和曲率的表达式：

{\left[\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right]}_k={\left[\begin{array}{ccc}{Q}_{\mathrm{11}}& Q_{\mathrm{12}}& Q_{\mathrm{16}}\\ Q_{\mathrm{12}}& Q_{\mathrm{22}}& Q_{\mathrm{26}}\\ Q_{\mathrm{16}}& Q_{\mathrm{26}}& Q_{\mathrm{66}}\end{array}\right]}_k\left[\left[\begin{array}{c}{\epsilon }_x^{\mathrm{0}}\\ {\epsilon }_y^{\mathrm{0}}\\ {\gamma }_{xy}^{\mathrm{0}}\end{array}\right]+z\left[\begin{array}{c}{\kappa }_x\\ {\kappa }_y\\ {\kappa }_{xy}\end{array}\right]\right]（3）

式中，Q_ij 为式（1）中刚度矩阵 Q 的展开。

将复合材料层合板截面内应力沿层合板厚度方向积分可得层合板单位宽度上的合力，如式（4）所示：

\left[\begin{array}{c}{N}_x\\ N_y\\ N_{xy}\end{array}\right]={\int }_{-t/\mathrm{2}}^{t/\mathrm{2}}\left[\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right]\mathrm{d}z={\displaystyle \sum _{k=\mathrm{1}}^N}{\int }_{z_{k-\mathrm{1}}}^{z_k}\left[\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right]\mathrm{d}z（4）

式中，N_x、N_y、N_{xy}为单位长度下的合力。

同理，将复合材料层合板截面内应力与厚度方向的坐标z的乘积沿厚度方向积分可得层合板单位宽度上的合力矩，如式（5）所示：

\left[\begin{array}{c}{M}_x\\ M_y\\ M_{xy}\end{array}\right]={\int }_{-t/\mathrm{2}}^{t/\mathrm{2}}\left[\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right]z\mathrm{d}z={\displaystyle \sum _{k=\mathrm{1}}^N}{\int }_{z_{k-\mathrm{1}}}^{z_k}\left[\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right]z\mathrm{d}z（5）

式中，M_x、M_y、M_{xy}为单位长度下的合力矩。

最终合并式（3）～（5）并化简得到描述复合材料层合板力及力矩与应变及曲率关系的ABD刚度矩阵：

\left[\begin{array}{c}\begin{array}{c}{N}_x\\ N_y\\ N_{xy}\end{array}\\ \begin{array}{c}{M}_x\\ M_y\\ M_{xy}\end{array}\end{array}\right]=\left[\begin{array}{cc}\begin{array}{ccc}{A}_{\mathrm{11}}& A_{\mathrm{12}}& A_{\mathrm{16}}\\ A_{\mathrm{12}}& A_{\mathrm{22}}& A_{\mathrm{26}}\\ A_{\mathrm{16}}& A_{\mathrm{26}}& A_{\mathrm{66}}\end{array}& \begin{array}{ccc}{B}_{\mathrm{11}}& B_{\mathrm{12}}& B_{\mathrm{16}}\\ B_{\mathrm{12}}& B_{\mathrm{22}}& B_{\mathrm{26}}\\ B_{\mathrm{16}}& B_{\mathrm{26}}& B_{\mathrm{66}}\end{array}\\ \begin{array}{ccc}{B}_{\mathrm{11}}& B_{\mathrm{12}}& B_{\mathrm{16}}\\ B_{\mathrm{12}}& B_{\mathrm{22}}& B_{\mathrm{26}}\\ B_{\mathrm{16}}& B_{\mathrm{26}}& B_{\mathrm{66}}\end{array}& \begin{array}{ccc}{D}_{\mathrm{11}}& D_{\mathrm{12}}& D_{\mathrm{16}}\\ D_{\mathrm{12}}& D_{\mathrm{22}}& D_{\mathrm{26}}\\ D_{\mathrm{16}}& D_{\mathrm{26}}& D_{\mathrm{66}}\end{array}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{\epsilon }_x^{\mathrm{0}}\\ {\epsilon }_y^{\mathrm{0}}\\ {\gamma }_{xy}^{\mathrm{0}}\end{array}\\ \begin{array}{c}{\kappa }_x\\ {\kappa }_y\\ {\kappa }_{xy}\end{array}\end{array}\right]（6）

式中，A_{ij}为拉伸刚度矩阵系数，B_{ij}为耦合刚度矩阵系数，D_{ij}为弯曲刚度矩阵系数。A_{ij}、B_{ij}和D_{ij}的计算分别如式（7）~（9）所示：

A_{ij}={\displaystyle \sum _{k=\mathrm{1}}^N}(Q_{ij}{)}_k(z_k-z_{k-\mathrm{1}})={\displaystyle \sum _{k=\mathrm{1}}^N}(Q_{ij}{)}_k{t}_k（7）

B_{ij}=\frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \sum _{k=\mathrm{1}}^N}(Q_{ij}{)}_k\left(z_k^{\mathrm{2}}-z_{k-\mathrm{1}}^{\mathrm{2}}\right)={\displaystyle \sum _{k=\mathrm{1}}^N}(Q_{ij}{)}_k{t}_k{\overline{z}}_k（8）

D_{ij}=\frac{\mathrm{1}}{\mathrm{3}}{\displaystyle \sum _{k=\mathrm{1}}^N}(Q_{ij}{)}_k\left(z_k^{\mathrm{3}}-z_{k-\mathrm{1}}^{\mathrm{3}}\right)={\displaystyle \sum _{k=\mathrm{1}}^N}(Q_{ij}{)}_k\left(t_k{\overline{z}}_k^{\mathrm{2}}+\frac{t_k^{\mathrm{3}}}{\mathrm{12}}\right)（9）

式（7）～（9）中，t_k为第k层的厚度，{\overline{z}}_k为第k层z方向中心的坐标值。

假设求解过程中应力单位为GPa，厚度单位为mm，则式（6）中A_{ij}单位为\mathrm{k}\mathrm{N}/\mathrm{m}\mathrm{m}，B_{ij}的单位为\mathrm{k}\mathrm{N}，D_{ij}的单位为\mathrm{k}\mathrm{N}\cdot \mathrm{m}\mathrm{m}。通过式（6）能够计算复合材料层合板的刚度特性，但更改铺层顺序、铺层角度后，需要重新进行复合材料ABD刚度矩阵推导，获取矩阵内元素数值，该过程较为烦琐和耗时。因此，本文开发了一种复合材料ABD刚度矩阵有限元快速获取方法，能够根据任意铺层顺序、铺层角度快速给出复合材料ABD刚度矩阵各元素数值，服务于工程应用且便于后续进行复合材料薄壁带簧结构数值模拟。

1.2 ABD刚度矩阵有限元计算

本文开发了一种适用于任意铺层形式、材料参数的复合材料层合板ABD刚度矩阵计算方法。以玻璃纤维及聚丙烯材料所制成的复合材料层合板为例，单层板材料属性如表1所示^[32]，单层板厚度为0.213 mm，铺层方式为[+{\mathrm{45}}^{\mathrm{o}}/-{\mathrm{45}}^{\mathrm{o}}/{\mathrm{0}}^{\mathrm{o}}/+{\mathrm{45}}^{\mathrm{o}}/-{\mathrm{45}}^{\mathrm{o}}]。表1中，E₁、E₂分别为主方向（纤维方向）和垂直于主方向（横向）的弹性模量，G₁₂为面内剪切弹性模量，v₁₂、v₂₃为泊松比，分别代表纵向应变引起的横向应变和横向应变引起的纵向应变。

采用有限元软件ABAQUS建立复合材料层合板代表性体积单元（RVE），施加周期性边界条件（PBC），之后通过静力分析依次对代表性体积单元施加单一方向的单位应变，提取支反力和支反力矩，并最终求得ABD刚度矩阵各元素数值。具体步骤如下。

首先，建立任意尺寸几何模型作为代表性体积单元，如图2所示，该代表性体积单元为边长l_x=l_y的正方形平面，将四边依次编号为1～4。按照表1中相关材料性能参数及反对称铺层方式，在商用有限元软件ABAQUS中定义平面壳结构的截面属性。进行网格划分时，确保相对边上的网格节点一一对应。

建立虚拟参考点A、虚拟参考点B，分别与代表性体积单元两组对边各相对应节点进行自由度方程关联，并在参考点上建立线性方程对相对边上的对应节点施加约束，完成周期性边界条件的施加。实现结构变形后，两组对边1和3、2和4形状完全一致，即代表性体积单元模型所代表的整体结构是一种空间中无限阵列、没有边界的理想结构。

随后，进行ABD刚度矩阵各元素求解。沿x方向向虚拟参考点A施加拉伸位移载荷，并约束虚拟参考点A、虚拟参考点B其他方向自由度，使代表性体积单元仅产生沿x方向的单位应变，即{\epsilon }_x^{\mathrm{0}}=\mathrm{1}，{\epsilon }_y^{\mathrm{0}}={\gamma }_{xy}^{\mathrm{0}}={\kappa }_x^{}={\kappa }_y^{}={\kappa }_{xy}^{}=\mathrm{0}，此时可得：

\left[\begin{array}{c}\begin{array}{c}{N}_x\\ N_y\\ N_{xy}\end{array}\\ \begin{array}{c}{M}_x\\ M_y\\ M_{xy}\end{array}\end{array}\right]=\left[\begin{array}{cc}\begin{array}{ccc}{A}_{\mathrm{11}}& A_{\mathrm{12}}& A_{\mathrm{16}}\\ A_{\mathrm{12}}& A_{\mathrm{22}}& A_{\mathrm{26}}\\ A_{\mathrm{16}}& A_{\mathrm{26}}& A_{\mathrm{66}}\end{array}& \begin{array}{ccc}{B}_{\mathrm{11}}& B_{\mathrm{12}}& B_{\mathrm{16}}\\ B_{\mathrm{12}}& B_{\mathrm{22}}& B_{\mathrm{26}}\\ B_{\mathrm{16}}& B_{\mathrm{26}}& B_{\mathrm{66}}\end{array}\\ \begin{array}{ccc}{B}_{\mathrm{11}}& B_{\mathrm{12}}& B_{\mathrm{16}}\\ B_{\mathrm{12}}& B_{\mathrm{22}}& B_{\mathrm{26}}\\ B_{\mathrm{16}}& B_{\mathrm{26}}& B_{\mathrm{66}}\end{array}& \begin{array}{ccc}{D}_{\mathrm{11}}& D_{\mathrm{12}}& D_{\mathrm{16}}\\ D_{\mathrm{12}}& D_{\mathrm{22}}& D_{\mathrm{26}}\\ D_{\mathrm{16}}& D_{\mathrm{26}}& D_{\mathrm{66}}\end{array}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}\mathrm{1}\\ \mathrm{0}\\ \mathrm{0}\end{array}\\ \begin{array}{l}\mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\end{array}\end{array}\right]（10）

由式（10）可知，在该工况下N_x=A_{\mathrm{11}}，N_y=A_{\mathrm{12}}，N_{xy}=A_{\mathrm{16}}，M_x=B_{\mathrm{11}}，M_y=B_{\mathrm{12}}，M_{xy}=B_{\mathrm{16}}。之后通过求解静力分析，提取虚拟参考点A、虚拟参考点B在x、y方向的支反力和支反力矩，并根据式（11）计算该代表性体积单元所对应的单位长度下合力N_x\mathrm{、}{N}_y\mathrm{、}{N}_{xy}和合力矩M_x\mathrm{、}{M}_y\mathrm{、}{M}_{xy}：

\left\{\begin{array}{l}{N}_x=\frac{F_x^{\mathrm{A}}}{l_y},N_y=\frac{F_y^{\mathrm{B}}}{l_x},N_{xy}=\frac{F_y^{\mathrm{A}}}{l_y};\\ M_x=\frac{M_y^{\mathrm{B}}}{l_y},M_y=\frac{M_x^{\mathrm{A}}}{l_x},M_{xy}=\frac{M_x^{\mathrm{A}}}{l_y}\end{array}\right.（11）

将以上单位长度下合力和合力矩的值代入式（10），可求得ABD刚度矩阵第1列元素的值。

重复此求解过程，依次令y方向拉伸应变{\epsilon }_y^{\mathrm{0}}=\mathrm{1}、x-y方向剪切应变{\epsilon }_{xy}^{\mathrm{0}}=\mathrm{1}、x方向曲率{\kappa }_x=\mathrm{1}、y方向曲率{\kappa }_y=\mathrm{1}以及x-y方向曲率{\kappa }_{xy}=\mathrm{1}，可获得ABD刚度矩阵的第2～6列中各元素的值。图3为不同方向单位应变加载下代表性体积单元变形示意图。

最终，求得整个ABD刚度矩阵中每一项的数值，由于数值计算获取的部分元素值相比其他元素值小7～8个数量级，因此令这些小量约等于0，进而得到复合材料薄壁壳结构ABD刚度矩阵的有限元等效计算结果为：

\left[\begin{array}{cccccc}\mathrm{13.885}& \mathrm{6.004}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{0.546}\\ \mathrm{6.004}& \mathrm{8.758}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{0.546}\\ \mathrm{0}& \mathrm{0}& \mathrm{6.194}& -\mathrm{0.546}& -\mathrm{0.546}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -\mathrm{0.546}& \mathrm{0.971}& \mathrm{0.673}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -\mathrm{0.546}& \mathrm{0.673}& \mathrm{0.952}& \mathrm{0}\\ -\mathrm{0.546}& -\mathrm{0.546}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0.692}\end{array}\right]（12）

此外，将表1中给出的材料性能参数及铺层角度代入式（6）中，可得到复合材料层合板ABD刚度矩阵的理论解：

\left[\begin{array}{cccccc}\mathrm{13.90}& \mathrm{6.00}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{0.55}\\ \mathrm{6.00}& \mathrm{8.80}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{0.55}\\ \mathrm{0}& \mathrm{0}& \mathrm{6.19}& -\mathrm{0.55}& -\mathrm{0.55}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -\mathrm{0.55}& \mathrm{0.97}& \mathrm{0.67}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -\mathrm{0.55}& \mathrm{0.67}& \mathrm{0.95}& \mathrm{0}\\ -\mathrm{0.55}& -\mathrm{0.55}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0.69}\end{array}\right]（13）

通过有限元软件仿真分析计算复合材料层合板均匀化刚度性能并与经典层合板理论推导结果对比，即式（12）与式（13）对比，可以看到二者误差较小，矩阵内各元素值最大偏差为0.73%，验证了本文所开发数值计算获取ABD刚度矩阵方法的正确性，能推广到任意铺层形式和材料参数的复合材料层合板ABD刚度矩阵快速计算。

2 复合材料薄壁带簧弯曲折叠数值模拟

2.1 计算模型

复合材料薄壁带簧是一种可弯曲折叠的薄壁结构，被广泛应用于空间可展开结构中。首先，在商业有限元软件ABAQUS中建立长度为252 mm、宽度为40 mm的复合材料薄壁带簧模型，复合材料属性及铺层信息见表1，利用上述开发的复合材料层合板ABD刚度矩阵有限元计算方法获取描述该复合材料薄壁带簧刚度特性的ABD刚度矩阵，并采用壳单元（S8R5）进行有限元模型构建，模型示意图及边界条件如图4所示。参考点A、B分别与带簧两边进行运动耦合绑定。参考点A除绕x轴弯曲外其余自由度均0，即位移U_x^{\mathrm{A}}=U_y^{\mathrm{A}}=U_z^{\mathrm{A}}=\mathrm{0}，转动自由度U{R}_y^{\mathrm{A}}=U{R}_z^{\mathrm{A}}=\mathrm{0}。参考点B能沿z轴平动及绕x轴弯曲，其余自由度均为0，即U_x^{\mathrm{B}}=U_y^{\mathrm{B}}=\mathrm{0}，转动自由度U{R}_y^{\mathrm{B}}=U{R}_z^{\mathrm{B}}=\mathrm{0}。参考点A、B绕x轴的转动自由度通过虚拟参考点C相关联以实现加载过程中复合材料薄壁带簧两侧均匀弯曲，即U{R}_x^{\mathrm{A}}-U{R}_x^{\mathrm{B}}=U{R}_x^{\mathrm{C}}。虚拟参考点C的边界条件为U_x^{\mathrm{C}}=U_y^{\mathrm{C}}=U_z^{\mathrm{C}}=\mathrm{0}，U{R}_y^{\mathrm{C}}=U{R}_z^{\mathrm{C}}=\mathrm{0}。逐步加载转角度位移载荷U{R}_x^{\mathrm{B}}=\mathrm{0}，使复合材料薄壁带簧绕x轴的弯曲转角θ_x 在一定的步长时间内从0°增加到170°（或从0°减小到-170°）。

2.2 求解器选取

商用有限元软件ABAQUS中各求解器为非线性力学问题的数值计算工具，其主要作用是求解有限元模型在给定载荷和边界条件下的非线性力学行为。针对复合材料薄壁结构弯曲折叠数值模拟问题，主要有隐式和显式两种求解方法。ABAQUS/Standard分析模块采用隐式求解方法对线性、非线性平衡方程进行求解，包含Static General静力通用求解器和Dynamic Implicit隐式动力学求解器等。Static General静力通用求解器的核心思想是使用Newton-Raphson迭代来不断接近结构的平衡解，但该求解器在每一次迭代中必须获得切线刚度矩阵并对整个平衡方程完成求解，且切线刚度矩阵有时难以用公式形式表达，需采用数值计算方法进行求解，导致计算成本较高，不适合用于复杂的大型有限元计算。Dynamic Implicit隐式动力学求解器中使用的算法主要为加入可调数值阻尼的广义纽马克法（Newmark），即HHT法。此方法可以通过减小分析步大小，同时引入能量耗散项等方式解决收敛困难等问题。

隐式求解方法每完成一次迭代的计算成本非常高，且面对大变形等强非线性问题容易出现收敛困难等问题。为此，显式动力学Dynamic Explicit分析模块采用中心差分法对运动方程进行显式时间积分，将一个增量步的动力学分析结果作为计算下一个增量步的初始动力学条件。增量步在初始时提供了满足动力学平衡条件的加速度，并在时间上“显式地”前推速度和位移。显式求解方法最显著的特点是不涉及整体切向刚度矩阵求解，其为前推模型状态，无需迭代和收敛准则，可以较好地模拟一些复杂接触和极度不连续的情况。但该求解方法的时间增量步必须足够小，以保证数值模拟结果准确性，如图5所示。

2.3 计算结果

针对本计算模型，分别使用3种求解器计算得到复合材料薄壁带簧在弯曲折叠过程中单位长度力矩M_x对应的带簧结构中心位置处曲率{\kappa }_x、{\kappa }_y、{\kappa }_{xy}的值如图6（a）所示。

将由Dynamic Explicit显式动力学求解方法计算得到的结果与实验测试结果^[32]进行对比，结果如图6（b）所示。从图6（b）可以看出，3种求解器对复合材料薄壁带簧的弯曲折叠数值模拟均能够得到较好的有限元计算结果。进一步分析，基于隐式动力学的求解方法虽能成功处理复合材料薄壁带簧大变形下的数值计算问题，但其求解过程需要多次迭代计算，整个计算周期较长。基于显式动力学的求解方法避免了反复迭代计算，但存在计算过程中为保证求解精度所用时间增量步较小，从而导致计算成本相对较大等不足。然而，综合来看，显式动力学求解方法从计算原理上来看，更适用于大变形、复杂接触等强非线性问题求解。

3 结 论

针对复合材料层合板结构，本文开发了一种基于有限元方法的复合材料层合板ABD刚度矩阵快速获取方法。此方法可服务于工程应用，实现任意铺层顺序、铺层角度下复合材料ABD刚度矩阵各元素数值的快速准确获取，且与基于经典层合板理论计算结果一致。在此基础上，使用商用有限元软件ABAQUS中3种有限元求解器分别进行复合材料薄壁带簧大变形弯曲折叠数值计算，对比了隐式求解方法与显式求解方法的计算结果及计算精度。结果表明，Static General，Dynamic Implicit两种隐式求解方法和Dynamic Explicit显式求解方法均能准确实现复合材料薄壁结构大变形折叠数值模拟，然而Dynamic Explicit显式求解方法由于无需考虑迭代和收敛准则，避免了大变形等非线性因素下计算收敛困难等问题，对复合材料薄壁带簧大变形弯曲折叠等过程的数值模拟更具优势。

参考文献

[1] Dixit A, Mali H S.Modeling techniques for predicting the mechanical properties of woven-fabric textile composites:A Review[J].Mechanics of Composite Materials,2013,49(1):1‒20. doi:10.1007/s11029-013-9316-8

[2] Ma Xiaofei, Li Tuanjie, Ma Jingya,et al.Recent advances in space-deployable structures in China[J].Engineering,2022,17:207‒219. doi:10.1016/j.eng.2022.04.013

[3] Liu Tianwei, Bai Jiangbo, Fantuzzi N,et al.Thin-walled deployable composite structures:A review[J].Progress in Aerospace Sciences,2024,146:100985. doi:10.1016/j.paerosci.2024.100985

[4] Yang Chao, Wang Bing, Zhong Shuncong,et al.On tailoring deployable mechanism of a bistable composite tape-spring structure[J].Composites Communications,2022,32:101171. doi:10.1016/j.coco.2022.101171

[5] Wang Bing, Seffen K A, Guest S D.Shape of a bistable composite tape-spring in folding[C]//Proceedings of the AIAA Scitech 2019 Forum.San Diego:AIAA,2019:AIAA2019‒2276. doi:10.2514/6.2019-2276

[6] Fernandez J M, Viquerat A, Lappas V J,et al.Bistable over the whole length (BOWL) CFRP booms for solar sails[M]//Advances in Solar Sailing.Berlin:Springer,2014:609‒628. doi:10.1007/978-3-642-34907-2_38

[7] Fernandez J M, Visagie L, Schenk M,et al.Design and development of a gossamer sail system for deorbiting in low earth orbit[J].Acta Astronautica,2014,103:204‒225. doi:10.1016/j.actaastro.2014.06.018

[8] Lim J H, Kim Y B, Jang I S,et al.Parabolic deployable mesh antenna with a hingeless system of superelastic SMA ribs and composite tape springs[J].Acta Astronautica,2022,200:149‒162. doi:10.1016/j.actaastro.2022.08.007

[9] Wu Zhangming, Li Hao.Novel one-layer-cable-net antenna with high surface accuracy[J].Journal of Spacecraft and Rockets,2022,59(4):1179‒1192. doi:10.2514/1.a35225

[10] Duan Baoyan.Large spaceborne deployable antennas (LSDAs)—A comprehensive summary[J].Chinese Journal of Electronics,2020,29(1):1‒15. doi:10.1049/cje.2019.09.001

[11] Hilgemann E, Rajguru A, Sauder J F,et al.Development of a 6-meter crossed dipole antenna deployment mechanism for CubeSat applications[C]//Proceedings of the AIAA SCITECH 2024 Forum.Orlando:AIAA,2024:AIAA2024‒1430. doi:10.2514/6.2024-1430

[12] Fischer G, Panchenko M, Macher W,et al.Calibration of the JUICE RWI antennas by numerical simulation[J].Radio Science,2021,56(11):e2021RS007309. doi:10.1029/2021rs007309

[13] Mobrem M, Adams D S.Deployment analysis of the lenticular jointed antennas onboard the Mars express spacecraft[J].Journal of Spacecraft and Rockets,2009,46(2):394‒402. doi:10.2514/1.36890

[14] Ishikawa T, Chou T W.Elastic behavior of woven hybrid composites[J].Journal of Composite Materials,1982,16(1):2‒19. doi:10.1177/002199838201600101

[15] Ishikawa T, Chou T W.Stiffness and strength behaviour of woven fabric composites[J].Journal of Materials Science,1982,17(11):3211‒3220. doi:10.1007/bf01203485

[16] Angioni S L, Meo M, Foreman A.A critical review of homogenization methods for 2D woven composites[J].Journal of Reinforced Plastics and Composites,2011,30(22):1895‒1906. doi:10.1177/0731684411408158

[17] Gommers B, Verpoest I, Van Houtte P.The Mori-Tanaka method applied to textile composite materials[J].Acta Materialia,1998,46(6):2223‒2235. doi:10.1016/s1359-6454(97)00296-6

[18] An Ning, Jia Qilong, Jin Hao,et al.Multiscale modeling of viscoelastic behavior of unidirectional composite laminates and deployable structures[J].Materials & Design,2022,219:110754. doi:10.1016/j.matdes.2022.110754

[19] Omairey S L, Dunning P D, Sriramula S.Development of an ABAQUS plugin tool for periodic RVE homogenisation[J].Engineering with Computers,2019,35(2):567‒577. doi:10.1007/s00366-018-0616-4

[20] Daniel I M, Ishai O, Daniel I M,et al.Engineering mechanics of composite materials[M].New York:Oxford university press,1994.

[21] Quek S C, Waas A M, Shahwan K W,et al.Analysis of 2D triaxial flat braided textile composites[J].International Journal of Mechanical Sciences,2003,45(6/7):1077‒1096. doi:10.1016/j.ijmecsci.2003.09.003

[22] Mallikarachchi H M Y C, Pellegrino S.Quasi-static folding and deployment of ultrathin composite tape-spring hinges[J].Journal of Spacecraft and Rockets,2011,48(1):187‒198. doi:10.2514/1.47321

[23] Mallikarachchi H M Y C, Pellegrino S.Deployment dynamics of ultrathin composite booms with tape-spring hinges[J].Journal of Spacecraft and Rockets,2014,51(2):604‒613. doi:10.2514/1.a32401

[24] Yan Zhongxi, Wu Minger.Analysis on the peak moment of cylindrical tape spring under pure bending[J].Engineering Mechanics,2024,41(9):225‒233.

[25] 闫中曦,吴明儿.圆弧形带状弹簧纯弯曲峰值力矩分析[J].工程力学,2024,41(9):225‒233.

[26] Leclerc C, Wilson L L, Bessa M A,et al.Characterization of ultra-thin composite triangular rollable and collapsible booms[C]//Proceedings of the 4th AIAA Spacecraft Structures Conference.Grapevine:AIAA,2017:AIAA2017‒0172. doi:10.2514/6.2017-0172

[27] Leclerc C, Pellegrino S.Nonlinear elastic buckling of ultra-thin coilable booms[J].International Journal of Solids and Structures,2020,203:46‒56. doi:10.1016/j.ijsolstr.2020.06.042

[28] Seffen K A, Wang B, Guest S D.Folded orthotropic tape-springs[J].Journal of the Mechanics and Physics of Solids,2019,123:138‒148. doi:10.1016/j.jmps.2018.09.017

[29] Ye Hongling, Zhang Yang, Yang Qingsheng,et al.Quasi-static analysis and multi-objective optimization for tape spring hinge[J].Structural and Multidisciplinary Optimization,2019,60(6):2417‒2430. doi:10.1007/s00158-019-02331-6

[30] Ye Hongling, Li Bowen, Shi Xuesong, Zhang Yang.Quasi-static folding mechanical behavior analysis and optimization design for composite tube hinge[J].Journal of Mechanical Engineering,2020,56(5):172‒180.

[31] 叶红玲,李博文,史雪松,张洋.复合材料可展管铰链准静态折叠力学行为分析与优化设计[J].机械工程学报,2020,56(5):172‒180.

[32] Jin Hao, Jia Qilong, An Ning,et al.Surrogate modeling accelerated shape optimization of deployable composite tape-spring hinges[J].AIAA Journal,2022,60(10):5942‒5953. doi:10.2514/1.j061668

[33] Jin Hao, An Ning, Jia Qilong,et al.Optimizing deployment dynamics of composite tape-spring hinges[J].Thin-Walled Structures,2024,198:111738. doi:10.1016/j.tws.2024.111738

[34] Iqbal K, Pellegrino S.Bi-stable composite shells[C]//Proceedings of the 41st Structures,Structural Dynamics,and Materials Conference and Exhibit.Atlanta:AIAA,2000:AIAA2000‒1385. doi:10.2514/6.2000-1385

基金资助

空间先进机构与智能飞行器教育部重点实验室开放课题（序号3）