耦合VOF的非线性k-ε模型三维溃坝数值模拟

第３８卷第５期　２０１７年１Ｏ月　水道港　口　

Ｊｏｕｒｎａｌ　ｏｆ　Ｗａｔｅｒｗａｙ　ａｎｄ　Ｈａｒｂｏｒ　ＶｏＩ．３８　Ｎｏ．５　

Ｏｅｔ．２０１７　

耦合Ｖ０Ｆ的非线性　一　模型三维溃坝数值模拟　丁伟业，金生，艾丛芳　（大连理工大学海岸和近海工程国家重点实验室，大连１１６０２４）　

摘要：溃坝是一种伴随有巨大破坏力的灾害性的水流现象。三维溃坝的数值模拟能够为溃坝水流的　研究分析提供科学的依据。针对瞬间全溃模型水流的复杂性，考虑水库下游障碍物的滞水及阻水作　用，利用开源程序ＯｐｅｎＦＯＡＭ建立了耦合ＶＯＦ对流方程的三维非线性（ＮＬ）　一　紊流数学模型求解雷　诺时均纳维斯托克斯方程。采用有限体积法进行离散，利用ＰＩＳＯ算法进行数值求解。通过模拟下游为　三角形障碍物、矩形障碍物以及９０。弯道的溃坝实验，并将模拟结果与实验结果进行对比分析，初步证　明了该模型的准确性和鲁棒性。　关键词：ＶＯＦ法；ＮＬｋ—ｓ模型；溃坝；ＯｐｅｎＦＯＡＭ；三维数值模拟　中图分类号：ＴＶ　１２２　文献标识码：Ａ　文章编号：１００５—８４４３（２０１７）０５—０４８９—０６　

２０１６年１月１４日１１时许，湖北省恩施龙凤镇三河村白庙突发堰塞湖溃坝事故。２０１６年２月３日，伊　拉克摩苏尔水坝面临即将溃坝的风险，成千上万的人处在危险之中，工作人员对水坝进行加固。溃坝是指　坝体溃决，是一种低频率高危害的社会致灾因素，会对下游地区人民生命财产安全造成灾难性的破坏。因　此进行溃坝洪水数值计算与模拟分析具有重要的理论与现实意义¨　。研究溃坝问题的方法主要有三种：理　论分析、模型试验和数值模拟。三种研究方法各有其优缺点：当处理十分简单的模型问题时理论分析方法　能够得到较为精确的解析解，而对于复杂问题并不能得到其控制方程的解析解；模型试验方法是研究水流　流动机理、分析水流流动现象、探讨并获得水流流动新概念的主要手段，然而模型试验受经费、时间及观测　精度的限制，并且难以摆脱原型和模型之间不完全相似的困扰；数值模拟方法能够解决复杂的水流流动问　题，同时不受经费和时间限制　。　近年来随着计算技术的不断发展，数值模拟越来越多的应用于溃坝以及溃坝水流对下游复杂地形冲刷　的研究当中。虽然１一Ｄ和水深２一Ｄ溃坝水流模型能够模拟溃坝水流的演进过程，但是由于对控制方程进　行了如静水压强分布假设，以及假设无明显的垂向加速度和自由表面曲率等，使得这些模型在溃坝的初始　阶段、溃坝波前端和靠近结构物水流等方面均不适用了　ｊ。垂向２Ｄ和３Ｄ模型能够更精细地模拟溃坝水　流的复杂流态及细节，尤其对于自由水面变化强烈的溃坝水流，３Ｄ模型所得到的丰富信息为研究局部水流　特性、溃口发展机理以及水工建筑物结构设计等提供了依据＿５　Ｊ。对于自由水面的追踪捕捉，Ｌｕｃｙ　ＬＢ＿６］，　Ｇｉｎｇｏｌｄ和Ｍｏｎａｇｈａｎ　提出了利用流体粒子运动代替计算网格来求解控制方程的光滑粒子运动方法　（ＳＰＨ）。近些年三维非静压数值模型在自由表面的问题研究中具有较快的发展，从求解控制丰城数值方法　分类主要有：显示投影法　；半隐分步法　；全隐式法¨　。由Ｈｉｒｔ和Ｎｉｃｈｏｌｓ口纠提出的ＶＯＦ方法则更为研　究者们所熟知。体积分数Ｆ介于０和１之间，当Ｆ等于０或１时，则单元内没有交界面。本文利用开源计算　软件ＯｐｅｎＦＯＡＭ研究３Ｄ溃坝流模型。利用ＶＯＦ方法耦合ＮＬｋ一　紊流模型求解雷诺时均Ｎ—ｓ方程。采　用六面体网格单元建立网格模型，有限体积法进行离散，ＰＩＳＯ算法进行求解。每一个ＰＩＳＯ算法的内迭代步　

收稿日期：２０１７—０３—２７；修回日期：２０１７—０５—０５　基金项目：国家自然科学基金项目（５１３０９０５２；５１４７９０２２）；中央高校基本科研业务费专项基金（ＤＵＴ１５ＬＫ０１）；辽宁省教育　厅重点实验室基础研究项Ｅｌ（ＬＺ２０１５０１２）．　作者简介：丁伟业（１９８８一），男，黑龙江省人，博士研究生，主要从事流体自由表面研究。　Ｂｉｏｇｒａｐｈｙ：ＤＩＮＧ　Ｗｅｉ—ｙｅ（１９８８一），ｍａｌｅ，ｄｏｃｔｏｒ　ｓｔｕｄｅｎｔ．　４９０　水道港　口　第３８卷第５期　中的压力梯度通过压力泊松方程来求解，并用Ｒｈｉｅ—Ｃｈｏｗ插值原理来消除压力波　。通过与实验案例进　行对比验证来证明模型的可靠性。　

１　数学模型　溃坝波的运动主要是受重力作用，水流和气体有明显的分界面，因此可以用水气两相流的分层模型　ＶＯＦ　来进行模拟自由液面。自由水面的运动通过流体体积函数的对流输运方程来捕捉，可以较为精细　地描述溃坝洪水中的水面变化，克服了静压假定和刚盖假定对变化剧烈的自由水面的限制和导致的压力场　失真。　ＶＯＦ方法是将每个单元中一项流体体积与单元体积之比定义为流体体积函数，即　ｆ１　液体　（　，ｔ）＝｛０＜Ｏｔ＜１　气液交界面　（１）　【０　气体　并且单元内的物理特性为　Ｐ＝１９１　＋（１一ａ）ｐ　（２）　＝　１　＋（１一　）　（３）　其中ｚ与ｇ分别为气相和液相。则流体体积分数　满足对流输运方程　

＋Ｖ　“　］＋Ｖ［ｕ，　（１一　）］＝０　（４）　

且　满足控制条件　～＿－１，　为计算流体总数。式（４）中的Ｖ［ＵｒＯ／（１一　）］为限制交界面人工压缩项，　

仅适用于薄交界面区域，在此区域外的影响很小。　为可压缩交界面流体相对流速。　模型的基本微分方程包括连续性方程、动量方程、以及＾　一　方程，分别表示如下。　连续性方程　Ｖ・　＝０　（５）　雷诺时均动量方程　

＋Ｖ・［ｐｕｔｔ　］＝一Ｖｐ　一ｇ・　Ｖｐ＋Ｖ・［肛Ｖ　ｕ＋　舭］＋ｏｒｒｋ　Ｖ　（６）　式中：　＝（　，　，Ｗ）为笛卡尔坐标系下流速；ｐ　＝ｐ—Ｐｇ・　为修正压力；Ｐ＝ｐ（　）　＝　（　）由式（２）、式（３）得　到；　和Ｉｉ｝　＝７（　）分别为气液交界面处的表面张力和曲率。ｒ　为非线性雷诺应力张量　

ＴＮＬ＝ｑ　－（　＋　）　（　／）＝　十戈　Ｊ　南｛ｃ　［７　ｕ。Ｖ　＋（Ｖ　Ｍ‘Ｖ　）　］＋ｃ祀［Ｖ　（Ｖ“）　】＋ｃ　【（Ｖ　Ｍ）　・Ｖ　】）（８）　式中：Ｃ　１＝一４．０，Ｃｒ２＝１３．０，Ｇ　＝－２．０和Ａ２＝１　０００均为常数。ｒ／＝ＳＫ／ｅ由Ｙａｋｈｏｔ　ｅｔ　ａ１．（１９９２）获得，且　ｓ为流体的平均应变率，有　Ｓ＝　丽，Ｓ＝÷（ｖ　＋（ｖ　）　）　（９）　

ｋ方程　警＋Ｖ’（ｐｕｋ）＝Ｖ【Ｐ（　＋丝ｏｒｋ１／Ｖ｜ｉ｝】＋ｐＰ　（１ｏ）　

方程　警＋７　＝Ｖ［Ｐ（　＋　）Ｖ小ｐＣ・　ｚ譬　…）　丁伟业，等耦合ＶＯＦ的非线性　一　模型三维溃坝数值模拟４９１　其中　譬　

丽　式中：　＝０．９和Ａ．＝１．２５为常数，且　＝￣２ｋ／ｅ，且　

＝　，　＝　

１（Ｖ　ｚｚ一（Ｖ　ｕ）　）　

ＮＬｋ—　中生成项　同样由非线性描述表达，有　Ｐ　＝ｐ（　，Ｓ：Ｖ“一ｒ舭：Ｖ　Ｉｔ）　式Ｉ１ｌ：经验常数取值为Ｃ，＝１．４４，Ｃ２＝１．９２，　＝０．７７，Ｏ－　＝１　

２模型验证　２．Ｉ　下游为三角形障碍物溃坝模拟　模型为长３８　ｍ、宽１　ｍ的矩形河道。河道上游为长１５．５　ＨＩ，初始水深　为０．７５　ｍ的水库。障碍物位于河道起点下游２５．５　１１１位置，障碍物与水库　问为干燥的河床，障碍物以下河床水深为０．１５　ｒｎ，具体尺寸如图１所示。　图中测点Ｇ４、Ｇ１０、Ｇ１３、Ｇ２０分别位于河道中部距离水库下游４　ｍ、１０　ｍ、　１３　Ｉｌｌ和２０　Ｉｌ　’　。对于水平网格：障碍物上游选择３　ｅｍ，障碍物附近选　择ｌ　Ｃｌｌ，障碍物下游选择２　ｅｌｌ。在宽和高方向均采用１　ｃｉｎ网格。　图２为各测点测得水位实验值与计算值的　。。　对比结果。由图可知，Ｃ４点整体吻合较好，但是　ｏ　ｅ　汁算模拟的反射波与实验值相比略微滞后；Ｇ１０　：：　点计算模拟的最大水位略高于测量水深；Ｇ１３点　：　计算模拟上游溃坝来流时间与实验结果相比稍　有提前，且在反射波过后的ｔ＝３０～３４　ｓ，计算值　：　与测量值差异较大；计算模拟溃坝水流抵达（；２０　点时有冲刷，【大１而计算模拟的水位低于实验结　暑ｏ　ｓ　果。虽然各测点的模拟值与实验值均存在差　异，但模拟值　实验测量值整体吻合还是较　好的。　２．２　下游为矩形障碍物溃坝模拟　模型为长３．２２　Ｉｌｌ、宽１　ｎｌ的矩形河道，河道　上游水库长１．２３　Ｉｌ，水库初始水深ｈ　为０．５５　

（１２）　（１３）　

（１４）　（１５）　

ｕｎ１　Ｌ：ｎｌ　图１模型三维尺寸图　Ｆｉｇ．１　Ｔｈｒｅｅ—ｄｉｍｅｎｓｉｏｎａｌ　ｓｉｚｅ　ｏｆ　ｍｏｄｅｌ　

ｍ，下游河床干燥无水。矩形障碍物尺寸为０．４０３×０．１６１×０．１６１，位于水库下游１．１７　ｍ。图中Ｈ４和Ｈ２为　水位测点，ｐｌ、ｐ３、ｐ５和ｐ７分别为压力测点　。　。采用２　ｃ・ｎｌ×２　ｔｉｌ×ｌ　ＯＩｌ网格划分模型。　采用　一　和ＮＬｋ一　两种紊流方程对模型２进行数值模拟。图３为　Ｈ４与Ｈ２水位测点数值模拟与实验测量和Ｋ．Ｍ．Ｔ．Ｋｌｅｅｆｓｍａｎ（２００５）　结果的对比。图４为ｔ＝０．５　Ｓ，０．８　Ｓ，１．４　Ｓ时的流域特征量。图４一ａ为　流速矢量压力云图，图４一ｂ为相应时刻的河道中心剖面体积分数图，图４　

一ｃ为３Ｄ流场示图。ｔ＝０．５　Ｓ时，溃坝波抵达障碍物，水流由中心向两侧　及Ｌ方绕流经过障碍物，因此在障碍物前端形成高压　，Ｈ２点测得的水　深为０．０８　ｎ　。ｔ＝０．８　ｓ时，溃坝波仍然作用于障碍物上，障碍物后方形成　

．　１　２３　一　Ｉ　ｌ７一Ｏ　．　ｕｎｉｔ：ｍ　、ｇ姒　

图３模型三维尺寸图　Ｆｉｇ．３　Ｔｈｒｅｅ—ｄｉｍｅｎｓｉｏｎａｌ　ｓｉｚｅ　ｏｆ　ｍｏｄｅｌ