劉家峽洮河口排沙洞水下巖塞爆通后多相流運(yùn)動(dòng)數(shù)值模擬

徐小武，薛立梅，李俊富，吳希華

（1.中水東北勘測(cè)設(shè)計(jì)研究有限責(zé)任公司，吉林 長(zhǎng)春 130061；2.水利部松遼水利委員會(huì)，吉林 長(zhǎng)春 130021）

1 概 述

劉家峽水電站位于甘肅省永靖縣的黃河干流上，主壩為混凝土重力壩，最大壩高147 m，總庫(kù)容57.01億m3，總裝機(jī)容量1390 MW，為一等大（1）型工程。

洮河水少沙多，在距劉家峽大壩上游1.5 km的右岸匯入黃河干流。洮河口段死庫(kù)容于1987年淤滿，在河口附近黃河干流形成沙坎，且淤積面逐年抬高，現(xiàn)有排沙設(shè)施已不能解決洮河泥沙淤積并向壩前推移的問題，給電站的安全運(yùn)行和渡汛造成了嚴(yán)重危害。

為解決泥沙淤積對(duì)大壩安全運(yùn)行帶來的嚴(yán)重危害和對(duì)大壩安全渡汛構(gòu)成的威脅及電站機(jī)組磨損等問題，在劉家峽水電站左岸增建排沙洞工程，排沙洞的進(jìn)水口采用水下巖塞爆破方案。

巖塞爆破口位于洮河出口，黃河左岸；排沙洞全長(zhǎng)1486 m，巖塞進(jìn)口底板高程為1664.53 m；巖塞內(nèi)口為圓形，內(nèi)徑10 m，外口近似橢圓，外口尺寸約20.3 m×27.84 m，巖塞最小厚度12.3 m，巖塞進(jìn)口軸線與水平面夾角45°，塞體方量約2474 m3。

2 計(jì)算目的和內(nèi)容

排沙洞進(jìn)口段采用水下巖塞爆破施工。在排沙洞巖塞口爆破區(qū)范圍內(nèi)，淤積層厚度變化較大，11～58 m；另外，淤積層內(nèi)有522 m3的水下人工堆渣，渣塊塊徑一般5～10 cm，最大塊徑100 cm。厚淤泥層的存在對(duì)巖塞口爆破會(huì)帶來不利影響；另外，在排沙洞運(yùn)行過程中，若渣塊進(jìn)入洞內(nèi)，也會(huì)影響排沙洞的正常運(yùn)行，成為工程的隱患。為了進(jìn)一步掌握排沙洞水下巖塞爆破后的水沙運(yùn)動(dòng)規(guī)律，有必要開展排沙洞水下巖塞爆通后多相流運(yùn)動(dòng)數(shù)值模擬計(jì)算，為排沙洞水下巖塞爆破方案設(shè)計(jì)提供參考。

為了解水下巖塞爆破后厚覆蓋層對(duì)排沙洞內(nèi)水沙運(yùn)動(dòng)的影響程度，本文基于當(dāng)前多相流計(jì)算技術(shù)的發(fā)展水平，假定在排沙洞水下巖塞爆破后，爆破巖塞所產(chǎn)生的巖渣瞬間全部堆滿預(yù)挖的集渣坑，忽略巖塞巖渣和集渣坑的影響。采用Euler多相流模型，以排沙洞巖塞口區(qū)及排沙洞為研究對(duì)象，建立三維排沙洞多相流非恒定多相流計(jì)算模型；計(jì)算正常蓄水位、淤泥層不同顆粒組成等條件下，排沙洞水下巖塞爆破后排沙洞中多相流運(yùn)動(dòng)過程；對(duì)比分析不同時(shí)間過程排沙洞內(nèi)各相體積率及其流動(dòng)速度、壓力分布等。

3 計(jì)算模型原理和方法

3.1 多相流模型

排沙洞巖塞口區(qū)覆蓋層較深且夾雜人工堆渣，在排沙洞巖塞爆通后,巖塞口區(qū)的淤泥將夾帶人工堆渣流入排沙洞。夾帶堆渣的淤泥在排沙洞內(nèi)的流動(dòng)是一個(gè)多相流問題，需要建立排沙洞多相流數(shù)學(xué)模型來研究分析排沙洞內(nèi)淤泥的流動(dòng)規(guī)律。目前，處理多相流多采用歐拉模型的計(jì)算方法。

歐拉模型是一種非常復(fù)雜的多相流模型，它包含有n個(gè)動(dòng)量方程和連續(xù)方程，來求解每一相；壓力項(xiàng)和各界面交換系數(shù)是耦合在一起的。耦合的方式則依賴于所含相的情況，顆粒流(流—固)的處理與非顆粒流(流—流)是不同的。對(duì)于顆粒流，可應(yīng)用分子運(yùn)動(dòng)理論來求得流動(dòng)特性。不同相之間的動(dòng)量交換也依賴于混合物的類別。

3.2 流體流動(dòng)的基本控制方程

流體流動(dòng)受物理守恒定律的支配，要滿足質(zhì)量守恒方程、動(dòng)量守恒方程。對(duì)于定常流動(dòng)，密度為常數(shù)的情況下，質(zhì)量守恒方程可表述為：

對(duì)于不可壓縮常黏度的流體，在不考慮流體黏性的情況下，歐拉方程可表述為：

3.3 湍流模型

湍流是一種高度復(fù)雜的三維非穩(wěn)態(tài)、帶旋轉(zhuǎn)的不規(guī)則流動(dòng)。湍流中流動(dòng)的各個(gè)物理量參數(shù)，如流速、壓力、溫度等都隨時(shí)間和空間發(fā)生隨機(jī)變化，特點(diǎn)是物理量的脈動(dòng)，非穩(wěn)態(tài)的N-S方程對(duì)湍流運(yùn)動(dòng)仍是適用的。湍流模型[1]分為標(biāo)準(zhǔn)k-ε湍流模型、RNGk-ε模型、Realizablek-ε模型3種，這些模型均是針對(duì)湍流發(fā)展非常充分的湍流流動(dòng)來創(chuàng)建的。

3.4 計(jì)算模型的求解參數(shù)

計(jì)算模型的參數(shù)選取對(duì)于數(shù)學(xué)模型的計(jì)算結(jié)果至關(guān)重要。正確設(shè)置求解的控制方程、壓力速度耦合方法、松弛因子、離散格式等參數(shù)，有利于提高模型計(jì)算的仿真精度。

1）求解的控制方程

在求解參數(shù)設(shè)置中，可選擇的方程包括流動(dòng)方程、湍流方程、能量方程、體積分?jǐn)?shù)方程等。在求解過程中，有時(shí)為了得到收斂的解，等一些簡(jiǎn)單的方程收斂后，再開啟復(fù)雜的方程一起計(jì)算。

2）壓力速度耦合方法選擇

在基于壓力求解器中，一般有SIMPLE，SIMPLEC，PISO以及Coupled等4種壓力速度耦合的方法。定常狀態(tài)計(jì)算一般使用SIMPLE或SIMPLEC方法，但對(duì)于許多問題如果使用SIMPLEC可能會(huì)得到更好的結(jié)果。

3）亞松弛因子

亞松弛因子是基于壓力求解器所使用的加速收斂參數(shù)，用于控制每個(gè)迭代步內(nèi)所計(jì)算的場(chǎng)變量的更新。除耦合方程之外的所有方程，包括耦合隱式求解器中的非耦合方程（如湍流方程），均有與之相關(guān)的亞松弛因子。一般，壓力、動(dòng)量、k和ε的亞松弛因子分別為0.2，0.5，0.5和0.5。

4 計(jì)算模型建立

4.1 計(jì)算范圍及計(jì)算網(wǎng)格

按照設(shè)計(jì)要求，模型計(jì)算水位選用正常蓄水位1735 m；庫(kù)區(qū)水平范圍以巖塞進(jìn)口中心點(diǎn)為中心、按75 m爆破影響半徑范圍選?。挥?jì)算假定在排沙洞水下巖塞爆破后，爆破巖塞所產(chǎn)生的巖渣瞬間全部堆滿預(yù)挖的集渣坑，忽略巖塞巖渣和集渣坑的影響。

選取的計(jì)算庫(kù)區(qū)尺寸為99 m×150 m×78.82 m(長(zhǎng)×寬×高)；巖塞為圓臺(tái)型，上口直徑為22.3 m、下口直徑為10 m、高為12.3 m，圓臺(tái)軸線傾角為45°，巖塞進(jìn)口中心點(diǎn)高程為1673.18 m；巖塞的下面接直徑為10 m、高為3 m的圓柱；圓柱后面接一直徑為10 m、中心角為51°的圓環(huán)；圓環(huán)下游接一坡度為6°的逆坡段圓柱，圓柱直徑10 m，長(zhǎng)度116.86 m；最后與逆坡段圓柱連接的是一段長(zhǎng)度為1286.417 m、直徑為10 m的順坡圓柱，坡角為1.488°，圓柱出口高程為1631.6 m。參照黃河劉家峽洮河口排沙洞工程巖塞地質(zhì)勘測(cè)工程地質(zhì)勘查報(bào)告及剖面圖（H-H，I-I），模型水庫(kù)計(jì)算域分為上、中、下三層：上層為空氣層，高度為5 m（高程從1735.0～1740.0 m）；中層為庫(kù)水層，高度為41.5 m（高程從1693.5～1735.0 m）；下層為淤泥層，淤泥層高度為 32.32 m（高程 1661.18～1693.5 m）、含水量為 33.9%、孔隙比為1.114。

4.2 建立求解模型

1）模型求解器設(shè)置：①開啟標(biāo)準(zhǔn)-湍流模型；②選擇k-epsilon作為 Model；③在 k-epsilon Model下，保留其默認(rèn)選擇Standard不變；④保持Near-Wall Treatment下的默認(rèn)設(shè)置Standard Wall Functions不變。

2）設(shè)置多相流模型，開啟歐拉模型選項(xiàng)。

3）設(shè)置用于計(jì)算相間動(dòng)量傳遞的拖曳規(guī)律，即將各相間相互影響的拖拽作用選擇為schiller-naumann。

4）設(shè)置重力加速度。

4.3 離散格式

對(duì)流項(xiàng)可選擇不同的離散格式：

1）一般，當(dāng)使用基于壓力求解器時(shí)，所有方程中的對(duì)流項(xiàng)均用一階精度格式進(jìn)行離散。

2）當(dāng)使用基于密度求解器時(shí)，流動(dòng)方程使用二階精度格式，其他方程使用一階精度格式進(jìn)行離散。

3）此外，當(dāng)使用分離式求解器時(shí)，可為壓力選擇插值方式。

4.4 流體材料

模型計(jì)算材料分為空氣、水和沙三類，相應(yīng)的計(jì)算參數(shù)見表1。

4.5 計(jì)算相

為了研究水下巖塞爆破后淤泥顆粒組成及堆渣含量對(duì)排沙洞中水沙運(yùn)動(dòng)規(guī)律的影響，該多相流模型的計(jì)算相大體分為空氣、水和沙粒三類。針對(duì)工程特點(diǎn)并結(jié)合多相流模型的要求，模型中將空氣作為主相，其他類均設(shè)為副相。由于淤泥層孔隙比為1.114，所以淤泥層中水的體積率為0.527。

表1 計(jì)算參數(shù)表

4.6 計(jì)算方案

按淤泥層中不同沙粒的組成，模型計(jì)算分4個(gè)方案，見表2。

表2 計(jì)算方案表

4.7 計(jì)算邊界

邊界條件的正確設(shè)置是多相流計(jì)算的關(guān)鍵。巖塞爆通后，巖塞口區(qū)淤泥將流向排沙洞，同時(shí)巖塞口區(qū)外的淤泥和水又會(huì)流向巖塞口區(qū)，使得巖塞口區(qū)的淤泥和水得到補(bǔ)充。

由于庫(kù)區(qū)邊界處各相的流速或流量均未知，為了便于模型計(jì)算，同時(shí)又能真實(shí)模擬水沙運(yùn)動(dòng)現(xiàn)象，特將模型水庫(kù)計(jì)算域上游側(cè)設(shè)置為壓力進(jìn)口邊界，從水庫(kù)由上到下依次為空氣、庫(kù)水和淤泥3種介質(zhì)，按照上述水庫(kù)計(jì)算域的3個(gè)分層，進(jìn)口邊界相應(yīng)分為空氣壓力進(jìn)口、庫(kù)水壓力進(jìn)口和淤泥壓力進(jìn)口三部分；排沙洞出口設(shè)置為壓力出口邊界，壓力按自由出流邊界控制。

4.8 堆渣不沖流速

對(duì)巖塊來說，最常用的不沖流速公式為伊茲巴士公式：

式中：V為不沖平均流速，m/s；D為粒徑，m；g為重力加速度，m/s2；γs，γ 分別為巖塊和水的容重；C 為反映巖塊穩(wěn)定狀況的無量綱系數(shù)，取0.86時(shí)塊石是處于滑動(dòng)臨界失穩(wěn)狀態(tài)，取1.2時(shí)塊石是處于傾覆失穩(wěn)臨界狀態(tài)；K一般為5～7 m1/2/s。

0.6m 堆渣，γs為 2.73 g/cm3，當(dāng)C取 0.86～1.2 時(shí)，由（3）式計(jì)算的不沖流速為3.88～5.42 m/s。說明當(dāng)流速大于3.88 m/s時(shí)，0.6 m堆渣處于滑動(dòng)狀態(tài)；當(dāng)流速大于5.42 m/s時(shí)，0.6 m堆渣處于滾動(dòng)狀態(tài)。20 mm的沙粒，相應(yīng)的不沖流速為0.71～0.99 m/s。

5 計(jì)算成果分析

本文采用Euler多相流模型，建立三維排沙洞水沙非恒定多相流計(jì)算模型，計(jì)算劉家峽排沙洞正常蓄水位、淤泥層不同顆粒組成等條件下，水下巖塞爆破后排沙洞中多相流運(yùn)動(dòng)過程，計(jì)算成果包括不同時(shí)間過程排沙洞內(nèi)各相體積率及其流動(dòng)速度、壓力分布等。

1）對(duì)于不考慮巖塞口區(qū)淤泥層的情況，參照有壓管道流速計(jì)算公式計(jì)算的巖塞爆破后排沙洞內(nèi)水的運(yùn)動(dòng)速度為25.08 m/s。

2）方案1（3相流，不含堆渣)各相流速過程線見圖1（a）。爆破后160 s，水相和淤泥相均流出排沙洞出口，淤泥（1 mm沙粒）相在逆坡段中心的流速為12.7 m/s，是排沙洞內(nèi)的最小流速。

圖1 排沙洞各相流速過程線圖

3）方案2（3相流）各相流速過程線見圖1（b）。爆破后175 s，水相和堆渣相均流出排沙洞出口，但在逆坡段中心處，堆渣（0.6 m堆渣）相的流速為5.13 m/s，是排沙洞內(nèi)的最小流速，小于0.6 m堆渣5.42 m/s的不沖流速，存在0.6 m堆渣停止流動(dòng)的可能。

4）方案3（4相流，不含堆渣）各相流速過程線見圖1（c）。爆破后173 s，水相流至排沙洞全長(zhǎng)97.66%的地方，淤泥相流至排沙洞全長(zhǎng)97.41%的地方，堆渣（20 mm沙粒）相在逆坡段中心的最小流速為7.28 m/s遠(yuǎn)遠(yuǎn)大于20 mm沙粒0.99 m/s的不沖流速。

5）方案4（4相流，淤泥含0.6 m堆渣）各相流速過程線見圖1（d）。巖塞爆破后198 s，在逆坡段中心，水相、堆渣(0.6 m堆渣)相及淤泥(1 mm沙粒)相均達(dá)到最小流速，依次為 5.38，5.08，5.16 m/s，均小于 0.6 m 堆渣 5.42 m/s的不沖流速，存在0.6 m堆渣在逆坡段停止運(yùn)動(dòng)的可能。

[1]計(jì)算流體動(dòng)力學(xué)分析[M].北京：清華大學(xué)出版社，2004.