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

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

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

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

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

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

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

1 數(shù)值計算模型

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

1.1 CFD計算模型

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

圖1 水錘試驗臺示意圖

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

為了實現(xiàn)對發(fā)射系統(tǒng)的CFD模擬， 需要對管道試驗段的求解區(qū)域進行網(wǎng)格剖分。 相比于非結(jié)構(gòu)化網(wǎng)格， 結(jié)構(gòu)化網(wǎng)格劃分方法可在減少網(wǎng)格數(shù)的同時，大幅提高網(wǎng)格的質(zhì)量，尤其可以有針對性地控制特定區(qū)域的節(jié)點分布， 提高計算的收斂性和光滑性。此外，為了避免網(wǎng)格劃分對計算結(jié)果產(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 求解方法和設(shè)置

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

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

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

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

由于文獻［10］ 基于CFX求解器進行數(shù)值模擬， 為了探尋模擬水錘事件的最優(yōu)方法， 采用Fluent求解器進行上述兩種方法的求解。 求解的過程中，湍流模型選用對強逆壓梯度的流動具有很好的計算精度的SST k-ω；離散方程的求解采用壓力速度耦合的Coupled算法。閥門網(wǎng)格的轉(zhuǎn)動采用動網(wǎng)格方法實現(xiàn)， 通過通過用戶自定義函數(shù)（UDF）控制閥門勻速關(guān)閉，實現(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計算模型

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

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

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

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

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

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

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

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

式中 L——管長，m。

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

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

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

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

3 管道關(guān)閥水錘過程中壓力、 流速變化和內(nèi)流場分析

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

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

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

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

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

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

4 水錘壓力載荷下管道結(jié)構(gòu)力學特性

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

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

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

5 結(jié)論

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

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

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

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