熱通量計算方法在開敞海域溫排水模擬中的應用——以田灣核電站為例

鄭曉琴，堵盤軍，徐麗麗，陳昞睿，費岳軍

（國家海洋局東海預報中心，上海 200081）

1 引言

溫排水攜熱量進入海洋后主要有兩個方面的熱量交換[1]，一是高溫的溫排水與環(huán)境水體之間由于熱湍流產(chǎn)生的海水內(nèi)部的水平和垂直方向的熱量交換；二是在海氣界面處與大氣的熱量交換，包括太陽的短波輻射、海面有效回輻射、海水蒸發(fā)或凝結(jié)引起的潛熱通量、海氣溫差引起的感熱通量。由于有這兩方面熱通量的存在，首先溫差所產(chǎn)生的熱浮力效應及熱湍流輸送，使溫排水引起的溫升場具有顯著的三維結(jié)構(gòu)特征，需要建立三維數(shù)學模型才能更加準確地模擬溫排水引起的溫升場[2-4]，同時現(xiàn)有的數(shù)學模型要考慮海氣熱通量才能更加準確地模擬溫排水造成的溫升場[5]。

田灣核電站一期工程（1、2號機組）自2007年投入運行以來，給海洲灣及其附近海域帶來了大量的溫排水，影響了周圍的環(huán)境。對于田灣核電站海域的溫排水研究，大多基于遙感監(jiān)測方面[6]。但是，由于遙感近海水溫調(diào)查的應用過程中受限于遙感測量的深度和計算溫升時溫度本底值的選取[7]。因此，本文應用三維FVCOM數(shù)值模式，采用《物理海洋學》中的海表面熱通量計算方法，建立了田灣核電站附近海域溫排水數(shù)值模型，對溫排水造成的溫升進行了模擬。通過科學、準確定量化預測電廠溫排水對海域溫升范圍、溫度垂直與水平結(jié)構(gòu)，為海灣水體污染、生物資源、環(huán)境容量、生態(tài)服務損失等方面評估提供數(shù)據(jù)基礎，也是溫排水海域管理的重要技術(shù)支撐。

2 研究海域概況

2.1 電廠概況

圖1 江蘇岸外M2分潮等值線圖[8]

圖2 冬季實測大潮表層潮流玫瑰圖

田灣核電站廠址位于江蘇省連云港市后云臺山南側(cè)，距連云港港口約5 km。廠址所在海域?qū)冱S海北部的海洲灣，海域開闊，其南面為黃海灘地，海灘較為平坦，水深變化緩慢，廠址附近水深僅為5 m以下。取水口水域在洋山島南，水域比較開闊，有利于海水置換交流；排水口位于連云港高公島附近，海域水深僅為2 m，附近潮流為往復流，溫排水隨漲落潮流來回擺動。

田灣核電站海水冷卻水系統(tǒng)是采用以黃海為水源和最終熱阱的直流供水系統(tǒng)。海水取水工程采用輸水隧洞和引水明渠的組合取水方案。循環(huán)冷卻水系統(tǒng)的排水暗溝由循環(huán)水泵房出口經(jīng)虹吸井接至雙孔排水暗渠，最終通過東護岸處的循環(huán)冷卻水排水口入海。核電站一期工程兩臺機組的取水量為102.00 m3/s，取排水溫差為10℃，即排水水溫等于取水水溫加10℃。

2.2 潮流特征

田灣核電站所屬的連云港海域受黃海旋轉(zhuǎn)潮波和東海前進潮波的交會作用（見圖1），多年平均潮差為3.69 m，受旋轉(zhuǎn)流和蘇北沿岸流的作用，該海域潮流特征較復雜。根據(jù)2010年冬季在田灣核電站及其附近海域觀測資料的分析得到（見圖2），該海域各測站大潮期間實測流速基本小于100 cm/s，最大流速出現(xiàn)在測流表層的TW06測站，約80 cm/s左右，為旋轉(zhuǎn)流。位于排水口附近的TW32和TW53受到岸線影響，為往復流，流速較小，僅為40 cm/s左右。從表1可以看出，落潮流速普遍比漲潮流速小，該海域為漲潮優(yōu)勢流。

3 溫排水數(shù)值模型

3.1 三維無結(jié)構(gòu)網(wǎng)格水動力數(shù)值模式FVCOM簡介

FVCOM是一套基于有限體積法的三維無結(jié)構(gòu)網(wǎng)格數(shù)值模式[9-10]。它的網(wǎng)格采用來自有限元法的無結(jié)構(gòu)網(wǎng)格，易于擬合岸界和局部加密，而其采用的有限體積法則有離散簡潔、守恒性好等優(yōu)點。使用FVCOM對田灣核電站海域計算，網(wǎng)格可方便地在電廠附近局部加密，使得細致刻畫溫升的擴散范圍成為可能。模式分內(nèi)外模，內(nèi)模時間步長較大，計算三維溫度、鹽度和流速，外模時間步長較小，計算垂向平均二維流速和水位。時間上采用4階-榮格庫塔法。模式整合了Mellor-Yamada的2.5階湍流閉合模型。物質(zhì)輸運平流項（包括溫度、鹽度的平流項）則采用MPDATA方法求解。具體控制方程詳見FVCOM用戶手冊。

表1 研究海域冬季大潮期間各測站實測最大流速單位/（cm/s；/°）

3.2 熱通量計算方法

海表面和水體內(nèi)部存在熱通量交換，所以需要確定海表面熱通量的計算方案，水體內(nèi)部的熱通量可由海表面輻射推得。如果海表面熱通量系數(shù)采用常數(shù)或時間序列插值，那將不能即時反應水溫局部變化對熱通量變化的影響，無法體現(xiàn)海表高溫水體對大氣的加速散熱。FVCOM程序默認的熱通量計算方法是給出常數(shù)或時間序列插值，所以在此替換了原有的熱通量計算方法，新的方法參考了葉樂安和李鳳岐的《物理海洋學》中介紹的熱通量計算方法[11]。

海表面熱通量為：

Qθ=Qs-Qb-Qe-Qh(w/m2)，即單位為焦耳每秒每平方米。其中，Qθ海水吸收的凈熱能,Qs通過海面進入水體的輻射能,Qb有效回輻射,Qe海水蒸發(fā)消耗的熱能（潛熱）,Qh水氣之間的顯熱交換的熱能（感熱）。其中在計算Qs時，晴空時到達海面的總輻射Qs0通過隨時間變化的日地訂正距離和太陽高度角計算得到，考慮了晝夜和季節(jié)變化。

熱通量計算方法中，用到了相對濕度、氣溫等參數(shù)，這些參數(shù)根據(jù)田灣核電站附近海域歷史實測資料設定[12]。廠址所在區(qū)域為南溫帶半濕潤季風氣候區(qū)，年平均氣溫為14.3℃，以8月份最高，平均為26.8℃，以1月份最低，為1.1℃。年平均濕度為69%，以7月份最高，平均為83%，12月份最低，平均為61%。根據(jù)這些資料，模式熱通量系數(shù)中氣溫冬季取2，夏季取27；相對濕度冬季取65%，夏季取85%。

3.3 模式設置及試驗配置

模式垂向分為均勻的6層。網(wǎng)格在田灣電廠海域附近作了局部加密（見圖3，4），最小網(wǎng)格的邊長約20 m。

模式的外模時間步長設為0.1 s，內(nèi)模時間步長設為1 s。外海開邊界給出了16個分潮（M2、S2、N2、K2、K1、O1、P1、Q1、MU2、NU2、T2、L2、2N2、J1、M1、OO1）的調(diào)和常數(shù)以計算開邊界水位。

在電廠排水口位置（見圖5）設置徑流以考慮底層溫排水，同時在取水口設置與排水等量的底層取水。由于連云港海域溫鹽水平變化較小，這里溫鹽初始場取為均值。選取冬季和夏季作為特征季節(jié)進行有、無溫排放兩種情況進行模擬計算，將計算后溫度場直接相減，得到海域電廠溫排水排放后造成的海域溫升。具體參數(shù)設置詳見表2。

圖3 模式的網(wǎng)格圖

圖4 田灣核電站海域的網(wǎng)格

表2 數(shù)值試驗基本設置表

3.4 模型驗證

采用2010年12月21—22日（大潮）期間的海流觀測數(shù)據(jù)進行模式的驗證。圖5分別給出了三個測站的表、底流速流向的計算值和實測值的比對情況，可見模型計算值與實測值趨勢基本一致。此外，利用車牛山測站（2010年12月15—18日和2010年12月21—24日）、連云港測站（2010年12月15—25日）的潮位資料進行水位驗證（見圖6—8）。通過大小潮實測資料的驗證，模型計算值與實測值基本一致。

實測水溫資料采用2011年7月2—3日（大潮）和7月7—8日（小潮）tw32和tw53測站水溫的實測數(shù)據(jù)。驗證結(jié)果顯示：模擬的溫度在潮周期里面均有高低變化，模式計算的溫度與實際溫度趨勢一致，表明熱通量的計算是有效的（見圖9）。

4 溫升模擬結(jié)果

利用驗證后的三維無結(jié)構(gòu)網(wǎng)格FVCOM數(shù)值模式，分別對冬季和夏季進行了模擬，得到典型時刻的溫升分布圖。

圖5 測站的大潮表底層流速流向?qū)Ρ龋▎挝?d）

圖6 水位驗證點位置

4.1 冬季

圖10是計算海域冬季大潮溫升分布圖。可以發(fā)現(xiàn)，漲急時刻，高溫水從排水口分別向南北兩支擴散，南支的溫升較高，4℃溫升等值線以逆時針向北延伸，中心溫升值可達6—8℃，小于實際的排水溫升。由于關(guān)注海域水深較淺，僅為5 m以下，表底層的形態(tài)相似，但表層影響范圍大于底層。溫升分布與流場特征緊密相關(guān)，根據(jù)流場分布圖（見圖11）看出，漲潮期間，漲潮流自東向西涌入本海區(qū)，潮流有沿近岸自東南向西北方向運動的趨勢；在排水口處有一明顯流速較大的出流區(qū)，但排水口周圍區(qū)域流速較小，僅為20—30 cm/s，且流場成逆時針向往北流。這表明從排水口出來的高溫水首先沿岸向南流，其后隨著漲潮流的頂托作用，溫升場向東北方向移動，呈逆時針旋轉(zhuǎn)的擴散場，然而在排水口東南向有一相對溫升低值，可能該區(qū)域是受轉(zhuǎn)流影響。

圖7 車牛山測站水位驗證

圖8 連云港測站水位驗證

圖9 測站大小潮期間表層溫度驗證結(jié)果（點為實測，線為計算）

從落急流場來看，落潮流呈自北向南運動，并逐漸調(diào)整為西北至東南向流動，落潮流速沿岸向外逐漸增大，因此整個溫升區(qū)域在落潮流的帶動下，到達很遠，其中2℃等溫線幾乎與岸平行。在排水口的東北側(cè)仍留有漲潮時殘留的高溫水。

經(jīng)過溫升距離計算和面積統(tǒng)計，大潮漲急時刻的溫升面積最大，1℃溫升面積為54 km2，溫排水擴散距離可達15.5 km，表層面積大于底層，但是擴散距離表底層相差不大。小潮落急時刻的擴散面積最小，1℃溫升面積為7.6 km2，擴散距離為9.6 km。表3、表4給出了大潮平均的溫升距離和面積統(tǒng)計。

表3 冬季大潮平均的溫升距離統(tǒng)計表(單位/km)

表4 冬季大潮平均的溫升面積統(tǒng)計表（單位/km2）

圖10 冬季大潮期間漲落急溫升平面分布

圖11 冬季大潮期間漲落急流場平面分布

圖12 夏季大潮期間漲落急溫升平面分布

圖13 夏季大潮期間漲落急流場平面分布

4.2 夏季

在排水口區(qū)域高溫水分別向南北擴展，漲急時刻（見圖12）南部的溫度等值線以逆時針向東北方向擴展，中心溫升值超過6℃，從溫升距離和面積統(tǒng)計（見表5、表6）來看，表層的擴散范圍明顯大于底層，海洲灣內(nèi)受溫排放影響，水溫維持高值。與冬季形態(tài)相似的是，從排水口出來的溫升場有一個逆時針旋轉(zhuǎn)的擴散；落潮時整個高溫水沿岸向兩邊流動，南邊沿岸的擴展范圍比較遠（見圖13）。

從大潮平均溫升平面圖和面積統(tǒng)計表來看，各溫升值的擴散距離表底層差距不大，但是溫升面積相差較大，表層溫升面積遠大于中底層。其中表層1℃溫升值面積可達36.54 km2，2℃溫升值面積可達15.62 km2，3℃溫升值面積可達7.48 km2，4℃溫升值面積可達3.87 km2。

可以看出，夏季由于本底溫度較高，熱交換均勻。整體的溫升擴散距離和面積比冬季來得大。

表5 夏季大潮平均的溫升距離統(tǒng)計表(單位/km)

表6 夏季大潮平均的溫升面積統(tǒng)計表（單位/km2）

5 結(jié)論

通過采用三維FVCOM數(shù)值模式，考慮海表面熱通量計算方法，本文對田灣核電站一期工程的溫排水造成的溫升在冬季和夏季分別進行了模擬，主要得到以下結(jié)論：

（1）利用海氣熱界面平衡方程，綜合考慮太陽輻射、海面有效回輻射、潛熱通量和感熱通量，在模型中嵌入塊體公式計算感熱和潛熱通量，對在溫排水影響下田灣的溫度場進行了數(shù)值模擬。模型結(jié)果表明：考慮了海氣界面的熱通量的模型結(jié)果能比較真實反映海域溫度；

（2）從漲落潮模擬結(jié)果得出：漲潮期間，從排水口出來的高溫水首先沿岸向南流，其后隨著漲潮流的頂托作用，溫升場向東北方向移動。受旋轉(zhuǎn)流的影響，從排水口排出的高溫水有一個呈逆時針旋轉(zhuǎn)的擴散場。落潮期間，整個高溫水沿岸向南擴展范圍比較遠；

（3）冬、夏季各溫升值的擴散最遠距離表底層差距不大，但是溫升面積相差較大，表層溫升面積遠大于中底。溫升面積最大值出現(xiàn)在大潮的漲急時刻，最小值出現(xiàn)在小潮的落急時刻；

（4）在同等溫排水量的情況下，夏季溫升面積和溫升距離大于冬季，這是由于夏季本底溫度較高，熱交換均勻的原因。

[1]余宙文.海洋溫排水環(huán)境影響評估的若干問題[A].電廠溫排水環(huán)境影響專題研討會論文集[C].北京:國家環(huán)境保護總局環(huán)境工程評估中心,2006.

[2]曹穎,朱軍政.基于FVCOM模式的溫排水三維數(shù)值模擬研究[J].水動力研究與進展,2009,24(4):432-433.

[3]苗慶生,周良明,鄧兆青.象山港電廠溫排水的實測和數(shù)值模擬研究[J].海岸工程,2010,29(4):4-5.

[4]周巧菊.大亞灣海域溫排水三維數(shù)值模擬[J].海洋湖沼通報,2007,4:37-46.

[5]張長寬,姚靜,陶建峰.濱海火-核電廠溫排放數(shù)值模擬研究進展[J].2010,30(3):85-86.

[6]孫戀君,王鳳英,朱曉翔.田灣核電站溫排水環(huán)境影響遙感調(diào)查[J].中國輻射衛(wèi)生,2011,20(3):330-332.

[7]李靜晶,張永興,姜子英,等.濱海核電廠熱污染評價中幾個問題的探討[J].核安全,2012,1:30-34.

[8]任美鍔.江蘇省海岸帶和海涂資源綜合調(diào)查[M].北京:海洋出版社,1986.

[9]Chen C S,Liu H,Beardsley R C.An unstructured grid,finite volume,three-dimensional,primitive equations ocean model:application to coastal ocean and estuaries[J].Journal of Atmospheric and Oceanic Technology,2003,20:159-186.

[10]Chen C S,Zhu J R,Ralph E,et al.Prognostic modeling studies of the Keweenaw Current in Lake Superior,Part I:formation and evolution[J].Journal of Physical Oceanography,2001,31:379-395.

[11]葉安樂,李鳳歧.物理海洋學[M].青島:中國海洋大學出版社,1992.

[12]劉正暉,等.田灣核電站一期工程環(huán)境影響報告書[R].江蘇核電有限公司,2005.