王希晨，張 健，許聽雨，俞曉東，陳 勝

(1.南水北調(diào)東線江蘇水源有限責任公司，江蘇 南京 210019； 2.河海大學水利水電學院,江蘇 南京 210098)

為提高大型城市防洪排澇能力，針對中心城區(qū)地面硬化率高、建筑密度大、地下市政管線錯綜復雜等特點，美國芝加哥、日本東京、英國泰晤士、法國巴黎以及我國廣州、上海、成都、重慶等城市已開始規(guī)劃建設深層隧道排水調(diào)蓄系統(tǒng)工程(簡稱深隧工程)[1-13]。主隧系統(tǒng)是整個深隧工程中的重要組成部分，由主隧及豎井組成，豎井與主隧水氣相通，相互作用，豎井入流影響主隧流態(tài)，主隧排氣影響豎井入流，主隧從明流到滿流過程中易產(chǎn)生涌浪、滯留氣團、彌合水錘等諸多復雜水力學現(xiàn)象。近年來，國內(nèi)外學者對管(隧)道系統(tǒng)水氣兩相流進行了深入研究：Zhou等[14]對有壓輸水管道中含有兩個氣團的瞬變壓力進行了數(shù)值模擬研究；Trindade等[15]通過數(shù)值模擬和模型試驗對有壓輸水管道考慮空氣增壓的充水過程進行了研究；Catano-Lopera等[16]采用數(shù)值模擬方法對芝加哥深隧系統(tǒng)工程中的瞬變流和水氣兩相流進行了研究，提出不穩(wěn)定工況和水氣相互作用是間歇泉發(fā)生的重要原因；王才歡等[17]通過某水利工程導流隧洞水力模型試驗，指出在一定水流條件下，在較長的洞段及較大的流量區(qū)間，會出現(xiàn)明滿交替流、明流沖擊波、水流折沖等引發(fā)洞壁有害振動或空蝕破壞的不良流態(tài)；郭永鑫等[18]結合南水北調(diào)中線工程北京西環(huán)暗涵充水試驗，建立了長輸水管道充水過程的水氣兩相流數(shù)值模型；吳建華等[19]針對深隧工程中豎井的高落差和大流量等特點提出了一種大轉角階梯泄水道豎井結構，采用物理模型試驗方法研究了流態(tài)、時均動水壓力、消能和摻氣特性等水力特性；王斌等[20]通過物理模型試驗分析了矩形截面折板豎井內(nèi)的流態(tài)及消能特性,提出了改善豎井流態(tài)的整治措施；李璐等[21]采用數(shù)值模型與物理模型相結合的方法對深隧工程的折板式豎井的消能特性進行了研究；夏海等[22]對跌落式入流豎井內(nèi)部流場進行了數(shù)值模擬，研究了不同豎井體型和不同入口流速工況下豎井內(nèi)部壓力場、速度場的變化規(guī)律及其消能率；徐劍喬[23]采用數(shù)值模擬方法對武漢大東湖深隧工程入流豎井進行了建模分析,對比了渦流式豎井和折板式豎井在不同工況下的消能和排氣效果；Wang等[24]針對深隧系統(tǒng)流量大范圍變化工況，采用數(shù)值模型和物理模型相結合的手段，對某對稱結構的新型折板式豎井的水力特性進行了研究。

綜上可知，當前國內(nèi)外學者對管(隧)道系統(tǒng)水氣兩相流研究主要集中在水利工程過水隧洞中的明滿流水力特性方面，對于城市深隧工程的水氣兩相流研究大多集中在單體消能豎井等局部設施，對于大排水流量、深埋地下的深隧工程而言，其包含豎井和相連主隧的主隧系統(tǒng)的水力特性研究仍相當有限。目前，國內(nèi)外對于管道中的明滿流問題研究主要采取一維數(shù)值模擬方法[25-29]，該方法對主隧系統(tǒng)進水過程的模擬難以準確描述主隧系統(tǒng)水氣兩相流復雜的三元結構，以及主隧與豎井之間水氣相互作用的具體流動現(xiàn)象?；诹黧w動力學計算(computational fluid dynamics, CFD)方法[30]的三維數(shù)值模擬在形象地反映局部復雜流態(tài)方面有顯著優(yōu)勢，目前已在深層隧道調(diào)蓄系統(tǒng)的消能豎井等局部流場的數(shù)值模擬中得到了應用[21-24]。該方法對深隧全系統(tǒng)進行模擬，對于研究整個深隧系統(tǒng)水氣運動特性具有重要意義，但三維CFD模擬消耗大量計算資源，現(xiàn)有硬件水平難以滿足深隧全系統(tǒng)三維模擬需求。

為對主隧調(diào)蓄過程中的水力學問題進行研究，筆者所在課題組搭建了包含兩豎井及相連主隧的主隧系統(tǒng)物理模型。但受到試驗條件、試驗方法、量測手段的限制，物理模型方法對于主隧系統(tǒng)全流場的水氣運動現(xiàn)象捕捉有限。筆者采用CFD三維數(shù)值模擬與物理模型試驗相結合的方式，建立包含兩豎井和相連主隧段的主隧系統(tǒng)三維數(shù)值模型，對大流量工況下的進水過程進行三維動態(tài)數(shù)值模擬，旨在揭示主隧系統(tǒng)在大進水流量下的水氣運動及壓力變化規(guī)律，為主隧系統(tǒng)工程結構設計和工況控制提供科學依據(jù)。

1 模型建立與驗證

1.1 數(shù)值模型的建立

建立深隧工程主隧系統(tǒng)三維數(shù)值模型，對其在大流量進水工況下的水氣運動特性進行模擬，并通過筆者所在課題組搭建的深隧系統(tǒng)物理模型試驗平臺對模擬結果進行驗證。數(shù)值模型與物理模型的體型結構尺寸相同，模型與原型比尺為1∶30。

模型包含兩端豎井及一段與之相連的主隧，主隧管徑0.3 m，水平全長51.79 m，坡度0.001，兩豎井間高差0.051 79 m，主隧布置線路如圖1所示。模型采用了如圖2所示的簡化豎井結構，豎井高1 m，由中心內(nèi)井及外井構成，內(nèi)、外井直徑分別為0.4 m和1 m，內(nèi)井底部連接一直徑0.4 m的進水管，頂部開設一直徑0.1 m的排氣孔；內(nèi)井壁面上、下端各開設4個使內(nèi)外井相連通的高0.08 m、間距0.1 m方孔，水流由進水管進入內(nèi)井后從下端方孔進入外井，外井氣體從上端方孔進入內(nèi)井后從豎井頂部排氣孔排出。模型中對豎井采用內(nèi)井底部進流的結構形式，模擬雨水跌落后經(jīng)豎井內(nèi)部消能工充分消能后，在豎井底部匯集并流入主隧的進流狀態(tài)，旨在未確定豎井具體結構形式的情況下有效減小豎井水流對主隧排氣作用的影響，獲得較為良好的主隧進流條件。

圖1 主隧系統(tǒng)模型布置線路

圖2 主隧系統(tǒng)三維數(shù)值模型概化豎井結構

對主隧系統(tǒng)三維數(shù)值模型進行計算網(wǎng)格劃分：豎井內(nèi)井、外井主體部分和主隧段采取結構式網(wǎng)格，外井與主隧端口連接部分采取非結構式網(wǎng)格；對豎井下半部分和主隧兩端進行網(wǎng)格加密。采用經(jīng)網(wǎng)格無關性驗證確定的單元格尺寸進行網(wǎng)格劃分，總網(wǎng)格數(shù)約為295.2萬個。

該主隧系統(tǒng)三維模型包含2個進口邊界和2個出口邊界，其中主隧兩端豎井的進水口給定具有恒定流量的滿管入流，豎井的排氣口設為大氣壓邊界，豎井內(nèi)、外壁和主隧管壁均設為壁面邊界，粗糙度忽略不計?？紤]到模型兩端豎井連通外界大氣，進水過程為明流到明滿流，無瞬變工況，且在水位高于管道端口上方0.2 m時就停止進流，滯留在管道內(nèi)的氣體壓縮程度較小，將氣體視為不可壓縮氣體且無溫度變化。采用VOF方法(流體體積法)對水氣兩相交界面進行追蹤，定義水為主相、氣體為第二相，通過求解質(zhì)量方程每一相的體積分數(shù)來對水氣交界面進行模擬：體積分數(shù)為1則為水；體積分數(shù)為0則視作空氣；體積分數(shù)為0.5視為水氣交界面。由于模型體型復雜且存在彎曲流動，故選擇對彎曲流動具有較高精度的Realizablek-ε湍流模型?；贔LUENT軟件對主隧系統(tǒng)三維模型進行動態(tài)求解。由于模型網(wǎng)格數(shù)大，三維動態(tài)計算需要消耗大量計算資源及計算時間，選取兩典型大流量極端工況進行模擬，計算時間步長設為0.005 s，當豎井1水位超出主隧上表面約0.2 m時停止進流。工況1豎井1進水流量為0.008 7 m3/s，豎井2進水流量為0.013 9 m3/s。工況2兩端豎井進水流量均為0.0162 3 m3/s。

1.2 物理模型試驗驗證

為實時監(jiān)測主隧內(nèi)水氣兩相流沿程壓力變化過程，在主隧底部和頂部沿程布置了一系列測壓點：從豎井一側主隧端口開始，于主隧內(nèi)靠近底部位置沿管線方向每隔10 m布置1個測壓點(記作P1～P6)，其中P1點位于距主隧端口界面0.1 m處。此外，在P1測壓點對應的管頂位置布置P7測壓點，在靠近豎井外壁、主隧端口界面頂部向上0.1 m處布置P8測壓點，在側壓點上布置壓力傳感器，進行同步采集壓力數(shù)據(jù)，見圖3。主隧系統(tǒng)三維數(shù)值模型與物理模型相同位置布置8個測壓點，計算過程中進行實時監(jiān)測壓力，每0.05 s提取1個壓力值，得到壓力隨時間變化的曲線。

圖3 主隧系統(tǒng)三維模型測壓點布置

將主隧系統(tǒng)在工況1下的三維數(shù)值模擬與物理模型試驗結果進行比對，對主隧系統(tǒng)三維數(shù)值模型的可靠性進行驗證。工況1進水條件下，物理模型試驗觀測和三維數(shù)值模擬的進水過程基本相似，主隧系統(tǒng)三維數(shù)值模擬及物理模型試驗部分測壓點上壓力隨時間變化曲線如圖4所示。三維數(shù)值計算結果與物理試驗測量值與在壓力變化趨勢和壓力值上均較為接近，僅在曲線糙率、壓力變化時間點等處存在微小誤差。造成誤差的原因可能包括：三維數(shù)值模型由算法、網(wǎng)格質(zhì)量產(chǎn)生的計算誤差；數(shù)值模擬與物理模型在入流工況控制上產(chǎn)生的誤差；物理模型對壓力數(shù)據(jù)量測和采集精度產(chǎn)生的誤差；將氣體視為不可壓縮氣體帶來的誤差。三維數(shù)值模型結果與試驗量測的壓力曲線變化趨勢上基本接近且誤差在合理范圍內(nèi)，因此建立的主隧系統(tǒng)三維數(shù)值模型是可靠的，可用于研究主隧系統(tǒng)在大進水流量工況下水氣運動及壓力變化規(guī)律。

圖4 工況1主隧系統(tǒng)三維數(shù)值模擬與物理模型試驗部分測壓點壓力變化對比

2 水力特性分析

在前述模型建立的基礎上，對主隧系統(tǒng)兩種工況下的動態(tài)模擬結果進行分析，對進水過程、水氣運動及壓力變化過程進行分析總結。

2.1 進水過程

通過計算對比分析可知，工況1和工況2入流條件下主隧系統(tǒng)進水過程相似，可歸納為4個階段：①第一階段：水流由兩端豎井開始進入主隧。②第二階段：兩端入流在主隧中部附近沖撞產(chǎn)生一對反射液面波，沿主隧中部向兩端推進過程中又產(chǎn)生較弱液面波，多股液面波相互碰撞、疊加，主隧水位迅速上升。③第三階段：主隧下游、上游側端口先后被水流封堵，主隧進流受阻，兩端豎井水位超過主隧頂部并持續(xù)上升，大量氣體滯留于主隧上部，僅少量氣體從兩端口頂部排出。④第四階段：當上游豎井水位超出主隧頂部并繼續(xù)上漲時對進流進行控制；進流停止后，主隧內(nèi)氣體逐漸排出，兩端豎井水位逐步回落。

2.2 水氣運動特性

分別截取兩端豎井包括相連主隧局部段沿主隧入流方向截取的中心平面以及端口附近主隧橫截面，對各截面上進行壓力、速度矢量和水面線流場綜合分析，見圖5～8，進水過程的4個階段水氣運動特性總結如下。

圖5 進水第二階段工況1水氣運動特性

圖6 進水第二階段工況2水氣運動特性

圖7 進水第三階段工況1水氣運動特性

圖8 進水第三階段工況2水氣運動特性

a.第一階段：水流進入主隧初期，流場整體壓力和速度分布相對均勻，兩豎井排氣口速度相近，主隧兩端口處速度場均存在微弱擾動——水流從主隧下部進入，氣體從主隧上部排出。工況1中，豎井2一側主隧端口水位略高于豎井1一側主隧端口水位，而工況2由于兩端進水流量相同，主隧兩側端口水位差異并不明顯。

b.第二階段：當兩股水流在主隧中部碰撞并產(chǎn)生第一對反射波時，反射波向主隧兩端推進，使得主隧端口上部氣流排氣速度增大，氣體沖擊兩端豎井內(nèi)壁。該階段兩種工況水氣運動特征略有差異，工況2較工況1水氣擾動更為劇烈。

工況1：隨著主隧水位上升，豎井2側主隧端口排氣面積減小且低于豎井1側排氣面積，因而該側主隧上部排出氣流速度增大，并進一步影響到兩豎井頂部排氣口速度，使得豎井1排氣口速度明顯減小，豎井2排氣口速度顯著增大。當?shù)谝还煞瓷洳ū平Q井2時，主隧上部氣流對豎井內(nèi)壁沖擊速度增大，沖擊面積減小，主隧上部出現(xiàn)局部低壓，豎井2頂部排氣速度達到最大；而豎井1端口處主隧上部氣流速度減小，豎井1頂部排氣口排氣進一步減弱。當豎井2端口處產(chǎn)生反向液面波并向主隧內(nèi)部傳遞時，主隧內(nèi)上部氣流速度迅速降低，而主隧端口排氣方向開始朝上變化，與內(nèi)壁間存在一局部氣團；同時豎井2頂部排氣速度開始減小，而豎井1頂部排氣速度開始增加。

工況2：受到坡度影響，液面波向豎井2側傳遞較快，因而主隧上部氣流從豎井2側快速排出，豎井2頂部排氣速度增大。當?shù)谝还煞瓷洳ū平Q井2時，主隧端口排出氣流對外井產(chǎn)生渦旋擾動。液面波抵達豎井2，對內(nèi)井壁面產(chǎn)生強烈沖擊，外井及主隧端口水位大幅上升，排氣面積減小，豎井2頂部排氣減弱，豎井1頂部排氣加大。

c.第三階段：豎井2處主隧端口被封堵，少量氣流從端口頂部小區(qū)域并貼向豎井外井壁向上排出，豎井2外井水氣作用明顯，水面受到氣流擾動出現(xiàn)局部壅高，豎井2頂部排氣明顯減?。恢魉韮?nèi)氣流向豎井2側排氣受阻，開始大量排向豎井1一側，豎井1側主隧端口排氣速度劇烈增加，對豎井1內(nèi)井壁形成強烈沖擊，豎井1頂部排氣速度也顯著增大。豎井1處主隧端口很快也被水流封堵，豎井1頂端排氣開始減小，對內(nèi)井上部氣流造成強烈擾動并產(chǎn)生氣旋；豎井2處排氣開始增加。兩端豎井水位開始超過主隧頂部并持續(xù)上升，兩端主隧端口氣流向上排出水面，水面受擾動均產(chǎn)生局部壅高，兩端豎井頂部排氣開始趨于接近。豎井1水位超過主隧頂部約0.2 m，此時主隧內(nèi)部壓力明顯增加，內(nèi)井均存在不同程度的氣旋擾動，由于主隧存在一定坡度，豎井2側主隧端口排出的氣流方向較豎井1更貼向豎井外井壁面。

d.第四階段：停止進流，兩豎井頂部幾乎停止向外排氣，主隧內(nèi)氣體持續(xù)從兩側主隧端口進入兩豎井內(nèi)部，豎井水位逐漸回落，主隧壓力降低。由于坡降因素，豎井1較豎井2水位降低較快，壓力分布也較低。

綜上可知：主隧系統(tǒng)兩種工況下進水過程各階段水氣運動特性相似，但由于工況2較工況1進水流量大，其氣體擾動、排氣速度均較工況1種更為劇烈。兩種極端工況下，主隧氣流從兩主隧端口排出時，會對豎井內(nèi)井壁產(chǎn)生沖擊，易產(chǎn)生空蝕；主隧端口處水氣擾動明顯，當下游側和上游側主隧端口先后被水流封堵時，兩端排氣變化劇烈。主隧端口被水流封堵后，氣體難以順暢排出且進流受阻，大量氣體滯留在主隧頂部，使得主隧的蓄水能力大大減弱。進流較大的工況下主隧端口更易被封堵，滯氣現(xiàn)象更為嚴重。

2.3 壓力變化過程

兩種工況下，進水過程4個階段主隧沿程壓力隨時間變化趨勢相似，見圖9、圖10,其壓力變化過程如下。

圖10 工況2主隧系統(tǒng)三維數(shù)值模型測壓點壓力變化

a.第一階段(工況1，0～45 s；工況2，0～35 s)：水流從兩端口進入主隧，主隧壓力平穩(wěn)上升。

b.第二階段(工況1，45～130 s；工況2，35～90 s)：水流在主隧中部匯集碰撞形成液面波，多股液面波沿主隧在兩豎井間碰撞、疊加，主隧沿程各壓力隨液面波的經(jīng)過出現(xiàn)跳躍式上升。

c.第三階段(工況1，130～170 s；工況2，90～115 s)：水流開始封堵主隧端口，主隧進流受阻，各壓力曲線出現(xiàn)大幅上升，主隧沿程各測壓點壓力值趨于接近；隨著豎井水位超出主隧頂部持續(xù)上升，豎井兩測壓點(P7、P8)壓力也開始大幅增加。

d.第四階段(工況1，170～240 s；工況2，115～200 s)：進流停止，豎井水位回落，壓力緩慢降低。各壓力曲線下降過程中存在小波動，豎井水面回落過程中出現(xiàn)小幅震蕩，影響到主隧內(nèi)部進水排氣，進而對主隧沿程壓力產(chǎn)生波動影響。當豎井水位超過主隧上升階段，主隧壓力迅速上升，但因對進流進行了有效控制，整個充水過程中各測壓點最大壓力值均未超過5 kPa(相當于原型的176 kPa)。

3 結 論

a.主隧系統(tǒng)在大流量下的進水過程可歸納為4個階段：第一階段——水流由兩端豎井開始進入主隧；第二階段——兩端進流在主隧中部附近碰撞產(chǎn)生一對反射液面波，沿主隧向兩端推進過程中，多股液面波相互碰撞疊加，主隧水位迅速上升；第三階段——主隧下游、上游端口先后被水流封堵，兩端豎井水位超過主隧頂高并持續(xù)上升，大量氣體滯留于主隧上部，僅少量氣體從兩側主隧端口頂部排出；第階四段——停止進流，豎井水位波動式逐漸回落。

b.進水流量主要影響主隧系統(tǒng)端口封堵程度、滯氣現(xiàn)象和進水各階段進程時間。流量越大，主隧端口的水流封堵、滯氣現(xiàn)象越嚴重，進水過程各階段進程時間大幅縮短，主隧端口氣體擾動也更為劇烈。

c.依據(jù)兩種大流量工況下主隧系統(tǒng)進水過程水氣運動特性分析，可對主隧系統(tǒng)優(yōu)化設計提出以下建議：①優(yōu)化豎井消能結構以減小主隧進流的摻氣量，增強進流的平穩(wěn)性，提高豎井排氣性；②加固豎井內(nèi)井下部壁面，以增強結構穩(wěn)定性，提升抗空蝕、抗沖擊能力；③改善主隧與豎井的銜接結構，以增強主隧端口過流排氣能力，減少進水對主隧端口的局部封堵，減小主隧內(nèi)部滯氣；④對進水流量進行有效控制，減緩進水流量，減短峰值流量進流時間，當水位高于主隧頂部時停止進流。

d.筆者針對主隧系統(tǒng)物理模型建立三維數(shù)學模型并開展CFD模擬研究，揭示了主隧系統(tǒng)在大流量進水工況下水氣運動及壓力變化規(guī)律。但主隧系統(tǒng)三維數(shù)值模型采用了簡化豎井結構方案，入流條件較為理想化。而實際工程中，經(jīng)豎井消能傳輸后，進入主隧兩端的水流摻氣量和能量較大，主隧端口處的擾動也會更加嚴重，對入流和排氣均會造成阻礙，特別是大流量下主隧內(nèi)滯氣現(xiàn)象會嚴重影響系統(tǒng)結構安全。因此，主隧系統(tǒng)三維數(shù)值模型主隧兩端入流豎井應依據(jù)工程明確的豎井具體設計方案進行建立，以模擬真實入流情況下主隧系統(tǒng)入流過程，為實際工程設計提供具體優(yōu)化方案。