基于无结构网格有限体积法的Roe隐格式探讨
- 格式:pdf
- 大小:513.98 KB
- 文档页数:6
溃漫堤洪水多维耦联数值模型及应用苑希民;王亚东;田福昌【摘要】针对山洪沟防洪标准偏低、遭遇暴雨洪水时溃漫堤风险较大且较难准确预测的问题,开展洪水多维耦联数值模型研究.基于圣维南方程与VOF法的标准k-e 双方程原理,融合Abbott六点隐格式法、非结构化网格Roe格式的单元中心显式有限体积法与结构化网格PISO算法优势,系统构建溃漫堤洪水多维耦联数值仿真模型.以贺兰山汝箕沟为研究对象,模拟其遭遇百年一遏洪水时,沟道一维、溃口三维与淹没区域二维的多维水流动态耦合演进过程及影响情况.结果表明,溃漫堤洪水多维运动过程与流态模拟效果较好,所建多维洪水耦联数值模型计算精度较高,对多流态洪水动态耦合精确计算与淹没风险准确评估具有重要意义.【期刊名称】《天津大学学报》【年(卷),期】2018(051)007【总页数】9页(P675-683)【关键词】溃漫堤洪水;多维耦联;动态耦合;数值模型;汝箕沟【作者】苑希民;王亚东;田福昌【作者单位】天津大学水利工程仿真与安全国家重点实验室,天津300350;天津大学水利工程仿真与安全国家重点实验室,天津300350;天津大学水利工程仿真与安全国家重点实验室,天津300350【正文语种】中文【中图分类】TV131.2随着自然环境日趋恶劣,极端天气常有发生,致使超常洪水较以往更加频繁.堤防在防洪减灾方面起到了至关重要的作用,但当遭遇超标准洪水时,漫溃堤风险加大.水流数值模拟技术具有计算精度高、水力特征信息提取便利等明显优势,可为防汛指挥决策提供可靠的技术支撑.近年来,各国学者对洪水演进数值模拟技术进行了广泛深入的研究.在一、二维水动力模型及其耦合研究方面,Blade等[1]基于数值通量耦合方法,分析模拟了天然河道水流运动情况;张靖雨等[2]采用一、二维水动力耦合模型对溃堤洪水与内涝洪水进行叠加分析,模型能够合理反映保护区洪涝风险分布特征,结果较为可靠合理;苑希民等[3]将全二维气相色谱理论与全二维水动力模型结合,为洪水演进模拟提供了新思路,并将该模型成功应用于黄河内蒙段河道与左右岸灌区洪水风险动态耦合模拟;施露等[4]运用Mike Flood建立了水文学模型、一维和二维水动力耦合模型,模拟分析了中小河流洪涝风险情况,研究结果具有重要价值.在二、三维水动力模型及其耦合研究方面,练继建[5]基于VOF法建立了溃坝水流运动三维数值模型,并对河道丁坝附近溃坝波传播的特性及其内部结构进行了模拟研究;王桂芬[6]提出了二、三维嵌套数学模型概念,并应用于天津新港抛泥区的潮流计算;Zounemat等[7]也对二、三维耦合水动力模型进行了研究,其耦合方式为水体表面采用二维浅水模型计算,而深水水体采用三维水动力模型计算.综上所述,国内外学者对溃堤或漫堤洪水一、二维耦合模拟研究较多,少数学者对二、三维水动力耦合模型进行了理论与应用研究,而对于溃漫堤洪水多维耦联计算模型的研究却少见报导.实际堤防溃决洪水分流是极其复杂的三维非恒定流运动过程,堤防及溃口边界对水流流态与分流过程影响很大,通常所建一、二维水动力耦合模型仅能计算河道与溃口断面平均水力要素,难以精确获取纵向水流变化规律与分流特征.因此,本文拟融合不同维度水动力数值模型计算优势,基于河道实测断面与DEM数字高程,采用Abbott六点隐格式法、Roe格式离散法与PISO算法,研究建立河道一维、平面区域二维与溃口分流三维的洪水多维耦联数值仿真模型,较为准确地计算沟道一维洪水水面线、淹没区二维洪水淹没特征与溃口分流三维洪水运动过程,并将其应用于贺兰山汝箕沟遭遇百年一遇超标准洪水时溃漫堤洪水演进动态变化的数值模拟.河道一维水动力学模型采用描述明渠非恒定流的圣维南方程,如式(1)和式(2)所示[8].式中:A为过水面积;q为区间流量;Q为断面过水流量;g为重力加速度;t为时间;x为沿程距离;a,为动量修正系数;Z为河道断面水位;K为流量模数.采用稳定性与精准度较高的Abbott六点隐式格式离散上述方程.二维平面淹没区域,垂向洪水流速等水力要素的变化远小于水平方向,故采用二维浅水方程对其进行描述计算.具体如式(3)~式(5)所示[3].式中:H为水深;Z为水位;q为源汇项流量;M与N分别为x、y方向上的垂向平均单宽流量;u、v分别为垂向平均流速在x、y方向的分量;n为曼宁系数;g 为重力加速度.方程不考虑科氏力和风力影响.在保证计算区内水量和动量守恒的基础上,采用Roe格式的单元中心显式有限体积法对控制方程进行离散求解,并通过干湿边界处理技术、线状地物与网格耦联方法,优化淹没区二维水动力模型.1)干湿边界处理优化为提高模型计算效率及其稳定性,采用干湿边界理论判断网格淹没状态,通过设置干水深及湿水深,确定洪水移动边界.当网格单元计算水深值小于干水深时,单元通量为0;当计算水深值大于干水深而小于湿水深时,网格单元只计算质量通量,无动量通量;而当其大于湿水深时,网格单元质量通量与动量通量均进行计算处理.具体关系如图1所示.2)线状地物与网格耦联方法为在保护区内准确定位线状地物,将实测地理参考点线性连接,以其为内边界创建网格,可使堤防、道路等特殊边界与非结构化网格密切耦联.该方式可避免线状地物与网格界面相交,使得参与计算折线与实际线状地物重合,利于模型稳定运行并提高计算精度.图2(a)为常规线状地物耦联方式,图2(b)为改进后的线状地物耦联示意图.洪水经由河道通过溃口进入堤外保护区域时,水流方向发生急剧改变,堤防受到溃堤水流的强烈冲击,水流流态复杂,三维特征十分明显.一维水动力模型计算溃口分流洪水时,只能得到纵向上无梯度分布的河道及溃口断面平均水力要素,与实际情况差异较大.故本文采用具有较强界面捕捉能力的VOF方法的标准k-e双方程模型对溃口洪水分流过程进行模拟计算[9].基本控制方程如式(6)和式(7)所示.式中,cμ为经验值;p为修正压力.流体自由表面处理方程如式(8)~式(10)所示.输运方程为式中:F为流体体积函数,F=1表示控制体内充满流体,F=0表示无流体;ui为xi方向上的流速.k-e紊流方程为模型常数式中:μ为分子动力黏性系数;μt为紊流黏性系数;k为紊动动能;G为紊动能生成率;ε为紊动耗散率;σk、σe分别为k、ε的紊流普朗特数;C1z、C2z为经验常数.在时间与空间上分别采用Crank-Nicholson格式和LUD格式离散求解控制方程.为使其能更好地满足动量方程和连续性方程,采用PISO算法对方程进行修正.具体见文献[9].河流溃漫堤洪水耦合计算关系如图3所示.1.4.1 一、二维水动力耦合模型当河道水位超过堤顶高程时,发生洪水漫溢现象,可将漫溢洪水作为宽顶堰流处理[8],建立一、二维水动力耦合模型,采用侧向连接方式实现河道与堤外保护区域的动态耦联,利用Villemonte堰流公式计算洪水漫溢流量(如式(11)所示),模拟漫溢洪水在堤外平面区域的动态演进过程及淹没特征.式中:q为漫堤流量;C1和C2=为堰流系数;W为堰顶宽度;Hu为上游水位高程;Hd为下游水位高程;Hw为堤顶高程.1.4.2 一、三维水动力耦合模型为减小河道下游水位顶托影响,采用河床具有正向坡的断面作为一、三维水动力模型的耦合衔接面,边界条件设定如下.1)上游入流边界条件以流量过程作为河道上游入流边界条件.为保证一、三维水动力耦合模型计算结果的准确性,考虑将衔接面河道水位作为模型动态耦合衔接条件,如式(12)、式(13)所示.由经验公式确定k和ε,如式(14)和式(15)所示.式中:Q1为衔接面上游流量;Q3为衔接面下游流量;h1是一维河道水位;h3为三维河道水位;V3是三维河道入口水流流速;l为一、三维耦合衔接面长度.2)下游出流边界条件以河道末端断面水位-流量关系作为下游出流边界条件.为确保模型计算稳定,连同溃口上、下游一定范围内的河道及淹没区作为三维水动力模型建模范围.在此可选择溃口下游地形平坦、水流流态较为平稳的淹没区域作为二、三维水动力耦合模型的过渡衔接面,如图4所示.1.4.3 二、三维水动力耦合模型在计算二、三维水动力耦合模型衔接通量时,三维模型沿水深方向具有多层计算节点,而二维模型仅计算水平方向速度分布,两者相差较大.计算二维模型对三维模型的反馈影响,需获得三维洪水纵向水力要素的变化规律,以此来补充其垂向分层边界条件.为确保流速和水位等变量信息在二、三维模型耦合衔接面上的交互传递[10](如图5所示),在压强及流量相等的基础上,对三维洪水纵向速度分布进行平均处理,如式(16)和式(17)所示.对二、三维耦合模型过渡衔接处回水反馈影响,按式(18)和式(19)对平均流速进行垂线化分布处理.式中:m为流速分布指数;H为水深;z为距离河床高程;U、V分别为x、y方向基于水深平均的流速;u(z)、v(z)分别为沿河床高程各层节点的x、y方向流速,水流流经二、三维耦合过渡衔接区域后,垂向无速度梯度.为验证多维耦联模型的准确性,本文选用Soares等[11]于2002年做的物理试验对耦联模型进行验证.试验模型由水库段和渠道段组成,其中渠道宽0.495,m,渠底高于水库底部0.33,m,距溃口下游3.92,m处有一90°直角弯,渠道的曼宁糙率系数为0.011,末端为自由出流;水库宽2.44,m、长2.39,m,初始水位0.58,m,高出渠道底部0.25,m.该试验主要研究水库坝体瞬间全溃时水流在渠道中的演进过程.在此分别建立多维耦联模型和一、二维耦合模型,如图6所示.首先,从上游至下游依次建立二维(A区)、三维(C区)和一维(D区)模型.当水流流经直角弯时,水体流态空间分布较为复杂,故C区域进行三维模拟.其次,因该耦联模型只有二、三维耦合和一、三维耦合,缺少对多维耦联模型中的一、二维耦合模型部分的验证,在此单独建立一、二维耦合模型,如图6(b)所示.模型处理:①采用非结构化三角形网格剖分二维模拟区,最大网格面积0.009,m2;②二、三维过渡段区(B区)采用0.02,m×0.02,m×0.1,m,C区采用0.02,m×0.02,m×0.01,m结构化正交网格进行剖分;③一维模拟计算区,渠道断面距离0.02~0.04,m,在直角弯进行插值加密.模型上边界条件为给定的初始水深0.58,m,下边界为出口端水位-流量关系.为使得水库水面水平,在渠道处设一控制闸门,当其上游水位等于0.58,m时,闸门突然提起;模型计算步长为0.01,s,二、三维模型最小步长为0.1,μs,模拟时长t=10,s.图7为水深验证结果,提取模型在t=5,s和t=7,s时的模拟水深值,将其与实测数据进行对比分析.由图可知,二维模拟计算值与实测值最为接近,因二维模拟区水流受直角弯影响较小,流态较为平稳,故该段模拟结果较为准确;而三维模拟水深总体较小,主要原因有:①涡黏性系数各向同性假定,即模型自身限制;②因模型尺寸总体较小,水体表面张力影响较大,而数学模型并未考虑表面张力作用,致使模型计算结果误差较大;对于一维渠道,由模拟结果可知,水深变化总体趋势基本一致,其中该段上游模拟值与实测值较为接近.结果表明,本文所建多维耦联数值计算模型模拟结果较为准确可靠.本文选择贺兰山汝箕沟出山口至长胜墩拦洪库段作为研究对象(如图8所示),该河段地处贺兰山东麓中段地区石嘴山市,为宁夏两大暴雨中心之一.区域地势总体为西南高东北低,洪水自出山口流出后,经行洪沟道汇入长胜墩拦洪库.行洪沟道长约11.01,km,左右岸防洪堤防顶宽约5.00,m,沟道平均比降为3.29%,防洪标准为20年一遇.研究区域内主要线状地物有G110、姚汝公路和石银高速公路等,因姚汝公路桥上游转弯处河道弯度较大、水流对右侧凹岸淘刷作用较强,溃堤风险偏高,且堤防一旦溃决,将对溃口下游居民及大量农田造成严重影响,故将此处设定为溃口位置.根据河道洪水、溃口分流洪水、堤防漫溢洪水与淹没区洪水流动各自表现出的不同维度特征,建立溃漫堤多维洪水动态耦联数值计算模型.为减小河道一维与溃口三维水流动态交互影响,于溃口断面上游200,m处设置一、三维耦合面域,并通过二、三维耦合过渡衔接区,实现溃口三维水流与淹没区二维洪水的动态耦合,模拟河道洪水水面线、溃口三维洪水分流过程、溃漫堤洪水运动特征及其淹没风险.研究区网格剖分及模型衔接关系,如图9所示.汝箕沟沟道一维水动力模型计算河段长度11.01,km,设置控制断面192个.根据《石嘴山市汝箕沟大武口段治理工程初步设计报告》内容论述,设定汝箕沟沟道综合糙率为0.034.为提高模型计算稳定性及运算效率,设定汝箕沟一维水动力模型计算步长为1,s.选择2009年汝箕沟水文站汛期实测洪水过程作为典型洪水,并采用同倍比放大法求得其百年一遇设计洪水流量过程(如图10所示),作为一、三维水动力耦合模型的上游入流边界条件.采用非结构化三角形网格剖分堤外平面保护区域,区域总面积318.68,km2,网格面积控制在100~5,000,m2之间,剖分网格数量17.50×104个,计算节点数8.92×104个.计算区域多为农田耕地,房屋村庄面积较小且零散分布(见图11),因缺乏实测糙率资料,在此依据以往类似区域模型的计算经验[12],可采用不同地理条件的糙率值与面积权重之积,确定计算区域综合糙率值为0.06.为提高模型计算精度与运算效率,设置干水深值为0.005,m,湿水深值为0.1,m,在满足临界值CFL≤1的条件下,设置计算时间步长为0.01~10,s,模型可根据计算效率与稳定性大小,实现计算步长的自我动态调整.针对姚汝公路桥上游转弯处溃口上、下游200,m河道及溃口外延50,m淹没区,采用矩形网格剖分技术构建研究区域三维地形,网格边长为0.60~1.50,m,剖分网格数量74.47×104个,整体建立溃口分流过程三维仿真计算模型.当溃口处河道断面水位等于设计水位1,106.10,m时,堤防瞬间溃决,并参考《洪水风险图编制技术细则(试行)附录》中的溃口宽度经验公式(如式(20)),溃口处河宽77.00,m,溃口宽度计算值为59.98,m,模型设定溃口宽度为60,m,溃口底高程为1,105.87,m.通过模拟溃口三维洪水分流及二、三维过渡衔接区域水流运动过程,实现多维洪水流态的动态耦合与过渡,准确计算溃漫堤分流洪水淹没风险分布特征.式中:Bb为溃口宽度,m;B为河宽,m.根据汝箕沟溃漫堤多维洪水耦联数值模型计算结果分析可知:当汝箕沟遭遇百年一遇洪水时,出山口下游850~2,100,m之间沟段发生洪水漫溢现象,漫溢洪水最高水位超过堤顶高程0.51,m.漫溢洪水流量过程如图12所示,漫溢分流呈现2个峰值,最大分流量为132,m3/s,漫溢总水量为17.70×104,m3,上游来水洪峰过后,河道流量迅速减小,堤防漫溢洪水流量急剧降低,且由于G110及下游沟道堤防阻水壅水作用明显,使得发生漫溢洪水倒灌回流现象,即漫溢流量出现短暂负值,随后水流平稳漫溢流量为零.汝箕沟溃口位于出山口下游3,150,m处,溃口分流流量过程如图13所示.堤防瞬间溃决,溃口分流流量急剧增加,起溃时分流流量为59,m3/s,分流洪峰流量为139,m3/s,溃堤分流总水量为68.35×104 m3.汝箕沟洪水峰值时刻溃口处河道横断面与区域三维流速分布分别如图14和15所示.由图14分析可知:溃口附近流速较大,溃口外侧流速逐渐减小,水面流速低于底部流速;溃口出流水面收缩处流速最大,下游发生水跃现象,水面波动剧烈,流速减小,水深增加,水面逐渐趋于平稳.由图15分析可知:堤防对洪水破堤分流干扰作用明显,表现为侧向出流,且在溃口边缘部位产生漩涡,溃口处流速急剧变化且明显高于其余部位.图16为研究区地形,汝箕沟溃漫堤洪水演进不同时刻淹没水深分布情况如图17所示,与图16对比分析可知:洪水淹没范围不断向地势低洼处扩展,其中以石银高速公路为例,因涵洞位置地形较低,路面较高,使得石银高速靠近溃口一侧发生积水现象,该处水深值较大,阻水作用明显.当洪水演进2,h,漫溢洪水淹没风险较大,主要影响区域为北庄子;当洪水演进6,h,主影响区域包括北庄子、崇岗村、金达煤业、崇岗镇民生服务中心和崇岗村砖厂;当洪水演进12,h,结合图11可知,溃口分流洪水淹没大量农田耕地,影响范围扩展至农三队.按不同淹没水深等级(见表1)统计分析溃漫堤洪水淹没区域面积与受影响道路长度,如表1所示.洪水淹没总面积为6.70,km2,其中溃堤洪水淹没面积为5.84,km2,漫堤洪水淹没面积为0.86,km2,淹没水深小于0.5,m的风险区域面积为6.58,km2,占淹没总面积的98.21%,;受影响道路总长度为4.00,km,其中积水深小于0.5,m的道路长度为3.16,km,占受影响道路总长度的79.00%,.本文基于圣维南方程与VOF法的标准k-ε双方程,建立了溃漫堤洪水多维耦联计算模型.以贺兰山汝箕沟为研究对象,分析模拟了其遭遇百年一遇超标准洪水情况下,沟道一维洪水演进、溃口三维分流以及平面区域二维动态淹没的多维洪水耦合联动过程.结果表明,溃漫堤多维洪水演进过程模拟效果较好,流态较为真实,计算结果基本合理可靠.所建模型在水面追踪与多维洪水动态耦联模拟等方面具有明显优势,横纵向水力要素模拟与多空间尺度处理能力较高,具有很好的推广应用前景.【相关文献】[1] Blade,Gmez-Valentn M,Dolz J,et al. Integration of 1D and 2D finite volume schemes for computations of water flow in natural channels[J]. Advances in Water Resources,2012,42:17-29.[2]张靖雨,徐佳,袁先江,等. 外洪溃堤一、二维耦合模型与内涝模型叠加在防洪保护区内的应用探讨[J]. 水利水电技术,2017,48(5):87-94.Zhang Jingyu,Xu Jia,Yuan Xianjiang,et al. Discussion on application of superpositionof 1-D and 2-D coupling model for dike failure due to outside flood and water-logging model to flood control protection area[J]. Water Resources and Hydropower Engineering,2017,48(5):87-94(in Chinese).[3]苑希民,田福昌,王丽娜. 漫溃堤联算全二维水动力模型及应用[J]. 水科学进展,2015,26(1):83-90.Yuan Ximin,Tian Fuchang,Wang Lina. Comprehensive two-dimensional associate hydrodynamic models for overflow and levee-breach flood and its application[J]. Advances in Water Science,2015,26(1):83-90(in Chinese).[4]施露,董增川,付晓花,等. Mike Flood在中小河流洪涝风险分析中的应用[J]. 河海大学学报:自然科学版,2017,45(4):351-358.Shi Lu,Dong Zengchuan,Fu Xiaohua,et al. Application of Mike flood and waterlogging risks of medium and small rivers[J]. Journal of Hohai University:Natural Science,2017,45(4):351-358(in Chinese).[5]练继建. 溃坝水流在复杂河道中传播的三维数值模拟[J]. 水利学报,2007,38(10):1151-1157.Lian Jijian. Three-dimensional numerical simulation of dam-break flow in complicated river sections[J]. Shuili Xuebao,2007,38(10):1151-1157(in Chinese).[6]王桂芬. 二、三维潮流数学模型嵌套连接技术[J]. 交通与计算机,1988(3):39-40. Wang Guifen. Second-dimensional and three-dimensional trend mathematical model nested connection technology [J]. Journal of Transport Information and Safety,1988(3):39-40(in Chinese).[7] Zounemat K M,Sabbagh-Yazdi S R. Coupling of two-and three-dimensional hydrodynamic numerical models for simulating wind-induced currents in deep basins[J]. Computers and Fluids,2010,39:994-1011.[8]苑希民,薛文宇,冯国娜,等. 溃堤洪水分析的一、二维水动力耦合模型及应用[J]. 水利水电科技进展,2016,36(4):53-58.Yuan Ximin,Xue Wenyu,Feng Guona,et al. A couple one-and two-dimensional hydrodynamic model for analysis of levee-breach flood and its application[J]. Advances in Science and Technology of Water Resources,2016,36(4):53-58(in Chinese).[9]李超. 长距离调水工程高填方渠道溃堤洪水演进情景仿真[D]. 天津:天津大学建筑工程学院,2014.Li Chao. Scenario Simulation on Dike Breach Flood of High Embankment Channel in Long Distance Water Diversion Project[D]. Tianjin:School of Civil Engineering,Tianjin University,2014(in Chinese).[10]黄玉新,张宁川. 二、三维耦合水动力模型研究I:模型的建立[J]. 水道港口,2013,34(4):304-310.Huang Yuxin,Zhang Ningchuan. Research on a couple 2D-3D hydrodynamic model I:Model establishment [J]. Journal of Waterway and Harbor,2013,34(4):304-310(in Chinese).[11] Soares Frazao S,Zech Y. Dam break in channels with 90° bend[J]. Journal of Hydraulic Engineering,2002,128(11):956-968.[12]贺新娟,徐国宾,苑希民. 溃口宽度和倒灌、退水对防洪保护区洪水演进的影响[J]. 水资源与水工程学报,2016,27(2):142-147.He Xinjuan,Xu Guobin,Yuan Ximin. Effect of levee-breach width and water intrusion and recession on flood routine in flood protected zone[J]. Journal of Water Resources & Water Engineering,2016,27(2):142-147(in Chinese).。
非结构化网格下求解二维浅水方程的和谐Roe格式
吕彪;金生;艾丛芳
【期刊名称】《水利水运工程学报》
【年(卷),期】2010(000)002
【摘要】为保证底坡源项和重力梯度项的平衡离散,采用把底坡源项分解为两不同部分单独处理的方法,在非结构化网格上建立了求解带复杂地形的二维浅水方程数值模型.采用Roe格式的计算界面通量,隐式求解摩擦力源项以增加格式的稳定性,并给出计算格式在非结构化网格上满足和谐性条件的证明.通过实例验证了此格式是和谐的,并具有良好的间断捕捉能力和稳定性.
【总页数】6页(P39-44)
【作者】吕彪;金生;艾丛芳
【作者单位】大连理工大学土木水利学院海岸和近海工程国家重点实验室,辽宁,大连,116024;大连理工大学土木水利学院海岸和近海工程国家重点实验室,辽宁,大连,116024;大连理工大学土木水利学院海岸和近海工程国家重点实验室,辽宁,大连,116024
【正文语种】中文
【中图分类】TV131.4
【相关文献】
1.求解二维浅水方程的高分辨率完全松弛格式 [J], 胡彦梅;陈建忠
2.基于修正Roe格式的有限体积法求解二维浅水方程 [J], 刘刚;金生
3.三角形网格下求解二维浅水方程的KFVS格式 [J], 潘存鸿;徐昆
4.求解具有复杂地形二维浅水方程的修正HLL格式 [J], 艾丛芳;金生
5.三角形网格下求解二维浅水方程的和谐Godunov格式 [J], 潘存鸿
因版权原因,仅展示原文概要,查看原文内容请购买。
非结构网格的有限体积法研究张海军李雪(长安大学理学院陕西·西安710064)摘要基于双曲守恒律方程,详细阐述了非结构三角形网格的有限体积方法。
在该方法中,通过运用数值流函数来近似计算线积分,并且详细介绍了常见的三种数值流函数。
时间的离散采用三阶TVD Runge-Kutta方法。
关键词有限体积法非结构网格数值流函数龙格-库塔法中图分类号:V211文献标识码:A0引言有限体积方法可以认为是有限元法和有限差分法的结合,所以有限体积法吸取了有限元网格剖分灵活的优点,克服了差分法网格适应性差的缺点。
二十世纪八十年代以来,由于非结构网格的发展,有限体积法取得了很大的进步,为双曲型守恒律方程的发展提供了很大的空间。
因此,非结构网格下的有限体积法已经成为数值模拟复杂、高速流动的重要方法。
1格式的构造1.1空间离散对于二维标量双曲守恒律方程对计算区域采用规则的三角形网格剖分,以三角形单元A上积分得:A的网格平均值,是三角形的第k条边,是第k条边对应的外法向量。
定义,且定义外法向量的模长为对应边的长度,用中矩形公式近似上式中的积分可得:其中表示第k A的第k条边共边的三角形对应边的重构函数,为数值流函数。
1.2数值流函数的近似数值流函数是近似,常用的有三种形式:第一种形式为算术平均形式:在有限体积法的构造中,人们习惯对时间和空间分别进行处理,时间方向的离散一般采用文献[3]中的TVD Runge-Kutta方法。
2结语在非结构的三角形网格下,详细描述了双曲守恒律方程的有限体积方法,通过运用数值流函数来近似计算积分。
时间的离散用三阶TVD Runge-Kutta方法表达式。
参考文献[1]赵延生.非结构网格的ENO有限体积方法研究[D].长沙:国防科学技术大学研究生院,2004,1-52.[2]朱华君.二维浅水波方程的高阶有限体积格式[D].长沙:国防科学技术大学研究生院,2006,1-54.[3]Fjordholm U,Mishra S,Tadmor E.Energy preserving and energy stable sche-mes for the shallow water equations[C].2009,93-139.数|学|研|究—科教导刊(电子版)·2017年第24期/8月(下)—153。
非结构网格数值计算格式的研究及其在环境水力学中的应用的开题报告标题:非结构网格数值计算格式的研究及其在环境水力学中的应用一、研究背景和意义随着计算机技术的不断发展和水力学研究的推进,非结构网格数值计算格式在环境水力学领域中的应用越来越广泛。
与传统的结构网格数值计算相比,非结构网格数值计算格式具有更高的灵活性和适应性,能够更准确地描述复杂的水动力学现象。
因此,非结构网格数值计算格式的研究已经成为环境水力学领域内的一个重要研究方向。
二、研究内容和方法本研究主要探究非结构网格数值计算格式在环境水力学中的应用。
具体研究内容包括以下几个方面:1.非结构网格的生成方法以及其在水力学中的应用;2.非结构网格数值计算格式的基本原理和常用算法;3.基于有限体积法和有限元法的非结构网格数值计算格式在环境水力学中的应用及其数值效果分析;4.非结构网格数值计算格式的优化方法和相关问题的研究。
研究方法主要采用理论分析和数值模拟相结合的方法。
首先,通过文献综述和理论分析,探究非结构网格数值计算格式的基本原理和常用算法;然后,采用有限体积法和有限元法的非结构网格数值计算格式在环境水力学中进行数值模拟,并通过实例分析来检验数值方法的有效性和准确性;最后,通过对现有算法的分析和对问题的深入理解,深入研究非结构网格数值计算格式的优化方法和相关问题。
三、预期成果和意义研究成果主要包括以下几个方面:1.论文发表:通过此次研究,论文发表至国内外权威学术期刊。
2.参数选取:确定非结构网格数值计算格式在环境水力学中的最佳参数选取。
3.数值模拟:基于有限体积法和有限元法的非结构网格数值计算格式在环境水力学中的应用及其数值效果分析。
此次研究的意义在于提高非结构网格数值计算格式在环境水力学中的应用水平,为环境水利工程的相关设计和研究提供理论基础和技术支持。
基于无结构化网格浅水方程的隐式解法唐岳灏【摘要】In order to improve the numerical stability of shallow water equation calculated by Finite Volume Method, by emplo-ying Roe's approximate Riemann solution to calculate the interface flux and TVD-MUSCL Format to reconstruct the conservation variable, a high efficient implicit computation scheme is derived. On the basis of the unstructured grids, this format improves the computation accuracy to grade 2. It computes the velocity gradient by the area weight and satisfies the stationary hydraulic pres-sure equilibrium by handling the bed slope term. In order to use the implicit scheme for the time integration, the full resolution form of Jacobian matrix is analytically derived, which was solved by Newton-Raphson algorithm iteratively. By the comparison with various numerical studies on dam-breaking cases, this computation method is proved to be stable, compatible and efficient with the capability of accurately capturing the shock wave in dam-breaking problems.%为提高有限体积法计算浅水方程的数值稳定性,采用Roe方法近似Riemann解计算界面通量,利用TVD-MUSCL格式对守恒变量进行重构,推导并建立了高效的隐式计算格式。
基于Roe格式黎曼近似解的二维FVM模型汪梅华;张铭;柳杨;乌景秀【摘要】提出一种基于黎曼近似解Godunov格式的二维FVM模型求解口门区二维通航水力特性,并采用Roe格式计算界面通量.采用水面坡度代表源项中压力项的作用,有利于复杂地形条件计算的稳定性;通过对斜底单元干湿特性的合理划分,确保计算单元的水量和动量平衡及数值计算精度.利用该模型计算分析了衢江梯级塔底枢纽上下游口门区及引航段发电及泄洪条件下的通航水力特性,结果表明Roe格式的有限体积方法计算复杂地形条件下的水流流场稳定性好、计算精度高,为合理制定衢江梯级枢纽调度规则提供数据支持.【期刊名称】《水利水运工程学报》【年(卷),期】2016(000)003【总页数】8页(P27-34)【关键词】山溪型航道;有限体积法;Roe格式;二维水力计算模型【作者】汪梅华;张铭;柳杨;乌景秀【作者单位】衢州市港航管理局,浙江衢州324002;南京水利科学研究院,江苏南京210029;南京水利科学研究院,江苏南京210029;南京水利科学研究院,江苏南京210029【正文语种】中文【中图分类】TV135.4金沙江、乌江、西江、衢江等为我国重要通航水运要道,一大批通航枢纽已陆续建成并投入运行[1-3]。
西部航道具有显著的山溪型航道特征,河道地形复杂,水力特性敏感,电站调峰或大坝泄洪对航运的影响十分突出,枢纽瞬时下泄流量的快速改变恶化枢纽附近口门区、上下游引航道、连接段及下游航道通航水力学条件[4-6],对枢纽下游航运带来安全隐患。
本文提出一种基于黎曼近似解Godunov格式的二维FVM模型求解口门区二维通航水力特性,并采用Roe格式计算界面通量。
该模型源项中采用水面坡度代表压力项的作用,避免了对底坡项的复杂处理,有利于复杂地形条件下计算的稳定性;通过对斜底单元干湿特性的合理划分,实现水位和水深合理转换;基于物理通量严格守恒,确保计算单元的水量和动量平衡及数值计算精度[9-11]。
浅水方程无结构网格有限体积法研究进展
于守兵
【期刊名称】《水运工程》
【年(卷),期】2009(000)005
【摘要】回顾了近些年国内外基于无结构网格的有限体积法模拟浅水流动的发展,着重论述了无结构网格的有限体积法的应用、高性能的法向数值通量的计算格式、物理量的重构和浅水流动静水压力项与底坡项的平衡问题.讨论了无结构网格有限体积法隐式求解、减小局部一维化误差和构造具有明确物理意义的消除虚假流动的计算方法等有待于进一步研究的方向.
【总页数】6页(P15-20)
【作者】于守兵
【作者单位】南京水利科学研究院,江苏,南京,210024;南京水利科学研究院水文水资源与水利工程科学国家重点实验室,江苏,南京,210024
【正文语种】中文
【中图分类】O352
【相关文献】
1.基于非结构网格的不可压N-S方程多矩有限体积法 [J], Bin Xie;Peng Jin;Feng Xiao
2.基于无结构网格的有限体积法和颗粒方法相结合求解湍流的联合PDF方程 [J], 葛海文;朱旻明;胡国辉;李艺;董刚;王海峰;陈义良
3.二维无结构网格浅水波方程的高阶非振荡有限体积法 [J], 袁美;龙巧云;宋松和
4.二维浅水波方程的非结构网格ENO型有限体积法 [J], 朱华君;宋松和
5.基于无结构网格单元中心有限体积法的二维对流扩散方程离散 [J], 耿艳芬;王志力;陆永军
因版权原因,仅展示原文概要,查看原文内容请购买。
基于非结构网格有限体积法的Roe 隐格式初探杨彬河海大学,南京 (210098)E-mail :tuohuangzhe119@摘 要:本文采用非结构网格有限体积法,尝试建立适合河流、湖泊和近海浅水水流流动的Roe 隐式离散格式的数学模型。
对于隐式离散方程组采用线性化的处理方法以后,用对称的Gauss-Seidel(GS)法进行迭代求解,提高了求解效率。
鉴于地形的处理对数学模型稳定性和合理性的影响,将斜底模型的概念应用于模型的地形处理当中。
运用模型对几种标准算例进行验证模拟,并将模拟结果与相同条件下建立的Roe 显式离散格式模型的计算结果进行分析比较,通过计算表明本模型可以有效的缩短模拟计算时间,并且可以较好的处理间断解和各种复杂流态的过渡。
关键词:二维浅水方程,无结构网格,Roe 隐格式,对称GS 法,地形处理 中图分类号:TV131.21. 引言在实际计算不规则天然地形上的浅水流动时,为了贴切地模拟陆地边界,整个计算区域的网格尺寸常相差较大,在利用显式格式时,时间步长往往受到所有网格稳定性限制中最小者的控制,极大地降低了处理效率,特别是在计算历时较长的数值模拟中,计算机的实际耗时往往让人难以接受,从而使数值模拟受到限制。
鉴于隐式格式的高效率和较好的稳定性,以及近年来非结构网格技术的发展,有限体积法的广泛应用,许多高性能的隐式新算法在计算空气动力学领域中脱颖而出[1,2,3]。
然而在计算水力学中,目前利用无结构网格有限体积法,基本上是用的显式格式。
因此借鉴空气动力学中高性能的隐式格式,在浅水动力学中建立与本领域相适应的高性能的隐式格式不但有利于恒定流或解的时间变率小、流动过程历时长等一些实际问题的数值模拟,而且也能够有效的处理间断解。
本文尝试在二维浅水方程的基础上建立基于无结构网格有限体积法的Roe 隐式离散格式数学模型,利用此模型对几种标准算例进行验证模拟,并将计算结果与相同条件下建立的Roe 显式离散格式模型的计算结果进行分析比较。
2. 控制方程基于圣维南方程的基本假设,本模型采用浅水方程组[4]:S yGx E t U F t U =∂∂+∂∂+∂∂=⋅∇+∂∂ (1) [][][];2/,,;,2/,;,,2222TTTgh hv huv hu G huvgh hu hu E hv hu h U +=+==式中:=+=f S S S 0[][];,,0,,000TfyfxTyxghSghSghSghS−−+u ,v 分别为x ,y 方向的水深平均流速; t 为时间;h 为水深;g 为重力加速度;y x S S 00, 分别为x ,y 方向的底坡;fy fx S S ,分别为x ,y 方向的摩阻坡度,在本文中摩阻坡度可用曼宁公式计算,即:=fx S ,/3/4222h v u u n +3/4222/h v u v n S fy +=, n 为曼宁系数。
3. 空间离散对区域布置三角形网格,h 、u 、v 等变量定义在网格中心,采用网格中心格式(CC),控制体即为网格单元。
在每个三角形网格单元上利用有限体积法对浅水方程组进行积分,将方程组(1)写成积分形式:Ω=Ω⋅∇+Ω∂∂∫∫∫∫∫∫ΩΩΩSd d F d tU)( (2)其中:Ω为控制体积的区域。
运用高斯格林公式的原理将(2)式进一步离散成如下形式: ∫Ω=++∂∂i i iS A dl G E tUA )sin cos (θθ (3) 令θθsin cos G E F n +=,方程(3)转化为:∑=+=+∆−mj i i j n i n n i n ii S A L U F t U U A 11)( (4) 其中:i A 为控制体的面积,m 为控制体的边数,j L 为第j 条边的边长。
利用欧拉方程的旋转不变性,将二维问题转换为沿各边法向分别求解的一维问题: 定义旋转矩阵⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡−=θθθθcos sin 0sin cos 0001T ,其逆矩阵⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡−=−θθθθcos sin 0sin cos 00011T 由关系()()(U F TU F U TF n ==可以得到)()(1U F TU F n −=,代入方程(4)得:∑=−+=+∆−mj i i j n i n i i S A L U F T t U U A 111)( (5)其中:令U θθsin cos v u += ,V θθcos sin v u +−=,所以[]T hV hU h U ,,=,=)(U F[]ThUVgh hUhU ,2/,22+。
U ,V 分别表示单元边法向流速和切向流速。
对于(5)式,法向数值通量)(U F 的精确求解是问题的关键,对于一阶法向数值通量的求解本文采用了Roe 逆风格式。
其表达式为:{})(),(~)()(21),()(i j j i j i j i U U U U A U F U F U U F U F −⋅−+== (6) 其中:A ~是关于j i U U ,的某种意义下的Jacobian 平均矩阵,也称Roe 平均矩阵。
4. 时间离散将半离散化的方程(5)写成如下的隐式格式:∑=+−+=+∆−m j i i j n n i n i i S A L UF T t U U A 1111)( (7)将方程(7)进一步化为如下形式:{}j n mj i i j n n m j n i i L U F T S A L U F T U F T t U A )()()(111111∑∑=−−+=−−=−+∆∆ (8) 其中:n i n i n i U U U −=∆+1。
利用一阶近似,则[3]:{),(),(),()()(1111111111++−−++−−+−=−=−n j n i n j n i n j n i n n U U F T U U F T U U F T U F T U F T}{}),(),(),(11111n j n i n j n i n j n i U U F T U U F T U U F T −+−+−−+− (9)对式(9)的两项残值进行线性化,则:n j jn i i n n U U U F T U U U F T U F T U F T ∆∂∂+∆∂∂=−−−−+−)()()()(11111 (10) 对于式(6)在这里我们可以写成如下形式:{})()()(21),()(i j ij j i j i U U U F U F U U F U F −⋅−+==ρ (11) 其中:ij ρ为通量Jacobian 平均矩阵的谱半径[5]。
将式(11)代入式(10)中,则:{}n i ij i i n n U I U F T U F T UF T ∆⋅+∂∂=−−−+−ρ)(21)()(1111{}n j ij jj U I U F T ∆⋅−∂∂+−ρ)(211 (12)将式(12)代入式(8)中,并进行整理得到:{}i i m j j n mj n j j ij j j n i m j j ij i S A L U F T U L I U U F T U L t A −−∆⋅−∂∂−=∆+∆∑∑∑==−−=11111)()(21)21(ρρ (13) 因每一个控制体都是闭合的,所以:0)(2111=∂∂−=∑j ii m j L U F T 方程(13)中右边第二项,可以看作是显格式的净输出通量之和,仿照Roe 的黎曼间断问题的解法,并对间断处进行熵修正,将式(6)可以改写成如下形式[6]:{}∑=−+==31~)~(~)()(21),()(k k k k j i j i U F U F U U F U F γλψα (14)其中:kk γλ~~、分别为A ~的特征值、对应右特征向量;k 为特征值个数;ψ为熵修正函数。
本文利用(14)式计算显格式项的净输出通量之和。
5. 地形处理本模型中地形的处理采用斜底模型的概念,地面高程布置在各个节点上,控制体单元各边中点和单元中心的高程用线性插值得到,由于本模型采用三角形单元,所以节点都在同一平面上。
利用单元内部流量守恒确定单元各边法向流速和切向流速。
⎩⎨⎧−=−=c i c i i iZ h Z h ηη,所以⎩⎨⎧+−=+=ic c c i ic c c i h h v u V h h v u U /)cos sin (/)sin cos (θθθθ 其中:i η表示单元平均水位,c i Z Z ,分别表示单元边中点和单元中心的底高程,c i h h ,分别表示单元边中点和单元中心的水深,i i V U ,表示单元边中点处的法向和切向流速。
参考文[]7中推导的结论,将方程(2)中底坡项的积分推导成下面的形式:=Ω∫∫Ωd S 0Tm j j i j i m j L gh L gh ⎥⎦⎤⎢⎣⎡⋅⋅⋅⋅∑∑==1221sin 21,cos 21,0θθ,将此项代入到(13)中。
通过以上步骤将非线性隐式离散方程组化成线性方程组,并对底坡项进行处理,对于摩阻项采用迭代计算的最新值代入,以提高计算稳定性。
对各个隐元进行计算时,用求解线性方程组的Gauss-Seidel(GS)法,逐个隐元求解。
此时相邻控制体时段末的解设为已知,而该隐元的解又用于后面的隐元计算中。
为了消除控制体编号顺序的影响,这里采用对称扫描的方法:先按编号从小到大的顺序计算一次,再按编号从大到小的顺序计算一次。
按照上面的对称GS法,对隐元进行反复迭代计算,直到收敛。
6.验证计算与讨论6.1圆形溃坝一个边长为50m的正方形区域,在区域中心有一半径为11m,水深10m的圆形坝,外区域水深2m,假定坝体在某一时刻瞬时溃决。
计算网格采用三角形网格,t∆为0.1S,x∆为1m,单元网格数为5642个。
图1,图2为本模型计算溃坝后0.69S的自由水面分布图和等值线图。
与前人的模拟结果相吻合[7],验证了本模型在处理对称间断流方面的能力。
h/mx/my/m图 1 圆形溃坝计算水面分布图(t=0.69S) 图 2 圆形溃坝水位等值线图(t=0.69S)Fig. 1Computed water surface of the circular Fig. 2Water depth contours of the circular dam-break ( t = 0. 69 s) dam-break ( t = 0. 69 s)6.2二维矩形部分溃坝一长、宽各200m矩形计算区域中间设有一坝,,坝上游水位10m,下游水位5m,某一时刻,堤坝在y为95到170m之间瞬时溃决。
计算网格采用三角形网格,t∆为0.1S,x∆为5m,网格单元数3758个。