周 睿,程永光,吳家陽
(武漢大學(xué) 水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室,湖北 武漢 430072)
水庫水溫時(shí)空變化預(yù)測是水電工程生態(tài)環(huán)境評價(jià)和改善研究領(lǐng)域中重要的研究方向,借助數(shù)學(xué)工具和數(shù)學(xué)方法來解決類似問題最經(jīng)濟(jì)可行[1]。預(yù)測水庫水溫,應(yīng)解決水面大幅變動模擬及其計(jì)算的問題,可以通過改進(jìn)浸沒邊界法加以實(shí)現(xiàn)。對于水庫水溫類型的判別,我國現(xiàn)行水庫環(huán)境影響評價(jià)中普遍采用經(jīng)驗(yàn)公式法[2]:徑流-庫容比數(shù)法和密度弗勞德數(shù)法。也有研究者提出其他一些方法,如蔡為武[3]考慮到水庫調(diào)節(jié)性能、年內(nèi)泄水狀況、泄水孔口相對位置而提出的水庫水溫分層法;郄志紅等[4]提出的一種判別水庫水溫分層模式的人工神經(jīng)網(wǎng)絡(luò)法;劉金祿[5]根據(jù)水庫水溫垂直結(jié)構(gòu)分布的模糊性特點(diǎn),應(yīng)用多目標(biāo)模糊模式識別和回歸分析理論,研究了水庫水溫垂直分布類型的模糊識別回歸預(yù)測法。國內(nèi)研究者還提出了許多經(jīng)驗(yàn)性水溫估算方法[6-8]。
目前國內(nèi)水庫水溫水質(zhì)預(yù)測的難點(diǎn)在于對速度場和溫度場耦合作用、水庫形態(tài)等因素的處理,模型通用性不足。而新近發(fā)展起來的格子Boltzmann方法(Lattice Boltzmann Method,LBM)[9]是一種先進(jìn)的流體力學(xué)數(shù)值模擬方法,具有算法簡單、并行性好、邊界條件易處理和壓力直接計(jì)算等優(yōu)點(diǎn)[10]。但是將LBM法用于模擬水庫水溫的研究比較有限,且模型預(yù)測也存在較大誤差。因此為了完善水溫預(yù)測方法,本文利用簡便易行、計(jì)算效率高的浸沒邊界法,開展水庫水溫的時(shí)空變化規(guī)律研究,模擬真實(shí)水庫自由液面大幅度變動[11],并通過簡單算例來證實(shí)其可靠性。
首先引入階梯函數(shù)即Heaviside函數(shù)H,規(guī)定在一種流體(比如液體)中H為1,在氣體中H為0,自由液面則標(biāo)記H為0.5。為將Heaviside函數(shù)與Dirac函數(shù)聯(lián)系起來,將它表示為一維Dirac函數(shù)乘積的積分,
例如在二維條件下積分區(qū)域是由周線S包圍的A(t)面積。
很顯然,當(dāng)(x,y)落在面積A(t)內(nèi)時(shí),Heaviside函數(shù)取值為1而在其他地方取值為0。為了計(jì)算Heaviside函數(shù)的梯度,設(shè)定兩類變量,以有上標(biāo)的變量為第二類變量。首先考慮梯度是針對第一類變量計(jì)算,而積分是針對第二類變量計(jì)算,這樣梯度符號可以和積分號交換順序。同時(shí)Dirac對于兩類變量是反互易的,所以梯度符號也可以針對第二類變量計(jì)算,即
本模型考慮的是同一種不可壓縮流體,密度為常數(shù),則全場密度可用Heaviside函數(shù)表示
式中:ρ1和ρ0分別為H取1和0處的密度。而密度的梯度則可以表示為
在此▽ρ=ρ1-ρ0。
流體的運(yùn)動假定由N-S方程描述,在考慮表面張力的情況下,有
上述方程適用于全流場(含密度場)。β為所涉及問題的維數(shù);κ為自由液面曲率;σ為表面張力系數(shù);δβ是β維的Dirac函數(shù),由一維Dirac函數(shù)的乘積構(gòu)成;n為自由液面的法向量。式(5)最右端積分項(xiàng)的積分區(qū)間為自由液面。模型盡管可以考慮流體黏性系數(shù)的不連續(xù)性,但依然假定流體黏度是常數(shù),即
自由面邊界條件獲取途徑為:在自由液面附近取一小塊圓柱體,圓柱體上表面在氣相中,而下表面在液相中,分別計(jì)算作用在圓柱體表面上的面力及圓柱體內(nèi)部的體積力。根據(jù)達(dá)朗貝爾原理,考慮到體積力中慣性力作用,所以各力應(yīng)該是平衡的。將圓柱體體積縮小并逐漸趨于零,由于體積力是較面積力的高階小量,故可忽略所有體積力項(xiàng)。將方程在所考察點(diǎn)附近沿自由液面法線方向的平行和垂直方向分解,可得應(yīng)力張量的垂直分量在自由液面兩側(cè)連續(xù),而平行于法線方向的分量則應(yīng)該滿足
式中:[]表示相關(guān)變量的階躍值。速度矢量在自由液面兩側(cè)也是連續(xù)的。
求解攜帶外力項(xiàng)的不可壓縮流體的動量方程和連續(xù)性方程。如果底層網(wǎng)格是結(jié)構(gòu)化正交網(wǎng)格,則采用projection method。首先重構(gòu)全場密度:
式中:n為時(shí)間步數(shù)。
動量方程在固定的結(jié)構(gòu)化正交網(wǎng)格上加以求解,而表面張力可以通過變動水面上的網(wǎng)格來計(jì)算,所以需要將表面張力擴(kuò)散至變動水面周圍的流體上。因界面由Dirac函數(shù)定義,所以這種擴(kuò)散必然包含多維Dirac函數(shù)離散形式的構(gòu)造。這種構(gòu)造可從多方面進(jìn)行,會有不同結(jié)果,但其最終形式必須滿足一定的要求[12],如:
為更新自由液面位置,在界面網(wǎng)格上定義真實(shí)的界面移動速度矢量。由于自由液面由流體組成且速度場在界面兩側(cè)連續(xù),所以液面按照當(dāng)?shù)亓黧w速度運(yùn)動,即滿足無滑移邊界條件。為計(jì)算定義在界面網(wǎng)格上的速度,利用多維Dirac函數(shù)由流體速度插值得到界面移動速度,即
式中:小寫和大寫字母分別代表定義在流體結(jié)構(gòu)化正交網(wǎng)格和界面網(wǎng)格上的變量。得到自由液面速度后,其位置更新可按下式進(jìn)行:
式(11)的離散格式很多。離散格式關(guān)乎到浸沒邊界法的穩(wěn)定性。采用恰當(dāng)?shù)碾x散格式能帶來無條件穩(wěn)定的優(yōu)良性質(zhì)??刹捎蔑@式向前差分格式。
流體的物性系數(shù),如密度、黏度等,在一種相的流體中取值為常數(shù),而在界面處存在階躍值,需要在每一時(shí)刻利用自由液面位置,來確定全場的密度或黏度分布。
自由液面處密度場的梯度已知,而其他地方梯度值均為零,進(jìn)而根據(jù)邊界密度逐步更新全場密度值。但當(dāng)界面卷曲較嚴(yán)重時(shí),該方法不再采用,可以換用更為一般化的方法。由于式(4)解決了計(jì)算全場的密度梯度問題,所以可在其基礎(chǔ)上繼續(xù)取散度得到密度場的拉普拉斯算符:
式(12)右端由傳統(tǒng)的中心差分格式計(jì)算。求解由此得到的泊松方程,從而得到全場密度。在三維條件下,為加速計(jì)算過程,考慮到密度梯度只是在界面附近才取非零值,所以將泊松方程的計(jì)算區(qū)域局限在界面附近,而邊界均采用狄利克雷邊界條件。采用格子Boltzmann方法來求解泊松方程,具體過程參考文獻(xiàn)[13]。
表面張力計(jì)算參考拉普拉斯公式:
式中:f為表面張力面密度;σ為表面張力系數(shù);R1和R2分別為曲面內(nèi)過曲面上一點(diǎn)的曲線最大和最小曲率,即主曲率。表面張力的方向嚴(yán)格沿曲面上一點(diǎn)的法線方向。在離散網(wǎng)格上利用式(13),采用三角形網(wǎng)格來離散自由液面:
式中:k(xl)為坐標(biāo)為xl的點(diǎn) Nl的平均主曲率的2倍,而定義在點(diǎn)Nl上的面積是泰森多邊形的面積:
式中:點(diǎn) q 是與點(diǎn) l 共享一個(gè)三角形的另外兩個(gè)節(jié)點(diǎn);α和β為這些節(jié)點(diǎn)在三角形內(nèi)的內(nèi)角度數(shù)。從而法向平均主曲率可按下式計(jì)算:
法向平均主曲率乘以表面張力系數(shù)即可得到張力的面密度。采用平面三角形來離散曲面,即網(wǎng)格只能通過加密才能逐漸逼近真實(shí)的曲面;而在實(shí)際中無法還原曲面,也可采用曲面三角形來離散液面,按照以下方法來計(jì)算曲面上某一點(diǎn)的平均法向主曲率:
則作用在一小塊三角形上的表面張力為:
在自由液面的模擬中,液面與計(jì)算區(qū)域的邊界相交,從上面計(jì)算某一拉格朗日點(diǎn)的表面張力來看,該點(diǎn)需要完全被其他點(diǎn)所包圍,而這一條件對于處在邊界上的拉格朗日點(diǎn)來說是無法滿足的??梢栽谶吔缤饷嬖傩略鲆慌殴?jié)點(diǎn),這排節(jié)點(diǎn)并不會像內(nèi)部節(jié)點(diǎn)一樣參與拉氏力的計(jì)算及其擴(kuò)散和位置的更新,其坐標(biāo)根據(jù)具體所設(shè)置的邊界條件人為地給定,比如對于滑移固壁、對稱邊界和周期邊界來說,液面在邊界處的法向量應(yīng)該與邊界平行。
模擬中速度場的局部差異性較大。隨著計(jì)算的進(jìn)行,網(wǎng)格質(zhì)量開始下降,為使網(wǎng)格能夠分布均勻,且三角形網(wǎng)格的邊長能夠小于預(yù)先設(shè)定好的上限值,存在下面3種情形需要處理(圖1)。
預(yù)先設(shè)定三角形邊長的取值范圍,當(dāng)其邊長小于下限時(shí),利用圖1(a)所示方法來添加三角形;大于上限時(shí),利用圖1(b)所示方法刪除三角形;特殊情形需重新劃分三角形(圖1(c))。上述網(wǎng)格重構(gòu)可借助C語言中鏈表容器實(shí)現(xiàn)。
圖 1 液面三角形網(wǎng)格重構(gòu)Fig. 1 Reconstruction of liquid surface triangle meshes
為驗(yàn)證上述模型的可行性和準(zhǔn)確性,可嘗試模擬包含自由液面和熱傳導(dǎo)的流場。為便于表述,以下物理量均表述為格子單位,時(shí)間以計(jì)算時(shí)步作為單位。計(jì)算區(qū)域仍然為正立方體,均勻離散成60×60×60的網(wǎng)格,流體運(yùn)動黏度υ=0.1,熱擴(kuò)散系數(shù)D=0.1,自由液面表面張力系數(shù)σ=0.7。如圖2所示,初始時(shí)刻自由液面處于平衡位置,其上部為氣體,密度ρg=0.1,其下方為液體,密度ρf=1.0,全場流體處于靜止?fàn)顟B(tài),溫度統(tǒng)一賦值為T0=1.0。計(jì)算區(qū)域的頂部為壓力邊界和溫度邊界,四周為周期邊界,底部為壁面。在底面中央開一半徑為10的小孔,從里面以豎向速度u=0.2注入溫度為2.0T0的熱流體,當(dāng)計(jì)算進(jìn)行到t=1 100后,將此小孔改為壓力出口邊界。預(yù)期在熱流體的注入下,自由液面逐漸上升的同時(shí)會伴隨一定波動;而當(dāng)將速度入口改為壓力出口后,液體將在內(nèi)外壓差作用下開始泄流,而自由液面則會逐漸下降。
為清楚顯示自由液面形態(tài),單獨(dú)將液面與通過計(jì)算區(qū)域中心且垂直于底面的平面交線繪制在圖3中,給出的是幾個(gè)不同時(shí)刻自由液面的位置和形態(tài)。
模擬前期,自由液面初始位置靠近底面,注入的液體來不及擴(kuò)散,分布極不均勻;而在注入液體的推動下液面開始上漲,因而在中央部位向上鼓起。液面曲率越大,表面張力就越大,并且在距離小孔較近的中央部位液面上漲速度較四周快,液面向上鼓起成山包狀(圖3(a)),最終表面張力和液面中央部位的上下壓差平衡,流體在表面張力作用下開始減速,并向四周擴(kuò)散,液面持續(xù)上漲,之前形成的突起逐漸消失(圖3(b))。由于立方體四周設(shè)置為周期邊界,在邊界上流速只可能沿豎向,所以速度場在遇到側(cè)邊界之前再次轉(zhuǎn)向垂直,進(jìn)而導(dǎo)致液面周圍部分開始上漲。該突起的消失并非一次性單調(diào)完成,而是以在液面總體向上平移的基礎(chǔ)上疊加一振幅逐漸衰減的波動形式完成。不久突起完全消失,液面基本上成一平面(圖3(c)),當(dāng)液面上漲到一定高度后,底面注入的液體在撞擊到液面之前已經(jīng)擴(kuò)散得較為充分,在計(jì)算區(qū)域的水平截面上垂向速度分布得較為均勻,因此,液面不同部位的上漲速度差別不大,所以才會呈現(xiàn)出平面狀(圖3(c)。液面在停止注入液體的t=1 100時(shí)達(dá)到最高點(diǎn)(圖3(d)),在將速度入口邊界修改為壓力出口邊界后,液體開始自小孔流出,液面逐漸下降(圖3(e)),下降前期液面依然成一平面,降至距離出口較近時(shí),液面中央部位開始向下凹陷,形成一漏斗(圖3(f))。
為分析液面的波動特征,將測點(diǎn)布置在液面正中央,并統(tǒng)計(jì)了該點(diǎn)縱坐標(biāo)變化情況,結(jié)果如圖4所示。在注水前期t<100時(shí),液面中央上升速度非??欤@是突起形成階段。之后曲線在100<t<200時(shí)趨于水平,這是由于液面突起達(dá)到最高點(diǎn)后向下回彈所致。這段曲線體現(xiàn)該波動疊加在液面總體上漲趨勢,液面上升趨勢和突起向下的波動疊加在一起相互抵消,所以曲線成水平狀。波動過程一直持續(xù),波動幅值同時(shí)也按照某種規(guī)律衰減。
圖 2 初始時(shí)刻液面位置和密度分布Fig. 2 Liquid level position and density distribution at initial time
圖 3 不同時(shí)刻的液面位置Fig. 3 Liquid level position at different times
圖 4 液面中央測點(diǎn)縱坐標(biāo)隨時(shí)間變化規(guī)律Fig. 4 Variation in vertical coordinates with time at central measuring point of liquid level
為定量分析模擬準(zhǔn)確性,圖5給出了液體所占據(jù)的體積隨時(shí)間的變化。由于在注水期間入口速度均勻分布,流體不可壓縮,根據(jù)連續(xù)性定理,液面以下液體體積的增長速度應(yīng)該等于入口流量。同理在放水階段液面下降的速度也應(yīng)該等于出口的放水流量。從圖5中可以看出,雖然在注水階段液體體積的增長曲線并非嚴(yán)格成一條直線,但與理論直線的差距處于可接受范圍。在放水階段由于給定的是壓力邊界,出口速度并非均勻分布,所以理論曲線也并非是一條直線,得到的結(jié)果與理論解在總體趨勢一致的基礎(chǔ)上存在一定偏差。究其上述偏差的來源,可能一方面來自自由液面的離散處理,由于采用的是三角形平面面元來離散液面,實(shí)際上的液面應(yīng)該是光滑曲面;另一方面來自自由液面與計(jì)算區(qū)域邊界相交面的處理。前文已經(jīng)闡述了如何處理液面與物理邊界相交的問題,但實(shí)際上卻無法完全做到。
為了比較不同擴(kuò)散系數(shù)對模擬結(jié)果的影響,在保持其他參數(shù)不變條件下改變擴(kuò)散系數(shù),圖6為對應(yīng)于D=0.005和0.5的溫度場。
圖 5 液體體積的時(shí)間變化規(guī)律Fig. 5 Temporal variation law of liquid volume
圖 6 擴(kuò)散系數(shù)D=0.005和0.5時(shí)不同時(shí)刻的溫度場云圖Fig. 6 Temperature field cloud diagram at different times under diffusion coefficients D=0. 005 and D=0.5
擴(kuò)散系數(shù)表示擴(kuò)散強(qiáng)度與對流強(qiáng)度的比值,取值越小說明對流作用越明顯。比如對于液面上升至頂端的t=1 100時(shí)刻,當(dāng)擴(kuò)散系數(shù)取值較小時(shí),D=0.005,溫度場受到對流控制,其分布僅局限在計(jì)算區(qū)域中央的狹長地帶,并不像當(dāng)D=0.5時(shí)分布得比較均勻。為定量分析溫度場分布隨時(shí)間以及擴(kuò)散系數(shù)的變化,圖7給出了不同擴(kuò)散系數(shù)下計(jì)算區(qū)域中心測點(diǎn)溫度的變化過程曲線。
溫度變化來自于對流以及擴(kuò)散兩方面。液面的存在阻止了注入液體進(jìn)一步向氣體方向侵入,使溫度無法僅通過對流盡快從液體區(qū)域向氣體區(qū)域遷移;但是擴(kuò)散作用則不受流場分布的影響,僅取決于溫度場的局部差異。對于擴(kuò)散系數(shù)取值較大的情況(圖7),溫度分布主要受到擴(kuò)散作用的控制,即使自由液面還沒有上漲到中心測點(diǎn),溫度場不需要依靠對流機(jī)制仍然可以憑借擴(kuò)散影響到中央測點(diǎn),所以在圖7中可以看到測點(diǎn)溫度很快就有所上升。而對于擴(kuò)散系數(shù)取值較小的情況,擴(kuò)散作用影響有限,溫度只能等到液面上漲到測點(diǎn)位置后憑借對流作用逐漸上升,所以測點(diǎn)溫度在t=420液面抵達(dá)時(shí)才開始上升,而此前保持在初始值左右。不管抵達(dá)的速度快慢如何,測點(diǎn)溫度一旦抵達(dá)其最大值附近便不再繼續(xù)增加,而是在其最大值附近波動,波動原因是自由液面在上升過程中的幾何波動。對于泄水過程,孔口由速度入口邊界修改為壓力出口邊界后,因考慮的是不可壓縮流體,壓力波的傳播速度為無窮大,所以立刻就可以觀察到降壓波。該降壓波疊加在之前液面波動所導(dǎo)致的波動之上。對于擴(kuò)散系數(shù)比較小的情況,降壓波均出現(xiàn)在溫度下降起始點(diǎn)之前,且擴(kuò)散系數(shù)越小,溫度下降起始點(diǎn)出現(xiàn)的時(shí)間越遲,這是因?yàn)闇囟认陆抵饕芸赜跀U(kuò)散作用。在本次模擬中,頂面邊界的溫度是給定的,并且其取值小于注入液體的溫度,所以擴(kuò)散系數(shù)越小,測點(diǎn)溫度受到來自于頂面邊界條件影響所需要的時(shí)間也就越長,而對于擴(kuò)散系數(shù)較大的情形,擴(kuò)散作用明顯,測點(diǎn)溫度下降與降壓波幾乎同時(shí)產(chǎn)生。
圖 7 不同熱擴(kuò)散系數(shù)下區(qū)域中心測點(diǎn)溫度變化Fig. 7 Temperature change curves of central measuring point of calculating region under different thermal diffusion coefficients
利用浸沒邊界法,在格子Boltzmann雙分布函數(shù)模型基礎(chǔ)上添加了能夠考慮自由液面大幅波動的模塊。通過簡單算例驗(yàn)證了該方法的可行性和準(zhǔn)確性,可在一定精度要求下較好地重現(xiàn)自由液面變動過程,為該方法在今后的應(yīng)用奠定了基礎(chǔ)。