k-kl模型在風(fēng)生混合流中的應(yīng)用

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

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

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

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

1 k-kl模型

k-kl模型以紊動(dòng)能k和紊動(dòng)能通量kl建立對(duì)流擴(kuò)散方程:

式中:νt——垂向紊動(dòng)系數(shù);u——流速矢量;P,B,ε——紊動(dòng)能產(chǎn)生項(xiàng)、浮力項(xiàng)和紊動(dòng)耗散項(xiàng);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ù)關(guān)聯(lián);t,z——時(shí)間和垂向坐標(biāo).

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

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

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

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

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

2 算例簡(jiǎn)介

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

式中:Dm——混合層深度;u*s——表面摩阻流速,在試驗(yàn)中取10-2m/s;t——施加風(fēng)應(yīng)力的時(shí)間.

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

3 結(jié)果分析

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

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

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

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

3.2 壁函數(shù)的影響

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

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

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

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

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

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

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

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

圖2 不同穩(wěn)定函數(shù)混合層邊界位置隨時(shí)間變化關(guān)系(水面在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的影響

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

結(jié)合式(5),(6)可推導(dǎo)出c3和里查德森數(shù)Ri存在如下對(duì)應(yīng)關(guān)系

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

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

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

4 結(jié) 語

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

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

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

[1]范聰慧.多因素對(duì)海洋上混合層深度影響的數(shù)值模擬[D].北京:中國(guó)科學(xué)院,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]龔政.長(zhǎng)江口三維斜壓流場(chǎng)及鹽度場(chǎng)數(shù)值模擬[D].南京:河海大學(xué),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.