管管间接触挤压变形的数值模拟方法研究①

刘巨保 陈 健 姚利明

(东北石油大学机械科学与工程学院)

在石油石化工程和装备中，存在一类管束结构，如井下油管与套管、管壳式换热器中的换热管、海洋钻井隔水管等。这些管束结构在外力作用下，常发生接触挤压或振动碰撞而产生变形磨损甚至破坏，危及到管束的安全运行。此类问题属于接触变形问题，且往往伴随材料非线性和(或)几何非线性，难以得到准确的解析解[1～5]。因此对管束结构接触变形的全过程进行三重非线性仿真分析是非常必要的。

在以往的管束接触和碰撞问题数值计算中，一般不考虑接触的局部变形，采用横截面不变形的梁单元建立数值分析模型[6，7]，由于模型的限制无法考虑局部变形的问题。而对于采用实体单元建模的接触碰撞问题研究中，更多的是关注接触力、冲击载荷和能量损失，而未考察结构自身的具体变形规律情况[8，9]。在材料成型领域，多采用无间隙接触的挤压工艺对产品进行成型[10，11]，这在数值模拟计算的收敛上降低了计算难度。

笔者以平行管束为研究对象使用实体单元与梁单元结合建立数值模型，采用增广拉格朗日算法，引入了管挤压变形产生的大位移几何非线性和材料非线性效应，对有间隙的管管接触变形进行非线性仿真分析。

1 管管间接触挤压变形分析力学模型

图1所示为两根平行管束，将上管束定义为1管，下管束定义为2管。1管为简支约束，2管两端全约束，在1管轴线中点处施加一个集中载荷，使它变形并与两管发生接触，进而挤压两管也产生变形。两管长度均为L，间隙为g，两管外径相同均为D，内径不同，1管内径为d1，2管内径为d2，1管材料弹性模量为E，惯性矩为I。图2为管管接触的3种状态，在第2种管管刚接触的状态下，由于接触区域较小，可假设为点接触，且接触位置位于L/2处。

图1 管管间接触挤压力学模型

图2 管管接触的3种状态

由材料力学叠加原理建立平衡方程为：

ωF-ωC=g

(1)

(2)

(3)

(4)

其中ωF，ωC分别为外载荷F与接触力Fc对1管挠度的贡献的大小。联立上述方程可以解得此时的接触力：

(5)

继续增大外载荷接触区域也会增大，出现第3种状态，2管横截面发生大变形，此时上述假设不再成立，所得接触力计算结果也不具参考性，应采用数值分析方法进行计算。

2 管管间接触挤压变形分析的数值方法

2.1 数值模型

对于梁单元的有限元分析不能考虑到梁横截面上、内部各点的应力位移，需要实体单元建立模型。但如果全部采用实体单元，由于模型几何尺寸较大使得有限元模型的节点数过多，造成计算时间过长。笔者首先使用梁单元试算上述例题获得大致的接触位置，之后使用实体单元代替接触管和目标管接触段的梁单元并将实体单元网格细化，其他位置仍采用梁单元，进行混合建模。采用节点接触绑定连接两种单元，以保证两种类型单元在力和位移上的正确传递。

图3为数值模型，将实体单元段的接触管与目标管轴线方向的单元长度设置为3mm，将圆周分成18份，将目标管横截面分成3份，接触管分成1份。图4为模型中的接触对，将2管与下方平面设置为无间隙接触，将1管与2管的外表面设置为含间隙接触，由于挤压变形将造成2管内壁发生自接触，应对2管内壁上下设置为含间隙接触。

图3 数值模型示意图

图4 数值模型接触对示意图

其中接触管(忽略其横截面变形)材料定义为弹性模量E1=200GPa的弹性材料，将目标管材料定义为弹性模量E2=71GPa、屈服强度为280MPa、切线模量为500MPa的双线性弹塑性材料。应用Von Mises屈服准则和随动强化模型。通过多个载荷步加载，一个载荷步卸载。

2.2 数值接触算法

对于接触算法即将有限元法的约束变分原理应用于接触问题的求解，即分别采用拉格朗日乘子法、罚函数法和增广拉格朗日乘子法。对于包含接触界面的接触问题，泛函可以表示为：

Π=Πu+ΠCL

(6)

其中Πu是原问题中不包括接触约束条件的总位能，ΠCL是不同算法引入约束条件的附加泛函。增广拉格朗日乘子法中取：

(7)

其中Λ′=diag(a1′,a2′,a3′)；a1′、a2′和a3′为罚系数，为设定常数；λ′={λ1′λ2′λ3′}T，g′={g1′g2′g3′}T。

增广拉格朗日乘子法是为了找到精确的拉格朗日乘子(即接触力)，而对罚函数法进行一系列修正迭代。与罚函数法相比，增广拉格朗日乘子法的优点是容易得到良态条件，对接触刚度的敏感性较小。缺点是，可能需要更多的迭代，特别是变形后网格变得太扭曲。

由于笔者研究的管管接触不会出现扭曲畸形网格，且该问题经常因为接触界面耦合关系限制罚参数的取值，容易引起振荡。所以笔者采用对接触刚度敏感性小的增广拉格朗日乘子法进行有限元管管接触问题的计算方法研究[10]。

2.3 非线性方程的收敛计算

由于本例题考虑了3种非线性效应，除网格的划分、接触对的正确建立外，难点在于非线性方程组的收敛，如何快速并准确地求解非线性方程成为了本题的关键。

由于非线性问题需要采用迭代逼近技术进行求解，通常采用牛顿-拉普森法计算，其收敛速度较快，但和载荷子步、非线性设置有关。除此之外还需要对载荷子步的收敛判定方法和准则进行定义。图5给出了静力学非线性方程求解应设置的选项。

图5 静力学非线性方程求解设置选项

2.3.1载荷子步相关设置

由牛顿拉普森迭代法通过计算各载荷(位移)增量的结果，进而得出极限载荷下的计算结果，步长的设置对计算能否收敛影响较大。图6是通过试算得到的载荷作用点的力位移曲线，第1个载荷步为0～20kN，最小时间增量为0.001，最大时间增量为0.005。第1个载荷步为20～40kN，最小时间增量为0.005，最大时间增量为0.050。并将自动时间步设置为程序自行选择。线性方程求解器设置为程序自行选择。对于收敛的判定，例题采用力收敛准则，残差的计算方法采用二范数，收敛容差为0.001。如果子步达到某一物理限制，程序将视此步长过大，进而自动将子步二分，例题将这一物理限制定义为子步最大迭代次数设置为25次。由于例题中结构产生了塑性应变，为了更好地控制时间步长上的缩减，增加另一个物理限制即最大塑性应变增量限制，将其值设置为0.1。

图6 载荷作用点力位移曲线

2.3.2非线性设置

ANSYS提供了4种切向刚度矩阵修改方法，并提供了一系列的非线性设置来保证非线性求解的收敛。例题中对于各子步切向刚度矩阵的修改采用程序自行选择。非线性设置选项包括自由度预测、自适应下降和线性搜索。其中自由度预测适用于无位移突变的非线性分析，可加速收敛。由图6可以看出力位移曲线非线性响应相对平滑，故例题中激活自由度预测。自适应下降与线性搜索均采用程序自行选择。

2.3.3实常数的设置

文中的算例将接触刚度系数设置为0.01，穿透容差系数设置为0.1，收敛速度较快且得到的计算结果精确度较高。

2.3.4弧长法的使用

由图6可以看出在极限载荷段，载荷变化较大但位移变化较小，出现力位移曲线趋于水平的现象，这可能造成切向刚度矩阵接近奇异矩阵(0刚度)，造成迭代不易收敛，需要更小的子步步长，而过小的步长会造成计算时间过长。本文算例在极限载荷阶段采用弧长法计算，最大弧长半径为0.15，最小弧长半径为0.03，收敛速度较快。图7为弧长法迭代原理图。

图7 弧长法迭代过程

2.4 数值模型计算的合理性验证

图8为极限载荷下管管接触的变形图。将混合单元建立模型的接触力数值解与小载荷状态下的解析解和梁单元状态下的数值解进行对比，得出表1数据，由此可以看出题模型建立的正确性。

图8 极限载荷下管束变形

N

由表可以看出，3种计算结果吻合较好，从而证明了本文模型建立的正确性。

3 横截面变形情况

3.1 卸载前后管横截面变形

在加载和卸载后，管束不同位置的横截面在不同大小的外载荷作用下和卸载后表现出不同变形情况。图9将横截面上下端点距离设为h，将δ=(D-h)/D定义为截面的变形率，用来衡量横截面的变形程度。

图9 横截面端点位置

由于本文例题中外载荷施加位置为接触管轴线中点处，故横截面变形成对称分布，取距轴线中点距离为0.00、0.03、0.06、0.09、0.12、0.15、0.18、0.21、0.24m的截面作为研究对象。

本组研究根据内固定稳定性提取了3组钢板模型在不同工况条件下的应变能指标，如表1所示。其中，在2工况条件下，从应变能和计算的2种刚度来看，FP整体刚度要高于其他2组模型。而RP在轴向压缩工况下，应变能（结构柔度指标）相比SP降低了21.4%，相比FP仅仅提高了7.2%；轴向刚度则比SP提高了21.29%，比FP刚度仅仅降低了6.8%。在扭转工况条件下，RPDE应变能相比SP降低了16.28%，相比FP则增加了13.5%；在扭转刚度方面，RP比SP提升了19.5%，而相比FP则仅下降了12.0%。由此可见，RP在两组工况条件下较之SP实现了固定刚度上的显著提升。

图10、11为不同载荷下同一位置横截面的变形率，可以看出随外载荷的增大，各位置横截面接触变形率与塑性变形率逐渐增大，且两者变化趋势相似，其中在外载荷为10～35kN过程中，随着载荷的增大横截面的变形率变化较大，之后横截面变形幅度趋于平缓。

图10 不同载荷下横截面的接触变形率

图11 不同载荷下横截面的塑性变形率

图12、13为不同位置横截面的变形率，可以看出接触变形率与塑性变形率两者变化趋势相似，距离中心越远的横截面变形率越小，且变形率与距中心位置距离几乎成线性。

图12 不同横截面的接触变形率

图13 不同横截面的塑性变形率

3.2 卸载前后管横截面变形率对比

在加载中接触变形较大的截面，卸载后会产生较大的塑性变形，图14、15给出了各截面在不同载荷作用下的接触变形率与塑性变形率的差值，即管横截面受压卸载后的残余变形。

图14 不同载荷下横截面的变形率差值

图15 不同载荷下横截面的变形率差值

可以看出，0.00、0.03、0.06、0.09m处截面变化趋势相似，可以看出变形率的差值随外载荷的增大，呈幅值逐渐减小的正弦式变化，且最大差值不超过3%。0.12、0.15、0.18、0.21m处截面变化趋势相似，可以看出变形率的差值随外载荷的增大，呈幅值逐渐减小的正弦式变化，相比于靠近中心位置的截面，变化周期变小。且最大差值不超过2%。0.24m靠近接触变形边缘处的截面，出现了塑性变形率较接触变形率大的现象。

3.3 卸载前后管横截面各点位移

以管中心位置为坐标原点，建立极坐标系。图16、17分别为2管0m位置横截面上的点，在不同大小外载荷作用下卸载前后的位移情况，可见在不同外载荷作用下点的位移趋势是相同的，由图16可以看出管发生接触变形时横截面下端点位移最小，上端点位移最大，0～75°左右各点位移逐渐减小，至75°左右处点位移达到极小值，75～105°左右点位移逐渐增大，至105°左右处点位移达到极大值，之后逐渐减小，至180°处达到最小值。由图17可以看出管发生塑性变形时横截面上端点位移最大，由于卸载后截面的回弹，管横截面下端点出现较大位移，0～75°左右各点位移逐渐减小，至75°左右处点位移达到极小值，75～90°左右点位移逐渐增大达到极大值，90～120°左右处点位移减小，横截面点位移最小处位于120°左右，120～180°点位移逐渐增大。

图16 0m处横截面接触变形点位移

图17 0m处截面塑性变形点位移

4 接触力与接触区域

4.1 接触力与最大接触应力变化情况

在加载过程中，随着外载荷的增大，管管间接触力会随之改变，当载荷达到一定值时，2管内壁将发生自接触。由图18、19可以看出管管间的接触力的合力与2管(目标管)自接触的接触力合力随外载荷的增大几乎成线性增长。管管间的最大接触应力的大小较为稳定，2管自接触的最大接触应力，增长幅度较大。最大接触应力发生的位置，随载荷增大而变化，图20中标出了最大应力的大致位置，在加载接触发生的初期位于接触区域中心处，随外载荷增大位于管轴向边缘附近，继续增大外载荷位于管圆周方向边缘处，在加载后期位于接触区域中心处。

图18 接触力合力变化

图19 最大接触应力变化

4.2 接触区域变化情况

如图20的接触力与接触区域所示，以a1代表接触区域长度，以a2代表分离区域长度，以b1代表接触区域宽度，以b2代表分离区域的宽度。则a1-a2为实际接触区域的长度，b1-b2为实际接触区域的宽度。

图20 接触力与接触区域

图21、22反映了接触区域、分离区域和实际接触区域长度、宽度的变化情况，接触区域长度随外载荷增大逐渐增大，但增大幅度逐渐减小。接触区域宽度随载荷增大逐渐增大，但增大幅度逐渐减小。

图21 区域长度

图22 区域宽度

5 结论

5.1笔者以管结构作为研究对象，采用数值方法对弹塑性管束接触变形的全过程进行模拟，这其中伴随着3种非线性效应。在加载中采用多个载荷步加载，在非线性程度较强段，采用小步长加载。在极限载荷段采用弧长法代替传统的牛顿拉普森迭代法。

5.2在管横截面变形的研究中得出，各位置横截面在不同载荷加载的接触变形率与卸载后的塑性变形率变化趋势相似，且各截面接触变形率与塑性变形率的差值呈幅值减小的正弦式变化。变形率差值最大不超过3%。

5.3在管管接触力和接触区域研究方面得出，管管间接触力与外载荷呈线性关系，最大接触应力变化不大。2管自接触接触力与外载荷呈线性关系，最大接触应力变化较大。接触区域长度与宽度随外载荷增大逐渐增大，但增大幅度逐渐减小。

[1] 刘巨保，王秀文，岳欠杯.接触油管屈曲变形与内外层杆管接触分析[J].力学与实践，2011，33(3)，25～29.

[2] 颜惠庚，郁翠菊.换热器液压胀管残余接触压力的工程图算法[J].化工机械，2001，28(4)：211～214.

[3] 于洪杰，钱才富.液压胀接接头密封性能的力学表征[J].化工机械，2010，37(6)：758～762.

[4] 刘冰，綦耀光，韩军.同心双管柱受力分析及变形计算[J].石油机械，2013，41(1)：64～67.

[5] 陈银强，桂春，王先元.外来物对蒸汽发生器传热管微动磨损的分析研究[J].核动力工程，2011，32(1)：21～25.

[6] 付茂青，刘巨保.管束横向接触与碰撞的计算方法研究[D].大庆：东北石油大学，2015.

[7] 董世民，张万胜，张红，等.定向井有杆抽油系统杆管分布接触压力的研究[J].工程力学，2011，28(10)：179～184.

[8] 剧锦三，杨蔚彪，蒋秀根.刚体撞击弹塑性直杆时冲击荷载之数值解[J].工程力学，2007，24(6)：49～53．

[9] 肖会芳, 邵毅敏, 周晓君. 非连续粗糙多界面接触变形和能量损耗特性研究[J].振动与冲击，2012，31(6)：83～89.

[10] 王晓军，何兆坤，苗承鹏.C10100铜合金管挤压成形过程的有限元模拟[J].锻压技术，2015，40(6)：144～149.

[11] 沈群，吴志林，袁人枢，等.镁合金扩管挤压过程的有限元数值模拟[J].热加工工艺，2013，40(7)：82～85.

[12] 刘巨保，罗敏.有限单元法及应用[M].北京：中国电力出版社，2013：181～192.

[13] 苏春峰，艾延廷，娄小宝，等.接触非线性仿真中接触刚度因子选取的方法研究[J].沈阳航空航天大学学报，2009，26(3)：5～9.

[14] 刘金梅，周国强，韩国有.弧长法在服役石油井架非线性全过程仿真中的应用研究[J].应用力学学报，2012，29(2)：229～233.

[15] Feng Y T，Peri D，Owen D R J.A New Criterion for Determination of Initial Loading Parameter in Arc-length Methods[J].Computers & Structures，1996，58(3)：479～485.

[16] 王瑁成.有限单元法及应用[M].北京：清华大学出版社，2003：556～558.

[17] 王新敏，李义强，许宏伟.结构分析单元与应用[M].北京：人民交通出版社，2011：498～502.