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

劉巨保 陳 健 姚利明

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

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

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

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

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

圖1所示為兩根平行管束，將上管束定義為1管，下管束定義為2管。1管為簡(jiǎn)支約束，2管兩端全約束，在1管軸線(xiàn)中點(diǎn)處施加一個(gè)集中載荷，使它變形并與兩管發(fā)生接觸，進(jìn)而擠壓兩管也產(chǎn)生變形。兩管長(zhǎng)度均為L(zhǎng)，間隙為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對(duì)1管撓度的貢獻(xiàn)的大小。聯(lián)立上述方程可以解得此時(shí)的接觸力：

(5)

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

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

2.1 數(shù)值模型

對(duì)于梁?jiǎn)卧挠邢拊治霾荒芸紤]到梁橫截面上、內(nèi)部各點(diǎn)的應(yīng)力位移，需要實(shí)體單元建立模型。但如果全部采用實(shí)體單元，由于模型幾何尺寸較大使得有限元模型的節(jié)點(diǎn)數(shù)過(guò)多，造成計(jì)算時(shí)間過(guò)長(zhǎng)。筆者首先使用梁?jiǎn)卧囁闵鲜隼}獲得大致的接觸位置，之后使用實(shí)體單元代替接觸管和目標(biāo)管接觸段的梁?jiǎn)卧?shí)體單元網(wǎng)格細(xì)化，其他位置仍采用梁?jiǎn)卧?，進(jìn)行混合建模。采用節(jié)點(diǎn)接觸綁定連接兩種單元，以保證兩種類(lèi)型單元在力和位移上的正確傳遞。

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

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

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

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

2.2 數(shù)值接觸算法

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

Π=Πu+ΠCL

(6)

其中Πu是原問(wèn)題中不包括接觸約束條件的總位能，ΠCL是不同算法引入約束條件的附加泛函。增廣拉格朗日乘子法中?。?/p>

(7)

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

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

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

2.3 非線(xiàn)性方程的收斂計(jì)算

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

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

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

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

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

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

2.3.2非線(xiàn)性設(shè)置

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

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

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

2.3.4弧長(zhǎng)法的使用

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

圖7 弧長(zhǎng)法迭代過(guò)程

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

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

圖8 極限載荷下管束變形

N

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

3 橫截面變形情況

3.1 卸載前后管橫截面變形

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

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

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

本組研究根據(jù)內(nèi)固定穩(wěn)定性提取了3組鋼板模型在不同工況條件下的應(yīng)變能指標(biāo)，如表1所示。其中，在2工況條件下，從應(yīng)變能和計(jì)算的2種剛度來(lái)看，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%。由此可見(jiàn)，RP在兩組工況條件下較之SP實(shí)現(xiàn)了固定剛度上的顯著提升。

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

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

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

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

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

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

3.2 卸載前后管橫截面變形率對(duì)比

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

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

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

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

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

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

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

圖18 接觸力合力變化

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

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

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

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

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

圖21 區(qū)域長(zhǎng)度

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

5 結(jié)論

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[14] 劉金梅，周?chē)?guó)強(qiáng)，韓國(guó)有.弧長(zhǎng)法在服役石油井架非線(xiàn)性全過(guò)程仿真中的應(yīng)用研究[J].應(yīng)用力學(xué)學(xué)報(bào)，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.