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

(大連理工大學(xué) 海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連 116024)

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

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

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

1 數(shù)值模型

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

▽·u=0

(1)

(2)

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

流體密度的變化是一個(gè)對(duì)流-擴(kuò)散的過程。相應(yīng)的方程為

(3)

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

2 模型驗(yàn)證

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

圖1 計(jì)算域示意圖Fig.1 Schematic illustration of lock-exchange

表1為5T時(shí)刻不同計(jì)算域離散尺度的異重流特征值收斂性分析。根據(jù)Berntsen[14]描述，xf=(x-Lx/2)/H為ρ0+△ρ/2處最小的x值。大密度流體運(yùn)動(dòng)最前端的位置，Umax,min為水平速度的最大值和最小值，Wmax,min為垂向速度的最大值和最小值。在收斂性分析算例中粘度為1E-6 m2/s。由表中結(jié)果可知，模擬結(jié)果的穩(wěn)定性與網(wǎng)格精度相關(guān)。隨著網(wǎng)格精度的增加，采用MS和SS兩種離散方法得到的計(jì)算結(jié)果之間的誤差逐漸減小。當(dāng)網(wǎng)格精度達(dá)到1 600×200時(shí)，兩種離散方法的計(jì)算結(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í)本文的數(shù)值模擬結(jié)果與Berntsen的結(jié)果對(duì)比Tab.2 Comparisons in the simulations with v=1E-6 m2/s between the present results and the results from Berntsen

表2為5T時(shí)刻，粘度為1E-6 m2/s情況下，本文采用的模型與Berntsen[14]模擬結(jié)果的對(duì)比。從表2中可知，MS與SS獲得的Umin，Umax和Wmin存在明顯的差值。這表明動(dòng)量離散格式對(duì)模擬結(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ù)值模擬異重流時(shí)空演化Fig.2 Temporal and spatial variations for the lock-exchange calculated using numerical model

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

3 異重流特征值分析

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

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

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

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