基于半数值方法的限弯器剪弯刚度分析

丁乐声，陈金龙，陈 潇，郭 豪，张 聪

(1. 南方海洋科学与工程广东省实验室(湛江)，广东 湛江 524006; 2. 大连理工大学 宁波研究院，浙江 宁波 315700; 3. 江苏亨通海洋光网系统有限公司，江苏 苏州 215000; 4. 华中科技大学 船舶与海洋工程学院， 湖北 武汉 430074)

柔性管缆是海洋渔业养殖[1-2]、油气和风电[3]等海洋资源开发的关键装备之一。其在多种载荷共同作用下易产生过大曲率，从而导致弯曲失效[4]。因此美国石油协会(API)发表过相关规范[5]，其中指出柔性管缆易发生弯曲失效的位置需要安装限弯器。并且需对不同工况的限弯器进行整体分析[6]，以确定其抵抗变形能力符合使用要求。

限弯器是由多组同尺寸变径半圆筒式子结构相互嵌套而成，当其子结构受力将相互锁合，从而限制弯曲变形进一步发生[7]。在限弯器锁合过程中，其锁合结构具有接触非线性，进而导致剪弯刚度分析困难。针对该类问题，部分研究人员进行了相关研究。其中理论方面，Souza和Ramos[8]建立了剪弯荷载下连续异形梁结构受弯的理论模型，但无法应用于限弯器这一包含接触非线性的梁结构。张聪[9]建立了限弯器剪弯刚度的整体数值分析方法。Noh和Ha[10]运用ABAQUS的接触分析模块[11]进行了限弯器与缆的多体接触分析，并进行限弯器结构拓扑优化设计。

综上所述，目前限弯器这一包含接触非线性梁结构的剪弯刚度分析计算都基于整体数值方法，无法直接给出限弯器截面尺寸对剪弯刚度的影响，导致限弯器剪弯刚度设计多基于经验。为此提出剪弯刚度半数值分析方法以便通过公式直接分析部分参数与剪弯刚度关系，并满足限弯器剪弯刚度迭代设计的快速分析需求。先基于悬臂梁方程，建立包含锁合结构扭转刚度的剪弯刚度理论模型；然后通过构建考虑接触边界的锁合结构数值模型，得到其扭转刚度曲线；最后结合理论分析和数值计算求得考虑接触非线性的限弯器剪弯刚度曲线，并同传统的整体数值方法相比较，验证该半数值分析方法的有效性。

1 剪弯刚度理论分析

如上所述，限弯器是典型的非线性短梁结构。传统的限弯器剪弯刚度采用整体数值分析，先建立三维悬臂梁有限元模型计算该结构在剪力作用下的端部挠度；然后使用悬臂梁挠度方程计算出剪弯刚度。如图1所示，限弯器的整体结构可分为子梁结构和锁合结构；在弯曲过程中子梁结构主要发生挠曲变形，锁合结构主要发生扭转变形[12]，其中锁合结构发生的扭转会导致梁结构端部发生垂直于轴线方向的位移进而使结构挠度增加；因此限弯器的剪弯刚度受子梁结构弯曲刚度和锁合结构扭转刚度共同影响。

进行理论分析时将限弯器结构简化为由两部分组成的端部受剪力P的周期性悬臂短梁结构，其中直线为子梁结构、圆点为锁合结构，如图2所示；限弯器端部的总挠度w由所有子梁结构挠度w1与所有锁合结构扭转所对应挠度w2组成，即w=w1+w2。子梁结构是典型的直梁，可引入小变形与直线梁假设；同时考虑其截面剪切变形，由Timoshenko梁理论[13-14]进行悬臂梁结构挠度计算，如式(1)所示。值得注意的是限弯器结构中的锁合结构会使整体产生转角，导致在引入直线梁假设后基于长度计算力矩时产生误差，但当其圆心角较小时该影响可忽略[15]。

图2 限弯器结构简化悬臂梁模型Fig. 2 Simplified cantilever beam model with bend restrictor

(1)

式中：w1为所有子梁结构挠度，P为端部剪力，E为材料弹性模量，n为子梁结构节数，l为限弯器有效节长，I为梁截面惯性矩，c为剪切系数，Mr为剪切模量，A为梁截面面积。若忽略梁截面剪切变形，则式(1)退化为基于Euler-Bernoulli梁理论的悬臂梁挠度方程(2)：

(2)

(3)

(4)

n节子梁结构需要n-1节锁合结构连接，进一步求得有锁合结构扭转所对应挠度w2。

(5)

(6)

同理，将式(2)、(5)联立，并也代入总挠度公式和等效剪弯刚度方程，可得到基于Euler-Bernoulli梁理论的忽略梁截面剪切变形的剪弯刚度K计算方法。

(7)

对比可知，方程(6)仅分母比方程(7)多一项，该项为短梁截面的剪切应变造成，并恒为正值；因此基于Timoshenko梁理论的计算方程(6)更适用于短梁结构，并且结果更为保守。同时也需要注意，以上理论模型是引入了直线梁假设和小变形假设建立的，限弯器本身的曲率与变形将导致该模型存在误差。

2 限弯器剪弯刚度计算

2.1 算例几何与材料参数

限弯器在进行测试和使用时常因节省成本而选择较短的组装长度，以5节典型限弯器子结构为例，子结构形状如图3所示。使用过程中，限弯器子结构头尾相互嵌套于海缆外部，以提高局部剪弯刚度。构件的几何参数以及材料力学性能参数在表1中给出，并以5节限弯器子结构嵌套为例进行计算。

图3 限弯器算例尺寸Fig. 3 Bend restrictor calculation example

表1 限弯器结构参数Tab. 1 Bend restrictor structure parameter

2.2 半数值分析

根据理论模型，剪弯刚度求解需要锁合结构扭转角，但限弯器中的每个锁合结构包含两个接触边界，使其计算较为复杂。因此使用半数值方法进行剪弯刚度分析，具体的流程[16]如图4所示，先基于有限元软件ABAQUS建立限弯器锁合结构局部三维模型，并计算提取不同力矩下锁合结构扭转角度；同时基于理论分析中的式(1)或(2)使用不同梁理论计算子梁结构弯曲挠度；再进一步通过式(3)分析不同力矩下的锁合结构扭转刚度，通过式(5)分析锁合结构所对应的挠度，进而算出限弯器端部挠度；最后再将其与各参数代入式(6)或(7)基于Matlab软件计算限弯器剪弯刚度。

图4 剪弯刚度半数值分析流程Fig. 4 Shear bend stiffness analysis process

限弯器锁合结构可简化为平面对称结构，因此建立半体三维模型[2-3, 9]进行数值分析以减少计算量。模型在锁合结构末端设置固定约束，并在半体模型轴对称面设置对称约束；在进行锁合结构局部数值分析时使用扭矩来等效5节限弯器端部受剪力造成的力矩，分析中限弯器端部受到的最大剪力为20 kN，则等效的最大扭矩为21.6 kN·m；又由于采用半体三维模型，实际该模型端部z轴方向施加最大扭矩的一半；具体如图5所示。

图5 锁合结构相互作用及边界条件设置Fig. 5 Interaction and boundary condition setting of locking structure

该模型有两个接触对，法向行为选择硬接触算法，切向行为则选择罚函数算法。网格划分单元选择分析耗时短的C3D8R单元，共3 980个。该模型每迭代100次的计算时间为96 s。通过有限元分析方法不仅可以直接计算锁合结构扭转变形进而计算弯曲刚度，还可以直接读取力矩荷载达到21.6 kN·m时锁合结构各接触面的接触应力云图，如图6所示。从图6中可知，锁合结构接触面较小且接触压力分布不均匀[2]，这都表明该处有较强的接触非线性。

图6 锁合结构各部分接触压应力云图Fig. 6 Contact compressive stress cloud diagram of each part of locking structure

从分析结果中可提取不同力矩的锁合结构转角，进一步可基于数值分析得到力矩—转角关系以及式(3)计算限弯器锁合结构非线性扭转刚度，如图7所示。图7(a)力矩与转角的关系曲线与切线相比，具有明显的非线性；从几个给定力矩条件下的母头端部后侧接触压力云图可以看出在力矩较小时接触面积随着力矩快速增加；而力矩超过10 kN·m后接触面趋于稳定，随着力矩的增加接触非线性减弱，力的传递主要体现在接触压力的增加。图7(b)中扭转刚度的变化趋势进一步说明锁合结构接触非线性在加载初期较为明显，随着荷载增加非线性逐渐减弱。

图7 限弯器锁合结构转角与刚度Fig. 7 Torsion angle and rigidity of the bend restrictor locking structure

在半数值分析中为考虑锁合结构的扭转刚度非线性，采用四阶多项式对局部数值分析的结果进行拟合，得到R平方系数[17](确定系数)为0.997的扭转刚度函数G(M)；其拟合结果在图7(b)中给出。最后根据计算流程，将扭转刚度函数G(M)与计算方程(1)、(5)或(2)、(5)分别联立得到端部挠度，进一步代入等效剪弯刚度方程即可获得基于式(6)与(7)两种理论的限弯器剪弯刚度。

2.3 整体数值分析

传统方法中限弯器的剪弯刚度采用整体数值分析[18-19]；利用结构对称性建立上述算例的半体三维模型[5]，选择与半数值分析方法相同的悬臂梁剪力加载，在根部施加完全固定约束，在对称面设置对称约束，并对半体模型端部施加10 kN的集中力，进而模拟出实际结构根部21.6 kN·m的力矩，如图8所示。

图8 限弯器相互作用及边界条件设置Fig. 8 Interaction and boundary condition setting of bend restrictor

模型共有8个接触对，接触算法以及单元属性与限弯器锁合结构数值分析模型相同；模型共包含28 390个C3D8R单元，是限弯器锁合结构模型的7倍，其中锁合结构处的单元尺寸和分布与锁合结构局部数值分析模型相同。整体数值分析可得到限弯器剪弯刚度用来验证半数值分析方法。并且整体数值分析能计算根部力矩和端部转角关系，如图9所示，进一步说明了数值方法可以有效分析接触非线性问题。

图9 限弯器根部力矩—端部转角曲线Fig. 9 Curve of bend restrictor root moment-end rotation angle

通过整体数值分析还可以进一步分析小变形假设与直线梁假设对该结构的影响。限弯器装配体最大对数应变和挠度分别如图10、图11所示。从图10中可以看出该结构整体应变较低，应变峰值出现在靠近悬臂梁末端的子梁结构与锁合结构交界的应力集中处，且均小于5%；从图11中可以看出限弯器装配体结构最大挠度为140.5 mm，比其整体结构长度1 174 mm低一个数量级。应变所代表的局部变形与挠度代表的整体变形远小于结构原始尺寸[20]，因此理论分析引入小变形假设所造成的误差是有限的。限弯器装配体中子梁结构是典型的沿中性层对称的直梁，直线梁假设导致的误差主要由分析整体结构力矩时产生，但该结构对应的圆心角为12.8°，而当圆心角小于22°30′时曲线梁荷载可以近似的按直梁计算[15]，因此理论分析引入直线梁假设所造成的误差也是有限的。

图10 限弯器对数应变云图Fig. 10 Logarithmic strain cloud diagram of bend restrictor

图11 限弯器挠度云图Fig. 11 Deflection cloud diagram of bend restrictor

3 各算法对比

3.1 计算效率对比

限弯器设计过程中常进行迭代计算，以获得较优的设计[10, 21]。因此计算效率是限弯器分析方法的关键指标。上文所述半数值分析与整体数值分析均在i3-8350k CPU中进行，其耗时于表2给出。可知Euler-Bernoulli梁理论与Timoshenko梁理论2种计算方法效率极其接近，而半数值分析效率远高于整体数值分析。

表2 各算法计算效率对比Tab. 2 Calculation efficiency comparision of various algorithms

3.2 计算结果对比

剪弯刚度是限弯器的关键指标，传统计算采用整体数值方法；同时其也可基于上文所述的半数值分析流程计算。通过不同的计算方法得到了相应的限弯器剪弯刚度曲线如图12所示，可知基于Euler-Bernoulli梁理论与Timoshenko梁理论的两组不同半数值计算结果与传统数值方法的平均相对误差均小于6%，并且两种半数值方法均可以较好地反映出限弯器锁合结构的接触非线性，Timoshenko梁理论对应的算法由于考虑了短梁的剪切变形更加保守。

图12 限弯器剪弯刚度Fig. 12 Shear bend stiffness of bend restrictor

半数值方法的误差主要由理论模型中引入的小变形假设和直线梁假设引起，进一步计算基于不同梁理论的两种半数值方法与传统数值方法的误差，于表3中给出；证明了Timoshenko梁理论对应的算法比Euler-Bernoulli梁理论准确，且Timoshenko梁理论的平均相对误差为3.398%。

表3 两种理论算法误差对比Tab. 3 Comparison of two theoretical calculation errors

4 结 语

基于Timoshenko梁理论和Euler-Bernoulli梁理论剪力限弯器剪弯刚度半数值分析方法，并采用该方法和传统的整体数值方法分别进行了剪弯刚度分析并比较，得到如下结论：

1) 在限弯器剪弯刚度分析时半数值方法的计算效率远高于传统数值方法；并且基于不同梁理论半数值方法的平均相对误差均小于6%。

2) 限弯器的锁合结构处受到接触边界的影响呈现接触非线性，这对限弯器锁合结构刚度造成很大的影响；并且其非线性在加载初期较为明显，随着荷载增加逐渐减弱。

3) 限弯器的子梁结构是典型的短梁结构，因此基于Timoshenko梁理论的半数值分析更适用于限弯器剪弯刚度分析，结果误差更小且更保守。

综上所述，基于半数值方法可以高效地分析限弯器剪弯刚度，其结果可与数值方法相印证，为相关分析与设计提供参考。但是本文中所使用的梁理论模型需引入直线梁假设和小变形假设，仍对分析结果造成一定误差，可以在未来进一步研究。