近海區域二階紊流封閉模型的比較研究

張卓，宋志堯，孔俊

（1.河海大學 水文水資源與水利工程科學國家重點實驗室 江蘇 南京 210098；2.虛擬地理環境重點實驗室（教育部） 南京師范大學，江蘇 南京 210097）

近海區域二階紊流封閉模型的比較研究

張卓1，宋志堯2，孔俊1

（1.河海大學 水文水資源與水利工程科學國家重點實驗室 江蘇 南京 210098；2.虛擬地理環境重點實驗室（教育部） 南京師范大學，江蘇 南京 210097）

將雷諾應力方程不同封閉方式歸結為穩定函數的差異，并以此作為區分標準，總結了6種國內外采用過的穩定函數并以統一的形式列出，分析不同紊流模型穩定函數的分布和影響因素，并且通過實驗資料進行了比較分析，發現（1）2000年后提出的3種穩定函數得到的臨界理查德森數大于2000年前的3種，更符合實際情況，這與其在壓強應變相關項和浮力相關項中增加了旋轉效應和重分布效應有關；（2）盡管這3種模型考慮的因素更加全面，但在某些特定場合下，如邊界層，與實驗的吻合程度反不如前3種較早的穩定函數，反映了封閉模型本身存在的缺陷。最后，提出了二階紊流封閉模型還待解決的問題。

二階紊流封閉模型；垂向紊動擴散；臨界理查德森數；

引 言

紊流模型發展于上世紀70年代，最初主要用于氣體動力研究中，而真正能直接運用到近海水流的紊流模型卻不多，Launder (1975) 和Rodi (1976) 對單方程，雙方程紊流模型有系統的闡述。直到1982年海洋模式POM (Princeton Ocean Model) 的提出，里面附帶了一個能處理垂向紊動擴散的模型，這就是最早的二階紊流封閉模型——MY (Mellor, Yamada) 模型。從此，二階紊流封閉模型開始逐步取代0階模型和經驗公式，成為三維模型中計算垂向紊動擴散的主流算法。近些年來，隨著人們對紊流機理的認識加深及計算機性能的增強，二階紊流封閉模型又有了新的發展。由于它在模擬邊界混合層中的良好表現，在海洋數值計算中得到了廣泛的應用，成為解決海洋中溫度，鹽度垂向擴散問題最主要的方式之一[1]。

近年來，隨著計算機性能的提高，大渦模擬 (LES) 和直接數值模擬 (DNS) 已開始應用于紊流的研究，但這些精細程度遠高于紊流模型的方法在工程實際問題中卻應用寥寥，僅限于處理一些邊界簡單的問題。這主要是因為：（1）在工程應用中，人們更關心的是紊流對平均流場及其輸運效果帶來的影響，而對于紊流細節一般并不關心。（2）從目前看，二階紊流封閉模型達到了精度和計算效率的完美統一，在滿足工程要求精度的條件下，計算量遠小于大渦模擬和直接模擬。（3）紊流模型在上世紀90年代后又有了新的發展[2]，考慮的問題更加全面，使模擬結果更接近實際，并可以統一的形式表達，就使人們可以將新的模型和原有的模型進行比較，程序編制容易，人們也樂于實行。因此，二階紊流封閉模型仍是最具實用價值的方法。

本文將雷諾應力方程的不同封閉方式歸結為穩定函數的差異，并以此作為區分標準，總結了6種國內外采用過的穩定函數并以統一的形式列出，比較在不同影響因素下的穩定函數取值差異并與前人的實驗結果和經驗公式進行對比。通過這些比較，可以為如何選用和改進二階封閉模型提供參考，最后，總結全文并對二階紊流封閉模型的問題及解決方案提出自己的看法。

1 二階紊流封閉模型

1.1 雷諾應力方程

式中：Pij為剪切產生項，Gij為浮力產生項，Φij為壓強應變項，εij為耗散項，具體形式如下：

式中：表示流速時均值；ui和p表示脈動流速和壓強；ρ0和ν分別是參考密度和運動粘滯系數；b可

表達為?gρ'/ρ0;xi表示坐標；δij是克羅內克符號；Dij表示雷諾應力擴散項。將雷諾應力輸運方程里的

j取i, 可以得到常見的紊動能方程，上式中＜.＞表示時均。浮力通量輸運方程式為：

由式（1），式（6）和式（11）構成整個二階封閉模型的控制方程?？梢钥吹剑死字Z應力與浮力通量以外，還有一些二階甚至三階量需要采用近似或模型進行封閉，這些量包括：壓強應變項Φij和壓強

下面將分別介紹如何封閉和簡化這些項。

1.2 壓強應變相關項Φij和壓強浮力相關項的封閉

Φij和對雷諾應力的分布有著重要的影響，對它的建模是整個二階封閉模型的關鍵，一般二階封閉模型的分類也基于此。紊流模擬工作者的任務就是通過各種假設和近似，將它們簡化成為工程應用中可以接受的形式。

從1975年Launder開始，陸續有學者提出各自的模型來封閉Φij和根據文獻[2]，這些模型可以

寫成通用的形式：

表示無量綱各項異性張量；第二項到第四項表示應變（包括變形和旋轉）引起能量重分配，代表中的非線形效應，一般都被忽略。 最后一項表示浮力的貢獻。式 (13) 類似，只是第四項表示浮力貢獻，最后一項表示浮力自相關量的貢獻。限于篇幅，式 (12) 和式 (13) 的具體展開形式不再贅述，詳見參考文獻[2]。

Mellor, Yamada (MY) 和Kantha (KC) 模型忽略了浮力及旋轉對雷諾應力的壓強重分配效應，而本文列舉的其他模型對這些作用均予以不同程度的考慮。這些模型都可以表示為式 (12) 和式 (13) 的形式，只是每項的系數c1～ c5, cb1～ cb5不同。

1.3 kb,εij 和的求解

(11) 得到，但一般用一個簡化的方法，即假定式 (11) 達到平衡態，等號左端為0，可以得到下式關系：r被看作常數，各個模型的取值有所不同，但差別也不大，在0.47 ～ 0.82之間。

值得一提的是，假定紊動的耗散都通過小渦來完成，而小渦一般認為是各向同性的，因此認為紊動的耗散項為各向同性即

一般假定浮力脈動和流速脈動沒有必然的相關性，所以

1.4 顯式代數應力模型（EASM）

由于雷諾應力方程過于復雜，Rodi等首先提出了代數應力模型的設想[5]。主要思路是設法將應力的的微分方程簡化為代數表達式，同時保持紊流各項異性的特點。采用Rodi第一近似，假設應力＜uiuj＞的對流擴散線性正比于紊動能k的對流與擴散，即：

將式 (18)、式 (19) 代入 (1) 和 (6) 便可得到代數應力方程的形式，擴散項Dij和的 作用相當于用Dk來按比例近似?；蛘哒f，將和的作用統一進了紊動能擴散項中，此時雷諾應力及通量方程中的未知量只剩紊動能k和耗散項ε，可在雙方程模型中離散求解。根據式 (14)，式 (19) 等號右邊第二項等于0。

因為這里的紊流模型主要是用于近海區域，屬于“淺水型”，水平尺度遠大于垂向尺度，兩者比值一般達到或超過104。通過尺度分析，可以忽略雷諾應力及浮力通量的水平對流和擴散，在剪切產生項中可忽略對水平向的導數，只留對垂向的導數，這樣便大大簡化了方程。

為便于表達，首先定義下面關系：

從形式上, 式 (20) 很像忽略掉對水平向的導數后的Boussnesq假設，但又有區別，Boussnesq假設中

的νt和是各項同性的變量，而在代數應力模型中是二階張量。將式 (20) 代入代數應力方程，可以將νt

和表示成為：式中，n0, n1, n2, nb0, nb1, nb2, d0, d1, d2, d3, d4, d5都是常數，其值依賴于在壓力應變項中采用的不同模型。表1給出了不同模型的取值。

和被稱為穩定函數，是和的顯函數，前者表示剪切效應，是引起不穩定的因素，后者表

示浮力效應，是引起穩定（分層）的因素。

從式（21）可以看到，垂向渦粘系數的分布取決于兩個要素：一是穩定函數，主要由對壓強應變項建立的模型決定；二是紊動參數k和ε，主要由雙方程模型決定。而雙方程模型的分歧主要在第二個方程中及其采用不同的物理量，包括k-kl模型、k-ε模型和k-ω模型等。本文的研究主要集中于第一個方面，即不同模型的穩定函數分布。

表 1 六種穩定函數的系數Tab.1 Coefficients for the six stability functions

2 穩定函數的比較分析

在工程上較常用的穩定函數主要是準平衡態形式，這里的平衡指的是局部平衡，即紊動能的產生和耗散正好抵消時的狀態。本節對6種穩定函數在這種情形下的分布進行比較。這6種穩定函數的系數分布列于表1中。

2.1 準平衡態下的穩定函數

準平衡態是指局部區域的紊動能產生等于耗散，寫成式子就是：

式中：P=0.5Pii,G=0.5Gii分別代表由剪切和浮力引起的紊動產生，忽略掉水平梯度項，將式 ( 20 ) 代入式 (24)，可寫成：

通過式 (23) 和式 (21)，可以將式 (24) 寫成如下形式：

2.2 六種模型下的穩定函數分布比較

根據式 (23)αM的定義，舍去式 (29) 中沒有物理意義的解，可以得到αM～Ri關系曲線（圖1）。從圖1可以看到，隨著Ri的增加，αM呈增加的趨勢，并且Ri有個極限值Ric，當Ri→Ric時，αM→+∞。這表明Ri的取值是有界限的。不同的模型這個界限值不一樣，表2給出六種穩定函數Ric的取值。結合圖1可以看到，早期的模型Ric都比較小，在0.3以下，而2000年以后的模型Ric值都接近1。

從圖2可以看到，在準平衡態下，穩定函數cμ和cμ'都是Ri的函數。各模型分布趨勢一致即隨Ri數增大而減小, 但由于不同模型對壓強應變項的不同假設，因此在數量上存在明顯差異。這些差異主要表現在：(1) 2000年前提出的模型 (MY, KC, LD) 穩定函數變化較陡，而2000年后的模型 (CA, CB, CB) 穩定函數變化相對較緩；（2）Ri的取值范圍不同，MY, KC, LD最大可取值到0.19 ～ 0.28，CA, CB, CH最大可取到0.84 ～ 0.97。這是由臨界理查得森數Ric決定的（表2）。在實際水流中，Ric的值也并不唯一，如：Miles等曾于1961利用線性平衡理論，推導出Ric=0.25；Abarbanel等于1984年用非線性理論導出Ric=1；近年來的大渦模擬和直接模擬都表明紊流的Ri等于甚至大于1都存在[6]。因此，CA, CB, CH的結果更符合實際些。

圖 1 不同模型的αM ～ Ri的關系曲線Fig.1 Curve of αM ～ Ri for different stability functions

圖 2 準平衡狀態下穩定函數和Ri的變化關系Fig.2 Quasi-equilibrium versions of stability functions displayed as functions of Ri

2.3 理論值和實驗值的比較

很少有實驗資料是直接關于穩定函數分布的，因此，需要通過間接的方法以及在一些特殊情況下的取值來進行驗證和比較。

首先探討對數層中cμ的取值。通常在定常明渠中水流的流速呈對數型分布，而在非均勻水流中一般也認為在近底區域的流速呈對數分布。在對數層，紊動完全通過層與層之間的剪切作用產生，并且達到局部平衡即P=ε，將式 (21)、式 (25) 代入可得到：

把式 (30) 入式 (31) 就可以求得，一般令此時的將六種模型的值列于表2。

將紊流模型計算得到的結果與經驗公式以及Webster和Gerz的實驗數據[3]做了對比，如圖3所示?？梢钥吹剑S著Ri數的增加，σt逐漸增加。MY, KC, CB的結果和經驗公式（34）極其接近，只不過MY和 KC受Ric所限，僅僅能表示Ri小于0.25的部分。而CA的結果和實驗數據的相關性則更好一些。

最后，將穩定函數用Monin-Obukhov[7]的形式表示

圖 3 Prandtl數與Ri的關系曲線 (包括6種紊流模型，經驗公式（34）和Webster以及Gers的實驗結果[7])Fig.3 Curve between Prandtl number and Ri ( including 6 turbulence models, equation (34) and measurements by Webster and Gers

圖 4 Monin-Obukhov形式下的理論解和實驗數據比較Fig.4 Comparison of theory and experiment in form of Monin-Obukhov

將6種模型與Businger的實驗值（圖4）的比較發現，MY, KC, LD除了φM～?曲線在非穩定的部分（?＜0）偏小以外，其余部分與實驗值吻合較好。說明MY模型在壓強應變項中忽略掉浮力和旋轉的影響，至少在Monin-Obukhov相似理論適用的邊界層區域是符合實際情況的。至于考慮了浮力和旋轉作用的CA, CB, CH與Businger的實驗值的吻合程度為何反不如MY模型，這可能與這些模型的系數率定以及模型本身的缺陷有關。CA,CB,CH壓力重分布的假設很可能還有不合理的地方，至少在邊界層區域，式 (12)，(13) 需要進行修正。而對于準平衡態的穩定函數，當Ri ＞ Ric時，Cμ→0,相當于此時沒有紊動擴散存在，這又與實際情況不一致，因為紊動能可以通過周圍水體將紊動傳播過來，因此可以考慮使用非平衡態穩定函數。

3 結 論

本文總結了二階封閉模型從雷諾應力方程簡化為顯式代數應力模型的整個過程，發現影響垂向渦粘系數和擴散系數分布的要素有兩個：一是穩定函數，主要由對壓強應變項建立的模型決定；二是紊動參數k和ε，取決于雙方程模型。在本文采用6種模型對穩定函數的分布情況進行的比較研究中發現，2000年后提出的模型CA, CB, CH的臨界理查德森Ric比較早的MY, KC, LD的大一些，說明前者允許的紊流范圍比后者大一些，這是符合實際情況的。而采用Monin-Obukhov的邊界層相似理論的結果與實驗數據比較時卻發現， MY, KC, LD與實驗數據吻合更好些，反映了CA, CB, CH模型的參數率定以及在雷諾應力方程簡化中本身存在的一些問題,需進一步在實際應用中檢驗。

[1]Mellor G L, Yamada T.Development of a turbulence closure model for geophysical fluid problems [J].Rev Geophys Space Phys, 1982, 20 (4): 851-875.

[2]Umlauf L, Burchard H.Second-order turbulence closure model for geophysical boundary layers.A review of recent work [J].Continental shelf Research 2005, 25: 795-827.

[3]Canuto V M, Howard A, Cheng Y, et al.Ocean turbulence.Part Ⅰ: one-point closure model.Momentum and heat vertical diffusivities [J].Journal of Physical Oceanography, 2001, 31: 1 413-1 426.

[4]Cheng Y, Canuto V M.Howard A M.An improved model for the turbulent PBL [J].Journal of the Atmospheric Science, 2002, 59: 1 550-1 565.

[5]Rodi W.A new algebraic relation for calculating the Reynolds stresses [J].Zeitschrift fur Angewandte Mathematic und Mechanik, 1976, 56: T219-T221.

[6]Wang D, Large W G, McWillimas J C.Large eddy simulation of the equatorial ocean boundary layer: Diurnal cycle, eddy viscosity and horizontal rotation [J].J Geophys Res, 1996, 101 (C2): 3 649-3 662.

[7]Monin A S, Yaglom A M.Statistical Fluid Mechanics: Mechanics of Turbulence [M].1971: 1-769.

Comparison of two-order closure model in coastal areas

ZHANG Zhuo1, SONG Zhi-yao2, KONG Jun1

(1.State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai University, Nanjin, 210098, China; 2.Key Laboratory of Virtual Geographic Environment (Ministry of Education), Nanjing Normal University, Nanjing 210097, China)

With advance of knowledge on turbulence and computer capability in recent years, two-order closure model has made new process.Different two-order closure models lead to different stability functions.Three earlier models (before 2000) and three later ones (after 2000) have been concluded in a unified form and then have been compared in the paper. After that some experimental data are compared with the model results.It is found that: (1) the later models (after 2000) obtain more realistic Richardson Number than the earlier models (before 2000) ; (2) in some particular tests, such as boundary layer, performance of earlier models is even superior to the later’s, which reflects the faults of the later ones. Finally, problems to be solved are put forward.

two-order closure model; vertical turbulent mixing; critical Richardson Number

P731.21; P735

A

1001-6932 (2010) 01-0012-10

2009-01-06；

2009-06-18

水文水資源與水利工程科學國家重點實驗室開放基金（2007491211）

張卓（1981－），男，浙江海寧人，博士研究生，主要從事港口，海岸及近海工程研究。電子郵箱：mercury1214@hhu.edu.cn