算法网格计算

收稿日期:作者简介:杨彬(%$'!
—),男,江苏淮安人,硕士,主要从事计算流体力学研究(非结构网格上浅水方程的!
"#$%$隐式算法杨彬,汪德(河海大学环境科学与工程学院,江苏南京!
%""$')摘要:针对浅水数值模拟的特点,建立了可以在非结构网格上求解浅水方程的)*+,-,隐式算法(对复杂地形的底坡项进行了修正,并利用淹没节点法处理移动边界问题(利用算例对此隐式算法进行了模拟验证,并与.
/0显式算法进行了比较(结果表明,隐式算法可以有效地缩短计算时间,具有良好的计算和收敛效果(长江口实际潮流模拟结果也证明了此隐式算法处理实际浅水问题的能力(关键词:非结构网格;浅水方程;)*+,-,隐式算法;复杂地形;移动边界中图分类号:12%&%文献标识码:3文章编号:%"""!
%$'"(!
""')"4!
"4'&!
"5在计算浅水动力学领域中,有结构网格有时难以准确描述计算域边界和复杂地形,而非结构网格能够解决复杂边界问题(但是,在非结构网格上采用隐式算法的确比较困难,因为很多快速收敛的高效算法都基于有结构网格(比较典型的隐式算法如交替方向隐格式(367)法及近似因子分解法(38)等,较适用于有结构的矩形或曲线网格[%](由于非结构网格舍去了网格节点的结构性限制,很多结构网格计算中的高效算法都无法适用(在非结构网格的计算中,一些常用的隐式算法往往要进行复杂的矩阵计算,需要大量的存储空间(这对于大型的数值计算问题来说,可能无法完成计算(由9:;0/0@+ABB0@)隐式算法已成功地应用于非结构网格[&!
5]隐式算法通过对隐式算子的近似因式分解,完全消除了对矩阵的存储,从而使构造的),*算子具有最大程度的对角占优,解的收敛速度和程度有了很大提高(本文首先建立了可以在非结构网格上求解浅水方程的)*+,-,隐式算法,然后利用算例对此)*+,-,隐式算法和.
/0显式算法进行了模拟验证,最后通过对长江口潮位过程的模拟证明了本文所建)*+,-,隐式算法在非结构网格上求解浅水方程的可行性(&控制方程计算模型采用二维浅水方程组!
!
!
!
"!
"!
#"!
#!
$%$(%)其中!
%(&,&',&()1"%(&',&'!
")&!
*!
,&'()1#%(&',&'(,&(!
")&!
*!
)1$%(",)&+"#,)&+"$)1H#,,)&+H$)1式中:+"#,+"$———#,$方向的底坡;+H#,+H$———#,$方向的摩阻坡度;',(———#,$方向的水深平均流速;!
———时间;&———水深;)———重力加速度-分解将空间半离散化的方程写成隐式格式:.
/!
(0"%)/,!
(0)/"!
""12%%%,%3(#!
(0"%))42%.
/$(!
(0)/)(!
)其中#!
%(&,&5,&6)1式中:%I%———欧拉旋转矩阵的逆矩阵;5J'F//HM/N:E*=EO0@,FE0=F0(&LP/(49A>(!
""'轴的夹角(逆时针)!
方程(!
)左边第!
项为""#时间步净输出通量之和!
隐式算法是将式(!
)的隐式部分线性化后,做近似因式分解处理的!
具体方法如下:在方程(!
)的两边同时减去"次计算净输出通量之和:#$!
!
(")$!
%&!
'()#"*#+("!
("&#))*"*#+("!
(")[]'()#"*#+("!
(")),(())#$!
!
(")$!
%&!
'()#"*#+("!
("&#))*"*#+("!
(")[](")$(,)其中!
!
(")$)!
("&#)$*!
(")$*+-(")'()#"*#+("!
(")),(称为"次计算的残值项!
式())右边第!
项一阶数值通量的求解仿照*-+的黎曼间断问题的解法,并对间断处进行熵修正[.
]!
+("!
))+("!
$,"!
())#!
+("!
$)&+("!
()*!
).
)##!
.
"(##.
)$$[].
(/)式中:##.
,$$.
———%$的特征值、对应右特征向量;.
———特征值个数;"———熵修正函数!
对方程(,)左边第!
项采用一阶近似,则[0]"*#+("!
("&#))*"*#+("!
(")))"*#+("!
("&#)$,"!
("&#)()*"*#+("!
(")$,"!
(")())"*#+("!
("&#)$,"!
("&#)()*"*#+("!
(")$,"!
("&#)([])&"*#+("!
(")$,"!
("&#)()*"*#+("!
(")$,"!
(")([])(.
)"*#+("!
("&#))*"*#+("!
(")))"*#&+("!
)&"!
$!
"!
(")("!
)&"!
(!
"!
(")((0)$%&'('因式分解时,一阶数值通量写成如下形式:+("!
))+("!
$,"!
())#!
+("!
$)&+("!
()*%$(("!
(*"!
$[])(1)式中%$(为通量234-5637平均矩阵的谱半径[1]!
将式(1)代入式(0),得"*#+("!
("&#))*"*#+("!
(")))"*##!
&+("!
$)&"!
$&%$([]%!
"!
(")!
&+("!
()&"!
(*%$([]%!
"!
(")((8)将式(8)代入式(,),整理后得#$!
%&#!
!
'()#%$(,()(!
!
(")$)*#!
!
'()#"*#&+("!
()&"!
(*%$([]%,(!
"!
(")(&*+-(")$(#9)因每一个控制体都是闭合的,所以#!
!
'()#"*#&+("!
$)&"!
$,()9!
令/)#$!
%&#!
!
'()#%$(,()(,,与0分别对应计算单元编号小和大的邻居单元,那么经过$%&'('因式分解的完整迭代形式如下:向前扫描(即按编号从小到大进行计算)!
!
'$)&*#*+-$*#!
!
((,"*#&+("!
()"!
(*%$([]%,(!
"!
'{}((##)向后扫描(即按编号从大到小进行计算)$*#!
&*#!
((0"*#&+("!
()"!
(*%$([]%,(!
"!
((#!
)最后,更新数值解,得第""#时间层上的守恒变量值!
("&#)$)!
(")$&!
!
$!
地形和动边界的处理在浅水模型中,地形处理效果直接影响到数学模型的稳定性和合理性!
本模型地形处理时采用了斜底模型的概念,地面高程布置在各个节点上,控制体单元各边中点和单元中心的高程通过线性插值得到!
因为本模型采用三角形单元,所以节点都在同一平面上!
利用单元内部流量守恒确定单元各边法向流速和切向流,1,河海大学学报(自然科学版)第).
卷速!
对于底坡项的修正,将方程(!
)中底坡项的积分表示成!
!
"#$"$$!
%&!
'($%&'","#$"$$!
%&!
'($'()()"*并将其纳入源项后进行显式处理[+]!
移动边界的处理是浅水模拟中的难题!
本文采用淹没节点法[$"],即根据单元各节点的淹没情况对单元和各界面的平均水深进行修正!
采用修正后的平均水深作为单元是否露滩的标准,采用修正后的界面平均水深进行界面数值通量的计算,并通过求解修正单元平均水深达到间接求解实际平均水深的目的!
!
算例为了验证,-.
/0/隐式算法在收敛速度方面的优越性以及本模型所采用的地形处理方法的正确性,利用图"计算水位0,-12-3-24,-.
/0/隐式算法和1&2显式算法对缓、急流过渡流动算例进行模拟计算,并对二者的残值收敛历史及收敛时的计算机总耗时进行比较3由于显式算法的时间步长受到计算稳定性要求的限制,显式算法的时间步长取满足稳定性要求的时间步长的最大值,而,-.
/0/隐式算法在全场可以采用比较大的时间步长!
设矩形水槽的长为!
45、宽为$"5,地形条件为)6"578{)5*#(+*+%)!
,"}其中)5""9!
5+%"$"5#""9"40&缓流!
边界条件:上游流量条件5@ABA49;-)-格式!
0(1)总计算时步2次稳定步2次计算耗时2'缓流激波过渡流动缓流激波过渡流动缓流激波过渡流动缓流激波过渡流动1&2显,-.
/0/隐"9"$"9"!
4949>"9$"9$94!
4@9;"9;"9;9]中的实验数据相当吻合3从显、隐式算法的收敛曲线(图!
)可以看出,,-.
/0/隐式算法的收敛速度明显优于1&2显式算法,而且随着时间步长的加大,收敛速度也在提高3但从收敛曲线4>;第;期杨彬,等非结构网格上浅水方程的,-.
/0/隐式算法和计算耗时的比较结果可以看出,随着时间步长的加大,残值收敛曲线不再有明显变化,计算耗时也不再明显减少!
这说明,对"#$%&%隐式算法而言,时间步长也不是越大越好,时间步长也要根据模拟精度要求和计算耗时来选取!
!
长江口潮位验证长江口是一个中等强度的潮汐河口,全年均为涨、落潮双向往复流!
通过对长江口水流现象的模拟,可以很好地检验数学模型的性能!
本文利用"#$%&%模型模拟长江口'(('!
('!
')*+,:((至'(('!
(-!
(-*+,:((共,'.
的潮流过程!
开边界均采用水位边界条件,全场采用'/的初始水位进行迭代计算!
采用三角形网格,计算网格单元为)0()个,节点为011+个!
最大空间步长约为+(2/,最小空间步长为-((/,时间步长为'0(3!
'(('!
(-!
('*++:((口门大潮涨急和'(('!
(-!
(-*1:((口门大潮落急'个典型瞬时潮流场模拟结果如图-所示,'(('!
(-!
(+*(:((至'(('!
(-!
('*'0:((共计0).
的高桥二站和连兴港两观测站的潮位验证结果如图0所示!
图"计算表层潮流场01$.
*(21+)*+,*)0$2/$/03.
4,-1)+.
*)5.
-图6潮位验证结果#$%&6784.
-9.
()/(+2:;,0.
(0$()**.
9.
*4从模拟计算结果可以看出,相对于456显式算法,"#$%&%隐式算法虽然每时间步存在迭代计算,但时间步长可以取456显式算法的+(倍以上,所以总的模拟时间大大减少!
两观测站的潮位验证结果表明,计算潮位与实测潮位吻合较好[+(]!
这说明本文所建立的"#$%&%隐式算法能够较好地应用于含复杂地形和动边界的实际工程问题!
;>G:HFI&DH>$E"%=>$J>%E>I"$#%DH$#EH#$>J$:D"I&$K:H$&$LH"%H$>&'J&D3&=:>$BCH"Q>DE";;>[!
].
7M77!
"#$%&',()),(():N)(!
N)*.
[-]王波兰.
二维非结构网格2#'>$方程的@ABC8C算法[R].
南京:南京航空航天大学,+,,N.
[*]王立强.
基于动态非结构重叠网格的有限体积法研究[R].
南京:南京航空航天大学,+,,O.
[S]于普兵.
二维浅水水流数值模拟技术研究[R].
南京:南京水利科学研究院,+,,-.
[)]朱德军.
处理二维浅水流动中动边界问题的淹没节点法[R].
北京:清华大学,+,,-.
[(,]周杰.
河口黏性泥沙输运数值模拟:以长江口为例[R].
南京:河海大学,+,,O.
00.
1)-)23.
14,5',*631.
14,5,535&1-3'13-)78-47&9:;%;>I"$HG"B:;>%D:"%&'DF&''"GG&H>$>T#&H:"%E";J$:DG&DD>%H>.
7D'"HF"G&D>;I"$H"&H;>%H.
UF>D#K;>$J>B%">;>HF"(C3P)G&D#D>H">&'G:HFHF>;DE"%E>$%:%JHF>;"=:%JK"#%&$L.
UG">V&;DG>$>E&'E#'&H>&%HF>$>D#'HDG>$>E";G:HFHF"D>E";KLHF>>VDEF>;>.
MH:DDF"G%HF&HHF:DDEF>;>J$>&H'L$>#E>DE&'E#'&H:"%H:;>&%F&DJ"">II:E:>%EL&%E"%=>$J>%E>.
UF>$>D#'HD"I&%#;>$:E&'D:;#'&H:"%"IHF>W&%JHX>4:=>$2DH#&$LI#$HF>$>;"%DH$&H>HF&HHF>DEF>;>:D&EE#$&H>&%>II:E:>%H.
@)A0,-7&:#%DH$#EH#$>J$:;DF&''"GG&H>$>T#&H:"%D;:;;>;E";VH"$'&K.
E%Y&E&>;:EYDF"GP>>H:%J.
&DZNN[:Z(-/-.
(本刊编辑部供稿)*S/第/期杨彬,等非结构网格上浅水方程的@ABC8C隐式算法