張 杰，梁博豐，潘 波，沈坤蓉，陳林燕

(1.西南石油大學 機電工程學院，成都 610500;2.西南石油大學 南充校區(qū)，四川南充 637001)

0 引言

隨著國家能源結(jié)構(gòu)調(diào)整，天然氣輸送管道建設得到了迅猛發(fā)展。由于中國地質(zhì)條件復雜，且管網(wǎng)規(guī)?；?、密集化等特點日益突出，長輸天然氣管道不可避免地穿越山地、丘陵、沙漠、平原、水網(wǎng)等地帶，容易受多種地質(zhì)災害影響，如滑坡、崩塌、水毀、沉降等，其中滑坡是影響管道安全的常見地質(zhì)災害之一。滑坡作用下天然氣管道易發(fā)生破壞，造成天然氣泄漏，進而引發(fā)火災、燃爆等事故，甚至造成災難性后果[1-2]。因此，研究滑坡作用下的管道力學響應對其安全評價、維修防護至關重要。

根據(jù)滑坡與管道的位置關系，可將其分為橫向滑坡和縱向滑坡，其中橫向滑坡對管道危害最為嚴重，圖1為管道穿越橫向滑坡區(qū)域示意圖。當橫向滑坡作用時，埋地管道將受到滑坡側(cè)向力，其與管道鋪設后的形態(tài)有著密切關系[3]。按引起滑動的力學性質(zhì)可以將橫向滑坡分為牽引式和推移式兩類：牽引式滑坡對管道部無被動土壓力，管道呈懸空狀態(tài)；而推移式滑坡在發(fā)生整體滑動的情況下，下部巖土層對上部巖土層有一定的支撐作用[4]。

圖1 管道穿越橫向滑坡區(qū)示意

目前，對管道力學響應問題的研究主要是基于彈性地基梁模型開展的。如張土喬等[5]建立了垂直載荷作用下的管道縱向受力簡化模型，研究了影響管道危險區(qū)應力的主要因素；鄧道明等[6]分析了埋地管道在無內(nèi)壓和有內(nèi)壓作用下的內(nèi)力和變形，推導出了管環(huán)彎矩和撓度的一系列計算式；袁祥忠等[7]建立了牽引式滑坡作用下管道橫向受力的分段模型，采用分段模型，利用微元思想，計算出各管段的解析表達式；Zhang等[8]對滑坡區(qū)管道力學行為進行了數(shù)值仿真，探究多種參數(shù)對其力學影響，較為真實地反映出管道在牽引式滑坡作用下的形變及受力情況；Li等[3]建立了海底滑坡下推移式管道受力的分段模型，將海底管道分為4段，列出相應平衡方程，應用迭代法得出了各段解析表達式，但該模型待定系數(shù)過多，邊界條件及連續(xù)性條件復雜，求解過程較繁瑣。此外，大部分學者對處于埋設狀態(tài)管段未作較詳細分析，將懸跨段視為簡支梁，忽略了非滑坡段管道變形對滑坡段管道變形的影響，使得結(jié)果誤差偏大。

為此，本文提出一套改進的解析模型來描述牽引式滑坡下天然氣管道力學響應問題。參照懸鏈線理論和大變形梁理論[9-10]，得出管道在牽引式滑坡作用下的最終狀態(tài)。根據(jù)邊界條件，提出一種迭代的方法計算待定系數(shù)，最終得到相應的解析表達式，并在此基礎上進行參數(shù)分析，以指導工程實踐。

1 管道力學模型

由彈性地基梁理論[11]可知，當管道發(fā)生彎曲時，管道下方的土體也會發(fā)生微小形變，以抵抗外載作用，如圖2所示。由懸鏈線特征可知，管道與土體在點A1處分離，其邊界A1處的曲率與A1點以后有限范圍(A1A2)內(nèi)管道曲線曲率的符號相同，均為正號，即A0A1A2段呈現(xiàn)出凸曲線。A2點以后，由于滑坡作用，根據(jù)大變形梁理論可知管道呈現(xiàn)出凹曲線，A2點是管道曲線曲率由正變負的反曲點，該點曲率為零。管道在A3點處的側(cè)向位移最大，曲線斜率為零。假設管道在幾何及受力上都對稱于y軸。

圖2 滑坡下管道受力變形示意

對于管道下面土體凹陷形成的懸空段，其橫向力q按式(1)計算[12]，不考慮管道內(nèi)介質(zhì)質(zhì)量。

(1)

式中D——管道外徑，m；

γt——覆土容重，kg/m3；

h——覆土高度(管道中心距填土表面的深度)，m；

φ——土壤內(nèi)的摩擦角,(°)。

為簡化假設，取橫向滑坡力q為常數(shù)。

管道側(cè)向運動時，土體響應采用理想線彈塑性模型[13]，由此可得圖3所示的雙線性p-y關系曲線。土體抗力p以k的增長率線性增長，直到達到最大值，當其達到其最大抗力p0.1時，對應的管道位移為0.1D[14]，所以在0.1D以內(nèi)的土體抗力增長率為k=10p0.1/D。

圖3 載荷位移關系曲線

2 控制方程

由于模型的對稱性，選取管道右半段作為研究對象。該管段根據(jù)其受力情況可分為3部分，如圖4所示。

圖4 滑坡下管道受力變形力學模型

2.1 管段Ⅰ

管段Ⅰ處于滑床段，自管道側(cè)向位移最大處A3點至滑床結(jié)束點A5。該管段側(cè)向位移最大，對整體管道的安全非常重要。根據(jù)彈性地基梁理論[15-16]，通過y軸方向平衡可得到控制方程：

(2)

式中E——管道彈性模量；

I——管道慣性矩；

yⅠ——管道從A3到A5的側(cè)向位移；

T——未知拉力；

w——管道質(zhì)量;

x——A4點的x軸坐標。

由式(2)可得：

(3)

式中C1,C2,C3,C4為未知系數(shù)，其中:

(4)

從而可得彎矩：

MⅠ=-EIy″Ⅰ

(5)

然后得到彈性范圍內(nèi)由管道軸力和彎矩造成的管道應變：

(6)

2.2 管段Ⅱ

該管段側(cè)向土體抗力隨著管道位移而發(fā)生線性變化，其控制方程為:

(7)

式(7)的解為：

+e-βx×(C7cosγx+C8sinγx)

(x2≤x≤x3)

(8)

式中C5,C6,C7,C8為未知參數(shù)，且:

(9)

(10)

同理,可得管段Ⅱ的彎矩和應變：

MⅡ=-EIy″Ⅱ

(11)

(12)

需指出，當x≥x3時，管道側(cè)向位移會非常小，利用A6處邊界條件可知，C5和C6都非常小，為簡化計算(此處對計算精度的影響極小，可以忽略不計)，令C5和C6均為零。

2.3 管段Ⅲ

該管段側(cè)向位移非常小，可以忽略。因此，可認為該管段僅受軸向土體抗力作用。軸向土體抗力可以通過將管道法向接觸力wn和摩擦系數(shù)相乘得到。摩擦系數(shù)與土體性質(zhì)密切相關。Bruton等[14]指出摩擦系數(shù)介于0.15～1.5之間。

遠離滑坡區(qū)管道的軸向拉力逐漸減小。這部分管道末端的軸向拉力為零。其長度可以通過A6點拉力和單位管道長度受到的土體抗力得到。

3 力學模型求解

3.1 待定系數(shù)分析

所建立力學模型中共有C1，C2，C3，C4，C7，C8，T，x1，x2，x3等10個未知量。

由A5點連續(xù)性條件可知：

yⅠ∣x=x2=0.1D

(13)

yⅡ∣x=x2=0.1D

(14)

y′Ⅰ∣x=x2=y′Ⅱ∣x=x2

(15)

y″Ⅰ∣x=x2=y″Ⅱ∣x=x2

(16)

y?Ⅰ∣x=x2=y?Ⅱ∣x=x2

(17)

由A3及A4點的邊界條件可知：

y′Ⅰ∣x=0=0

(18)

y″Ⅰ∣x=x1=0

(19)

y?Ⅰ∣x=0=0

(20)

其中:

(21)

(22)

y?Ⅰ=α3(C3eαx+C4e-αx)

(23)

y′Ⅱ=e-βx[C8(γcosγx-βsinγx)-C7(γsinγx

+βcosγx)]

(24)

y″Ⅱ=e-βx{C8[(β2-γ2)sinγx-2βγcosγx]

+C7[(β2-γ2)cosγx+2βγsinγx]}

(25)

y?Ⅱ=e-βx{C8[β(3γ2-β2)sinγx+γ(3β2-γ2)

×cosγx]+C7[-γ(3β2-γ2]sinγx

+β(3γ2-β2)cosγx)]}

(26)

當滑坡作用在整個懸跨段時，管道在拉力作用下伸長量和管道幾何變形Δl1相等。Δl1可以通過幾何積分得到：

(27)

其中，x3是管道側(cè)向位移非常小的位置。x3右側(cè)管道軸向拉力T在軸向摩擦力f作用下逐漸減小，所以管道伸長量可以通過下式[17]得到：

(28)

式中A——管道截面面積。

從而有：

(29)

3.2 待定系數(shù)求解

由式(13)～(20)，(29)組成方程組，包含9個方程，10個未知數(shù)，且方程形式復雜。因此，應用傳統(tǒng)方法無法求解。根據(jù)其中參數(shù)特征，采用如下計算方法。

由式(10)可知，若要使γ存在，則有：

(30)

(31)

當x2到其最大值時，管道從O到x2滑坡力等于從O到x3土體抗力[4]，由此可得：

qx2=p(x3-x2)

(32)

解得:

(33)

將式(21)～(26)代入式(13)～(20)得到以下方程式(考慮到滑坡力相對于管道質(zhì)量很大，此處省略管道質(zhì)量)。

cos(γx2)C7+sin(γx2)C8=0.1Deβx2

(34)

(35)

C2+αeαx2C3+αe-αx2C4+e-βx2(γsinγx2

(36)

α2eαx2C3-α2e-αx2C4-e-βx2[(β2-γ2)cosγx2

+2βγsinγx2]C7-e-βx2[(β2-γ2)sinγx2-2βγcosγx2]

(37)

C2+αC3+αC4=0

(38)

e-βx2{C8[β(3γ2-β2)sinγx2+γ(3β2-γ2)cosγx2]

+C7[-γ(3β2-γ2)sinγx2+β(3γ2-β2)cosγx2]}

=α3(C3eαx2+C4e-αx2)

(39)

(40)

C3+C4=0

(41)

觀察式(34)～(41)，可以將C1，C2，C3，C4，C7，C8，x1用x2,T和x3表示。其中x2為給定的滑坡寬度，同時假定一個x3和T，將其代入式(34),(35),(36),(37),(38),(40),(41)中，然后不斷修正x3和T，首先需要使其近似滿足式(17)的連續(xù)性條件，即∣y?Ⅰ(x2)-y?Ⅱ(x2)∣<δ1，其中δ1是一個定義的小量，為得到較高的精度，這里推薦1×10-6。由前文分析可知，將x3和T得初值設為x3max和T的最大值(由式(31)得到)。然后，將計算得到的x3和T以及其他參數(shù)代入到式(42)計算每個迭代過程中的軸向摩擦力f。

由式(29)得：

-2EAx3-2Tx3]

(42)

每次計算得到的摩擦力都需驗證，以滿足不等式∣f-f0∣<δ2。其中，δ2是一個定義的小量，本文推薦δ2=0.1 N/m。如果達不到δ1的精度，就需要首先對T進行修正，重新將x3和新的T代入求解；如果摩擦力的不等式不滿足，就需要對軸向拉力x3進行修正，再次將新的x3和上一輪的T代入求解。當不等式滿足了，最后得到的軸向拉力及其他參量即可達到標準，能夠進行結(jié)果輸出。需要說明的是，A6是管道側(cè)向位移較小的位置，可取x3為較大值，以保證足夠的精度。上述算法的流程如圖5[18]所示。

圖5 算法流程

4 算例分析

為開展案例分析，表1列出了管道力學參數(shù),表2列出了解析模型算法中需要輸入的初值及誤差值。

表1 管道力學參數(shù)

表2 解析模型待定系數(shù)求解參數(shù)

使用MATLAB編程，由于方程組系數(shù)矩陣存在數(shù)量級差距較大的元素，本文對系數(shù)矩陣進行截斷奇異值處理，可求得管段Ⅰ和管段Ⅱ的解析表達式[19]如下：

yⅠ=-5×10-5x2+0.7469+5.56×10-9

×(e0.13x+e-0.13x)

(43)

yⅡ=103×e-0.0845x[8.0173cos(0.014x)

+2.403sin(0.014x)]

(44)

4.1 滑坡寬度

圖6示出滑坡力為700 N/m時不同滑坡寬度的管道位移曲線。當滑坡寬度為140 m時，管道最大位移為1.006 8 m。隨著滑坡位移增加，管道最大位移也增加。當滑坡寬度達到240 m時，最大位移為2.022 4 m。

當滑坡力為700 N/m時，不同滑坡寬度下管道總應力分布如圖7所示。在彈性范圍內(nèi)，管道最大總應力隨著滑坡寬度的增加而增加。當滑坡寬度為180 m時，最大應力為256.21 MPa；當滑坡寬度為240 m時，最大應力為290.9 MPa，相差9.69%。由圖可知，應力曲線的峰值是尖銳的，沒有平滑過渡區(qū)，即管道中沒有塑性變形，原因是其最大應力為290.9 MPa，遠低于X80的屈服極限555 MPa；管道在滑坡出應力突變，從計算模型的角度考慮，這是由于滑坡力在滑坡寬度處從700 N/m突變?yōu)?造成的。

圖6 不同滑坡寬度下的管道位移曲線

圖7 不同滑坡寬度下的總應力分布

4.2 滑坡力

當滑坡寬度為180 m時，不同滑坡力作用下管道位移如圖8所示。隨著滑坡力增加，管道最大位移不斷增大。當滑坡力為500 N/m時，管道最大位移為1.195 5 m；當滑坡力為1 000 N/m時，管道最大位移為2.299 m。最大位移增加量較為接近，即在500～1 000 N/m的范圍內(nèi)，管道最大位移隨著滑坡力的增加而呈近似線性增加。

圖9示出當滑坡寬度為180 m時，管道最大彎曲應變隨著滑坡力的增大而增大的變化曲線。當滑坡力為500 N/m時，管道最大彎曲應變?yōu)?.257 7×10-4；當滑坡力為1 000 N/m時，管道最大彎曲應變?yōu)?.231 1×10-4。不同滑坡力下管道最大彎曲應變發(fā)生在滑坡寬度約180 m處，因而在滑坡末端是發(fā)生最大變形的地方，其原因是該處管道變形最大，會出現(xiàn)應力集中現(xiàn)象。因此，對于易發(fā)生滑坡區(qū)域管道，應在滑床段與非滑床段接觸區(qū)域處采取保護措施。

圖8 不同滑坡力下的管道位移曲線

圖9 不同滑坡力下的管道應變變化曲線

4.3 管道壁厚

管道內(nèi)徑為0.57 m，通過外徑變化來改變管道壁厚，對壁厚為0.02，0.035，0.05，0.065，0.08，0.095 m的管道進行研究。當滑坡力為700 N/m，滑坡寬度為180 m時，不同管道壁厚對管道位移影響如圖10所示，可以看出，管道最大位移隨著管道壁厚增加而減小。當管道外徑為0.61 m時，管道最大位移為2.299 m;當管道外徑為0.76 m時，管道最大位移為1.195 5 m。

隨著管道壁厚增加，管道最大位移不會無限減小。在方程有解的前提下，即在一定管道壁厚范圍內(nèi)，當管道壁厚增加到一定值時，管道最大位移減小量逐漸降低，最大位移逐漸接近一定值。圖11示出由最大位移組成的曲線呈微弱凹曲的形態(tài)，這是由于隨著壁厚增加，管道剛度增加，其更容易抵抗外力并減少自身變形造成的。

圖10 不同外徑管道位移曲線

圖11 不同外徑管道最大位移變化曲線

當滑坡力為700 N/m，滑坡寬度為180 m時，不同外徑管道彎曲應變分布如圖12所示。

圖12 不同管道外直徑下的管道彎曲應變曲線

由圖12可看出，隨著管道厚度增加，管道最大彎曲應變呈下降趨勢。當管道外徑為0.61 m時，最大彎曲應變?yōu)?.31×10-4；當管道外徑為0.76 m時，管道最大彎曲應變?yōu)?.902×10-4。這種現(xiàn)象在管道最大彎曲應變中更明顯，其原因與上述相同。管道最大彎曲應變增長率的絕對值降低；同時，管道最大彎曲應變出現(xiàn)在滑坡末端，因此應特別注意確保該段管道安全運行。

5 結(jié)論

(1)本文考慮了懸鏈線和大變形梁的特點，詳細分析了牽引式滑坡段管道的形態(tài)特征，保證所建立模型的準確性。該模型的邊界條件增加，連續(xù)性條件減少，方程形式簡單，需要求解的待定系數(shù)從10減少到7。結(jié)合本文提出的迭代算法，降低了整個解析求解難度，提高了求解效率。

(2)探討了滑坡寬度，滑坡力和管道厚度等參數(shù)對管道力學響應影響。分析表明，滑坡寬度的增大加劇了管道變形，使管道更容易破壞。因此，在鋪設管道時，應采取措施防止滑坡段擴張，并對滑床與非滑移段接觸部位的管道采取保護措施，同時通過增加管道壁厚也可降低其失效概率。