張 卓,宋志堯,孔 俊,張 鳳

(1.河海大學港口海岸與近海工程學院,江蘇 南京 210098;2.南京師范大學虛擬地理環(huán)境教育部重點試驗室,江蘇南京 210097)

紊流模型發(fā)展于20世紀70年代,當時主要用于氣體動力研究中,近十多年來,包括k-kl,k-ε模型開始在海洋工程中得到廣泛應用.

k-kl模型是在海洋工程中較為常用的雙方程紊流模型,在模擬大氣海洋邊界層及全球環(huán)流方面有較強的優(yōu)越性[1].自從首個k-kl模型Meller-Yamada(MY)模型于20世紀末在河口、近海及海洋領域得到廣泛使用以來[2],k-kl模型開始被許多學者關注并且提出了各自的改進方法,包括對雙方程(主要是kl方程)系數(shù)的調整,壁函數(shù)的改進[3],壓強應變項封閉模型的改進[4]等,這些改進方法都在各自的應用領域內取得了令人滿意的結果,但缺乏同一種情況下進行多因素綜合比較的例子,這就給k-kl模型在實際工程中的應用帶來了不便.因此,本文主要考慮在風生混合流情形下,不同方法對紊流模型模擬結果的影響,為在實際類似情形下,如河口的鹽水入侵,污染物遷移等情況下,合理地選取、調整和改進k-kl模型提供參考和借鑒.

1 k-kl模型

k-kl模型以紊動能k和紊動能通量kl建立對流擴散方程:

式中:νt——垂向紊動系數(shù);u——流速矢量;P,B,ε——紊動能產(chǎn)生項、浮力項和紊動耗散項;Fwall——壁函數(shù);σk,σl,c1,c2,c3——常數(shù),分別取1.96,1.96,0.9,0.5,其中c3的取值和穩(wěn)定函數(shù)以及穩(wěn)態(tài)里查德森數(shù)關聯(lián);t,z——時間和垂向坐標.

底面邊界條件和表面邊界條件一般采用對數(shù)層邊界,即

式中:u*——摩阻流速;cμ 0——穩(wěn)定函數(shù)在對數(shù)層下的常值;κ——卡門常數(shù),取0.4.

由式(1)～(4)可以數(shù)值求解得到k和kl的分布,再通過式(5),(6)便可以得到垂向紊動系數(shù)和垂向擴散系數(shù)的分布,最后可以利用這些系數(shù)來求解N-S方程及溫度擴散方程.

式中 :cμ,c′μ——穩(wěn)定函數(shù) ;l——紊動尺度 ;ν′t——垂向擴散系數(shù).

通過眾多文獻中k-kl模型的比較[5-7],發(fā)現(xiàn)這些模型主要有3個方面不同:(a)穩(wěn)定函數(shù)cμ.從式(5),(6)可以看到,cμ大小直接影響擴散系數(shù)的分布.(b)壁函數(shù)Fwall.作為kl對流擴散方程中源項的一部分,會對kl的分布產(chǎn)生影響.(c)式(2)中的c3取值.在對c1和c2取值大同小異的情況下,對同一種穩(wěn)定函數(shù)c3的取值其實由穩(wěn)態(tài)里查德森數(shù)Rist唯一確定.

2 算例簡介

Kato等[8]曾通過風生密度混合流試驗來研究混合層的變化情況.在試驗中,由常數(shù)風應力作用在穩(wěn)定分層的水流表面,引起混合層自上而下逐漸發(fā)展.水深取得較大,因此可以忽略底面邊界層的影響,看成無限水深的情況.Prince[9]根據(jù)Kato等的試驗數(shù)據(jù),總結了經(jīng)驗公式:

式中:Dm——混合層深度;u*s——表面摩阻流速,在試驗中取10-2m/s;t——施加風應力的時間.

本文的數(shù)學模型根據(jù)Kato等的試驗建立,水流模型采用ELCIRC斜壓模型[10],計算區(qū)域為5000m×100m,深度為50m,水平網(wǎng)格尺度為100m×20m,南北閉邊界,東西采用定水位開邊界,水流可自由出入.垂向網(wǎng)格采用自下而上逐漸加密的方式,總共分50層.底面阻力設0以消除底面邊界層的影響.初始表面溫度12.5℃,底面溫度10℃,相應的;數(shù)模中混合層深度采用文獻[11]的定義,即從水面開始向下所有k＞10-5m2/s2的區(qū)域都看作混合層.其余參數(shù)如下:網(wǎng)格結點96,單元數(shù)75,表面切應力τs=0.1Pa,時間步長120s,參考密度 ρ0=1000kg/m3.

3 結果分析

3.1 穩(wěn)定函數(shù)的影響

海洋數(shù)值模型必須考慮水流的斜壓及分層效應,通常河道中將cμ直接等于cμ 0的方法顯然已不能滿足需要,因此從Meller開始,陸續(xù)提出了MY,KC,LD,CA,CB,CH等6種穩(wěn)定函數(shù),其中CA,CB,CH都是2000年后提出的.這6種穩(wěn)定函數(shù)都是在雷諾應力方程中簡化而來,差異主要源于對壓強應變相關量的不同封閉模型.而Burchard等[11-12]將這6種穩(wěn)定函數(shù)表達成統(tǒng)一的形式,在準平衡態(tài)下,又可表達為理查德森數(shù)Ri的函數(shù):

將6種穩(wěn)定函數(shù)的曲線繪制于圖1.

本文關心的主要是穩(wěn)定分層區(qū)域即Ri＞0.可以看到,隨著Ri的增大,cμ和c′μ減小,表示了分層作用越強或剪切作用越弱,垂向擴散會因此得到削弱的物理機制.穩(wěn)定函數(shù)對垂向紊動擴散系數(shù)的影響體現(xiàn)在2個方面:一是直接體現(xiàn)在計算中,如式(5)和(6);二是通過影響c3的取值從而決定kl的分布.將壁函數(shù)統(tǒng)一取拋物改進型,穩(wěn)態(tài)里查德森數(shù)Rist取0.19.圖2表示6種穩(wěn)定函數(shù)的混合層邊界位置隨時間變化關系.可以看到,MY,KC,LD的結果比較接近,三者之中以KC的結果最接近經(jīng)驗公式(7),但三者得到的混合層擴散速率都顯得快了些,CA,CB,CH的結果比較接近,總體來講,好于MY,KC,LD,這主要是因為CA,CB,CH在壓強應變項中加入了浮力項、各項異性產(chǎn)生項以及渦重分布項,考慮因素更全面的緣故.

3.2 壁函數(shù)的影響

在近壁區(qū)域,水流處于對數(shù)層區(qū),忽略浮力項且P=ε,則式(2)可轉化為

圖1 準平衡狀態(tài)下6種穩(wěn)定函數(shù)和Ri的關系Fig.1 Relationship between 6 kinds of stability functions and Richardson number Ri under quasi-equilibrium condition

其中 ε=c3μ0k3/2l-1,νt可以通過式(5)得到,只是對數(shù)層下cμ=cμ0,并將k近似為常數(shù),代入式(9),整理后得

如果l=κz,那么可以得到 σl的表達式

可以看到要使 σl始終為正,必須使c2Fwall＞c1.

比較常用的壁函數(shù)有3種[3]:

式中卡門常數(shù)k=0.4,db和ds分別是距離底面和表面的長度.

穩(wěn)定函數(shù)采用CH模型,穩(wěn)態(tài)里查德森數(shù)取0.25,得到3種壁函數(shù)條件下的混合層邊界位置隨時間變化曲線(圖3).可以看到,拋物改進型較拋物線型和對稱型的混合層擴散速率小一些,更接近于經(jīng)驗公式,這主要是因為拋物改進型較拋物線型的值在邊界附近相對較小,σl的值較大,kl的擴散速率下降的緣故.

圖2 不同穩(wěn)定函數(shù)混合層邊界位置隨時間變化關系(水面在50m處)Fig.2 Development of mixed layer depth under different stability functions(water level at a depth of 50m)

3.3 穩(wěn)態(tài)里查德森數(shù)Rist的影響

考慮平衡時狀況,即紊動產(chǎn)生等于耗散,將式(1),(2)簡化后,代入耗散項與紊動尺度之間的關系

結合式(5),(6)可推導出c3和里查德森數(shù)Ri存在如下對應關系

準平衡態(tài)下,穩(wěn)定函數(shù)可表達為理查德森數(shù)Ri的函數(shù),c3僅由Ri確定.此時的式(1),(2)都處于平衡態(tài),Ri稱為穩(wěn)態(tài)里查德森數(shù),用Rist表示.試驗表明,Rist的值在0.2～0.25之間.這就意味c3的取值不僅僅要使結果和經(jīng)驗公式(7)吻合,而且要使Rist的取值在合理范圍,否則便失去了模型本身的物理意義[6].在這點上,MY模型顯然有不足(Rist最大只能取到0.1939),圖4表示CH模型在不同Rist下的混合層擴散曲線.從圖4可以看到,增加Rist會增大同時刻的混合層深度,也就是使垂向紊動加強.為分析產(chǎn)生這種情況的原因,可以首先看某個時刻水層中Ri的分布.表層動量顯然最大,然后向下遞減,因此Ri在表層最小,然后向下逐層增大,一直增大到Rist,此時紊動能k的產(chǎn)生和耗散正好達到平衡,k即保持常數(shù)值.再往下,隨著混摻長度的增加,k的耗散大于產(chǎn)生,紊動逐漸減小,直至消失(消失是指風應力產(chǎn)生的垂向紊動小到某個界限下,如10-5,水體本身的紊動依然存在.對于一個較小的Rist,Ri很快就能達到,因此紊動很早就被抑制;而一個相對較大的Rist,必須在更深的水層中才能達到,所以紊動也要在更深層才被抑制.

圖3 不同壁函數(shù)混合層邊界位置隨時間變化關系(水面在50m處)Fig.3 Development of mixed layer depth under different wall functions(water level at a depth of 50m)

圖4 不同 Rist下混合層邊界位置隨時間變化關系(水面在50m處)Fig.4 Development of mixed layer depth under different Rist(water level at a depth of 50m)

4 結 語

a.穩(wěn)定函數(shù)的選擇會影響k-kl模型的計算結果,這種影響不僅體現(xiàn)在計算紊動擴散系數(shù)中,而且還體現(xiàn)在對kl分布的影響上,比較而言CA,CB,CH的結果稍好于MY,KC,LD.

b.壁函數(shù)采用拋物改進型能使結果更接近實際,而對稱型和拋物線型壁函數(shù)都高估了邊界附近kl的擴散速率.

c.對同一種穩(wěn)定函數(shù)而言,如果k-kl又采用了同一種壁函數(shù),那么影響混合層的主要是Rist.當Rist增大時,混合層也會加快發(fā)展,這是由Rist本身的性質決定的.

[1]范聰慧.多因素對海洋上混合層深度影響的數(shù)值模擬[D].北京:中國科學院,2007.

[2]MELLER G L,YAMADA T.Development of a turbulence closure models for geophysical fluid problems[J].Review of Geophysics,1982,20:851-875.

[3]WARNER J C,SHERWOOD C R.Performance of four turbulence closure models implemented using a generic length scale method[J].Ocean Modelling,2005,8:81-113.

[4]CANUTO V M,HOWARD A,CHENG Y,et al.Ocean turbulenceⅠ:one-point closure model:momentum and heat vertical diffusivities[J].J Phys Oceanogr,2001,31:1413-1426.

[5]GALPERIN B,KANTHA H,HASSID S.A equasi-equilibrum turbulent energy model for geophysical flows[J].J AtmosSci,1988,45:55-62.

[6]龔政.長江口三維斜壓流場及鹽度場數(shù)值模擬[D].南京:河海大學,2002.

[7]BURCHARD H,PETERSON O.Models of turbulence in the marine environment:a comparative study of two-equation turbulence models[J].JMar Syst,1999,21:29-53.

[8]KATO H,PHILLIPS O M.On the penetration of a turbulent layer into stratified fluid[J].J Fluid Mech,1969,37:643-655.

[9]PRICE J E.On the scaling of stress-driven entrainment experiment[J].J FluidMech,1979,90:509-529.

[10]ZHANG Y,BAPTISTA A M.A cross-scale model for 3D baroclinic circulation in estuary-plume-shelf systems:Ⅰ:formulation and skill assessment[J].Continental Shelf Research,2005,25:935-972.

[11]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.

[12]CHENG Y,CANUTO V M,HOWARD A M.An improved model for the turbulent PBL[J].Journal of the Atmospheric Sciences,2002,59:1550-1565.