σ坐標系下淺水方程水平擴散項守恒型式

于守兵

(1.南京水利科學研究院,江蘇南京 210024;2.水文水資源與水利工程學科國家重點實驗室,江蘇南京

210024)

σ坐標系下淺水方程水平擴散項守恒型式

于守兵1,2

(1.南京水利科學研究院,江蘇南京 210024;2.水文水資源與水利工程學科國家重點實驗室,江蘇南京

210024)

三維淺水模型控制方程經σ坐標變換后,其水平擴散項還存在一變換因子,故不能被常見的擴散模型表示成守恒型式。對σ坐標系下淺水方程水平擴散項進行處理,得到守恒型式的擴散模型,并采用垂向二維靜水鹽度擴散算例進行驗證,結果表明守恒型模型能夠保持初始的鹽度分布。

σ坐標系;淺水方程;水平擴散項;鹽度擴散

三維水流模型中遇到的一個突出問題是復雜的地形變化。Philips在1957年提出自動邊界擬合的垂向σ坐標變換,后被引入到海洋模型并得到很快發展。然而,水平擴散項經σ坐標變換后其表達式變得十分復雜,給數值離散帶來不便,并帶來很大的離散誤差。

很多學者對此進行了研究,并提出處理擴散項的全模型、Mellor模型[1]和Huang模型[2-3],以及不同的插值方法[4-5]等。這些模型都是針對擴散項本身的,而淺水方程擴散項在σ坐標變換后,還有一變換因子,采用上述模型得到的擴散項不能表示成守恒型式。故在有限體積(FVM)離散時作為源項處理,不能保持擴散通量在整體區域上的守恒性。因此,有必要尋找守恒型的擴散模型。

1 σ坐標系中水平擴散項模型

1.1 淺水方程中的水平擴散項

引入算子:

則直角坐標系下各微分在σ坐標系下表示為

直角坐標系下某一變量Φ在x和y方向的擴散項為

式中:DΦ為擴散系數,對于水平流速Du=Dv=vt。

引入σ坐標系下垂直方向的速度協變量ω:

則得到三維淺水模型控制方程為

為進一步說明問題,將σ坐標系下淺水方程中水平擴散項統一表示為如下型式

1.2 常見的水平擴散項模型

目前的擴散模型主要是針對型如式(2)的擴散項,由于x和y方向的處理方式相同,下面以x方向為例進行說明。

根據式(1),可以得到x方向擴散項的表達式為

該式直接由σ坐標變換得到,未作任何處理,稱為全模型。

Mellor等[1]注意到全模型計算結果與實際存在很大差異,對式進行簡化處理,提出Mellor模型:

這種表達式簡單,計算易實現。但是,Mellor模型處理σ變換中擴散項的二階混合偏導項時省略的因素較多,在模擬諸如無初始運動,溫度僅為垂直方向的函數且水平擴散系數為常數的流動時,會出現溫度界面的熱傳導這一虛假擴散。

Huang等[2]對σ變換變換產生的二階混合偏導項進行簡化處理,提出Huang模型:

該模型保留了擴散項中重要的部分,略去一些次要項,很好地防止了由于σ變換而導致的虛假擴散的發生。但是,在地形坡度較大時,計算誤差仍較大。不能滿足Hanney[6]提出的“靜水一致性”條件。

Huang等[3]提出一種方法,水平梯度仍在原直角坐標系中處理,而離散變量由σ坐標系下變量插值確定。該方法體現了問題的本質,計算精度較高。張景新等[5]在插值方法上對其進行改進,并提高該方法的使用范圍。

然而,σ坐標系下淺水方程中水平擴散項(3)還存在一變換因子φσ,也即水深H。那么,考慮φσ時采用上述常見的三種擴散模型對應的擴散項分別為(仍以x方向為例):

可以看出,這三種擴散項的表達式中每一項都不是關于變量Φ的某種關系式的散度,也即不是守恒型式。由于在FVM離散時通常將不能表示成散度的項歸入源項,這樣采用上述擴散模型得到σ坐標系下淺水方程中的水平擴散項只能作為源項處理,從而不能在理論上保證擴散通量在整個計算區域上的守恒。

2 淺水方程中水平擴散項的守恒型式

為得到守恒型的擴散模型,對淺水方程中x方向水平擴散項進行變換:

根據引入的算子:

那么,可得到x方向擴散項表達式

將其中直角坐標系下偏導數用σ坐標系表示,則式(5)變為

同理,可得到y方向擴散項表達式:

那么,由式(5)和式(6)可得到淺水方程中水平擴散項(3)的表達式:

式(7)與上述三種擴散模型表達式最大不同之處在于其中每一項都表示成變量Φ的某種關系式的散度,也即守恒形式。這樣,在FVM離散時就能夠以擴散通量處理,從而在理論上保證擴散通量在整個計算域上的守恒。從表達式組成來看,守恒模型采用了兩次直角坐標下偏導數坐標系下的形式,簡潔程度僅次于最簡單的Mellor模型,因此具有較高的計算效率。

3 數值試驗

3.1 靜水鹽度擴散算例

用于評價σ坐標系下擴散模型的算例是二維垂直平面上具有傾斜底床的封閉靜止水域中的鹽度擴散問題。計算域水平方向長6 km,最大水深260m,最淺處80m,底床坡度為0.03(圖1)。初始鹽度分布以指數型式給定:

式中:Hmax為最大水深,Hmax=260;cmax=34 ppt為計算域鹽度最大含量。

圖1 靜水鹽度擴散算例求解域與兩個固定觀測點垂向位置Fig.1 Computation area and vertical locationsof two fixed observatory pointsof stillwater salinity diffusion test

直角坐標系下垂向二維靜水鹽度擴散方程為

式中:Dx為鹽度的水平擴散系數,Dx=5m2/s。

由初始鹽度分布(8)可知水平方向的鹽度梯度等于零,故由式(9)計算得到的鹽度分布將不隨時間而變。如果在σ坐標系下計算得到的濃度場與初始濃度場有所不同,則一定是由所采用的擴散模型引起的。故通過對該濃度擴散問題的計算,可以檢驗σ坐標系下擴散模型的可行性。

3.2 模型計算結果分析

計算時采用三維模型模擬本算例,根據式(7)得到σ坐標系下的靜水濃度擴散方程為:

與垂向二維方程不同的是,三維模型中還包含y方向的水平擴散項,擴散系數取與x方向相同,也即Dx=Dy=5m2/s。

設定計算域寬1 km,采用開源程序EasyMesh以邊長為200m的三角形進行網格剖分,共得到364個單元(圖2)。除若干單元外,大部分單元均勻性良好。垂向均勻分為7層(圖1)。表層、底層以及固壁處邊界條件為擴散通量為零。

圖2 求解域平面三角形剖分及兩個固定觀測點平面位置Fig.2 Plan triangle dissection of computation area and locationsof two fixed observatory pointsof stillwater salinity diffusion test

圖3 鹽度等值線分布Fig.3 Salinity contour

圖3 (a)給出初始鹽度等值線分布,圖3(b)給出采用守恒型模型計算30天后得到的等值線分布。可以看出,守恒型模型得到的鹽度等值線分布與初始分布基本吻合,尤其是垂向位于中部的25 ppt等值線分布變化非常小,其它等值線只是在模型的控制邊界處有所差異。對于固定的兩個觀測點A和B(垂向位置見圖1,水平位置見圖2),鹽度隨時間的變化曲線如圖4所示,可以看出兩個點的鹽度變化很小,A點濃度由30.42 ppt變為30.43 ppt(圖4(a)),B點濃度由23.09 ppt變為23.04 ppt(圖4(b))。

圖4 固定點鹽度變化Fig.4 Salinity change at fixed observatory points

4 結 語

由于淺水方程σ坐標變換后水平擴散項還存在一變換因子,目前的擴散模型不能將其表示成守恒型式,從而不能保證整個計算域上擴散通量的守恒。推導了淺水方程水平擴散項的守恒型式,并采用垂向二維鹽度擴散算例進行驗證,結果表明該擴散模型能夠保持鹽度的初始分布狀態。

[1] Mellor GL,Blumberg A F.Modeling vertical and horizontal diffusivitieswith the sigma coordinate system[J].MonthlyWeather Review,1985,113(8):1379-1383.

[2] Huang W R,Spaulding M.Modeling horizontal diffusion with sigma coordinate system[J].Journal of Hydraulic Engineering,1996,122(6):349-352.

[3] Huang W R,Spaulding M.Reducing horizontal diffusion errors inσ-coordinate coastaloceanmodelswith a second order lagrangian-interpolation finite-difference scheme[J].Ocean Engineering,2003,29:495-512.

[4] Stelling GS,Kester AM.On the approximation of horizontal gradients in sigma coordinates for bathymetrywith steep bottom slopes[J].International Journal for Numerical Methods in Fluids,1994,18:915-935.

[5] 張景新,劉 樺.σ坐標系下水平擴散項的有限差分計算[J].力學季刊,2006,27(3):377-386.

[6] Haney R L.On the pressure gradient force over steep topography in sigma coordinate oceanmodels[J].J.Phys.Oceanography,1991,21:610-619.

Conservativemodel of horizontal diffusion term of shallow water equations inσcoordinate

YU Shou-bing1,2
(1.Nanjing Hydraulic Research Institute,Nanjing 210024,China;2.State Key Laboratoryof Hydrology-Water Resourcesand Hydraulic Engineering,Nanjing 210024,China)

For shallow water equations inσcoordinate,the horizontal diffusion term has an factor,which leads to non-conservative form exploiting now diffusionmodels.The paper dealswith the diffusion term and deducesa conservative diffusionmodel,which is tested by the classical example of salinity diffusion in stillwater.The result shows that the conservative diffusionmodel is capableof preserving the initial salinity status.

coordinate;shallow water equations;horizontal diffusion term;salinity diffusion

TV131.2

A

1005-9865(2010)03-0123-05

2009-11-30

南京水利科學研究院博士學位論文基金資助項目(Yy20803)

于守兵(1980-),男,河南襄城人,博士,從事河口海岸及近海工程方面的研究。E-mail:yushoubing@163.com