洪水作用下穿河懸空輸氣管道的動力響應研究

符煒東 田宇航 譚亮

(1.四川大學水力學及山區河流開發保護國家重點實驗室，四川 成都 610000；2. 四川省廣元水文水資源勘測局，四川 廣元 628000）

穿河管段的安全性問題是能源領域非常典型的工程問題。針對穿河管道，人們起初關注的是管道的合理埋深，在保證管道安全的基礎上力求減小工程量，從而達到節省工程投資的目的。為此，國內外學者總結了一些在一定條件下合適的計算沖刷深度的公式：如1981年，俞樂群借用橋渡建筑物附近一般沖刷深度的公式來計算管道穿越河流處的河床沖刷深度[1]；2011年Azamathulla 將人工神經網絡方法帶入到管道沖刷計算中[2]。有學者針對一些實際工程作出安全性分析：1998 年，黃金池提出了一些穿河管段保護措施[3]；2002年，先智偉分析了水下穿越事故的主要特征，指出了穿越事故主要原因，并提出相應防護措施[4]。隨著計算機技術的發展，人們開始利用各類軟件對管道的荷載進行研究：2012 年，徐濤龍等對穿越管道臨界懸空高度進行動態數值模擬[5]；2013 年，姚安林等使用數值模擬的方法，對不同壁厚穿河管段的臨界懸空高度進行了研究[6]；

這些學者大多關注管道的埋深、管道懸空壁厚、管道的懸空長度、管道的應力特性等單一參數，而很少關注洪水與管道之間產生流固耦合致使管道發生形變以及管道的應力特征改變。對跨河管道進行數值模擬時，通常將河道簡化為矩形水槽，忽略了河床的地貌變化及河流自由水面的模擬，致使水流的流態與實際間存在差異，因而對于實際工程問題，數值模擬得到結果的準確性難以保證。為此，本文使用管道所在河道橫斷面數據建模，基于Workbench 軟件對某山區穿河管工程進行流固耦合仿真分析，得出在一定的沖刷深度下，不同流速的洪水及管道壁厚所對應的懸空管段的應力特性及管道破壞模式。

1 計算模型

1.1 控制方程

對于一般的可壓縮牛頓流體，守恒定律通過如下控制方程描述控制方程包括質量守恒方程、動量守恒方程、能量守恒方程，流體運動方程的具體表達式如下：

質量守恒方程

式中:t為時間,ρf為流體密度,υ為流體速度矢量。

動量守恒方程

式中：t為時間,ρf為流體密度,υ為流體速度矢量，ff為體積力矢量。

能量守恒方程

式中：t為時間,ρf為流體密度,υ為流體速度矢量，ff為體積力矢量，qf為單位體積熱量損失，E為單位質量內能。

固體運動控制方程

式中：ρs為固體密度?為固體域當地加速度矢量，?σs為柯西應力張量，fs為體積力矢量。

流固耦合系統

式中：τf、τs為流體、固體應力，nf、ns為流體、固體單位方向向量，df、ds為流體、固體位移。

1.2 基本原理

流固耦合問題分為有單向流固耦合和雙向流固耦合。由于懸空管道形變尺度相較于流場很小，可以忽略不計，故采用單向耦合法[8]。本文依據某實際穿河工程管道所在橫斷面數據建立流體三維模型與管道三維模型，將流體網格導入Fluent 進行不同流速下的仿真計算，提取管道所受阻力與升力，再將Fluent的計算結果作為管道荷載導入workbench里的單項流固耦合模塊(FSI)得到管道在沖擊下的形變及應力特征。

1.3 三維模型的建立

1.3.1 工程概況

某輸氣管道全長136.5公里，管徑1000毫米，管道沿線地形地貌、地質構造復雜，全線均為中低山和丘陵地貌區，山巒起伏、溝壑縱橫，管道多次穿越大小河流[7]。該管道穿越一河流時，管道順斜坡敷設，管道直徑1000mm，埋深1m，采用大開挖方式，回填土較松散，地表水土流失嚴重，在地表形成沖溝，少數沖溝已深達1-2米。

1.3.2 流體三維模型的建立

本文采用某實際穿河工程管道所在河道橫斷面數據建模，管道管徑1m，壁厚14mm，距進口斷面50 米，距出口斷面150米，管道中心點高程177.48m，埋深1m，采用穿越的方式跨過河流，與河流成90°角。假設若干年后，河床整體由于水流沖刷下切2.5m變為懸空管道。

斷面數據如表1.1，沖刷后的管道與河床示意圖如圖1.2。研究的洪水高程為178.98m，考慮自由液面波動，使用水平線181m與沖刷后的河床斷面數據建立流體三維模型，如圖1.3。

表1.1 斷面沖刷前后數據表Tab.2.1 Data sheet before and after cross section scouring

圖1.2相對位置示意圖Fig.1.2Relative position diagram

圖1.3流體三維模型圖Fig.1.3 3D model of fluid

1.3.3 管道三維模型建立

管道模型根據5 種不同的工況的壁厚分別進行建模，管道兩端截止于河床邊壁，懸空長度為40m左右。管道三維模型見圖1.4

圖1.4管道三維模型圖Fig.1.4 3D model of pipeline

2 網格劃分與邊界條件

2.1 網格劃分

2.1.1 流體網格劃分

將流體三維模型導入Icem劃分網格，對管道中間網格部分進行加密，進口段和出口段采用結構化網格劃分，管道中間段由于管道距離河床底太近，且河床底部斷面形狀不規則，故采用非結構化網格劃分。流體網格見圖2.1、圖2.2。

圖2.1 流體網格圖Fig.2.1 fluid mesh

圖2.2 流體網格局部放大圖Fig.2.2 Fluid partially enlarged view of the grid

2.1.2 固體網格劃分

固體部分根據5種不同工況的壁厚的管道進行建模并導入Workbench中進行網格劃分，管道網格劃分見圖2.3、圖2.4。

圖2.3 管道網格圖Fig.2.3Pipelinemesh

圖2.4 管道網格局部放大圖Fig.2.4 Pipeline grid enlarged view

2.2 邊界條件

2.2.1 流體邊界條件

流體求解器類型選擇pressure-based，速度選擇absolute，時間選擇transient，考慮重力，由于建立的模型中流體包括水和空氣，計算模型選擇多相流（Volumeof Fluid），k-epsilon（2eqn），迭代中采用SIMPLE 算法。流體模型入口邊界條件為速度入口，上邊界為空氣出口，下邊界為wall，出口邊界為壓力出口，耦合面為管道外壁面（wall）。時間步長采用0.01s，分別模擬流速4m/s、5m/s、6m/s、7m/s、8m/s的流場。

2.2.2 固體邊界條件

大多數情況下由于土壤的剛度等因素影響，穿河管道兩端既不是兩端固定也不是純鉸支的情況，一般鑒于兩者之間[9]，本例采用一端固定另一端簡支的約束條件使結果更符合實際情況。

穿河管道采用西氣東輸工程常用的L485 管線鋼，密度7900kg/m3，彈性模量210Gpa，切線模量13.5Gpa，泊松比0.3，屈服應力485Mpa[10]。外徑1m，分別模擬11mm、14mm、17mm、20mm、23mm壁厚下管道的應力特征及形變。

3 結果分析

3.1 流場分析

從下圖3.1 組分分布圖可以看出當流速為4m/s 時，由于管道阻隔，洪水自入口斷面到管道附近水位一直在壅高，到管道處洪水水位高于入口斷面設置水位，在管道后發生明顯的跌水現象，水位迅速降低，在管道下游，水位略微抬高后又小幅下降，最終流態穩定。

圖3.1 組分分布圖Fig.3.1 Component distribution diagram

從下圖3.2流速矢量圖可以看出，在洪水行進過程中，洪水到達管道附近時，由于管道阻隔，水流流態比較紊亂。從進口到接近管道過程中，水流速度從4m/s 在逐漸減小，洪水越過管道后流速增加，管道上下兩端附近流速較高，最大為7.97m/s，而后流速逐漸減小后趨于穩定，比較圖3.1 與3.2 可以發現，水深變化規律與流速變化規律是協調統一的，即流速增加，水深減小，流速降低，水深增加。

圖3.2 流速矢量圖Fig.3.2 Velocity vector

從下圖3.3壓力分布圖可以看出管道迎水面中點受到的壓強最大為31KPa，管道上端壓力最小，產生負壓，管道底部壓力較上端大。迎水面中點受到的水流沖擊最大，故所受壓強最大。頂部水深較小而流速較大，根據伯努利能量方程可知，管道頂部出現負壓是合理的[11]。

圖3.3 壓力分布圖Fig.3.3 Pressure profile

流體沖擊管道時，會與管道產生相互作用，作用力在來流方向的分力為阻力（FD），與來流方向垂直方向的分力為升力（FL）。阻力可分為摩擦阻力和壓差阻力，其中摩擦阻力是黏性力作用的直接結果，壓差阻力是邊界層分離的結果，管道屬于鈍體，壓差阻力占阻力的絕大部分[12]。升力是與來流方向垂直方向上產生的壓力差，本例中由于河床與管壁距離太近，全斷面平均間隙比小于0.5（管道底部與河床表面的間隙，該間隙與管道直徑的比值定義為間隙比），屬于高雷諾數下的近壁圓柱繞流，根據相關文獻，壁面阻礙了下部剪切層與外側流動間的相互作用，使得下側旋渦得不到充分發展，并且間隙比越小，與壁面接觸的漩渦越靠前，壁面對漩渦脫落的抑制作用越大，致使升力的時均值為正數[13-14]。流場穩定后，升力大于阻力保持在257879N 左右，與壓強分布圖及文獻描述相同，驗證了數值模擬的可靠性。

以上為來流流速為4m/s 時的流場分析，通過改變流速，得到不同流速下的管道所受升力及阻力如下表3.4。從表中可以得知，管道所受的升力FL大于所受到的阻力FD，兩者均為威脅管道安全的主要因素。

表3.4 速度、升力及阻力關系表Tab.3.4 Speed,lift and resistance relationship table

從表3.4可以得出水流流速與升力、阻力關系圖（如圖4.5）

圖3.5 流速與升力、阻力關系圖Fig.3.5 Relationship between flow velocity,lift and resistance

從圖4.5可以得知，水流速度增大，管道所受的升力與阻力也增大，且隨著水流速度的增加，管道所受的升力與阻力增大的越來越快，管道失穩的風險也越來越大。

3.2 管道應力分析

分別將4m/s、5m/s、6m/s、7m/s、8m/s入口速度生成的荷載結果導入Workbench 中，再分別模擬11mm、14mm、17mm、20mm、23mm 壁厚懸空管道的應力及形變。以4m/s 時，工程中所使用的14mm壁厚L485管道為例。

如圖3.6 所示，管道的最大有效應力產生在管道左端底部為487.88MPa，管道使用L485管線鋼，管線鋼為塑性材料，適用第三強度理論，即不論材料處于何種應力狀態，只要最大有效應力超過材料屈服應力材料就屈服[15]。L485 管線鋼屈服應力為485Mpa，可知該實際穿河管道工程采用的14mm壁厚的L485管道在河床沖刷下切2.5m后，流速4m/s時難以抵擋洪水的沖擊而失穩。比較各工況下的管道最大有效應力產生處均為管道左端底部，說明失穩破壞最先發生在管道左端底部。

圖3.6 管道有效應力分布圖Fig.3.6 Distribution of effective stress in pipelines

如下圖3.7 所示，管道的最大形變發生在距左側管道25m處。

圖3.7 管道形變分布圖Fig.3.7Pipeline deformation distribution

以上為流速4m/s、管壁14mm時的工況，通過改變流速大小及管壁壁厚可以得出，該工程懸空管道的有效最大應力總是發生在管道底部左端與土壤相交處附近，這與河床斷面形狀以及流體形態有關。

不同流速下懸空管道的應力和形變如下表3.8

表3.8 不同水流速度及壁厚下管道最大有效應力及最大形變表Tab.3.8 maximum effective stress and maximum displacement under different water velocity and wall thickness

從表3.8 數據可的管道最大有效應力與管道壁厚的關系（如圖3.9所示），管道最大形變與管道壁厚的關系（如圖3.10所示）水流速度與管道最大有效應力的關系（如圖3.11），水流速度與管道最大形變的關系（如圖3.12所示）

圖3.9壁厚與最大有效應力關系圖Fig.3.9 Wall thickness and maximum effective stress diagram

圖3.10壁厚與最大形變關系圖Fig.3.10 Wall thickness and maximum deformation diagram

由4.9 及4.10 可以看出，管道的最大有效應力及最大形變隨著管道壁厚的減小而增大，且隨著管道壁厚減小，最大有效應力增大的越來越快；

圖3.11流速與最大有效應力關系圖Fig.3.11 The maximum effective stress and flow diagram

圖3.12流速與最大形變關系圖Fig.3.12 Flow rate and the maximum deformation diagram

由3.11及3.12可以看出，管道的最大有效應力及最大形變隨著水流速度的增大而增大，且隨著管道流速增大，最大有效應力增大的越來越快。

4 結語

本文依托某實際穿河管道工程，基于workbench平臺，采用流固耦合的數值計算方法，分析了不同流速下管道周圍的流場特性及管道本身的動力響應特征，研究結果表明：管道懸空后，上游段出現壅水現象，水位升高，管道處形成跌水，水位迅速下降；管道迎水面中線所受壓力最大，背水面次之；本文模擬工況下的升力FL均大于阻力FD。不同工況下的最大有效應力位置均為管道左端底部與岸邊交接處，根據不同工況下，壁厚、最大有效應力、最大形變及流速的關系可以得出，流速越大，壁厚越小，最大有效應力及形變越大且增長越快，管道越容易屈服破壞。