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

楼主^#

更多发布于：2005-05-31 15:21

韩龙喜
河海大学水文水资源及环境学院
摘要：对二维浅水方程守恒形式，采用控制体积法离散，可得流速与水深的耦合关系.据此，本文在SIMPLE类程式的基础上，提出一种新的、具有更快收敛速度的流速校正公式.作为测试算例，对由两种不同底坡组成的棱柱形明渠水流进行了数值模拟，取得了预期的效果.
关键词：二维浅水方程，守恒形式，SIMPLE类程式，流速校正
作者简介：韩龙喜(1964-)，男，江苏邗江人，河海大学副教授，博士.
1 SIMPLE类程式概述
　　自SIMPLE算法[1]问世以来，该方法被广泛应用于不可压流体的数值模拟.在此基础上，相继出现了SIMPLER程式[2]及SIMPLEC程式[3]等改进形式.上述程式，当速度场与压力场成线性关系(如有压流)时，流速校正公式因所有的系数独立于压力，而精确反映了流速校正与压力校正的耦合关系.但在求解明渠水流时，由于流速场与水深场成平方关系而非线性关系，现行的流速校正公式对流速与水深的耦合关系产生了近似.
　　本文采用控制体积法，由守恒形式的二维浅水方程导出差分方程，得到流速与水深的耦合关系.据此关系对SIMPLE类程式的流速校正公式进行改进.
2 守恒形式的二维浅水运动方程
　　忽略哥氏力及风成应力，二维非恒定浅水运动方程的守恒形式为
<TABLE cellSpacing=0 cellPadding=0 width="100%" border=0>

<TR>
<TD width="81%">
ht+(uh)x+(vh)y=0</TD>
<TD width="19%">
(1)</TD></TR>
<TR>
<TD align=middle width="81%">
(uh)t+(u·uh)x+(u·vh)y=-g(h2/2)x-ghzx-gn2<IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/1602.jpg" align=middle border=0>+(hεux)x+(hεuy)y</TD>
<TD align=right width="19%">
(2)</TD></TR>
<TR>
<TD align=middle width="81%">
uh)t+(u·uh)x+(u·vh)y=-g(h2/2)x-ghzx-gn2<IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/1603.jpg" align=middle border=0>+(hεux)x+(hεuy)y</TD>
<TD align=right width="19%">
(3)</TD></TR></TABLE>
式中：h为水深，z为水位，u、v为x、y方向的流速，n为糙率，ε为紊动扩散系数.
3 方程的求解
3.1 方程的离散 采用控制体积法离散方程，h,u,v分布采用交错网格技术，假设水深在相邻节点间呈线性变化，得h,u,v差分方程：
<TABLE cellSpacing=0 cellPadding=0 width="100%" border=0>

<TR>
<TD width="82%">
[(hP-h0P)/Δt]ΔxΔy+(hp+hE)ueΔy-(hP+hW)uwΔy+(hP+hN)unΔx-(hP+hS)usΔx=0</TD>
<TD align=right width="18%">
(4)</TD></TR>
<TR>
<TD width="82%">
aPuP=aEuE+aWuW+aNuN+aSuS+a0</TD>
<TD align=right width="18%">
(5)</TD></TR></TABLE>
式中
<IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/16-1.jpg">
式中:<IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/16-2.jpg">为v在u的控制体内的平均值.同理可得v的离散方程
<TABLE cellSpacing=0 cellPadding=0 width="100%" border=0>

<TR>
<TD width="85%">
bPvP=bEvE+bWvW+bNvN+bSvS+b0</TD>
<TD width="15%">
(6)</TD></TR></TABLE>
式中:bE,bW,bN,bS与a相似，
<TABLE cellSpacing=0 cellPadding=0 width="100%" border=0>

<TR>
<TD width="87%">
b0=b0Pv0P+g/2△x(h2s-h2n)+g/2△x(hs-hn)(zs-zn)
bP=bE+bW+bN+bS+b0P-SP(v)△x△y</TD>
<TD width="13%">
(7)</TD></TR></TABLE>
3.2 流速校正公式
3.2.1 流速与水深的耦合 采用欠松弛技术，引入松弛系数a(0<α<1)，式(5)，式(6)用下式表示
aPuP=α(aEuE+aWuW+aNuN+aSuS+a0)+(1-α)aP<IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/16-3.jpg" align=middle border=0>P
bPuP=α(bEuE+bWuW+bNuN+bSuS+a0)+(1-α)bP<IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/16-4.jpg" align=middle border=0>P
<IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/16-3.jpg" align=middle border=0>P,<IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/16-4.jpg" align=middle border=0>P表示u及v在前一次迭代值.以上二式可进一步写成
<TABLE cellSpacing=0 cellPadding=0 width="100%" border=0>

<TR>
<TD width="82%">
aPuP=aΣanbunb+ac+ah(h2w-h2e)+az(hw+he)</TD>
<TD align=right width="18%">
(8)</TD></TR>
<TR>
<TD width="82%">
ac=(1-a)aP<IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/16-3.jpg" align=middle border=0>P+aa0pu0P</TD>
<TD align=right width="18%">
(9)</TD></TR>
<TR>
<TD width="82%">
ah=a(g/2)△y</TD>
<TD align=right width="18%">
(10)</TD></TR>
<TR>
<TD width="82%">
az=a(g/2)△y(zw+ze)</TD>
<TD align=right width="18%">
(11)</TD></TR></TABLE>
同理可得
<TABLE cellSpacing=0 cellPadding=0 width="100%" border=0>

<TR>
<TD width="84%">
bPuP=aΣbnbvnb+bc+bh(h2w-h2e)+bz(hs+hn)</TD>
<TD align=right width="16%">
(12)</TD></TR>
<TR>
<TD width="84%">
bc=(1-a)bP<IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/16-3.jpg" align=middle border=0>P+ab0pv0P</TD>
<TD align=right width="16%">
(13)</TD></TR>
<TR>
<TD width="84%">
bh=a(g/2)△y</TD>
<TD align=right width="16%">
(14)</TD></TR>
<TR>
<TD width="84%">
bz=a(g/2)△y(zs-zn)</TD>
<TD align=right width="16%">
(15)</TD></TR></TABLE>
由式(8)，式(12)可知，流速分量(u,v)与水深平方及水深本身成正比，本文拟从此关系出发，导出SIMPLE类计算程式速度校正公式的改进形式.
3.2.2 流速校正公式的改进 给定水深场的猜测值h*，由下式求得相应的流速场
<TABLE cellSpacing=0 cellPadding=0 width="100%" border=0>

<TR>
<TD width="86%">
<IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/1616.jpg"></TD>
<TD align=right width="14%">
(16)</TD></TR>
<TR>
<TD width="86%">
<IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/1617.jpg"></TD>
<TD align=right width="14%">
(17)</TD></TR></TABLE>
式中:上标星号代表此值由h*求得.通常h*,u*,v*不满足连续方程，因此需对水深h*不断进行改进.设正确的水深场h与猜测的水深场h*之差为h′,则有
<TABLE cellSpacing=0 cellPadding=0 width="100%" border=0>

<TR>
<TD width="83%">
h=h*+h′</TD>
<TD width="17%">
(18)</TD></TR></TABLE>
h′为水深校正.正确的速度场u,v与猜测的u*,v*之差为u′,v′，则有
<TABLE cellSpacing=0 cellPadding=0 width="100%" border=0>

<TR>
<TD width="85%">
u=u*+u′, v=v*+v′</TD>
<TD width="15%">
(19)</TD></TR></TABLE>
假设aP=a*P,anb=a*nb,将方程(8)减去方程(16)，(12)减去(17)得
<TABLE cellSpacing=0 cellPadding=0 width="100%" border=0>

<TR>
<TD width="92%">
<IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/1620.jpg"></TD>
<TD width="8%">
(20)</TD></TR></TABLE>
显然，上式精确反映了流速校正与水深校正之间的本构关系.
　　若采用SIMPLE程式，引入松弛系数，通常采用以下流速校正公式[4]
<TABLE cellSpacing=0 cellPadding=0 width="100%" border=0>

<TR>
<TD width="86%">
uP=u*P+du(h′w-h′e)</TD>
<TD align=right width="14%">
(21)</TD></TR>
<TR>
<TD width="86%">
vP=v*P+dv(h′s-h′n)</TD>
<TD align=right width="14%">
(22)</TD></TR></TABLE>
du=ah(hw+he)/aP dv=bh(hs+hn)/bP
　　显然，式(21)、(22)对速度校正与水深校正的关系式(20)作了近似.本文以SIMPLE程式为例，导出反映流速-水深耦合关系的流速校正公式的改进形式.
　　以u方程为例，将式(18)代入式(20)得
<TABLE cellSpacing=0 cellPadding=0 width="100%" border=0>

<TR>
<TD width="85%">
aPu′P=ah(2h*wh′w+h'2w-2h*ehe-h2e)+az(h′w+h′e)</TD>
<TD width="15%">
(23)</TD></TR></TABLE>
忽略二阶小量水深校正的平方，得流速u的校正方程
<TABLE height=35 cellSpacing=0 cellPadding=0 width="100%" border=0>

<TR>
<TD width="85%" height=17>
uP=u*P+duwh′w-dueh′e</TD>
<TD align=right width="15%" height=17>
(24)</TD></TR>
<TR>
<TD width="85%" height=18>
<IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/1625.jpg"></TD>
<TD align=right width="15%" height=18>
(25)</TD></TR></TABLE>
due表示：在u的控制体内，u的校正随边界e上的水深校正的变化关系，duw类同.同理，可导得流速v的校正公式
<TABLE cellSpacing=0 cellPadding=0 width="100%" border=0>

<TR>
<TD width="84%">
vP=v*P+dvwh′s-dvnh′n</TD>
<TD align=right width="16%">
(26)</TD></TR>
<TR>
<TD width="84%">
<IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/1627.jpg"></TD>
<TD align=right width="16%">
(27)</TD></TR></TABLE>
　　从物理机制的角度出发，分析流速校正公式系数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)代入连续方程，得水深校正方程
<TABLE cellSpacing=0 cellPadding=0 width="100%" border=0>

<TR>
<TD width="87%">
CPh′p=CEh′E+CWh′W+CNh′N+CSh′S</TD>
<TD width="13%">
(28)</TD></TR></TABLE>
式中
<TABLE cellSpacing=0 cellPadding=0 width="100%" border=0>

<TR>
<TD width="86%">
<IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/1629.jpg"></TD>
<TD align=right width="14%">
(29)</TD></TR>
<TR>
<TD width="86%">
<IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/1630.jpg"></TD>
<TD align=right width="14%">
(30)</TD></TR>
<TR>
<TD width="86%">
<IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/1631.jpg"></TD>
<TD align=right width="14%">
(31)</TD></TR></TABLE>
式中: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.
<TABLE cellSpacing=0 cellPadding=0 width="100%" border=0>

<TR>
<TD align=middle width="50%"><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/1601.gif"> 图1 渠底形状及计算水深示意</TD>
<TD align=middle width="50%"><IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/1602.gif"> 图2 流速矢量</TD></TR></TABLE>
<TABLE cellSpacing=0 cellPadding=0 width="100%" border=0>

<TR>
<TD width="67%">
　　本文还采用常用的流速校正公式进行计算，经过1786次迭代得到收敛解，结果与前算例相近.两种不同方法质量源残差随迭代次数的变化对比见图3.?
5 结论
　　(1)从二维浅水方程的守恒形式导得的差分方程，较从非守恒形式控制方程导出的差分方程，更精确地反映了水流运动遵循的物理机制.?</TD>
<TD width="33%">
<IMG src="http://www.cws.net.cn/Journal/slxb/200007/images/1603.gif"> 图3 质量源残差变化曲线</TD></TR></TABLE>
　　(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

举报回复

xmxinping 路人甲注册日期2005-09-26 发帖数2 QQ 铜币104枚威望0点贡献值0点银元0个加关注写私信	1楼^# 发布于：2005-09-26 18:52 <P><img src="images/post/smile/dvbbs/em06.gif" /></P>
	举报回复(0) 喜欢(0) 评分

xmxinping 路人甲注册日期2005-09-26 发帖数2 QQ 铜币104枚威望0点贡献值0点银元0个加关注写私信	2楼^# 发布于：2005-09-26 18:54 <P>请问版主有这篇文章的程序吗？xinpingy@hotmail.com.</P>
	举报回复(0) 喜欢(0) 评分

发帖回复

« 返回列表

您需要登录后才可以回帖，登录或者注册

返回顶部

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

最新喜欢：