浓度及温度场问题三维数值模拟及应用

，

(长江科学院 河流研究所，武汉 430010)

1 国内外研究现状

随着全球范围内工农业的快速发展，越来越多的污染物被排放到河流、湖泊和大海，从而导致其污染的发生。如何解决污染问题受到了人们的重视。

20世纪80年代以来，国内外学者相继建立物理和数学模型，对水体中污染物的对流和扩散规律进行了研究，并有一些数学模型被成功应用到工程领域。如Demuren[1]进行了明渠的侧向排放实验；Demuren和Rodi[2]采用数值模拟方法研究了污染物在连续弯曲水槽中的运动规律；董耀华[3]采用平面二维数学模型研究了河道水流、温度及浓度输运的规律，并采用该模型对实验室空腔水槽温度场、含沙量场、长江口南支河段温排水及盐水入侵进行了研究；林彬良等[4-5]建立了基于河口和海湾的平面二维模型对英国的Humber湾进行了污染物浓度场的计算；Sladkvich等[6]进行了平面二维浅水水生态环境模拟；Rueda-Valdivia[7]建立了三维湖泊模型，并对Oaks湖进行水温模拟；沈永明等[8]建立了考虑盐度场、温度场、浓度场等相互耦合作用的三维湍流模型，并对大连湾的水环境进行了模拟；许金电等[9]通过建立σ坐标系下的Navier-Stoker方程及三维物质扩散模型，对浅海三维污染物扩散规律进行了模拟和研究；崔丹等[10]建立了基于三角形网格有限体积法的三维水流运动与温度输运数学模型，模拟了温排水的时空分布及取水口的温升过程线。

随着对电力的大量需求，相继修建了众多热电厂，因而电厂温排水带来的环境及生态问题也越来越严重。

美国学者从20世纪70年代就开始研究温水排放引起的流场和水质变化：Beer等对Michigan湖的温水排放进行了环境影响研究，1975年Louisville大学开始预测温水排放的温度分布规律。此后如何预测和控制热(核)电厂热污染的问题，引起了世界范围内广泛的关注[11-12]。

我国学者从20世纪80年代开始相关研究，到现在己经有了比较长足的发展。倪浩清等[13]建立了大尺度水域深度平均的湍流水环境模式，并应用在温排水的模拟;2002年陈华[14]利用平面二维ADI(Alternating Direction Implicit)数学模型对后石电厂海区做了流场和温度场的数值模拟；郝瑞霞等[15]采用平面二维水流和温度场数学模型，对拟建的深圳前湾电厂冷却水工程进行了模拟研究；朱军政[16]选用平面二维水动力和物质输运数值模型对强潮海湾进行了流场、温度场的计算，并分析了温排水温升包络面积的影响因素、取水口温升变化的特点，同时分析了网格尺寸、综合海面散热系数对温升的影响。在三维温度场研究方面，杨纯等[17]采用简化的三维模型，预测辽西(高岭)电厂温排水受纳水域的流场和温度场；黄平[18]提出了一种准三维的温排水数值模型；1997年王丽霞等[19]建立了三维热扩散预测模型，对黄岛电厂温水排放进行了模拟计算；2004年郝瑞霞等[20]采用浮力修正的湍流模型，准三维离散型边界拟合坐标变换下的控制体积法，进行了滨海电厂冷却水工程的潮汐水流和热传输的数值模拟；2006年王一航等[21]基于POM模式，对北海近岸的庄河电厂附近海域进行了潮汐潮流和温排水三维数值模拟;2008年张细兵等[22]采用结构和非结构网格，建立了三维温度场和浓度场数学模型，并应用于天然河道浓度及温度场问题中。

综合以上研究进展可知，国内在温排水和污染物浓度场数学模型方面主要还是平面二维模型，采用三维模型求解的较少，而受平面二维模型本身局限性影响，在模拟热源或污染源近区时不能反映温度或浓度场的分层现象，会产生较大的误差，因此采用三维模型来研究是很有必要的。

2 模型的建立

2.1 基本方程

天然情况下，计算区域的边界都是不规则的，对于不规则边界，阶梯形的逼近方法在网格数较少时过于粗糙，网格数增加又会导致计算经济性恶化，采用坐标变换，使得流体控制方程相应地复杂化，也大大增加了计算量。为此，本文引入空度(流体可流通区域在整个积分区域中所占的比例)概念处理不规则边界及河床变形。在网格化了的积分区域中，各个网格可以有不同的空度。引入空度和VOF之后的方程组如下所述。

水流连续方程及动量方程：

(1)

(2)

(3)

其中：k，ε分别为紊动能及其扩散率；δij为Kronecker delta符号(当i=j时，δij=1，否则δij=0)；νt可通过RNG模型[23]求出。

(4)

(5)

RNG模型中各参数取值如表1和式(6)所示。表中σk,σε,C2ε。

表1 RNG紊流模型参数表

(6)

引入空度后的VOF方程：

(7)

式中，F为流体体积百分数。

浓度场及温度场方程：

(8)

(9)

式中：νe=ν+νt；C为浓度；ω为沉速；I为源汇项；T为温度；Cp为水的热容系数；HT为热交换项。

2.2 方程的离散及求解

采用有限体积法(FVM)离散方程组，为避免压力场的震荡，网格剖分采用交错网格法。为保证数值解的稳定性，按照负坡线形化原则对各方程源项进行处理。离散方程组的求解采用SIMPLE方法。

2.3 边界条件

(1) 进出口边界：进出口边界由已知条件给出。

(2) 固壁边界：采用对数律壁面函数处理，具体可参考文献[24]。

(3) 自由表面边界：自由面采用VOF法计算。自由表面位置通过求解液相体积百分数传输方程确定，即式(7)，这样通过求解每个单元的液相体积百分数，就可确定自由表面的位置。

(4) 为避免普通的差分格式带来的数值扩散导致界面模糊现象，本文在求解时采用Hirt和Nichols的“给予-接受”方法，即DA法。

(5) 水体表面散热系数为

(10)

(6) 水体密度为

(11)

式中：t为水面温度(℃)；Uz为风速(m/s)。

3 模型验证及应用

3.1 例1——水流模型检验

Muneta和Shimizu于1994年在实验室水槽中进行了丁坝绕流试验[15]，并详细测量了流场及流速分布。试验流量为1.87 L/s，水槽宽为0.4 m，槽底坡度为1‰，平均水深0.07 m，平均流速约0.067 m/s。丁坝为长0.2 m、宽0.04 m的长方体，为非淹没丁坝。

计算区域选在丁坝轴线上游1.5 m与下游1.5 m之间的3.0 m×0.4 m×0.1 m长方体范围内，如图1所示。网格节点数为80×42×15，采用非均匀划分，底面、侧面、丁坝附近及自由水面附近网格局部加密。

图1 计算区域平面布置

图2为C-C横断面的横向流速计算与实测的比较图，图3分别为A-A纵剖面、C-C横断面及D-D横断面流速矢量图。

注：图中虚直线为速度等于零的位置

图3 A-A剖面、C-C横断面、D-D横断面流速矢量图

从图2可见，计算和实测结果较吻合。图3(a)可见，A-A纵剖面上在丁坝附近有一顺时针方向的涡；从图3(b)、图3(c)可见本文计算和实测及文献[12]计算的流场很一致，整个横断面上有2个涡存在，除y=30～40 cm范围内，有一个靠近水面的逆时针方向较大的涡外，还有在y=20～25 cm范围内，靠近底面处有一个逆时针方向的小涡，Muneta和Shimizu在实验中观察到了这个涡，而这个小涡是水流绕过丁坝后三维特性的重要特征。

3.2 例2——浓度场模型检验

假定有一断面呈矩形，深、宽均为300 m，长为2 800 m的顺直水道，水道仅在纵向存在流速，该流速处处相同且恒定不变，在该水道的上游断面中心，存在一稳定点源连续排放，水中的扩散速率各向同性且为常数，此时有理论解

假定M=400 g/s，u=0.2 m/s，D=0.1 m2/s，Δx=100 m,Δy=Δz=100 m,Δt=300 s,点源所在网格的源项为M/(Δx·Δy·Δz)。

图4给出了计算和理论值的通过排放点的x-y剖面的浓度等值线的对比情况，从图上看，吻合较好。只是在距离排放口比较近的地方，小浓度的计算值有所差异。

图4 计算与理论值浓度等值线对比图

3.3 实3——温度场模型检验

所选计算区域为一模型水库，是美国陆军兵工团于1980—1981年为验证水库水温动力学模型而做的模型实验。该模型水库长为24.39 m，宽度从入口处的0.3 m逐渐增大到出口的0.91 m。水库由2部分组成，上部分长6.1 m，库底水平，宽度由0.3 m线性变化到0.91 m；第2部分长为18.29 m，宽度保持不变为0.91 m，底部下降了0.61 m。水库几何尺寸见图5。

图5 水库尺寸示意图

计算条件：初始时刻水库库区内充满21.44℃的热水，16.67℃冷水从入口进入水库，入流流量为0.000 33 m3/s，入口断面上半部分为挡板，冷水从挡板以下进入，出口于出口断面0.15 m处，出流流量与入流流量相等。

数值模拟中采用结构化网格，网格节点数为80×10×25，最大网格尺寸0.45 m，最小网格尺寸0.005 m，平均网格尺寸为0.025 m。数值求解过程中假定初始值，先求解水动力学方程，随后求解能量方程和紊动方程，将计算结果作为初场进行迭代计算，直到达到收敛精度，再进行下一时间步的计算。

图6 计算与实测水库中心纵向流速分布

图6给出了模型水库中心纵向流速分布，图7给出了4 min和50 min时数值模拟的水库中温度场的分布。冷水进入水库后，其密度比热水大，基本上沿库底流动，热水在浮力作用下不断上升，在冷热水接触界面附近区域很快形成回流区。沿库底流动冷水的形态类似于层流，冷水上层为回流区，在紊动扩散作用下冷热水相掺混，出现明显的温度分层。50 min后，回流区充满整个水库，温度等值线上移。

图7 4 min和50 min时水库中冷热水速度场及温度场的分布

3.4 实4——温度场和浓度场实际检验

为检验模型在天然情况下的适应情况，选取长江中游的宜都河段进行计算。拟建宜都电厂位于长江宜都河段下段右岸弯顶部位，枝城水文站上游约4 km处。

(1) 水流条件：本次计算选取多年冬季平均流量进行计算。根据宜昌水文站1950—2003年流量进行统计，并考虑清江流量的汇入，多年冬季平均流量约为5 000 m3/s，计算河段出口水位根据枝城站水位推算。

(2) 工程条件：排水流量为13.8 m3/s，温升15℃。氯离子浓度0.3 mg/m3。

(3) 其它条件，水温10.9℃，气温7.9℃，相对湿度88%，可见度70%。

注：图中的纵、横坐标为大地坐标，单位均为m。

图8为表层与底层温升等值线平面分布。由图8中可见，温排水影响主要位于排水口下游的近岸区域，沿程温升不断减小，且越靠近排水口温升减小梯度越大。表层与底层温升图对比来看，表层图中高温升等值线范围明显大于底层，而低温升等值线差别不大，反映了近区存在温升分层，而远区经充分掺混后分层不明显的一般规律。计算结果统计表明，1℃温升线沿岸扩散长度为550 m，离岸扩散宽度为100 m。两图中均较好反映了温升在排水口附近沿垂线分层现象。

图9(a)为离子浓度横断面图，图9(b)为纵剖面图。图10中均反映了离子浓度在排水口附近沿垂线减小，在离排水口下游100 m以外，浓度沿垂线已接近于均匀分布。

图9 横断面、纵剖面离子浓度分布

4 结 语

在总结国内外研究的基础上，建立了河道浓度和温度场的三维数学模型，同时进行了概化模型计算和天然河道浓度及温度场计算与验证。

计算结果表明，建立的模型能较好模拟浓度及水温沿垂向的分层现象。由实测资料的验证表明，模型计算结果具有较高的精度，可在电厂温排水影响及污染物排放问题研究中推广应用，同时也可应用于湖泊、水库、河口海岸等水域的计算。

参考文献：

[1] DEMUREN A O, ROID W. Side Discharge in Open Channels: Mathematic Model[J]. Journal of Hydraulic Engineering, 1983, 109(12):1707-1723.

[2] DEMUREN A O, ROID W. Calculation of Flow and Pollutant Dispersion in Meandering Channels [J]. Journal of Fluid Mechanics, 1986, 172(11):65-72.

[3] 董耀华.河道温度及浓度场的平面二维数值模拟[J].长江科学院院报,1997,14(2):23-26. (DONG Yao-hua. Two-Dimensional Numerical Simulation of Temperature or Concentration Field in River Channels[J]. Journal of Yangtze River Scientific Research Institute, 1997,14(2):23-26. (in Chinese))

[4] LIN B L, CHANDLER S N. A Depth-integrated 2D Costal and Estuarine Model with Conformal Boundary-fitted Mesh Generation [J]. International Journal for Numerical Method in Fluids, 1996,23(8):819-846.

[5] LIN B L, SHIONO K. Numerical Modeling of Solute Transport in Compound Channel Flows [J]. Journal of Hydraulic Research, 1995, 33(6): 773-788.

[6] SLADKEVICH M, MILITEEV A N, RUBIN H,etal. Simulation of Transport Phenomena in Shallow Aquatic Environment [J]. Journal of Hydraulic Engineering, ASCE, 2000, 126(2):123-136.

[7] RUEDA-VALDIVIA F J. A Three-dimensional Hydrodynamic and Transport Model for Lake Environments[D]. USA: University of California, 2001.

[8] 沈永明,郑永红,吴修广.近岸海域污染物迁移转化的三维水质动力学模型[J]. 自然科学进展,2004,14(6):694-699. (SHEN Yong-ming, ZHENG Yong-hong, WU Xiu-guang. The 3D Water Quality Numerical Modeling about Pollutant Transport at Offshore Waters[J]. Progress in Natural Science, 2004,14(6):694-699. (in Chinese))

[9] 许金电,江毓武.港湾三维水动力和污染物扩散数值模型[J]. 台湾海峡,2003,22(1):85-91. (XU Jin-dian, JIANG Yu-wu. A Three-dimensional Coastal Hydrodynamic and Pollutant Diffusion Model[J]. Journal of Oceanography in Taiwan Strait, 2003,22(1):85-91. (in Chinese))

[10] 崔 丹,金 峰.近岸海域温排水的三维数值模拟[J]. 长江科学院院报,2010,27(10):55-59. (CUI Dan, JIN Feng. Three-dimensional Numerical Simulation on Thermal Discharge in Coastal Area[J]. Journal of Yangtze River Scientific Research Institute, 2010,27(10):55-59. (in Chinese))

[11] BOULET F. Modeling of Heated Water Discharges on the French Coast of the English Channel[C]∥ Transport Models for Inland and Coastal Waters. New York: Academic Press, 1981.

[12] HAMRICK J M, MILLS W B. Analysis of Water Temperature in Conowingo Pond as Influenced by the Peach Bottom Atomic Power Plant Thermal Discharge[J]. Environmental Science and Policy, 2000, 3(Sup.1): 197-209.

[13] 倪浩清,沈永明,陈惠泉.深度平均的k-ε紊流全场模型及其验证[J].水利学报,1994, (11):8-17. (NI Hao-qing, SHEN Yong-ming, CHEN Hui-quan. The Depth-averaged k-ε Turbulence Model for the Entire Water Environment and Its Verification[J]. Shuili Xuebao, 1994, (11):8-17. (in Chinese))

[14] 陈 华.后石电厂温排水的数学模型研究[D]. 厦门：厦门大学,2002. (CHEN Hua. Numerical Model of Thermal Discharge from Houshi Power Plant[D]. Xiamen: Xiamen University, 2002. (in Chinese))

[15] 郝瑞霞,齐 伟,李海香,等.潮汐水域流速场和温度场的数值模拟研究[J].太原理工大学学报,2005,36(3):235-237,245. (HAO Rui-xia, QI Wei, LI Hai-xiang,etal. Application of Matpower in Calculating Spot Price of Power System with Optimal Power Flow[J]. Journal of Taiyuan University of Technology, 2005,36(3):235-237,245. (in Chinese))

[16] 朱军政.强潮海湾温排水的数值模拟[J].水动力学研究与进展(A辑),2007,22(1):17-23. (ZHU Jun-zheng. Numerical Simulation of 2-D Model of Cooling Water Draining into the Bay with Strong Tide[J]. Chinese Journal of Hydrodynamics, 2007,22(1):17-23. (in Chinese))

[17] 杨 纯,傅龙泉,黄东涛,等.电厂温排水受纳水域传热与流动的数值计算[J].上海电力学院学报,1995, 11(2):1-9. (YANG Chun, FU Long-quan, HUANG Dong-tao,etal. Numerical Analysis of Heat Transfer and Fluid Flow in the Power Plant Thermal Discharge Waters[J]. Journal of Shanghai Institute of Electric Power, 1995, 11(2):1-9. (in Chinese))

[18] 黄 平.汕头港水域温排水热扩散的三维数值模拟[J].海洋环境科学,1996,(1):59-65. (HUANG Ping. Three-dimensional Numerical Simulation of the Heat Diffusion from the Cooling Water in Shantou Harbour[J]. Marine Environmental Science, 1996,(1):59-65. (in Chinese))

[19] 王丽霞,孙英兰，田 晖.热扩散预测方法研究概况II-热扩散的研究现状[J].海洋科学,1997，(6)：16-17. (WANG Li-xia, SUN Ying-lan, TIAN Hui. Survey of Predicting Method of Coastal Thermal Diffusion Ⅱ:Review of Thermal Diffusion[J]. Marine Sciences, 1997,(6):16-17. (in Chinese))

[20] 郝瑞霞，韩新生.潮汐水域电厂温排水的水流和热传输准三维数值模拟[J].水利学报，2004,(8):66-70. (HAO Rui-xia, HAN Xin-sheng. Quasi 3-D Numerical Model for Flow and Heat Transport of Power Plant Cooling Water Discharged into Coastal Area[J]. Shuili Xuebao, 2004,(8):66-70. (in Chinese))

[21] 王一航,魏泽勋,王永刚,等.潮汐潮流三维数值模拟在庄河电厂温排水问题中的应用[J].海洋通报,2006,25(1):8-15. (WANG Yi-hang, WEI Ze-xun, WANG Yong-gang,etal. Application of Tidal Numerical Model to Thermal Plume Dispersion from Power Plant at Zhuanghe[J]. Marine Science Bulletin, 2006,25(1):8-15. (in Chinese))

[22] 张细兵,崔占峰.浓度及温度场问题三维数值模拟报告(中央级公益性科研院所基本科研项目)[R].武汉：长江科学院，2008. (ZHANG Xi-bing, CUI Zhan-feng. Report on 3-D Numerical Simulation of Temperature and Concentration Field[R]. Wuhan: Yangtze River Scientific Research Institute, 2008. (in Chinese))

[23] SPEZIALE C G, THANGMA S. Analysis of a RNG Based Turbulence Model for Separated Flows[J]. International Journal of Engineering Science, 1992, 30(10): 1379-1388.

[24] 陶文铨.数值传热学[M].西安：西安交通大学出版社, 1988. (TAO Wen-quan. Numerical Heat Transfer Theory[M]. Xi’an: Xi’an Jiaotong University Press, 1988. (in Chinese))