非正交曲线坐标下二维水流计算的 SIMPLEC

2003年2月
水 利 学 报SH UI LI X UE BAO 第2期
收稿日期:2002201204
基金项目:国家杰出青年科学基金项目(50125924);高等学校博士学科点专项科研基金项目(2000014112);辽宁省自然科学基
金项目(2001101073)
作者简介:吴修广(1974-),男,山东阳谷人,博士生,从事环境水力学研究。

文章编号:055929350(2003)022*******非正交曲线坐标下二维水流计算的SIMPLEC 算法
吴修广1,沈永明1,郑永红2,王平义
3(11大连理工大学海岸和近海工程国家重点实验室　辽宁大连　116023;21中国科学院广州能源研究所,广东广州　510070;
31重庆交通学院河海工程系　重庆　400074)摘要:本文采用Laplace 方程坐标变换方法生成正交曲线网格,并对浅水流动的控制方程进行坐标变换,方程离散时采用B 型交错网格。

利用“水位扫描法”结合壁面函数法来处理移动边界,用SI MP LEC 算法解非正交曲线坐标下的k -ε双方程紊流模型,修正了由网格的非正交性引起的误差。

通过对美国C olorado 洲Fall River 的资料进行流场验证,计算结果与实测资料基本符合,显示了本模型在不规则水域计算中的实用价值。

关键词:坐标变换;k -ε紊流模型;水位扫描法;壁面函数;SI MP LEC 算法
中图分类号:T V13114文献标识码:A
随着经济发展和社会进步,水利工程建设的步伐也在进一步加快,其中港航建设、大坝建设中的泥沙问题以及近来倍受世人关注的水污染问题已经成为制约水利发展的瓶颈问题,弄清河流、湖泊、海洋中水动力因素,是解决以上问题的重要基础。

近年来,数学模型已逐步取代物理模型实验成为研究水流的重要手段,而浅水流动模型是处理大区域流场的一种非常有效的模型。

它属于非线性方程组,在目前只能用数值方法求解,因此,有必要研究一种简单、高效的方法来求解浅水流动问题。

自Patankar 和S palding [1]
发展了SI MP LE 算法以来,该方法被广泛应用于不可压缩流场的数值模拟,而且该方法还得到了进一步的发展,主要有SI MP LER 算法
[2]、SI MP LEC 算法[3]、SI MP LEX 算法[4]和SI M 2P LET 算法[5]等。

这些模型均成功地应用于速度—压力耦合的流场计算,深度平均的浅水流动模型是在静压假定下导出的,一般流体模型中的速度—压力耦合也就转换成浅水流动模型中的速度—水深耦合[6]。

天然河流、海湾的边界曲折、地形复杂,采用坐标变换是解决问题的途径之一。

目前多数N -S 方程的坐标变换中,流程全部采用逆变分量,这样就增加了方程的复杂程度。

于是忽略掉方程中的非正交项,利用正交变换下的方程进行数值求解[7,8]。

对于具有复杂边界的海湾及弯曲的河流,坐标变换中很难保证每个点都正交,特别是边界附近。

水位变化是水力计算中难点之一,在目前的紊流数学模型中,多简单的利用“冻结法”,这样做将失去对边界出流动模拟的准确性。

本文研究中,采用正交曲线坐标变换生成数值网格,而数值计算中采用非正交曲线坐标下的k -ε双方程紊流模型,这样可以自动修正网格生成中的非正交项。

流速除对流项中采用逆变分量,在其余各项中均采用原始分量,这样使得方程书写简单,有利于将各方程写成通用形式,编写的程序变得更规范。

作者受Jian Y e 同位网格[9]
的启发,对普通交错网格做了修改,即采用B 型交错网格,使得u ,v ,k ,ε的计算布置在一个节点上,有利于节省计算程序代码,使程序书写更加规范。

引入动边界扫描技术,结合紊流模型的壁面函数法,使壁面随着真实边界而变化。

数值求解时,采用控制体积法离散方程,运用SI MP LEC 算法,使计算的流场更符合实际流场。

1　数值网格
本文对计算区域用Laplace方程实施坐标变换,生成正交的贴体网格,控制方程:
52x 5ξ2+52x
5η2=0
52y 5ξ2+52y
5η2=0
(1)
方程组(1)可采用S OR方法求解,然后用正交曲线网格边界正交化处理的三次样条插值边界滑动法[10],处理其边界处的网格,以提高其正交性。

与其它方法相比,此方法虽在一定程度上提高了边界处网格的正交性,但难以保证完全正交,如果再把曲线坐标系下N-S方程中的非正交项去掉,由斜交网格引起的计算误差是不容忽视的。

为了修正此项误差,本文采用非正交曲线坐标下的平面二维k-ε双方程紊流模型,模型能自动修正网格生成中的非正交项。

2　数学模型
笛卡儿坐标下深度平均的k-ε双方程紊流模型的通用微分方程:
5
5t(HρΦ)+5
5x(HρuΦ)+
5
5y(HρvΦ)=
5
5x(HΓΦ
5Φ
5x)+
5
5y(HΓΦ
5Φ
5y)+SΦ(2)
方程(2)转换到曲线坐标(ξ、η)下,仅在对流项中使用流速的逆变分量,而在其它项中使用原始变量,这样既简化了方程,又使所有方程仍可写为曲线坐标下的通用方程,模型的微分方程可写为如下通用形式:
5
5t(HρΦ)+1
J
5
5ξ(HρUΦ)+
1
J
5
5η(HρVΦ)=
1
J
5
5ξ(
HΓΦ
J
(αΦξ-βΦη)+
1
J
5
5η(
HΓΦ
J
(-βΦξ+γΦη))+SΦ(ξ,η)(3)
式中:Φ为所求问题的因变量;U和V分别为直角坐标下流速u和v的逆变分量,仅在对流项中出现;ΓΦ为扩散系数;SΦ为源项。

当Φ表示某一特定量时,ΓΦ和SΦ对此特定量有特定的意义和表达式,这时方程(3)亦赋予特定的意义,模型的控制方程组如表1所示。

表1　模型控制方程组
方程ΦΓΦSΦ
连续100
ξ———动量uμ
eff
S u
η———动量νμ
eff
Sν
紊流动能k μ
eff
σ
k
H(G0k+G kν-ρε)
紊流动能耗率εμ
eff
σεH
ε
k
(Cε
1
G0k-Cε
2
ρε)+Gεν
S u=-1
J
ρgΗ(y
η
5z s
5ξ-yξ
5z s
5η)+
1
J
5
5ξ
Hμeff
J
y2η
5u
5ξ-yξyη
5u
5η
+1
J
5
5ξ
Hμeff
J
yξxη
5v
5η-xηyη
5v
5ξ+
1
J
5
5η
Hμeff
J
xξyη
5v
5ξ-xξyξ
5v
5η
+1
J
5
5η
Hμeff
J
y2ξ
5u
5η-yξyη
5u
5ξ-τbx
S V=-1
J
ρgΗ(x
ξ
5z s
5η-xη
5z s
5ξ)+
1
J
5
5ξ
Hμeff
J
x2η
5v
5ξ-xξxη
5v
5η
+1
J
5
5ξ
Hμeff
J
xξyη
5u
5η-xηyη
5u
5ξ+
1
J
5
5η
Hμeff
J
yξxη
5u
5ξ-xξyξ
5u
5η
+1
J
5
5η
Hμeff
J
x2ξ
5v
5η-xξxη
5v
5ξ-τby
U=uyη-vxη, V=vxξ-uyξ,
α=x2
η
+y2η,　β=xξxη+yξyη,　γ=x2ξ+y2ξ,　J=xξyη-yξxη
G0k=2μeff 5u
5ξ
yη
J
-
5u
5η
yξ
J
2
+2μeff-
5v
5ξ
xη
J
+
5v
5η
xξ
J
2
+
μ
eff 5v
5ξ
yη
J
-
5v
5η
yξ
J
+-
5u
5ξ
xη
J
+
5u
5η
xξ
J
2
G kV=c kρU33
H
　Gε
V=
cερU43
H2
　U3=c
f
(u2+v2)　H=z
s
-z b
c k=1
c f
　cε=
316Cε
2
c0175f
Cμ　μeff=ρ(v+v t)　v t=Cμ
k2
ε
τ
bx =gn
2u u2+v2
H1Π3
　τ
by
=gn
2v u2+v2
H1Π3
　c
f
=01003
式中:u和v分别为笛卡儿坐标下x和y方向的深度平均的流速分量;k和ε分别为深度平均的紊动
动能和紊动动能耗散率;H为水深;z
s 为水位;z
b
为河床高程;ν为分子运动粘性系数;ν
t
为漩涡
运动粘性系数;τ
bx 和τ
by
分别为x和y方向上的底摩擦力μ
eff
为有效粘性系数;Cμ、Cε
1
、Cε
2
、σ
k
和σε
为经验常数,其中Cμ为0109,Cε
1为1144,Cε
2
为1192,σ
k
为110,σε为113。

3　模型离散
如图1所示,本文采用B型交错网格,即将u,v,k,ε及各标量布置在网格中心,将水深相关项H、z
s、z b项布置在网格的左下方节点上。

这样布置有以下两个优点:①有利于边界条件的引入;
②在编写程序u,v,k,ε各方程共用一套系数,可节省程序代码,使程序简单明了。

本文用控制容积法,采用混合格式来导出离散方程。

其边界条件为:
在进口边界,所有边界条件都按本征条件给出,即:u=u
,v=0,k=k0,ε=ε0;
在出口边界,给定出口处的水位z
s
,令v=0,其余变量都取法向导数为零,即:
5u 5n=5k
5n=
5ε
5n=0。

图1　B型交错网格布置图2　计算平面内水位变化引起的动边界
在水陆边界,在计算河流、河口、海湾的非恒定流场时,水位不是恒定不变的,而且每一时间步长的水位可能都不一样,这就带来了水陆边界的移动,如果每一时间步长都进行坐标变换而生成新的网格是非常不经济的。

于是本文在前人“冻结法”的基础上,提出了适合k-ε双方程紊流模型的
“移动边界的壁面函数法”。

如图2、图3,水位变化后,使边界处的网格干出,原来的岸边界12、34变成了新的岸边界PQ 、ST ,而壁面函数仍然布置在12、34岸边界处;水域内岛屿由于水位的变化而出露,其岛屿内部按“冻结法”来处理,其边界ABC D 周围也需用壁面函数法来处理。

本文采用“水位扫描法”来判断新的水陆边界:①从边界12向34扫描,如水深不为零,则令其为起始点;②从边界34向12扫描,如水深不为零,则令其为终点;③然后扫描水域内部,确定岛屿的边界ABC D 。

之后对水陆边界实施壁面函数法。

图3　水位变化引起的动边界断面在壁面边界,k -ε双方程紊流模型是高Rey 2
nolds 数模型,只适用于充分发展的紊流,而在近
壁处,粘性效应起主导作用。

另外,由于近壁处
物理量变化很快,需要非常密的网格,才能真实
地模拟实际流动。

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

4　SI MP LEC 算法求解模型
SI MP LE 算法的基本步骤如下:(1)计算坐标变换的相关系数;(2)根据进、出口边界的水位,
给定全场水位,计算全场初始水深H 3;(3)解动量方程求u n ,v n ;(4)解k -ε紊流模型求k 、
ε;(5)解水深校正方程求H ′,并修正水深H =H 3+αh H ′;(6)修正流速u 、v 及有效粘性系数μeff ,然
后重复(2)～(5)步直至收敛。

模型求解中,为了利于非线性迭代的收敛,计算中采用亚松弛技术。

计算中采用ADI 技术和T D 2M A 算法,以连续方程的误差小于给定的值作为判断收敛的依据。

在恒定流计算中,以相邻时间步之间全场各点流速误差的最大绝对值小于给定的值作为达到定常状态的依据。

5　模型应用
为了使本文数学模型具有代表意义,本文采用天然连续弯道河流———美国C olorado 洲Fall River
部
图4　Fall River
河床地形图5　物理平面内数值网格
分河段的实测水流资料进行验证。

该河段长120m ,平均河宽815m ,其河床地形资料如图4所示。

计算流量为411m 3/s ,上游水位2185m ,下游水位2161m 。

用坐标变换的方法将物理平面转换到计算平面。

如图5所示,整个河段可划分为91×20个网格,网格长度平均014m ,宽度平均0109m 。

本模型计算所选用的糙率平均为0162,并在计算中自动调整
[13]。

本文属于恒定流计算,所采用的精度为10
-5,计算经过1500步收敛,得全河段流场图,如图6所示。

为了进一步对本模型进行验证,作者选用具有实测资料的6个断面的流速和水位与计算值进行
图7　各断面计算与实测流速比较
图8　各断面计算与实测水位比较
比较。

由于本河段属于冲积性多沙河流,计算中糙率以及河底高程很难调整到与实际情况完全相符,使得计算中存在一定的误差,而且在一些断面上个别点的计算结果与实测值存在着较大的离散,但总体看来计算结果能够满足工程需要。

流速比较见图7,水位比较见图8。

6　结论
本文建立了非正交曲线坐标下的k-ε双方程紊流模型,可以对具有不规则边界的水域,如海湾、天然连续弯道水流的进行较好的模拟。

模型中采用了一系列针对天然水流的处理技巧:特殊的交错网格、与“水拉扫描法”结合的壁面函数法的动边界处理、保留坐标变换后模型方程中的非正交项等等,均提高了模型的实用性,显示了本模型向水库、河流、海湾等水域的污染预报模型发展的前景。

参　考　文　献:
[1]　Patankar S V,S palding D B.Calculation Procedure for Heat,Mass and M omentum T rans fer in32D Flows[J].Int.
J.Heat Mass T rans fer,1972(15):1787-1806.
[2]　Patankar S V.Numerical heat trans fer and fluid flow[M].New Y ork:McG raw2Hill,1980.
[3]　Van D oormaal J P.Raithby GD.Enhancement of SI MP LE method for predicting incompressible fluid flows[J].Num2
er Heat T rans fer,1984(7):147-163.
[4]　Raithby G D,Schneider G E.E lliptic systems:finite difference methodsⅡ[M].New Y ork:John Wiley&S ons,
1981.
[5]　Sheng Y,Shoukri M,Sheng G,W ood P.A m odification to the SI MP LE method for buoyancy2driven flows[J].Nu2
mer Heat T rans fer,Part B,1998,33(1):65-78.
[6]　Jian G uo Zhou.Velocity2Depth C oupling in Shallow2Water Flows[J].J.Hydr.Engrg.,ASCE,1995,121(10):
717-724.
[7]　王船海,程文辉.河道二维非恒定流场计算方法研究[J].水利学报,1991(1):10-17.
[8]　魏文礼,金忠青.复杂边界河道流速场的数值模拟[J].水利学报,1994,(11):26-30.
[9]　Y e J,McC orquodale J A.Depth2averaged hydrodynamic m odel in curvilinear collocated grid[J].J.Hydr.Engrg.,
ASCE,1997,123(5):380-388.
[10]　吴修广.河道平面二维水流数值模拟与阻力特性研究[D],重庆:重庆交通学院硕士论文,2001.
[11]　Launder B E,S palding D B.The Numerical C omputation of Turbulent Flows[J].C omp.Meth.In Appl.and En2
grg.1974(3):269-289.
[12]　沈永明.油—水两相湍流的数值模拟[D].成都:成都科技大学博士论文,1991.
[13]　吴修广,王平义.山区河流二维阻力研究[J].重庆交通学院学报,2001,20(3):102-105.
(下转第37页)
5:212
[4]　孙双科,等.下游高水位条件下导流洞改建为旋流竖井式泄洪洞的水力学问题研究.水利学报,2000,
(10):22-27
[5]　董兴林,等.高水头大泄量旋涡竖井式泄洪洞的设计研究[J].水利学报,2000,(11):27-33
Changing rule of the encircle flow hollow diameter of horizontal vortex in
an internal energy dissipating tunnel
NI U Zheng2ming1,S UN Jing1,ZH ANG Ming2yuan2
(1.Xi’an University o f Technology,Xi’an　710048,China; 2.Xi’an Jiao Tong University,Xi’an　710048,China)
Abstract:The authors studied and analyzed the changing rule of the encircle flow hollow diameter(d0) of a horizontal v ortex in an internal energy dissipation tunnel with a inlet shaft by analyzing the interrelat2 ed in fluence factors based on experimental data.The results show that the changing rule of d0has an ob2 vious character of zoning.The main interrelated in fluencing factors are different in different flow zones.
A general function of synthetically dimensionless variables and em pirical formula were established accord2
ing to the experimental results,which will be a g ood reference in hydraulic characteristic study,analysis and design of this type of tunnel.
K ey w ords:Inlet shaft;encircle flow hollow diameter;internal energy dissipation tunnel;horizontal v ortex;dimension analysis
(上接第30页)
22D flow SIMP LEC algorithm in non2orthogonal curvilinear coordinates
W U X iu2guang1,SHE N Y ong2ming1,ZHE NG Y ong2hong2,W ANG Ping2yi2
(1.State K ey Laboratory o f Coastal and Offshore Engineering,Dalian University o f Technology,Dalian　116024,China;
21Guangzhou Institute o f Energy Conversion,Chinese Academy o f Science,Guangzhou　510070,China;
3.Department o f River and Ocean Engineering,Chongqing Jiaotong University,Chongqing　400074,China.)
Abstract:In this paper,the Laplace equation was adopted in coodrdinate trans formation to generate or2 thog onal curvilinear grid.The shallow2water equations were trans formed under the curvilinear coordinate system.During discretization,B type staggered grid was introduced to sim plify the program.Water lev2 el scanning method combined with wall2function was used to treat the m oving boundary.The turbulence m odel of double k2εequations in non2orthog onal curvilinear coordinate was s olved through SI MP LEC alg o2 rithm,which m odified the error derived from the non2orthog onal grid.The m odel was validated by the measured data in the Fall River in C olorado State,US.The calculation results were in g ood accordance with the measurement.
K ey w ords:coordinate trans formation;k-εturbulence m odel;water level scanning method;wall2 function;SI MP LEC alg orithm。