埋地管道關閥水錘數(shù)值模擬及結構響應流固耦合分析

虞學軍 鄭舟斌 熊偉東 金建波 張文奇

（1.舟山市特種設備檢測研究院；2.杭州市特種設備檢測研究院；3.浙江大學能源工程學院）

管道廣泛應用于城市供水/暖、化工及核電等行業(yè)中，對于長距離管道,為了布設需求，常常以埋地的方式進行布管。 但有時管道部件運行狀態(tài)的突然改變（如下游閥門的突然關閉）會引起管道內流體流動狀態(tài)的瞬間變化， 產(chǎn)生流體沖擊力，并在介質彈性作用下持續(xù)較長時間。 這種管道流體動態(tài)過渡的過程， 會導致系統(tǒng)部件振動，嚴重時會引起管道的破壞和斷裂，這種現(xiàn)象稱為水錘［1～4］。 水錘引起的埋地管道損壞，使埋地管道的檢測頻率和成本大幅提高。 因此，研究水錘造成的沖擊對管道的可靠性分析和事故預測至關重要。

水錘導致埋地管道振動甚至損壞，是一種典型的流體流動現(xiàn)象引起的結構振動。 需要在研究水錘的水動力學特性的基礎上，再進行管道結構力學的研究。

對于水錘的水動力學特性的試驗研究，Ruus E和Karney B B在考慮管道常數(shù)和摩擦參數(shù)的基礎上， 探討止回閥關閉后引起的水錘波特性，確定了水錘影響下泵不同位置壓頭的升高和降低的規(guī)律［5］。在理論研究的方法方面，早期主要基于特征線法進行。 如Tian W X等通過特征線法研究并聯(lián)泵交替啟動過程中的閥門水錘現(xiàn)象，并提出采用阻尼矩陣降低止回閥的關閉速度的方法，以減輕閥門水錘的潛在危害［6］。 由于特征線法主要針對管道系統(tǒng)研究，結果只可得到關鍵節(jié)點的動態(tài)數(shù)據(jù)，無法從內流場角度進行分析，也就無法進一步將內流場的壓力等數(shù)據(jù)作為載荷進一步研究結構力學的流固耦合。Al-Khomairi A M通過實驗證明基于特征線法計算瞬態(tài)管道流動存在誤差［7］。

基于管道水錘的水動力學特性的研究，Keramat A等率先進行管道水錘的流固耦合特性研究，但采用的方法是特征線法-有限元法，只能計算管道的一維變形［8］。 隨著計算流體動力學（CFD）的發(fā)展和計算機技術的進步。 CFD方法相比傳統(tǒng)的一維方法，其結果具有更高的準確性且更加詳細。 Wu D Z等采用Fluent軟件對管道系統(tǒng)的閥門快速開啟過程進行模擬， 模擬結果與PIV實驗結果有良好的吻合度，得到過渡過程中內流暢的瞬態(tài)演化過程［9］。 將CFD應用在水錘的研究中，尤其是將CFD得到的壓力等流場數(shù)據(jù)作為載荷，并通過計算結構動力學（CSD）方法研究管道結構在水錘影響下的響應特性，可以得到傳統(tǒng)一維方法難以計算的管道水錘流固耦合特性。

綜上所述， 探索一種將CFD和CSD耦合的計算方法，以準確預測埋地管道的關閥水錘水動力學特性，并得到管道結構的響應特性，在得到管道水錘的流固耦合機理的同時，還能對水錘引起管道破壞的位置進行預測，減少埋地管道檢測的成本。 因此，以Bergant經(jīng)典管道關閥水錘試驗［2，3］為基礎，通過多種CFD方法探究最優(yōu)數(shù)值解法，再采用流固耦合的方法， 借由CSD軟件分析埋地管道在水錘壓力波載荷下的響應特性，揭示管道結構的振動機理。

1 數(shù)值計算模型

通過流固耦合方法分析管道關閥水錘的水動力學特性和結構響應特性， 需要分別建立CFD計算模型和CSD計算模型。

1.1 CFD計算模型

澳大利亞研究委員會（ARC）曾經(jīng)針對管道的水錘進行了試驗研究。 該水錘試驗主要由Bergant A完成。 試驗臺由一根長37.2 m，內徑22 mm，壁厚1.6 mm的管道作為試驗段（圖1），兩端連接有球閥和水箱，球閥的下游設有段管道和水箱，兩水箱間連有壓力調節(jié)器，并由計算機控制水箱的壓力。 此外，試驗臺還包括安裝在管道試驗段的壓力傳感器和數(shù)據(jù)采集儀。 Bergant A完成了多組試驗，通過關閥使管道試驗段中產(chǎn)生水錘波，并測量水錘波的特性。 其中，管道入口壓力為0.32 MPa， 流場初始速度為0.1 m/s被武漢理工大學的楊成選取，并通過ANSYS CFX計算軟件進行數(shù)值模擬，并分析了多種初始流速和關閥速度對水錘特性的影響，但其計算值與試驗值吻合性較差［10］，而且沒有得到壓力波的衰減規(guī)律。 為此，筆者采用相同的工況，在上述研究的基礎上繼續(xù)進行數(shù)值方法的研究，以探索水錘數(shù)值模擬的最優(yōu)方法。

圖1 水錘試驗臺示意圖

1.1.1 網(wǎng)格的收斂性分析

為了實現(xiàn)對發(fā)射系統(tǒng)的CFD模擬， 需要對管道試驗段的求解區(qū)域進行網(wǎng)格剖分。 相比于非結構化網(wǎng)格， 結構化網(wǎng)格劃分方法可在減少網(wǎng)格數(shù)的同時，大幅提高網(wǎng)格的質量，尤其可以有針對性地控制特定區(qū)域的節(jié)點分布， 提高計算的收斂性和光滑性。此外，為了避免網(wǎng)格劃分對計算結果產(chǎn)生影響， 采用Roache提出的網(wǎng)格收斂系數(shù)法（GCI方法） 對網(wǎng)格收斂性進行檢驗。 GCI方法的指標——GCI值是模擬值偏離漸進數(shù)值的度量。選擇兩種網(wǎng)格劃分方案，即粗糙網(wǎng)格（網(wǎng)格數(shù)為96 065）和加密網(wǎng)格（網(wǎng)格數(shù)為704 267），如圖2所示，加密網(wǎng)格方案中網(wǎng)格的3個方向的尺寸，均為粗糙網(wǎng)格劃分方案的一半。 經(jīng)計算，加密網(wǎng)格方案的GCI值小于1%，滿足網(wǎng)格收斂性要求。 因此，采用加密網(wǎng)格方案進行數(shù)值模擬。

圖2 網(wǎng)格劃分方案示意圖

1.1.2 求解方法和設置

為了模擬關閥引起的水錘波，需要選擇合適的模擬閥門關閉過程的方法——控制邊界法和模擬關閥法。

控制邊界法，即設定連接閥門的管道出口為速度出口邊界條件，并且在一定時間內衰減為0 m/s，這也是文獻［10］選擇的方法。 在選擇的工況中，速度出口的邊界條件初始為0.1 m/s，并且在0.009 s內線性衰減為0 m/s；管道入口處的邊界條件為壓力入口，恒壓0.32 MPa。

模擬關閥法， 即額外建立球閥的流體域，通過旋轉球閥的流體域實現(xiàn)關閥過程的模擬，該方法的流體域設置如圖3所示， 球閥的兩個交界面分別與管道的上下游貼合，此時為球閥的關閉狀態(tài)，管道為不流通狀態(tài)。 在選擇的工況中，球閥初始為全開，并且在0.009 s內勻速關閉；管道入口處的邊界條件同前一種。

圖3 模擬關閥數(shù)值方法的流體域

由于文獻［10］ 基于CFX求解器進行數(shù)值模擬， 為了探尋模擬水錘事件的最優(yōu)方法， 采用Fluent求解器進行上述兩種方法的求解。 求解的過程中，湍流模型選用對強逆壓梯度的流動具有很好的計算精度的SST k-ω；離散方程的求解采用壓力速度耦合的Coupled算法。閥門網(wǎng)格的轉動采用動網(wǎng)格方法實現(xiàn)， 通過通過用戶自定義函數(shù)（UDF）控制閥門勻速關閉，實現(xiàn)閥門由全開度至閉合狀態(tài)。 在計算中，考慮流體的可壓縮性，通過UDF控制密度隨壓力變化來反映水的可壓縮性變化規(guī)律，其表達式為：

式中 K——水的體積彈性模量，2.18GPa；

p——水的實時壓力，Pa；

pop——標準狀況下大氣壓力，101 325 Pa；

ρ——水的密度，kg/m3；

ρref——標準狀態(tài)下水的密度，998.2 kg/m3。

單個時間步長為0.000 5 s， 單個時間步迭代20次；總的計算時長為0.5 s，共1 000個時間步。

1.2 CSD計算模型

對在水錘波作用下的管道結構響應特性的分析，采用與試驗相同的管道，但是為了貼合實際工況，材料選用不銹鋼以替代銅。 同樣，對管道進行結構化網(wǎng)格劃分，網(wǎng)格如圖4所示。

圖4 CSD計算模型網(wǎng)格

計算模型的載荷包括重力載荷和隨時間變化的內部流體壓力載荷。 為了加載時刻變化著的壓力載荷，先選用ANSYS的Transient Structural模塊進行結構動力學計算，該模塊可以充分考慮載荷變化。再通過ANSYS的ACT擴充功能，提取CFD計算得到的管道壁面的實時載荷，同時提取載荷的大小、時間和網(wǎng)格位置信息，并作為載荷施加給管道結構。

采用彈性約束的方式模擬土壤對管道的作用，彈性約束以管道的外壁面為邊界，基礎的剛度為使管道結構產(chǎn)生單位法向偏移所需要的壓力。 不同地基土的抗壓剛度系數(shù)各不相同，相差數(shù)倍或十幾倍，因此以典型的碎石土（硬質土壤）為管道埋地環(huán)境，其抗壓剛度為70 kN/m3，并作為約束。

2 數(shù)值模擬與試驗結果對比

水錘波引起的壓力變化是水錘最重要的特性，因此選用水錘事件中的壓力變化，作為衡量數(shù)值模擬與真實情況貼合程度的指標。 分別總結Bergant A試驗、文獻［10］和上述基于Fluent求解結果中的閥門上游點和管道軸向中點的壓力隨時間變化的規(guī)律，其中，閥門點處壓力變化規(guī)律如圖5所示。

圖5 閥門點處壓力變化規(guī)律

從水錘波的周期入手，試驗測得的水錘波周期約為0.135 s，CFX-控制邊界法為0.095 s，F(xiàn)luent的兩種解法同為0.105 s。 按照水錘波理論，其壓力波周期T和波速a的關系式為：

式中 L——管長，m。

由式（2）計算得到的水錘波速分別為：試驗值1 102 m/s；CFX-控制邊界法1 566 m/s；Fluent控制邊界和模擬關閥法1 417 m/s。 可見，F(xiàn)luent模擬水錘波的波速和周期更接近試驗值，模擬值和試驗值的偏差緣由為：水錘波會在試驗中管道上游的水箱中繼續(xù)傳遞， 而模擬中忽略了水箱，受水箱影響的試驗數(shù)據(jù)會使計算得到的波速略低；水錘發(fā)生時會伴隨著一定的空化現(xiàn)象，水錘波在水蒸氣中的傳播速度遠低于液態(tài)水，而模擬忽略了空化現(xiàn)象，算得的波速會偏高。

在水錘波壓力幅值方面， 基于Fluent的兩種方法計算結果幾乎相同， 且在第1周期中呈短時間內在最高值維持不變的規(guī)律，這與試驗結果相同，而且是CFX求解器沒能成功仿真的現(xiàn)象。對于壓力波波峰，試驗值為0.45 MPa（水頭46.01 m）；CFX求解的結果為0.49 MPa（水頭50.10 m），與試驗值偏差為8.89%， 且出現(xiàn)在計算的最后一個周期，在0.9 s時刻［10］；Fluent求解的結果為0.48 MPa（水頭49.07 m），出現(xiàn)在水錘波的第1個周期，與試驗值的偏差為6.67%。 可見，F(xiàn)luent的求解結果更精確。此外，CFX求解器沒能得到水錘波壓力波峰隨周期遞減的規(guī)律， 而Fluent求解結果中水錘波波峰遞減規(guī)律與試驗結果相同。

管道中點壓力的變化規(guī)律如圖6所示。 同樣的， 基于Fluent的模擬關閥法和控制邊界法得到的結果相同， 且在幅值和周期上相較CFX求解器更為準確。 此外，基于Fluent求解水錘的優(yōu)勢在于可以很好地模擬水錘波壓力衰減規(guī)律。

圖6 管道中點壓力變化規(guī)律

3 管道關閥水錘過程中壓力、 流速變化和內流場分析

由上述分析可知， 通過Fluent求解器得到的水錘波壓力變化規(guī)律更趨近于試驗值， 因此，采用Fluent模擬關閥法的求解結果分析管道關閥水錘過程中壓力、流速的變化和內流場變化。 除閥門處和管道中點外，同時監(jiān)測了管道入口和距管道入口5 m處的兩個測點。 各測點壓力變化規(guī)律如圖7所示。 由圖7可以看出，水錘波的壓力峰值隨距閥門的距離逐漸降低。 入口處壓力由邊界條件確定，始終保持在恒定初壓，因此與試驗的真實情況存在偏差，這也是引起數(shù)值方法產(chǎn)生誤差的重要因素。

圖7 管道各測點壓力變化規(guī)律

各測點速度變化規(guī)律如圖8所示。 對于閥門處，由于閥門關閉，閥門處流體流速迅速降低至0 m/s，期間有輕微波動。在管道的其他位置，流體的流速隨壓力波在±0.1 m/s的幅值范圍內波動。但在第1個周期內， 各測點的速度變化規(guī)律并不同步，而是存在一定的相位差，這意味著在管道的特定位置，流體流速梯度較大，流體相互擠壓，甚至產(chǎn)生相向流動。

圖8 管道各測點速度變化規(guī)律

圖9為0.127 s時閥門上游區(qū)域流場的速度云圖和矢量圖。 由圖9可看出，因閥門關閉，上游區(qū)域的流速降低至0 m/s，且流體隨水錘壓力波反向流動，但更上游的流體由于慣性，保持原先的正向流動，故在該區(qū)域產(chǎn)生了相向流動的現(xiàn)象。

圖9 管道內流體相向流動的區(qū)域

4 水錘壓力載荷下管道結構力學特性

圖10是管道結構在水錘載荷作用下，振動加速度和應力隨時間的變化規(guī)律。 對于管壁的變形加速度，壓力波的波峰和波谷均會引起管壁的振動，且幅值達到了一個重力加速度。 在同一個水錘波周期內， 變形加速度會出現(xiàn)兩次波峰和波谷。 所以，在水錘的影響下，管道會劇烈振動，與土壤的沖擊也會使得結構損壞。 對于管道結構的應力，可以發(fā)現(xiàn)管道內的初壓引起了管壁一定的壓力，當水錘波壓力波峰出現(xiàn)時由于作用在壁面的壓力陡然增大， 應力的最大值也隨即增大，當壓力值達到最大時，應力值也達到最大，升高幅度達到了46.4%，這種應力的急劇增大，極易引起結構的破壞。 而水錘壓力波波谷到達時，由于管道內部分位置的壓力降低，使得應力的最小值降低。 但由于壓力波在管道內時刻傳遞，這使得不同位置的應力不同， 既存在最小值也存在最大值。

管壁應力以上述規(guī)律呈現(xiàn)明顯的周期性，且隨著周期增大應力的變化減小。 圖11為4個不同時刻應力最大值和最小值出現(xiàn)的位置，4個時刻分別為圖10所示的t1～t4。 由圖11可以看出，水錘壓力波波峰時，應力最大出現(xiàn)在閥門處，且最小值在入口處管段；當壓力波波谷時，應力最大值在管道入口，而最小值在管道中段。 而且，下一個周期也呈現(xiàn)出相同的規(guī)律。

圖11 管道應力最大值和最小值位置

5 結論

5.1 相對于CFX，F(xiàn)luent求解水錘壓力波的幅值和變化規(guī)律更為準確，且控制邊界法和模擬關閥法得到的結果相同。 為了得到更為準確的周期信息，需要建立包括水箱在內的全流體域。

5.2 在水錘波的影響下，管道內流體流動方向發(fā)生改變，但不同位置的流速改變不同步，在特定位置會出現(xiàn)流體相向流動，產(chǎn)生擠壓。

5.3 水錘壓力波會引起管道劇烈振動，且一個壓力波周期內，會導致兩個振動周期。

5.4 水錘壓力波波峰會引起管道結構的應力最大值顯著增大， 增幅達46.4%， 極易引起結構破壞。