楊延濤

(北京大學 工學院，湍流與復雜系統(tǒng)國家重點實驗室，北京 100871)

0 引 言

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

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

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

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

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

1.1 控制方程

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

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

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

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

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

1.2 數(shù)值方法

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

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

1.3 模擬參數(shù)范圍

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

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

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

2 主要結果與分析

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

2.1 溫鹽傳輸特性

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

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

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

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

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

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

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

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

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

2.2 鹽指結構

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

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

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

圖8 典型對流渦結構對應的垂向速度場云圖（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 結 論

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

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

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

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