虞學(xué)軍 鄭舟斌 熊偉東 金建波 張文奇

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

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

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

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

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

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

1 數(shù)值計(jì)算模型

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

1.1 CFD計(jì)算模型

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

圖1 水錘試驗(yàn)臺(tái)示意圖

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

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

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

1.1.2 求解方法和設(shè)置

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

控制邊界法，即設(shè)定連接閥門的管道出口為速度出口邊界條件，并且在一定時(shí)間內(nèi)衰減為0 m/s，這也是文獻(xiàn)［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所示， 球閥的兩個(gè)交界面分別與管道的上下游貼合，此時(shí)為球閥的關(guān)閉狀態(tài)，管道為不流通狀態(tài)。 在選擇的工況中，球閥初始為全開，并且在0.009 s內(nèi)勻速關(guān)閉；管道入口處的邊界條件同前一種。

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

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

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

p——水的實(shí)時(shí)壓力，Pa；

pop——標(biāo)準(zhǔn)狀況下大氣壓力，101 325 Pa；

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

ρref——標(biāo)準(zhǔn)狀態(tài)下水的密度，998.2 kg/m3。

單個(gè)時(shí)間步長為0.000 5 s， 單個(gè)時(shí)間步迭代20次；總的計(jì)算時(shí)長為0.5 s，共1 000個(gè)時(shí)間步。

1.2 CSD計(jì)算模型

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

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

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

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

2 數(shù)值模擬與試驗(yàn)結(jié)果對(duì)比

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

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

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

式中 L——管長，m。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖11 管道應(yīng)力最大值和最小值位置

5 結(jié)論

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

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

5.3 水錘壓力波會(huì)引起管道劇烈振動(dòng)，且一個(gè)壓力波周期內(nèi)，會(huì)導(dǎo)致兩個(gè)振動(dòng)周期。

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