基于浸沒(méi)邊界法的水庫(kù)變動(dòng)水面模擬及驗(yàn)證

周 睿，程永光，吳家陽(yáng)

(武漢大學(xué) 水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430072)

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

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

1 自由液面處理算法介紹

1.1 控制方程

首先引入階梯函數(shù)即Heaviside函數(shù)H，規(guī)定在一種流體(比如液體)中H為1，在氣體中H為0，自由液面則標(biāo)記H為0.5。為將Heaviside函數(shù)與Dirac函數(shù)聯(lián)系起來(lái)，將它表示為一維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)的變量為第二類變量。首先考慮梯度是針對(duì)第一類變量計(jì)算，而積分是針對(duì)第二類變量計(jì)算，這樣梯度符號(hào)可以和積分號(hào)交換順序。同時(shí)Dirac對(duì)于兩類變量是反互易的，所以梯度符號(hào)也可以針對(duì)第二類變量計(jì)算，即

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

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

在此▽?duì)?ρ1-ρ0。

流體的運(yùn)動(dòng)假定由N-S方程描述，在考慮表面張力的情況下，有

上述方程適用于全流場(chǎng)(含密度場(chǎng))。β為所涉及問(wèn)題的維數(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ù)的。

1.2 N-S方程求解

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

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

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

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

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

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

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

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

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

表面張力計(jì)算參考拉普拉斯公式：

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

式中：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ù)即可得到張力的面密度。采用平面三角形來(lái)離散曲面，即網(wǎng)格只能通過(guò)加密才能逐漸逼近真實(shí)的曲面；而在實(shí)際中無(wú)法還原曲面，也可采用曲面三角形來(lái)離散液面，按照以下方法來(lái)計(jì)算曲面上某一點(diǎn)的平均法向主曲率：

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

在自由液面的模擬中，液面與計(jì)算區(qū)域的邊界相交，從上面計(jì)算某一拉格朗日點(diǎn)的表面張力來(lái)看，該點(diǎn)需要完全被其他點(diǎn)所包圍，而這一條件對(duì)于處在邊界上的拉格朗日點(diǎn)來(lái)說(shuō)是無(wú)法滿足的?？梢栽谶吔缤饷嬖傩略鲆慌殴?jié)點(diǎn)，這排節(jié)點(diǎn)并不會(huì)像內(nèi)部節(jié)點(diǎn)一樣參與拉氏力的計(jì)算及其擴(kuò)散和位置的更新，其坐標(biāo)根據(jù)具體所設(shè)置的邊界條件人為地給定，比如對(duì)于滑移固壁、對(duì)稱邊界和周期邊界來(lái)說(shuō)，液面在邊界處的法向量應(yīng)該與邊界平行。

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

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

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

2 水箱底部注水液面變動(dòng)算例

2.1 計(jì)算條件說(shuō)明

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

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

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

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

為分析液面的波動(dòng)特征，將測(cè)點(diǎn)布置在液面正中央，并統(tǒng)計(jì)了該點(diǎn)縱坐標(biāo)變化情況，結(jié)果如圖4所示。在注水前期t＜100時(shí)，液面中央上升速度非?？欤@是突起形成階段。之后曲線在100＜t＜200時(shí)趨于水平，這是由于液面突起達(dá)到最高點(diǎn)后向下回彈所致。這段曲線體現(xiàn)該波動(dòng)疊加在液面總體上漲趨勢(shì)，液面上升趨勢(shì)和突起向下的波動(dòng)疊加在一起相互抵消，所以曲線成水平狀。波動(dòng)過(guò)程一直持續(xù)，波動(dòng)幅值同時(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 液面中央測(cè)點(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ù)性定理，液面以下液體體積的增長(zhǎng)速度應(yīng)該等于入口流量。同理在放水階段液面下降的速度也應(yīng)該等于出口的放水流量。從圖5中可以看出，雖然在注水階段液體體積的增長(zhǎng)曲線并非嚴(yán)格成一條直線，但與理論直線的差距處于可接受范圍。在放水階段由于給定的是壓力邊界，出口速度并非均勻分布，所以理論曲線也并非是一條直線，得到的結(jié)果與理論解在總體趨勢(shì)一致的基礎(chǔ)上存在一定偏差。究其上述偏差的來(lái)源，可能一方面來(lái)自自由液面的離散處理，由于采用的是三角形平面面元來(lái)離散液面，實(shí)際上的液面應(yīng)該是光滑曲面；另一方面來(lái)自自由液面與計(jì)算區(qū)域邊界相交面的處理。前文已經(jīng)闡述了如何處理液面與物理邊界相交的問(wèn)題，但實(shí)際上卻無(wú)法完全做到。

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

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

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

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

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

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

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

3 結(jié) 語(yǔ)

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