周 睿,程永光,吳家陽
(武漢大學 水資源與水電工程科學國家重點實驗室,湖北 武漢 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],并通過簡單算例來證實其可靠性。
首先引入階梯函數(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ù)的。
求解攜帶外力項的不可壓縮流體的動量方程和連續(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
為驗證上述模型的可行性和準確性,可嘗試模擬包含自由液面和熱傳導的流場。為便于表述,以下物理量均表述為格子單位,時間以計算時步作為單位。計算區(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)外壓差作用下開始泄流,而自由液面則會逐漸下降。
為清楚顯示自由液面形態(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)闡述了如何處理液面與物理邊界相交的問題,但實際上卻無法完全做到。
為了比較不同擴散系數(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
利用浸沒邊界法,在格子Boltzmann雙分布函數(shù)模型基礎(chǔ)上添加了能夠考慮自由液面大幅波動的模塊。通過簡單算例驗證了該方法的可行性和準確性,可在一定精度要求下較好地重現(xiàn)自由液面變動過程,為該方法在今后的應(yīng)用奠定了基礎(chǔ)。