复合材料飞轮的径向强度过低是制约飞轮转速提高的主要因素，在飞轮缠绕过程中对纤维束施加足够的张紧力是解决这一问题的有效的办法之。通过纤维束张紧力缠绕成型工艺，可使得成型的复合材料飞轮径向初应力处于数值足够大的压应力状态；其关键是建立纤维张归力与缠绕后飞轮的初应力之间的定量关系。Springer将纤维束连续缠绕简化为多层复合材料圆筒的外压叠加模型，但该模型实质上忽略了缠绕过程中芯轴的变形对均匀压力的影响，从而造成较大的计算误差。本文作者将纤维束张紧力缠绕看成多层复合材料薄环连续过盈装配的过程，并在此基础上建立了计算纤维张紧力与复合材料飞轮初应力的平面应力型全弹性轴对称模型。
然而实际的力学问题是三维的具有轴向效应；为了分析张紧力缠绕过程中的轴向效应，必须建立三维分析模型。本文中采用ANSYS软件对张紧力缠绕复合材料飞轮的初应力进行三维数值分析，以期望对前述建立的平面应力型全弹性模型作一比较和验证，并给出三给轴向效应的表征方法。
1 三维有限元分析模型
三维模型和平面应力模型的共同点在于其物理思想，都把张紧力缠绕简化为多层复合材料薄环连续过盈装配过程，且均采用了当前缠绕层的拉应力等于张紧力除以纤维层截面积的假设（σθ=F/A）；区别在于三维模型的控制方程、边界条件和连接条件均经平面应力模型高一维。
对于图1所示的三维轴对称和轴向对称过盈装配模型，图中复合芯轴的定义相同，包括金属芯轴（或模具芯轴）已缠绕的k-1层复合材料环；该模型基于柱坐标形式给出的轴对称的控制方程详见文献，其边界条件为

显然对于上术三维边值问题不能给出解析解，必须用数值分析的方法才能得到三维过盈装配的全场应力和变形分布。图2给出了缠绕层数为N，轴向长度为L的纤维束张紧力缠绕飞轮的整体模型的有限元网络，由于结构的轴对称和轴向对称，计算模型仅取轴向的1/2；当N=100，L=40mm时，模型由82987个轴对称单元和63719个节点组成。在金属芯轴和复合材料薄环的连接部分，采用了三角形单元与四边形单元过渡（如图3的局部放大所示），这样可以减小模型单元和节点的规模，提高计算效率。边界条件、装配界面处连接条件分别同式（1）～（3），三维过盈装配界面处的连接条件采用ANSYS的面-面接触算法实现，ANSYS的面-面接触算法是经过多步循环迭代使得主面（单元定义为CONTA171）和从面（单元定义为TARGE169）上对应节点的径向位移值之和等于过盈量（式（2）），同时满足径向应力平衡条件（式（3）），从而实现主从面对应节点的重合，完成过盈装配计算。

图4给出了利用该三维轴对称模型求解的具体过程示意图。在缠绕第1层（图4（a））时，金属芯轴和第1层复合材料薄环之间的地盈量为δ（1）利用而-面接触算法求解第1次过盈装配，得到一次装配初应力及变形，A点变形至A′点。当缠绕第2层（图4（b））时，在A′点外侧建立过盈量为δ（2）的复合材料薄环，将由金属芯轴和第1层薄环组成的复合芯轴看成初应力为零的初始状态，并利用面-面接触算法求解第2次过盈装配，得到二次装配初应力及变形，B点变形至B′点。依此类推，即可计算得到N次装配初应力及变形（如图4（c））；后利用叠加原理求得N次装配完成后的飞轮的总体初应力，具体的计算流程见图5。

由于在三维分析中不能预先给出地盈量与张紧力之间的解析关系，即不能由张紧力直接得到过盈理，只能由预先设定的过盈量来实施三维的初应力分析，故本文中，将平面应力模型计算得到的真实过盈量作为三维有限元模型的地盈量；后依初应力分析的结果并按
бθ=F/A反过来得到张紧力。
三维分析的有关材料参数和几何参数见表1。

2 计算结果与分析
针对表1的材料参数，先进行了4个算例分析。这4个算例的缠绕层数均为100层，轴向长度L分别为40、80、120、160mm；预先设定的过盈量如图6所示，它们是在张紧力恒为10N、缠绕层数为100层时，由平面应力模型得到的真实过盈量。

图7给出了对L=400mm的算例计算得到的纤维束张紧力缠绕复合材料飞轮（实心）的环向初应力以及径向初应力分布云图。从图7可以看出，环向初应力均处理拉应力状态，外层处的环向初应力大，且在轴向方向分布比较均匀，与缠绕过程中施加的恒定张紧力吻合；飞轮的径向初应力均处于压应力状态，且与金属芯轴连接的里层的初始径向压应力大。
 为了分析三维模型的轴向效应，在上述4个算例中分别定义如图8所示的3个截面：z=0(Section),z=L/4(Section2),z=L/2(Section3)。以下将给出这4个算例的在这3个截面上的初应力分布，还将给出相应的平面应力模型的解析结果。
2．1．1 不同截面的径向初应力分布
图8（a）～图8(d)分别给出了上述4个算例的3个截面上的径向初应力以及平面应力模型的结果。从图8中可以看出，三维模型的初始径向压应力（数值）均略小于平面应力模型的，且这种差距随着L增大则略有增大；截面1和截面2上径向初始压应力分布几乎重合，截面3的初始径向压应力分布与前两者相比有所差别，这种差别在外层十分微小，随r的减小而明显增加（初始径向压应力数值减小）。
2．1．2 不同截面的环向初应力分布
图9(a) ～图9(d)分别给出了上述4个算例的3个截面的环向初应力分布曲线以及平面模型的计算结果。从图9可以看出，三维模型计算的环向应力均略小于平面应力模型的，且这种差距随着L的增大而缓慢增友。截面1与截面2的环向应力分布几乎重合，面截面3的环向应力分布与前两者有一定偏差。

．1．3 三维轴向效应的表征方法
图10给出了轴向长度L=400mm的算例的外层复合材料薄环的环向初应力沿轴向分布曲线。从图10可以看出，个层的环向初应力沿z的分布曲线从中面（z=0）到端面（z=1/2）的大部分区域（以下称为内部区域）可近似看成与z轴的平行直线，仅是在靠近端面的边界区域有一定的跳动。由此内部区域环向初应力的平均值бθ（63.113MPa）并按бθ=F/A即可得到F=9.97N，这与用平面应力模型为计算过盈量δ（k）而设定的张紧力F=10N非常接近，相对误差仅为－0.3%。这说明平面应力模型给出的过盈量与张紧力之间的解析关系在三维情况下仍有很高精度，且本文作者提出的利用平面应力模型为三维计算提供设定过盈量，后再通过应力分析确定张紧力F的方法也是有效的。
图10所示的环向初应力在轴向边界区域的跳

动是由于三维模型的轴向效应引起的。为阐述的方便，将环向初应力跳动的区域称为轴向效应的影响区l，并用相对影响长度K=2*l/来表示该影响区的范围；用相对幅值H=（бθmax－бθ）来表示影响区бθ的跳动，其中，бθmax为图中大环向初应力值。
本文中还对不同的缠绕层数（1，10，50）以及不同的轴向长度（40,80,120,160mm）组合的12个算例进行了三维分析，其中设定的过盈量均由图6给出。对于本文中的所有算例，按照前述方法分别计算相对影响长度K以及相对幅值H等。
图11（a）给出了相对影响长度K与缠绕层数和轴向长度之间的关系。从图中可以看出，当缠绕层数一定，随着模型轴向长度的增大，相对影响长度逐渐变小；当模型轴向长度一定，随着缠绕层数的增大，相对影响长度逐渐变大。图11(b)给出的是相对幅值H与缠绕层数和轴向长度之间的关系，从图中可以看出，相对幅值随着轴向长度和缠绕层数的增大而缓慢增大，但绝对数值均较小。

3 结 论
（1） 利用平面应力全弹性模型的思想（即将纤维束张紧力缠绕看成多层复合材料薄环连续过盈装配的过程），建立了三维纤维张紧力缠绕复合材料飞轮初应力分析模型，并给出了基于面-面接触算法求解张紧力缠绕复合材料飞轮初应力的三维数值方法。算例分析表明，该三维数值模型和方法是合理有效的，同时利用平面应力模型为三维计算提供设定过盈量，后再通过应力分析确定张紧力F的方法也是有效的。该三维数值分析方法可直接用于多环过盈装配的初应力分析。
（2） 三维数值分析得到的环向初应力以及径向初始压应力（数值）均略低于平面应力模型的，且这种差距随着飞轮轴向长度增加而缓慢增大。
（3） 三维计算中的相对影响长度随缠绕飞轮的轴向长度增大而变小，随缠绕层数增加而变大；相对幅值随缠绕飞轮轴向长度及缠绕层数增大而缓慢增大，但绝对数值较小。
（4） 三维分析证实了平面应力模型关于复合材料飞轮张紧力缠绕的初应力分析有足够的精度，可作为复合材料飞轮初步设计的基础
（5） 在三维数值分析模型中，层内沿轴向方向设定的过盈量相同，但算出的张紧力在轴向边界区域有一定跳动；而实际缠绕过程中，一层缠绕的张紧力均相同，故该模型尚不能完全反应出同一层等张紧力的缠绕过程；这也是该模型尚需进一步改进的地方。