《1 前言》

1 前言

湍流模型在求解工程中的复杂水流现象做出了重要贡献, 工程中最常用的就是k-ε双方程湍流模型, Launder, Rodi, Demuren, Ye Jian等做了很多有意义的工作 [1,2,3,4]。在很多情况下, k-ε湍流模型能够得出比较符合实际的结果, 但由于使用了Boussinesq的线性本构关系及标量的粘性系数概念, 无法考虑离心力与浮力在不同方向的作用, 因而用于模拟环流、强漩流与浮力分层流等各向异性的湍流时, 误差很大, 甚至定性失真。而在各相异性的湍流模拟中, 代数应力湍流模式精度较高, 而复杂程度又远小于应力通量湍流模型, 所以在近年来得到较快发展, 如Gibson、倪浩清、Ross和Rahman等均给出了较成功的算例 [5,6,7,8], 但这些算例都是针对规则边界的流场及标量场, 对于复杂边界的天然河流及海湾, 则受到很大限制。曲线坐标下环境水动力模型是研究污染物在天然水体扩散输移规律及在水体中的浓度分布的重要工具。

笔者建立了曲线坐标下平面二维水流计算和污染物扩散输移的代数应力湍流模型, 并用于实验室连续弯道污染物扩散输移的数值模拟, 计算了流场及岸边和中心污染物排放的浓度分布, 对该模型的计算值与笔者所建k-ε模型的计算值以及实测值进行了比较 [9], 结果表明代数应力湍流数学模型明显优于基于各向同性假设条件下的k-ε模型。

《2 数学模型》

2 数学模型

《2.1通用方程》

2.1通用方程

笛卡儿坐标下深度平均的代数应力湍流模型的通用微分方程为

t(ΗρΦ)+x(ΗρuΦ)+y(ΗρvΦ)=x(ΗΓΦΦx)+y(ΗΓΦΦy)+SΦ(1)

方程式 (1) 转换到非正交曲线坐标 (ξ, η) 下 [10], 仅在对流项中使用流速的逆变分量, 而在其他项中使用原始变量, 这样既简化了方程, 又使所有方程仍可写为曲线坐标下的通用方程, 模型的通用微分方程可写为如下通用形式:

t(ΗρΦ)+1Jξ(ΗρUΦ)+1Jη(ΗρVΦ)=1Jξ(ΗJΓΦαΦξ-ΗJΓΦβΦη)+1Jη(ΗJΓΦγΦη-ΗJΓΦβΦξ)+SΦ(ξη)(2)

式中Φ为所求问题的因变量;UV分别为Cartesian坐标下x, y流速u, v的逆变分量, 仅在对流项中出现;ΓΦα, ΓΦγ, ΓΦβ为扩散系数;SΦ (ξ, η) 为源项。模型的控制方程组如表1所示。

表1 通用方程中各变量

Table 1 Variable of general equations

 

《表1》


方程
Φ ΓΦα ΓΦγ SΦ (ξ, η)

连续
1 0 0 0

x-动量
U y2ηΓux+x2ηΓuy x2ξΓuy+y2ξΓux Su

y-动量
v y2ηΓvx+x2ηΓvy x2ξΓvy+yξ2Γvx Sv

浓度
C y2ηΓCx+x2ηΓCy x2ξΓCy+y2ξΓCx SC

湍流动能
K (y2ηΓkx+x2ηΓky) /σk (x2ξΓky+y2ξΓkx) /σk Sk

耗散率
ε (yη2Γεx+xη2Γεy) /σε (x2ξΓεy+y2ξΓεx) /σε Sε

 

 

表中:

Su=-1J[ρgΗ(yηzsξ-yξzsη)]-τxb-yηJξ(ρΗC11k)+yξJη(ρΗC11k)+xξJη[Η(μ+μty)J(yηvξ-yξvη)]-xηJξ[Η(μ+μty)J(yηvξ-yξvη)]-1Jξ{ΗΓuβJuη}-1Jη{ΗΓuβJuξ}(3)Sv=-1J[ρgΗ(xξzsη-xηzsξ)]-τyb-xξJη(ρΗC12k)+xηJξ(ρΗC12k)+yηJξ[Η(μ+μtx)J(xξuη-xηuξ)]-yξJη[Η(μ+μtx)J(xξuη-xηuξ)]-1Jξ(ΗΓvβJvη)-1Jη(ΗΓvβJvξ)(4)Sk=yηJξ[ΗΓkλJ(xξkη-xηkξ)]-yξJη[ΗΓkλJ(xξkη-xηkξ)]+xξJη[ΗΓkλJ(yηkξ-yξkη)]-xηJξ[ΗΓkλJ(yηkξ-yξkη)]-1Jξ(ΗΓkβJkη)-1Jη(ΗΓkβJkξ)+Ηρ(Ρk-Ρkv-ε)(5)Sε=yηJξ[ΗΓελJ(xξεη-xηεξ)]-yξJη[ΗΓελJ(xξεη-xηεξ)]+xξJη[ΗΓελJ(yηεξ-yξεη)]-xηJξ[ΗΓελJ(yηεξ-yξεη)]-1Jξ(ΗΓεβJεη)-1Jη(ΗΓεβJεξ)+Ηρ[εk(Cε1Ρk-Cε2ε)-Ρεv)](6)SC=SC0+1Jξ{1Cφ1kεuv¯[2CξxηyηJ-Cη(xξyηJ+xηyξJ)]}+1Jη{1Cφ1kεuv¯[2CηxξyξJ-Cξ(xηyξJ+xξyηJ)]}+1Jξ{(1-Cφ2)[uφ¯(uξyη2J-uηyξyηJ)+vφ¯(uηxξyηJ-uξxηyηJ)]}-1Jη{(1-Cφ2)[uφ¯(uξyξyηJ-uηyξ2J)+vφ¯(uηxξyξJ-uξxηyξJ)]}-1Jξ{(1-Cφ2)[uφ¯(vξxηyηJ-vηxηyξJ)+vφ¯(vηxξyηJ-vξxη2J)]}+1Jη{(1-Cφ2)[uφ¯(vξxξyηJ-vηxξyξJ)+vφ¯(vηxξ2J-vξxξxηJ)]}(7)ΓΦα=ΓΦxyη2+ΓΦyxη2,ΓΦγ=ΓΦyxξ2+ΓΦxyξ2,ΓΦβ=ΓΦxyξyη+ΓΦxxξxη(Φ=u,v,C,k,ε)Γux=2μ+2μt0+ρDu,Γuy=μ+μtx+ρDu,Γvx=μ+2μty+ρDv,Γvy=2μ+2μt0+ρDu,ΓCx=μ+ρuu¯k/Cφ1ε,ΓCy=μ+ρvv¯k/Cφ1ε,Γkx=μ+ρCkuu¯k/ε,Γky=μ+ρCkvv¯k/ε,Γεx=μ+ρCεuu¯k/ε,Γεy=μ+ρCεvv¯k/ε,Γkλ=ρCkuv¯k/ε,Γελ=ρCεuv¯k/ε,μt0=ρCμ0k2/ε,U=uyη-vxη,V=vxξ-uxξ,α=xη2+yη2,β=xξxη+yξyη,γ=xξ2+yξ2,J=xξyη-yξxη,μty=ρCμ0k2ε[1-2(1-C2)C1kε(uξyηJ-uηyξJ)]μtx=ρCμ0k2ε[1-2(1-C2)C1kε(vηyξJ-vξyηJ)]Ρk=-ρJ[uu¯(uξyη-uηyξ)+uv¯(uηxξ-uξxη)+uv¯(vξyη-vηyξ)+vv¯(vηxξ-vξxη)]Ρkv=caρU*3/Η,Ρεv=cbρU*4/Η2,U*=(cf(u2+v2))1/2,Η=zs-zb,Du=Dv=0.2ΗU*,ca=cf-1/2cb=3.6Cε2cf-0.75Cμ1/2,cf=0.003,τbx=gn2u(u2+v2)1/2/Η1/3,τby=gn2v(u2+v2)1/2/Η1/3,

其中, SC0为污染物源强, kε分别为深度平均的湍流动能和湍流动能耗散率, H为水深, zs为水位, zb为河床高程, μ为分子动力粘性系数, -uu¯-vv¯,-uv¯分别代表u, v及交叉方向的Reynolds应力, -uφ¯-vφ¯为方向的通量。

《2.2雷诺应力》

2.2雷诺应力

曲线坐标系下雷诺应力表达式:

-uu¯=2Cμ0k2ε(uξyηJ-uηyξJ)-C11k(8)

-vv¯=2Cμ0k2ε(vηxξJ-vξxηJ)-C12k(9)-uv¯=Cμ0k2ε[1-2(1-C2)C1kε(vηxξJ-vξxηJ)](uηxξJ-uξxηJ)+Cμ0k2ε[1-2(1-C2)C1kε(uξxηJ-uηyξJ)](vξyηJ-vηyξJ)(10)uφ¯=1Cφ1kε{(1-Cφ2)[uφ¯(uξyηJ-uηyξJ)+vφ¯(-uξxηJ+uηxξJ)]-uu¯(CξyηJ-CηyξJ)-uv¯(-CξxηJ+CηxξJ)}(11)vφ¯=1Cφ1kε{(1-Cφ2)[uφ¯(vξyηJ-uηyξJ)+vφ¯(-vξxηJ+uηxξJ)]-uv¯(-CξyηJ-CηyξJ)-vv¯(-CξxηJ+CηxξJ)}(12)Cμ0=13(1-C2)(1-C2+2C1)/{C12+(1-C2)2k2ε2[(uηxξJ-uξxηJ)-(vξyηJ-vηyξJ)]2}(13)C11=13{(1-C2+2C1)C1+C13(uηxξJ-uξxηJ)[(uηxξJ-uξxηJ)-(vξyηJ-vηyξJ)]Κ¯2ε¯2}/{C12+(1-C2)2k2ε2[(uηxξJ-uξxηJ)-(vξyηJ-vηyξJ)]2}(14)

C12=13{(1-C2+2C1)C1-C13(vξyηJ-vηyξJ)[(uηxξJ-uξxηJ)-(vξyηJ-vηyξJ)]k2ε2}/

{C12+(1-C2)2k2ε2[(uηxξJ-uξxηJ)-(vξyηJ-vηyξJ)]2}(15)C13=2(1-C2)2(1-C2+2C1)/C1(16)

式中 Ck=0.24, C1=2.2, C2=0.55, Cε=0.15, Cε1=1.44, Cε2=1.92, Cφ1=3.0, Cφ2=0.5, σk=1.0, σk=1.3为湍流经验常数。

《3 模型离散和求解》

3 模型离散和求解

模型的离散采用控制体积法, 控制方程的系数采用乘方格式, 模型求解采用SIMPLEC算法。针对代数应力湍流模型中湍流动能源项容易出现负值而产生发散的问题, 采用文献[11]中的处理办法。

《3.1边界条件》

3.1边界条件

在进口边界, 所有边界条件都按本征条件给出, 即u=u0, v=v0, k=k0, ε=ε0

在出口边界, 给定水位zs, 并设各变量法向导数为零, 即zsn=un=vn=kn=εn=0

在壁面边界, 采用壁面函数法 [12], 即用半分析的方法得到的解来近似由壁面到湍流核心区之间的流速、湍流动能、湍流动能耗散率的分布规律, 将壁面的影响 (如壁面应力) 附加到差分方程中 (应先将有关的边界系数置为零) 。

《3.2SIMPLEC算法》

3.2SIMPLEC算法

Step 1 计算坐标变换的相关系数;

Step 2 根据进、出口边界的水位, 给定全场水位, 计算全场初始水深H*;

Step 3 解动量方程, 求un, vn;

Step 4 解代数应力湍流模型, 求k, ε;

Step 5 解水深校正方程, 求H′, 并修正水深H=H*+αhH′和流速场u, v;

Step 6 求解各个方向的雷诺应力-uu¯-uv¯-vv¯-uφ¯-vφ¯及各个方程、各个方向的扩散系数;

Step 7 重复Step 3至Step 6直至收敛;

Step 8 求解浓度C方程, 直至得到稳定解。

模型求解中, 为了利于非线性迭代的收敛, 计算中采用亚松弛技术, ADI技术和TDMA算法, 以连续方程的残余质量小于一给定值 (10-5) 作为判断收敛的依据。

《4 模型应用》

4 模型应用

模型采用Chang 在实验室连续弯道中进行的一系列水流和污染物扩散输移的试验结果作为验证资料 [13]。弯道断面为矩形并具有光滑底面, 其形状与具体尺寸如图1所示。在试验中, 污染物扩散输移按三种情况进行, 即中心排放, 岸边排放 (左岸排放和右岸排放) , 在第一个弯道进口处排放。水深H=0.115 m, 进口流量Q=0.09486 m3/s, 污染物初始浓度为零。模型将物理区域划分为166×45个网格, 时间步长取Δt=0.2 s, 糙率取0.015。计算至1 200步, 计算收敛值收敛至规定值。选了4个断面浓度的模型计算值与k-ε模型计算值以及实测值进行比较 [9], 见图2。

《图1》

图1 实验室连续弯道示意图

图1 实验室连续弯道示意图  

Fig.1 Sketch of meandering channel in lab

《图2》

图2 断面浓度比较

图2 断面浓度比较  

Fig.2 Comparison of section concentration

从图2的比较, 可以看出本模型计算的浓度在横向的扩散明显大于基于各向同性假定的k-ε模型的计算值, 调整Schmidt数σC可以加大浓度的横向扩散, 但在没有实测值的情况下, σC的调整是没有根据的。连续弯道中污染物浓度的横向扩散还受到弯道二次流的加强, 所以, 平面二维模型在计算曲率比较大的弯曲河流时, 缺陷比较明显。

《5 结论》

5 结论

建立了曲线坐标下平面二维污染物扩散输移的代数应力湍流模型, 以对具有不规则边界水域的湍流各向异性特征进行模拟。将模型用于实验室连续弯道的污染物扩散输移的数值模拟, 分别计算了岸边和中心污染物排放的浓度分布, 对模型的计算值与笔者所建k-ε湍流模型的计算值以及实测值进行了比较 [9], 结果表明代数应力湍流数学模型明显优于基于各向同性假设条件下的k-ε模型。由于平面二维模型在计算曲率比较大的弯曲河流时具有比较大的缺陷, 需要发展曲线坐标下三维弯曲河流水流和污染物扩散输移的湍流模型。