管管間接觸擠壓變形的數值模擬方法研究①

劉巨保 陳 健 姚利明

(東北石油大學機械科學與工程學院)

在石油石化工程和裝備中，存在一類管束結構，如井下油管與套管、管殼式換熱器中的換熱管、海洋鉆井隔水管等。這些管束結構在外力作用下，常發生接觸擠壓或振動碰撞而產生變形磨損甚至破壞，危及到管束的安全運行。此類問題屬于接觸變形問題，且往往伴隨材料非線性和(或)幾何非線性，難以得到準確的解析解[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.