0 引　言

海洋立管作为连接海底与海面的桥梁，是海洋油气开采系统的重要组成部分。作业时立管附着在海洋平台上，但随着开采深度的增加，传统金属海洋立管会增加海洋平台的建造成本。为此，纤维增强复合材料凭借强度高、质量轻以及耐腐蚀等优点成为深水立管开发中替代金属材料的首要选择。但纤维复合材料的各向异性特点，导致结构设计复杂且变量众多。传统纤维复合海洋立管铺层优化设计以经典层合板理论为基础，假设复合材料由各层之间完全粘合的多个单层板叠合而成^[1]。计算时将每个片层视为均匀的正交异性材料，其应力应变关系遵循广义胡克定律。设计时，由设计者根据经验给出初始参数进行试算，再依据计算结果调整设计，如此往复得出最优结果。这种工作方式不但对设计者的经验要求较高，而且工作量庞大且效率较低。因此，借助智能算法开展纤维复合海洋立管铺层优化具有重要意义。

目前，国内外已有许多学者借助智能算法完成海洋立管优化研究。由于传统基于导数的优化方法在求解高维问题时容易陷入局部最优解，于是，适合求解多变量优化问题的元启发式算法受到研究学者的青睐。余杨等^[2]应用NSGA-Ⅱ算法对缓坡型立管截面、线型和浮筒进行多目标优化。黄松等^[3]利用神经网络完成了动态缓波型立管的浮力块数量优化，降低了浮力块数量。刘昊等^{[4 − 6]}采用NSGA-Ⅱ算法对钢悬链立管纤维铺层、纤维层的厚度以及铺设角度进行了优化。王竹^[7]采用改进的粒子群算法（GA-PSO）完成了立管铺层与接头的优化，降低了立管的结构重量。Silva等^[8]利用梯度算法和遗传算法对复合材料钢悬链立管的冒口进行纤维厚度和铺层方向优化，节省了结构重量。现有研究主要集中在纤维复合海洋立管的构型优化设计上，而关于纤维复合海洋立管铺层优化的研究相对较少。同时，现有成果中的优化变量设计较为有限，尚未涵盖所有可能的优化方向。在这些设计中，纤维复合海洋立管的铺层结构通常采用[轴向纤维层（0°）/角度增强纤维层（±$ \theta $°）/环向纤维层（90°）]的结构形式。本文拟在已有研究的基础上，针对该结构设定9个研究变量，采用改进的粒子群算法，在Matlab和Ansys中对立管铺层进行优化。以减轻结构单位重量为目标，设计一组纤维复合海洋立管铺层设计参数，为研究人员提供实用的优化方法和可供参考的结构数据。

1 纤维复合海洋立管有限元模型 1.1 材料选择

纤维复合海洋立管结构如图1所示。其中，$ \theta $为管道环向方向；$ r $为管道径向方向；$ \phi $为角度增强纤维缠绕角度；1为纤维方向；2为垂直纤维方向；3为厚度方向。

内管选用抗腐蚀和重量更轻的聚醚醚酮（PEEK）。纤维材料选用强度高、重量轻且在海水中力学性能不会受到影响的高强度碳纤维AS4。根据基体的材料特性选用聚醚醚酮（PEEK）作为基体。PEEK材料力学性能如表1所示。王春光等^[9]结合理论模型计算和前期研究数据对纤维材料和基体进行等效处理得到单层AS4-PEEK纤维复合材料的弹性常数和长期强度，具体数值如表2所示。

1.2 有限元模型及验证

根据美国船级社ABS标准^[10]内压值155.25 MPa。由于是设计考虑的最大内压值，所以只需建立局部模型即可，可以减少大量计算时间。局部有限元模型长度为3 m，内径为0.25 m，立管环向设置80个单元，轴向设置150个单元。立管两端采用固定约束，同时考虑管内流体产生的内压端部效应。内管采用均质实体单元建模，纤维复合层采用分层实体单元建模，两者均为Solid186单元，通过关键项Keyword（3）来定义单元的均质与分层属性。有限元模型如图2所示。

为了验证所建立模型的正确性，与文献[1]中的设计结果进行了对比。文献[1]设计参数如下：内管选用铝合金，厚度为2 mm；纤维复合层选用碳纤维-环氧树脂；铺层结构为[0°₄/(±53.5°)₅/90°₄]，厚度比为1.58∶1.86∶1.62。验证结果如图3所示，可见当前模型与文献相比，安全系数偏高。考虑到以下2个因素：1）建模方式，文献采用Ansys Workbench中的ACP模块进行建模，当前模型采用Ansys经典模块通过APDL命令流建立；2）网格划分，文献中采用环轴比80∶30的划分形式，与其相比当前模型网格划分更精细。综上，2个模型安全系数的变化趋势相同数值相差不大，2种模型具有较强一致性，如图3所示。

2 改进粒子群算法

粒子群算法（Particle Swarm Optimization, PSO）思想来源于模拟自然界鸟群的觅食行为。每只鸟都被视为一个没有质量和体积的粒子，每个粒子代表优化问题的一组候选解。粒子间通过信息交互逐步向食物靠拢，粒子与食物的距离表示为适应度函数值，适应度值越小代表当前候选解的质量越好。每次迭代完成后粒子群通过位置和速度更新公式，更新个体历史最优解和全局最优解。最终经过多轮迭代获得最优解。

传统的粒子群算法存在2个缺陷：1）不能很好地平衡局部开发能力和全局探索能力；2）容易过早收敛陷入到局部最优解，特别是在搜索空间复杂、存在多个局部最优解的问题中。在实际工程中主要有以下几种改进措施^{[11 − 15]}：1）参数改进，其中包括线性、非线性和自适应方法；2）种群初始化改进；3）混合其他算法；4）调整速度和位置更新公式。粒子群算法陷入局部最优解示意图，如图4所示。

2.1 基于LHS的种群位置初始化

为了使粒子的初始分布更具有随机性和均匀性，将拉丁超立方采样（Latin Hypercube Sampling，LHS）引入粒子群算法中，增加粒子在空间分布中的均匀性和随机性。LHS是一种用于多维空间中的均匀采样方法，其基本思想是将每个维度均匀分割成若干区间，然后从每个区间中随机选择一个样本点，例如在第$ i $个区间随机选择的点为：$ {x}_{i,j}$～$ U\left(\left(i-1\right)/n,i/n\right) $。随后将上述随机点进行随机打乱，确保每个维度上的投影都是均匀分布的。LHS的数学表达式如下：

$ {{\boldsymbol X}}=\left\{{x}_{ij}\right\} , i=\mathrm{1,2},\cdots \cdots ,n ; j=\mathrm{1,2},\cdots \cdots ,d。$

(1)

式中：$ {x}_{ij} $为第$ i $个样本在第$ j $个维度的值。

2.2 自适应惯性权重及学习因子的非线性动态调整策略

针对粒子群算法局部开发与全局搜索不平衡的问题，将信息熵引入算法中，通过计算粒子熵值判断种群多样性，粒子根据信息熵值自适应动态调整惯性权重和学习因子，以平衡全局搜索和局部开发能力。在粒子群算法中，对于离散随机优化变量$ X $，其熵$ H\left(X\right) $定义为：

$ H\left(X\right)=-\sum _{i=1}^{N}{P}_{i}\mathrm{log}{P}_{i} 。$

(2)

式中：N为种群规模；$ {P}_{i} $为第$ i $个粒子的位置概率。

基于信息熵自适应参数非线性动态调整公式如下：

$ \omega\left(t\right)=\omega\mathrm{_{min}}+\left(\omega\mathrm{_{max}}-\omega\mathrm{_{min}}\right)\cdot\left(1-H/H\mathrm{_{max}}\right)。$

(3)

$ c_1\left(t\right)=c_{1_\mathrm{min}}+\left(c_{1_\mathrm{max}}-c_{1_\mathrm{min}}\right)\cdot\left(H/H_{\mathrm{max}}\right)，$

(4)

$ c_2\left(t\right)=c_{2_\mathrm{max}}-\left(c_{2_\mathrm{max}}-c_{2_\mathrm{min}}\right)\cdot\left(H/H_{\mathrm{max}}\right)。$

(5)

式中：$ H $为当前迭代次数下粒子信息熵；$ H\mathrm{_{max}} $为最大信息熵。$ \omega\mathrm{_{max}} $和$ \omega\mathrm{_{min}} $分别为惯性权重的最大值和最小值，本文取0.9和0.3。$ c_{\mathrm{1,2}_\mathrm{max}} $和$ c_{\mathrm{1,2}_\mathrm{min}} $为学习因子的最大值和最小值，本文取2和0.5。

2.3 “靠拢、变异及重分布”机制

在高维优化过程中发现，部分变量在经过几次迭代后迅速收敛至某一固定值，加快算法收敛至局部最优解。为解决这一问题，本文通过“靠拢、变异及重分布”机制对粒子运动进行扰动，旨在保持种群多样性，防止部分维度过早收敛影响算法性能。

计算每个粒子与全局最优粒子之间的距离，并以平均距离为参考，按照式(6)引导处于平均距离范围内的粒子向全局最优粒子靠拢；对于其他粒子，则根据变异概率对部分维度进行随机扰动使其发生变异，变异后重新计算其与全局最优粒子之间的距离，并依据式(7)向全局最优粒子移动；对于变异后仍被淘汰的粒子，随机初始化其位置，并重新分布在种群中。

靠拢位置更新公式：

$ D_i^{\mathrm{new}}=\alpha\cdot D_i+\left(1-\alpha\right)\cdot D_{\mathrm{gbest}}。$

(6)

式中：$ D_i^{\mathrm{new}} $为第$ i $个粒子的新位置；$ {D}_{i} $为第$ i $个粒子当前的位置；$ D\mathrm{_{gbest}} $为当前全局最优位置；$ \alpha $为靠拢系数，取值为$ \left[0，1\right] $，本文设为0.9。

变异位置更新公式：

$ D_i^{\mathrm{mutated}}=D_i+\delta\cdot{\boldsymbol{N}}\left(0,\sigma^2\right)。$

(7)

式中：$ D_i^{\mathrm{mutated}} $为第$ i $个粒子的新位置；$ \delta $为变异幅度因子；$ {{\boldsymbol N}}\left(0,{\sigma }^{2}\right) $为一个正态分布随机变量。

3 优化问题 3.1 问题描述

根据美国船级社ABS^[10]的规定，纤维复合海洋立管应进行局部设计确保其局部承载能力。局部设计阶段，假定立管处于静水环境下，承受荷载主要来自于管内流体压力、海水静水压力以及平台拉伸荷载等。其中，管内流体产生的内压荷载是主要工作荷载。因此，本文以内压荷载为工况，以减轻结构单位重量为目标，优化得到合理的铺层设计参数。由于纤维复合海洋立管由内管和纤维复合层组成，为防止立管失效，以内管等效应力$ \sigma $和纤维复合层在最大应力破坏准则下的最小安全系数$ sf_\mathrm{min} $为约束条件，其中安全系数为许用应力与实际应力的比值；并将内管厚度$ t $、轴向纤维单层厚度及层数$ {t}_{0} $、$ {l}_{0} $，角度增强纤维单层厚度及层数$ {t}_{\theta } $、$ {l}_{\theta } $，环向纤维单层厚度及层数$ {t}_{90} $、$ {l}_{90} $，角度增强纤维层增强角度$ \theta $和层叠顺序$ n $作为设计变量进行优化。

纤维复合海洋立管铺设优化问题描述为：

$ \begin{split} &\mathrm{find}\;r,t_0,t_{\theta},t_{90},\theta,l_0,l_{\theta},l_{90},n \\ &\mathrm{min}\;W \end{split} $

$ \mathrm{s.t.}\left\{\begin{aligned} &\mathop {{{t}}}\limits_ -\leqslant t\leqslant \overline{t}，\\ &\mathop {{{t}_{0}}}\limits_ - \leqslant {t}_{0}\leqslant \overline{{t}_{0}}，\\ &\mathop {{{{t}_{\theta }}}}\limits_ -\leqslant {t}_{\theta }\leqslant \overline{{t}_{\theta }}，\\ &\mathop {{{{t}_{90}}}}\limits_ -\leqslant {t}_{90}\leqslant \overline{{t}_{90}}，\\ &\mathop {{{\theta }}}\limits_ -\leqslant \theta \leqslant \overline{\theta }，\\ &\mathop {{{{l}_{0}}}}\limits_ -\leqslant {l}_{0}\leqslant \overline{{l}_{0}}，\\ &\mathop {{{{l}_{\theta }}}}\limits_ -\leqslant {l}_{\theta }\leqslant \overline{{l}_{\theta }}，\\ &\mathop {{{{l}_{90}}}}\limits_ -\leqslant {l}_{90}\leqslant \overline{{l}_{90}}，\\ &\mathop {{{n}}}\limits_ -\leqslant n\leqslant \overline{n}，\\ &sf\mathrm{min}\leqslant 1，\\ &\sigma \leqslant {\sigma }_{s}。\end{aligned}\right. $

(8)

式中：$ {\sigma }_{s} $为内管屈服应力；下划线“_”表示设计变量的下界；上划线“^—”表示设计变量的上界。

目标函数为：

$ W={\rho }_{1}\times {\text{π}} \times \left({R}_{1}^{2}-{R}_{0}^{2}\right)+{\rho }_{2}\times {\text{π}} \times \left({R}_{2}^{2}-{R}_{1}^{2}\right)。$

(9)

式中：$ {R}_{0} $为管道内径；$ {R}_{1} $为内管外半径；$ {R}_{2} $为管道外半径；$ {\rho }_{1} $为内管材料密度；$ {\rho }_{2} $为纤维复合材料密度；$ W $为结构的单位重量。

3.2 约束处理

利用罚函数法将状态变量与目标函数相联系，把约束问题转化为无约束问题。约束值转换为惩罚项，与目标函数组合在一起构成惩罚函数，目标函数值寻优转换为寻找罚函数值最小值，使得PSO算法能够有效地逼近满足约束条件的解。惩罚项构造借助文献[16]思想，采取分段非线性的构造方式。罚函数构造如下：

$ {F}=OBJ\times \left(1+P\right) \text{；} P=\left\{\begin{aligned} &0,P\leqslant 0，\\ &\sqrt{P},0 < P\leqslant 1\\ &{P}^{2},P > 1。\end{aligned}\right. , $

(10)

$ \boldsymbol{P}=\kappa\times\left(\left(1-\mathit{sf}_\mathrm{min}\right)+\frac{\sigma-\sigma_s}{\sigma_s}\right)。$

(11)

式中：$ {F} $为惩罚函数；$ OBJ $为目标函数；$ {{\boldsymbol P}} $为惩罚项；$ \kappa $为惩罚系数，本文取10；$ sf_\mathrm{min} $为纤维层最小安全系数；$ \sigma $为内管等效应力；$ {\sigma }_{s} $为内管屈服应力。

3.3 优化流程

Matlab与Ansys联合仿真步骤如下：

步骤1　 在Matlab主程序中设定粒子群算法的基本参数，并定义设计变量的上下界，利用LHS生成初始设计变量输出到文件中。

步骤2　Ansys写入变量数据，建立立管有限元模型（FEM），在内压工况下进行静力分析，将内管等效应力和纤维复合层最小安全系数输出到文件中。

步骤3　Matlab读入设计变量，根据式(9)计算目标函数$ OBJ $；读入状态变量，根据式(11)计算惩罚项$ {{\boldsymbol P}} $，最终依据式(10)计算适应度函数F，并在粒子群主程序中得到该粒子适应度值。每个粒子计算完成之后返回步骤2。

步骤4　种群中所有粒子计算完成后，获取全局最优粒子的适应度，并根据最大迭代次数Maxiter和目标精度判断是否满足终止条件。如果不满足则更新粒子速度和位置，生成新的设计变量并输出到文件，返回步骤2；否则，输出最优设计变量和目标函数值。

值得注意的是，纤维层层数和堆叠顺序均为整数型变量。因此，在优化中当Matalb写出数据时均已做整数化处理。优化流程图如图5所示。

4 优化结果

优化变量[$ t $，$ {t}_{0} $，$ {t}_{\theta } $，$ {t}_{90} $，$ \theta $，$ {l}_{0} $，$ {l}_{\theta } $，$ {l}_{90} $，$ n $]的上下界设计值为U_b=[9，3，3，3，70，5，16，5，6];L_b=[5，1，1，1，30，1，6，1，1]。其中，堆叠顺序$ n $取$ n $₁=[0°/(±$ \theta $°)/90°]；$ n $₂=[0°/90°/(±$ \theta $°)]；$ n $₃=[(±$ \theta $°)/0°/90°]；$ n $₄=[(±$ \theta $°)/90°/0°]；$ n $₅=[90°/0°/(±$ \theta $°)]；$ n $₆=[90°/(±$ \theta $°)/0°]。粒子群算法种群数量设置为30，迭代次数为100。

为测试改进粒子群算法的性能，将遗传算法（Genetic Algorithm，GA）和传统粒子群算法（PSO）作为对比。遗传算法采用Matlab中的遗传算法工具箱。3种算法的对比结果如图6所示。可知，PSO算法迅速收敛后迭代曲线不再波动，虽然收敛速度最快，但算法陷入局部最优解导致过早收敛；GA算法虽然能够在过早收敛后跳出局部最优解，但其收敛的适应度值较高；改进后的PSO算法在前期收敛速度较快，但随着迭代的进行，后期收敛速度逐渐变慢。并在该过程中成功跳出了几次局部最优解，最终收敛到一个较低的适应度值。改进PSO优化结果如表3所示。

优化后的结构单位重量较初始重量减少了41.25 kg/m³，减重52.63%；较文献[1]减轻了2.78 kg/m³，下降7%。实际工程中深水立管长达千米，若以立管长度1500 m计算，优化后单根立管重量分别减轻61.88 t和4.17 t，对减小平台荷载具有重要贡献。依据设计结果，对纤维复合层安全系数进一步分析（见图7），可得以下结论：

1） 纤维方向上轴向纤维复合层安全系数较大；在垂直纤维方向上，轴向纤维层和环向纤维层的安全系数较小，而角度增强纤维层的安全系数较大；在剪切方向上的安全系数变化则与垂直纤维方向刚好相反。

2） 纤维方向、垂直纤维方向和剪切方向的最小安全系数分别为1.26（LAYER 3）、1.00（LAYER 13）和1.68（LAYER 3）。显然，平面内的横向应力（垂直纤维方向的应力）是最关键的应力，决定了结构的设计结果。

如图8所示，对设计变量和状态变量在优化过程中的详细数据进行了提取和分析，旨在为后续研究提供数据参考。分析数据可得以下结论：

1） 由图8(a)可知，内管厚度的可行解搜索范围区域趋于下边界，增强角度高频区间置于50°~70°，堆叠顺序组合1较其他组合在种群中有更高的被选择概率。

2） 由图8(b)可知，轴向纤维单层厚度和角度增强纤维单层厚度的优秀解在2 mm左右时，而环向纤维则低一些，其区间在[1,1.5]。

3） 由图8(c)可知，轴向纤维和环向纤维的高频搜索空间均倾向于下边界。

4） 由图8(d)可知，纤维层最小安全系数是影响设计结果的关键约束。

图9为在Ansys中得到的内管等效应力云图、第3层纤维在纤维方向和剪切方向上的应力云图，以及第13层在垂直纤维方向上的应力云图。可以清晰地观察到立管在内压作用下的变形形式：立管横截面向外扩展，壁厚减薄。应力云图显示，纤维方向、剪切方向以及内管等效应力的最大值均出现在管道端部，而垂直纤维方向的应力分布与之相反。此外，剪切方向的应力明显小于其他方向的应力。

5 结　语

针对深水纤维复合海洋立管内压下局部优化设计问题的特点，考虑了9种设计变量，以纤维复合层的最小安全系数和内管的等效应力作为状态变量，以结构总重量作为目标函数，借助改进粒子群算法，通过Matlab与Ansys联合仿真完成优化设计。可得以下结论：

1） 改进后的粒子群算法性能显著提升。通过LHS方法生成初始化粒子，并借助“信息熵”调整粒子群参数，以及采用“靠拢、变异和重分布”机制干扰粒子运动，增强了全局搜索与局部开发之间的平衡能力。算法在避免陷入局部最优解方面表现出较强的能力。

2） 内管厚度的可行解搜索范围趋近于变量的下边界；当增强角度的搜索空间处于50°～69°时，结构性能表现最佳；轴向纤维层的设计厚度大于环向纤维层；[0°/(±$ \theta $°)/90°]形式的堆叠顺序相比其他铺层方式展现出更优的性能。

3） 在内压作用下，立管纤维方向的最大应力与剪切方向的最大应力发生在同一层；纤维方向上的最大应力分布位置与垂直纤维方向的应力分布位置相反。

4） 优化结果获得了更轻的设计方案，提供了一组可供设计人员参考的纤维复合材料立管结构参数，具有更低的重量。

然而，纤维复合材料海洋立管的铺层设计是一项复杂的工作，除需进行局部荷载下的局部设计外，还需进行整体荷载下的整体计算。