默认头像
gis
管理员
管理员
  • 注册日期2003-07-16
  • 发帖数15951
  • QQ
  • 铜币25345枚
  • 威望15368点
  • 贡献值0点
  • 银元0个
  • GIS帝国居民
  • 帝国沙发管家
  • GIS帝国明星
  • GIS帝国铁杆
阅读:3600回复:2

求解守恒形式的二维浅水方程的SIMPLE类程式

楼主#
更多 发布于:2005-05-31 15:21

韩龙喜

河海大学水文水资源及环境学院

摘 要:对二维浅水方程守恒形式,采用控制体积法离散,可得流速与水深的耦合关系.据此,本文在SIMPLE类程式的基础上,提出一种新的、具有更快收敛速度的流速校正公式.作为测试算例,对由两种不同底坡组成的棱柱形明渠水流进行了数值模拟,取得了预期的效果.

关键词:二维浅水方程,守恒形式,SIMPLE类程式,流速校正

作者简介:韩龙喜(1964-),男,江苏邗江人,河海大学副教授,博士.

1 SIMPLE类程式概述

  自SIMPLE算法[1]问世以来,该方法被广泛应用于不可压流体的数值模拟.在此基础上,相继出现了SIMPLER程式[2]及SIMPLEC程式[3]等改进形式.上述程式,当速度场与压力场成线性关系(如有压流)时,流速校正公式因所有的系数独立于压力,而精确反映了流速校正与压力校正的耦合关系.但在求解明渠水流时,由于流速场与水深场成平方关系而非线性关系,现行的流速校正公式对流速与水深的耦合关系产生了近似.

  本文采用控制体积法,由守恒形式的二维浅水方程导出差分方程,得到流速与水深的耦合关系.据此关系对SIMPLE类程式的流速校正公式进行改进.

2 守恒形式的二维浅水运动方程

  忽略哥氏力及风成应力,二维非恒定浅水运动方程的守恒形式为

ht+(uh)x+(vh)y=0

(1)

(uh)t+(u·uh)x+(u·vh)y=-g(h2/2)x-ghzx-gn2+(hεux)x+(hεuy)y

(2)

uh)t+(u·uh)x+(u·vh)y=-g(h2/2)x-ghzx-gn2+(hεux)x+(hεuy)y

(3)

式中:h为水深,z为水位,u、v为x、y方向的流速,n为糙率,ε为紊动扩散系数.

3 方程的求解

3.1 方程的离散  采用控制体积法离散方程,h,u,v分布采用交错网格技术,假设水深在相邻节点间呈线性变化,得h,u,v差分方程:

[(hP-h0P)/Δt]ΔxΔy+(hp+hE)ueΔy-(hP+hW)uwΔy+(hP+hN)unΔx-(hP+hS)usΔx=0

(4)

aPuP=aEuE+aWuW+aNuN+aSuS+a0

(5)

式中

式中:为v在u的控制体内的平均值.同理可得v的离散方程

bPvP=bEvE+bWvW+bNvN+bSvS+b0

(6)

式中:bE,bW,bN,bS与a相似,

b0=b0Pv0P+g/2△x(h2s-h2n)+g/2△x(hs-hn)(zs-zn)

bP=bE+bW+bN+bS+b0P-SP(v)△x△y

(7)

3.2 流速校正公式

3.2.1 流速与水深的耦合   采用欠松弛技术,引入松弛系数a(0<α<1),式(5),式(6)用下式表示

aPuP=α(aEuE+aWuW+aNuN+aSuS+a0)+(1-α)aPP

bPuP=α(bEuE+bWuW+bNuN+bSuS+a0)+(1-α)bPP

P,P表示u及v在前一次迭代值.以上二式可进一步写成

aPuP=aΣanbunb+ac+ah(h2w-h2e)+az(hw+he)

(8)

ac=(1-a)aPP+aa0pu0P

(9)

ah=a(g/2)△y

(10)

az=a(g/2)△y(zw+ze)

(11)

同理可得

bPuP=aΣbnbvnb+bc+bh(h2w-h2e)+bz(hs+hn)

(12)

bc=(1-a)bPP+ab0pv0P

(13)

bh=a(g/2)△y

(14)

bz=a(g/2)△y(zs-zn)

(15)

由式(8),式(12)可知,流速分量(u,v)与水深平方及水深本身成正比,本文拟从此关系出发,导出SIMPLE类计算程式速度校正公式的改进形式.

3.2.2 流速校正公式的改进 给定水深场的猜测值h*,由下式求得相应的流速场

(16)

(17)

式中:上标星号代表此值由h*求得.通常h*,u*,v*不满足连续方程,因此需对水深h*不断进行改进.设正确的水深场h与猜测的水深场h*之差为h′,则有

h=h*+h′

(18)

h′为水深校正.正确的速度场u,v与猜测的u*,v*之差为u′,v′,则有

u=u*+u′, v=v*+v′

(19)

假设aP=a*P,anb=a*nb,将方程(8)减去方程(16),(12)减去(17)得

(20)

显然,上式精确反映了流速校正与水深校正之间的本构关系.

  若采用SIMPLE程式,引入松弛系数,通常采用以下流速校正公式[4]

uP=u*P+du(h′w-h′e)

(21)

vP=v*P+dv(h′s-h′n)

(22)

du=ah(hw+he)/aP       dv=bh(hs+hn)/bP

  显然,式(21)、(22)对速度校正与水深校正的关系式(20)作了近似.本文以SIMPLE程式为例,导出反映流速-水深耦合关系的流速校正公式的改进形式.

  以u方程为例,将式(18)代入式(20)得

aPu′P=ah(2h*wh′w+h'2w-2h*ehe-h2e)+az(h′w+h′e)

(23)

忽略二阶小量水深校正的平方,得流速u的校正方程

uP=u*P+duwh′w-dueh′e

(24)

(25)

due表示:在u的控制体内,u的校正随边界e上的水深校正的变化关系,duw类同.同理,可导得流速v的校正公式

vP=v*P+dvwh′s-dvnh′n

(26)

(27)

  从物理机制的角度出发,分析流速校正公式系数d的表达式可知,水深较正通过两方面 影响流速校正:①ah(bh)——水深校正引起水深变化,进一步引起水压力的变化,从而引起动量的变化;②az(bz)——水深校正引起控制体内水体体积(重力)的变化,由于某方向渠底底坡的存在,引起底坡对控制体反作用力在该方向分量的变化,而最终引起动量的变化.特殊地,底坡为平坡时,az(bz)=0,此影响为0.若忽略②的影响,并以u(或v)控制 体中心的水深(h*w+h*e)/2(或h*s+h*n)/2分别代替控制体界面上的h*w,h*e(或h*s,h*n),则方程(24)及(26)退化成方程(21)~(22).

  从上述分析可知,改进的流速校正公式,从流速与水深的耦合关系出发,更精确地反映了水流运动的物理机制,比常用的流速校正公式更趋合理,具有更高的效率.用同样的方法可得SIMPLEC程式的流速校正公式.

3.2.3 水深校正方程 采用改进的流速校正公式,把式(24)及式(26)代入 连续方程,得水深校正方程

CPh′p=CEh′E+CWh′W+CNh′N+CSh′S

(28)

式中

(29)

(30)

(31)

式中:d(e)uE表示在水深hP的控制体交界面e上的流速u的校正值随E点水深校正的变化关系,其余类同.比较方程(29)与(30)可知,系数CP并不严格等于相邻系数之和,但分析表达式可知,其间差异源于d的表达式的不同,忽略此细微的差异,令CP等于相邻系数之和,可保证高斯-赛德尔迭代法求解水深校正方程的收敛.

4 算例

  图1所示明渠由两段不同的底坡组成,负坡底坡i1=-0.01;顺坡底坡i2=0.02.两段明渠长同为28m,宽度同为1.12m,渠底糙率取n=0.02.上边界采用流量边界条件Q=0.2m3/s.由水力学公式求得,临界水深为0.1571m,顺坡段正常水深为0.1106m.水深在负坡上发生A2型水面曲线;在顺坡下游为均匀急流;在顺坡初始段水流从缓流向急流过渡,发生S2型水面曲线,此水流包含了缓流、急流及结流向急流的过渡流态,从上游至下游水深变化急剧.因此对此水流的模拟计算,可作为本模型测试的合适算例.

  对SIMPLEC程式采用本文提出的流速校正公式计算,Δx=0.56,Δy=0.14,取质量源 残差R=8×10-4为收敛判据.计算时,流速松弛系数取0.7;水深松弛系数,在前两 百次迭代取0.3,后取0.9.经过1689次迭代得到收敛解,在顺坡与负坡的交界面上,计算水深为0.1549m,接近临界水深值;在顺坡下游出口断面附近计算水深为0.1112m,流速为1.629m/s,入口断面水深为0.5052m.水深沿渠道纵向变化见图1,流速矢量分布见图2.


图1 渠底形状及计算水深示意

图2 流速矢量

  本文还采用常用的流速校正公式进行计算,经过1786次迭代得到收敛解,结果与前算例相近.两种不同方法质量源残差随迭代次数的变化对比见图3.?

5 结论

  (1)从二维浅水方程的守恒形式导得的差分方程,较从非守恒形式控制方程导出的差分方程,更精确地反映了水流运动遵循的物理机制.?


图3 质量源残差变化曲线

  (2)本文从守恒形式导出的流速校正公式,更准确地反映了流速校正对水位校正的依赖关系,算例表明:当水深变化较大时,具有更快的收敛速度.

参 考 文 献

[1] Patankar S V and Spalding D B.A calculation proced ure for heat,mass and momentum transfer in three-dimensional flows[J].Int.J.H eat Mass Transfer,1972,(15):1787~1806.

[2] Patankar S V.A calcuation procedure for two-dimensional e lliptic situation[J].Numerical Heat Transfer,1981,(4):409~425.

[3] Van Doormaal J P and Raithby G D.Enhancement of the SIMPL E method for prediction incompressible fluid flows[J].Numerical Heat Transfer, 1984,7(2):147~163.

[4] Jian Guo Zhou.Velocity-Depth Coupling in Shallow-Water F lows[J].J.Hydr Engrg.Div.,ASCE,1995.718-723.

喜欢0 评分0
GIS麦田守望者,期待与您交流。
默认头像
路人甲
路人甲
  • 注册日期2005-09-26
  • 发帖数2
  • QQ
  • 铜币104枚
  • 威望0点
  • 贡献值0点
  • 银元0个
1楼#
发布于:2005-09-26 18:54

请问版主有这篇文章的程序吗?xinpingy@hotmail.com.

举报 回复(0) 喜欢(0)     评分
默认头像
路人甲
路人甲
  • 注册日期2005-09-26
  • 发帖数2
  • QQ
  • 铜币104枚
  • 威望0点
  • 贡献值0点
  • 银元0个
2楼#
发布于:2005-09-26 18:52

举报 回复(0) 喜欢(0)     评分
默认头像

返回顶部