太湖風(fēng)生流及污染物輸移擴(kuò)散數(shù)值模擬

舒葉華，高晨晨

(1.上海市水利工程設(shè)計(jì)研究院有限公司，上海 200061； 2.上海灘涂海岸工程技術(shù)研究中心，上海 200061)

湖流從成因上可分為重力流、密度流和風(fēng)生流。平原地區(qū)的淺水湖泊，重力流和密度流微乎其微，風(fēng)生流是湖流的主要形式[1-3]。風(fēng)生流是由湖面風(fēng)場對湖泊表層水體的剪切、拖曳作用形成的流動，能顯著影響淺水湖泊底泥再懸浮與污染物輸移擴(kuò)散，進(jìn)而對湖泊水質(zhì)產(chǎn)生重要影響[4-8]，因此研究淺水湖泊風(fēng)生流及其驅(qū)動下的污染物輸移擴(kuò)散規(guī)律具有重要的科學(xué)意義。

太湖是我國五大淡水湖之一，水域面積寬廣，但平均水深不到2 m，是典型的大型淺水湖泊。迄今為止，已有很多學(xué)者借助野外觀測或數(shù)值模擬等手段對太湖風(fēng)生流開展了研究，取得了不少有價(jià)值的成果。梁瑞駒等[2]發(fā)現(xiàn)太湖風(fēng)生流沿垂向和水平向均有較大變化，流向沿垂向和水平向也并不一致。逄勇等[9]研究了湖底地形和不規(guī)則湖區(qū)邊界對風(fēng)生流場的影響。胡維平等[10]基于三維數(shù)值模擬揭示了太湖風(fēng)生流存在垂直切變，表層流向與風(fēng)切應(yīng)力方向偏角不大，底層湖流為補(bǔ)償流。周杰源等[3]提出太湖夏季環(huán)流產(chǎn)生的主要?jiǎng)恿σ蜃訛轱L(fēng)，風(fēng)的大小直接影響平面環(huán)流流速的大小，風(fēng)向結(jié)合地形基本決定了環(huán)流的形態(tài)。王建威等[11]基于原位觀測數(shù)據(jù)，分析了太湖流場垂向分布特征，發(fā)現(xiàn)表層水體的流向與風(fēng)向基本一致，其余各層流場與風(fēng)場的對應(yīng)關(guān)系并不明顯，甚至出現(xiàn)完全反向的情況。Li等[12]利用現(xiàn)場觀測的風(fēng)速和湖流資料，研究了風(fēng)場對太湖不同水深處流場的影響。Luo等[13]建立了太湖三維數(shù)值模型，發(fā)現(xiàn)在盛行風(fēng)作用下湖區(qū)存在2個(gè)順時(shí)針環(huán)流和2個(gè)逆時(shí)針環(huán)流。Liu等[14]基于二維水流-污染物耦合數(shù)值模型，研究了風(fēng)場驅(qū)動下的太湖環(huán)流及污染物遷移過程。Liu等[15]認(rèn)為風(fēng)向和風(fēng)速決定了太湖平面環(huán)流的方向、強(qiáng)度和位置，并采用示蹤劑方法研究了定常風(fēng)作用下不同位置的點(diǎn)源污染物泄露對太湖水質(zhì)的影響。

以往對于太湖風(fēng)生流及污染物輸移擴(kuò)散的數(shù)值模擬還存在一些不足，如模型區(qū)域采用矩形網(wǎng)格進(jìn)行離散，無法很好地貼合太湖復(fù)雜的岸邊界(邊界呈現(xiàn)出鋸齒狀)；網(wǎng)格精度較低，很多甚至是千米量級等，這些不足必然會影響到水流結(jié)構(gòu)的計(jì)算精度。本文通過建立高精度的三維數(shù)值模型對太湖風(fēng)生流及污染物輸移擴(kuò)散進(jìn)行精細(xì)化的數(shù)值模擬，以期為太湖水污染治理和水資源調(diào)度提供參考。

1 研究區(qū)域概況

太湖位于長江三角洲南緣(北緯30°55′40″～31°32′58″、東經(jīng)119°52′32″～120°36′10″)，湖泊面積2 427.8 km2，實(shí)際水面面積2 338.1 km2，湖岸線總長405 km。西側(cè)與南側(cè)湖岸線平順、呈圓弧狀，東側(cè)與北側(cè)湖岸線曲折、灣岬相間。太湖是典型的淺水湖泊，平均水深僅1.89 m，最大水深不超過3 m。湖底地形較為平坦，湖中無深槽，也無大片的湖灘，1 m等深線靠近湖岸。環(huán)太湖現(xiàn)有出入湖河道200余條，其中入湖河道主要分布在湖區(qū)西部和西南部，出湖河道主要位于湖區(qū)東部。

本文結(jié)合河道流量特征、位置分布等信息，將環(huán)太湖口門概化為望亭、銅坑閘、胥口、瓜涇口、太浦閘、長兜港、長興港、大港河、城東港、浯溪橋、百瀆口、雅浦港、直湖港、梁溪河等14個(gè)，編號 1～14，具體布置見圖1。

根據(jù)太湖小雷山氣象站(圖1)2007年12月至2008年11月的逐時(shí)風(fēng)資料，太湖湖面各月平均風(fēng)速主要集中在4～6 m/s，年平均風(fēng)速約為5 m/s。由圖2可知，太湖地區(qū)盛行東南風(fēng)(SE)和西北風(fēng)(NW)，頻率分別為12.5%和12.3%，其次為西北偏北風(fēng)(NNW)，頻率為11.7%。西南風(fēng)(SW)、西南偏西風(fēng)(WSW)和西風(fēng)(W)發(fā)生的頻率很低，均不超過3%。

圖1 太湖分區(qū)及環(huán)湖主要口門分布

圖2 太湖湖面風(fēng)向玫瑰圖

2 模型建立與驗(yàn)證

2.1 控制方程

模型采用三維不可壓縮雷諾平均Navier-Stokes方程作為控制方程，服從Boussinesq假定和靜水壓力假定。采用交替方向隱式迭代法(ADI算法)對控制方程進(jìn)行差分離散求解。其連續(xù)性方程和兩個(gè)水平向動量方程為

(1)

(2)

(3)

式中：vx、vy、vz分別為x、y、z方向的速度分量，m/s；vxs、vys分別為x、y方向的點(diǎn)源速度分量，m/s;η為自由水面高程，m；h為全水深，m；t為時(shí)間，s；f為柯氏力參數(shù)，s-1；g為重力加速度，m/s2；ρ為水體密度，kg/m3；ρ0為水體參考密度，kg/m3；sxx、sxy、syx、syy為輻射應(yīng)力張量分量，kg/(m·s2)；νt為垂向紊動黏性系數(shù)，m2/s；pa為自由表面大氣壓，Pa；Fx、Fy分別為x、y方向的水平應(yīng)力；S為源匯項(xiàng)。

2.2 模型區(qū)域和網(wǎng)格劃分

模型計(jì)算范圍為整個(gè)太湖湖區(qū)，具體計(jì)算區(qū)域與網(wǎng)格劃分如圖3所示。由于東側(cè)與北側(cè)湖岸線較為曲折，為使網(wǎng)格更好地貼合岸線，采用三角形的非結(jié)構(gòu)化網(wǎng)格對模型區(qū)域進(jìn)行水平離散，其中西側(cè)與南側(cè)岸線平順處網(wǎng)格相對較大，最大網(wǎng)格邊長 1 200 m；對岬角及島嶼附近區(qū)域進(jìn)行了網(wǎng)格加密處理，最小網(wǎng)格邊長約60 m。平面上共劃分了23 827個(gè)網(wǎng)格單元，節(jié)點(diǎn)總數(shù)為12 448個(gè)?？紤]到太湖為淺水湖泊，為使模型垂向上也具有較高的分辨率，采用Sigma分層法將模型垂向均勻劃分為5層，即每一層的垂向網(wǎng)格邊長都是該位置處總水深的1/5。Sigma分層法最大的優(yōu)點(diǎn)是垂向分層數(shù)不受水深影響，能夠保證在淺水區(qū)也具有較高分辨率，此外湖底和等Sigma面重合，可以真實(shí)地反映湖底地形。

圖3 計(jì)算區(qū)域與網(wǎng)格

2.3 模型參數(shù)設(shè)置

模型綜合考慮重力、風(fēng)力、底床摩擦力和柯氏力的作用，其中表面風(fēng)力τs、底床摩擦力τb和柯氏力f的計(jì)算公式分別為

τs=ρa(bǔ)cd|uw|uw

(4)

(5)

f=2ωsinφ

(6)

由于環(huán)太湖各口門的空間尺度相對于整個(gè)太湖來說很小，因此模型計(jì)算時(shí)，將太湖視為一個(gè)四周封閉的水體，各主要出入湖河道的流量分別按點(diǎn)匯和點(diǎn)源處理。

2.4 模型驗(yàn)證

2.4.1水位驗(yàn)證

采用太湖西山站和吳溇站的實(shí)測水位資料對模型進(jìn)行驗(yàn)證，實(shí)測資料的時(shí)間跨度為1個(gè)月。各水位站位置見圖1，水位驗(yàn)證過程如圖4所示。驗(yàn)證結(jié)果表明，水位計(jì)算值與實(shí)測值基本吻合，模型精度滿足要求。

圖4 計(jì)算水位與實(shí)測水位的對比

2.4.2風(fēng)生流驗(yàn)證

南京地理與湖泊研究所曾基于多年的現(xiàn)場觀測資料，總結(jié)了太湖在北風(fēng)、西風(fēng)、南風(fēng)及東南風(fēng)4種風(fēng)場作用下形成的環(huán)流方向與平面形態(tài)?，F(xiàn)以東南風(fēng)作用下的風(fēng)生流為模擬對象，模型模擬時(shí)風(fēng)速選擇太湖地區(qū)的平均風(fēng)速5 m/s，用上述建立的三維水流模型模擬太湖風(fēng)生流流場，并將模擬得到的垂向平均流場與南京地理與湖泊研究所的觀測成果做對比，結(jié)果如圖5所示。可以看出，模型模擬的垂向平均流場與南京地理與湖泊研究所實(shí)際觀測的流場在流態(tài)上較為一致，表明模型能夠較好地模擬太湖的水流運(yùn)動和流場結(jié)構(gòu)。

(a) 太湖風(fēng)生流觀測成果

3 結(jié)果與分析

3.1 盛行風(fēng)作用下的太湖風(fēng)生流特征

太湖流域夏季盛行東南風(fēng)，冬季盛行西北風(fēng)[16-17]，因此模型中選取東南和西北兩個(gè)盛行風(fēng)向的定常風(fēng)進(jìn)行數(shù)值試驗(yàn)，風(fēng)速值取太湖地區(qū)的平均風(fēng)速5 m/s。圖5(b)和圖6為5 m/s東南風(fēng)作用下的風(fēng)生流流場模擬結(jié)果。需要指出的是，由于風(fēng)向相反，5 m/s西北風(fēng)作用下的垂向平均流場與各分層流場呈現(xiàn)跟東南風(fēng)完全相反的水流形態(tài)和結(jié)構(gòu)特征。

(a) 表層

在5 m/s東南風(fēng)作用下達(dá)到穩(wěn)定狀態(tài)的垂向平均流場中(圖5(b))，湖區(qū)西南部存在以大雷山為中心的大范圍順時(shí)針環(huán)流，該環(huán)流途經(jīng)大雷山東部向東南直至太湖南岸，此后沿太湖西南側(cè)岸線折向西北，至蘭右山后再轉(zhuǎn)至大雷山北面形成閉合順時(shí)針環(huán)流。太湖湖心區(qū)出現(xiàn)一正(順時(shí)針)兩反(逆時(shí)針)3個(gè)環(huán)流，其中順時(shí)針環(huán)流以平臺山為中心循環(huán)流動，兩個(gè)逆時(shí)針環(huán)流分別位于順時(shí)針環(huán)流南側(cè)和東北側(cè)。此外，在東太湖、胥口灣、貢湖灣以及梅梁灣、竺山灣等湖灣區(qū)域出現(xiàn)幾個(gè)大小不等、方向不同的小尺度閉合環(huán)流。上述結(jié)果與王謙謙[18]和胡維平等[19]的研究結(jié)果基本一致。

圖6為5 m/s東南風(fēng)作用下達(dá)到穩(wěn)定狀態(tài)的分層流場?？梢钥闯?，表層由于受到湖面風(fēng)場的直接剪切和拖曳作用，風(fēng)生流流速較大，大致在0.005～0.100 m/s之間，其中西南沿岸以及岬角處流速較大，湖心區(qū)流速相對較小，風(fēng)生流流向基本上與風(fēng)向一致。底層除沿岸淺水區(qū)域以及湖灣內(nèi)有局部環(huán)流存在外，其余區(qū)域流向與表層大致相反，表現(xiàn)為很明顯的補(bǔ)償流。受底部摩擦力影響，底層流速較表層流速小，在0.005～0.060 m/s之間。中層流態(tài)與底層基本類似，但流速比底層略小。這是由于從湖面到湖底依次為風(fēng)應(yīng)力影響控制區(qū)、壓強(qiáng)梯度力影響占優(yōu)區(qū)以及底摩擦力影響控制區(qū)。表層水體處于風(fēng)應(yīng)力影響控制區(qū)，受湖面風(fēng)應(yīng)力直接影響，流速最大；往下流速逐漸減小，至風(fēng)應(yīng)力影響與壓強(qiáng)梯度力影響相平衡時(shí)，流速達(dá)到極小值；再往下過渡到壓強(qiáng)梯度力影響占優(yōu)區(qū)，流速逐漸增大，至壓強(qiáng)梯度力與底摩擦力影響平衡區(qū)出現(xiàn)流速的極大值；底層水體處于底摩擦力影響控制區(qū)，流速隨水深的增加而逐漸減小[19]。

對比圖5(b)和圖6可知，垂向平均流場和分層流場存在顯著差異，因此在用平面二維模型研究淺水湖泊風(fēng)生流時(shí)，應(yīng)充分認(rèn)識到其在揭示流場垂直結(jié)構(gòu)方面的不足。

3.2 太湖風(fēng)生流驅(qū)動下的污染物輸移擴(kuò)散

污染物輸移擴(kuò)散規(guī)律研究是太湖水質(zhì)水環(huán)境研究中的一項(xiàng)重要內(nèi)容。考慮到太湖污染物輸移擴(kuò)散的復(fù)雜性，本文通過在一具有代表意義的入湖口釋放保守示蹤劑，模擬了太湖在5 m/s東南風(fēng)作用下的污染物輸移擴(kuò)散情況，模擬總時(shí)長為90 d。示蹤劑釋放點(diǎn)設(shè)置在“引江濟(jì)太”工程[20-22]的望虞河入湖口處(圖1中的概化口門1)，引水流量設(shè)置為 100 m3/s，示蹤劑釋放速率為5 kg/s，釋放時(shí)間為模擬時(shí)間90 d中的前30 d。

3.2.1污染物質(zhì)量濃度的垂向分布特性

雖然在定常風(fēng)作用下，太湖流場呈現(xiàn)出明顯的分層流動現(xiàn)象，但由于太湖為淺水湖泊，其垂向尺度遠(yuǎn)小于水平尺度，在垂向擴(kuò)散作用下，污染物質(zhì)量濃度沿水深方向近似均勻分布。圖7為5 m/s東南風(fēng)作用30 d后太湖表、中、底層的示蹤劑質(zhì)量濃度分布?？梢钥闯觯瑢τ谔@類淺水湖泊，由于垂向擴(kuò)散比較強(qiáng)烈，污染物質(zhì)量濃度并未呈現(xiàn)明顯的分層現(xiàn)象，與楊益洪[23]對西湖污染物濃度垂向分布特性的研究結(jié)果吻合。

(a) 表層

3.2.2污染物輸移及平面分布特征

前述分析表明，風(fēng)生流驅(qū)動下太湖污染物質(zhì)量濃度在垂向上并未呈現(xiàn)明顯的分層分布，因此以下在對太湖污染物輸移及平面分布特征進(jìn)行分析時(shí)示蹤劑質(zhì)量濃度場均取垂向平均的質(zhì)量濃度場。

5 m/s東南風(fēng)作用下，示蹤劑在不同時(shí)刻的質(zhì)量濃度分布情況如圖8所示?？梢钥闯?，10 d時(shí)，自望虞河入湖的示蹤劑在太湖東北部沿岸流作用下在貢湖灣北岸形成高濃度“污染帶”。此后，由于受到梅梁灣內(nèi)部兩個(gè)小尺度環(huán)流的牽引影響，示蹤劑一部分進(jìn)入梅梁灣，另一部分則直接越過梅梁灣口進(jìn)入竺山灣。至50 d時(shí)，在以平臺山為中心的順時(shí)針環(huán)流帶動下，示蹤劑進(jìn)入到湖心區(qū)，并最終在西南部順時(shí)針大環(huán)流作用下對南太湖水質(zhì)產(chǎn)生影響。從模擬開始后的第30天直至模擬結(jié)束，太湖北部湖灣區(qū)的“污染物”質(zhì)量濃度始終較高，這是由于北部的湖灣區(qū)存在多個(gè)小尺度環(huán)流，給污染物的聚集提供了水動力條件。這也從一定程度上解釋了夏季太湖藍(lán)藻水華主要集中暴發(fā)在梅梁灣、竺山灣等北部湖灣區(qū)的原因。

4 結(jié) 論

a. 太湖風(fēng)生流流態(tài)較為復(fù)雜，且流場存在明顯的垂向分層。穩(wěn)定后的風(fēng)生流場表層流速較大，流向基本上同風(fēng)向一致；底層流速較小，除沿岸淺水區(qū)域以及湖灣內(nèi)有局部環(huán)流存在外，其余區(qū)域流向與表層大致相反，表現(xiàn)出補(bǔ)償流的特征。

b. 西北風(fēng)作用下的太湖風(fēng)生流流場呈現(xiàn)跟東南風(fēng)作用下完全相反的水流形態(tài)和結(jié)構(gòu)特征。此外，風(fēng)場驅(qū)動下的太湖垂向平均流場和分層流場存在顯著差異，在用平面二維模型研究淺水湖泊風(fēng)生流時(shí)，應(yīng)充分認(rèn)識到其在揭示流場垂直結(jié)構(gòu)方面的不足。

c. 風(fēng)場驅(qū)動下的太湖污染物輸移規(guī)律體現(xiàn)出同其風(fēng)生流特征的高度一致性。太湖為大型淺水湖泊，其垂向尺度遠(yuǎn)小于水平尺度，在垂向擴(kuò)散作用下，污染物質(zhì)量濃度沿水深方向近似呈均勻分布，并未呈現(xiàn)出垂向分層特征。

d. 太湖風(fēng)生流流速較小，水動力條件較差，再加上在盛行風(fēng)向的定常風(fēng)作用下，湖區(qū)出現(xiàn)眾多大小不一的環(huán)流結(jié)構(gòu)，尤其是北部的湖灣區(qū)存在多個(gè)閉合環(huán)流，具備積聚污染物的水動力條件。這也是太湖藍(lán)藻水華主要集中暴發(fā)在梅梁灣、竺山灣等北部湖灣區(qū)的一個(gè)重要原因。

(a) 10 d