材料非线性对复合材料柱壳双稳态特性的影响

Impact of Material Nonlinearity on the Bistable Characteristics of Composite Cylindrical Shell

摘要

双稳态复合材料柱壳是一种新型可展开结构，已广泛应用于可折叠机翼、能量收集器及自适应结构等.当双稳态结构用于复杂工作环境时，结构材料属性将发生性能变化，进而影响结构的双稳态特性.结合理论研究和数值研究，分析材料属性对复合材料柱壳双稳态特性的影响.建立反对称铺设复合材料柱壳理论模型，导出结构变形过程中应变能的解析表达式，并分析材料的纵向弹性模量、横向弹性模量、剪切模量和泊松比等特征常数对柱壳结构应变能、主曲率和扭曲率的影响.结果表明，材料的纵向弹性模量和横向弹性模量的变化对双稳态柱壳结构的单位面积应变能和第二稳态主曲率均有明显影响，另外剪切模量G₁₂从10 GPa到2 GPa，单位面积应变能减少约72.98%，泊松比对第二稳态性能几乎没有影响.

Abstract

As a unique type of deployable structure， bistable composite material structures have extensive applications in various fields such as foldable wings， energy harvesters， and adaptive structures. When these bistable structures are employed in complex environments， the changes in material properties have a significant impact on their bistable characteristics. By combining theoretical and numerical investigations， this study examines the impact of material properties on the bistability of composite cylindrical shells. A theoretical model is developed for antisymmetrically laminated composite cylindrical shells， and analytical expressions for the strain energy of the shell during deformation are derived. Moreover， the effects of characteristic constants including longitudinal modulus of elasticity， transverse modulus of elasticity， shear modulus， and Poisson’s ratio on strain energy， principal curvatures， and torsional rates are analyzed for the shell structure. The results indicate that variations in the longitudinal and transverse modulus of elasticity of the material significantly affect the strain energy per unit area and the principal curvature of the second stable state in bistable cylindrical shells. Specifically， when the shear modulus G₁₂ decreases from 10 GPa to 2 GPa， the strain energy per unit area is reduced by approximately 72.98%.The Poisson’s ratio has almost no impact on the performance of the second stable state.

Graphical abstract

关键词

双稳态 / 复合材料 / 柱壳 / 材料属性 / 曲率

新型双稳态复合材料柱壳结构具有质轻、高强、空间利用率高等优点，在航空航天、能量收集等领域具有广泛应用前景^［1-3］.

双稳态柱壳结构具有伸展和卷拢两种几何构型的稳定状态，在外荷载、压电片、温度场以及记忆合金等激励下，结构可在两种稳定状态间相互跳变，且不需要持续的能量输入即可保持在某个特定稳定状态^［4-5］.通常柱壳结构在稳定状态改变时会处于不稳定的高应变能状态，此时需要外力或约束装置限制结构释放应变能^［6］. 但是约束装置增加了结构的复杂度. 双稳态柱壳结构在两种形态下均能保持自稳定，利用这种特性可有效简化设备的设计和建造.

1996年，Daton-Lovett^［7］研究发现反对称铺设复合材料层合柱壳具有双稳态特性. 依据简化模型研究柱壳的稳定性，结果表明反对称复合材料柱壳结构具有双稳态或多稳态特征，即结构总势能有两个或多个局部最小值. Iqbal等^［8-9］基于经典层压理论提出了双稳态复合材料柱壳的简单线弹性理论模型，给出了壳体总应变能随横向曲率、纵向曲率和铺设角变化的解析表达式. Guest和Pellegrino^［10］考虑扭曲变形的影响，基于Iqbal理论建立双稳态壳的本构方程，并提出双参数模型. 依据该模型，可计算得到柱壳构型变化时应变能的解析表达式，并可获得柱壳的稳定状态判别准则.

树脂基碳纤维复合材料易受环境温度和湿度的影响，材料特性会发生变化，进而导致结构的双稳态特性发生改变. Hyer^［11］研究了碳纤维T300/5208环氧树脂复合材料的弹性模量、泊松比和热膨胀系数等随温度的变化规律，并通过数据拟合得出了这些材料参数与温度之间的非线性关系式. Barker等^［12］研究了100~170 ℃范围内碳纤维复合材料的弹性模量、剪切模量和泊松比. 结果表明，树脂较其他复合材料组分更易受温度影响. 潘豪等^［13］研究了温度和湿度影响下反对称铺设复合材料圆柱壳的双稳态构型，预测了圆柱壳结构的第一稳态扭曲率、第二稳态主曲率和扭曲率. Chillara等^［14］研究了温度对纤维增强非对称铺设复合材料层压板稳态构型的影响. Moore等^［15］应用有限元法研究了非对称正交铺设复合材料柱壳力学性能的温度依赖性，分析表明热膨胀系数对结构的固化形状有重要影响. Zhang等^［16］采用试验和数值模拟相结合的方法研究了水分对复合材料柱壳双稳态特性的影响，并应用载荷-位移曲线以及跳变载荷表示柱壳在两个稳定状态之间的跳变过程.

研究反对称铺设复合材料柱壳结构的双稳态特性时通常仅考虑常温条件，但在复杂环境下材料参数将发生变化，进而对结构的稳态构型和稳态跳变产生影响.本文首先分析了反对称铺设复合材料柱壳的本构方程，推导结构形变过程中应变能的解析表达式，并结合典型算例探讨了材料的纵向弹性模量、横向弹性模量、剪切模量以及泊松比等对柱壳结构双稳态特性的影响. 应用数值方法模拟了反对称铺设复合材料柱壳的双稳态行为，得到柱壳在第二稳态时应变能和主曲率的表达式.

1 理论预测

1.1 层合壳的广义本构方程

复合材料柱壳结构通常在较高温度下经过固化制备得到，故结构的使用温度需低于固化温度.另外，树脂基体较增强纤维更容易受湿度的影响.因此，复合材料柱壳结构对湿热环境非常敏感，需重点探讨.

当复合材料结构仅受温湿作用时，单层柱壳的湿热应变为：

\left[\begin{array}{c}{\epsilon }_x^{\mathrm{e}}\\ {\epsilon }_y^{\mathrm{e}}\\ {\gamma }_{xy}^{\mathrm{e}}\end{array}\right]=\left[\begin{array}{cc}{\alpha }_x& {\beta }_x\\ {\alpha }_y& {\beta }_y\\ {\alpha }_{xy}& {\beta }_{xy}\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }T\\ C\end{array}\right]（1）

式中：ΔT为温度变化量；C为湿度.

单层壳的应力表达式为：

\left[\begin{array}{c}{\sigma }_x^{\mathrm{e}}\\ {\sigma }_y^{\mathrm{e}}\\ {\tau }_{xy}^{\mathrm{e}}\end{array}\right]=\mathit{Q}\left[\begin{array}{cc}{\alpha }_x& {\beta }_x\\ {\alpha }_y& {\beta }_y\\ {\alpha }_{xy}& {\beta }_{xy}\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }T\\ C\end{array}\right]（2）

式中：{\sigma }_x^{\mathrm{e}}、{\sigma }_y^{\mathrm{e}}、{\tau }_{xy}^{\mathrm{e}}分别为单层壳的x、y方向的正应力和x-y平面上的切应力；{\alpha }_x、{\alpha }_y、{\alpha }_{xy}和{\beta }_x、{\beta }_y、{\beta }_{xy}分别为x、y和x-y平面的热膨胀系数和湿膨胀系数，且在材料1、2主方向，温度和湿度变化不引起剪应变，因此{\alpha }_{\mathrm{12}}和{\beta }_{\mathrm{12}}为0；{\gamma }_{\mathrm{12}}^{\mathrm{e}}=\mathrm{0}； Q 为折减刚度矩阵.

\mathit{Q}=\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]（3）

将单层柱壳的应力-应变方程沿柱壳厚度方法积分，则复合材料柱壳的广义本构方程为：

\left[\begin{array}{c}{N}^{\mathrm{e}}\\ {\mathit{M}}^{\mathrm{e}}\end{array}\right]=\left[\begin{array}{cc}\mathit{A}& \mathit{B}\\ \mathit{B}& \mathit{D}\end{array}\right]\left[\begin{array}{c}{\epsilon }^{\mathrm{e}}\\ k^{\mathrm{e}}\end{array}\right]（4）

式中：N^e为温湿度场梯度所产生的等效力；M^e为温湿度场产生的力矩；{\epsilon }_{}^{\mathrm{e}}为湿热荷载引起的应变；k_{}^{\mathrm{e}}为湿热载荷引起的曲率； A、B、D 分别为柱壳的拉伸刚度矩阵、耦合刚度矩阵、弯曲刚度矩阵. 反对称复合材料柱壳的内力应变关系为：

\left[\begin{array}{c}{N}_x\\ N_y\\ N_{xy}\\ M_x\\ M_y\\ M_{xy}\end{array}\right]=\left[\begin{array}{cccccc}{A}_{\mathrm{11}}& A_{\mathrm{12}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& B_{\mathrm{16}}\\ A_{\mathrm{21}}& A_{\mathrm{22}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& B_{\mathrm{26}}\\ \mathrm{0}& \mathrm{0}& A_{\mathrm{66}}& B_{\mathrm{61}}& B_{\mathrm{62}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& B_{\mathrm{16}}& D_{\mathrm{11}}& D_{\mathrm{12}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& B_{\mathrm{26}}& D_{\mathrm{21}}& D_{\mathrm{22}}& \mathrm{0}\\ B_{\mathrm{61}}& B_{\mathrm{62}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& D_{\mathrm{66}}\end{array}\right]\left[\begin{array}{c}{\epsilon }_x^{\mathrm{0}}\\ {\epsilon }_y^{\mathrm{0}}\\ {\gamma }_{xy}^{\mathrm{0}}\\ k_x\\ k_y\\ k_{xy}\end{array}\right]（5）

A₁₆=A₂₆=0和D₁₆=D₂₆=0，表明反对称层合薄柱壳无拉剪耦合效应和弯扭耦合效应；B₁₆≠0和B₂₆≠0，表明柱壳存在拉扭耦合效应，但变形相对较小，且拉扭耦合效应会随着铺层数的增加而逐步减小；B₁₁=B₁₂=B₂₂=0和B₆₆=0，表明柱壳没有拉弯耦合效应和剪扭耦合效应.

柱壳的应变和曲率为：

\left[\begin{array}{c}{\epsilon }^{\mathrm{e}}\\ k^{\mathrm{e}}\end{array}\right]={\left[\begin{array}{cc}\mathit{A}& \mathit{B}\\ \mathit{B}& \mathit{D}\end{array}\right]}^{-\mathrm{1}}\left[\begin{array}{c}{\mathit{N}}^{\mathrm{e}}\\ {\mathit{M}}^{\mathrm{e}}\end{array}\right]（6）

可知，层合柱壳的吸湿变形与热变形在形式上基本相同，热膨胀系数等效于湿膨胀系数，温差等效于吸水浓度. 在湿热环境下，柱壳的本构方程包括外载荷和湿热环境影响产生的应变和曲率，其表达式为：

\left[\begin{array}{c}\mathit{N}\\ \mathit{M}\end{array}\right]=\left[\begin{array}{c}{\mathit{N}}^{\mathrm{M}}\\ {\mathit{M}}^{\mathrm{M}}\end{array}\right]+\left[\begin{array}{c}{\mathit{N}}^{\mathrm{e}}\\ {\mathit{M}}^{\mathrm{e}}\end{array}\right]=\left[\begin{array}{cc}\mathit{A}& \mathit{B}\\ \mathit{B}& \mathit{D}\end{array}\right]\left[\begin{array}{c}\epsilon \\ k\end{array}\right]=

\left[\begin{array}{cc}\mathit{A}& \mathit{B}\\ \mathit{B}& \mathit{D}\end{array}\right]\left(\left[\begin{array}{c}{\epsilon }^{\mathrm{M}}\\ k^{\mathrm{M}}\end{array}\right]+\left[\begin{array}{c}{\epsilon }^{\mathrm{e}}\\ k^{\mathrm{e}}\end{array}\right]\right)（7）

中性面的应变和曲率为：

\left[\begin{array}{c}\epsilon \\ k\end{array}\right]={\left[\begin{array}{cc}\mathit{A}& \mathit{B}\\ \mathit{B}& \mathit{D}\end{array}\right]}^{-\mathrm{1}}\left[\begin{array}{c}\mathit{N}\\ \mathit{M}\end{array}\right]（8）

1.2 薄柱壳结构的应变能

湿热作用会引起膨胀变形，影响柱壳结构的拉伸应变能. 依据最小势能原理，可得到反对称铺设层合柱壳结构的应变和曲率. 柱壳结构的弯曲应变能u_{\mathrm{b}}和拉伸应变能u_{\mathrm{s}}分别为：

u_{\mathrm{b}}=\frac{\mathrm{1}}{\mathrm{2}}\left(M_x{k}_x+M_y{k}_y+M_{xy}{k}_{xy}\right)（9）

u_{\mathrm{s}}=\frac{\mathrm{1}}{\mathrm{2}}\left(N_x{\epsilon }_x+N_y{\epsilon }_y+N_{xy}{\gamma }_{xy}\right)（10）

即

u_{\mathrm{b}}=\frac{k_x}{\mathrm{2}}(B_{\mathrm{16}}{\gamma }_{xy}+D_{\mathrm{11}}{k}_x^{}+\mathrm{2}{D}_{\mathrm{12}}{k}_y)+

\frac{k_y}{\mathrm{2}}(D_{\mathrm{22}}{k}_y^{}+B_{\mathrm{26}}{\gamma }_{xy})+\frac{k_{xy}}{\mathrm{2}}(B_{\mathrm{16}}{\epsilon }_x+B_{\mathrm{26}}{\epsilon }_y+D_{\mathrm{66}}{k}_{xy}^{})（11）

u_{\mathrm{s}}=\frac{{\epsilon }_x}{\mathrm{2}}(A_{\mathrm{11}}{\epsilon }_x^{}+\mathrm{2}{A}_{\mathrm{12}}{\epsilon }_y+B_{\mathrm{16}}{k}_{xy})+

\frac{{\epsilon }_y}{\mathrm{2}}(A_{\mathrm{22}}{\epsilon }_y^{}+B_{\mathrm{26}}{k}_{xy})+\frac{{\gamma }_{xy}}{\mathrm{2}}(A_{\mathrm{66}}{\gamma }_{xy}^{}+B_{\mathrm{16}}{k}_x+B_{\mathrm{26}}{k}_y)（12）

式中：{\gamma }_{xy}={\gamma }_{xy}^{\mathrm{0}}+z{k}_{xy}+{\gamma }_{xy}^{\mathrm{e}}.

令{\mathit{A}}^{\mathrm{\text{'}}}={\mathit{A}}^{-\mathrm{1}}和{\mathit{B}}^{\mathrm{\text{'}}}=-{\mathit{A}}^{-\mathrm{1}}\mathit{B}，则柱壳结构在第二稳态时的扭率为：

{\gamma }_{xy}=B_{\mathrm{16}}^{\mathrm{\text{'}}}{k}_x+B_{\mathrm{26}}^{\mathrm{\text{'}}}{k}_y+{\gamma }_{xy}^{\mathrm{e}}（13）

将式（13）代入式（11）和式（12），可得总的单位面积应变能为：

u=u_{\mathrm{b}}+u_{\mathrm{s}}=

\frac{\mathrm{1}}{\mathrm{2}}[(B_{\mathrm{16}}{B}_{\mathrm{16}}^{\mathrm{\text{'}}}+D_{\mathrm{11}})k_x^{\mathrm{2}}+(B_{\mathrm{26}}{B}_{\mathrm{26}}^{\mathrm{\text{'}}}+D_{\mathrm{22}})k_y^{\mathrm{2}}]+

\frac{k_{xy}}{\mathrm{2}}(B_{\mathrm{16}}{\epsilon }_x+B_{\mathrm{26}}{\epsilon }_y+\frac{D_{\mathrm{66}}}{\mathrm{2}}{k}_{xy}^{\mathrm{2}})+

\mathrm{2}{k}_x{k}_y(B_{\mathrm{16}}{B}_{\mathrm{26}}^{\mathrm{\text{'}}}+D_{\mathrm{12}})+\frac{{\epsilon }_x}{\mathrm{2}}(A_{\mathrm{11}}{\epsilon }_x^{}+\mathrm{2}{A}_{\mathrm{12}}{\epsilon }_y)+

\frac{\mathrm{1}}{\mathrm{2}}[A_{\mathrm{22}}{\epsilon }_y^{\mathrm{2}}+A_{\mathrm{66}}({\gamma }_{xy}^{\mathrm{e}}{)}^{\mathrm{2}}]（14）

k_x2、k_y2和k_xy2分别为柱壳结构第二稳态的x方向、y方向曲率与扭曲率，湿热环境下柱壳结构的应变和曲率为：

\left\{\begin{array}{c}{\epsilon }_x={\epsilon }_x^{\mathrm{0}}+{\epsilon }_x^{\mathrm{e}}\\ {\epsilon }_y={\epsilon }_y^{\mathrm{0}}+{\epsilon }_y^{\mathrm{e}}\end{array}\right.{,}_{}\left\{\begin{array}{l}{k}_x=k_{x\mathrm{2}}-k_x^{\mathrm{e}}\\ k_y=k_{y\mathrm{2}}-\mathrm{1}/R-k_y^{\mathrm{e}}\\ k_{xy}=k_{xy\mathrm{2}}-k_{xy}^{\mathrm{e}}\end{array}\right.（15）

柱壳变形时不考虑横向应变，即{\epsilon }_y^{\mathrm{0}}\approx \mathrm{0}.当柱壳结构处于第二稳定状态时应变能极小，且此时横向曲率k_y2≈0. 当应变能满足du/dk_x2=0和du/dk_xy2=0时，可得到层合柱壳第二稳态的曲率.

\left\{\begin{array}{l}{k}_{x\mathrm{2}}=k_x^{\mathrm{e}}+\frac{B_{\mathrm{16}}{B}_{\mathrm{62}}^{\mathrm{\text{'}}}+B_{\mathrm{26}}{B}_{\mathrm{61}}^{\mathrm{\text{'}}}+\mathrm{2}{D}_{\mathrm{12}}}{\mathrm{2}\left(B_{\mathrm{16}}{B}_{\mathrm{61}}^{\mathrm{\text{'}}}+D_{\mathrm{11}}\right)}\left(\frac{\mathrm{1}}{R}+k_y^{\mathrm{e}}\right)\\ k_{xy\mathrm{2}}=k_{xy}^{\mathrm{e}}-\frac{B_{\mathrm{16}}{\epsilon }_x^{\mathrm{e}}+B_{\mathrm{26}}{\epsilon }_y^{\mathrm{e}}}{D_{\mathrm{66}}}\end{array}\right.（16）

1.3 材料参数分析

采用理论分析与数值仿真相结合的方法研究材料属性对复合材料柱壳双稳态特性的影响.由于壳体的材料性能受温度、湿度等因素影响，针对材料的4个独立参数（纵向弹性模量、横向弹性模量、剪切模量和泊松比），研究各材料参数对结构双稳态特性的影响. 选用碳纤维环氧树脂（T700/Epoxy）复合材料，铺设方式为5层45°反对称铺设，结构的几何参数以及常温下材料参数分别如表1和表2所示. 文中：L为柱壳的长度，mm；t为单层壳厚度，mm；ν₁₂为泊松比，不受温湿度影响；u_a为柱壳结构的单位面积应变能，N·mm/mm²；k_x2为主曲率，mm^-1；T为温度，℃.根据Hyer^［11］研究所得碳纤维环氧树脂材料参数与温度的关系式，通过最小二乘法可导出材料参数随环境变化的非线性关系式：

\begin{array}{l}{E}_{\mathrm{1}}=\mathrm{108.1}+\mathrm{0.023}\mathrm{5}T-\mathrm{0.001}\mathrm{62}{T}^{\mathrm{2}}\\ E_{\mathrm{2}}=\mathrm{32.73}{\mathrm{e}}^{-\mathrm{0.111}T}+\mathrm{3.746}{\mathrm{e}}^{-\mathrm{0.003}\mathrm{219}T}\\ G_{\mathrm{12}}=\mathrm{21.53}-\mathrm{12.43}{T}^{\mathrm{0.091}\mathrm{88}}\\ {\alpha }_{\mathrm{1}}=(-\mathrm{2.2}+\mathrm{0.002}\mathrm{5}T+\mathrm{0.000}\mathrm{125}{T}^{\mathrm{2}}){\mathrm{10}}^{-\mathrm{6}}/℃\\ {\alpha }_{\mathrm{2}}=(\mathrm{68.283}\mathrm{3}-\mathrm{0.197}\mathrm{5}T){\mathrm{10}}^{-\mathrm{6}}/℃\\ {\beta }_{\mathrm{2}}=\mathrm{0.005}{W}^{-\mathrm{1}}\end{array}

其中，W表示含水量的质量百分比.材料的纵向弹性模量E₁、横向弹性模量E₂、剪切模量G₁₂的材料参数随温度变化的非线性关系曲线分别如图1、图2和 图3所示.

1.3.1 纵向弹性模量

依次取复合材料的纵向弹性模量E₁为108、103、98、93、88 GPa，不改变其他材料参数，则单位面积应变能u_a随θ的变化情况如图4所示.

由图4可知，柱壳结构在θ=π/2时应变能极小. 通过计算得到反对称铺设双稳态柱壳在第二稳态时的单位面积应变能u_a2和主曲率k_x2，如表3所示.

纵向弹性模量E₁从108 GPa降低到88 GPa，单位面积应变能u_a2减少约3.43%，第二稳态主曲率k_x2减少约6.74%，反对称铺设层合柱壳的第二稳态应变能与主曲率逐级递减. 主要原因为：材料的弹性模量越高，结构抵抗弹性变形的能力越强，即刚度越大，双稳态结构的第二稳态卷曲半径越小，卷拢效果越明显.

1.3.2 横向弹性模量

依次取复合材料的横向弹性模量E₂为10、8、6、4、2 GPa，不改变其他材料参数，则单位面积应变能u_a随θ的变化情况，如图5所示.

由图5可知，柱壳结构在θ=π/2时应变能极小，对应柱壳的第二稳态，此时单位面积应变能u_a2和主曲率k_x2如表4所示.

横向弹性模量E₂从10 GPa降低到2 GPa，单位面积应变能u_a2减少约2.07%，第二稳态主曲率k_x2减少约4.44%，柱壳的第二稳态应变能与主曲率减少幅度较纵向弹性模量E₁小. 主要原因为：材料的横向弹性模量E₂远小于纵向弹性模量E₁，取值范围更小，虽与E₁变化趋势相同，但变化幅度较小.

1.3.3 剪切模量

依次取复合材料的剪切模量G₁₂为10、8、6、4、2 GPa，不改变其他材料参数，则第二稳态时单位面积应变能u_a随θ的变化曲线，如图6所示.

由图6可知，柱壳结构在θ=π/2时应变能极小，计算可得柱壳在第二稳态时的单位面积应变能u_a2和主曲率k_x2，如表5所示.

剪切模量G₁₂从10 GPa降低到2 GPa，单位面积应变能u_a2减少约72.98%，第二稳态主曲率k_x2增加约73.85%，薄柱壳的第二稳态应变能逐渐减少，第二稳态主曲率逐级递增. 主要原因为：材料剪切模量越大，表示结构抵抗剪切变形的能力越强，即刚度越大，越不易发生剪切变形. 对于双稳态结构，即第二稳态卷曲半径越小，卷拢效果越明显.

1.3.4 泊松比

依次取复合材料的泊松比{\nu }_{\mathrm{12}}为0.45、0.40、0.35、0.30、0.25，不改变其他材料参数，则单位面积应变能u_a随θ的变化情况，如图7所示.

由图7可知，不同泊松比时，柱壳结构的应变能随θ变化的曲线很接近，且都在θ = π/2处应变能极小. 通过计算可得双稳态柱壳在第二稳态时的单位面积应变能u_a2和主曲率k_x2，如表6所示.

泊松比{\nu }_{\mathrm{12}}从0.45降低到0.25，单位面积应变能u_a2减少约0.42%，第二稳态主曲率k_x2减少约1.24%. 反对称铺设复合材料柱壳的第二稳态应变能和主曲率均线性减少，但减小幅度较小，即泊松比对复合材料柱壳的主曲率和应变能几乎没有影响.

2 数值模拟

数值模拟可得到结构各区域的应力和应变数据，并可模拟难以试验的复杂环境，且仿真过程和结果易于实现可视化. 本研究涉及几何非线性和材料非线性，ABAQUS擅长处理非线性问题，可以得到双稳态柱壳的曲率和扭曲率.

2.1 有限元模型

利用ABAQUS软件建立壳结构的有限元模型.首先根据壳体的几何参数定义壳体的截面，然后使用拉伸命令生成构件. 为了准确、快速地定义边界条件、荷载和预定义场，一般将壳结构分割为四个部分，如图8所示.

壳体几何模型部件建立完成后，通过属性模块赋予其材料属性，并定义薄壳结构的截面参数，将截面属性分配给壳模型.

在网格模块中对壳体进行网格划分. 采用S4R减缩积分壳单元作为单元类型. 网格尺寸设置为 3 mm，通过“布种”命令为壳体各边指定单元数，参数设置完成后生成的网格模型如图9所示.

为分析壳体结构的稳态跳变过程，设计3个位移驱动的分析步骤，并设置相应的边界条件与荷载. 在初始分析步中，对壳体中心点施加全约束，固定所有自由度，确保结构初始状态稳定［图10（a）］，该步骤全程激活. 在分析步1中，通过对壳体两条直边施加向下的位移荷载，使壳体展平［图10（b）］. 随后，在分析步2中，对壳体两条曲边的中点施加向上的位移荷载，诱导壳体卷曲［图10（c）］. 最后，在分析步3中移除所有的外荷载，仅保留边界条件约束，使壳体达到自平衡的第二稳定状态［图10（d）］. 每个分析步仅在当前步骤激活.

2.2 模拟结果分析

通过有限元模拟预测材料参数对反对称铺设复合材料柱壳单位面积应变能u_{\mathrm{f}}和主曲率k_{x\mathrm{2}}^{\mathrm{f}}的影响，并将理论结果与数值结果进行对比分析.

2.2.1 纵向弹性模量对比分析

依次取复合材料的纵向弹性模量E₁为108、103、98、93、88 GPa，其他参数保持不变，薄柱壳结构第二稳态的特征参数随纵向弹性模量E₁的变化曲线如 图11、图12所示.

由图11和图12可知，E₁在88~108 GPa范围内，复合材料柱壳在第二稳态时，单位面积应变能和主曲率的理论结果与有限元结果变化趋势相同. 单位面积应变能u的平均误差值约为5.27%，主曲率k_x2的平均误差值约为-7.28%，二者相差不大，吻合程度好.

2.2.2 横向弹性模量对比分析

依次取复合材料的横向弹性模量E₂为10、8、6、4、2 GPa，其他参数保持不变，柱壳结构的第二稳态特征参数随横向弹性模量E₂的变化曲线如图13、 图14所示.

由图13和图14可知，E₂在2~10 GPa范围内，复合材料柱壳在第二稳态时，单位面积应变能和主曲率的理论结果与有限元结果变化趋势基本相同，应变能u的平均误差值约为4.3%，主曲率k_x2的平均误差值约为-6.82%，二者吻合较好.

2.2.3 剪切模量对比分析

依次取复合材料的剪切模量G₁₂为10、8、6、4、 2 GPa，其他参数保持不变，柱壳第二稳态特征参数随剪切模量G₁₂的变化曲线，如图15、图16所示.

由图15和图16可知，G₁₂在2~10 GPa范围内，复合材料柱壳第二稳态单位面积应变能和主曲率的理论结果与有限元结果的变化趋势相同，均随着G₁₂的增大而增大，单位面积应变能u的平均误差值约为5.91%，主曲率k_x2的平均误差值约为-8.67%，二者吻合较好.

2.2.4 泊松比对比分析

依次取复合材料的泊松比ν₁₂为0.45、0.40、0.35、0.30、0.25，其他参数保持不变，ν₁₂柱壳结构第二稳态特征参数随泊松比ν₁₂的变化曲线如图17、图18所示.

由图17和图18可知，ν₁₂在0.25~0.45范围内，柱壳第二稳态单位面积应变能和主曲率的理论结果与有限元结果随泊松比的变化趋势相同，均随ν₁₂的增大而增大. 单位面积应变能u的平均误差值约为11.00%，主曲率k_x2的平均误差值约为-7.35%，二者吻合较好.

总体而言，数值结果与理论结果变化趋势基本相同，但数值略有差异，主要原因是理论分析时简化了运算过程，忽略了稳态跳变时高应变能对结构的影响.

3 结 论

复合材料柱壳结构在受温湿环境影响时，结构自身的材料属性会发生非线性变化. 本文建立了复合材料柱壳理论模型，导出结构变形过程中应变能的解析表达式，利用最小势能原理研究柱壳的第二稳态，并探讨材料的4个特征常数对反对称铺设复合材料层合柱壳第二稳态应变能、主曲率和扭曲率的影响. 主要结论如下：

1）纵向弹性模量E₁从108 GPa降低到88 GPa，单位面积应变能减少约3.43%，第二稳态主曲率减少约6.74%，表明材料刚度减小，卷拢效果减弱.

2）横向弹性模量E₂从10 GPa降低到2 GPa，单位面积应变能减少约2.07%，第二稳态主曲率减少约4.44%. 由于E₂远小于E₁，E₂对结构的影响相对较小，但二者变化趋势基本一致.

3）剪切模量G₁₂从10 GPa降低到2GPa，单位面积应变能减少约72.98%，第二稳态主曲率增加约73.85%，表明随着剪切模量的减小，结构更易发生变形，卷曲效果更明显.

4）泊松比{\nu }_{\mathrm{12}}从0.45降低到0.25，单位面积应变能减少约0.42%，第二稳态主曲率减少约1.24%. 泊松比的影响相对较小，对柱壳的第二稳态主曲率和单位面积应变能几乎没有影响.

综上，材料的纵向弹性模量和剪切模量对双稳态薄柱壳的双稳态性能影响显著，但横向弹性模量和泊松比的影响相对较小.

参考文献

[1] LU Z Q， SHAO D， FANG Z W，et al ．Integrated vibration isolation and energy harvesting via a bistable piezo-composite plate［J］．Journal of Vibration and Control，2019，26（9/10）：779-789．

[2] 陈炳彬，张征，鲁聪达，等 ．复合材料层合结构在防覆冰/除冰系统中的应用［J］．中国机械工程，2019，30（7）：771-776．

[3] CHEN B B， ZHANG Z， LU C D，et al ．Applications of composite laminated structures in anti-icing and de-icing systems［J］．China Mechanical Engineering，2019，30（7）：771-776．（in Chinese）

[4] MALLOL P， MAO H N， TIBERT G ．Experiments and simulations of the deployment of a bistable composite boom［J］．Journal of Spacecraft and Rockets，2018，55（2）：292-302．

[5] DANO M L， HYER M W. Thermally-induced deformation behavior of unsymmetric laminates［J］. International Journal of Solids and Structures， 1998， 35（17）： 2101-2120．

[6] PORTELA P， CAMANHO P， WEAVER P，et al ．Analysis of morphing，multi stable structures actuated by piezoelectric patches［J］．Computers & Structures， 2008， 86（3/4/5）：347-356．

[7] WU Y P. Study on bi-stable behaviors of isotropic shell structures and numerical simulation［J］. Advanced Materials Research，2010， 168/169/170： 341-344．

[8] DATON-LOVETT A. An extendible member： US6217975［P］. 2001-04-17.

[9] IQBAL K， PELLEGRINO S， DATON-LOVETT A. Bi-stable composite slit tubes［M］//IUTAM-IASS symposium on deployable structures：theory and applications. Dordrecht：Springer Netherlands， 2000： 153-162．

[10] IQBAL K， PELLEGRINO S ．Bi-stable composite shells［C］//41st Structures，Structural Dynamics，and Materials Conference and Exhibit. Atlanta， GA， USA： AIAA， 2000： AIAA2000-1385．

[11] GUEST S D， PELLEGRINO S. Analytical models for bistable cylindrical shells［J］. Proceedings of the Royal Society A：Mathematical，Physical and Engineering Sciences，2006，462（2067）： 839-854．

[12] HYER M W. The room-temperature shapes of four-layer unsymmetric cross-ply laminates［J］. Journal of Composite Materials， 1982， 16（4）： 318-340．

[13] BARKER A J， VANGERKO H. Temperature dependence of elastic constants of CFRP［J］. Composites， 1983， 14（1）： 52-56．

[14] 潘豪， 柴国钟， 张征，等 ．温湿环境反对称铺设圆柱壳结构的双稳态模型与变形调控研究［J］．机电工程，2017，34（11）：1235-1242．

[15] PAN H， CHAI G Z， ZHANG Z，et al. Hygrothermal effect on bistable model of anti-symmetric cylindrical shells and deformation control［J］. Journal of Mechanical & Electrical Engineering，2017， 34（11）： 1235-1242．（in Chinese）

[16] CHILLARA V S C， DAPINO M J. Mechanically-prestressed bistable composite laminates with weakly coupled equilibrium shapes［J］. Omposites Part B：Engineering， 2017， 111： 251-260．

[17] MOORE M， ZIAEI-RAD S， SALEHI H. Thermal response and stability characteristics of bistable composite laminates by considering temperature dependent material properties and resin layers［J］. Applied Composite Materials， 2013， 20（1）： 87-106．

[18] ZHANG Z， PAN H， WU H L，et al ．Hygroscopic influence on bistable characteristics of antisymmetric composite cylindrical shells：an experimental study［J］. Journal of Composite Materials，2018， 52（26）： 3565-3577．

基金资助

陕西省自然科学基础研究计划资助项目（2023-JC-YB-435）