符煒東 田宇航 譚亮

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

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

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

1 計算模型

1.1 控制方程

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

質(zhì)量守恒方程

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

動量守恒方程

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

能量守恒方程

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

固體運動控制方程

式中：ρs為固體密度?為固體域當(dāng)?shù)丶铀俣仁噶浚?σs為柯西應(yīng)力張量，fs為體積力矢量。

流固耦合系統(tǒng)

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

1.2 基本原理

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

1.3 三維模型的建立

1.3.1 工程概況

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

1.3.2 流體三維模型的建立

本文采用某實際穿河工程管道所在河道橫斷面數(shù)據(jù)建模，管道管徑1m，壁厚14mm，距進口斷面50 米，距出口斷面150米，管道中心點高程177.48m，埋深1m，采用穿越的方式跨過河流，與河流成90°角。假設(shè)若干年后，河床整體由于水流沖刷下切2.5m變?yōu)閼铱展艿馈?/p>

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

表1.1 斷面沖刷前后數(shù)據(jù)表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 管道三維模型建立

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

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

2 網(wǎng)格劃分與邊界條件

2.1 網(wǎng)格劃分

2.1.1 流體網(wǎng)格劃分

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

圖2.1 流體網(wǎng)格圖Fig.2.1 fluid mesh

圖2.2 流體網(wǎng)格局部放大圖Fig.2.2 Fluid partially enlarged view of the grid

2.1.2 固體網(wǎng)格劃分

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

圖2.3 管道網(wǎng)格圖Fig.2.3Pipelinemesh

圖2.4 管道網(wǎng)格局部放大圖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 固體邊界條件

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

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

3 結(jié)果分析

3.1 流場分析

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

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

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

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

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

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

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

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

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

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

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

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

3.2 管道應(yīng)力分析

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

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

圖3.6 管道有效應(yīng)力分布圖Fig.3.6 Distribution of effective stress in pipelines

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

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

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

不同流速下懸空管道的應(yīng)力和形變?nèi)缦卤?.8

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

從表3.8 數(shù)據(jù)可的管道最大有效應(yīng)力與管道壁厚的關(guān)系（如圖3.9所示），管道最大形變與管道壁厚的關(guān)系（如圖3.10所示）水流速度與管道最大有效應(yīng)力的關(guān)系（如圖3.11），水流速度與管道最大形變的關(guān)系（如圖3.12所示）

圖3.9壁厚與最大有效應(yīng)力關(guān)系圖Fig.3.9 Wall thickness and maximum effective stress diagram

圖3.10壁厚與最大形變關(guān)系圖Fig.3.10 Wall thickness and maximum deformation diagram

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

圖3.11流速與最大有效應(yīng)力關(guān)系圖Fig.3.11 The maximum effective stress and flow diagram

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

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

4 結(jié)語

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