基于OpenFOAM的異重流數(shù)值模擬分析

丁偉業(yè)，金 生，林金波

(大連理工大學 海岸和近海工程國家重點實驗室，大連 116024)

異重流是指密度相差不大的兩種可以相混的流體在適宜條件下因密度差異而產(chǎn)生的相對運動。在運動過程中，各種流體基本能夠保持原有的面貌，不因交界面上存在的紊動摻混作用而發(fā)生全局性的混合現(xiàn)象。1892年Forel[1]在萊曼湖中首次觀察到了攜沙水流在注入湖水中時并未與湖水混合，而是潛入到湖水底部形成異重流的現(xiàn)象。作為河流注入?yún)R水體三種方式之一的異重流，1953年由Bates[2]在研究三角洲時給出定義。Mulder和Syvitski[3]于1995年重新修訂了Bates關(guān)于異重流的定義。他們將異重流定義為一種負浮力的流體。按照產(chǎn)生的原因，異重流可以分為溫差異重流(如熱電廠冷卻池)、鹽水異重流(如河口鹽水楔)以及渾水異重流(如挖入式港池，水庫等)。對于河口地區(qū)，上游淡水泄入海洋中，與隨潮上溯的咸水混合形成異重流。異重流在其流動的過程中會攜帶大量的泥沙，隨著泥沙濃度的增加，泥沙顆粒構(gòu)成絮凝團,引起沉降，最終形成淤積。因為港池、船閘、渠道中的流體常常是靜止的，具有形成異重流的條件，容易因為異重流的入侵形成淤積進而阻礙航行。因此，關(guān)于異重流的研究對于水道港口的正常運行，乃至河口海岸演變有著重要的實際意義。

Keulegan[4]采用實驗的方式，獲得了一系列的關(guān)于異重流鋒頭速度特性的曲線。Benjanmin[5]指出，當上、下層流體達到渠道半深時，上、下層流體分別以不同的速度前進。在研究異重流流動的過程中Benjanmin發(fā)現(xiàn)，能量損失的最大值約為10%，能量損失較少。Shin[6]在完善Benjamin的理論的基礎(chǔ)上指出，當Re較大時，異重流的能量耗散可以忽略不計。Hartel[7-9]采用直接數(shù)值模擬(DNS)，更細致準確的觀察了異重流鋒頭的形狀與特性。Birman[10]在研究異重流鋒頭空間結(jié)構(gòu)的過程中觀察到了Kelvin-Helmholtz (KH)不穩(wěn)定，與Hartel所觀察到的一樣，KH涌浪都是成對出現(xiàn)的。Cantero[11]通過二維和三維數(shù)值模擬對異重流的鋒頭結(jié)構(gòu)與速度等相關(guān)特性進行了研究。他指出通常高Re情況下，二維異重流的模擬中界面處翻滾的渦更為強烈。

為了更清楚的了解異重流的結(jié)構(gòu)以及相關(guān)特性，本文利用開源計算軟件OpenFOAM對異重流的進行數(shù)值模擬。利用有限體積法對控制方程進行離散，PISO算法對方程進行求解。每一個PISO算法的內(nèi)迭代步中的壓力梯度通過壓力泊松方程來求解[12]。

1 數(shù)值模型

本文研究的異重流是由不同密度的流體相對運動引起的。相應的數(shù)值模型的基本微分方程包括連續(xù)方程，雷諾時均方程以及輸運方程。

▽·u=0

(1)

(2)

其中：t為時間；u=(u,v)為笛卡爾坐標系下的速度矢量場；ν為動力粘度；ρ為流體密度；p*=p-ρg·x為修正壓力；g為重力加速度；x=(x,y)為笛卡爾坐標。

流體密度的變化是一個對流-擴散的過程。相應的方程為

(3)

式中：Dm為分子耗散系數(shù)；Sct=νt/Dk為紊流施密特數(shù)；νt為渦粘度；Dk為渦耗散。

2 模型驗證

如圖1所示，在平底渠道中兩種不同密度的流體被閘門分隔。當閘門被抽取后，密度大的流體向下流動，密度小的流體向上流動，此過程被稱作異重流。許多研究人員花費大量的精力通過實驗[13]或者數(shù)值模擬[9, 14-16]來研究異重流發(fā)生過程中的細節(jié)。作為經(jīng)典算例，異重流特別適合用來驗證非靜壓海洋模型。在異重流發(fā)生的過程中，包含了切應力驅(qū)動的流體混合和內(nèi)波。流體粘性的大小以及模型的精度均會對異重流產(chǎn)生明顯的影響。

圖1 計算域示意圖Fig.1 Schematic illustration of lock-exchange

表1為5T時刻不同計算域離散尺度的異重流特征值收斂性分析。根據(jù)Berntsen[14]描述，xf=(x-Lx/2)/H為ρ0+△ρ/2處最小的x值。大密度流體運動最前端的位置，Umax,min為水平速度的最大值和最小值，Wmax,min為垂向速度的最大值和最小值。在收斂性分析算例中粘度為1E-6 m2/s。由表中結(jié)果可知，模擬結(jié)果的穩(wěn)定性與網(wǎng)格精度相關(guān)。隨著網(wǎng)格精度的增加，采用MS和SS兩種離散方法得到的計算結(jié)果之間的誤差逐漸減小。當網(wǎng)格精度達到1 600×200時，兩種離散方法的計算結(jié)果最為相近。在后文的模擬過程中，均采用1 600×200網(wǎng)格精度。

表1 收斂性分析， γ= 1E-6 m2/sTab.1 Convergence analysis with 1E-6 m2/s

表2 v= 1E-6 m2/s時本文的數(shù)值模擬結(jié)果與Berntsen的結(jié)果對比Tab.2 Comparisons in the simulations with v=1E-6 m2/s between the present results and the results from Berntsen

表2為5T時刻，粘度為1E-6 m2/s情況下，本文采用的模型與Berntsen[14]模擬結(jié)果的對比。從表2中可知，MS與SS獲得的Umin，Umax和Wmin存在明顯的差值。這表明動量離散格式對模擬結(jié)果有著重要的影響。MS模擬的結(jié)果與BOM模型的結(jié)果相近，SS模擬的結(jié)果則與MITgcm模型給出的結(jié)果相近。MS與BOM各特質(zhì)值之間的誤差分別為0%，1.95%，2.98%，0.72%和1.86%。SS與MITgcm各特質(zhì)值之間的誤差分別為2.7%，1.68%，0.37%，0.96%和1.16%。

圖2 數(shù)值模擬異重流時空演化Fig.2 Temporal and spatial variations for the lock-exchange calculated using numerical model

圖2為v= 1E-6 m2/s情況下，采用MS方法模擬的異重流流場演化過程。從圖2中可以觀察到，當在異重流傳播過程中，其兩種流體的頭部均形成James[17]所述的膨脹波。當異重流演化一段時間后，輕流體前緣的膨脹波仍然穩(wěn)定，然而在膨脹波影響下重流體前緣產(chǎn)生了耗散。如圖2-b所示，在膨脹波后方，有內(nèi)涌浪形成。從圖2-a～2-c中可以看出，界面處的兩種流體從層流過渡成為紊流。兩種流體以不同的速度傳播。流體界面因剪切速度變得不穩(wěn)定，進而不斷生成渦。對于由不同密度的連續(xù)流體形成的異重流而言，不同速度的流體在界面處的不穩(wěn)定現(xiàn)象可以用Kelvin-Helmholtz (KH)不穩(wěn)定來描述。KH不穩(wěn)定狀態(tài)下生成的渦稱為KH渦。從圖2-b中可以觀察到，在鋒頭后方的第一個KH渦更為明顯。隨著KH渦逐漸摻混到重流體，兩種流體界面處形成混合渦，如圖2-d所示。在前進的異重流鋒頭后方，有新的KH渦生成。此時，重流體鋒頭的高度近似為H。

3 異重流特征值分析

圖3和圖4分別為5T和10T時刻采用v= 5.6E-8 m2/s, 1E-6 m2/s, 2E-6 m2/s和5E-6 m2/s四種粘度的密度分布。圖中變換后的無量綱坐標為x1=(x-Lx/2)/H和z1=(z-H)/H。由圖3展示的結(jié)果可見，SS的耗散性小于MS。這表明動量方程的對流格式對密度的分布有重要的影響。當粘性為5.6E-8 m2/s和1E-6 m2/s的時候，MS與SS的結(jié)果具有明顯的差別。粘度為5.6E-8 m2/s時，在SS模擬的結(jié)果中可以捕捉到更多小尺度的渦，這與BOM模型獲得的結(jié)果相一致。然而，當v=1E-6 m2/s時，可以觀察到相對小尺度的中等渦，這表明MS得到的結(jié)果與MITgcm模型獲得的結(jié)果幾乎完全一致。在Ai[16]和Hartel[9]都能觀察到類似的結(jié)果。由此可知，本文所采用的模型能夠精確的模擬異重流。當使用粘度2E-6 m2/s和5E-6 m2/s時，MS和SS在5T時刻得到的密度分布之間的差別已經(jīng)非常小了，并且均與BOM模型得到的結(jié)果相近。

注:左側(cè)欄為MS方法的結(jié)果,右側(cè)欄為SS方法的結(jié)果。從上到下的粘度為5.6E-8m2/s,1E-6m2/s,2E-6m2/s和5E-6m2/s。圖3 5T時刻密度分布Fig.3Resultsofdensitydistributionsat5T

相比于5T時刻，10T時刻異重流發(fā)展的更為充分。通過圖4可以觀察到，當v= 2E-6 m2/s和5E-6 m2/s時，流場內(nèi)的渦關(guān)于點(0,0) 近似對稱分布。隨著粘度的減小，密度場的渦分布逐漸變得復雜。MS得到的密度場的耗散性強于SS得到的密度場。在采用SS方法計算的密度場中能夠觀察到更多小尺度的復雜渦。

注:左側(cè)欄為MS方法的結(jié)果,右側(cè)欄為SS方法的結(jié)果。從上到下的粘度為5.6E-8m2/s,1E-6m2/s,2E-6m2/s和5E-6m2/s。圖4 10T時刻密度分布Fig.4Resultsofdensitydistributionsat10T

表3 前緣速度Uf，F(xiàn)roude數(shù)Fr和Reynolds數(shù)Re本文的數(shù)值模擬結(jié)果與Berntsen的結(jié)果對比Tab.3 Comparisons of front speed Uf, Froude number Fr and Reynolds number Re between the present results and the results from Berntsen

4 結(jié)論

本文利用開源CFD軟件OpenFOAM對異重流進行了數(shù)值模擬，并對其演化過程中的結(jié)構(gòu)變化以及特征值進行了分析。異重流在傳播過程中，兩種流體前緣均會形成膨脹波。相對于上層的輕流體，膨脹波會對重流體前緣的耗散產(chǎn)生更明顯的影響。在膨脹波后方會有內(nèi)涌浪生成。由于上、下兩層流體以不同的速度傳播，流體界面會產(chǎn)生KH不穩(wěn)定現(xiàn)象，進而形成KH渦。隨著異重流的演化，KH渦不斷與重流體混摻，進而生成混合渦。當異重流穩(wěn)定傳播時，其前緣高度為水槽高度的一半。通過對不同粘性的異重流模擬發(fā)現(xiàn)，流體粘性越大，異重流演化過程中的耗散現(xiàn)象越明顯。在數(shù)值模擬過程中，動量方程的離散格式對模擬的結(jié)果有著顯著的影響。當流體粘性較小時，采用MS方法模擬的異重流耗散現(xiàn)象比SS方法得到的結(jié)果更為明顯。通過對異重流特征值的分析可知，MS方法與BOM模型得到的結(jié)果更為接近，SS方法與MITgcm模型得到的結(jié)果擬合的更好。