平面二维水流泥沙数值模拟

第13卷第6期

2002年11月 水科学进展ADVANCESINWATERSCIENCEVol 13,No 6 Nov.,2002

平面二维水流泥沙数值模拟

张细兵,殷瑞兰

(长江科学院河流研究所,湖北武汉 430010)

摘要:基于三角形网格划分,采用有限元方法,建立了河道平面二维水流泥沙数学模型。采用

质量集中的简化处理和预估校正的时间推进算法,较好地解决了有限元计算存储量和计算速度

问题。以空腔流和突扩段两种情况为例对模型进行了检验计算,结果表明,模型能较好地模拟

河道水流泥沙运动及河床变形情况,且计算稳定性好、速度快、精度较高。

关 键 词:水流泥沙;平面二维;有限元;数学模型

中图分类号:TV142 文献标识码:A 文章编号:1001-6791(2002)06-665-05

随着计算机技术的发展和广泛应用,河道水流泥沙数学模型得到了迅速发展,目前已形成多种计算方法。其中,有限元法可采用无结构化网格,能够较好地适应不规则的河道边界,因此有限元方法是一种理想的数值计算方法。

自20世纪70年代以来,有限元法在流体计算中有了长足的发展。日本学者Kawahara针对隐式有限元法计算量和存储量大的缺点提出了 选择参数质量集中二步显示有限元法 ,并成功地运用此方法来解决河口、海湾等方面的问题。接着,美国密西西比大学的王书益教授将该方法运用于河口海湾的泥沙输移问题也获得了成功。在前人研究基础上,本文首先采用网格自动剖分系统对计算域进行三角形网格划分,推导得到了基于三角形单元的有限元方程。采用了质量集中的简化处理方法,化方程系数矩阵为三对角阵,采用了预估校正的时间推进算法,迭代求解方程组,从而大大减少了存储量和计算量,使得有限元法能很好地运用于河道水流泥沙问题的数值模拟。[2][1]

1 基本方程及定解条件

平面二维水流方程

Zs (HU) (HV)++=022x ++U+V=-g-+ t(1)(2)

收稿日期:2001-10-22;修订日期:2002-01-10

基金项目:国家自然科学基金和长江水利委员会联合资助重大项目(50099620);长江科学院基金项目作者简介:张细兵(1976-),男,湖北鄂州人,长江科学院河流研究所助理工程师,主要从事水力学及

河流泥沙动力学研究。

666水科学进展第13卷

22y +U+V=-g- + t+ x y(3)

悬移质泥沙扩散方程

+U+V= +-(S-S)s*H x y

Zb = (S-S*)22(4)(5)河床变形方程

式中 U、V分别为垂线平均流速在x、y方向上的分量;Zs、Zb和H分别为水位、河底高程和水深;g为重力加速度; 为水的密度; y方向上的分量:( x、y分别为底部切应力在x、x; y)=U+V(U;V);C为谢才系数; t为水流紊动粘性系数,由零方程紊流模型确定;S和S*CH

分别为垂线平均含沙量和挟沙力; s为泥沙紊动扩散系数; 为泥沙沉速; 为床沙干容重; 为悬移质泥沙恢复饱和系数。

水流挟沙力和分组挟沙力级配采用李义天方法进行计算。

数值计算定解条件为:进口给定流量、含沙量及级配,出口给定水位。对固壁边界,取法向流速为0。初始时刻给出地形、水位、流速和含沙量等物理量的初始值。[3]

2 有限元方程

为进行有限元计算,首先需对计算区域进行网格划分。模型基于波前推进法思想,建立了河道三角形网格自动剖分系统。即通过编制通用的程序,只需输入少量信息(边界线及地形),便可自动生成充满整个计算域的三角形网格,并配以自动绘图及屏幕显示。该系统能自动获得网格的节点坐标、高程及单元关联信息,从而解决了有限元计算前处理工作量大的困难。

整个计算域共划分为NE个单元,NP个节点,节点整体编号为i,i=1,2, ,NP。对水沙方程(1)~(4),用Calerkin加权余量法进行积分,将单元形函数

经整理得到如下有限元方程:

Aij

AijdZsjdt=-D1ij(HU)j-D2ij(HV)jdUj xj=-BijUj-gD1ijZsj-Aij- t(Cij+Eij)Ujdt(6)(7)

(8)

(9)[5][5][4]代入积分方程,dVj yjAijdt=-BijVj-gD2ijZsj-Aij- t(Cij+Eij)VjAij

其中Aij=

B2ijkdSj=-BijSj- Aij[(S-S*)/H]js(Cij+Eij)Sj- dt = dxdy , j

ki

j idxdy, Bij=B1ijkUk+B2ijkVk, B1ijk=Cij= k i jdxdy, i j i jdxdy,+

第6期张细兵、殷瑞兰:平面二维水流泥沙数值模拟667

D1ij= i jdxdy, D2ij= i jdxdy, Eij= j j idx- idy 为形函数;i、j、k为节点整体编号。

3 数值解法

为便于书写,采用通用变量P来代替式(6)~(9)中的变量Zs、U、V和S,用FP表示方程中等号右边项,并采用了质量集中的简化处理方法,即用对角矩阵A ij来代替Aij,则上述方程可化成如下统一的形式为

dPj=FPi(10)dt

对式(10)的求解,模型采用了预估校正的时间推进算法,即用二阶显式Adams公式进行预A ij报,用隐式梯形公式进行校正,来构造Adams二阶PC公式。离散后的方程不需联接,计算过程稳定性好,时间步长可取得较长,同时还有效避免了数值振荡。

n+1*n tnn-1预报 A ijPj=A ijPj+(3FPi-FPi)2

校正 A ijPj

n+1*n+1n+1**nnn+1*=A ijPj+2(FPi+FPi)[6](11)(12)n+1**为提高计算精度,可进行迭代Pi|< ,则令Pi=Pi

n+1n+1**:给定误差 ,对于所有的1 i NP,若|Pin+1*-,否则令Pi=(Pin+1**+Pin+1*)/2;转到第二步校正,继续迭代,直到满足精度要求为止。当求出Hi、Ui、Vi和Si后,代入河床变形方程便可求得冲淤变形后的河床高程Zbi,这样便完成了一个时段的计算。

以上角标n为时段;*为预报值;**为校正值;无*标记的为该时段的计算结果。

由水位变化引起的自由水面边界变动问题,模型进行了动边界处理,即在每个 t时段,对所有节点按水深区分为水域与陆域节点,并给陆域节点一很小富余水深(如Hmin=0 001m)。4 计算实例

天然河道中的盲肠河段和突然放宽段往往容易形成回流,对回流区的水流运动及泥沙淤积特性进行研究是很有意义的。为检验模型可靠性,本文选用两种情况下的概化模型来进行计算。

空腔流流场计算:主槽为一顺直矩形渠道,中部有一空腔,只进行恒定流流场计算。初始条件:流速为0,水深1m。边界条件:进口流量Q=200m s,出口水深H=1m。主要参数:n=0 02, t=10s。图1为流场示意图。图中水流从空腔下游侧壁进入回流区,并形成一顺时针回流,整个流场水流平顺。图2为流速等值线图,由图可知,空腔回流区水流运动具有不对称性和不封闭性。从定性上看,计算结果合理。

突扩段河床冲淤计算:横断面为矩形,底坡0 1 ,进行水流及泥沙冲淤计算。初始条件:流速为0,水深1m。边界条件:进口流量Q=200m s,含沙量S=2 5kg m;出口H=1m。主要参数:n=0 02, t=10s,挟沙系数k=0 2,悬沙及河床组成均划分为8组。图3为计算稳定后流场示意图,图中水流平顺,并在突然放宽处形成一顺时针回流。图4为计算35000个时段(4d)后33[7]3

668水科学进展第13卷

图1 空腔流流场示意图 图2 空腔流流速等值线

Fig 1Flowfieldofthecavityrecirculation Fig 2Velocityisoplethsofthecavity

circumfluence

图3 突扩段流场示意 图4 河床冲淤分布

Fig 3Flowfieldofsingle-sudden-expansionchannel Fig 4Depositionanderosionofriver-bed

河床冲淤分布图。从纵向看,自上而下,河床由冲变淤。从局部看,河道突然放宽断面附近,右岸靠近突变处冲刷较左岸剧烈,因为右岸突然放宽造成跌水,使附近流速增大所致;回流区由于流速减小,淤积明显;回流区下游小范围内表现为冲刷,主要由于上游挟带的泥沙在回流区大幅度淤积后,使其下游含沙量减小,故而表现为冲刷。另将数模计算结果与曾经做过的突扩段水槽试验结果进行对比分析得知,两者冲淤分布规律及相对冲淤幅度基本一致。由此说明该模型能较好地反映突扩段河床的冲淤变形规律。

5 结 论

有限元法不仅具有稳定性好、灵活性强和精度高的优点,而且网格划分灵活,可对局部区域任意加密,因此能较好地模拟不规则河道边界,很适合于天然河道水流泥沙问题的数值模拟。

有限元法建立的平面二维水流泥沙数学模型,采用了质量集中的处理方法及预估校正的时间推进算法,较好地解决了计算存储量和速度问题,且计算过程稳定,计算速度较快。

本文重点研究模型的可行性和可靠性,因此以空腔流和突扩段两种概化模型为例进行了计算(天然实测资料的验证计算可参阅文献[8])。计算结果表明,模型能较好地模拟河道的水流泥沙运动及河床变形情况。

致谢:数模研制过程中,得到了长江科学院黄煜龄教授级高级工程师、董耀华高级工程师和武汉大学韦直林教授的细心指导,在此一并致谢!

第6期张细兵、殷瑞兰:平面二维水流泥沙数值模拟669参考文献:

[1]KawaharaM,HiranoH,TsubotuK,etal.Selectivelumpingfiniteelementmethodforshallowwaterflow[J].Interna-

tionalJournalforNumericalMethodsinFluids,1980,(2):89-112

[2]WangSY.Methodologyonmodelingoflooseboundaryhydraulicsystems[A].2ndIntSymp OnRiverSedimentation

[C].1982.219-239

[3]李义天,等.河道平面二维泥沙数学模型研究及应用[A].第三届全国泥沙基本理论研究学术讨论会论文集

[C].北京:中国建材工业出版社,1995

[4]张细兵.河道有限元网格自动剖分方法研究[J].长江科学院院报,2002,(3):19-21.

[5]JJ康 纳,等.流体流动的有限元方法[M].北京:科学出版社,1981.

[6]韦直林.河道水流泥沙问题的一种有限元解法[J].武汉水利电力学院学报,1990.

[7]董耀华.空腔回流区水沙特性的计算分析[J].泥沙研究,1999,(2):34-39

[8]张细兵.长江中游监利河段平面二维水沙数模研制及验证报告[R].长江科学院,2000

Planar2-Dflowandsedimentmathematicalmodeling

ZHANGX-ibing,YINRu-ilan

(YangtzeRiverScientificInstitute,Wuhan430010,China)

Abstract:Basedontriangulargrids,theplanar2-Dflowandsedimentmathematicalmodelisestablishedbythefiniteelementmethod Sincethemassconcentrationandthepre-estimatedandverifiediterationareusedtosolvetheequations,thequestionsofmemoryunitsandcalculationvelocityaresettledquitewell Simulationcalculationiscarriedoutonthecavityrecirculationandsingle-sudden-expansionchannel Theresultsshowthatflowfieldandsedimentdistributionaresimulatedsuccessfullyinthemodelandthecalcula-tionsarestable,fastandaccurate

Keywords:flowandsediment;planar2-D;finiteelementmethod;mathematicalmodel

TheprojectissupportedbyNationalNaturalScienceFundofChinaandYangtzeWaterResourcesCommission(No.50099620),andbytheFundofYangtzeRiverScientificInstitute.

本文来源：https://www.bwwdw.com/article/7234.html