劉巨保 陳 健 姚利明

(東北石油大學(xué)機(jī)械科學(xué)與工程學(xué)院)

在石油石化工程和裝備中，存在一類管束結(jié)構(gòu)，如井下油管與套管、管殼式換熱器中的換熱管、海洋鉆井隔水管等。這些管束結(jié)構(gòu)在外力作用下，常發(fā)生接觸擠壓或振動碰撞而產(chǎn)生變形磨損甚至破壞，危及到管束的安全運(yùn)行。此類問題屬于接觸變形問題，且往往伴隨材料非線性和(或)幾何非線性，難以得到準(zhǔn)確的解析解[1～5]。因此對管束結(jié)構(gòu)接觸變形的全過程進(jìn)行三重非線性仿真分析是非常必要的。

在以往的管束接觸和碰撞問題數(shù)值計算中，一般不考慮接觸的局部變形，采用橫截面不變形的梁單元建立數(shù)值分析模型[6，7]，由于模型的限制無法考慮局部變形的問題。而對于采用實(shí)體單元建模的接觸碰撞問題研究中，更多的是關(guān)注接觸力、沖擊載荷和能量損失，而未考察結(jié)構(gòu)自身的具體變形規(guī)律情況[8，9]。在材料成型領(lǐng)域，多采用無間隙接觸的擠壓工藝對產(chǎn)品進(jìn)行成型[10，11]，這在數(shù)值模擬計算的收斂上降低了計算難度。

筆者以平行管束為研究對象使用實(shí)體單元與梁單元結(jié)合建立數(shù)值模型，采用增廣拉格朗日算法，引入了管擠壓變形產(chǎn)生的大位移幾何非線性和材料非線性效應(yīng)，對有間隙的管管接觸變形進(jìn)行非線性仿真分析。

1 管管間接觸擠壓變形分析力學(xué)模型

圖1所示為兩根平行管束，將上管束定義為1管，下管束定義為2管。1管為簡支約束，2管兩端全約束，在1管軸線中點(diǎn)處施加一個集中載荷，使它變形并與兩管發(fā)生接觸，進(jìn)而擠壓兩管也產(chǎn)生變形。兩管長度均為L，間隙為g，兩管外徑相同均為D，內(nèi)徑不同，1管內(nèi)徑為d1，2管內(nèi)徑為d2，1管材料彈性模量為E，慣性矩為I。圖2為管管接觸的3種狀態(tài)，在第2種管管剛接觸的狀態(tài)下，由于接觸區(qū)域較小，可假設(shè)為點(diǎn)接觸，且接觸位置位于L/2處。

圖1 管管間接觸擠壓力學(xué)模型

圖2 管管接觸的3種狀態(tài)

由材料力學(xué)疊加原理建立平衡方程為：

ωF-ωC=g

(1)

(2)

(3)

(4)

其中ωF，ωC分別為外載荷F與接觸力Fc對1管撓度的貢獻(xiàn)的大小。聯(lián)立上述方程可以解得此時的接觸力：

(5)

繼續(xù)增大外載荷接觸區(qū)域也會增大，出現(xiàn)第3種狀態(tài)，2管橫截面發(fā)生大變形，此時上述假設(shè)不再成立，所得接觸力計算結(jié)果也不具參考性，應(yīng)采用數(shù)值分析方法進(jìn)行計算。

2 管管間接觸擠壓變形分析的數(shù)值方法

2.1 數(shù)值模型

對于梁單元的有限元分析不能考慮到梁橫截面上、內(nèi)部各點(diǎn)的應(yīng)力位移，需要實(shí)體單元建立模型。但如果全部采用實(shí)體單元，由于模型幾何尺寸較大使得有限元模型的節(jié)點(diǎn)數(shù)過多，造成計算時間過長。筆者首先使用梁單元試算上述例題獲得大致的接觸位置，之后使用實(shí)體單元代替接觸管和目標(biāo)管接觸段的梁單元并將實(shí)體單元網(wǎng)格細(xì)化，其他位置仍采用梁單元，進(jìn)行混合建模。采用節(jié)點(diǎn)接觸綁定連接兩種單元，以保證兩種類型單元在力和位移上的正確傳遞。

圖3為數(shù)值模型，將實(shí)體單元段的接觸管與目標(biāo)管軸線方向的單元長度設(shè)置為3mm，將圓周分成18份，將目標(biāo)管橫截面分成3份，接觸管分成1份。圖4為模型中的接觸對，將2管與下方平面設(shè)置為無間隙接觸，將1管與2管的外表面設(shè)置為含間隙接觸，由于擠壓變形將造成2管內(nèi)壁發(fā)生自接觸，應(yīng)對2管內(nèi)壁上下設(shè)置為含間隙接觸。

圖3 數(shù)值模型示意圖

圖4 數(shù)值模型接觸對示意圖

其中接觸管(忽略其橫截面變形)材料定義為彈性模量E1=200GPa的彈性材料，將目標(biāo)管材料定義為彈性模量E2=71GPa、屈服強(qiáng)度為280MPa、切線模量為500MPa的雙線性彈塑性材料。應(yīng)用Von Mises屈服準(zhǔn)則和隨動強(qiáng)化模型。通過多個載荷步加載，一個載荷步卸載。

2.2 數(shù)值接觸算法

對于接觸算法即將有限元法的約束變分原理應(yīng)用于接觸問題的求解，即分別采用拉格朗日乘子法、罰函數(shù)法和增廣拉格朗日乘子法。對于包含接觸界面的接觸問題，泛函可以表示為：

Π=Πu+ΠCL

(6)

其中Πu是原問題中不包括接觸約束條件的總位能，ΠCL是不同算法引入約束條件的附加泛函。增廣拉格朗日乘子法中取：

(7)

其中Λ′=diag(a1′,a2′,a3′)；a1′、a2′和a3′為罰系數(shù)，為設(shè)定常數(shù)；λ′={λ1′λ2′λ3′}T，g′={g1′g2′g3′}T。

增廣拉格朗日乘子法是為了找到精確的拉格朗日乘子(即接觸力)，而對罰函數(shù)法進(jìn)行一系列修正迭代。與罰函數(shù)法相比，增廣拉格朗日乘子法的優(yōu)點(diǎn)是容易得到良態(tài)條件，對接觸剛度的敏感性較小。缺點(diǎn)是，可能需要更多的迭代，特別是變形后網(wǎng)格變得太扭曲。

由于筆者研究的管管接觸不會出現(xiàn)扭曲畸形網(wǎng)格，且該問題經(jīng)常因?yàn)榻佑|界面耦合關(guān)系限制罰參數(shù)的取值，容易引起振蕩。所以筆者采用對接觸剛度敏感性小的增廣拉格朗日乘子法進(jìn)行有限元管管接觸問題的計算方法研究[10]。

2.3 非線性方程的收斂計算

由于本例題考慮了3種非線性效應(yīng)，除網(wǎng)格的劃分、接觸對的正確建立外，難點(diǎn)在于非線性方程組的收斂，如何快速并準(zhǔn)確地求解非線性方程成為了本題的關(guān)鍵。

由于非線性問題需要采用迭代逼近技術(shù)進(jìn)行求解，通常采用牛頓-拉普森法計算，其收斂速度較快，但和載荷子步、非線性設(shè)置有關(guān)。除此之外還需要對載荷子步的收斂判定方法和準(zhǔn)則進(jìn)行定義。圖5給出了靜力學(xué)非線性方程求解應(yīng)設(shè)置的選項(xiàng)。

圖5 靜力學(xué)非線性方程求解設(shè)置選項(xiàng)

2.3.1載荷子步相關(guān)設(shè)置

由牛頓拉普森迭代法通過計算各載荷(位移)增量的結(jié)果，進(jìn)而得出極限載荷下的計算結(jié)果，步長的設(shè)置對計算能否收斂影響較大。圖6是通過試算得到的載荷作用點(diǎn)的力位移曲線，第1個載荷步為0～20kN，最小時間增量為0.001，最大時間增量為0.005。第1個載荷步為20～40kN，最小時間增量為0.005，最大時間增量為0.050。并將自動時間步設(shè)置為程序自行選擇。線性方程求解器設(shè)置為程序自行選擇。對于收斂的判定，例題采用力收斂準(zhǔn)則，殘差的計算方法采用二范數(shù)，收斂容差為0.001。如果子步達(dá)到某一物理限制，程序?qū)⒁暣瞬介L過大，進(jìn)而自動將子步二分，例題將這一物理限制定義為子步最大迭代次數(shù)設(shè)置為25次。由于例題中結(jié)構(gòu)產(chǎn)生了塑性應(yīng)變，為了更好地控制時間步長上的縮減，增加另一個物理限制即最大塑性應(yīng)變增量限制，將其值設(shè)置為0.1。

圖6 載荷作用點(diǎn)力位移曲線

2.3.2非線性設(shè)置

ANSYS提供了4種切向剛度矩陣修改方法，并提供了一系列的非線性設(shè)置來保證非線性求解的收斂。例題中對于各子步切向剛度矩陣的修改采用程序自行選擇。非線性設(shè)置選項(xiàng)包括自由度預(yù)測、自適應(yīng)下降和線性搜索。其中自由度預(yù)測適用于無位移突變的非線性分析，可加速收斂。由圖6可以看出力位移曲線非線性響應(yīng)相對平滑，故例題中激活自由度預(yù)測。自適應(yīng)下降與線性搜索均采用程序自行選擇。

2.3.3實(shí)常數(shù)的設(shè)置

文中的算例將接觸剛度系數(shù)設(shè)置為0.01，穿透容差系數(shù)設(shè)置為0.1，收斂速度較快且得到的計算結(jié)果精確度較高。

2.3.4弧長法的使用

由圖6可以看出在極限載荷段，載荷變化較大但位移變化較小，出現(xiàn)力位移曲線趨于水平的現(xiàn)象，這可能造成切向剛度矩陣接近奇異矩陣(0剛度)，造成迭代不易收斂，需要更小的子步步長，而過小的步長會造成計算時間過長。本文算例在極限載荷階段采用弧長法計算，最大弧長半徑為0.15，最小弧長半徑為0.03，收斂速度較快。圖7為弧長法迭代原理圖。

圖7 弧長法迭代過程

2.4 數(shù)值模型計算的合理性驗(yàn)證

圖8為極限載荷下管管接觸的變形圖。將混合單元建立模型的接觸力數(shù)值解與小載荷狀態(tài)下的解析解和梁單元狀態(tài)下的數(shù)值解進(jìn)行對比，得出表1數(shù)據(jù)，由此可以看出題模型建立的正確性。

圖8 極限載荷下管束變形

N

由表可以看出，3種計算結(jié)果吻合較好，從而證明了本文模型建立的正確性。

3 橫截面變形情況

3.1 卸載前后管橫截面變形

在加載和卸載后，管束不同位置的橫截面在不同大小的外載荷作用下和卸載后表現(xiàn)出不同變形情況。圖9將橫截面上下端點(diǎn)距離設(shè)為h，將δ=(D-h)/D定義為截面的變形率，用來衡量橫截面的變形程度。

圖9 橫截面端點(diǎn)位置

由于本文例題中外載荷施加位置為接觸管軸線中點(diǎn)處，故橫截面變形成對稱分布，取距軸線中點(diǎn)距離為0.00、0.03、0.06、0.09、0.12、0.15、0.18、0.21、0.24m的截面作為研究對象。

本組研究根據(jù)內(nèi)固定穩(wěn)定性提取了3組鋼板模型在不同工況條件下的應(yīng)變能指標(biāo)，如表1所示。其中，在2工況條件下，從應(yīng)變能和計算的2種剛度來看，F(xiàn)P整體剛度要高于其他2組模型。而RP在軸向壓縮工況下，應(yīng)變能（結(jié)構(gòu)柔度指標(biāo)）相比SP降低了21.4%，相比FP僅僅提高了7.2%；軸向剛度則比SP提高了21.29%，比FP剛度僅僅降低了6.8%。在扭轉(zhuǎn)工況條件下，RPDE應(yīng)變能相比SP降低了16.28%，相比FP則增加了13.5%；在扭轉(zhuǎn)剛度方面，RP比SP提升了19.5%，而相比FP則僅下降了12.0%。由此可見，RP在兩組工況條件下較之SP實(shí)現(xiàn)了固定剛度上的顯著提升。

圖10、11為不同載荷下同一位置橫截面的變形率，可以看出隨外載荷的增大，各位置橫截面接觸變形率與塑性變形率逐漸增大，且兩者變化趨勢相似，其中在外載荷為10～35kN過程中，隨著載荷的增大橫截面的變形率變化較大，之后橫截面變形幅度趨于平緩。

圖10 不同載荷下橫截面的接觸變形率

圖11 不同載荷下橫截面的塑性變形率

圖12、13為不同位置橫截面的變形率，可以看出接觸變形率與塑性變形率兩者變化趨勢相似，距離中心越遠(yuǎn)的橫截面變形率越小，且變形率與距中心位置距離幾乎成線性。

圖12 不同橫截面的接觸變形率

圖13 不同橫截面的塑性變形率

3.2 卸載前后管橫截面變形率對比

在加載中接觸變形較大的截面，卸載后會產(chǎn)生較大的塑性變形，圖14、15給出了各截面在不同載荷作用下的接觸變形率與塑性變形率的差值，即管橫截面受壓卸載后的殘余變形。

圖14 不同載荷下橫截面的變形率差值

圖15 不同載荷下橫截面的變形率差值

可以看出，0.00、0.03、0.06、0.09m處截面變化趨勢相似，可以看出變形率的差值隨外載荷的增大，呈幅值逐漸減小的正弦式變化，且最大差值不超過3%。0.12、0.15、0.18、0.21m處截面變化趨勢相似，可以看出變形率的差值隨外載荷的增大，呈幅值逐漸減小的正弦式變化，相比于靠近中心位置的截面，變化周期變小。且最大差值不超過2%。0.24m靠近接觸變形邊緣處的截面，出現(xiàn)了塑性變形率較接觸變形率大的現(xiàn)象。

3.3 卸載前后管橫截面各點(diǎn)位移

以管中心位置為坐標(biāo)原點(diǎn)，建立極坐標(biāo)系。圖16、17分別為2管0m位置橫截面上的點(diǎn)，在不同大小外載荷作用下卸載前后的位移情況，可見在不同外載荷作用下點(diǎn)的位移趨勢是相同的，由圖16可以看出管發(fā)生接觸變形時橫截面下端點(diǎn)位移最小，上端點(diǎn)位移最大，0～75°左右各點(diǎn)位移逐漸減小，至75°左右處點(diǎn)位移達(dá)到極小值，75～105°左右點(diǎn)位移逐漸增大，至105°左右處點(diǎn)位移達(dá)到極大值，之后逐漸減小，至180°處達(dá)到最小值。由圖17可以看出管發(fā)生塑性變形時橫截面上端點(diǎn)位移最大，由于卸載后截面的回彈，管橫截面下端點(diǎn)出現(xiàn)較大位移，0～75°左右各點(diǎn)位移逐漸減小，至75°左右處點(diǎn)位移達(dá)到極小值，75～90°左右點(diǎn)位移逐漸增大達(dá)到極大值，90～120°左右處點(diǎn)位移減小，橫截面點(diǎn)位移最小處位于120°左右，120～180°點(diǎn)位移逐漸增大。

圖16 0m處橫截面接觸變形點(diǎn)位移

圖17 0m處截面塑性變形點(diǎn)位移

4 接觸力與接觸區(qū)域

4.1 接觸力與最大接觸應(yīng)力變化情況

在加載過程中，隨著外載荷的增大，管管間接觸力會隨之改變，當(dāng)載荷達(dá)到一定值時，2管內(nèi)壁將發(fā)生自接觸。由圖18、19可以看出管管間的接觸力的合力與2管(目標(biāo)管)自接觸的接觸力合力隨外載荷的增大幾乎成線性增長。管管間的最大接觸應(yīng)力的大小較為穩(wěn)定，2管自接觸的最大接觸應(yīng)力，增長幅度較大。最大接觸應(yīng)力發(fā)生的位置，隨載荷增大而變化，圖20中標(biāo)出了最大應(yīng)力的大致位置，在加載接觸發(fā)生的初期位于接觸區(qū)域中心處，隨外載荷增大位于管軸向邊緣附近，繼續(xù)增大外載荷位于管圓周方向邊緣處，在加載后期位于接觸區(qū)域中心處。

圖18 接觸力合力變化

圖19 最大接觸應(yīng)力變化

4.2 接觸區(qū)域變化情況

如圖20的接觸力與接觸區(qū)域所示，以a1代表接觸區(qū)域長度，以a2代表分離區(qū)域長度，以b1代表接觸區(qū)域?qū)挾龋詁2代表分離區(qū)域的寬度。則a1-a2為實(shí)際接觸區(qū)域的長度，b1-b2為實(shí)際接觸區(qū)域的寬度。

圖20 接觸力與接觸區(qū)域

圖21、22反映了接觸區(qū)域、分離區(qū)域和實(shí)際接觸區(qū)域長度、寬度的變化情況，接觸區(qū)域長度隨外載荷增大逐漸增大，但增大幅度逐漸減小。接觸區(qū)域?qū)挾入S載荷增大逐漸增大，但增大幅度逐漸減小。

圖21 區(qū)域長度

圖22 區(qū)域?qū)挾?/p>

5 結(jié)論

5.1筆者以管結(jié)構(gòu)作為研究對象，采用數(shù)值方法對彈塑性管束接觸變形的全過程進(jìn)行模擬，這其中伴隨著3種非線性效應(yīng)。在加載中采用多個載荷步加載，在非線性程度較強(qiáng)段，采用小步長加載。在極限載荷段采用弧長法代替?zhèn)鹘y(tǒng)的牛頓拉普森迭代法。

5.2在管橫截面變形的研究中得出，各位置橫截面在不同載荷加載的接觸變形率與卸載后的塑性變形率變化趨勢相似，且各截面接觸變形率與塑性變形率的差值呈幅值減小的正弦式變化。變形率差值最大不超過3%。

5.3在管管接觸力和接觸區(qū)域研究方面得出，管管間接觸力與外載荷呈線性關(guān)系，最大接觸應(yīng)力變化不大。2管自接觸接觸力與外載荷呈線性關(guān)系，最大接觸應(yīng)力變化較大。接觸區(qū)域長度與寬度隨外載荷增大逐漸增大，但增大幅度逐漸減小。

[1] 劉巨保，王秀文，岳欠杯.接觸油管屈曲變形與內(nèi)外層桿管接觸分析[J].力學(xué)與實(shí)踐，2011，33(3)，25～29.

[2] 顏惠庚，郁翠菊.換熱器液壓脹管殘余接觸壓力的工程圖算法[J].化工機(jī)械，2001，28(4)：211～214.

[3] 于洪杰，錢才富.液壓脹接接頭密封性能的力學(xué)表征[J].化工機(jī)械，2010，37(6)：758～762.

[4] 劉冰，綦耀光，韓軍.同心雙管柱受力分析及變形計算[J].石油機(jī)械，2013，41(1)：64～67.

[5] 陳銀強(qiáng)，桂春，王先元.外來物對蒸汽發(fā)生器傳熱管微動磨損的分析研究[J].核動力工程，2011，32(1)：21～25.

[6] 付茂青，劉巨保.管束橫向接觸與碰撞的計算方法研究[D].大慶：東北石油大學(xué)，2015.

[7] 董世民，張萬勝，張紅，等.定向井有桿抽油系統(tǒng)桿管分布接觸壓力的研究[J].工程力學(xué)，2011，28(10)：179～184.

[8] 劇錦三，楊蔚彪，蔣秀根.剛體撞擊彈塑性直桿時沖擊荷載之?dāng)?shù)值解[J].工程力學(xué)，2007，24(6)：49～53．

[9] 肖會芳, 邵毅敏, 周曉君. 非連續(xù)粗糙多界面接觸變形和能量損耗特性研究[J].振動與沖擊，2012，31(6)：83～89.

[10] 王曉軍，何兆坤，苗承鵬.C10100銅合金管擠壓成形過程的有限元模擬[J].鍛壓技術(shù)，2015，40(6)：144～149.

[11] 沈群，吳志林，袁人樞，等.鎂合金擴(kuò)管擠壓過程的有限元數(shù)值模擬[J].熱加工工藝，2013，40(7)：82～85.

[12] 劉巨保，羅敏.有限單元法及應(yīng)用[M].北京：中國電力出版社，2013：181～192.

[13] 蘇春峰，艾延廷，婁小寶，等.接觸非線性仿真中接觸剛度因子選取的方法研究[J].沈陽航空航天大學(xué)學(xué)報，2009，26(3)：5～9.

[14] 劉金梅，周國強(qiáng)，韓國有.弧長法在服役石油井架非線性全過程仿真中的應(yīng)用研究[J].應(yīng)用力學(xué)學(xué)報，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] 王瑁成.有限單元法及應(yīng)用[M].北京：清華大學(xué)出版社，2003：556～558.

[17] 王新敏，李義強(qiáng)，許宏偉.結(jié)構(gòu)分析單元與應(yīng)用[M].北京：人民交通出版社，2011：498～502.