海湾电厂三维斜压水流和温排水数值模拟

姚 静，张长宽，陶建峰

（河海大学港口海岸与近海工程学院水文水资源与水利工程科学国家重点实验室，南京 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