張卓，宋志堯，孔俊

（1.河海大學(xué) 水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室 江蘇 南京 210098；2.虛擬地理環(huán)境重點(diǎn)實(shí)驗(yàn)室（教育部） 南京師范大學(xué)，江蘇 南京 210097）

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

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

（1.河海大學(xué) 水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室 江蘇 南京 210098；2.虛擬地理環(huán)境重點(diǎn)實(shí)驗(yàn)室（教育部） 南京師范大學(xué)，江蘇 南京 210097）

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

二階紊流封閉模型；垂向紊動(dòng)擴(kuò)散；臨界理查德森數(shù)；

引 言

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

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

本文將雷諾應(yīng)力方程的不同封閉方式歸結(jié)為穩(wěn)定函數(shù)的差異，并以此作為區(qū)分標(biāo)準(zhǔn)，總結(jié)了6種國內(nèi)外采用過的穩(wěn)定函數(shù)并以統(tǒng)一的形式列出，比較在不同影響因素下的穩(wěn)定函數(shù)取值差異并與前人的實(shí)驗(yàn)結(jié)果和經(jīng)驗(yàn)公式進(jìn)行對比。通過這些比較，可以為如何選用和改進(jìn)二階封閉模型提供參考，最后，總結(jié)全文并對二階紊流封閉模型的問題及解決方案提出自己的看法。

1 二階紊流封閉模型

1.1 雷諾應(yīng)力方程

式中：Pij為剪切產(chǎn)生項(xiàng)，Gij為浮力產(chǎn)生項(xiàng)，Φij為壓強(qiáng)應(yīng)變項(xiàng)，εij為耗散項(xiàng)，具體形式如下：

式中：表示流速時(shí)均值；ui和p表示脈動(dòng)流速和壓強(qiáng)；ρ0和ν分別是參考密度和運(yùn)動(dòng)粘滯系數(shù)；b可

表達(dá)為?gρ'/ρ0;xi表示坐標(biāo)；δij是克羅內(nèi)克符號；Dij表示雷諾應(yīng)力擴(kuò)散項(xiàng)。將雷諾應(yīng)力輸運(yùn)方程里的

j取i, 可以得到常見的紊動(dòng)能方程，上式中＜.＞表示時(shí)均。浮力通量輸運(yùn)方程式為：

由式（1），式（6）和式（11）構(gòu)成整個(gè)二階封閉模型的控制方程。可以看到，除了雷諾應(yīng)力與浮力通量以外，還有一些二階甚至三階量需要采用近似或模型進(jìn)行封閉，這些量包括：壓強(qiáng)應(yīng)變項(xiàng)Φij和壓強(qiáng)

下面將分別介紹如何封閉和簡化這些項(xiàng)。

1.2 壓強(qiáng)應(yīng)變相關(guān)項(xiàng)Φij和壓強(qiáng)浮力相關(guān)項(xiàng)的封閉

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

從1975年Launder開始，陸續(xù)有學(xué)者提出各自的模型來封閉Φij和根據(jù)文獻(xiàn)[2]，這些模型可以

寫成通用的形式：

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

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

1.3 kb,εij 和的求解

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

值得一提的是，假定紊動(dòng)的耗散都通過小渦來完成，而小渦一般認(rèn)為是各向同性的，因此認(rèn)為紊動(dòng)的耗散項(xiàng)為各向同性即

一般假定浮力脈動(dòng)和流速脈動(dòng)沒有必然的相關(guān)性，所以

1.4 顯式代數(shù)應(yīng)力模型（EASM）

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

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

因?yàn)檫@里的紊流模型主要是用于近海區(qū)域，屬于“淺水型”，水平尺度遠(yuǎn)大于垂向尺度，兩者比值一般達(dá)到或超過104。通過尺度分析，可以忽略雷諾應(yīng)力及浮力通量的水平對流和擴(kuò)散，在剪切產(chǎn)生項(xiàng)中可忽略對水平向的導(dǎo)數(shù)，只留對垂向的導(dǎo)數(shù)，這樣便大大簡化了方程。

為便于表達(dá)，首先定義下面關(guān)系：

從形式上, 式 (20) 很像忽略掉對水平向的導(dǎo)數(shù)后的Boussnesq假設(shè)，但又有區(qū)別，Boussnesq假設(shè)中

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

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

和被稱為穩(wěn)定函數(shù)，是和的顯函數(shù)，前者表示剪切效應(yīng)，是引起不穩(wěn)定的因素，后者表

示浮力效應(yīng)，是引起穩(wěn)定（分層）的因素。

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

表 1 六種穩(wěn)定函數(shù)的系數(shù)Tab.1 Coefficients for the six stability functions

2 穩(wěn)定函數(shù)的比較分析

在工程上較常用的穩(wěn)定函數(shù)主要是準(zhǔn)平衡態(tài)形式，這里的平衡指的是局部平衡，即紊動(dòng)能的產(chǎn)生和耗散正好抵消時(shí)的狀態(tài)。本節(jié)對6種穩(wěn)定函數(shù)在這種情形下的分布進(jìn)行比較。這6種穩(wěn)定函數(shù)的系數(shù)分布列于表1中。

2.1 準(zhǔn)平衡態(tài)下的穩(wěn)定函數(shù)

準(zhǔn)平衡態(tài)是指局部區(qū)域的紊動(dòng)能產(chǎn)生等于耗散，寫成式子就是：

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

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

2.2 六種模型下的穩(wěn)定函數(shù)分布比較

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

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

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

圖 2 準(zhǔn)平衡狀態(tài)下穩(wěn)定函數(shù)和Ri的變化關(guān)系Fig.2 Quasi-equilibrium versions of stability functions displayed as functions of Ri

2.3 理論值和實(shí)驗(yàn)值的比較

很少有實(shí)驗(yàn)資料是直接關(guān)于穩(wěn)定函數(shù)分布的，因此，需要通過間接的方法以及在一些特殊情況下的取值來進(jìn)行驗(yàn)證和比較。

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

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

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

最后，將穩(wěn)定函數(shù)用Monin-Obukhov[7]的形式表示

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

圖 4 Monin-Obukhov形式下的理論解和實(shí)驗(yàn)數(shù)據(jù)比較Fig.4 Comparison of theory and experiment in form of Monin-Obukhov

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

3 結(jié) 論

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

[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

水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室開放基金（2007491211）

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