0 引　言

复合材料以其结构轻质、强度高、优质的防蚀特性，在工程领域得到广泛应用。船舶制造业为追求轻量化，逐步采用复合材料建造船舶^{[1 − 3]}。然而，全复合材料船舶建造成本高，并且会降低整体强度，因此，将复合材料应用于上层建筑，主船体仍采用钢材，这样较为合理。这使得上层建筑与主船体连接所用的钢-复合材料接头成为研究重点。材料的连接方式中，机械连接方式如铆接会引起应力集中，导致结构损伤，而胶接可使应力分布均匀和减轻结构重量^[4]，故在复合材料结构中占据重要地位^{[5 − 7]}，巴沙木是一种原产于南美洲的热带木材，它以其优良的低密度、抗震、隔音、隔热、抗疲劳、生长周期快等优点著称，尽管现如今尽管有大量的竞争性夹层芯材（例如聚合物泡沫、金属泡沫、蜂窝）, 但巴沙木仍然是夹层结构施工的重要材料^{[8 − 9]}。

因此，研究钢-复合材料胶接接头的力学性能及其结构优化，对提高其承载能力至关重要。对船用巴沙木夹芯钢-复合材料接头虽有较多研究，但有限元仿真难以准确模拟其失效点与失效模式，主要因巴沙木材料力学特性复杂。本文通过Abaqus软件中的VUMAT子程序编写巴沙木失效准则，准确模拟该接头的力学特性与失效形式，之后利用Abaqus与Python联立创建建模计算脚本，凭借ISIGHT软件里集成的实验设计（DOE）模块中的最优拉丁超立方算法，捕捉不同参数下结构的力学性能设计区域并输出样本点，根据样本点建立响应面模型，最后以质量最小、力学性能最优为目标，在多目标优化问题里，NCGA 遗传算法展现出了有效的优势，故以此来确定最优解。

1 巴沙木芯钢-复合材料胶接接头

巴沙木夹芯钢-复合材料接头^[9]主要被应用于船舶甲板延伸部分以及甲板与上层建筑的连接之中。该接头的设计为了契合特定的应用需求，需要维持一侧的光滑性^[8]。接头的左侧为甲板钢板，右侧则为上层建筑舷侧或甲板外延。此结构由复合材料外蒙皮、胶层、巴沙木夹芯以及钢板共同构成（见图1），复合材料外蒙皮由4层玻璃纤维复合材料组成，其铺层方式为[0°/90°]₂，整体采用胶接连接方式。

2 模型材料参数与求解方法 2.1 巴沙木材料参数及失效准则

木材作为一种具有高度不均匀性和各向异性的材料，在拉伸、剪切等负载条件下分别呈现出脆性特征、塑形延性损伤，给数值模拟仿真带来了严峻的挑战。本文基于Sandhass^[10]和Linde^[11]的工作成果，凭借 Yamada-Sun 屈服准则对木材屈服的状况予以判定。基于此准则，开发一种新的木材累积损伤演化模型，特别是针对木材在受拉和受压条件下的脆性破坏与延性破坏的差异进行了深入的考虑。巴沙木的径向和切向参数相同。每个损伤变量都遵循特定的线性软化规则。表达式如下：

$ {{d}}\left( k \right) = 1 - \frac{1}{k} ，$

(1)

$ d(k) = 1 - \frac{1}{{f_{\max }^2 - 2{G_f}E}}\left( {f_{\max }^2 - \frac{{2{G_f}E}}{k}} \right)。$

(2)

式中：$ k $、$ f_{\max }^{} $、E和$ {G_f} $分别为状态变量、强度、弹性常数和临界能量释放率，状态变量$ k $用来跟踪整个加载历程，在每个时间增量点，其表达式可以表示为：

$ {k^t} = \max \left\{ {\left. {1,{{\max }_{increment = 0,t}}\left\{ {\left. {{F^{increment}}} \right\}} \right.} \right\}} \right.。$

(3)

$ {F^{increment}} $为表中失效标准的数量，基于有效应力的概念，损伤材料的本构关系可表示为：

$ \varepsilon = {{\boldsymbol{C}}^{el}}{\sigma ^{ef}} = {{\boldsymbol{C}}^{el}}M\sigma = {{\boldsymbol{C}}^{dam}}\sigma 。$

(4)

式中：$ {{\boldsymbol{C}}^{dam}} $为损伤柔度矩阵。而在VUMAT子程序的实现过程中，需要通过对$ {{\boldsymbol{C}}^{dam}} $求逆矩阵来获得损伤刚度矩阵$ {{\boldsymbol{D}}^{dam}} $。

Yamada-Sun材料失效准则及巴沙木材料参数如表1和表2所示，下标0为沿木材纹理（即树木生长）方向，这个方向的材料弹性模量以及抗拉抗压强度是最高的。而下标90为木材的横截面方向，其中90、R和90、T分别为径向和切向。此外，下标t和c分别为拉伸和压缩；v为剪切。

2.2 玻璃纤维复合材料参数及失效准则

复合材料具各向异性和多相复合体特性，破坏过程复杂，失效位置和方式与各向应力相关。因此，选择合适的失效准则至关重要。复合材料的常用评估标准有最大应力、最大应变、蔡-希尔、霍夫曼、蔡-吴和Hashin等，Hashin准则能准确反映复合材料的破坏特性，尤其是层间分裂现象，被广泛采用。本文采用 Hashin3D失效准则来判断复合材料蒙皮失效^[12]。关于Hashin3D失效准则及刚度退化理论公式可参考文献[12]。

如表3所示：$ {X_t} $、$ {X_c} $、$ {Y_t} $和 $ {Y_c} $分别为单层板材在纵向上的拉伸与压缩强度，以及在横向上的拉伸与压缩强度；S₁₂、S₁₃和S₂₃分别为纵向与2个横向方向上的剪切强度。当以上各式的值达到1时，则材料开始失效。

2.3 胶层材料及钢材参数

本文中胶层的脱粘效应通过内聚区域模型（Cohesive zone models, CZM）^[14,15]来模拟。CZM采用的是双线性牵引-分离律，如图2所示，其材料特性表现为线性弹性和线性软化阶段。其中，$ \sigma $、$ {\sigma _{\max }} $分别为胶合面胶合应力与最大胶合应力，$ {\sigma _{\max }} $为82.5 MPa；$ {\delta _0} $、$ {\delta _{\max }} $分别为初始分离应变与最大分离应变，$ {\delta _{\max }} $=0.0027 mm。由于篇幅限制，关于胶层失效及B-K准则详请参考文献[15]。

本文参考N. Kharghani等^[9]文章中所使用的AH36钢材参数（见表4）进行数值模拟。

3 有限元仿真及形状优化 3.1 有限元模型及数值仿真验证

通过有限元软件Abaqus，进行对该接头在三点弯曲条件下极限承载能力的数值模拟。钢板、复合材料蒙皮以及巴沙木夹芯均采用了八节点六面体线性减缩积分单元（C3D8R），如图3所示。参考 Kharghani等^[9]的研究，胶层的厚度设定为1 mm，采用八节点三维粘结单元（COH3D8）。在此模型中，约束条件为：模型左右两端均设为简支边界条件。左端固定3个方向的位移，即在x、y、z轴上均受限制，右端仅在接头厚度方向和垂直于截面方向施加位移约束。此外，在距左端275 mm的上方，施加多点梁约束（MPC），并将此处的中点设为控制点。为实现模拟，对该控制点施加平滑过渡的位移载荷，通过这种载荷条件来实现模拟三点弯曲载荷。

最终得出结果如图4所示，其中，失效云图中的具体失效单元位置与实验中失效位置相吻合，由于Kharghani等^[6]数值仿真并未给出具体失效点，故结果与他们的实验结果对比，如图5和图6所示.试验极限失效平均载荷为9823.55 N；计算极限失效载荷为11073.4 N；极限失效载荷误差约为11.2；试验刚度平均值为942.76 N/mm；计算刚度值为1.01$ \times $10³ N/mm；误差约为6.7%。试验失效位移为13.8 mm；计算失效位移为12.1 mm；误差为12.3%。验证模拟计算精准度较高，具有参考价值。

在结构即将失效时，钢板最大应力为393.4 MPa，说明钢板已发生塑性形变，在力-位移曲线中的具体表现为斜率下降。故试件的初始软化是由钢板的局部屈服所引发的。

直至结构完全失效，复合材料蒙皮也并未发生明显失效。结构的失效是发生在巴沙木芯处，故对巴沙木芯进行应力分析，其即将失效的应力云图如图7所示。巴沙木芯简易受力图如图8所示。

在木芯与钢板蒙皮交界的尖端处应力较大，中部和两端已有单元失效（失效单元自动删除）。在木芯与钢板末端的交界面，尤其是钢板的折角处，出现明显的应力集中现象。另外，在结构的右侧简支上端，出现了应力较大的区域。

3.2 失效原因分析

在木芯与钢板、复合材料蒙皮交界尖端处的应力集中主要由于角度过小所引起，导致尖端处的单元过早失效。在木芯与钢板末端的交界面处，由于三点弯曲载荷，钢板传递的力导致木芯在钢板末端处发生了较大的应力集中。最后，右端简支上端应力大的原因主要是钢板向下传递的力，由于简支处的垂向自由度被固定，形成了一对力偶，从而在简支处产生了较大的剪切力。所以结构失效主要原因为2点：1）木芯尖端角度过小而导致的应力集中；2）在木芯与钢板末端的交界面处发生了大的应力集中。

3.3 钢板末端优化及不同结构转角分析

基于以上结果分析，钢板完全延伸将改善木芯的受力情况^[13]。故在保持结构右端厚度不变的情况下，建立钢板转角分别为15°、30°和45°的接头结构，依照与之前同样的数值模拟方法，对其进行对比，结果如图9所示。

其完全失效时失效云图如图10所示。完全延伸结构力-位移曲线如图11所示。界面优化后，15°结构的极限承载能力未提升，虽避免了界面折角的应力集中，但增加了钢板对木芯的垂向力，导致结构更早发生纵向剪切断裂失效。木芯角度增大，结构极限承载能力提升，这源于优化了钢板与木芯界面，增大了木芯角度，改善了其受力状况。从 45°木芯失效前的应力云图看，应力分布较均匀，但与原始结构的极限承载能力相比，差距不明显。由于结构右端厚度不变，增大木芯转角会减少钢材的材料占比，使结构整体刚度下降。

4 优化方法 4.1 DOE实验设计

基于Isight软件（DOE）组件对钢板完全延伸的巴沙木芯钢-复合材料接头进行参数分析，参数定义如图12所示。研究上下蒙皮的厚度（H₁）、钢板的厚度（H₂）、右端木芯的厚度（H₃）以及结构转角（H₄）4个参数对该结构的重量（mass）、极限承载力（RF_max）和刚度（Stiffness）的影响。

$ \begin{split} mass = & 100 \times ({H}_{2}\times 240\times {\rho }_{1} + ({H}_{3}/ \sin ({H}_{4}) - {H}_{2}/\sin ({H}_{4}))\times \\ & {H}_{2}\times {\rho }_{1}+(230\times {H}_{3}-{H}_{3}/ \tan ({H}_{4})\times {H}_{3}\times 0.5)\times {\rho }_{2}+ \\ & 470\times {H}_{1}\times {\rho }_{3}+({H}_{3}+2\times {H}_{1}-{H}_{2})/ \sin ({H}_{4})\times {H}_{1}\times \\ & {\rho }_{3}+(230-({H}_{3}+2\times {H}_{1}-{H}_{2})/ \tan ({H}_{4})。\\[-1pt]\end{split} $

(5)

式中：$ {\rho _1} $为钢材密度；$ {\rho _2} $为巴沙木密度；$ {\rho _3} $为玻璃纤维复合材料密度。

贡献图（见图13）中描绘了4个参数针对3个优化目标的贡献值情况。就重量（mass）而言，贡献最为显著的是钢板厚度（H₂），而结构转角（H₄）与质量呈负相关关系；极限承载力（RF_max）和刚度（Stiffness）与4个参数的影响关系大体一致。结构转角（H₄）与钢板厚度（H₄）对于二者的贡献值较大。

4.2 RSM响应面模型

鉴于响应面模型具备良好的可导性，优化时具有较快的收敛速度，能够拟合繁杂的响应关系，此文选用RSM模型来模拟巴沙木芯钢-复合材料接头的力学性能响应模型，其中上下蒙皮的厚度（H₁）、钢板的厚度（H₂）、右端木芯的厚度（H₃）以及结构转角（H₄）有关结构未退化刚度（Stiffness）和极限承载力（RF_max）的响应面图如图14所示。

针对误差分析，能够引入复相关系数（R²）来衡量响应面与样本点的契合程度，当R²值越发趋近于1时，意味着所构建的 RSM 响应面模型就越发可信。本文随意选用了25个样本点，算出与其相对应的极限承载力、未退化刚度以及质量响应值，对它们实施均方差归一化处理，呈现出了极限承载力和刚度实际值与预测值的误差分析图（见图15），样本点均匀地分布在预测值线的两边，而且较为贴近，两者的R²值分别为0.88934和0.85132。所以，该模型可以被认为是较为可靠的。

4.3 NCGA优化

在多目标优化中，各个优化目标之间或许会相互冲突，其优化解是一个Pareto最优解集，与之对应的目标函数空间的象即为Pareto前沿。本文采用了 NCGA 多目标遗传算法，该算法将问题的可能随机解构建成多个基因种群，并且借助排序分组来达成交叉，实现“相邻繁殖”的成效，从而使靠近 Pareto 前沿的解进行交叉繁殖的概率得以提升，加快收敛速度。考虑到巴沙木芯钢-复合材料接头应用于主甲板和上层建筑舷侧的连接部位，优化目的在于其轻质量以及安全可靠性，所以将钢板厚度（H₂）变量的范围设定为 5.5～7.5 mm，数学优化模型如下。

$ 约束条件\left\{\begin{array}{l}1.8\;{\rm{mm}} \leqslant {H}_{1} < 3.0\;{\rm{mm}}，\\ 5.5\;{\rm{mm}} \leqslant {H}_{2} < 7.5\;{\rm{mm}}，\\ 28.0\;{\rm{mm}} \leqslant{H}_{3} < 40.0\;{\rm{mm}}，\\ 20.0^{\circ} \leqslant{H}_{4} < 50.0^{\circ}，\end{array} \right.$

(6)

$ 优化目标\left\{\begin{array}{l}\mathrm{min}({\mathrm{mass}})，\\ \mathrm{max}({F}_{\mathrm{max}})，\\ \mathrm{max}(\text{stiffness}）。\end{array} \right.$

(7)

经过301次迭代优化，最优解为第212次迭代结果。根据Isight后台历史监控可得优化过程中响应目标质量、刚度、极限承载力随迭代变化过程如下图所示。优化结果显示，当H₁=2.56 mm，H₂=7.45 mm，H₃=34.84 mm，H₄=43.01°为最优参数组合，此时结构质量为2.2 kg，极限承载力为9380 N，未退化前的刚度为458.83 N/mm。

4.4 优化结果参数取整的有限元结果

基于工程实际，将最优解中的4种参数取整，H₁=2.6 mm，H₂=7.5 mm，H₃=35 mm，H₄=43°，进行有限元计算得到结果，力位移曲线见图16。结构质量为2.32 kg，极限承载力为8324.3 N，未退化刚度为392.15 N/mm。与10 mm钢板原始结构相比，质量下降了38.8%，极限承载力下降22.4%，而相比于7.5 mm厚钢板的原始结构，质量降幅为17.4%，极限承载力上升13.1%。

5 结　语

1）本文依据Yamada-Sun失效准则编写了巴沙木VUMAT子程序，构建了胶接巴沙木夹芯钢—复合材料连接接头模型，并从极限载荷、结构抗弯刚度和失效模式3个方面将计算结果与实验结果进行了对比，以此验证了子程序和数值仿真计算的准确性。

2）分析原始结构的失效原因，主要问题在于钢板与木芯连接处的转角过小，以及钢板末端界面曲率变化过大，这两者都导致了应力的集中。因此，对结构进行了优化。保持右端的高度不变，将钢板末端改为完全延伸的形式，并通过对15°、30°和45°的转角结构进行了仿真模拟。分析可知，适当增加木芯转角，可以改善木芯受力条件，从而提升结构的极限承载能力。

3）由于初始结构钢板厚度为10 mm，为了上层建筑质量从而降低船舶重心，将钢板优化范围设置为5.5～7.5 mm，最终得出最优结果。与10 mm初始结构相比，质量下降了38.8%，极限承载力下降22.4%，而相比于同钢板厚度原始形状结构，质量降幅为17.4%，极限承载力上升13.1%。质量与力学性能优化明显。