鹽指型雙擴(kuò)散湍流二維直接數(shù)值模擬研究

楊延濤

(北京大學(xué) 工學(xué)院，湍流與復(fù)雜系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室，北京 100871)

0 引 言

當(dāng)流體密度由于溫度等標(biāo)量場的空間不均勻分布出現(xiàn)差異時(shí)，流體在重力場的作用下會發(fā)生浮力驅(qū)動的對流湍流。對流湍流是眾多自然界和工程流動中熱量和物質(zhì)混合的關(guān)鍵物理機(jī)制。在地球海洋環(huán)境中，由于海水密度主要受溫度和鹽度的影響，溫度和鹽度的非均勻分布必然會導(dǎo)致浮力驅(qū)動的宏觀流動發(fā)生。這種由溫度和鹽度共同驅(qū)動的對流現(xiàn)象稱為雙擴(kuò)散對流湍流。特別地，由于溫度的分子擴(kuò)散速率比鹽度快大概兩個量級，雙擴(kuò)散湍流表現(xiàn)出單標(biāo)量場對流湍流所不具有的眾多獨(dú)特動力學(xué)過程和豐富的流動機(jī)理。

溫鹽驅(qū)動雙擴(kuò)散湍流廣泛分布于海洋上層海水中[1-2]。在熱帶和亞熱帶海域，由于日照加熱和表層海水蒸發(fā)，上層海水中特定深度范圍內(nèi)溫度和鹽度隨深度增加而下降，此時(shí)鹽度梯度驅(qū)動流動而溫度梯度致穩(wěn)流動；雙擴(kuò)散湍流處于鹽指流態(tài)[3]。在高緯度極地海域，由于表層海冰低溫和融化，上層海水中存在溫度和鹽度隨深度增加的區(qū)域，此時(shí)溫度梯度驅(qū)動流動而鹽度梯度致穩(wěn)流動；雙擴(kuò)散湍流處于擴(kuò)散流態(tài)[4]。無論哪種流態(tài)，雙擴(kuò)散湍流導(dǎo)致的溫鹽臺階形態(tài)都具有垂直尺度較小而水平相干性尺度極大的特征，并且對熱量和鹽度的垂直輸運(yùn)有重要的影響[3-5]。

對于鹽指型雙擴(kuò)散湍流，對觀測數(shù)據(jù)的分析表明其可能存在于超過30%的海洋水體中[1]，且能夠大大增強(qiáng)當(dāng)?shù)睾Ｓ虻拇怪被旌闲蔥3]。因此理解鹽指型雙擴(kuò)散湍流的流動結(jié)構(gòu)和輸運(yùn)規(guī)律，并提出參數(shù)化描述，具有非常重要的科學(xué)意義和應(yīng)用價(jià)值。自鹽指型雙擴(kuò)散湍流機(jī)制被提出以來[6]，相關(guān)研究工作非常豐富，包括實(shí)驗(yàn)、計(jì)算和理論工作[7-9]。雖然海洋雙擴(kuò)散湍流是溫度-鹽度體系，但是在實(shí)驗(yàn)中，由于實(shí)驗(yàn)條件的限制，也有采用如鹽度-糖度體系[10-11]、溫度-銅離子體系等[12-13]。數(shù)值模擬方面，由于鹽度場極小的分子擴(kuò)散率對計(jì)算網(wǎng)格分辨率提出非常高的要求，一般為節(jié)省計(jì)算資源會采用具有較小施密特?cái)?shù)的濃度場進(jìn)行模擬，而不直接采用與海水相同的物性參數(shù)[14-15]。

基于自主開發(fā)的高效直接數(shù)值模擬程序[16]，本課題組針對鹽指型雙擴(kuò)散湍流已開展了一系列數(shù)值模擬工作，包括真實(shí)海水物性參數(shù)模擬[17-18]，以及不同物性參數(shù)組合體系的模擬[19]。但是，之前的研究只針對固定密度比和少數(shù)不同瑞利數(shù)展開，關(guān)于不同物性參數(shù)組合的系統(tǒng)性研究還沒有進(jìn)行[19]。因此，本論文主要針對四種不同物性參數(shù)組合，通過變化瑞利數(shù)和密度比開展系統(tǒng)性的二維直接數(shù)值模擬工作，并在此基礎(chǔ)上考察流動結(jié)構(gòu)和傳輸特性的變化規(guī)律。

1 控制方程與數(shù)值方法

1.1 控制方程

本文考慮兩平行平板間的流體流動。平板間距為H，且重力方向垂直于平板。上下板分別給定恒定的溫度和鹽度值。本研究采用線性狀態(tài)方程ρ=ρ0(1 ?βT θ+βSs)，其中：ρ為密度，ρ0為參考狀態(tài)密度，θ為偏離參考狀態(tài)溫度，s為偏離參考狀態(tài)鹽度，βT為熱膨脹系數(shù)，βS為鹽度導(dǎo)致密度變化的線性系數(shù)。本文中腳標(biāo)T代表與溫度相關(guān)物理量，而腳標(biāo)S代表與鹽度相關(guān)物理量。則在Oberbeck-Boussinesq假設(shè)下不可壓縮流動動量方程和溫度、鹽度輸運(yùn)方程為：

其中，u為速度矢量，p為壓強(qiáng)，υ為運(yùn)動學(xué)黏性系數(shù)，g為重力加速度，κT為溫度場分子擴(kuò)散率，κS為鹽度場分子擴(kuò)散率。

溫度和鹽度用兩平板間的恒定溫度差ΔT和鹽度差ΔS無量綱化，速度用自由落體速度U= (g βSΔSH)1/2無量綱化，長度用板間距H無量綱化。則流動控制參數(shù)寫為：

本文還同時(shí)使用劉易斯數(shù)，其為普朗特?cái)?shù)和施密特?cái)?shù)的組合：

注意Le代表兩個組分分子擴(kuò)散率的比值，而密度比代表兩個組分差引起的密度差異的比值。對于鹽指型流態(tài)，密度比越大，表明致穩(wěn)流動的溫度梯度相對于驅(qū)動流動的鹽度梯度的強(qiáng)度越大。

1.2 數(shù)值方法

本文的流動模型中，上下兩個平板處速度場滿足無滑移和無穿透邊界條件，溫度場和鹽度場恒定。水平方向所有物理量滿足周期邊界條件。計(jì)算中選取區(qū)域水平寬度遠(yuǎn)大于流動結(jié)構(gòu)的特征水平尺度以保證周期邊界條件的有效性。由于計(jì)算中涉及大施密特?cái)?shù)和大瑞利數(shù)情況，本文采用二維模型，即只考慮一個水平方向。以往研究表明，對于鹽指型雙擴(kuò)散湍流二維和三維流動非常類似[20]。

直接數(shù)值模擬采用自主開發(fā)程序[16]。程序采用分步-投影法處理速度-壓強(qiáng)耦合。空間離散采用二階中心有限差分格式。時(shí)間方向采用二階Runge-Kutta格式，其中擴(kuò)散使用類似Crank-Nicolson格式的半隱式處理方法，非線性對流項(xiàng)采用類似Adams-Bashforth格式的顯式處理方法，浮力項(xiàng)處理方法與非線性項(xiàng)一致。針對大施密特?cái)?shù)情況，程序特別使用了雙分辨率方法。其中鹽度場在加密網(wǎng)格上計(jì)算，其他物理量在基礎(chǔ)網(wǎng)格上計(jì)算?；A(chǔ)網(wǎng)格和加密網(wǎng)格分別保證能夠解析速度場黏性Kolmogorov尺度和鹽度場Batchelor尺度，兩者的比值為Sc1/2。速度場從基礎(chǔ)網(wǎng)格向加密網(wǎng)格插值采用保證局部無散性的Hermite插值方法。利用該方法，能夠有效提高處理高施密特?cái)?shù)標(biāo)量場的模擬效率。注意由于溫度場擴(kuò)散較快，在本文參數(shù)范圍內(nèi)其特征尺度一般大于速度場，因此溫度場采用與速度場一致的網(wǎng)格分辨率。實(shí)際計(jì)算中考慮到并行效率，基礎(chǔ)網(wǎng)格不能太小，因此加密網(wǎng)格的加密系數(shù)要小于Sc1/2，一般單個方向在2～4倍左右。

1.3 模擬參數(shù)范圍

本研究選取四種不同物性參數(shù)組合，對應(yīng)的普朗特?cái)?shù)和施密特?cái)?shù)分別為（Pr,Sc） = （7, 21）、（7, 70）、（7, 700）和（700, 2100），其中（Pr,Sc） = （7, 700）對應(yīng)海洋真實(shí)物性參數(shù)，（700, 2100）對應(yīng)實(shí)驗(yàn)中使用的鹽度-糖度體系[10-11]，（7, 21）為直接數(shù)值模擬中常用的數(shù)值[14]。為系統(tǒng)研究Le的影響，此處額外計(jì)算組合（7,70）。

對于每個Pr和Sc的組合，計(jì)算三個不同密度比，Λ= 0.8、1.2、1.6。對于每個密度比，Ra= 1×106～1×1010，計(jì)算5個不同瑞利數(shù)。總計(jì)60個算例。圖1給出了Ra-Λ相平面上的模擬參數(shù)范圍。

圖1 Ra-Λ相平面上模擬的參數(shù)范圍Fig. 1 Simulated parameters on the Ra-Λ plane

2 主要結(jié)果與分析

本節(jié)首先討論溫鹽傳輸特性隨控制參數(shù)的變化行為以及相應(yīng)的冪律模型，然后討論流動結(jié)構(gòu)及相應(yīng)特征尺度的行為。

2.1 溫鹽傳輸特性

熱量和鹽度的傳輸率可由努塞爾特?cái)?shù)度量。本文定義努塞爾特?cái)?shù)為對流傳輸率與傳導(dǎo)傳輸率的比值；而湍流運(yùn)動強(qiáng)度可由均方根速度定義的雷諾數(shù)度量。這三個無量綱參數(shù)的數(shù)學(xué)定義如下：

其中，w為垂向速度，尖括號代表全場空間平均和時(shí)間平均；urms為速度矢量大小的全場時(shí)空均方根值。注意公式（10、11）所定義Nu只考慮對流運(yùn)動引起的傳輸率，而沒有考慮傳導(dǎo)引起的傳輸。圖2～圖4展示了上述溫鹽努塞爾特?cái)?shù)和雷諾數(shù)隨流動控制參數(shù)的變化行為?？梢钥闯?，對四種不同物性參數(shù)組合，三個主要輸運(yùn)參數(shù)都隨著Ra的增大而增加，且分別滿足類似的冪律規(guī)律。同時(shí)，密度比Λ增加導(dǎo)致努塞爾特?cái)?shù)和雷諾數(shù)減小；這是因?yàn)槊芏缺仍黾訉?yīng)致穩(wěn)溫度梯度的相對強(qiáng)度增加，從而導(dǎo)致流動和傳輸減弱。但是Λ的影響幅度隨著Le的增加而降低。尤其是鹽度傳輸率，當(dāng)Le= 100時(shí)NuS幾乎不隨Λ變化，見圖3（c）。特別地，對于不同物性參數(shù)組合，努塞爾特?cái)?shù)的取值與Le的大小相關(guān)，而與Pr和Sc的取值無關(guān)。如圖2（a、d）和圖3（a、d），其Le都為3，但是Pr和Sc各自相差100倍，溫度和鹽度努塞爾特?cái)?shù)大小非常相似。但是隨著Pr和Sc的增大，雖然Le保持不變，流動雷諾數(shù)顯著降低。

圖2 溫度努塞爾特?cái)?shù)隨流動參數(shù)變化趨勢Fig. 2 Heat Nusselt number versus control parameters

圖3 鹽度努塞爾特?cái)?shù)隨流動參數(shù)變化趨勢Fig. 3 Salt Nusselt number versus control parameters

對圖2～圖4數(shù)據(jù)的進(jìn)一步分析表明，兩個努賽爾特?cái)?shù)基本遵循Ra0.30的規(guī)律增加，而雷諾數(shù)近似按照Ra0.42冪律變化。由此可假設(shè)存在如下冪律標(biāo)度：

圖4 流動雷諾數(shù)隨流動參數(shù)變化趨勢Fig. 4 Reynolds number versus control parameters

其中Λ的冪指數(shù)取值與物性參數(shù)Pr和Sc相關(guān)，其數(shù)值可根據(jù)數(shù)據(jù)由最小二乘擬合給出。具體結(jié)果見表1?？梢钥吹剑瑢τ诤Ｑ笙嚓P(guān)的物性參數(shù)，即Pr= 7和Sc= 700，鹽度雷諾數(shù)對應(yīng)的冪律指數(shù)ξ非常接近于零，即與Λ相關(guān)性很小。但是對于其他物性參數(shù)組合，如鹽-糖體系等具有較小Le的情況，傳輸率受到Λ的顯著影響。圖5中展示了用相應(yīng)密度比冪律修正的響應(yīng)參數(shù)隨瑞利數(shù)的變化規(guī)律，可以看到對每組物性參數(shù)，數(shù)據(jù)均重合在一起，表明公式（13～15）的冪律模型確實(shí)可以描述響應(yīng)的數(shù)據(jù)。

圖5 用密度比冪律修正后的相應(yīng)參數(shù)隨瑞利數(shù)變化趨勢Fig. 5 Nusselt numbers and Reynolds numbers rescaled by the power laws of density ratio versus Rayleigh number

表1 不同Pr和Sc對應(yīng)的Λ冪律指數(shù)Table 1 Values of Λ exponents for different combinations of Pr and Sc

2.2 鹽指結(jié)構(gòu)

本節(jié)討論典型流動結(jié)構(gòu)和相關(guān)特征尺度。對于鹽指型雙擴(kuò)散湍流，主要的典型流動結(jié)構(gòu)為鹽指結(jié)構(gòu)，即沿垂直方向發(fā)展的細(xì)長對流胞結(jié)構(gòu)，其水平寬度可遠(yuǎn)小于豎直高度，如圖6所示。

圖6 典型鹽指結(jié)構(gòu)對應(yīng)的垂向速度場云圖（Pr = 7，Sc = 700，Ra = 1×1010，Λ = 1.2）Fig. 6 Typical finger structures illustrated by the contours of vertical velocity（ Pr = 7, Sc = 700,Ra = 1×1010, and Λ = 1.2）

為提取如圖6所示鹽指結(jié)構(gòu)的水平寬度，我們在流場中間高度為0.5H的區(qū)域內(nèi)計(jì)算垂向速度w的水平自相關(guān)系數(shù)C(δ)，并以C(δ)第一個零點(diǎn)位置對應(yīng)δ的兩倍作為鹽指寬度df。圖7中給出鹽指寬度df隨流動參數(shù)的變化趨勢。如圖所示，對于大部分算例，鹽指寬度均滿足冪律df～Ra?1/4，該冪律指數(shù)與線性穩(wěn)定性分析給出的結(jié)果一致[6]。但是，對于小Le情況，當(dāng)Ra足夠大時(shí)，鹽指寬度明顯偏離該冪律。如圖8所示，此時(shí)流場中主導(dǎo)結(jié)構(gòu)不再是鹽指結(jié)構(gòu)，而是較寬的對流渦結(jié)構(gòu)，因而不再滿足鹽指結(jié)構(gòu)的冪律規(guī)律。注意偏離冪律規(guī)律的算例發(fā)生在Λ= 0.8，即整體密度梯度處于不穩(wěn)定分層狀態(tài)，與經(jīng)典對流流動類似。對于最大的Le= 100，所有算例結(jié)果都滿足該冪律。當(dāng)流場主導(dǎo)結(jié)構(gòu)為鹽指時(shí)，鹽指寬度隨著Λ的增加而減小。

圖7 鹽指寬度隨流動參數(shù)的變化趨勢Fig. 7 Finger width versus control parameters

圖8 典型對流渦結(jié)構(gòu)對應(yīng)的垂向速度場云圖（Pr = 7，Sc = 21，Ra = 1×1010，Λ = 0.8）Fig. 8 Typical finger structures illustrated by the contours of vertical velocity（ Pr = 7, Sc = 21, Ra = 1×1010, and Λ = 0.8）

3 結(jié) 論

綜上所述，本文針對四種不同物性參數(shù)組合的雙擴(kuò)散系統(tǒng)進(jìn)行不同瑞利數(shù)和密度比的二維直接數(shù)值模擬研究。物性參數(shù)組合包括海洋真實(shí)物性參數(shù)，實(shí)驗(yàn)和數(shù)值模擬常用參數(shù)。模擬瑞利數(shù)范圍涵蓋4個數(shù)量級，同時(shí)密度比與海洋典型參數(shù)相當(dāng)。

結(jié)果表明，熱流、鹽流和流動雷諾數(shù)隨瑞利數(shù)的整體變化趨勢不隨物性參數(shù)的改變發(fā)生顯著變化，但其具體數(shù)值受密度比的影響。密度比的影響幅度與施密特?cái)?shù)和普朗特?cái)?shù)的比值，即劉易斯數(shù)有顯著的關(guān)系。對于小劉易斯數(shù)，密度比的增大導(dǎo)致努塞爾特?cái)?shù)和雷諾數(shù)的減小。隨著劉易斯數(shù)的增加，努塞爾特?cái)?shù)和雷諾數(shù)的減小幅度下降。對于海水對應(yīng)的劉易斯數(shù)近似為100的情況，溫度努塞爾特?cái)?shù)和雷諾數(shù)仍然隨密度比的增大出現(xiàn)明顯下降，但是鹽度努塞爾特?cái)?shù)幾乎不隨密度比的增大發(fā)生變化。基于數(shù)據(jù)，本研究給出了溫鹽努塞爾特?cái)?shù)和雷諾數(shù)隨瑞利數(shù)和密度比變化的冪律模型描述。并確定了幾種實(shí)驗(yàn)和計(jì)算中常用的雙擴(kuò)散系統(tǒng)的冪律指數(shù)。

對于本文中的大多數(shù)算例，流動主導(dǎo)結(jié)構(gòu)為沿垂向發(fā)展的細(xì)長鹽指結(jié)構(gòu)，其豎直高度遠(yuǎn)大于水平寬度。對于小劉易斯數(shù)和大瑞利數(shù)情況，主導(dǎo)結(jié)構(gòu)變?yōu)轭愃茖α鳒u結(jié)構(gòu)。當(dāng)主導(dǎo)結(jié)構(gòu)為鹽指時(shí)，鹽指寬度隨瑞利數(shù)的變化與線性穩(wěn)定性分析給出的結(jié)果類似。對于固定瑞利數(shù)，密度比的增加導(dǎo)致鹽指寬度減小。

本文給出的結(jié)果，如物性參數(shù)組合的影響、傳輸特性的冪律模型、主導(dǎo)結(jié)構(gòu)及尺度的變化趨勢，對理解海洋鹽指型界面?zhèn)鬏斕匦杂兄匾膬r(jià)值。同時(shí)，不同物性參數(shù)結(jié)果的對比，對解讀實(shí)驗(yàn)和模擬數(shù)據(jù)，以及將實(shí)驗(yàn)和模擬結(jié)果應(yīng)用于海洋參數(shù)時(shí)，有相應(yīng)的指導(dǎo)意義。