海灣電廠三維斜壓水流和溫排水數值模擬

姚 靜，張長寬，陶建峰

（河海大學港口海岸與近海工程學院水文水資源與水利工程科學國家重點實驗室，南京 210098）

隨著經濟建設的快速發展，海灣電廠的規模也逐步擴大。電廠所需的循環冷卻水直接排入海水中，改變了水體的物理、化學及生態環境，可能造成熱污染，因此有必要建立相應的模型對其影響進行預測。溫水與受納水體的摻混及傳熱過程本身是一個三維的過程，加之海灣地帶地形和水動力條件的復雜性，要了解溫排水垂向水流特性和溫升分布，必須進行三維模擬。已有學者將三維模型運用于溫排水的研究中［1-2］。在對流場及溫度場進行模擬時，大多都是先模擬水動力場，然后在此基礎上模擬溫度場，且忽略由于密度差異引起的密度梯度和斜壓效應，這一效應在水動力條件復雜的海灣地帶以及垂向密度分層明顯的溫水排放口附近水域尤為明顯［3］。因此很有必要建立三維斜壓水流和溫排水數學模型。

本文以象山灣內某電廠的溫排水為例，進行三維斜壓水流和溫排水數值模擬。先進行水動力場和溫度場計算，將溫度場計算結果代入狀態方程，實時修正水體密度，再將密度反饋到斜壓項，進行流場更新計算，進而得到新的溫度場，依此循環，實現水動力場和溫度場的耦合計算。在根據實測資料對模型進行充分驗證的基礎上，模擬了電廠溫排水擴散輸移過程中的三維溫升分布特征。

1 三維斜壓水流和溫排水數學模型

1.1 控制方程

模型水平方向采用正交曲線坐標系，垂直方向采用Phillips［4］于1957年提出的σ伸縮坐標變換，可得σ坐標下的控制方程。

（1）水動力方程

式中：ζ為潮位；D為全水深；u、v和ω分別為σ坐標系下ξ、η和σ方向的流速分量；h1和h2分別為坐標轉換系數；f為Coriolis力系數；g為重力加速度；ρ、ρ0分別為海水密度和海水參考密度；pa為海水自由表面大氣壓，正常天氣下作為常數處理；AM為動量水平渦粘系數；KM為動量垂向紊動渦粘系數；式（2）和式（3）等號右邊第三項即為斜壓項。

（2）溫排水方程

式中：T為溫升；AH為溫度水平擴散系數；KH為反映溫度垂向紊動混合的垂向擴散系數；R為源匯項，即取排水流量；Ks為海面綜合散熱系數；Cρ為水體的比熱。

（3）狀態方程

由于海水密度是由壓力、鹽度和溫度組成的一個非常復雜的非線性函數，密度根據式（5）［5］確定

式中：S為鹽度；p為壓力；ρ（S，T，0）為一個標準大氣壓（海壓為0）下的海水密度狀態方程；c由式（6）計算。

1.2 定解條件

初始條件以零啟動的形式給出，即全場3個流速分量給定為零；同時給出初始水位和初始溫升場，其值根據電廠附近水域實測資料插值給出。

邊界條件分為閉邊界條件和開邊界條件。流場岸邊界采用不可入邊界條件，底部邊界以壁函數法修正邊界鄰近計算節點上的變量值；溫度場閉邊界則給定為絕熱條件。潮位開邊界根據潮位站實測潮位過程外推給定；流場和溫度開邊界條件采用輻射條件確定。自由面熱力學邊界條件僅考慮水面散熱，忽略其他熱通量等邊界條件。取排水口邊界根據計算工況流量與熱量條件給定。

1.3 主要參數

模型水平渦粘系數AM和擴散系數AH采用Smagorinsky公式［6］計算。垂向渦粘系數KM和擴散系數KH采用Mellor-Yamada 2.5階紊流閉合模型［7］計算。根據全國通用公式［8］，考慮3.5 m/s的海面風速時，夏季平均水文氣象情況下的海面綜合散熱系數Ks=47 w/（m2℃）。

1.4 計算方法

針對限制時間步長的線性不穩定因子，采取相應措施，使模式的時間步長不受穩定性限制：（1）外重力波，采用Crank-Nicholson時間平均隱格式的θ參數法求解連續方程中的水平對流項和運動方程中的正壓梯度力項，以克服由外重力波所引起的CFL條件限制；（2）垂向粘性（擴散）項，對運動方程中的垂向粘性項和溫鹽方程中的垂向擴散項，采用半隱差分格式求解，以克服該兩項對穩定性的限制；（3）科氏力項，采用預報校正格式消除顯式處理運動方程中科氏力項引起的慣性不穩定性。通過以上處理，有效地提高了模式的時效性、穩定性和計算精度。

2 模型驗證與應用

2.1 模型設置

以象山灣內某電廠為例，研究三維溫排放模型及溫升影響。電廠排水口位于象山灣底部的強蛟東北方，西臨鐵港，東南為黃墩港（圖1）。象山灣為東北至西南向深入內陸的狹長型半封閉海灣，灣內主槽水深較大，一般為10～20 m，最大水深可達55 m，汊港及岸邊水域較淺，含有大片潮間灘地。電廠東側海域處于鐵港、黃墩港二港的分汊滯流區，灘涂發育，島礁較多，水深一般為5～10 m。該電廠三維溫排放模型東西向40 km，南北30 km，包含從象山灣口門的西澤潮位站至灣底的整個范圍。采用正交曲線網格進行剖分，水平網格總數為131×71個，最小網格尺寸194 m，最大356 m，平均240 m。垂向均勻分為5層，共6個層面。外海開邊界采用西澤潮位站實測潮位資料外推而得。

2.2 模型驗證

模型驗證采用2002年8月13日～30日連續18 d潮位以及一個連續潮汛大、中、小潮潮流的水文測驗資料。測站布設如圖1所示，強蛟和烏沙山為2個潮位測站，排水口附近的1#、2#為2條測流垂線，每條垂線分表層、0.2H、0.4H、0.6H、0.8H、底層6個測點流速值。根據給定的初邊值條件和計算參數，對象山灣潮流進行了驗證計算，潮位過程的驗證曲線見圖2，潮位計算值及相位均與實測值較吻合。圖3、圖4分別為測流點的大、小潮的表層、0.4H和底層流速、流向過程驗證曲線，可以看出，除了流速峰值略小，底層流向相位稍有偏差外，其余吻合良好。圖5為大潮漲、落急時刻表層和底層的流場圖，可以看到，漲、落潮流通道基本一致，都是沿深槽而行，流速分布也符合實測情況。中潮、小潮和大潮相比，除了流速量值上的差別之外，流場特性基本相同。

以上驗證結果說明，模型能較好地反映計算水域的三維水動力特性，在此基礎上耦合溫排水方程，可用來模擬溫排水的三維溫升分布。

2.3 溫升分布計算及分析

由于缺乏溫度場、鹽度場的實測資料，水體溫度的初始值采用海域的表層月平均值，為30℃，鹽度垂向分布較均勻，其初始值采用相應的月垂向平均值，為27‰。電廠溫排水的排放流量為85.5 m3/s，表層排放，連續均勻自由出流，出流溫升為8.6℃。將水動力與溫升耦合模擬一個月，溫升場基本達到穩定。考慮到中潮在潮差上具有平均意義上的代表性，選取中潮潮周期內計算的溫升場，分別繪制表層、底層的漲急、漲憩、落急、落憩四時刻的溫升等值線分布，分析溫升分布的三維特性（圖6）。

從圖6可以看到，受潮流影響，溫水排出后，在排放口附近形成了一個溫升帶。溫升帶分布范圍與漲落潮流主要流向相關。漲潮時，受潮流頂托作用，溫升帶主要分布在排水口以西，其中1℃以上的高溫帶沿岸線呈帶狀分布，漲憩時溫水舌最遠可伸至鐵港內部；落潮時，受潮流拖曳作用，溫升帶呈輻射狀分布在排水口以東，落憩時溫升帶變得狹長，離岸最遠。從垂向上看，由于溫差產生的浮力作用，表層排水口附近的溫度都在4℃以上，而底層溫度基本不超過2℃，表層明顯高于底層；而0～0.5℃的溫升范圍在垂向上相差不大，這是因為該溫升段主要分布在近岸的淺水區，水深不大的情況下，水體摻混均勻，分層效果不明顯。同時將該溫升分布計算結果與物理模型［9］試驗結果相比較，發現趨勢一致。

總體來說，溫排水影響范圍不大，溫升帶分布形態與潮流場密切相關，主要在以排水口為中心的地帶隨漲落潮擺動，最深探至鐵港分支這一帶，并未進入黃墩港。在排水口附近的，受浮力影響，溫升垂向分層較為明顯。

3 結論

將三維斜壓水流和溫排水數學模型用于象山灣內某電廠溫排放模擬，由于考慮了斜壓效應，計算密度時計及溫度變化，結果真實地反映了浮力作用下溫排水的三維運動特性，并且與物理模型試驗結果趨勢一致。海灣電廠溫升分布主要受漲落潮動力影響，形成以排水口為中心的溫水回蕩帶，在浮力作用下，排放口表層溫度高，底層溫度低。

［1］黃平.汕頭港水域溫排水熱擴散的三維數值模擬［J］.海洋環境科學，1996，15（1）：59-65.

HUANG P.Three-dimensional Numerical Simulation of the Heat Diffusion from the Cooling Water in Shantou Harbour［J］.Marine Environmental Science，1996，15（1）：59-65.

［2］王麗霞，孫英蘭，鄭連遠.三維熱擴散預測模型［J］.青島海洋大學學報，1998，28（1）：29-35.

WANG L X，SUN Y L，ZHENG L Y.A Three-Dimensional Prediction Method for Thermal Diffusion［J］.Journal of Ocean University of Qingdao，1998，28（1）：29-35.

［3］陶建峰.河口海岸三維斜壓水流數值模式研究［D］.南京：河海大學，2006.

［4］Phillips N A.A Coordinate System Having Some Special Advantages for Numerical Forecasting［J］.Journal of Meteorology，1957，14：184-185.

［5］Mellor G L.An Equation of State for Numerical Models of Oceans and Estuary［J］.Journal of Atmospheric and Oceanic Technology，1991，8：609-611.

［6］Smagorinsky J.General Circulation Experiments with the Primitive Equations I.The Basic Experiment［J］.Monthly Weather Review，1963，91：99-164.

［7］Mellor G L，Yamada T.Development of A Turbulence Closure Model for Geophysical Fluid Problems［J］.Reviews of Geophysics and Space Physics，1982，20（4）：851-875.

［8］陳惠泉，毛世民.水面蒸發系數全國通用公式的驗證［J］.水科學進展，1995，6（2）：116-120.

CHEN H Q，MAO S M.Calculation and Verification of an Universal Water Surface Evaporation Coefficient Formula［J］.Advances in Water Science，1995，6（2）：116-120.

［9］中國水利水電科學研究院.浙江國華寧海發電廠新建工程取排水口工程方案優化模型試驗研究報告［R］.北京：中國水利水電科學研究院，2003.YAO Jing，ZHANG Chang-kuan，TAO Jian-feng