海灣電廠三維斜壓水流和溫排水?dāng)?shù)值模擬

姚 靜，張長(zhǎng)寬，陶建峰

（河海大學(xué)港口海岸與近海工程學(xué)院水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，南京 210098）

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

本文以象山灣內(nèi)某電廠的溫排水為例，進(jìn)行三維斜壓水流和溫排水?dāng)?shù)值模擬。先進(jìn)行水動(dòng)力場(chǎng)和溫度場(chǎng)計(jì)算，將溫度場(chǎng)計(jì)算結(jié)果代入狀態(tài)方程，實(shí)時(shí)修正水體密度，再將密度反饋到斜壓項(xiàng)，進(jìn)行流場(chǎng)更新計(jì)算，進(jìn)而得到新的溫度場(chǎng)，依此循環(huán)，實(shí)現(xiàn)水動(dòng)力場(chǎng)和溫度場(chǎng)的耦合計(jì)算。在根據(jù)實(shí)測(cè)資料對(duì)模型進(jìn)行充分驗(yàn)證的基礎(chǔ)上，模擬了電廠溫排水?dāng)U散輸移過程中的三維溫升分布特征。

1 三維斜壓水流和溫排水?dāng)?shù)學(xué)模型

1.1 控制方程

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

（1）水動(dòng)力方程

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

（2）溫排水方程

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

（3）狀態(tài)方程

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

式中：S為鹽度；p為壓力；ρ（S，T，0）為一個(gè)標(biāo)準(zhǔn)大氣壓（海壓為0）下的海水密度狀態(tài)方程；c由式（6）計(jì)算。

1.2 定解條件

初始條件以零啟動(dòng)的形式給出，即全場(chǎng)3個(gè)流速分量給定為零；同時(shí)給出初始水位和初始溫升場(chǎng)，其值根據(jù)電廠附近水域?qū)崪y(cè)資料插值給出。

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

1.3 主要參數(shù)

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

1.4 計(jì)算方法

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

2 模型驗(yàn)證與應(yīng)用

2.1 模型設(shè)置

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

2.2 模型驗(yàn)證

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

以上驗(yàn)證結(jié)果說明，模型能較好地反映計(jì)算水域的三維水動(dòng)力特性，在此基礎(chǔ)上耦合溫排水方程，可用來(lái)模擬溫排水的三維溫升分布。

2.3 溫升分布計(jì)算及分析

由于缺乏溫度場(chǎng)、鹽度場(chǎng)的實(shí)測(cè)資料，水體溫度的初始值采用海域的表層月平均值，為30℃，鹽度垂向分布較均勻，其初始值采用相應(yīng)的月垂向平均值，為27‰。電廠溫排水的排放流量為85.5 m3/s，表層排放，連續(xù)均勻自由出流，出流溫升為8.6℃。將水動(dòng)力與溫升耦合模擬一個(gè)月，溫升場(chǎng)基本達(dá)到穩(wěn)定?？紤]到中潮在潮差上具有平均意義上的代表性，選取中潮潮周期內(nèi)計(jì)算的溫升場(chǎng)，分別繪制表層、底層的漲急、漲憩、落急、落憩四時(shí)刻的溫升等值線分布，分析溫升分布的三維特性（圖6）。

從圖6可以看到，受潮流影響，溫水排出后，在排放口附近形成了一個(gè)溫升帶。溫升帶分布范圍與漲落潮流主要流向相關(guān)。漲潮時(shí)，受潮流頂托作用，溫升帶主要分布在排水口以西，其中1℃以上的高溫帶沿岸線呈帶狀分布，漲憩時(shí)溫水舌最遠(yuǎn)可伸至鐵港內(nèi)部；落潮時(shí)，受潮流拖曳作用，溫升帶呈輻射狀分布在排水口以東，落憩時(shí)溫升帶變得狹長(zhǎng)，離岸最遠(yuǎn)。從垂向上看，由于溫差產(chǎn)生的浮力作用，表層排水口附近的溫度都在4℃以上，而底層溫度基本不超過2℃，表層明顯高于底層；而0～0.5℃的溫升范圍在垂向上相差不大，這是因?yàn)樵摐厣沃饕植荚诮兜臏\水區(qū)，水深不大的情況下，水體摻混均勻，分層效果不明顯。同時(shí)將該溫升分布計(jì)算結(jié)果與物理模型［9］試驗(yàn)結(jié)果相比較，發(fā)現(xiàn)趨勢(shì)一致。

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

3 結(jié)論

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

［1］黃平.汕頭港水域溫排水熱擴(kuò)散的三維數(shù)值模擬［J］.海洋環(huán)境科學(xué)，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］王麗霞，孫英蘭，鄭連遠(yuǎn).三維熱擴(kuò)散預(yù)測(cè)模型［J］.青島海洋大學(xué)學(xué)報(bào)，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］陶建峰.河口海岸三維斜壓水流數(shù)值模式研究［D］.南京：河海大學(xué)，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］陳惠泉，毛世民.水面蒸發(fā)系數(shù)全國(guó)通用公式的驗(yàn)證［J］.水科學(xué)進(jìn)展，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］中國(guó)水利水電科學(xué)研究院.浙江國(guó)華寧海發(fā)電廠新建工程取排水口工程方案優(yōu)化模型試驗(yàn)研究報(bào)告［R］.北京：中國(guó)水利水電科學(xué)研究院，2003.YAO Jing，ZHANG Chang-kuan，TAO Jian-feng