周 睿，程永光，吳家陽

(武漢大學 水資源與水電工程科學國家重點實驗室，湖北 武漢 430072)

水庫水溫時空變化預測是水電工程生態(tài)環(huán)境評價和改善研究領(lǐng)域中重要的研究方向，借助數(shù)學工具和數(shù)學方法來解決類似問題最經(jīng)濟可行[1]。預測水庫水溫，應(yīng)解決水面大幅變動模擬及其計算的問題，可以通過改進浸沒邊界法加以實現(xiàn)。對于水庫水溫類型的判別，我國現(xiàn)行水庫環(huán)境影響評價中普遍采用經(jīng)驗公式法[2]：徑流-庫容比數(shù)法和密度弗勞德數(shù)法。也有研究者提出其他一些方法，如蔡為武[3]考慮到水庫調(diào)節(jié)性能、年內(nèi)泄水狀況、泄水孔口相對位置而提出的水庫水溫分層法；郄志紅等[4]提出的一種判別水庫水溫分層模式的人工神經(jīng)網(wǎng)絡(luò)法；劉金祿[5]根據(jù)水庫水溫垂直結(jié)構(gòu)分布的模糊性特點，應(yīng)用多目標模糊模式識別和回歸分析理論，研究了水庫水溫垂直分布類型的模糊識別回歸預測法。國內(nèi)研究者還提出了許多經(jīng)驗性水溫估算方法[6-8]。

目前國內(nèi)水庫水溫水質(zhì)預測的難點在于對速度場和溫度場耦合作用、水庫形態(tài)等因素的處理，模型通用性不足。而新近發(fā)展起來的格子Boltzmann方法(Lattice Boltzmann Method，LBM)[9]是一種先進的流體力學數(shù)值模擬方法，具有算法簡單、并行性好、邊界條件易處理和壓力直接計算等優(yōu)點[10]。但是將LBM法用于模擬水庫水溫的研究比較有限，且模型預測也存在較大誤差。因此為了完善水溫預測方法，本文利用簡便易行、計算效率高的浸沒邊界法，開展水庫水溫的時空變化規(guī)律研究，模擬真實水庫自由液面大幅度變動[11]，并通過簡單算例來證實其可靠性。

1 自由液面處理算法介紹

1.1 控制方程

首先引入階梯函數(shù)即Heaviside函數(shù)H，規(guī)定在一種流體(比如液體)中H為1，在氣體中H為0，自由液面則標記H為0.5。為將Heaviside函數(shù)與Dirac函數(shù)聯(lián)系起來，將它表示為一維Dirac函數(shù)乘積的積分，

例如在二維條件下積分區(qū)域是由周線S包圍的A(t)面積。

很顯然，當(x，y)落在面積A(t)內(nèi)時，Heaviside函數(shù)取值為1而在其他地方取值為0。為了計算Heaviside函數(shù)的梯度，設(shè)定兩類變量，以有上標的變量為第二類變量。首先考慮梯度是針對第一類變量計算，而積分是針對第二類變量計算，這樣梯度符號可以和積分號交換順序。同時Dirac對于兩類變量是反互易的，所以梯度符號也可以針對第二類變量計算，即

本模型考慮的是同一種不可壓縮流體，密度為常數(shù)，則全場密度可用Heaviside函數(shù)表示

式中：ρ1和ρ0分別為H取1和0處的密度。而密度的梯度則可以表示為

在此▽ρ=ρ1-ρ0。

流體的運動假定由N-S方程描述，在考慮表面張力的情況下，有

上述方程適用于全流場(含密度場)。β為所涉及問題的維數(shù)；κ為自由液面曲率；σ為表面張力系數(shù)；δβ是β維的Dirac函數(shù)，由一維Dirac函數(shù)的乘積構(gòu)成；n為自由液面的法向量。式(5)最右端積分項的積分區(qū)間為自由液面。模型盡管可以考慮流體黏性系數(shù)的不連續(xù)性，但依然假定流體黏度是常數(shù)，即

自由面邊界條件獲取途徑為：在自由液面附近取一小塊圓柱體，圓柱體上表面在氣相中，而下表面在液相中，分別計算作用在圓柱體表面上的面力及圓柱體內(nèi)部的體積力。根據(jù)達朗貝爾原理，考慮到體積力中慣性力作用，所以各力應(yīng)該是平衡的。將圓柱體體積縮小并逐漸趨于零，由于體積力是較面積力的高階小量，故可忽略所有體積力項。將方程在所考察點附近沿自由液面法線方向的平行和垂直方向分解，可得應(yīng)力張量的垂直分量在自由液面兩側(cè)連續(xù)，而平行于法線方向的分量則應(yīng)該滿足

式中：[]表示相關(guān)變量的階躍值。速度矢量在自由液面兩側(cè)也是連續(xù)的。

1.2 N-S方程求解

求解攜帶外力項的不可壓縮流體的動量方程和連續(xù)性方程。如果底層網(wǎng)格是結(jié)構(gòu)化正交網(wǎng)格，則采用projection method。首先重構(gòu)全場密度：

式中：n為時間步數(shù)。

動量方程在固定的結(jié)構(gòu)化正交網(wǎng)格上加以求解，而表面張力可以通過變動水面上的網(wǎng)格來計算，所以需要將表面張力擴散至變動水面周圍的流體上。因界面由Dirac函數(shù)定義，所以這種擴散必然包含多維Dirac函數(shù)離散形式的構(gòu)造。這種構(gòu)造可從多方面進行，會有不同結(jié)果，但其最終形式必須滿足一定的要求[12]，如：

為更新自由液面位置，在界面網(wǎng)格上定義真實的界面移動速度矢量。由于自由液面由流體組成且速度場在界面兩側(cè)連續(xù)，所以液面按照當?shù)亓黧w速度運動，即滿足無滑移邊界條件。為計算定義在界面網(wǎng)格上的速度，利用多維Dirac函數(shù)由流體速度插值得到界面移動速度，即

式中：小寫和大寫字母分別代表定義在流體結(jié)構(gòu)化正交網(wǎng)格和界面網(wǎng)格上的變量。得到自由液面速度后，其位置更新可按下式進行：

式(11)的離散格式很多。離散格式關(guān)乎到浸沒邊界法的穩(wěn)定性。采用恰當?shù)碾x散格式能帶來無條件穩(wěn)定的優(yōu)良性質(zhì)?？刹捎蔑@式向前差分格式。

流體的物性系數(shù)，如密度、黏度等，在一種相的流體中取值為常數(shù)，而在界面處存在階躍值，需要在每一時刻利用自由液面位置，來確定全場的密度或黏度分布。

自由液面處密度場的梯度已知，而其他地方梯度值均為零，進而根據(jù)邊界密度逐步更新全場密度值。但當界面卷曲較嚴重時，該方法不再采用，可以換用更為一般化的方法。由于式(4)解決了計算全場的密度梯度問題，所以可在其基礎(chǔ)上繼續(xù)取散度得到密度場的拉普拉斯算符：

式(12)右端由傳統(tǒng)的中心差分格式計算。求解由此得到的泊松方程，從而得到全場密度。在三維條件下，為加速計算過程，考慮到密度梯度只是在界面附近才取非零值，所以將泊松方程的計算區(qū)域局限在界面附近，而邊界均采用狄利克雷邊界條件。采用格子Boltzmann方法來求解泊松方程，具體過程參考文獻[13]。

表面張力計算參考拉普拉斯公式：

式中：f為表面張力面密度；σ為表面張力系數(shù)；R1和R2分別為曲面內(nèi)過曲面上一點的曲線最大和最小曲率，即主曲率。表面張力的方向嚴格沿曲面上一點的法線方向。在離散網(wǎng)格上利用式(13)，采用三角形網(wǎng)格來離散自由液面：

式中：k(xl)為坐標為xl的點 Nl的平均主曲率的2倍，而定義在點Nl上的面積是泰森多邊形的面積：

式中：點 q 是與點 l 共享一個三角形的另外兩個節(jié)點；α和β為這些節(jié)點在三角形內(nèi)的內(nèi)角度數(shù)。從而法向平均主曲率可按下式計算：

法向平均主曲率乘以表面張力系數(shù)即可得到張力的面密度。采用平面三角形來離散曲面，即網(wǎng)格只能通過加密才能逐漸逼近真實的曲面；而在實際中無法還原曲面，也可采用曲面三角形來離散液面，按照以下方法來計算曲面上某一點的平均法向主曲率：

則作用在一小塊三角形上的表面張力為：

在自由液面的模擬中，液面與計算區(qū)域的邊界相交，從上面計算某一拉格朗日點的表面張力來看，該點需要完全被其他點所包圍，而這一條件對于處在邊界上的拉格朗日點來說是無法滿足的。可以在邊界外面再新增一排節(jié)點，這排節(jié)點并不會像內(nèi)部節(jié)點一樣參與拉氏力的計算及其擴散和位置的更新，其坐標根據(jù)具體所設(shè)置的邊界條件人為地給定，比如對于滑移固壁、對稱邊界和周期邊界來說，液面在邊界處的法向量應(yīng)該與邊界平行。

模擬中速度場的局部差異性較大。隨著計算的進行，網(wǎng)格質(zhì)量開始下降，為使網(wǎng)格能夠分布均勻，且三角形網(wǎng)格的邊長能夠小于預先設(shè)定好的上限值，存在下面3種情形需要處理(圖1)。

預先設(shè)定三角形邊長的取值范圍，當其邊長小于下限時，利用圖1(a)所示方法來添加三角形；大于上限時，利用圖1(b)所示方法刪除三角形；特殊情形需重新劃分三角形(圖1(c))。上述網(wǎng)格重構(gòu)可借助C語言中鏈表容器實現(xiàn)。

圖 1 液面三角形網(wǎng)格重構(gòu)Fig. 1 Reconstruction of liquid surface triangle meshes

2 水箱底部注水液面變動算例

2.1 計算條件說明

為驗證上述模型的可行性和準確性，可嘗試模擬包含自由液面和熱傳導的流場。為便于表述，以下物理量均表述為格子單位，時間以計算時步作為單位。計算區(qū)域仍然為正立方體，均勻離散成60×60×60的網(wǎng)格，流體運動黏度υ=0.1，熱擴散系數(shù)D=0.1，自由液面表面張力系數(shù)σ=0.7。如圖2所示，初始時刻自由液面處于平衡位置，其上部為氣體，密度ρg=0.1，其下方為液體，密度ρf=1.0，全場流體處于靜止狀態(tài)，溫度統(tǒng)一賦值為T0=1.0。計算區(qū)域的頂部為壓力邊界和溫度邊界，四周為周期邊界，底部為壁面。在底面中央開一半徑為10的小孔，從里面以豎向速度u=0.2注入溫度為2.0T0的熱流體，當計算進行到t=1 100后，將此小孔改為壓力出口邊界。預期在熱流體的注入下，自由液面逐漸上升的同時會伴隨一定波動；而當將速度入口改為壓力出口后，液體將在內(nèi)外壓差作用下開始泄流，而自由液面則會逐漸下降。

2.2 自由液面模擬結(jié)果及分析

為清楚顯示自由液面形態(tài)，單獨將液面與通過計算區(qū)域中心且垂直于底面的平面交線繪制在圖3中，給出的是幾個不同時刻自由液面的位置和形態(tài)。

模擬前期，自由液面初始位置靠近底面，注入的液體來不及擴散，分布極不均勻；而在注入液體的推動下液面開始上漲，因而在中央部位向上鼓起。液面曲率越大，表面張力就越大，并且在距離小孔較近的中央部位液面上漲速度較四周快，液面向上鼓起成山包狀(圖3(a))，最終表面張力和液面中央部位的上下壓差平衡，流體在表面張力作用下開始減速，并向四周擴散，液面持續(xù)上漲，之前形成的突起逐漸消失(圖3(b))。由于立方體四周設(shè)置為周期邊界，在邊界上流速只可能沿豎向，所以速度場在遇到側(cè)邊界之前再次轉(zhuǎn)向垂直，進而導致液面周圍部分開始上漲。該突起的消失并非一次性單調(diào)完成，而是以在液面總體向上平移的基礎(chǔ)上疊加一振幅逐漸衰減的波動形式完成。不久突起完全消失，液面基本上成一平面(圖3(c))，當液面上漲到一定高度后，底面注入的液體在撞擊到液面之前已經(jīng)擴散得較為充分，在計算區(qū)域的水平截面上垂向速度分布得較為均勻，因此，液面不同部位的上漲速度差別不大，所以才會呈現(xiàn)出平面狀(圖3(c)。液面在停止注入液體的t=1 100時達到最高點(圖3(d))，在將速度入口邊界修改為壓力出口邊界后，液體開始自小孔流出，液面逐漸下降(圖3(e))，下降前期液面依然成一平面，降至距離出口較近時，液面中央部位開始向下凹陷，形成一漏斗(圖3(f))。

為分析液面的波動特征，將測點布置在液面正中央，并統(tǒng)計了該點縱坐標變化情況，結(jié)果如圖4所示。在注水前期t＜100時，液面中央上升速度非常快，這是突起形成階段。之后曲線在100＜t＜200時趨于水平，這是由于液面突起達到最高點后向下回彈所致。這段曲線體現(xiàn)該波動疊加在液面總體上漲趨勢，液面上升趨勢和突起向下的波動疊加在一起相互抵消，所以曲線成水平狀。波動過程一直持續(xù)，波動幅值同時也按照某種規(guī)律衰減。

圖 2 初始時刻液面位置和密度分布Fig. 2 Liquid level position and density distribution at initial time

圖 3 不同時刻的液面位置Fig. 3 Liquid level position at different times

圖 4 液面中央測點縱坐標隨時間變化規(guī)律Fig. 4 Variation in vertical coordinates with time at central measuring point of liquid level

為定量分析模擬準確性，圖5給出了液體所占據(jù)的體積隨時間的變化。由于在注水期間入口速度均勻分布，流體不可壓縮，根據(jù)連續(xù)性定理，液面以下液體體積的增長速度應(yīng)該等于入口流量。同理在放水階段液面下降的速度也應(yīng)該等于出口的放水流量。從圖5中可以看出，雖然在注水階段液體體積的增長曲線并非嚴格成一條直線，但與理論直線的差距處于可接受范圍。在放水階段由于給定的是壓力邊界，出口速度并非均勻分布，所以理論曲線也并非是一條直線，得到的結(jié)果與理論解在總體趨勢一致的基礎(chǔ)上存在一定偏差。究其上述偏差的來源，可能一方面來自自由液面的離散處理，由于采用的是三角形平面面元來離散液面，實際上的液面應(yīng)該是光滑曲面；另一方面來自自由液面與計算區(qū)域邊界相交面的處理。前文已經(jīng)闡述了如何處理液面與物理邊界相交的問題，但實際上卻無法完全做到。

2.3 溫度場模擬結(jié)果及分析

為了比較不同擴散系數(shù)對模擬結(jié)果的影響，在保持其他參數(shù)不變條件下改變擴散系數(shù)，圖6為對應(yīng)于D=0.005和0.5的溫度場。

圖 5 液體體積的時間變化規(guī)律Fig. 5 Temporal variation law of liquid volume

圖 6 擴散系數(shù)D=0.005和0.5時不同時刻的溫度場云圖Fig. 6 Temperature field cloud diagram at different times under diffusion coefficients D=0. 005 and D=0.5

擴散系數(shù)表示擴散強度與對流強度的比值，取值越小說明對流作用越明顯。比如對于液面上升至頂端的t=1 100時刻，當擴散系數(shù)取值較小時，D=0.005，溫度場受到對流控制，其分布僅局限在計算區(qū)域中央的狹長地帶，并不像當D=0.5時分布得比較均勻。為定量分析溫度場分布隨時間以及擴散系數(shù)的變化，圖7給出了不同擴散系數(shù)下計算區(qū)域中心測點溫度的變化過程曲線。

溫度變化來自于對流以及擴散兩方面。液面的存在阻止了注入液體進一步向氣體方向侵入，使溫度無法僅通過對流盡快從液體區(qū)域向氣體區(qū)域遷移；但是擴散作用則不受流場分布的影響，僅取決于溫度場的局部差異。對于擴散系數(shù)取值較大的情況(圖7)，溫度分布主要受到擴散作用的控制，即使自由液面還沒有上漲到中心測點，溫度場不需要依靠對流機制仍然可以憑借擴散影響到中央測點，所以在圖7中可以看到測點溫度很快就有所上升。而對于擴散系數(shù)取值較小的情況，擴散作用影響有限，溫度只能等到液面上漲到測點位置后憑借對流作用逐漸上升，所以測點溫度在t=420液面抵達時才開始上升，而此前保持在初始值左右。不管抵達的速度快慢如何，測點溫度一旦抵達其最大值附近便不再繼續(xù)增加，而是在其最大值附近波動，波動原因是自由液面在上升過程中的幾何波動。對于泄水過程，孔口由速度入口邊界修改為壓力出口邊界后，因考慮的是不可壓縮流體，壓力波的傳播速度為無窮大，所以立刻就可以觀察到降壓波。該降壓波疊加在之前液面波動所導致的波動之上。對于擴散系數(shù)比較小的情況，降壓波均出現(xiàn)在溫度下降起始點之前，且擴散系數(shù)越小，溫度下降起始點出現(xiàn)的時間越遲，這是因為溫度下降主要受控于擴散作用。在本次模擬中，頂面邊界的溫度是給定的，并且其取值小于注入液體的溫度，所以擴散系數(shù)越小，測點溫度受到來自于頂面邊界條件影響所需要的時間也就越長，而對于擴散系數(shù)較大的情形，擴散作用明顯，測點溫度下降與降壓波幾乎同時產(chǎn)生。

圖 7 不同熱擴散系數(shù)下區(qū)域中心測點溫度變化Fig. 7 Temperature change curves of central measuring point of calculating region under different thermal diffusion coefficients

3 結(jié) 語

利用浸沒邊界法，在格子Boltzmann雙分布函數(shù)模型基礎(chǔ)上添加了能夠考慮自由液面大幅波動的模塊。通過簡單算例驗證了該方法的可行性和準確性，可在一定精度要求下較好地重現(xiàn)自由液面變動過程，為該方法在今后的應(yīng)用奠定了基礎(chǔ)。