溃坝水流数值模拟研究进展史宏达,刘　臻(中国海洋大学工程学院,山东青岛　266071)摘要:溃坝问题在水利工程的设计管理中具有重要地位,也是广大学者长期以来一直关注和研究的课题。

回顾和总结了国内外对溃坝水流演进问题的研究进展:介绍了溃坝水流的数学模型及解析解法存在的困难,进而讨论了数值解法的最新进展;论述了求解溃坝水流一维问题的有限差分法、近似黎曼解的G odunov 格式法、Boltzmann 法、K FVS 法和二维问题的T VD 格式法、间断有限元法、有限体积法、特征线法,并分析了各种方法的适用范围和优缺点,及讨论了限制函数的使用;介绍了利用自由水面追踪方法计算溃坝水流的研究进展,并根据目前存在的不足和实际工程的需要,提出了进一步研究的方向和发展趋势。

关　键　词:溃坝水流;数值模拟;研究进展中图分类号:T V1391231;G 353111 文献标识码:A 文章编号:100126791(2006)0120129207收稿日期:2005206220;修订日期:2005209221作者简介:史宏达(1967-),男,浙江宁波人,中国海洋大学副教授,主要从事工程水动力学研究。

E 2mail :hd 2shi @2631net在水利工程的设计和管理中预报坝体溃决这一灾害性水流现象十分重要。

溃坝计算是对水库和堤防的失事影响做出定量估算,并合理确定水库或堤坝防洪设计标准以及避险措施的有效手段。

溃坝计算的主要内容是算出溃坝坝址的流量和水位过程线,以及下游洪水演进过程中沿程各处的流量、水位、流速、波前和洪峰到达的时间等。

溃坝水流的构成复杂,通常包含激波,亚临界流,超临界流等区域。

通过数值解与试验数据比较,认为浅水方程能够较好的描述溃坝水流。

问题最终归结为求解控制水流运动的非恒定流拟线性双曲型偏微分方程组的有间断问题。

1　溃坝水流的控制方程对于一维溃坝问题,瞬间全溃引起的不稳定流动可视为一维流动,如果假定为静水压力分布和小底坡,则可用圣维南方程[1～4]描述。

对于二维溃坝问题,在静压假定和忽略风应力和柯氏力的条件下,描述溃坝洪水演进的二维控制方程为浅水方程[1～4]。

2　溃坝问题的理论解法尽管数值计算方法在模拟溃坝水流运动方面取得了一定成功,但其效果及优劣通常需要将其与解析法所得结果作比较来判断。

因此,解析解的重要意义是不容置疑的。

1957年,Stoker 将坝址流态分为连续波流、临界流和不连续流三个流态,推导了矩形河谷和下游有水但起始流速为零情况的瞬间全溃坝址峰顶流量公式[5]。

1982年,谢任之在吸收前人经验的基础上,综合连续波与间断波的解法,用抛物线概化河谷断面,去掉对下游水深和流速的控制条件,推出了“统一公式”,可用于各种情况的瞬间全溃的坝址峰顶流量计算,并给出了便于查用的表格[6]。

伍超从溃坝决口形状的任意性出发,定义了断面形态组合参数,提出了组合参数的分离方法,定义了溃坝特征数,探讨了相似性的解结构,建立了一个新的数学模型,反应了真实发生的复杂的溃坝决口的水力特性[7]。

1995年,谢任之对平底无阻力解进行了简化和延伸化的扩展,研究了无限水体的平底有阻力河床瞬间全溃的一阶和二阶渐近解,并提出了有限水体的平底有阻力的渐近解[8]。

第17卷第1期　2006年1月　水科学进展ADVANCES I N W ATER SCIE NCE V ol 117,N o 11　Jan.,2006　031水科学进展第17卷3　溃坝水流的数值模拟解法311　数值模拟溃坝水流的主要途径溃坝问题的特殊性是指其所对应的物理流场中存在间断波,该性质使数值研究溃坝水流具有特定的困难,多数算法常常失效。

溃坝水流数值模拟的主要途径有两条:激波拟合法(Shock2fitting)和激波捕捉法(Shock2cap2 turing)[9]。

激波拟合法是在光滑流动区对圣维南方程组求解,而在涌波两侧则通过间断条件将水流正确地衔接起来。

这类方法虽然精度较高,但其计算复杂,编制程序不便。

激波捕捉法的基本出发点是,若使用与守恒律微分方程组相容的守恒型差分格式,则所得差分解在间断两侧自动满足间断条件,因而不论解中是否存在间断,可以不加区别地统一进行计算,不必进行激波拟合的特殊处理。

捕捉法与拟合法的优缺点恰好相反。

目前常见的溃坝水流的模拟方法多是根据实际情况选择其中一种,或是结合使用。

312　控制方程的离散有限差分法(FDM)着眼于求解区域剖分节点上的函数值,方法简便、灵活,离散格式丰富多样,在收敛性、稳定性等理论研究方面也较为完善。

但由于计算中对求值节点要求规则分布,因而往往不能适应复杂的几何求解域。

此外,传统差分格式不能同时满足虚假振荡的抑制和高精度的要求,常常不是过分耗散就是数值振荡剧烈。

有限元方法(FE M)基于微分方程的弱解形式和广义变分原理,网格剖分能适应具有复杂几何形状的求解区域,在剖分单元上用形函数插值逼近求解,但处理大变形间断问题时会遇到困难。

有限体积法(FVM)在一定程度上吸收了FDM与FE M的长处,并克服了其缺点。

FVM从控制体的积分形式出发,对求解区域的剖分同FE M一样具有单元特征,能适应复杂的求解区域,离散方法具有差分方法的灵活性,并对间断解的适应性强。

随着研究的深入,学者们仍在探索溃坝问题数值解的新途径,或改进原有方法,或提出全新概念,并已在不同程度上取得成功。

313　一维溃坝问题(1)传统差分方法　在流体计算格式中,学者们首先想到的是中心差分格式,但中心差分是无条件不稳定的。

Lax和Friedrichs在中心差分格式中添加了一个二阶粘性项后,该格式才成为稳定格式,但耗散过强,致使数值解过分光滑,间断现象在很宽的范围内被展平[10,11]。

迎风格式在潮流计算中被广泛应用,该法认为:在急流中,只有上游信息传递至下游;在缓流中,上游与下游则相互影响。

因此,在该计算格式中,缓流部分采用根据上下游影响大小而定的加权中心差分,急流部分则采用向后差分。

该格式违反了熵条件,间断附近存在非物理解,故对间断的模拟是不准确的[10]。

MacC omack格式[12]分预测与校正两步实现:预测步用空间后差,校正步对平均解按通量前差校正;也可预测步采用前差,而校正步用后差;或者奇偶时间步轮流执行以上两种方案。

基于MacC ormack预测-校正技术的隐式数值格式可用于计算二维溃坝水流问题[13,14],但计算结果在间断附近有数值振荡,误差大且与Δt有关,Δt无限时不收敛于准确的恒定流解。

Lax2Wendroff格式是时间前差、空间中差(FT CS)格式,并附加粘性项。

该格式虽对间断有较高分辨率,但因为所含格式粘性不足,跨间断两侧取中心差分时与泰勒展开的前提不符,故在间断附近出现虚假的数值振荡(过冲或欠冲),需要加上适当人工耗散后才能保证计算的稳定[9,12]。

(2)近似黎曼解的G odunov型格式　溃坝问题可看作是满足熵条件的有间断黎曼问题。

G odunov格式是由Riemann问题得到启发而产生的[15],其基本思想是将各离散点上的值看作该值在离散点邻域内的平均值,即,将离散值看成某台阶函数。

于是,在离散点之间构成一系列间断,也形成了一系列Riemann问题,该间断在经过Δt时刻传播以后,各离散点上的值再次使用其邻域内的平均值,并重复进行相同过程。

R oe格式、Osher格式都可归结为G odunov型格式中,它们均含适应性格式粘性,不必附加人工粘性便可捕捉间断。

由于迭代求解黎曼问题时存在细节损失,故可利用具有足够精度且简便的近似解。

该格式在模拟溃坝水流时是和谐的,且分辨率较高[16],但一阶格式间断在一定程度上会被抹平。

高分辨率E NO [17]方法提出后,节点模板采用逐次扩展的方法选择,并在各阶差商绝对值极小的选择原则下实现了高分辨率和无振荡的效果。

一种基于E NO 模板选择技术的G odunov 型差分格式通过在时空控制体上的积分及自适应调节,将原方程的时间离散转换为空间离散,将其构造的高阶精度差分格式应用于一维溃坝水流问题[18],精度有所提高,但仍存在数值振荡。

(3)Lattice Boltzmann 方法　Lattice Boltzmann 方法与常见的从流体动力学微分方程出发的数值算法不同,它用离散的粒子在规则格子中的简单运动来模拟流体力学的复杂现象。

在离散的格子里,用虚拟的粒子对应某一状态的粒子分布函数。

对粒子分布函数进行相应计算,可得到质量和动量,从而分布函数的演变也就决定了流体运动的演变过程。

以这种方式建立的简化模型可以相当精确地逼近流体力学方程,进而对溃坝涌波进行模拟[19,20]。

若假定粒子运动碰撞的平均影响对分布函数的改变量与其偏离平衡分布的程度成正比,则可得到Lat 2tice BGK 格式[21],以此为基础建立的一维溃坝模型比文献[20]有一定的进步。

Boltzmann 方法严格满足熵条件,适于平行计算,易于延伸至多维,但目前仅局限于对缓流的模拟,而对急流的模拟却不够成功。

(4)流矢量分裂法　流矢量分裂法(FVS )方法是计算Euler 方程的常用方法,其基本思想是将Euler 方程中的流矢量分成正、负两部分,采用不同格式计算Euler 方程中的对流项。

通常的分裂方法是以方程的特征值或Riemann 问题的解为基础[22,23],而K FVS 方法则是假定可压缩气体分子的速度分布函数是局部平衡的Max well 分布函数,然后分别计算出沿正负两个方向运动的分子流矢量。

这种流矢量分裂方法的物理意义明确,可以看作是Boltzmann 方程的迎风差分格式[24]。

由于描述溃坝水流的浅水长波方程与可压缩无粘气体方程类似,通过对比,可得到模拟溃坝水流的K VFS 方法[25,26]。

该方法计算精度高,在间断点处基本上无振荡,但存在一定的数值耗散[10]。

314　二维溃坝问题(1)T VD 格式　近十多年来,由于气体动力学领域中精确计算激波的需要,在求解可压缩流体的欧拉方程组和纳维尔-司托克斯方程组中获得很大成功的T VD 格式[27]被广泛地应用于浅水动力学,特别是在溃坝问题的计算之中。

T VD 即T otal V oriation Diminishing ,所谓T otal V oriation 是指函数值上升总量和下降总量之和TV (u n )=Σi|u n i -u n i -1|。

为了避免间断两侧总变差的增加,其衰减格式应被设计为:在没有外力(非齐次项和边界条件)作用的条件下,随着时间的推移,总变差不断减小,使最后可能的最小值T V (u n +1)≤T V (u n )。

严格讲,该过程应称为“全变差非增”(T VNI ),只有在不等式严格成立时才称为T VD 。