三维自由面湍流场数值模拟及其在水利工程中的应用

  本文关键词：三维自由面湍流场数值模拟及其在水利工程中的应用，由笔耕文化传播整理发布。

河海大学 博士学位论文 三维自由面湍流场数值模拟及其在水利工程中的应用 姓名：王志东 申请学位级别：博士 专业：水力学及河流动力学 指导教师：汪德爟 20040322

河海大学博士论文

摘要
自然界中的流体流动问题一般都具有三维性和湍动性，对于最常见的流体一水则
同时还拥有自由表面，水利工程中绕闸墩的泄洪水流、船

舶工程中的船体绕流等均具 有上述流动的三种特性，这三种特性是计算流体力学（ＣＦＤ）领域非常重要的研究内 容和方向。
．

本文针对粘流场数值模拟技术的研究现状和发展动态作了比较全面的回顾与展 望，重点研究了网格生成技术、数值求解技术、湍流模型技术和动边界模拟技术，在 此基础上建立了模拟自由面湍流场的数学模型，并成功地应用Ｎ－维溢流坝和带闸墩
的三维溢流坝过坝水流的数值计算。

本文完成的具体工作和主要结论如下： １）以代数网格生成方法为基础提出了一种简单的、可独立于网格生成方法之外的边 界正交化技术；针对分区结构网格系统建立了分区交界面处的数据结构与计算模型： ２）利用有限体积方法在非正交同位网格系统中建立了ＳＩＭＰＬＥ求解算法，对非正交 网格系统中的边界条件、延迟修正技术及计算节点的梯度计算等专题进行了深入讨 论。通过对倾斜方腔驱动流、后台阶绕流、圆柱绕流等典型流场的数值模拟，表明本 文建立的非正交同位网格系统中的ＳＩＭＰＬＥ算法是合理的； ３）研究比较了各种湍流模型的特点及应用范围，以水翼绕流场为例完成了ＲＮＧｋ一￡ 模型与标准ｋ—ｓ模型的计算比较，结果表明ＲＮＧｋ—ｓ模型能够更好地模拟流场的湍
流特性：

４）研究了运动边界模拟技术中的ＶＯＦ方法，详细建立了一种精度更高的自由面重构 模型，通过对典型运动界面的数值模拟表明，该方法能够更加准确地确定运动界面的
形状和位置：

５）利用本文建立的非正交网格系统中自由面湍流场数学模型，对二维溢流坝过坝水 流进行了数值计算，分析了不同坝顶水头下的水面线位置、坝面压力及流量、流速等 水流特性，通过与典型实验资料的对比，水面线高度及坝面负压极值均与实验值具有
非常好的一致性。 ６）通过建立溢流坝和闸墩的三维整体模型，完成了对三维过坝水流粘流场的数值计

算，针对闸墩的型式及在坝面的布置计算了闸墩对水面线形状、坝面压力、流量系数 等设计参数的影响，并将不同墩型与布置形式的计算结果进行比较。结果表明，闸墩

ＩＩ

河海大学博士论文

的存在抬高了水面线的位置，其中在闸墩头部尤其明显：墩型对流量系数和坝面压力 影响较大．３Ａ型闸墩相对于２型闸墩具有更大的流量系数和更小的坝面压力：不同
的墩厚闸宽比ｆ／６对泄流能力也将产生显著的影响．随着墩厚闸宽比的增加，坝面压

力降低，而当ｔｌｂ＝０．２时溢流坝具有更大的流量系数。 ７）本文建立的三维自由面湍流场模型能够比较准确地模拟带闸墩溢流坝三维过坝水 流湍流场，可以针对不同的坝型、墩型、坝面布置形式以及不同的坝顶水头完成系列
水力计算，为泄水工程建筑物的设计提供可靠的分析依据。

关键词：数值模拟，网格，湍流模型，流体体积法，溢流坝，闸墩

ＩＩ

河海大学博士论文

ＡＢＳＴＲＡＣＴ
Ｔｈｒｅｅ ｄｉｍｅｎｓｉｏｎ ａｎｄ ｔｕｒｂｕｌｅｎｔ
ａｒｅ

ｃｏｍｍｏｎ ｃｈａｒａｃｔｅｒｉｓｔｉｃ ｏｆ ｆｌｕｉｄ ｆｌｏｗ

ｐｒｏｂｌｅｍ ｉｎ ｓｕｃｈ ｆｌｏｗ
ａｓ

ｎａｔｕｒｅ

ａｎｄ

ａｄｄｅｄ ｆｌｅｅ ｓｕｒｆａｃｅ ｆｏｒ ｔｈｅ ｂｅｓｔ ｆａｍｉｌｉａｒ ｆｌｕｉｄ－ｗａｔｅｒ，

ｆｌｏｏｄ ｄｉｓｃｈａｒｇｅ ｆｌｏｗ ａｒｏｕｎｄ ｆｒｕｓｔａ ｏｆ ｂｒａｋｅ ｉｎ ｗａｔｅｒｐｏｗｅｒ ｉｎ ｓｈｉｐ ｅｎｇｉｎｅｅｒｉｎｇ，ａｎｄ
ＳＯ

ｐｒｏｊｅｃｔ，

ａｒｏｕｎｄ ｈｕｌｌ

ｏｎ，ａｌｌ

ｈａｖｅ ｔｈｅ ｔｈｒｅｅ

ｃｈａｒａｃｔｅｒｉｓｔｉｃｓ ｗｈｉｃｈ ａｒｅ ｔｈｅ ｖｅｒｙ ｉｍｐｏｒｔａｎｔ ｓｔｕｄｙ ｃｏｎｔｅｎｔ ａｎｄ ａｓｐｅｃｔｓ ｉｎ
ＣＦＤ．
Ａ ｃｏｍｐｒｅｈｅｎｓｉｖｅ ｒｅｖｉｅｗ ａｎｄ ｅｘｐｅｃｔａｔｉｏｎ ａｂｏｕｔ ｒｅｓｅａｒｃｈｉｎｇ ｓｔａｔｕｓ ｉｎ ｑｕｏ

ａｎｄ ｄｅｖｅｌｏｐｍｅｎｔ ｔｒｅｎｄｓ ｏｆ ｖｉｓｃｏｕｓ ｆｌｏｗ ｎｕｍｅｒｉｃａｌ ｓｉｍｕｌａｔｉｏｎ ｔｅｃｈｎｉｑｕｅ ａｒｅ ｇｉｖｅｎ ｉｎ ｔｈｉｓ ｐａｐｅｒ

ａｎｄ ｇｒｉｄｓ ｇｅｎｅｒａｔｉｏｎ

ｍｅｔｈｏｄｓ，ｎｕｍｅｒｉｃａｌ
ｊｏｂ
ａ

ｃｏｍｐｕｔａｔｉｏｎ
ａｒｅ

ｍｅｔｈｏｄｓ，ｔｕｒｂｕｌｅｎｔ ｍｏｄｅｌ ａｎｄ ｍｏｖｉｎｇ ｂｏｕｎｄａｒｙ ｓｉｍｕｌａｔｉｏｎ ｍｅｔｈｏｄｓ ｅｍｐｈａｓｉｓ ｓｔｕｄｙ ｉｎ ｔｈｅ ｐａｐｅｒ．Ｏｎ ｔｈｅ ｂａｓｉｃ ｏｆ ａｂｏｖｅ

ｔｈｅ

ｍａｔｈｅｍａｔｉｃａｌ ｍｏｄｅｌ

ｏｆ ｓｉｍｕｌａｔｉｏｎ ｔｕｒｂｕｌｅｎｔ ｆｌｏｗ ｗｉｔｈ ｆｒｅｅ ｓｕｒｆａｃｅ ｉｓ ｅｓｔａｂｌｉｓｈｅｄ ａｎｄ ａｐｐｌｉｅｄ
ｔｏ

ｔｒｉｕｍｐｈａｎｔｌｙ

ｎｕｍｅｒｉｃａｌ ｃｏｍｐｕｔａｔｉｏｎ

ｏｎ

２Ｄ ｓｐｉｌｌｗａｙ

ａｎｄ

３Ｄ ｓｐｉｌｌｗａｙ ｗｉｔｈ ｆｒｕｓｔａ

ｏｆ ｂｒａｋｅ． Ｔｈｅ ｍａｔｅｒｉａｌ

ｊｏｂ ａｎｄ ｍａｉｎ

ｃｏｎｃｌｕｓｉｏｎ

ａｒｅ

ａｓ ｆｏｌｌｏｗｓ：

（１）Ａ ｓｉｍｐｌｅ ｂｏｕｎｄａｒｙ ｏｒｔｈｏｇｏｎａｌｉｚａｔｉｏｎ ｐｒｏｃｅｄｕｒｅ ｉｎｄｅｐｅｎｄｅｎｃｅ ｇｒｉｄｓ
ｇｅｎｅｒａｔｉｏｎ ｍｅｔｈｏｄ ｉｓ ｐｕｔ ｆｏｒｗａｒｄ
ｏｎ

ｔｈｅ ｂａｓｅ ｏｆ

ａｌｇｅｂｒａｉｃ鲥ｄ

ｇｅｎｅｒａｔｉｏｎ

ｍｅｔｈｏｄ；Ｄａｔａ
ａｉｍｉｎｇ
ａｔ

ｓｔｒｕｃｔｕｒｅ

ａｎｄ

ｃｏｍｐｕｔａｔｉｏｎａｌ ｍｏｄｅｌ

ｏｎ

ｉｎｔｅｒｆａｃｅ ａｒｅ ｅｓｔａｂｌｉｓｈｅｄ

ｂｌｏｃｋｓ ｏｆ ｓｔｒｕｃｔｕｒｅｄ ｇｒｉｄｓ． ｐｒｏｃｅｄｕｒｅ ｉｓ ｅｓｔａｂｌｉｓｈｅｄ ｉｎ ｎｏｎｏｒｔｈｏｇｏｎａｌ ｇｒｉｄｓ ｂｙ
ｏｎ
ｕｓｅ

（２）ＳＩＭＰＬＥ

ｏｆ
ａｓ

ｔｈｅ ｆｉｎｉｔｅ ｖｏｌｕｍｅ ｍｅｔｈｏｄ ａｎｄ ｉｎ－ｄｅｐｔｈ ｄｉｓｃｕｓｓｅｓ

ｓｐｅｃｉａｌ ｔｏｐｉｃ ｓｕｃｈ

ｂｏｕｎｄａｒｙ ｃｏｎｄｉｔｉｏｎ ｉｎ ｎｏｎｏｎｈｏｇｏｎａｌ ｇｒｉｄｓ，ｄｅｆｅｒｒｅｄ ｃｏｒｒｅｃｔｉｏｎ ｍｅｔｈｏｄ ｇｒａｄｓ ｃｏｍｐｕｔｅ
ｏｎ

ａｎｄ

ｃａｌｃｕｌａｔｉｏｎａｌ

ｎｏｄｅｓ，ａｎｄ

ＳＯ

ｏｎ．Ｔｈｅ

ＳＩＭＰＬＥ ｐｒｏｃｅｄｕｒｅ ｉｎ

ｎｏｎｏｒｔｈｏｇｏｎａｌ ｇｒｉｄｓ ｅｓｔａｂｌｉｓｈｅｄ ｉｎ ｔｈｅ ｐａｐｅｒ ｉｓ ｐｒｏｖｅｄ ｔｏ ｂｅ ｒｅａｓｏｎａｂｌｅ ｔｈｒｏｕｇｈ ｒｅｐｒｅｓｅｎｔａｔｉｖｅ ｎｕｍｅｒｉｃａｌ ｓｉｍｕｌａｔｉｏｎ ｓｕｃｈ ａｓ ｉｎｃｌｉｎｅ ｓｑｕａｒｅ ｃａｖｉｔｙ ｄｒｉｖｉｎｇ ｆｌｏｗ，ｆｌｏｗ ａｒｏｕｎｄ ｂａｃｋ ｓｉｄｅｓｔｅｐ，ｆｌｏｗ ａｒｏｕｎｄ ｃｙｌｉｎｄｅｒ．

（３）Ｃｈａｒａｃｔｅｒｉｓｔｉｃ ａｎｄ ｓｔｕｄｉｅｄ ａｎｄ ｃｏｍｐａｒｅｄ ａｎｄ

ａｐｐｌｉｃａｔｉｏｎ ｒａｎｇｅ ｏｆ ｄｉｆｆｅｒｅｎｔ ｔｕｒｂｕｌｅｎｔ ｍｏｄｅｌ ａｒｅ ｃｏｍｐｕｔｅ ｃｏｍｐａｒｅ ＲＮＧｋ－ｃ ｔｕｒｂｕｌｅｎｔ ｍｏｄｅｌ ｗｉｔｈ
ｃａｓｅ

ｋ－占ｍｏｄｅｌ ｉｓ ｅｓｔａｂｌｉｓｈｅｄ ｂｙ ｔｈｅ

ｏｆ ｖｉｓｃｏｕｓ ｆｌｏｗ ｔｈｒｏｕｇｈ ｈｙｄｒｏｆｏｉｌ
ｃａｌｌ

ａｎｄ

ｔ１１ｅ ｒｅｓｕｌｔ ｓｈｏｗ ｔｈａｔ ｂｅＲｅｒ ｃｈａｒａｃｔｅｒｉｓｔｉｃ ｔｕｒｂｕｌｅｎｔ ｆｌｏｗ
ＲＮＧｋ一占．

ｂｅ ｏｂｔａｉｎｅｄ ｂｙ

（４）ＶＯＦ ｍｅｔｈｏｄ ｉｎ ｍｏｖｉｎｇ ｂｏｕｎｄａｒｙ ｓｉｍｕｌｍｉｏｎ ｍｅｔｈｏｄ
ｉｖ

ｉｓ ｓｔｕｄｉｅｄ ａｎｄ

ａ

河海大学博士论文

ｈｉｇｈｅｒ ｐｒｅｃｉｓｉｏｎ ｆｒｅｅ ｓｕｒｆａｃｅ ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ ｍｏｄｅｌ ｉｓ ｅｓｔａｂｌｉｓｈｅｄ ｄｅｔａｉｌｅｄｌｙ． Ｔｈｅ ｍｅｔｈｏｄ ｉｓ ｐｒｏｖｅｄ ｍｏｔｉｏｎ ｉｎｔｅｒｆａｃｅ ｉｎｔｅｒｆａｃｅ．
ｔｏ

ｗｅｌｌ ａｎｄ ｔｒｕｌｙ ｃｏｎｆｉｒｍ ｔｈｅ ｓｈａｐｅ ａｎｄ ｌｏｃａｔｉｏｎ ｏｆ ｎｕｍｅｒｉｃａｌ ｓｉｍｕｌａｔｉｏｎ ｏｆ ｒｅｐｒｅｓｅｎｔａｔｉｖｅ ｍｏｔｉｏｎ

ｔｈｒｏｕｇｈ

（５）Ｃｏｍｐｕｔｅ ｏｆ ｆｌｏｗ
ｆｒｅｅ ｓｕｒｆａｃｅ、ｐｒｅｓｓｕｒｅ

ｏｖｅｒ

２Ｄ ｓｐｉｌｌｗａｙ ｉｓ ｄｏｎｅ ｂｙ

ＢＳｅ

ｏｆ

ｅｓｔａｂｌｉｓｈｅｄ ａｎｄ

ｆｒｅｅ

ｓｕｒｆａｃｅ ｔｕｒｂｕｌｅｎｔ ｍａｔｈｅｍａｔｉｃａｌ ｍｏｄｅｌ ｉｎ ｎｏｎｏｒｔｈｏｇｏｎａｌ ｇｒｉｄｓ
ｏｎ

ｌｏｃａｔｉｏｎ ｏｆ
ａｓ

ｄａｍ

ａｎｄ

ｆｌｏｗ ｃｈａｒａｃｔｅｒｉｓｔｉｃ

ｏｎ

ｔｈｅ ｓｕｃｈ

ｆｌｕｘ、 ｔｈｅ

ｖｅｌｏｃｉｔｙ ｏｆ ｆｌｏｗ ｏｆ ｄｉｆｆｅｒｅｎｔ ｗａｔｅｒ ｈｅｉｇｈｔ ｏｆ ｄａｍ ｐｅａｋ

ａｒｅ

ａｎａｌｙｚｅｄ ａｎｄ
ｏｎ

ｒｅｓｕｌｔ ｏｆ ｈｅｉｇｈｔ ｏｆ ｆｒｅｅ ｓｕｒｆａｃｅ ａｎｄ ｎｅｇａｔｉｖｅ ｅｘｔｒｅｍｕｍ ｏｆ ｐｒｅｓｓｕｒｅ ｓｕｒｆａｃｅ ａｒｅ ｗｅｌｌ ｉｎｏｓｃｕｌａｔｅ
ｔｏ

ｄａｍ

ｒｅｓｕｌｔ ｏｆ ｒｅｐｒｅｓｅｎｔａｔｉｖｅ ｅｘｐｅｒｉｍｅｎｔ．

（６）３Ｄ

ｉｎｔｅｇｅｒ ｍｏｄｅｌ ｏｆ ｓｐｉｌｌｗａｙ ａｎｄ ｆｒｕｓｔａ ｏｆ
ｏｖｅｒ

ｂｒａｋｅ ａｒｅ ｅｓｔａｂｌｉｓｈｅｄ ａｎｄ

ｎｕｍｅｒｉｃａｌ ｃｏｍｐｕｔａｔｉｏｎ ｏｆ ３Ｄ ｖｉｓｃｏｕｓ ｆｌｏｗ

ｓｐｉｌｌｗａｙ ｉｓ

ｃｏｍｐｌｅｔｅｄ．Ｔｈｅ
ｓｈａｐｅ ｏｆ ｆｒｅｅ

ｉｎｆｌｕｅｎｃｅ ｏｆ ｆｒｕｓｔａ ｏｆ ｂｒａｋｅ ｆｏｒ ｄｅｓｉｇｎ ｐａｒａｍｅｔｅｒｓ ｓｕｃｈ ｓｕｒｆａｃｅ，ｐｒｅｓｓｕｒｅ
ｏｎ

ａＳ

ｄａｍ，ｆｌｕｘ ｃｏｅｆｆｉｃｉｅｎｔ ｉｓ ｃｏｍｐｕｔｅｄ ａｉｍ ａｔ ｔｙｐｅ ｏｆ ｆｒｕｓｔａ ｏｆ
ｏｎ

ｂｒａｋｅ ａｎｄ

ｄｉｓｐｏｓａｌ

ｄａｍ

ａｎｄ ｃｏｍｐｕｔｅ ｒｅｓｕｌｔ ｏｆ ｄｉｆｆｅｒｅｎｔ ｔｙｐｅ ｏｆ ｆｒｕｓｔａ ａｎｄ
ｏｆ ｏｆ

ｄｉｓｐｏｓａｌ ｍｏｄｅｌ ｉｓ ｃｏｍｐａｒｅｄ．Ｔｈｅ ｒｅｓｕｌｔ ｓｈｏｗ ｔｈａｔ ｔｈｅ ｐｏｓｉｔｉｏｎ ｏｆ ｆｒｅｅ ｓｕｒｆａｃｅ ｉｓ

ｈｉｇｈｅｒ

ｂｅｃａｕｓｅ ｏｆ ｆｒｕｓｔａ ｏｆ

ｂｒａｋｅ；ｆｌｕｘ ｃｏｅｆｆｉｃｉｅｎｔ ｆｒｕｓｔａ

ｂｒａｋｅ，ｅｓｐｅｃｉａｌ ｏｂｖｉｏｕｓ ｉｎ ｔｈｅ ｈｅａｄ ｏｆ ｆｒｕｓｔａ ａｎｄ ｐｒｅｓｓｕｒｅ ｏｎ ｄａｍ ａｒｅ ｑｕｉｔｅ ｉｎｆｌｕｅｎｃｅｄ ｂｙ ｔｙｐｅ

ａｎｄ ３Ａ ｔｙｐｅ ｆｒｕｓｔａ ｏｆ ｂｒａｋｅ ｒｅｌａｔｉｖｅ ｔｏ ２ ｔｙｐｅ ｆｒｕｓｔａ ｏｆ ｂｒａｋｅ ｈａＳ ｂｉｇｇｅｒ ｆｌｕｘ ｃｏｅｆｆｉｃｉｅｎｔ ａｎｄ ｓｍａｌｌｅｒ ｐｒｅｓｓｕｒｅ ｏｎ ｄａｍ；Ｏｂｖｉｏｕｓ ｉｎｆｌｕｅｎｃｅ ｏｆ ａｂｉｌｉｔｙ ｏｆ
ｄｉｓｃｈａｒｇｅ ｆｌｏｗ ｆｏｒ ｄｉｆｆｅｒｅｎｔ ｒａｔｉｏ ｏｆ ｔｈｉｃｋｎｅｓｓ ｏｆ ｆｒｕｓｔａ

ａｎｄ

ｂｒｅａｄｔｈ ｏｆ ｂｒａｋｅ

ａｎｄ

ｐｒｅｓｓｕｒｅ

ｏｎ

ｄａｍ ｒｅｄｕｃｅｓ ｗｉｔｈ ｉｎｃｒｅａｓｉｎｇ ｒａｔｉｏｎ ｏｆ ｔｈｉｃｋｎｅｓｓ ｏｆ ｆｒｕｓｔａ

ａｎｄ

ｂｒｅａｄｔｈ ｏｆｂｒａｋｅ，ｗｈｉｌｅ ｂｉｇｇｅｒ ｆｌｕｘ ｃｏｅｆｆｉｃｉｅｎｔ ｉｓ ｏｂｔａｉｎｅｄ ｗｈｅｎ ｔ／ｂ＝０．２．

（７）Ｔｈｅ
ｂｒａｋｅ

３Ｄ ｆｒｅｅ ｓｕｒｆａｃｅ ｔｕｒｂｕｌｅｎｔ ｍｏｄｅｌ ｅｓｔａｂｌｉｓｈｅｄ ｉｎ ｔｈｅ ｐａｐｅｒ ｉｓ

ｃｏｍｐａｒａｔｉｖｅｌｙ ｗｅｌｌ

ａｎｄ

ｔｒｕｌｙ ｓｉｍｕｌａｔｅｓ ３Ｄ ｆｌｏｗ

ｏｖｅｒ

ｓｐｉｌｌｗａｙ ｗｉｔｈ ｆｒｕｓｔａ ｏｆ

ａｎｄ ｓｅｒｉｅｓ ｗａｔｅｒｐｏｗｅｒ ｃｏｍｐｕｔａｔｉｏｎ ｃａｎ ｂｅ ｃｏｍｐｌｅｔｅｄ ａｉｍ ａｔ ｄｉｆｆｅｒｅｎｔ ｔｙｐｅ ｏｆ ｄａｍ、ｔｙｐｅ ｏｆ ｆｒｕｓｔａ ａｎｄ ｄｉｓｐｏｓａｌ ｏｎ ｔｈｅ ｄａｍ ａｎｄ ｄｉｆｆｅｒｅｎｔ ｗａｔｅｒ ｈｅｉｇｈｔ
ｏｎ

ｐｅａｌ（ｏｆ ｄａｍ，ｗｈｉｃｈ ｗｉｌｌ ｓｕｐｐｌｉｅｓ ｒｅｌｉａｂｌｅ

ａｎａｌｙｓｉｓ

ｄａｔａ ｆｏｒ ｄｅｓｉｇｎ ｏｆ

ｗａｔｅｒｐｏｗｅｒ ｃｏｎｓｔｒｕｃｔｉｏｎ．

ｋｅｙｗｏｒｄ：ｎｕｍｅｒｉｃａｌ

ｓｉｍｕｌａｔｉｏｎ，酣ｄ，ｔｕｒｂｕｌｅｎｔ ｍｏｄｅｌ，Ｖｏｌｕｍｅ

ｏｆＦｌｕｉｄｓ，

ｓｐｉｌｌｗａｙ，ｆｒｕｓｔａ ｏｆ ｂｒａｋｅ

Ｖ

河海大学博士论文

．１上?

聃

——１一

日

本论文以水利工程中泄水工程的水力设计为背景，研究了自由面湍流场数值模拟
技术，并将之应用于带闸墩溢流坝三维过坝水流场的数值计算。

随着计算机科学和数值计算方法的发展，数值模拟在科学技术研究中具有更大的 优势和发展空间，一方面可以确定物理问题整体舶和局部的细致行为，另一方面在解 决实际工程问题的过程中便于进行方案优选，从而降低费用，缩短周期，此外还可以 替代一些昂贵的、危险的甚至难以实现的实验，因此数值模拟在整个科学技术进步中
起到非常重要的作用。
｝

然而，利用数值模拟技术解决实际工程问题必须要求数值模拟过程中的每一个环 节．如：数学模型、控制条件、数值算法等都是高效的ｊ先进的，并经过充分验证和 确认的，它们决定了数值模拟成果最终的可靠性和实用性。 在前人研究工作的基础上，本论文对数值模拟的关键技术：网格生成技术、数值 求解技术、湍流模型和动边界模拟等进行了较深入的研究，建立了非正交同位网格系 统中自由面湍流场数学模型，并成功地应用于二维溢流坝和带闸墩的三维溢流坝过坝
水流数值模拟。
．

通过本论文的工作，取得了以下研究成果： １）提出了一种简单的、可独立于网格生成方法之外的网格边界正交化方法； ２）利用分区结构网格技术和非正交同位网格系统中ＳＩＭＰＬＥ算法建立了求解复杂流 场的数学模型，提高了计算效率和计算精度； ３）详细建立了一种精度更高的自由面重构方法，能够更加准确地确定运动界面的形
状和位置：

４）利用以上模型首次计算了带闸墩溢流坝三维过坝水流粘流场。针对不同的坝顶水 头、闸墩的型式以及在坝面上的布置完成了系列水力计算，为泄水工程建筑物的水力 设计提供可靠的分析依据。

学位论文独创性声明：

本人所呈交的学位论文是我个人在导师指导下进行的研究工作 及取得的研究成果。尽我所知，除了文中特别加以标注和致谢的地方
外，论文中不包含其他人已经发表或撰写过的研究成果。与我一同工

作的同事对本研究所做的任何贡献均己在论文中做了明确的说明并 表示了谢意。如不实，本人负全部责任。

论文作者（签名）：
（注：手写亲笔签名）

￡色：垒
坠＆：坚

旦垒年

３

月

卫２日

学位论文使用授权说明： 河海大学、中国科学技术信息研究所、国家图书馆、中国学术
期刊（光盘版）电子杂志社有权保留本人所送交学位论文的复印件或

电子文档，可以采用影印、缩印或其他复制手段保存论文。本人电子
文档的内容和纸质论文的内容相一致。除在保密期内的保密论文外， 允许论文被查阅和借阅。论文全部或部分内容的公布（包括刊登）授

权河海大学研究生院办理。

论文作者（签名）：
（注：手写亲笔签名）

互盘垒

竺年

３月２工日

河海大学博士论文

符号表
松弛因子，控制网格边界正交化的程度
Ｐ，Ｑ
ｐ

阻尼参数，控制网格正交性的衰减程度
流体密度 流体运动的速度分量

Ｕ，（ｆ＝１，２，３）

Ｐ，
－“

粘性流体中的法向应力
粘性流体中的切向应力
流体的动力粘度

“

Ｒｅ＝ＵＬ／ｙ Ｕ
上

雷诺数

特征速度
特征长度 流体的运动粘度

＾
膏
占

涡粘系数 湍动能 湍动能耗散率
壁面切应力

Ｆ Ｈｓ

流体体积分数
溢流坝高 溢流坝顶设计水头 溢流坝顶水头
溢流坝流量

Ｈ１

日。
０
，竹

溢流坝流量系数
溢流坝面压力水头
闸墩厚度
闸孔宽度
．

＾
，

ｂ

ＩＸ

河海大学博士论文

第一章绪论
１．１数值模拟技术的工程研究意义
数值模拟、理论分析和模型试验是研究流体力学（水力学）的三种主要手段， ～般来说，这三种研究手段和方法必须相互配合，相互补充，相互促进，才能共同 推进流体力学学科的发展并解决各种工程实际问题。

理论分析方法是在研究流体运动规律的基础上提出各种简化流动模型，建立各
种控制方程，在一定的假使和条件下，经过一系列的解析推导和运算得到问题的解 析解，理论分析的许多方法仍是目前解决实际问题，主要是在初步设计阶段中常常 采用的方法，但理论分析法常常无法用于研究复杂的、以非线性为主的流动现象。

模型试验方法是研究流动机理、分析流动现象、探讨并获得流动新概念，推动
流体力学发展的主要手段，其不足是完成一个完整的试验过程需要解决一系列复杂 的技术问题，所需周期长、费用高。 数值模拟方法的特点是在用比模型试验所需花费少得多的情况下给出流场内部 细节的详细描述。并且只要数值模拟的数学提法（控制方程、边界条件、数值算法 等）是正确的，则可以在较广泛的流动参数（如雷诺数，空化数等）和物理设计参

数范围内较快地给出流场的定量结果，不受试验中固有约束条件的影响。因此随着 计算流体力学和计算机科学的迅速发展，数值模拟在科学技术研究中将具有更大的 优势和发展空间，主要体现在：①数值模拟在某种意义下比理论和实验研究对流体
运动过程的认识更加深刻、细致，不仅可以得到运动的结果，而且可以了解整体的

和局部的细致行为；②在解决实际工程问题的过程中，利用数值模拟可以选择最优 的设计方案或实验方案，从而减少实验次数，周期及费用：③数值模拟可以从理论 上探索流体运动新的现象和规律．而且可以替代一些昂贵的、危险的甚至难以实现
的实验，比如水坝坍塌造成的灾害预报、水下爆炸与核爆炸的过程和效应、大气运 动现象等。因此数值模拟在流体力学研究乃至整个科学技术进步中起到非常重要的
作用。

虽然数值模拟在流体力学研究中占有如此重要的地位，然而要使其一方面能够为
工程设计提供高质量、短周期、可靠的分析设计依据，另一方面为探索流动新概念、 新机理提供新途径，则要求数值模拟过程中的每～个环节，如：数学模型、控制条件、 数值算法等都是高效的、先进的，并经过充分验证和确认的，它们决定了数值模拟成 果最终的可靠性和实用性。

河海大学博士论文

１．２粘流场数值模拟关键技术研究现状 １．２．１网格生成技术
网格生成技术是计算流体力学（ＣＦＤ）发展的一个重要分支，在ＣＦＤ高度发展的

美国，网格生成所需的人力时间占计算任务全部人力时间的６０％，因此将ＣＦＤ技术进
一步推向工程应用，首先必须完善网格生成技术“。。对于复杂程度一般的计算区域边

界，采用结构网格系统完全能够满足工程要求的精度，由于算法相对简单，因而目前 仍被广泛采用。常用的三维结构网格生成方法大致可分为代数生成法“。、椭圆微分方 程生成法”’”和双曲微分方程生成法。”等三类。 代数生成法是利用已知的边界值进行中间插值来生成网格，其核心是确定插值函
数，晟常用的是“超限插值”（Ｔｒａｎｓｆｉｎｉｔｅ ｉｎｔｅｒｐｏｌａｔｉｏｎ），内部网格节点的分布 可以由边界上的伸缩函数来控制，通过一些专门的、不是十分复杂的技术能够实现网 格在边界的正交或近似正交“““，因此，对于几何形状不是十分复杂的流场，代数生

成方法不失为一种简单有效的方法。 椭圆型网格生成法是通过求解椭圆型微分方程来实现网格的生成，该方法的最
早代表为著名的ＴＴＭ方法”’，即通过求解Ｌａｐｌａｃｅ方程来生成区域网格，该方法的 优点是所生成的网格是光滑的，能够处理复杂的几何边界，同时若在方程右端增加

源项交成泊松方程则可以控制网格分布的疏密及倾斜度，其缺点是这样的源项很难 确定，且较难实现内部节点的控制。目前最常用的源项控制方法是根据边界正交性 和网格间距的要求壹接推导出右端项的表达式”’…，椭圆型网格生成方法是最常用
的一种单域贴体网格生成方法。 双曲型网格生成法是通过求解双益型微分方程来实现的，该方法的优点是能够保 证网格在边界处的正交性，容易控制两格的尺寸，且计算时间比椭圆方法少卜２个量

级，缺点是边界上的任何不连续都会传播到内部网格节点，外边界位置不能事先给定。 但随着对一些敏感问题，如边界条件、物面不连续和角点的处理等所开展的讨论和改
进。提高了该方法的适应能力，在近年来得到了比较广泛的应用”’”‘“。。

随着边界复杂程度的提高，生成单域贴体计算网格更加困难，同时网格质量不 能得到保证。影响最终的数值计算结果，为此人ｆ｜、］提出了分区结构网格、非结构网 格和结构／非结构混合网格技术。 分区结构网格系统是目前比较成熟的、并经常采用的构造复杂区域网格的方 法，分区结构网格方法的基本思想是根据边界的特点将整体流场划分为若干个子
区，对每一个子区分别建立网格，并在其中对流动控制方程求解。但在区域划分对

需遵循以下原则：尽量使每个子域的边界简单以方便结构网格生成；各子区的大小
也尽量相近以实现计算负荷的平衡，这一点对并行计算尤为重要。常用的分区网格

２

河海大学博士论文

方法有：对接网格方法和重叠网格方法。
国外分区网格生成技术于８０年代兴起并逐渐形成实用软件，具有代表性的有

ＥＡＧＬＥ…。和ＧＲＩＤＧＥＮ““。前者是以椭圆型网格生成理论为基础，兼有代数生成方法 （超限插值）和表面网格生成技术：后者是以代数网格生成方法为基础，功能不如 前者完善。到了９０年代，二者又补充了许多新的技术，如基于弧长的超限插值函 数，提高正交性的Ｈｅｒｍｉｔ三次型函数，良好的用户界面（６ＵＩ）和加强的与ＣＡＤ接 口能力等，从而发展成为功能很强的网格生成系统。其基本思想是依靠超限插值方 法生成表面和空间网格，然后由椭圆方法光顺，通过调整代数方法中的Ｈｅｒｍｉｔｅ超

限插值或椭圆方法中的源项控制函数来使网格线在交界面处的斜率具有连续性，以
保证交界面两侧网格的正交性。 非结构网格技术于８０年代末９０年代初得到迅速发展“““。，其基本思想是：由 于四面体是三维空间最简单的形状，任何空间区域都可以被四面体单元所填满，即

任何空间区域都可以被以四面体为单元的网格所划分。显然，非结构网格方法舍去
了网格节点的结构性限制，节点和单元的分布是任意的，能够较好地处理复杂边界。 非结构网格方法在其生成过程中都采用一定准则进行优化判定，因而能生成高质量 的网格，且很容易控制网格的大小和节点的疏密度，一旦在边界上指定网格分布， 在两个边界之间即可自动生成网格，因此非结构网格系统在近年来受到了很大的重

视，有了较大的发展，如文献［１６］针对二维粘流场构造出适应于非结构网格的有限 体积时间推进格式；文献［１７］利用点自动生成的Ｄｅｌａｕｎａｙ三角化方法构造非结构
网格，同时处理成自适应网格形式，能够计算带移动边界的瞬时流动问题等。非结 构网格的主要缺点是所需计算内存较大，计算效率较低。 结构／非结构混合网格系统则兼顾了二者的优点，对于几何形状复杂的三维问 题可以由近边界处的非结构网格和远离边界的结构网格构成，具有灵活模拟复杂边
界和计算效率高的特点。

目前分区结构网格的生成方法和并行计算方法、结构／非结构混合网格系统是计 算具有复杂几何形状流场的主要方法，其进一步的发展方向是减少工作量、实现自动

生成和自适应加密，具有良好的人机界面及可视化，强调更有效的数据结构等，最终
的目的是能够快速生成适用于粘流场计算的高质量计算网格。

今后，网格生成技术的研究重点应当是网格的专项技术研究，如：边界正交化技
术、分区结构网格及并行计算技术、网格的自适应加密技术以及动态网格技术等。

河海大学博士论文

１．２．２数值求解技术 数值求解技术实际上是一种离散近似的计算方法，各种各样的流体力学问题必须
从给定的微分方程或基本定律出发，建立物理上合理、数学上适定、能够在计算机上 进行计算的离散化数学模型。大多数数值求解技术的基本思想可归纳为：把原来在时 间和空间上连续的物理量的场（如速度场、温度场等），用有限个离散点上的值来代 替，然后按一定的方式建立关于这些值的代数方程组并求解，从而获得物理量场的近 似解。目前流场计算的主要数值求解方法有：有限差分法（ＦＤＭ）、有限元法（ＦＥＭ）、

有限分析法（删肘）、边界元法（ＢＥＭ）、有限体积法（ＦＶＭ）等，其中前四种数
值计算方法的特点及适应范围在众多的文献中均可查阅［１８，１９，２０，２ｌ，２２，２３］，下
面重点介绍有限体积法。

有限体积法由ＭｃＤｏｎａｌｄ在１９７１年首次用于求解二维欧拉方程，１９７２年被 Ｐａｔａｎｋｅｒ等人用于ＳＩＭＰＬＥ算法计算恒定不可压流。有限体积法的基本思想是：将计 算区域划分成一系列规则或不规则的单元或控制体积，并使每一个控制体积包含一个 网格节点，将待解的微分方程对每一个控制体积分得到一组离散方程，其中的未知数 是网格节点上因变量的值，首先计算出通过每个控制体边界沿法向输入（输出）流量 的动量通量，然后对每一个控制体分别进行流量及动量平衡计算，便可得到计算时刻

各个控制体内流动参数的平均值（平均流速、平均压强等）。因此，有限体积法实际 上是对推导原始微分方程所用控制体方法的回归，与ＦＤＭ和ＦＥＭ的数值逼近相比 其物理意义更加直接明确，因而在流场计算中获得了广泛的应用“。。’“。“。。
若边界通量的计算只使用前一个时刻的数值，为显式ＦＶＭ，反之，当涉及到下 一个时刻的数值时，则为隐式ＦＶＭ。此外还有半隐式ＦＶＭ“““。。

１．２．２．１数值方法对离散格式的要求
数值求解首先要对微分方程进行离散，然而不同的离散格式将带来完全不同的 结果，因为同一微分方程不同的离散格式并不等价，也不一定能保持原来微分方程 的性质，对渐变光滑的情况其结果可能相差不大，对迅变以致间断的情况其解则可

能相差很大。因此利用ＦＶＭ在建立和使用离散格式时需要满足以下条件才能保证 数值求解的精度和可靠性。①相容性：即当时空步长趋于零时截断误差（离散方程 与原始方程之差）为零；②守恒性：是指在每一个空间有限的控制体积内保持质量 和动量平衡，具有守恒性的数值格式在计算过程中不会产生守恒误差，并且计算精
度和稳定性较好，同时只有守恒格式才能处理计算空间中可能产生的间断面，即流

场中可能产生的激波：③逆风性：即空间导数离散时根据流速方向选择差分方向，
从而使粘性的耗散是合理的：④无虚假振荡：是指在计算间断（或解梯度大）时要

河海大学博士论文

求格式具有足够耗散，不出现虚假振荡，当然，无假振要求与计算精度的要求是相

互矛盾的，目前大多采用册格式计算有大梯度和间断的解；⑤高分辨率地捕捉间
断：是指要求人工粘性和格式粘性的强度与解的光滑性相适应，在光滑区粘性自动

减少以提高精度，在间断点附近或梯度大的区域又自动加强，以防假振并达高分辨 率，目前一些优良格式的间断分辩率可达二到三个网格：⑥计算稳定性：是指当数
值解在受到微小的扰动后不发生突变现象，目前很难通过数学推导给出实用的稳定 性条件，应用中几乎总是通过数值试验来确定时间步长，对时间步长起决定作用的

因素有时间离散格式、网格的规则性及边界条件的处理等；⑦解的收敛性：数值格
式的准确解与原微分方程的真解之差成为离散化误差，当时空步长趋于零时，离散 误差亦趋于零，则称格式具有收敛性，显然相容性只保证数值格式收敛于原微分问

题，更重要的是数值解收敛于真解；⑨精度：数值解的精度取决于算法、网格、边
界处理等诸多因素，除网格形状要求规则外，还要求网格光滑，步长渐变，否则要

按照步长比进行校正：⑨健全性及通用性：健全性是指能保证数值解不园所使用的 网格而异，例如可选择一个简单的明渠均匀流问题，所用格式在不同网格上的计算
结果必须一致，此外格式中不含专门为计算而引入的、需要人工加以调整的参数（如 人工粘性系数），目的是保证结果的客观性，有文献预言，随着格式其它性能的完 善，未来十年内对格式的主要要求为健全性；通用性则要求格式适用于各种流态（恒 定流／非恒定流，低速流／高速流，连续／间断）、任意几何形状，多种网格布置、多 种边界条件等尽可能多的情况。 能够满足上述大多数要求的格式可以称作为高性能计算格式”…。 为了达到上述要求，首先必须了解各种离散格式的优缺点及其改进方法。下面 对微分方程目前主要的空间和时间离散方法进行综述。

（１）空间离散化方法： ①中心格式：中心格式的原理是泰勒级数展开式，其最大优点是计算简单、精度
高，有算例表明，在加粘后其计算量仍只是二阶ＴＦＤ格式的］／５，因而成为计算流 体动力学的代表方法之一，常用来计算光滑流动；缺点是计算解在计算空间中会产

生假性振荡，同时稳定性较差，目前实际计算中最有代表性的中心型数值方法是
Ｊａｍｅｓｏｎ““等首先提出的。 ②逆风格式：逆风格式是根据流速方向来选择差分方向，优点是稳定性较好，如 二阶单侧逆风格式允许的ｆ为同样点数的对称格式的两倍，对急流和出流边界必须 或经常使用逆风格式，入流边界使用单侧逆风格式是不可取的，但是，一阶逆风的 耗散往往太强，精度低，二阶线性逆风格式在计算间断解时常产生和中心格式相似

的振荡，这就要求在经典逆风格式的基础上研究非线性的单调性保持或ＴＦＤ格式。

河海大学博士论文

文献［３２，３３］以中心格式和逆风格式为基础建立了通量修正格式，该格式既防止假
性振荡又可达到高阶精度，且格式较为健全。

③ＴＶＤ格式
ＴＶＤ（Ｔｏｔａｌ Ｖａｒｉａｔｉｏｎ

Ｄｉｍｉｎｉｓｈｉｎｇ）的含义是总变差减小，对于微分方程

的数值解“（一）来说，其总变差对应网格点的变差为：ＴＶ（ｕ（?））＝∑ｌｕ（ｘ，）－－Ｕ（Ｘ。）Ｉ
ｉ＝ｌ

如果对于不同时间步长的近似解有：ＴＶ（ｕ”、）≤ＴＶ（ｕ“） 则称总变差是减小的，即具有ＴＶＤ性质。

ＴＶＤ格式的建立是以方程罢＋罢：０的柯西问题为基本对象，在间断形成前， 拼盘
解的全变差Ｚｙ保持常数，以后ｒ矿下降，同时在解的递推过程中不产生新的极值，
原有极值也不增加，从而初值的单调性得以保持，据此可建立严格意义下的ＴＶＤ格
式。

建立二阶ＴＶＤ格式的途径有：修改通量（在粘性过大的～阶逆风格式中加入

受限制的反扩散通量）：通量限制（通过引入限制因子对反扩散项加以限制以提高
精度的阶数）：变量插值（在变量插值时对坡度进行限制）［１８，３０，３４Ｊ。

④ＮＮＤ格式 ＮＮＤ格式即为无波动、无自由参数的高分辨率格式，由张涵信所构造““，其 基本出发点是在Ｅｕｌｅｒ方程和Ｎ．ｓ方程的右端附加一个三阶导数项，并使该项的系
数在激波上游为正，在激波下游为负，以此来抑制差分解在激波上下游的非物理波 动，格式中不包含任何需要人工调整的参数，并具有二阶精度的ＴＩｅＤ性，对应不同 的时间推进形式，可以得到不同的显式全离散格式ＮＮＤ—ｌ到ＮＮＤ一５，并可构造

出隐式时间推进无条件稳定的ＮＮＤ格式，其中ＮＮＤ一２到ＮＮＤ一５格式对时间均 为二阶精度。因此ＮＮＤ格式以其准确可靠，程序简单，计算量小等特点在科学研 究及工程计算领域得到了广泛的应用并取得了良好的效果。 需要说明的是中心格式和逆风格式可以相互转化，中心格式添加人工粘性后可 解释为逆风格式，中心格式也可改写为逆风格式加上反扩散项，对后者加以限制即 构成ＴＶＤ格式。反之，逆风格式也写作中心格式加耗散项，对后者加以限制也构 成ＴＶＤ格式。因此两类格式的主要区别在于如何处理耗散，中心格式的人工粘性
和逆风或ＴＶＤ格式的粘性是相通的，故可以相互借鉴。

（１）时间离散化格式
时间坐标轴与空间坐标轴不同，它是单向的，只能是过去影响将来，时阊积分

河海大学博士论文

必须是采用相对ｒ。的单侧格式，且往往与空间离散化分别对待。时间积分主要分为 显式和隐式两类。显式是指可以由已知的ｆ。时的初值直接计算出ｔ。时刻的未知解，

不必求解方程组。其优点是每个时间步的计算工作量和存储量均较小，且程序简单， 缺点是时间步长受到稳定性要求的限制，当局部网格很密时，时间步长必须非常小
而使整个计算次数很大。隐式是指可同时由，。和，。时的解逼近空间导数和源项，

故通常需要求解非线性微分方程组。其优点是在线性分析时往往是无条件稳定的，
计算时间步长可以取较大的值从而提高总体效率。缺点是需要对每一个时间步长求 解方程组，因而计算工作量和内存占用量均很大，程序设计也比较复杂，主要用于 恒定流或解的时间变率小、流动时间历程长、高雷诺数流等情况。

应当注意，时间离散格式的选择对址的确定和计算量至关重要，址应在保证 稳定性的前提下由精度、单调性等综合考虑确定，通常需要通过数值试验来最终确
定。

１．２．２．２线性方程组的解法 线性方程组的解法有两类：直接法和迭代法。直接法可通过有限步得到准确解， 求解精度高，但占有内存大，仅适用于小型的方程组，而迭代法是从给定初值出发，
通过某种算法反复加以校正，直至收敛。迭代法主要用于求解大型多维问题的不对

称稀疏线性方程组，缺点是不一定能保证迭代收敛，且终止迭代的准则不明确。
常用的迭代法是Ｊａｃｏｂｉ法和Ｇａｕｓｓ．Ｓｅｉｄｅｌ（ＧＳ）法，在ＧＳ方法中为加速收敛一般

均引入松弛因子田，称为ＳＯＲ法（逐次超松弛法），若ｍ选取适当，可大大加速收 敛，∞的取值范围为［０，２］，若∞＞ｌ，则称超松弛法，ＳＯＲ法通常性能更好，因而
在计算流体力学中应用广泛。ＳＯＲ法的一种变型为ＳＳＯＲ（对称逐次超松弛法），

其收敛速度常远远超过ＳＯＲ，ＳＳＯＲ的每步迭代由两小步组成，一是向前ＳＯＲ，一 是向后ＳＯＲ，实际计算中需通过数值试验来选取甜。 １．２１３湍流模型技术
工程实际中存在的流动大部分都是湍流流动，湍流问题一直是计算流体力学中

的一个重要难题，湍流的特点是具有很强的随机性、非恒定性，并且无一例外地都 是三维问题。由于湍流现象的复杂性，湍流运动以及与之相联系的热和物质的输运

现象都很难描述，从而也极难进行理论预测。如果要通过计算来解决这一问题，就
不可避免地对湍流输运过程提出各种假设，采用一些经验性的结果和假定，把湍流 输运过程中的各种物理量与时均流场联系起来，这就是湍流模型研究的内容。所谓

湍流模型就是以雷诺时均运动方程与脉动运动方程为基础。根据现有的理论和经验 引进一系列假设而建立起来的一组描述湍流平均量的封闭方程组。湍流模型的选取

河海大学博士论文

不仅要考虑其模拟复杂流动的可靠度和精确度，而且还需要确定计算所需的费用和 处理问题所需的时间，后者对于工程应用尤为重要。通常评价湍流模型优劣的标准 是：能够适用于较多种类型的水流现象；具有足够的模拟精度；人力和计算机费用 适度；复杂程度适当。显然，那些适用范围大，精度高的湍流模型必然复杂，消耗

费用较大，而好的湍流模型正是在这些相互矛盾的要求之间达到某种理想的平衡。 目前的湍流模型主要有雷诺时均模型、大涡模拟模型（Ｌａｒｇｅ Ｓｉｍｕｌａｔｉｏｎ－ＬＥＳｌ和直接数值模拟湍流模型（Ｄｉｒｅｃｔ
Ｔｕｒｂｕｌｅｎｃｅ—ＤＮＳｎ。
Ｎｕｍｅｒｉｃａｌ
Ｅｄｄｙ Ｓｉｍｕｌａｔｉｏｎ

雷诺时均模型：该模型在总体上可划分为基于Ｂｏｕｓｓｉｎｅｓｑ假定的紊动粘性模型

和基于雷诺应力输运方程的二阶矩封闭模型，前者又可分为混合长度模型、一方程
模型和二方程（如ｋ—ｅ，ｋ—ｕ，ｋ－ｒ等）模型，后者则有代数应力模型（ＡＳＭ）和输
运方程模型（ＤＳＭ）两大类。

Ｂｏｕｓｓｉｎｅｓｑ假定的基本思想是：认为雷诺应力各向异性张量与平均场的应变率

张量成正比，即湍流应力的作用类似于层流应力，但其比例系数为一流动变量一涡
粘性系数，正是对涡粘性系数不同的计算方法导出了零方程、一方程、二方程等湍 流模型。不难想象，在涡粘性模型中，雷诺剪切应力的方向与平均场的应变率相同 这一特点可以说是涡粘性模型的主要缺陷之一，即雷诺应力与平均场应变率之间的

关系为线性关系，使得ｋ－ｅ模型无法再现流动中（如在方管中）出现的二次流“…， 这是由于复杂流动中，湍流场的变化总是滞后于平均场的变化。 在二阶矩封闭模型中，雷诺应力与平均速度之间的关系隐含在一组联立的偏微
分方程中，它们之间的相互关系包括一些复杂的物理机制，同时在雷诺应力方程中

的压力一应变率相关项为雷诺应力输运过程中提供了快慢两种不同的变化机理，在 一定程度上反映了湍流的非线性作用，然而通常形式的二阶矩封闭模型对雷诺应力 来说基本上是线性的，所包含的非线性作用非常小，还不能令人满意地用来解决复
杂的工程问题。

上述各类雷诺时均湍流模型的特点及应用可参见文献［１８，３３，３７，３８，３９，４０］，这 里就不再详细展开。

大涡模拟模型（ＬＥＳ）：通过不断研究人们认识到，在湍流流动中除存在许多随机
性很强的小尺度涡运动外，还有一些组织很好的大尺度涡结构，小尺度涡运动受流
动边界条件的影响很小，且近似是各相同性的，流动中大部分质量、动量或能量的

输运主要来自于大涡运动。因此，如果将比网格尺度大的大涡运动通过数值求解 Ｎ．Ｓ方程来计算，而将比网格尺度小的小涡运动利用湍流模型来模拟，使得通过模 型提供的计算在整个求解中只占很小的比重，降低最终结果对模型的依赖性，从而 达到提高计算精度的目的，这就是大涡模拟的基本思想。目前大涡模拟方面的工作

河海大学博士论文

开展不多，主要集中在气象学的流动问题或基础性研究及简单的流动问题，如文献 ［４ｌ，４２】利用大涡模型分别对方柱绕流和明渠流动进行了数值模拟。关于ＬＥＳ更详细
的介绍可参见文献［３８］

直接数值模拟模型（ＤＮｓＴ）：模型化总是存在缺陷驱使人们寻求更好的解决
端流问题的途径，由于Ｎ．Ｓ方程本身就是封闭的，不需要建立模型，由此提出了不 引入任何湍流模型，直接求解完整的三维非定常Ｎ．Ｓ方程，对湍流的瞬时运动进行 直接数值模拟，各种感兴趣的物理量都可通过对瞬时量的平均运算来取得，这就是 直接数值模拟湍流。该方法显然有许多优点，如精确、可提供流场的全部信息，能

实现实验室中很难做到的流动条件控制且比实验更为经济，但是直接数值模拟的计 算工作量却非常巨大，由于湍流运动所包含的单元（即小涡）比运动区域的尺度要
小得多，典型的数量级是流动区域尺度的ｌＯ。倍，而数值计算所需要的网格必须比

湍动单元的尺度更小，在三维情况下，为了覆盖流动区域至少需要ｌＯ。个网格点，

如此多的网格点和相应的计算次数是现代计算机的储存能力和运算速度所不能达 到的。目前国际上进行的湍流直接数值模拟还仅限于较低雷诺数（Ｒｅ＝２００）和简单
几何边界条件的问题，主要用于湍流的基础研究．如发现新结构，揭示新机理，提 供新概念，检验和改进湍流模型等。 显然，无论是大涡模拟还是直接数值模拟，主要的困难不仅在于计算枧的限制，

还有方法本身的问题，限于目前计算机的发展水平，短期内还不能够利用上述两种方 法对工程问题中的复杂流动进行数值模拟，因此对现有的湍流模型加以改进，提出更
加合理的新的湍流模型依然是解决工程实践中湍流问题的主要途径。 重整化群（Ｒｅｎｏｒｍａｌｉｚａｔｉｏｎ Ｇｒｏｕｐ。ＲＮＧ）是一种用于构筑许多物理现象模型的 通用方法，该方法起源于量子力学和高能物理中对基本粒子场的研究，自７０年代后

期开始把ＲＮＧ方法引入到湍流研究领域，Ｙａｋｈｏｔ和０ｒｓｚａｇ…。于１９８６年应用ＲＮＧ 方法建立了第一个湍流模型，简称为ＲＮＧｋ－ｅ。通过数值模拟表明它较之传统湍流 建模方法具有显著的优越性和发展潜力。自从基于ＲＮＧ方法的湍流模型提出以后，
许多学者应用该模型进行了数值模拟，Ｓｐｅｚｉａｌｅ Ｇ．Ｇ．”“和Ｙａｋｈｏｔ““利用该模型求

解了后台阶绕流中的湍流分离现象，蒋莉““将ＲＮＧ ｋ—ｅ湍流模型与壁面率相结合， 应用于９０。弯曲槽道湍流流动的数值模拟，结果表明，该湍流模式可以有效地模拟有 曲率影响的湍流流动，李玲”“利用该模型对钝体绕流的尾流场进行了数值模拟，结 果优于标准的ｋ—ｅ湍流模型。由于利用ＲＮＧ方法建立的湍流模型中不包含任何经验 常数和可调节的参数，其模型参数是利用ＲＮＧ理论精确推导出来的，因而是通用的：
该模型适用于各种雷诺数范围，包括层流、过渡流以及充分发展的湍流。此外，ＲＮＧ

模型能够较好地反映流动的各向异性，对于与时间相关的大尺度运动也能给出真实的
模拟。因此ＲＮＧ湍流模型在湍流计算中具有进一步的推广和应用的价值。

河海大学博士论文

１．２．４自由面模拟技术 自由液面问题是计算流体力学（ＣＦＤ）中非常重要的研究内容，自由液面模拟的
合理性直接关系到含自由液面问题数值模拟成果的可靠性及精度。由于自由液面形状 的复杂性和整体位置的不断变化，特别是当自由液面发生破碎丽使得该闯题成为强非 线性问题，所以精确计算和模拟自由液面的难度非常大，可以说，自由液面问题已经

成为计算流体力学领域的重要分支和前沿课题。目前在该领域的研究处于领先地位的 是日本的东京大学和日本船研所（ＳＲＩ），此外还有美国的ＬＡＷＯ大学和韩国的ＳＲＩ，，
我国的７０２研究所、上海交通大学、大连理工大学等单位对自由面模拟技术开展了一 些卓有成效的研究工作，但由于问题的复杂性，上述结果还远远没有达到令人满意的 程度ｃ４８－４９’５引。

下面对处理自由液面的一般方法及其应用作简单描述，重点介绍ＶＯＦ方法。 刚盖假定方法是将自由面看作是一个可移动的固体壁面，直接采用固体壁面的 “不可穿入”条件，仅适用大水体的宏观运动； 静压假定法”“则认为流体中的压强沿水深的分布服从静压规律，目前几乎所有
的大体积水流数值模拟都以此假定为基础；

高度函数法”“是将自由面处的水深ｈ定义为坐标（ｘ，ｙ，ｚ）和时间ｆ的函数，通过求 解高度函数ｈ满足的微分方程来确定自由液面的位置，该方法适用于非恒定自由表面
问题．但无法处理标高函数为多值函数的情况，如射流、波浪破碎等问题。

流线迭代法”“是直接从Ｎａｖｉｅｒ－Ｓｔｏｋｅｓ方程出发，迭代求解压力分布及流线位置。 文献［２３１采用此方法计算了溢流坝反弧曲线段的自由水面。 除上述几种方法外，目前处理自由液面问题的典型方法是欧拉模型中的ＭＡＣ法
和ＶＯＦ法。ＭＡＣ（Ｔｈｅ Ｍａｒｋｅｒ ａｎｄ Ｃｅｌｌ Ｍｅｔｈｏｄ）是设想在整个流场中均匀分布没有

体积和质量的微小颗粒一标记点，这些标记点以当地流体质点的速度随流体一起运 动，显然只要计算出标记点在空间的位置即可确定自由表面的形状和位置”‘…。该方 法的主要缺点是需要很大的计算机内存和计算工作量，并且对于非均匀流场在网格中 会出现虚假的密度很高或很低的标记点，造成自由液面形状的失真，而ＶＯＦ方法则 克服了上述不足，是目前研究自由液面问题的理想方法。
ＶＯＦ（Ｖｏｌｕｍｅ

ｏｆＦｌｕｉｄ）”“方法是定义这样一个体积分数Ｆ，在有流体的点取值

为ｌ，在无流体的点取值为０，这样，在一个网格单元内，Ｆ的平均值则代表单元内 流体所占的份额，若Ｆ的均值为１．则网格充满流体，若Ｆ的均值为０，则为空网格，

若Ｆ的均值在Ｏ和ｌ之阎，则网格包含自由液面，为自由面网格。由于ＶＯＦ法追踪
的是网格中的流体体积，而不是追踪流体质点的运动，因而具有容易实现、计算量小

和精度高等优点，并且可以处理自由面折叠、自由面入水等强非线性闷题，成为目前
Ｏ

河海大学博士论文

自由面运动模拟的主流方法。

在ＶＯＦ方法中，关键技术是建立运动界面的输运与重构模型，如文献［５７１研究了 运动界面追踪技术中的界面重构算法，并讨论了由此引起的数值误差问题，文献［５８１ 以ＶＯＦ理论为基础建立了三维空间自由面输运及重构的一阶模型，文献［５９１将有限体 积法和动网格技术应用到运动界面追踪，自由面的位置和形状利用运动边界条件和空
间守恒率确定等，总之７０Ｆ方法中有许多问题有待进一步深入研究。

Ｉ＇３本论文的工作意义 本论文的主要工作是研究流场数值模拟中的关键技术，并将之应用于水利工程中 带闸墩溢流坝过坝水流三维流场的水力计算。 ２１世纪人类面临的主要问题是生态保护型可再生能源的安全供电、用于抵抗饥 荒、贫穷和疾病的优质而充足的供水以及可满足人类基本需求的能源和食物，这些问 题的解决在很大程度上依靠设计开发新的水电工程来实现，水电工程对世界经济繁荣 和可持续发展起着至关重要的作用。在我国随着水电技术的迅速发展，高坝的建设日 益增多，其特点是水流落差大，流量大，功率大，多位于河谷狭窄地带，带来的一个 突出问题是泄洪消能、防冲及安全．可以说我国在该领域遇到问题的技术难度居世界
第一，因此需要更加重视大坝安全问题。

在泄水建筑物一溢流坝的水力设计中需要解决的问题不仅包括坝面下游反弧段 上的掺气、振动、空蚀防蚀及雾化等问题，同时需要解决坝面上游压力分布、水面线 的位置及形状、泄洪流量、闸墩对流量及水面线形状的影响等问题。目前主要的研究 手段是物理模拟（水工模型试验）和数值模拟．物理模拟有着公认的实际意义和科学 价值．它可以预演工程中非常复杂的三维水流现象，检验设计方案的可靠性和准确性， 但是它具有费用高、周期长、无法准确测取工程的某些湍流特征量和不便于进行多方 案优选分析等弱点，因而现在物理模拟中已有相当部分的工作为数值模拟所替代。数 值模拟便于进行方案比较和优选，同时研究高速水流、空化和空蚀现象已在相当程度 上可借助于湍流的数值模拟来实现。 在水利工程中，基于水深方向的二维模型在一定程度上能反映出坝面的基本流动 规律并解决了水利工程中的一些实际问题，如溢流坝反弧曲线的设计，因此在８０—９０ 年代得到了广泛的研究和应用。许协庆利用有限元方法计算了泄洪洞进口水流和溢流 坝反弧水流…。；顾兆勋等利用Ｎ—ｓ方程及相应的湍流模型模拟了二维非定常带自由 水面的湍流流动。“；王奇峰采用极坐标系下的水流控制方程和ｋ－ｉ湍流模型模拟了 反弧水流的水力特性”１；魏文礼利用正交曲线坐标系和曲率修正的ｋ—ｅ模型模拟了 二维溢流坝流动，并提出了一种抗蚀性能较优的反弧曲线形式“…。但是这些成果还 不能满足水利工程实践的需要，因为泄水工程中的水流具有明显的空间三维特性，二

河海大学博士论文

维模型无法揭示复杂流场内部的空间特性，必须借助于三维数学模型。目前国外在水
利工程中河口航道的三维水流数值模拟方面开展了较多的研究工作，提出了许多新模

型和新算法．国内学者对河口三维水动力学的研究也给予了很大的关注……“，开发 出一些有特色的三维潮流场计算软件。与此相比，在泄水工程的水动力学研究方面， 利用三维湍流模型全面揭示整个流场水力特性的研究工作开展不多，其中马福喜利用 修正的ｋ－ｅ湍流模型和流体体积分数法模拟自由表面计算了重力坝及其下游消力池 中的水流流场和一个三维水跃流场…。，金忠青利用ｋ—ｅ模型数值模拟了带有自由水 面的三维强紊动流场“…。对溢流坝坝型和平面布置（如：顶部溢流曲线、闸墩的形
状、位置及尺寸等）的综合三维数值模拟研究工作开展不足，而这些因素将直接关系

到工程的泄流能力、抗蚀体型的优化乃至整个大坝的安全。 正是由于泄水工程的水力设计对水利大坝工程的泄洪效率及安全营运起着至关
重要的作用，因此必须进一步加强对泄水工程复杂三维流场特性的研究，其中包括建 立具有足够精度的水动力学数学模型和提供高效、可靠的数值计算方法，在此基础上 则可以在较多的设计方案和大坝体型范围内快速而准确地给出整个流场的定量结果，

从而把大型水利工程的建设建立在更扎实的科学研究的基础之上。 １．４本论文的研究内容 １）对网格生成技术进行深入研究，详细给出了网格的代数生成方法：研究了边界正 交化技术，给出了一种简单的网格边界正交化方法；建立了非对接分区结构网格生
成方法，包括分区交界面上的数据结构形式和数值计算方法：

２）利用有限体积方法建立了非正交同位网格系统中的ＳＩＭＰＬＥ算法，详细给出了运 动方程的离散形式，数值计算采用二阶精度的通量修正格式，并完成了倾斜方腔驱 动流、后台阶绕流、圆柱绕流等流场的数值模拟： ３）研究比较了各种湍流模型的特点及应用范围，建立了ＲＮＧｋ一占湍流模型，并以水 翼粘流场为例完成了与标准ｋ一占模型的计算比较：

４）对运动界面模拟技术进行了比较深入的研究，针对ＶＯＦ方法中关键技术一自由面
重构方法详细给出了利用倾斜直线段模拟运动界面的方法，并对旋转流场中的马蹄
铁形界面运动、平移流场中圆形界面的运动、剪切流场中圆形界面的运动进行模拟：

５）利用建立的数学模型对ＷＥＳ型溢流坝过坝水流二维流场进行数值计算，给出了三 种不同坝顶永头下的自由水面线位置以及五种坝顶水头下的坝面压力分布曲线，并 与实验结果进行了比较。同时计算了流量系数及坝面流速分布。
６）研究了带闸墩溢流坝三维过坝水流粘流场，针对闸墩的型式及坝面的布置计算了

闸墩对水面线形状、坝面压力、流量系数等设计参数的影响，并将不同墩型与布置
形式的计算结果进行比较。

２

河海大学博士论文

第二章网格生成技术研究
在对流动问题进行数值计算之前，首先要解决如何将流动区域离散化的问题，目 前已经发展出许多种将物理区域进行离散以生成计算网格的方法，这些方法统称为网
格生成技术（ｇｒｉｄ
ｇｅｎｅｒａｔｉｏｎ

ｔｅｃｈｎｉｑｕｅｓ）。必须注意，网格生成技术在流场数值模拟过

程中占有非常重要的地位，网格生成质量的好坏直接关系到计算结果的精度，因此需
要引起足够的重视。

由于众多的网格生成方法各有特点，在第一章中己经介绍，关于一般的网格生成
方法可参阅文献［６７，６８，６９，７０，７１，７２１，因此本章主要针对网格生成方法中的专项技术进

行讨论，这些讨论有利于在数值计算之前帮助建立合理的计算网格系统。 ２．１有限体积方法中处理不规则区域的网格生成方法 工程中大多数的流动问题都具有一个不规则的复杂几何边界，如水利工程中的溢
流坝、泄洪道、水轮机蜗壳等过水建筑物，船舶工程中的船体及附体，天然河道中的 水流等问题都包含了非常复杂的几何边界，为此需要研究针对不规则区域的网格生成
方法，常用的方法有：

２．１．１采用阶梯形边界逼近真实边界 当流场边界为曲线或斜直线时，一种最简单的边界逼 近方法是采用阶梯形直线段逼近，如图２．１，这种方法的
特点是同一网格线上的网格点（或控制体）的个数是不等
＼
，＿

真实边界

ｌ

的，同时计算边界为相互垂直的锯齿状粗糙表面，当网格
较稀疏时会带来较大的计算误差，随着网格在曲线边界处
的不断细化，这种误差会逐渐减小。
／ ／

≯
，

由于这种网格构造简单，可以适应任何形状的边界， 若同时配合网格的局部加密技术，该方法仍然不失为一种
处理不规则边界的方法。

，

ｔ ｔ

计算边界
图２．１阶梯形边界逼近曲线边界

需要注意的是，当采用这种方法逼近真实边界时，如果把求解区域限制在计算边
界之内，则对每一个具体问题都需要单独开发计算程序，如果把计算区域扩充为一个

规则区域，则通过特殊的处理可以针对不同的问题使用同一个规则区域内的计算程 序，这是处理形状不是特别复杂的计算区域的有效方法，称之为区域扩充法（Ｄｏｍａｉｎ
ｅｘｔｅｎｓｉｏｎ

ｍｅｔｈｏｄ），具体细节可参见文献［７３］。

河海大学博士论文

２．１．２采用曲线坐标变换 在水利工程中经常采用的一种处理不规则边界的方法是引入广义曲线坐标进行 曲线坐标变换‘６８ｔ２３，７４，＇７５，７６】，即将原来的物理空间（ｘ，Ｙ，ｚ）变换到广义曲线坐标空间 （孝，刁，ｆ），同时保证物理空间的几何边界与曲线坐标系中的坐标轴重合或平行，从而 将（ｘ，Ｙ，ｚ）空间中不规则的物理区域变换到（｛，７７，ｆ）空间中形状规则的计算区域，如图 ２．２所示，坐标变换公式为
孝＝毒（ｘ，Ｙ，ｚ）’
ｒ／＝ｒ／（ｘ，Ｙ，：），

ｆ＝ｆ（ｘ，Ｙ，２）

因此如果在变换后的计算平面上进行数值计算，则网格的生成就十分简便。

物理空间

计算空间

图２．２曲线坐标变换

但是，由于数值计算是在变换后的广义曲线坐标系中进行的，因此描述流体运动 的控制方程也需要进行相应的变换，从而在方程中会出现一些附加项，使原本十分复
杂的控制方程变得更加复杂，同时这些附加项的物理意义不明确，而且在对附加项的 导数进行离散时增加了离散误差，这是曲线坐标变换方法的主要不足之处。

２．１．３边界贴体非正交网格 在计算复杂几何区域的流体运动时，经常使用边界贴体非正交网格系统【３”，其中
网格形式可以是结构网格、分区结构网格，也可以是非结构网格，这类网格的优点是

可以适应任意的几何形状，同时比正交曲线坐标变换更容易实现。由于非正交网格的 网格线是沿边界布置的，因而相对于阶梯形逼近曲线边界的情形更容易应用边界条 件，特别是当使用分区结构网格或非结构网格时，还可以沿流动方向布置网格以及在 物理量变化剧烈的区域进行网格加密，因此目前在众多的商用软件中均采用边界贴体
非正交网格系统。

当然，非正交网格也存在一些不足，如方程中由于网格非正交而产生的修正项增 加了方程求解的难度；非正交网格可能导致出现非物理解；此外变量在非正交网格系 统中的布置方式对求解精度和算法的有效性都将产生影响。

１４

河海大学博士论文

２．１．４分区结构网格
结构化网格是指网格系统中节点排列有序，相邻节点的关系明确；分区结构网格 是指把一个复杂的计算区域划分成若干个块，每～块内均采用结构化网格，但不同块

内的网格系统可以不同。分区结构网格的优点是能够用于几何形状复杂的流场计算， 同时可以使用结构网格各个系统中成熟的求解技术，适合于并行计算。便于提高计算
效率［７７，７８ ７９１。

２．１．５非结构两格 与结构化网格不同，在非结构化网格中节点的位置无法用一个固定的法则进行有 序地定义，但是它可以适应任何形状的几何边界，若在非结构网格系统中应用有限差 分方法和有限体积方法，则可以使这两种数值计算方法对不规则边界的适应性增强到 与有限元方法等同的程度，这是非结构网格的主要优点１１７１。非结构网格的缺点是网格
生成的工作量大，离散方程的求解速度较慢。
２．２

网格生成的代数方法一双边界方法 在众多的网格生成方法中，代数方法以其简单、通用性强而被广泛使用，若同时

配合网格的区域加密技术、边界正交化技术和分区网格技术，其适用的范围和精度将

被大大提高，因此首先介绍一种具有通用意义的代数网格生成方法一双边界方法Ｉｔ，９，舯１。
设在物理平面上有一不规则区域ＡＢＣＢ，其中ＡＢ，ＣＤ为两条不直接相联的边界， 首先选定这两条边界上的刁值，分别为ｒｉｂ，矾，于是这两条边界上的ｘ，Ｙ仅随ｆ变化，
即

ｋ＝ｘ＾（善）．ｙ。＝Ｙ．（舌）
Ｘ，＝ｚ，（考）， Ｙ。＝Ｙ，（孝）

其中：６，ｔ分别表示底边ＡＢ和顶边ＣＤ。 双边界方法提供了～种在两个边界之问内插的方法，从而确定内部节点的位置， 具体公式如下：

ｘ（善，Ｊ７）＝工ｒ（ｆ）＾（７７）＋ｘｚ（亭）也（Ｊ７）＋！ｉ！堕ｊ；；三塑２ｈ，（ｒｉ）＋‘！！丛ｊ尹＾４（叩）
彤，野）＝Ｙｌ（嬲（椰＋Ｙｚｇ）＾：（７）＋鱼生ｉ导盟缟（叩）＋鱼尘粤业‰（彳）
ｄｒ／

ｄ１７

其中：

河海大学博士论文

（１）ＸＩ

Ｙ－，ｘ２，Ｙ２为边界１和边界２上的离散点坐标：

（２）ｈｉ＝２矿一３叩２＋１，ｈ２＝一２ｒ／３＋３１／２，也＝ｒ／３—２ｒ／２＋ｒ／，ｈ４＝矿一矿

若将＾（ｆ＝１，２，３，４）表达式中的叩（２了乞击）用一维伸缩函数

．

刚叩＋（Ｉ－Ｐ丁１一、ｔｈ（Ｑ（１ｒ－ｒ／））］

（３）掣－－－Ｋ１皓）粤瞎），掣：ＫＩ瞎）掣
ｃＩ亡

来代替，则可实现网格密度的控制，上式中的Ｐ，Ｑ即为控制节点分布的参数。

ｄＣ

ｄ，７

ｄ‘

萼掣“：皓）掣，唑掣一础）訾
ｏＣ ０‘
ＯＣ

Ｏ‘

系数Ｋｌ，Ｋ２用于控制正交网格深入区域内部的程度。
说明：

（１）插值公式中的后两项导数项是为了保证区域边界处的正交性。

（２）当ｈ．＝２叩３—３叩２＋１，ｈ２＝－２ｒ／３＋３７／２时，可同时实现网格在两个边界处的加密
而当ｈ。＝１一，７，ｈ：＝ｒ／（线性情况）时，生成的网格比较均匀。 （３）当Ｐ＜１时，网格线向Ｊ＝ｌ的边界处加密 当Ｐ＞Ｉ时，网格线向Ｊ＝ＪＭ的边界处加密 （４）Ｑ用于控制善与刁的关系偏离线性的程度，对网格影响不大。

（５）若将红表达式中的叩（２孟矗）用下式替代，也能够实现网格的加密

其中：成为大于ｌ的常数，若玎接近于１，则可使网格在叩＝０和玎＝ｌ（即Ｊ＝ｌ
和Ｊ＝ＪＭ）处加密。

２．３网格边界的正交化技术 网格的正交性，特别是在物理区域边界附近的正交性，对于整个流场的计算精度 是至关重要的。实现网格正交化的方法有许多种，比如：在微分方程网格生成方法中，
改善正交性的途径是选取合适的控制函数，根据网格的正交条件ｇ．，＝０给出了控制

函数及控制方程，并实现了网格在边界上的近似正交”４”。：在网格生成过程中引入
Ｈｅｒｍｉｔｅ插值，通过求解一个泊松方程和一个拉氏方程来实现网格的正交化”。。 在上述方法中，网格的正交性算法隐含在整个网格生成过程当中，并且与网格的

河海大学博士论文

生成是同时实现的，此外由于控制函数是按正交性的要求而设定的，故无法同时实现
网格的加密。还有一种正交化的方法是通过移动边界点来实现的，但是当边界曲线的 曲率较大时容易引起网格的变形“…。 文献【６］介绍了一种网格边界正交化的代数方法，本节以该方法为基础，对其进

行了改进和完善，给出了更具有一般性，容易实现且独立于网格生成过程之外的代数 网格边界正交化方法，通过控制参数的设定很容易控制网格正交化的程度和范围。 ２．３．１正交化算法
．、ｊ

本算法的基本思想是给定边界点的分
一——

布，将区域内部网格节点向对应边界点的法 矢量方向移动，从而实现网格在边界区域附
近的正交化，下面以Ｊ＝１边界为例说明该正 交性算法。如图２．３，利用超限插值方法生
ｖ、ｍ

、、、ｊ．＾
。，．

一一。Ｊ－
（卜Ｉ，１）

，

、．蕊、～
，ｉ

成的区域内初始网格点为ｘ：，Ｙ：，首先求
出边界曲线Ｊ：ｌ在（ｉ，１）点处的法向矢量 月，则可以得到初始网格点在ｎ上的投形点
图２．３网格在边界处的正交化

《，Ｙ：。将初始网格点向投影点移动，则
可以实现网格在边界处的正交性，移动的程度可通过一个权重数∞来描述。

ｘ。＝（１一％）ｘ：＋％ｘ； ＹⅣ：（１一曲Ⅳ）ｙ：＋Ｄ口Ｙ；
显然，权重出，，应满足０≤∞，，≤ｌ 甜。＝０表示网格点没有移动 ∞。＝１表示网格点无阻尼地移动到对应边界点的法线上。

在算法的具体实现中需要解决两个问题： 第一个问题是网格点在对应边界点法线上投影点的计算。通过解析几何的方法推得投 影点的计算分式为：

ｘ：＝队，砘）＋吉ｈ＋吲他＋％，） ｙ：＝一÷ｘ：＋ｙ。＋古ｘ。
其中Ｉ．为边界线在点（ｉ，１）处的切线斜率

第二个问题是权重甜。的给定，本文采用文献［６］推荐的公式，并进行改进，给出了更

河海丈学博士论文

一般形式的权重计算公式。

ｏａ，，＝ｅｘｐ｛枷卜老）－２＿ｌ】卜，烈‰酬＿。一ｐ】。
其中：ｄ，ｊ＝√（ｘｇ—ｚ叫）２十（ｊ，Ⅱ一ｙ，｛）２
Ｐ，Ｑ称为阻尼参数，Ｐ值控制网格正交性在Ｊ方向上衰减程度，Ｐ值越大，则沿
Ｊ方向向区域内部传播的正交性衰减越快，Ｑ值控制网格正交性沿Ｉ方向的衰减程度， Ｑ值越大，则表明只有在边界的中心区域网格是正交的。 在解决了上述两个问题之后，即可实现网格的正交化。

２．３．２正交网格生成算例
为检验正交算法的可靠性，给出了正交网格生成的几个算例…，为便于比较，在

图２．４中给出了利用双边界法，但未采用正交化技术生成的初始网格。从图２．４、２。５ 中可以看出，使用正交化技术能够实现网格在边界处的完全正交，区域内部正交化的 程度可根据情况进行调整，周时网搐的加密与正交化可同时进行。

鐾蓁 。。冀熏
量垂肇薹蘩蕊浏
。一一 ０ ２０ ０．４０

三三三至至量耋善摹？’膏誓◆：曩

审器至毒害耋薹辇ｊ尊毫享三量 。”藿薹塞塞蠢蓑ｉ誊：摹霉 蔓芝三三主三三妄；｝■÷如ｉ、‘‘：一．ｊ

’∞一一…～——一一一一一…

～——

?一Ｏ¨一－一
１∞０．００

Ｏ瑚０．８０

４∞０．４０

～．一一
Ｏ∞

．

０∞

１∞

（１）

（２）

图２．４网格生成算例１一未使用正交化技术，２－－ｆｌ！ｈｑｉ交化技术

《三～
图２．５正交网格生成算例
培

Ｉ∞

河海大学博士论文

２．４分区结构网格技术 ２．４．１分区结构网格的一般特点
实际流动问题往往具有非常复杂的几何物理区域，传统的处理方法是通过曲线坐 标变换将复杂的物理区域转化为规则的计算区域，由此导致变换后的方程组更加复

杂，方程中各项的物理含义也模糊不清，计算速度和精度大大降低，目前流行的处理 方法是采用非结构网格技术和分区结构网格技术，由于非结构网搐的节点和单元分布 是任意的，因此具有非常优越的几何灵活性，可应用于任意形状流场的模拟，其主要 缺点是需要更多的计算机内存和ＣＰＵ时间，结构网格中已经成熟的算法在非结构网 格中无法应用，同时在复杂粘性流场中仍存在一些尚需解决的问题。分区结构网格的 基本思想则是把几何外形复杂的计算区域划分为若干个尽可能规则的子区域，在每个 子区域中独立生成网格并进行流场计算，其最大的优点是可以使用结构网格中非常有 效、可靠的求解技术，适合于并行计算，计算效率高，因此分区网格技术越来越受到
人们的重视。 分区网格系统一般分为相邻子区有重叠部分的重叠网格（ｏｖｅｒｌａｐｐｉｎｇ ｇｒｉｄｓ）１８
２＇驯

和相邻子区无重叠网格的拼接网格【３２】（ｐａｔｃｈｅｄ ｇｒｉｄｓ），分区拼按网格系统又分为相邻 予区网格线在交界面上对接和非对接两种情况，非对接情况比对接情况要复杂，但其 应用的空间和优势也更大，因为可以实行逐块加密而实现每一区的网格最优，同时允 许沿交界面移动网格，因此交界面处网格线非对接的分区结构网格拥有非结构网格的 几何灵活性和局部加密能力。同时可以利用可靠的结构网格求解技术，使计算过程更
加简单、高效．因而对于大多数流动问题都是非常有吸引力的。

２．４。２分区结构网格及其数据结构 分区结构网格的优点从图２．６中可以看出，如果采用单区网格技术，则不论是利 用代数生成方法还是解椭圆方程生成方法，都将会产生网格线在物面处非正交、网格 扭曲、尖点多、网格走向与流动方向不一致等缺点，从而严重影响计算精度及效果，
若将流场按照其几何外形和流动特点划分为图中所示的Ｉ，ｌＩ。ＩＩｌ区，则问题迎刃而
解。

为了保证在对分区结构网格数值计算时能够使用 单区结构网格求解器，在分区结构网格中需要采用两 种数据结构，即在每一区中的规则数据结构（ｆ√）和
分区交界面的数据结构，典型的分区交界面如图２．７。
一

在规则数据结构中，整个流场的网格节点和计算节点
按逐块和逐列的方法利用一维数组来标记和存放。

ｊ
９

河海大学博士论文

，＝Ⅳ。ｋ。‘＋（ｆ－１）Ｎ；＋Ｊ

ｋ＝１，．，ＮＢ，ｉ＝ｌ，．．，Ⅳ？，＿，＝１，．，Ｎｊ

其中：ＮＢ～总的分区数；Ⅳ？，川ｋ一第ｋ区网格在Ｉ?Ｊ方向的网格节点数
Ⅳ：＝Ｎ？×Ｎ：一第ｋ区网格的节点总数。
此外，在交界面数据结构中，需要对交界 面左右两侧控制体中心以及交界面上所有网格 节点进行序号标记。由于交界面上两侧子区的 网格线是非对称的，故同一个控制体积在交界
面上可能拥有几个控制面，如图２．７中典型的
ｒ，ｔｒ冉 工 月 月 卫

舡卉Ｊ

分区交界面，控制体积Ｌ在交界面上有３个控 制面．因此需要确认交界面上每一个网格面所 对应的左右两侧的控制体。关于分区结构网格

的数值计算方法在第三章中单独给出。
２．５小结

图２．７分区网格交界面

本章以通用形式的代数网格生成方法为基础研究了网格的边界正交化技术，提出 了一种简单的、能够独立于网格生成技术之外的边界正交化方法，给出了详细的数学

模型并编制了计算程序：通过网格生成算例表明，该方法能够大大改善网格在边界处 的正交性。此外对分区结构网格技术进行深入的讨论，并建立了非对接分区结构网格
的数据结构形式。

河海大学博士论文

第三章非正交同位网格系统中的ＳＩＭＰＬＥ算法
３．１交错网格系统与同位网格系统
数值计算的基本思想就是将求解区域划分成许多小的区域，即网格，在每一个网 格中定义求解变量并进行计算，网格中计算变量所在的点称之为计算点，计算点的布 置方式有多种，可以在网格节点，可以在网格中心点，也可以在网格线中心点，在图 ３．１中给出了三种情况计算节点对应的控制体。

——一～一一…一１一。一一＿’一
一一～～一＋÷—￥＾’＇｝～一一一Ｊ－２Ｌ２』‘“‘一—— ∽ｊ臻：麓
７

ｊ－：Ｅ， ７∥，■

＿■一ｊ—ｊ．

１『一再，７

～Ｍ

７，≯誊ｎ，≯ 。；：Ａ ｊ够

幽３．Ｉ计算节点及对应的控制体

（ａ）交错网格系统

（ｂ）同位网格系统

图３．２非正交网格系统中变量的布置形式

如果求解变量定义在不同的计算节点且拥有不同的控制体，则称为交错布置，对 应的网格系统称为交错网格（ｓｔａｇｇｅｒｅｄ ｇｒｉｄ），采用交错网格系统的主要目的是解 决Ｎ—ｓ方程中压力梯度项及连续方程的离散困难，同时避免出现波动的压力场和速度 场１４”，但是在进行编程计算时，需进行各物理变量的定位及复杂的插值运算。 如果求解域中所有的物理变量均定义在同一个计算点上且拥有相同的控制体，则 称为同位布置，对应的网格系统称为同位网格（ｃｏｌｏｃａｔｅｄ ｇｒｉｄ）。典型的交错网格
系统与同位网格系统见图３．２。

显然，只要能解决数值解的波动问题，同位网格在编程的简易性和通用性方面更 具有优势，因为所有的变量拥有同一个控制体，它们对应离散方程中的许多项是～样 的，计算量及存储量均减少，但是由于网格面上没有布置速度变量，在计算通过控制 面的通量时其速度需要由相邻两个计算节点上的速度插值得到，从而在程序中增加了

河海大学博士论文

更多的插值运算，尽管如此，同位网格系统在应用于复杂计算区域的非ｊ＿匕父嗍格和非

结构网格中更具有交错网格无法比拟的优判２ ５１。
３．２

Ｎａｖｉｅｒ－Ｓｔｏｋｅｓ方程及其在非正交同位网格系统中的离散方法

３．２．１积分形式的Ｎａｖｉｅｒ－Ｓｔｏｋｅｓ方程
应力形式的Ｎ．Ｓ方程的一般形式为【３３，８４ｌ：

判＋旦照趔：反一堡＋鱼
Ｏｔ ａｘ．

…苏．

（３．１）

ａｘ．

为便于理解，写出分量形式

掣ｏｔ埘?嘞＝六一警＋Ｃ等＋等］ 掣埘?删＝六一警＋（鲁＋鲁］
办
ｌ
ｄ．
’

＠ｚ，

ｄ，Ｊ

＠ｓ，

掣埘?咖＝正一警＋［等＋等］
其中：Ｐ，＝一Ｐ＋２∥警一姜脚ｖ矿为粘性流体中的法向应力； 戚
ｊ

慨ａ，

。＝∥（等＋斟“一烨一Ⅲ粘性流体中的切向鱿沩切向应
力所在面的法线方向，Ｊ表示应力本身的方向。

代入（３．１）式，并注意到研究对象是不可压缩的牛顿流体：

掣＋掣＝∥一若＋毒卜考］
西＝一ｇｒａｄ（ｐｇｈ）。
其中：ｈ＝ｚ一‰，ｚｏ为参考点的垂向坐标。

ｃ。．ｓ，

其中：，为重力加速度ｇ在三个坐标轴方向上的分量。显然，若ｐ＝ｃｏｎｓｔａｎｔ，则

则Ｐ＋触＝芦
从而可以将质量力项合并到压力项中处理，为方便起见，在以下各式中，将声仍 记作Ｐ，从而可以得到微分形式的Ｎ?Ｓ方程：

河海大学博士论文

——————卟一； 善十毒卜善ｊ 融，Ｉ’良，ｊ
８ｔ

ａ（∥，），ａｂ．“，）
瓠ｉ

（３．６）

苏。

将上述方程对控制体Ｑ积分，并利用Ｇａｕｓｓ定理即可得到积分形式的Ｎ．Ｓ方程：

昙￡艘地＋肛，矿．蕊＝一ｆＰｉ，?赢＋ｐ｝?蕊

（。。，）

其中：ｉ（ｆ－１，２，３）为三个坐标轴方向上的单位矢量，等式左端第一项为随时间变化的
非定常项，第二项为对流项，等式右端第一项为包含质量力的压力项，第二项为粘性 力项，又称扩散项。 积分形式的连续方程直接写出：

杪?蕊＝ｏ
利用有限体积法求解流场的原型Ｎ－Ｓ方程即为（３，７）、（３．８）式。

（３．８）

需要说明的是，关于有限体积方法的特点在文献ｆ３０，３３，５５】中有非常详细的介绍，在 这里就不详细展开，直接利用有限体积方法推导非正交同位网格系统中的离散模型。
３．２．２

Ｎａｖｉｅｒ－Ｓｔｏｋｅｓ方程的离散方法

３．２．２．１网格点的布置
本文采用非『Ｅ交同位网格系统，计算点布置在网格中心，对应的控制体即为网格
本身。

３．２．２．２动量方程的离散 以２Ｄ情况为例，推导非正交同位网格系统中动量方程离散形式，离散过程中涉
及到的计算节点、控制体、控制面等信息参见图３．３
ｚ．
～

～～

ｆ～
酽
●Ｅ

Ｃ

，ｉ

．

一～～—鑫＆兰
?

；
ｆ

ｓ

图３．３计算节点与控制体

河海大学博士论文

１）对流通量Ｉｐｕ，Ｖ?ｒｉｄｓ
仅以控制面中的ｅ面为例说明，其它控制面的形式完全相同，速度分量“，１１－－Ｉ殳 符号≯替代。同时为保证离散形式具有二阶精度，并避免产生波动解，采用中心格式
与上风格式联合的方法，具体方法在３．４中介绍。

①显式对流通量（中心差分格式）

ｋ。）＝…＝Ｌ却旷?ｉｄｓ＝Ｌ矽?ｈｄｓ。丸ｚｍｅ＇以
其中：

（３?９）

卅，＝￡ｐ矿－赢＝ｂ矿?二ｋ ｓ。＝ｐｐ。?以＋ｓ７丸）为通过ｅ面的质量通量；
ｓ。，Ｓｙ为控制面积Ｓ。沿Ｘ，Ｙ方向的投影面积： 以，≯，为速度在ｘ，Ｙ方向的分量：

以为控制面ｅ中心处ｅ点的速度，由线性插值确定，如图３．４所示。

吮吨＋（觌ｏ。－－。）＋㈣（ｙ。Ｍ）

锨埔㈣）＋㈣。＂（新一丑忙ｘ。，） ＋［（詈］。丸＋（考］，ｏ一九）］。。一
值的方法确定。

（３．１０）

其中，相邻计算节点连线与ｅ面的交点ｅ处的坐标、递厦及远厦梯度ｉ电辽线住艳

ｘ。．＝ｘＥ，，ｌ。＋斗（１－２ｃ） ｙ，＝ｙ，五＋¨（１一五。）

丸咆＋罢∽－－。ｐ等帆诫．）
砘”办（１－－五ｅ）＋罢（Ｘｅ－－Ｘｅ，）＋詈（圹虬）

（罢］。＝（罢］。以＋（罢］，‘（１吐）

ｃ。?Ⅲ

河海大学博士论文

其中：

丸＝Ｉ÷兰ｉ为线性插值矢量；计算节点上的速度梯度计算在后面单独给出。

图３．４控制面上的变量插值

②隐式对流通量（上风差分格式）
ｆ

ｋ。尸＝上厕蕊：卜办…
Ｉｍ。‘屯…

矿（ｙ?扦）。＞０

矿舻．’ｎｌ＜０
（３．１２）

＝ｍｉｎ【ｍ。，０ Ｉ毋Ｅ＋ｍａｘ（ｍ。，ｏ）≯Ｐ
＼ ／

则总的对流通量采用下式计算

ｃ。：ｋ。尸＋Ｉｃ）一一ｋ‘严ｒ

（３．１３）

其中【】中的项称为延迟修正项（ｄｅｆｆｅｒｅｄｃｏｒｒｅｃｔｉｏｎ），放在方程右端的源项之中， 即为离散方程中的Ｑ‘项，并利用上次的迭代值计算。采用延迟修正的原因是一方面 可使迭代数值解具有二阶精度，另一方面可以避免因单纯上风格式造成系数矩阵的非 对称性，从而导致迭代算法的发散。

２）扩散通量［∥挈ｉ淼
小Ｏ譬．

仍以控制面ｅ为例说明，动力粘性系数∥及速度”，分别用一股变量ｒ和矿替代，
于是扩散通量的一般形式为：

【Ｆｇｒａｄｑ）?ｒｉｄｓ

①显式扩散通量

河海大学博士论文

忙“广。＝ｆ，Ｆｇｒｎ彬－础＝（Ｆｇｒａｄｑ６?；Ｉ‘ｓ。

＝Ｌ伊＋－印ｆｆｓ’１

（３ １４）

上式中需要计算物理量在网格面上的梯度值（—兰）。，一般近似采用中心差分格式，即 佩
直接由控制体中心处的导数插值，只是需要对网格的非正交性影响而引入修正。

昏＊磬＝鳊ｒｅ ｒｅ—ｃ删烈隔吐，
出， 废，。
ｌ—１ ％一，尸｜

其中当网格在ｅ面处正交时，上式第二项等于０。

②隐式扩散通量

ｋ。尸＝￡（咧螂肛Ｌ（豢）。辽

２等¨¨邓。（删∥（禹以） 一Ｃ¨ｓｏ（ａ彘＃ａｘｍ＋雾知”卜蹦删删“ａ隔ｒＥ－ｒｐ吨）（３
其中：ｆ，２ｒａｄｄ，）“有网格中一Ｉ）处的梯度线性插值，

１５）

禹吨。与掣
则扩散通量

（，。。一＿ｙ。）ｉ＋（】‰一石。）＿７

扛乏面＝

Ｆ：

（Ｆ—ｍｒａ＋ｋ一广’一ｋ“严ｒ
ｃｏｒｒｅｃｔｉｏｎ）

（３．１６）

上式［】中的项称为延迟修正（ｄｅｆｆｅｒｅｄ

放在离散方程的源项之中，即

为∥项，并利用上一次的迭代值计算。
３）源项计算
首先给出压力源项的计算

Ｑ’＝一ｐ?淼＝一．ｆａｇｒａｄＰ?№

（３．１７）

河海大学博士论文

对于“勰酣＝一￡驷＝一（新△Ｑ
具甲：Ａｆ２为控制体体积，任意形状控制体体积的计算方法见３．４．４ 对于ｖ方程，情形类似。
非定常项的计算方法如下：

（昙胁碱＝害＠？＋１－－４Ｕ～一）
＝４；“ｉ；１一Ｑ＝．

煦耻等卅，＝害㈣？叫）
此外，在对流项与扩散项的计算中，分别出现了对流项修正和扩散项修正，也作
为离散方程的源项进行处理，即

Ｑ‘＝ｋ‘ｘｐ ７一ｋ。尸

ｇ“＝ｋ哼…一ｋ“严
以“方程中的ｅ面为例给出具体的表达式：

饼＝川。【丸无＋办（１一ｔ）＋Ｘ。－－Ｘｅ，ＯＸ 一ｍｉｎ（ｍ。，Ｏ）屯＋ｍａｘ（ｍ。，Ｏ）绋

）＋警（咒一以）】 删
。

Ｑ？２

Ｌ（罢ｓ‘＋警ｓ’）～等ｅ警缸。＋芳匈。）。

最后给出由于网格非正交而引起的修正项

Ｑｇ＝呱删Ｈ嚣高词】［１蝇】
３．２．３离散形式的Ｎａｖｉｅｒ－Ｓｔｏｋｅｓ方程
将（ｇ．１３）、（３．１６）式代入（３．７）式，得动量方程的离散形式：

４办＋∑铂＝ｇ
ｆ

（３．１８）

其中：Ｉ＝Ｅ，Ｗ，Ｎ，Ｓ为与计算点Ｐ相邻的计算点

河海大学博士论文

。

喊

咿

鹏一ｋ

一 鹏一‰ 一 丛‰
娜
ｍ

一

卟

髓一ｋ

ＡＰ＝√；一Ａｆ—Ａｗ～ＡⅣ一ＡＪ

Ｑ一：钟＋Ｑ；＋彰＋Ｑ；＋Ｑ８．该式右端分别为压力项、对流项修正项、扩散
项修正项、非定常项以及非正交网格修正项对源项的贡献。 ３．３非正交同位网格系统中的压力校正方程及其离散
由于隐式离散方法对时间步长没有严格的限制。故正如上节中介绍的那样，动量

方程采用隐式方法进行离散，为推导压力校正方程，我们将动量方程中的压力项从源
项中分离出来，以ｕ方程为例进行说明。

群砷＋∑群“，＝（Ｑ一饼）＋饼
，

（３．１９）

若动量方程的第ｍ次迭代值为“，，则

彳；“；’＋∑钟“？＋＝（Ｑ。一饼）＋钟
，

（３．２０）

”ｆ＝＿—一＋÷饼
（包一饼）一∑掣矿，
群
’

群一

（３．２１）
‘

’

记右端第一项为印同时咖一１１驷－－（孰垃
虿亥ＪＰ
△Ｑ，印、
（３．２２）

河海大学博士论文

显然，速度“：’满足动量方程，但不满足连续方程，为此需要对该速度进行修正，
以满足连续方程（３．８），在ＳＩＭＰＬＥ算法中，是通过修ｉＥ压，Ｊ７场来实现的。

咖夥一等嗜Ｏｐｍ）
其中：

（３ ２３）

“；＝”；＋＋“；为满足连续方程的速度，“ｊ为校正速度。
Ｐ”＝Ｐ＋Ｐ’为相应的压力场，Ｐ’为校正压强。

若将（３．２３）式代入连续方程塑竺卫：ｏ，可得压力泊松方程

砉【等（ｉＯｐｍ ｍ＝‘ＴＯ（ｐＨＴ＂）】，

（３．２ａ）

通过直接求解压力泊松方程（３．２４），即可得到满足连续方程的速度场Ⅳ？，但此 时的速度场及压力场并不满足动量方程，因此需要将迭代进行下去，上述的迭代过程，

即通过压力场校正来满足速度场，称为外部迭代，而内部迭代是指线性代数方程组的
内部求解迭代，这种直接求解压力泊松方程来确定压强的方法称为Ｐｒｏｊｅｃｔｉｏｎ
ｍｅｔｈｏｄ。

在实际计算中最常用的一种方法是压力校正法，（３．２３）式减去（３．２２）式可得
到校正压力Ｐ’与校正速度”’之间的关系

啦即筹（参
值之后控制面ｅ上的修正速度与修正压强梯度依然保持上述关系，即

（３．２５）

为了应用连续方程（３．８），需要知道控制面上的速度修正和压强修正，可以通过 插值的方法来确定。由于在修正速度的插值计算中包含了修正压强梯度的插值，故插

咖∥一等《）。
∑Ａ，ｕｌ
中，该项忽略不计。
即

（３ ２６）

式中ｑ＝一々，但由于在求解压力校正方程之前，珥未知，故在ＳＩＭＰＬＥ算法
虬一万‘玄’ｅ 吩一等ｃ拳
（３．２７）

河海大学博士论文

同理

吩一等（答

（３ ２８）

将（３．２７）、（３．２８）式代入离散型连续方程ｍ。＋ｍ。＋ｍ，＋ｍ。＝０，并注意到由动量方 程第ｍ次迭代解出的速度场并不一定满足连续方程，因此

捕：＋确：＋ｒｈ：＋廓：＋△确：＝０

其中△祝为迭代速度场对应的质量通量误差。

即（ｐｓｕ’）。一（ｐｓｕ’）。＋（ｐｓｕ’）。－（ｐｓｕ’），＋Ａｔｈ：＝０
在非正交网格系统中，修正质量通量计算采用下面的公式：

（３．２９）

咖（∥卜一（ｐＡｃ２ｓ）。（刁１。（拳

＝一器咖喊啊）－（倒ｐ从¨）】 鼽（拳ｚ南№ｈ卜（ｇｒａｄｐ’Ｈ瑶刊】
（ｇｒａｄｐｌ）。＝（ｇｒａｄｐ‘）Ｅ丑＋（ｇｒａｄｐ’）Ｐ（１～２０）

（３ ３０）

（ｇｒａ印’）。．（瑶一斥）＝【（譬Ｌ）。丑。＋（芒Ｌ），（１一九）】ｏ。一ｘ，）

“（婴）。丑。＋（孕），（１一五。）慨一ｙｐ） 们 却
应当注意，式（３．３０）的第二项是由于网格的非正交性而引入的，当网格非正交 性较弱时，该项影响不大。但当网格是强非正交的，则必须计入该项的影响，将另外

三个控制面上的质量通量代入（３．２９）式，并只进行一次修正，可得非正交网格中的 压力校正方程为：

爿；ｐ；＋∑彳，｜Ｄ；＝一△廊：＋Ｑ：

（３．３１）

舯班一（等）。而１
３０

河海大学博士论文

俨一（等）。瓦１丽
爿：＝＿（、ｐＡ雒ＤＳ）。百杀

班一ｃ等，，赤
爿ｆ＝一∑钟 ａ；＝∑（ｇｒａｄｐ。）；（不一‘），ｋ＝ｅ，ｗ，ｓ，ｎ，Ｋ＝Ｅ，∥，Ｓ，Ｎ
△痢：＝一（嘶：＋嘲：＋，ｆｌ：＋嘲：）”１为前一次迭代速度场对应的控制体质量通量。
至此，给出了非正交同位网格系统中ＳＩＭＰＬＥ算法中的动量方程和压力校正方程 的离散形式，下面将ＳＩＭＰＬＥ算法的计算步骤简述如下：

（１）将ｎ时刻的计算值“？，ｐ”作为初始值计算ｎ＋ｌ时刻的∥、ｐ∥；
（２）求解动量方程得到“？‘（该速度场不满足连续方程，要通过修正压力场来修正速
度场）

（３）求解压力校正方程得到Ｐ’

（４）修正速度及压力，得到满足连续方程的速度场和压力场 “●＝“？＋＋“’，Ｐ”＝Ｐ＋Ｐ。 （５）重复（２）．利用妒和Ｐ”作为对∥”，Ｐ”’的修正，计算直到所有的修正量足够
小。

（６）推进到下一个时间步长。
（２）～（５）称为外部迭代。

３．４非正交同位网格系统中的专题讨论 ３．４．１边界条件 ３．４．１．１壁面边界条件 壁面边界条件通常采用非滑移边界条件－（Ｄｉｒｉｃｈｌｅｔ条件），即流体速度等于壁面速 度，在有限体积方法中一般利用法向应力等于零的Ｎｅｕｍａｎｎ边界条件，以二维情况
的南边界为例进行说明，见图３．５

河海大学博士论文

Ｎ

—７一—７７＿，一，—！，——≥—＿———７—７—７——７———，——．ｒ ／／／／／７／，／一。，／，／／７／ ／／／／／
图３．５壁面边界条件

“＝２∥（罢）。＝０ 砂
扩散通量的贡献为

（３．３２）

由于壁面边界是不可穿透的，故对对流通量项没有贡献，只有对扩散通量的贡献。对

吩Ｉ，Ｔ、ａｄｓ＝ｊ砖㈣ｓ，（孰
Ｙｖ—ｙｓ

＝以Ｓ，盟 ：盟０，一“。）：４，一“，ｓＳ。ｔ
’

（３．３３）

其中Ａ，６，Ｓ。６分别为边界条件对离散方程中系数及源项的贡献。 对于３Ｄ情况，与２Ｄ情况相类似，下面以底面（Ｂｏｔｔｏｍ）为例．

对“方程：由于在底画上ｒ。＝娑＋＿Ｏｖ：０，所以 洲ｏ＇ｘ

‘４＝ｆｐ，＋ｋ这 ＝Ｉ‘ｒ，瓠

＝∥。￡（笔＋警产。
Ｍｊ，挚，

Ｍ＆警

河海大学博士论文

２警旷－以Ｌ‘－－“ｓ
对于ｖ方程，类似可得：

（３ ３４）

ｆ．。＝警Ｖ，一警Ｖｓ
其中，底面对Ｗ方程中的扩散通量没有贡献

（３．３５）

３．４．１．２动量方程中的入１３及出１３边界条件
由于入口边界处的速度是已知的，故对流通量可直接计算

ｍ，＝ｐＶｎｓ，＝ｐ吣，＋＂ｖｓ，Ｊ
如果采用隐式离散格式，则入口边界对对流通量的贡献为

（３．３６）

Ｆ。＝ｍｉｎ（ｍ，，ｏｈ。．＋ｍａｘ（ｍ，，ｏ－，

（３．３７）

其中：第一项为入口边界对ｕ方程源项的贡献，ｍａｘ■，，ｏ）为入１３边界对离散方程中系
数一．的贡献。入口边界对扩散通量的贡献可类似于一般情形扩散通量的计算，只是 速度梯度采用单向插值计算：
Ｆｄ

ｆＦｇｒａｄｕ‘赢＝ｒ（黔，＝号¨训

㈣。ｓ，

出口边界条件对对流及扩散通量的影响，其处理方法类似于入口边界，只是出口

边界要放在距离入口较远的下游，同时出口边界上的速度采用一阶上风近似，即 庐ｆ＝办，≯＝“，Ｖ，‘ｗ ３．４．１．３压力校正方程中的入口及出口边界条件
压力校正方程中的边界条件同样需要特别注意。若通过某边界的质量通量是给定 的，则压力校正方程中的质量通量修正等于０，即采用零梯度的Ｎｅｕｍａｎｎ条件．

在出口，若入口质量通量是给定的，对于定常运动，则出口边界上的速度可采用 外推法而获得（即“，＝＂，），但是需要对由外推法得到的出口速度进行修正，以保证
各个流场内的质量守恒（即出口质量通量与入口质量通量相等），出口速度修正方法

见图３．６。利用修正后的速度进行下一次的外部迭代，同时出口边界上的质量通量修
正值为零。

为保证解的唯一性，需要选取一个固定点的压强作为参考压强，所有其它点的压 强修正计算都要减去同一个参考点的压强，但是需要注意，在分区结构网格中，参考

点只能取一个，不能每区各取一个。

河海大学博士论文

图３．６出口边界的速度修正

３．４．２延迟修正技术（ｄｅｆｆｅｒｅｄ ｃｏｒｒｅｃｔｉｏｎ）
在有限体积方法中．网格控制面上的变量值都需要通过插值来确定，若利用高阶 插值的方法会导致计算量迅速增大，如对二维情况，利用四阶辛普生（Ｓｉｍｐｓｏｎ）公式， 则每个通量需要计算１５个节点上的变量值，并且每个控制体对应的代数方程包含２５ 个节点上的变量值，对于非线性方程组来说，其求解工作量是非常大的。

我们知道，利用前一步的迭代值显式计算高阶通量是简单的，而隐式的低阶通量 仅需要与计算节点最邻近的节点上的变量值，计算工作量小。若能将二者结合起来，
则将同时兼顾计算的精度与简单性，这种思想由Ｋｈｏｓｌａ
出【８卯，具体形式如下：
ａｎｄ

Ｒｕｂｉｎ于１９７４年首先提

Ｆ二＝Ｆ？９＋ｔＦ，ｉ—Ｆ？ｓｍｐｌ

ｌｏｌｄ

ｔ３．３９）

其中：Ｆ州为某种低阶隐式离散格式； ｒ州为某种高阶显式离散格式； 式中第二项由前一步的迭代值计算，与第一项的隐式部分相比是很小的，因而显 式处理不会影响计算的收敛性。这里将高阶与低阶离散格式联合的目、的有两个：一是
避免出现物理上不合理的波动解：二是改善某一迭代求解算法的收敛性［３２，３３１。

河海大学博士论文

３．４．３控制体中心的梯度计算
在对流通量，扩散通量及压力源项计算中 计算方法如下，以扩散通量为例：
Ｆｅｄ

都需要知道控制体中心Ｐ处的梯度

ｔ莩（黔。Ｊ

（３．４０）

其中交界面ｅ处的梯度可以通过控制体积中心点上的梯度插值得到

（考］。＝［考］，屯＋［等］，ｏ以，
控制体积中心处梯度可利用Ｇａｕｓｓ定理来确定

ｃ。．ｔ?，

Ｅ等ｍ肼ｈａ，ｓ
如缸． 垆‘一‘
△Ｑ△Ｑ

ｚ莩箬２面１融＋丸Ｓ。＇－＃ｉｗ吖一∥）
３．４．４三维任意形状网格控制体体积计算

（３．ｔ２）

任意形状的控制体积计算的原理是Ｇａｕｓｓ定理， 即将体积分转换成面积分
由于

击ｖ（ｘｉ）＝１，故控制体体积

△Ｑ＝肛＝ｐＶＧ彳妞＝Ｉｘ７ｄｓ ｚ∑ｔＳ。ｘ

（３．４３）

其中ｃ为控制体拥有的控制面，疋。为控制体的面积矢量在ｘ方向上的分量，即

控制面积Ｓ。在ｙｚ坐标面上的投影面积

ｓ。＝ｓ。。ｉ＋ｓ。ｙｊ＋墨。ｉ
下面以空间六面体为例，ｃ＝ｌ、２、３、４、５、６，工。可取ｃ控制

面上四个顶点Ｘ坐标的平均值计算，如图３．７，Ｓ。１为第ｃ个控
制面在ｙｚ坐标面上的投影面积
图３．７控制面投影

ｓ。Ｊ＝圭Ｉ舅一ｉ）×Ｅ—ｉＩ＝吉［（ｙ．一ｙ，）（Ｚ２－－Ｚ４）一（２＇Ｉ－－Ｚ３）（ｙ：一ｙ。）】
应为负号。

（３．４４）

将每个控制面上的Ｘｃ和Ｓ。‘计算出来之后，利用（３．４３）式即可计算出该控制体体积。 注意：若控制面外法线方向与Ｘ坐标轴方向的夹角大于９０。，则（３．４３）式中该项前

河海大学博士论文

３．５分区结构网格中ＳＩＭＰＬＥ算法的实现

动量方程的离散形式为：Ａｐ办＋∑Ａｋ≯ｋ＝Ｑ尸，女＝Ｅ，∥，Ｎ，Ｓ，在单区网格内部，
ｌ

每一个控制体拥有４个控制面，而交界面处的控制体则可能拥有多于４个的控制面， 故实际计算时每个网格面上的对流通量和扩散通量需单独处理，然后补充到交界面两 侧控制体的离散方程系数中去，如计算节点上对应的ｅ面为分界面，则取Ａ，＝０，在 离散方程的系数Ａ，中补充交界面的贡献一Ｒ，Ａ。，，Ａ马，同理，若控制体Ｒ：的∥面

为分界面，取Ａ，＝Ｏ，在对应的系数Ａ，中补充交界面的贡献Ａ。。 在计算对流通量时，需要知道物理量在控制面上的值，对于正交网格可直接有网
格面两侧的计算点插值得到。但在一般情况下，交界面上的网格线是非正交的，如图 ３．８所示，即，与，’不重合，因此交界面上每一个网格面的变量不能由两侧计算节点直 接插值得到，而需要进行修正，即

办叫＋（》“，一）
ｂＪｏｋ＾

?Ｎ

／
，

?

其中：

，

一一’一７Ｌ——’－

力＝≯．－＋≯。（１一五，－）

“？≮～
图３．８

一…ｉ。————；ｉ

学），＝（缸”尝Ｍ，１）
以上两式均由网格面两侧计算节点线 性插值确定。
非对接分区交界面

Ａ，＝；｝—鲁为线性插值矢量，
‘

ｈ一■ｌ

计算节点的梯度值采用下式计算

。耠挲
靠
ｙ

其中Ｓ：包括该控制体在分区交界面上的所有网格面。

上述格式可以保证交界面上的通量计算具有二阶精度。 在计算交界面上扩散通量的贡献时，需要知道交界面上的梯度矢量，若采用对流 通量的处理格式，则将出现高阶偏导数从而增加了计算工作量。常用的方法是直接采

河海大学博士论文

用中心差分格式，即由控制体中心的导数线性插值，同时对交界面处网格的非正交性 进行修正。

ｃ等小ｃ考）』２静焉一ｃ等，产ｃ音专“，，
其中当网格在交界面处正交时右端的第二项等于０。

非正交网格中的压力校正方程的离散形式为：爿；ｐ：＋∑彳，ｐ，，ｉ＝一△肌；＋鲜
其中需要注意的是对交界面上各个网格面上的质量通量进行修正，以，控制面为 例：

妒（棚＝＿（脚）『（扣挈，

一器（扣小础如ｄｐ＇Ｍ¨）］
其中：（！兰），ｚ＝＿！＝－ｉ［（ｐ：一ｐ：）一（ｇｒａｄｐ＋），（乓一五）】 靠（ｒＲ一丘妒
（ｇｒａｄｐ’），＝（ｇｒａｄｐ’）Ｒ＾＋（ｇｒａｃｌｐ’）￡（１—４）

（ｇｒａｄｐ‘），．（不一五）＝【冬）。丑＋（冬）。（１一删‰一吒） 盘 劣 ＋【冬）。＾＋（婴）。（１－丑）舰一ｙＬ） 卯
ｃｒｙ

应当注意，上式的第二项是由于网格的非正交性而引入的，当网格非正交性较弱时， 该项影响不大，但当网格是强非正交的，则必须计入该项的影响。 ３．６非正交同位网格系统中ＳＩＭＰＬＥ算法求解算例 ３．６．１倾斜方腔顶板驱动流 顶板移动的方腔内流场计算常被用来检验数值算法的稳定性和精度，为了便于比 较，本文首先计算了正交网格下方腔的内流场，网格数为８０×８０，迭代过程中速度 及压强的欠松弛矢量为０，８，０．８，０．３，给出了速度分布曲线，并与文献［３３】的结果进 行比较，非常一致，见图３．９。为了检验非正交网格系统下的ＳＩＭＰＬＥ算法，对方腔 壁倾斜角为４５。时的倾斜方腔内流场进行了数值模拟，网格数仍为８０×８０，见图

河海大学博士论文

３．１０，但迭代过程中速度及压强的欠松弛矢量为０．７。０．７，０．２。在图３．１ｌ中给出了雷

诺数Ｒｅ＝１００时的斜方腔中纵剖面及中横剖面上的速度分布曲线，其中Ｒｅ：堕，“为
，，

顶板移动速度，，为方腔长度，ｙ为流体的运动粘性系数，图３．１２，３．１３中分别给出 了Ｒｅ＝１００和Ｒｅ＝１０００时对应的流线矢量图和压力等值线分布图，可以看出，在斜方

腔右下方产生的次涡中心随着雷诺数的增加向上移动。
’∞］
：—４

ｙ州
ｏ∞一

，．一一一
彳

…ｊｔ忡ｍＩ

…
？

五
『『－
／．

。。矿＼
４∞＿／
、＇、＼

¨、ｊ

７

。。ｊ／
Ｄ∞ ａ∞
ａ∞

了
１∞

＼
０∞ 口棚ｏ∞
ｔＭ

图３．９

二ｉ一一一～．一７型”一…一一．一一Ⅲ
删

１∞

８－－９０。，Ｒｅ＝１０００，方腔中纵断面及中横断面上的速度分布曲线： 图中离散点为文献【２５】中的计算结果

图３．１０倾斜方腔中的非正交网格

河海大学博士论文

一一㈠
图３
，

一

ｌ １

Ｂ＝４５０，Ｒｅ＝１００，方腔中纵断面及中横断面上的速度分布曲线

一…一三！型．
０帅 ０∞
１

２０

●∞

¨ ∞ ¨ ∞

ｏ ０ １ １ ５ ２

（ａ）

Ｒｅ＝１００

，

¨ ％ ¨ 舱

ｏ ０ ２

（ｂ）

Ｒｅ＝１０００

图３．１２ １３＝４５。倾斜方腔内的流线矢量曲线

河海大学博士论文

图３ １３

Ｂ－＝４５。倾斜方腔内的流线矢量曲线

３．６．２后台阶绕流
后台阶绕流是典型的流动问题，如文献［８６，８７］分别对三维后台阶的层流流动及后

台阶绕流的出流边界条件进行讨论。为检验分区结构网格算法的可靠性，对计算流体 力学中的典型算例一后台阶绕流进行分区计算，台阶的尺寸为入口高度５．２ｍｍ，台阶 高度４．９ｒａｍ，台阶上游长度２００ｒａｍ，台阶下游长度５００ｒａｍ，计算结果给出了不同Ｒｅ下

４０

河海大学博士论文

的回流区长度的变化曲线，其中Ｒｅ中的特征速度和特征长度分别为台阶上游的平均 速度及台阶高度，同时给出了Ａｒｍａｌｙ
ｅｔ

ａ１．的试验曲线【７８】和文献【７９】的计算结果见图

３．１４，可以看出，在Ｒｅ＜４００时，计算结果的一致性很好，而当Ｒｅ＞４００时与试验曲 线偏离较大，其主要原因是试验数据是在三维流动的情况下获得的。

０——?－，－－－－－ｒ，－－－－———————－———，——————－－－—－———－—－————————一，——
ｌ

２

３

４

Ｒｅ

５

６

７

８

图３．１４二二维后台阶回流区长度计算结果比较

３．６．３非对称平板间圆柱绕流 当研究位于两个平板间的圆柱绕流时【Ｅ”，如图３．１５，重点关注的是柱面压力与尾 涡，因而需要对柱面附近和尾部区域使用较密的网格，此时使用单区网格是不合适的，

要么整体利用非结构网格，或者使用分区结构网格．本文采用六分区网格对其进行了
数值计算，计算条件为：圆柱直径为Ｄ，圆柱中心至上下平板的距离分别为２．０８Ｄ和 ２，ＯＤ，圆柱上游距离为８Ｄ，下游距离为２０Ｄ，来流速度为ｕ，雷诺数Ｒｅ＝ＵＤ／，，， 这里仅给出Ｒｅ＝１００和Ｒｅ＝１０００对应的压力等值线，速度等值线和流线图，见图３．１６。 从图中可以清楚地看到圆柱后尾涡的形成和脱落过程，随着Ｒｅ的增加，流线的波动 逐渐加剧。

≈≯ 誓１＿㈡－］ 一 一 ｜；纛一
，■≈

一

一

图３ １５平板间圆柱绕流

４Ｉ

河海大学博士论文

圈３ １６

Ｒｅ＝１００（ａ）和Ｒｅ＝１０００（ｂ）对应的压力，速度等值线与流线图

３．７小结 本章在非正交同位网格系统中推导出Ｎａｖｉｅｒ－Ｓｔｏｋｅｓ方程的离散模型、压力校正方 程及其离散模型，对非正交同位网格系统中的边界条件、延迟修正技术、控制体中心 的梯度以及任意形状控制体体积计算等进行专项讨论，建立了非正交同位网格系统中 的ＳＩＭＰＬＥ算法，并以此为基础建立了非对接分区结构网格数值求解模型。利用上述 算法编制计算程序，对倾斜方腔驱动流、后台阶绕流、非对称平板间圆柱绕流进行计 算，通过与实验结果的比较表明，本文建立的数学模型是可靠的和正确的。

河海大学博士论文

第四章
４．１湍流研究的工程背景和意义

湍流模型

大多数工程中的流体力学问题都是湍流问题，应当用湍流理论进行处理，但是由

于实际工程中的流动甚为复杂，目前的湍流理论和计算技术都未达到可以大量解决实
际工程中流体力学问题的水平，人们为了满足工程应用与设计方面的要求，往往顾及

不上湍流场的分辨率，而仅对湍流影响作简单化处理。有些涉及到流体内部结构的工 程流体力学问题，不用湍流理论不可能较好地得到解决，如环境工程、水利工程中的
高速水流和挟沙水流、空气动力学及船舶动力学中的绕流问题等都是如此，因而近

２０年来，国内外都在研究如何用湍流理论解决工程中的湍流问题。其主要内容包括： ①利用目前比较成熟的湍流理论解释一些特定的流动现象， ②应用不同的湍流模型预测工程中的流体力学问题。， ③以湍流理论为依据，提出一些解决工程问题的新理论和新方法。 上述研究内容紧密联系，相互促进，一方面湍流理论已被广泛地用于解决工程问
题，同时在工程应用中又促进了湍流理论本身的发展，如直接服务于工程的湍流模式
理论等。

根据目前湍流研究的趋势，除对湍流基本结构大力进行研究外，结合工程中的问 题进行工程湍流研究，己成为国际湍流研究的重要方面，如湍流模式、湍流控制、湍
流测量、湍流与传质传热等。

因此，以工程为背景开展湍流研究，对于国民经济建设和科学技术的发展都具有 十分重要的意义，为此，必须以发展和利用已有的湍流知识解决实际技术问题为目标， 在工程中广泛应用湍流理论，为解决工程湍流问题提供实践基础，同时积极促进湍流
理论及其相关技术的发展，正如周培元教授指出的那样需要“积极开展流体湍流运动 的实验、理论和应用研究”１３ ７１
ｏ

４．２湍流模式理论 由于湍流的复杂性，目前要想获得复杂的计算结果几乎是不太可能的，需要引入 各种假定，这些假定是对实际湍流流动的某种近似模拟。所谓湍流模式就是指将真实 的湍流运动模化为人为设计的模型，这种模型的平均行为应与实际湍流统计平均行为
基本一致。
。

雷诺时均模型在总体上可划分为基于Ｂｏｕｓｓｉｎｃｓｑ假定的湍动粘性模型和基于 雷诺应力输运方程的二阶矩封闭模型，前者又可分为混合长度模型、一方程模型和二 方程（如ｋ．￡，ｋ．ｍ，ｋ－ｒ等）模型，后者则有代数应力模型（ＡＳＭ）和输运方程模型（ＤＳＭ）
两大类。

河海大学博士论文

Ｂｏｔｔｓｓｉｎｅｓｑ假设的基本思想在于，认为湍流应力的作用类似于层流应力

ｆ＝∥罢，但其比例系数为～流动变量——涡粘系数以，即ｆ；∥．—ｄ＿ｕ。其中＾是
ｕｙ

８ｙ

与，完全不同的量，它不是流体的一种性质，而是决定于湍流的流动状态，因而几在 不同流动状态，甚至在同一流场中的不同点，其值都很不相同。根据对涡粘系数一的 近似计算方法的不同，即可构造出不同的湍流模型，如零方程，一方程、二方程等湍 流模式。需要注意，在涡粘性模式中，雷诺剪切应力的方向与平均场的应变率相同， 即涡粘系数是各向同性的标量，这一特点是Ｂｏｕｓｓｉｎｅｓｑ假设模式的主要缺陷之一。 与Ｂｏｕｓｓｉｎｅｓｑ假设模式不同，在二阶矩封闭模式中，雷诺应力与平均速度之间的 关系隐含在一组联立的偏微分方程中，在一定程度上反映了湍流的非线性作用。然而 通常形式的二阶矩封闭模式对雷诺应力来说基本上是线性的，所包含的非线性作用非 常小，还不能令人满意地用来解决复杂工程问题，可以说非线性湍流模式在理论和应
用上目前都已成为工程湍流计算研究中的核，Ｉｉ，问题之一【８９】。

４．３湍流运动的基本方程 在工程计算中，往往是流动参数的平均量对工程设计更具有重要作用。因此，流 体的运动方程都建立在雷诺平均的概念下，即对Ｎ．Ｓ方程取某种形式的平均，如时间
平均，空间平均、系统综合平均（指对相同初始条件和边界条件下大量实验结果求平
均）。

对不可压流体来说，平均场的动量方程（即雷诺方程ＲＡＮＳ＝Ｒｅｙｎｏｌｄｓ
Ｎａｖｉｅｒ－Ｓｔｏｋｅｓ

Ａｖｅｒａｇｅｄ

Ｅｑｕａｔｉｏｎｓ）的形式如下：
【ｑ．１ Ｊ （ｔ．１）
？。

鲁五筹一罴＋毒（，誊ｄｘ一瓦）坛 ａ。教ｋ础：融｜、．。… —二＋“Ｉ一一十一ｌｙ一一“．”｝＋Ｐ．

其中：“∽为雷诺应力，是方程中的未知量，Ｕ．和ＪＰ分别为平均速度分量和压力， 岛为体积力，，，为流体的动力粘度。当平均方式为时间平均时，－抛＇－；７。：０。如果流动 是非定常的，则不能取是平均。
不可压缩流体运动的连续方程为：

堕：０
敏．

（４．２）

由于在方程中出现了雷诺应力，因而方程（４．Ｉ）、（４．２）是不封闭的，即未知量

河海大学博士论文

的个数大于方程的个数，因此求解控制方程（４．１）、（４．２）的首要问题是确定雷诺应

力（一“，“，），在湍流模式理论中是通过引进低阶的关系或平均量来近似表示这些湍动 量的，由这些微分形式或代数关系表示的湍动量模拟关系式与控制方程（４．１）、（４．２）
共同组成封闭方程组。

需要注意的是，雷诺应力与粘性应力是不同的，粘性应力是流体分子运动与分子 相互作用的效应，而雷诺应力则是流体宏观的湍动效应，不是流体的物质本性，是流

体运动到一定阶段后的产物，从力的观点看，雷诺应力是迁移惯性力（非线性项）因 取时均而出现的，因此在实质上是属于惯性力的范畴，只有在统计意义下雷诺应力才
会出现。 在雷诺平均框架下的流动控制方程中，湍流的全部信息可以被认为是包含在二阶 速度相关项“。“，之中，若“，“，可以被忽略不计时，则方程（１）回归到层流的控制方 程。由此可见，层流可被看作是湍流的一个极限特例，在大部分工程湍流问题中，甜．“， 不仅不能忽视，而且它对流动混合，扩散等起决定性作用。 在建立“弘的封闭模式之前，有必要讨论雷诺应力输运的控制方程，该方程的具 体推导过程可参见文献［９０１，这里仅给出表达式：

警：导（，学一而一丛小丝西）
一—－＿ 一—－＿

＿—＿．

＿２ｖ当当 面 生％ 丽 堕奶 卜 ‰ｐ 赢一％ 巫魄
＋ ＋ ０冀．Ｏ％．

（４．３）

该方程等式左边为应力的对流输运项（ｑ），右边的第一个括号中的项（吒）代表
湍流的扩散特性，控制湍流的空间分布，由三个不同的物理机理组成：粘性扩散、三 阶速度相关和压力与速度的相互作用。第二个括号表示湍流的生成机制（只），即湍
流在通常情况下由雷诺应力和速度梯度的相互作用所产生。第三个括号为湍动压力与

应变率的相关项，最后一项为应力的耗散率（毛）。为了方便起见，（４．３）可写作

ｑ＝吒＋弓＋ｍ，一毛
通过对方程中的毛，中。和６１ｉ作适当模化即可得到不同的湍流模式。
４．４

（４?４）

ｔ一￡两方程模型。
对方程（４．４）中的ｆ√进行缩并，根据缩并原理及连续条件有中。＝０，可以得到

湍动能≈＝三瓦的控制方程：

河海人学博士论文

嬖：§ｙ婺一一ｋ＇ｕｔ一土丽卜瓦婺一ｙ挈晏 面２面咿面一 一石掣ｐ“１，嘶葛吖菇瓦
项为湍动扩散，右端第二项为湍动能发生项，最后一项为湍动能耗散项。

（４．５）

其中：左端为动能变化率，右端第一项为扩散项，【】中内第一项为分子扩散，后两

在ｋ一￡模型中，为了封闭湍流运动方程，还需要增加耗散率占的控制方程：

坐ｄｔ＝旦＆ｌ悟Ｌ罟磊一万］＿
２ｙ堕（堕ｔ堕兰）一２掣笪型一 ｚＨ等矗一
‘融｜、敏ｊ ａ］ｃ？ ａｘｌ瓠ｔ。 。ａ）ｃｔ瓠ｉ瓠ｌ
口ｆ

：隔２
（４．６）

式中宰为耗散能变化率，右端第一项为扩散（包括紊流扩散和分子扩散）项，第二
及第四项为产生能项，第三及第五项为破毁能项

方程（４．１）一（４．６）共１２个方程，若取Ｂｉ，Ｍ：“：，ｋ，ｓ共１２个量作为未知量?
则构成一个近似封闭的方程组，为了能够求解，还需要将方程组中其余的未知量与 １２个选定的量建立联系，因此还需要对方程（４．４）、（４．５）、（４．６）进行摸化，给出

季：水譬］掣＋，期扣¨。㈢阿却
‘（巴专妒），
ｃｓ个槲
（１个方程）
（４．７）

笔；＝毒［ｃ。（譬］善＋ｙ考］一占

一“

可

面ｉ

警＝言［ｃ。晤）毒＋ｙ爿一ｃｉ（妻］雨≤‘：（守ｃ，个方程，

韩只＝＿（而荨＋丽习～再等。

河海大学博士论文

Ｃｋ＝０．２２５，Ｃｌ＝２．２，Ｃ２＝ｏ．４，ｅ＝ｏ．１３，巳＝１．４５，ｅ，＝１．９２。
这些常数都是由实验或特例的计算给定。 方程（４．１）、（４．２）、（４．７）构成了可以求解的封闭方程组，需要指出的是，方

程（４．７）中出现的六个常数Ｇ，ｃｌ，ｃ：，ｃ。，ｔ．，ｃ。是在模化过程中采用近似的结果， 显然在模化中采用不同的假定会得到不同的常数。
在标准ｋ一占模型中，雷诺应力的计算并不采用它的摸化方程，而是直接采用 Ｂｏｕｓｓｉｎｅｓｑ假设，即雷诺应力与平均场应变率应有如下关系：

一瓦＝ｕ譬＋等）一；毛七
其中：涡粘系数＿在ｋ—ｓ模式中取

（４．ｓ）

ｖ，＝Ｃ。生

而湍流动能（膏＝寺蚝“，）和湍能的耗散率ｓ采用模化方程：

面Ｄｋ＝毒（（¨尝）毒）＋ｐｔ—ｓ 尝＝言ｃｃＶ＋≯毒，＋妻ｃｃ以一ｔ：ｓ，
其中仇为由于平均速度梯度引起的湍动能生成项，既：一ｐ丽挈

（４．９）

ｃａ．?。，

方程（４．９）（４．１０）与方程（４．１）（４．２）即构成了标准ｋ—ｓ湍流模型，在该模型中 的常数取值如下：

Ｉ
｜

ｃ， 。．。９

Ｃ引 １．４４

Ｃ；２ ｌ，９２

盯★

盯‘

ｌ，Ｏ

１．３

ｋ一￡两方程模型是目前应用比较广泛的一种湍流模型。在推演过程中，采 用了以下几项基本处理：（１）用湍能七反映特征速度；（２）用湍能耗散率占反映特征 长度尺度；（３）引进了■＝Ｃ。ｋ２／Ｆ的关系；（４）利用了Ｂｏｕｓｓｉｎｅｓｑ假定进行简化。 因此，ｋ一占模型具有以下优点：（１）通过求偏微分方程考虑湍流物理量的输运过程， 即通过求解偏微分方程确定湍动特征速度与平均场速度梯度的关系，而不是直接将两 者联系起来：（２）特征长度不是由经验确定，而是以耗散尺度作为特征长度，并由求 解相应的偏微分方程得到。由于脉动特征速度和特征长度是通过解相应的微分方程求 得，因而ｋ一占模型在一定程度上考虑了流动场中各点的湍能传递和流动的历史作用。
４７

河海大学博士论文

计算结果表明，它能比较好地用于某些复杂的流动，例如环流、渠道流、边壁射流和 自由湍流，甚至某些复杂的三维流等。但是它也有局限性，主要有以下几点：（１）采 用Ｂｏｕｓｓｉｎｅｓｑ假定，即采用了梯度型和ｕ各向同性的概念，因而使ｋ一占模型难以准

确地模拟剪切层中平均流动方向的改变对湍流场的影响；（２）引入了一系列的经验系 数，而这些系数都是在一定试验条件下得来的，因而也限制了模型的使用范围。 ４．５ＲＮＧ女一占模型 ＲＮＧｋ—ｓ模型由Ｙａｋｈｏｔ和Ｏｒｓｚａｇ于１９８６年应用重整化群（Ｒｅｎｏｒｍａｌｉｚａｔｉｏｎ
Ｇｒｏｕｐ）的方法导出，对于高雷诺数湍流，ＲＮＧｋ一占模型具有和标准ｋ一￡模型相同

的形式，只是在ｓ方程中出现了一个附加生成项，当流动快速畸变时，这一项将显著 增大。ＲＮＧｋ一占模型与标准ｋ一占的主要不同之处在于它们的模型系数不同，标准 ｋ一￡模型系数由经验给定，而ＲＮＧ＾一￡模型系数均由理论精确计算得出。此外
Ｙａｋｈｏｔ和Ｏｒｓｚａｇ认为ＲＮＧ ｋ一占模型不需要采用壁面率，可直接积分到壁面，而文

献【４６】作者在应用中发现，若采用ＲＮＧｋ—Ｅ模型直接积分到壁面，必须在壁面附近 划分足够细的网格，这必然增加计算工作量，尤其在三维工程湍流计算中巨大的网格 数量往往是难以接受的，因此类似于文献［４６］的作法，本论文采用考虑壁面率的
ＲＮＧ女一占模型，并将ＲＮＧｋ—ｓ湍流模型推广应用于三维湍流模拟。

ＲＮＧｋ—ｓ模型与标准的ｋ—ｓ模型具有相似的形式，该模型中的两个标量方程
ｋ方程与ｓ方程形式如下：

ｉＯｋ Ｍ婺：鲁（（ｙ＋互）婺）＋Ａ一占 ｏｌ
出Ｉ 盘Ｉ ｏｋ呶‘

（４．１１）

?０动６＇＋ｕ，毒＝云（（，＋与Ｏ＂。罢Ｏ＇Ｘｉ）一矗＋｛（ｅ．ｎ—ｅ：占）
程中的Ｒ项为ＲＮＧｋ一￡模型与标准ｋ～￡的主要区别，Ｒ＝

（４．１２）

其中：＾方程中Ｐ．为由于平均速度梯度引起的湍动能生成项，Ｐｋ：一ｐ丽挈，占方
。Ｏ＇Ｘ．

—百矿Ｔ’
巴矿（１一ｑｌｑ。）￡２

平均应变鄹㈣珊７＝厝孝船均湖搁瞄湍流时间尺趔Ｅ，

河海大学博士论文

湍流粘性系数∥。＝ｐ?Ｇ。七二／５，叩。＝４．３８是印在剪切流中的典型值， 口＝０．０１５。此外．与标准女一ｓ模型的不同之处还有方程中的常数不是用经验方法 确定的，而是采用ＲＮＧ理论推导出的精确值，具体数值如下：
２

ＣＦ：Ｏ．８４５，口ｋ＝０．７１７ ９，Ｏ－。＝０?７１７ ９，Ｃｃｌ－１．４２，ｃ。２

１?６８

与标准女一ｓ模型相比，ＲＮＧｋ一￡模型具有以下特点：

①ＲＮ６ｋ一￡模型的占方程中有一附加项Ｒ，该项代表平均应变率对￡的影响，可以 改善迅交流中的模拟精度： ②在ＲＮＧ≈一ｓ模型中包含了旋涡对湍流的影响，提高了对旋涡流的模拟精度； ③方程中的常数不是用经验方法确定的，而是利用ＲＮＧ理论推导出的精确值； ④ｋ一￡模型为高雷诺数模型，而ＲＮＧｋ一￡模型则同时适应于低雷诺数的情形，
此时需要在近壁面处作特别处理

综上所述，ＲＮＧ七一ｆ模型比标准七一占模型具有更高的精度和更广的应用领域。
４．６

ｋ—ｓ湍流模型中的边界条件
对ｔ一￡方程求解需要给出合理的边界条件，这些边界条件的应用对所有的方程

基本上是相类似的，只是对于固体壁面。问题可能变得比较困难。在壁面上取七＝０，
ｊ６

但￡≠０，由娑＝０替代，其中”为壁面的法向。
咖

对于高见流动问题，由于粘性底层很薄，无法在壁面处采用非常细密的网格模 拟，解决的方法是采用壁面函数（ｗａｌＩ Ｆｕｎｃｔｉｏｎｓ），该函数与速度分布的对数有关。湍
流边界层的速度分布在图４．１中给出，从图中可以看出，在对数区域的速度分布为：
，，

１

“＋：二土：二ｌｎｎ＋＋Ｂ
“．

（４。１３）

Ｆ

∞矿俗
ｍ

，

ｏ ２

５

１０

２０

５０

１００矿

图４．１湍流边界层速度分布

河海大学博士论文

斯∞与壁面平行嗍溉“。黼应力戤或称为摩擦斌铲层，
ｆ。为壁面切应力，Ｋ称为卡门常数（Ｋａｍａｎ ｃｏｎｓｔａｎｔ）（ｒ：ＯＡｌ），Ｂ表示粘性底层厚
度的实验常数（平板边界层中Ｂ 化。
ｚ

５．２），”为网格点到壁面的距离，ｎ＋为ｎ的无因次

ｎ＋：型
“

（４．１４）

由于通常假定流体是局部平衡的，即湍动能的产生与耗散是近似相等的，因此，

“，＝ｃ：４√．ｉ｝

（４．１５）

由（４。１４）（４．１５）两式可以推导出壁面上第一个网格点的速度与壁面切应力的关系

铲倒；＝斛４ｒ压赤
建立更加精确的方程。

（４．１６）

其中：Ｅ＝Ｐｏ，由于紧靠壁面的控制体的一个控制面在壁面上，因此在动量方程中

对于平行于壁面的控制体，其壁面上的切应力值是需要的，因此利用该边界条件可以

当采用壁面率边界条件时，女在壁面上的扩散通量取为ｏ，即婺：ｏ，关于耗
散率占的边界条件可以根据局部平衡的假定，即壁面处产生与耗散近似相等推导出 来。在壁面区，湍动能的产生项可由下式计算：

只一。鲁
ａＨ

（４．１７）

同时该区内切应力近似为常数，为了计算出控制体中心处的耗散率，通过对（１）式
的对数速度分布求导，可得：

舀Ｏ，：÷：掣
、

ｎ）８ｉ他～＇’ｊｏ Ｊ ’Ｐ槲， 肼ｐ

该式与（４．１６）式联立即可计算出耗散率。
当采用上述边界条件时，近壁面处的控制体中是不需要求解占方程的，而是由下

河海大学博士论文

式直接计算

。Ｐ

，：竺竖
砌口 必须注意，上述边界条件是使用只有当靠近壁面的第一个网格点处于对数区，即

月：＞３０时才有效。 此外，在入口边界处，ｋ，６通常未知，一般给膏赋一个小量值，如１０一云２，而￡的

值一般由式占。譬确定。
４．７算例：水翼流场模拟 利用非正交网格同位网格系统中的ＳＩＭＰＬＥ算法对ＮＡＣＡ００１５型二维水翼粘流 场进行数值计算，湍流模型分别采用标准ｋ一￡模型和ＲＮＧｋ—ｅ，雷诺数Ｒｅ＝１．９７ｅ～， 在图４．２中给出了不同攻角下升力系数曲线，并与文献［９１］的实验结果进行比较，在 小攻角下两种湍流模型的计算结果与实验值均比较一致，ＲＮＧｋ－￡模型略优于标准
ｋ一￡模型。
１．４ １．２ 甚１ ０．８ ０．６
０．４

０．２ ０
Ｏ

２

４

６

８

１，１４ ｇ（ｄｅａ）

１（１

１６

１８

２０

２２

图４．２
４．８小结

ＮＡＣＡ００１５水翼升力系数一攻角变化曲线

本章讨论了湍流模式理论研究的工程背景和意义，详细给出了标准ｋ一占模型和 ＲａＮＧｋ一￡模型，比较了两种模型的特点及应用范围，为更好地处理固壁边界，采用 与壁面率相结合的ＲＮＧｋ一占模型模拟湍流流动，通过对ＮＡＣＡ００１５型水翼湍流场的 模拟表明，ＲＮＧｋ—ｓ模型能够更好地模拟出流场的湍流特性。

河海大学博士论文

第五章运动界面数值模拟技术
５．１运动边界模拟技术简介 自由面模拟技术在水利工程、船舶工程、海洋工程、机械工程等领域有着深刻的 现实意义和应用价值，随着Ｈｉｒｔ等人提出用ＶＯＦ方法模拟自由面追踪技术，目前已 经发展和形成了许多实用软件，如ＶＯＦ．３Ｄ，ＮＡＳＡ．ＶＯＦ，ＳＯＬＡ．ＶＯＦ等，同时推广 到一般介质面的运动界面模拟，该理论已经成功应用于半导体工业、喷墨技术等现代
高新技术产业。据不完全统计，１９９９年度的Ｊ．Ｃｏｍｐｕｔ．Ｐｈｙｓ．杂志涉及运动界面追踪问

题的文章不少于３０篇，可见运动边界数值模拟的研究和分析方法受到越来越多的关 注，但是由于运动界面的形状和整体位置不断变化，所以精确模拟的难度非常大。 运动界面的模拟模型一般分为两大类：欧拉模型和拉格朗日模型。欧拉模型的特 征是坐标系或者固定或者按给定的方式移动以适应不断变化的计算域形状，其主要优 点是可以适应内部界面的大变形和大幅度移动而能保持一定的计算精度。利用欧拉模 型需要解决的主要问题是用什么形式的控制方程最方便，如何选取合适的数值离散方
法和最佳的内界面跟踪技术。拉格朗日模型的特征是坐标随流体雨运动，但每个网格 单元内包含的流体质点保持不变。该模型的优点是可以精确跟踪内界面并准确描绘出 运动界面的形状，便于使用界面的运动学和动力学边界条件；拉格朗日模型还可以在

运动界面附近进行网格加密，以精确模拟运动界面附近的流场内部结构，同时精细网 格区可随界面移动，大大提高了运动界面的模拟精度。拉格朗日模型的主要不足是随 着计算的进行，网格会出现重叠、缠绕等现象，随着网格越来越不规则，计算精度也
越来越低。 在欧拉模型中，模拟运动边界的典型方法有三种：即ＭＡＣ法，ＶＯＦ法，Ｌｅｖｅｌ
Ｓｅｔ

方法。ＭＡＣ法的基本原理是在计算区域中包含流体的网格内布设“标记点”（Ｍａｌ（ｅｒ），

标记点没有质量，不参与流场计算，只是用来标记流体质点的运动轨迹，随着流场计 算的进行，可获得每个标记点的速度，并由速度确定出标记点的最新位置，根据标记 点的最新位置即可确定哪些网格单元被流体充满，哪些网格单元被流空，从而确定运 动界面或自由液面在下一时刻的位置；ＶＯＦ法是Ｈｉｒｔ和Ｎｉｃｈｏｌｓ［９２１于１９８１年首先提

出，通过对主场的计算采用了Ｄｏｎｏｒ－－Ａｃｃｅｐｔｏｒ（施主一受主）的概念、Ｕｐｗｉｎｄ特征
思想和补偿效应，成功地模拟出溃坝和涌浪（ｂｒｏｋｅｎｄａｍ，ｂｒｅａｋｉｎｇｂｏｒｅ）自由面，避 开了对于自由面问题通常采用的高度函数法和计算量十分浩大的Ｍａｒｋｅｒ点的方案，
引入了一个流体体积函数方程，与流体运动基本方程联合求解，开辟了自由面计算的

新途径。尽管当初的格式、方法是比较粗糙的，但其基本思想至今一直受到重视，特 别是对运动界面的重构方法开创了一个新的研究方向：１９８８年，Ｏｓｈｅｒ［９３１等人提出了

河海大学博士论文

一种零等值面方法（Ｌｅｖｅｌ Ｓｅｔ方法），在许多复杂的界面追踪问题中获得成功，该方

法的主要思想是将运动界面定义为一个函数的零等值面（线），即妒（ｘ，，）＝０，然后始 终保持它是零等值面，并且在界面附近保持单调。
５．２

ＶＯＦ方法 ＶＯＦ方法的基本思想是定义一个函数Ｆ，在有流体的点取值为１，在无流体的点

取值为０，这样在一个单元内，Ｆ的平均值就可以代表单元内流体所占的份额，如果 Ｆ的均值为１，则单元内充满流体，如果单元内无流体，则Ｆ的值为零。而Ｆ值在０

和１之间则表示该单元必然包含自由液面，所以，ＶＯＦ法提供了与ＭＡＣ法类似的有 关内部界面的信息。由于ＶＯＦ方法只需对每一个单元提供一个内存就可以了，不像 ＭＡＣ法那样需要对标志点提供内存，从而节省的计算空间，提高了计算效率。 ５．２．１控制方程
流体体积分数Ｆ随时间的变化过程由下面的方程控制：

～ＯＦ＋“笪＋ｖ堡：０
ａｆ

缸

（５．１）

咖

对于不可压缩流体，控制方程的守恒形式为：

堡＋Ｏ（Ｆ／２）＋型：０
ａ

孤

（５．２）

砖

自由面模拟问题的控制方程一般由上面的流体体积函数方程与流体运动基本方
程构成。

５．２．２流体体积分数Ｆ的输运模型 以Ｈｉｒｔ和Ｎｉｃｈｏｌｓ建立的施主一受主模型（ｄｏｎｏｒ—ａｃｃｅｐｔｏｒ）为例说明Ｆ输运的基本思
想。

考虑网格的右侧边界，国时段内，通过该网格面的流体和空气的体积通量 Ｖ＝“田，其中”为网格面上的法向速度，其符号决定了网格是“施主网格”还是“受 主网格”，若“的符号为正，则网格面上游是施主网格，下游是受主网格。在一个时 间步长内通过网格面的Ｆ的通量等于６Ｆ与网格面积的乘积。８Ｆ由下式计算：

ｂ＇Ｆ＝ｍｉｎ眈ＤＭ＋ＣＦ，Ｒ衄Ｄ｝

（５．３）

其中第一项表示输运的流体的量，第二项表示施主网格所拥有的流体的量。

ＣＦ＝ｍａｘ｛（１．０一只Ｄ）ｌｙＩ一（１．ｏ—ＦＤ）缸Ｄ，０．０｝

（５．４）

河海大学博士论文

式中，下标Ａ代表受主网格，下标Ｄ代表施主网格，下标ＡＤ则表示可以是施

主网格，也可以是受主网格，根据边界面上流动的方向确定，比如，同一个网格相对
于不同网格边界的情况。（５．３）式中的Ｍｉｌｌ是为了防止通过界面的通量超过施主网格

所能提供的Ｆ的量，（５．４）式中的ｍａｘ是为了保证当空气的通量超过可能提供的量时，
需要补充流体通量，这个补充的量就是ＣＦ。

在图５．１中给出了施主一受主模型通量输运的基本形式，当流动通过一个垂直的单 元边界面时，在（ａ）中给出了施主网格与受主网格的定义：当下游网格为受主网格

即ＡＤ＝Ｄ时，则凹＝０，Ｆ的通量为Ｆ：Ｆｏｌｖｌ，其中Ｖ＝小６ｔ，Ｆ值用于确定网格
边界面通过流体的那部分断面面积，当然要求川＜Ａｘ以保证施主网格内的流体不会
流空，见（ｂ）：当ＡＤ＝Ａ时，则施主网格内的Ｆ值用来确定网格边界面上通过流体的
断面面积（ｃ），此时施主网格内的所有流体都将流失，即所有存储在虚线和边界面之 间的流体和空气均进入受主网格；在（ｄ１中表明，当ＡＤ＝Ａ，但通过网格面的流体通

量超过ＹｏＩＶＩ时，则需要补充额外的流体以满足需要的通量，其中虚线和边界面之间
的额外流体就等于ＣＦ。

（ｂ）

（ｃ） 图５．１
（ａ）

（ｄ）

施主．受主模型中Ｆ通量的计算方法举例

施主?受主模型网格的配置，（ｂ）、（ｃ）、（ｄ）中双阴影线区域表示实际传输的Ｆ的量

事实上，使用施主网格中的Ｆ值还是受主网格中的Ｆ值来确定边界面上通过流体 的断面面积完全取决于自由液面的平均定向，这可分为两种情况①当自由液西基本上

河海大学博士论文

沿其法向移动时，ＡＤ＝Ａ，②当自由液面基本上沿其切向移动时．ＡＤ＝Ｄ。此外如
果受主网格是空的，或者施主网格的上游是空的。则一律有ＡＤ＝Ａ，也就是说，施 主网格必须首先补充流体，然后才能向下游的空单元输运流体。 用通量值５Ｆ乘以流体通过边界的断面面积就可以得到输运流体的量，从施主网

格中减去该通量值添加到受主网格，然后对所有的网格边界重复上述步骤即可完成Ｆ
的输运计算。 在上述计算过程中可能会出现Ｆ值小于零和大于ｌ的情况，因此必须对Ｆ的 计算值归整到０和１之间，同时将由此引起的误差进行累计，以保证总的流体体积守

恒。此外由于在数值计算中产生的舍入误差使得Ｆ的值不可能正好是１或０，为此通 常设定一个很小的数如占，＝１０一，当Ｆ小于卸时，网格是空的，当Ｆ大于１一￡时，
网格是满的。

５。２。３运动界面的重构方法 Ｈｉｒｔ和Ｎｉｃｈｏｌｓ的贡献主要体现在：～是作者采用的施主一受主（Ｄｏｎｏｒ－Ａｃｃｅｐｔｏｒ） 方法给出了界面通量输运计算的初始格式，同时该方法不仅考虑本网格单元上的流体
体积分数，而且用到相邻网格上的流体体积分数，从而提高了计算精度：作者的另一 个重大贡献是提出了一种自由面的重构方法，因此掌握运动界面的重构方法必须从

Ｈ—Ｎ模型出发。

一赫
．÷一１Ｉ历黝嬲勉
图５．２施主．受主模型中的自由面重构方法

该模型重构运动界面的基本思想是：将流体自由面看作是局部的单值函数 ｒ（ｘ）或Ｘ（ｙ），计算中需要用到９个网格的信息，如图５．２。首先计算ｆ一１，ｆ，ｆ＋ｌ网格

列的ｒ值和Ｊ一１，，，Ｊ＋１网格列的ｚ，值，估算出每个网格上自由面的斜率值譬，譬， 靠聊
然后根据流体体积分数和斜率的大小确定网格（ｆ，Ｊ）上的自由面位置和方向，即

河海大学博士论文

Ｋ＝∑Ｆ，ｋＳｙＩ，

？＿ｉ一１，ｉ，ｉ＋１

（５ ５）

（氅，：ｏ暨斗
、ｄｒ“

（５ ６） 、…７

苏…＋２融，＋苏Ｈ

ｘ，＝∑ｒ，ｋ４ｖ＾，

，＝产１，Ｊ，ｊ＋ｌ

（５ ７）

（筝，：善婆兰业 Ｉ一｝
、咖。 ￥Ｈ＋２印ｆ＋咖『－ｌ

理论上讲，以上两个导数应该互为倒数，但数值计算结果则不然，作者的意图是比较
它们绝对值的大小，来决定自由面在该网格内是“水平”的还是“垂直”的。即如果

＝一

（５。，

ｌ：，．６， 、 。

引爿
则定义自由面是水平的；否则自由面是垂直的，如图５．２所示。

（５．９）

尽管Ｈ．Ｎ运动界面重构方法比较简单，但它却起到了抛砖引玉的开创性作用， 在这一思想的启迪下，自由面的重构技术得到了迅速的发展，比如：Ｃｈｏｆｉｎｌ９４】引入

ｌ，４和３／４单元以连接水平和垂直的自由面单元；Ｐｕｃｋｅｔｔ［９５１用离散的线段来表示自由
面，通过相邻网格的体积分数预测自由面的法向；Ｈｏｎｇｆａｎｇｗｅｎｌ９６ＪＮＮ斜线段来代替 自由面，其中斜线段方程中的系数由自由面网格周围８个网格的体积分数确定：

Ａｓｈｇｒｉｚ和Ｐｏｏ的ＦＬＡＩＲ方法ｐ。】则利用每两个水平相邻网格的体积分数求出斜线段来 模拟自由液面，但只适用于均匀网格情况。Ｒｕｄｍａｒｌｌ９”对上述两类自由面重构算法进
行了比较计算指出，利用斜线段来逼近自由液面比其它方法精度更高，这些方法的基

本思想都是首先确定网格内运动界面的具体位置及其运动方向，然后确定运动界面在 下一个时刻的新的位置，同时更加注意运动界面与网格线的交点的位置，因此得到的 运动界面也更加精确，普遍采用的格式是用网格内具有法向量的斜直线段Ｙ＝ａｘ＋ｂ 来逼近，也有用二次曲线Ｙ＝Ｏ．Ｘ２＋ｂｘ＋ｃ来构造运动界面１９８１，只是考虑得更加精细而
已。

本论文利用ＦＬＡＩＲ方法重构自由面，并对该方法进行了推广。
ＦＬＡＩＲ自由面重构方法的基本思想是在两个相邻的网格内构造一条倾斜的直线

段来近似该网格边界的界面，然后计算透过该网格边界的流体逶量，作为修正流体体
积分数的数值通量，由于界面构造仅涉及到两个网格，故需要根据网格内的流体体积

分数划分为许多种情况分别计算，各个网格内的流体体积分数可能有以下四种情况，
见图５．３。

河海大学博士论文

（Ｉ）施主网格是满网格，即瓦＝１ （２）施主网格是空网格，即屹＝０ （３）施主网格是半网格，受主网格也是半网格，即０＜ＦＤ＜１，０＜只＜１ （４）施主网格是半网格，受主网格满网格或半网格，即０＜％＜１，只＝０，只＝０

（ｃ）

ｉ－，？Ｔ Ｔ二

二

ｊｊｉ彩
（ｅ） （ｆ）

∥％／

图５．３相邻网格间的关系

情况（１）、（２）比较简单，情况（３）可以分为１６种，经过转换可以归结为４种，

见图５．４，即￡≥厶，且目标流体在运动界面以下。此时只要确定了直线的斜率及位
置即可计算出通过网格边界的流体输运量，如对情况（ａ），设左边是施主网格，则舀 时间内输运到相邻网格的体积量Ｖ＝”．西，即图中的斜线阴影区域的面积；如果右边 是施主网格，则应取右边网格内靠近网格边界对应区域面积，施主网格和受主网格由 网格边界上的速度方向确定。

河海大学博士论文

笏

］赫勿
（ｃ， （ｄ）

图５．４

ＦＬＡＩＲ模型运动界面重构方法

情况（４）最为复杂，以图５．５（ａ）为例，施主网格（ｆ√）是半网格，而受主网格（ｉ＋ｌ，，）

是空网格，可以用网格（ｆ，『）垂直方向的两个相邻网格分别按照第三种情况构造界面，
然后求通过网格边界流入受主网格的流体体积输运量，同样可以归结为图５．５（ｂ）所示 的四种情况：若受主网格是空网格，则取靠右边的宽度为“．西的区域面积作为输运

量，如果受主网格是满网格，则取左边的宽度为Ｉｔ．毋的区域面积作为输运量，即保
证流体分布靠近满网格一边。

Ｅ瓢 篓貔 遴～丝擞
（ａ）

（ｂ）

图５．５运动界面的四种重构方法

下面以情况（３）为例给出ＦＬＡＩＲ模型构造运动界面的方法，见图５．４，并将该方法

推广到非均匀网格，丘，兀分别表示左、右网格内的体积分数，见图５．６。

河海大学博士论文

Ｃａｓｅ １

Ｃａｓｅ ３

Ｃａｓｅ ４

图５．６
Ｃａｓｅｌ

运动界面结构形式

不妨假设两个相邻网格中的体积分数满足：Ｌ≥以，自由面认为是与网格边界相交 的直线段，故自由面可由Ｙ＝删＋６来描述，其中系数口，ｂ由已知的体积分数Ｚ，五确
定。因此可得到下面的方程组：

厶?血，缈＝袅，（ａｘ＋６协＝ｉ１口（ｘ｝ｘ２Ｉ）＋ｂ（ｘｉ－ｘ，＿Ｉ） 兀－Ａｘｉ＋ｌ缈＝璧＋１（ａｘ＋ｂ）ｄｘ＝去ｎ（毒ｌ一＃）＋６（鼍＋Ｊ—ｘｉ）
联立求解，得到

。：！亟二型心
ｂ＝厶．缈＋立二生Ｌ（厶一厂６）．缈
＾ｉ＋Ｉ—ｊ卜Ｉ
ｃａｓｅ

２

类似于ｃａｓｅｌ可列出如下方程

厶－Ａｘ。Ａｙ＝去Ⅱ（ｚ｝ｘ三１）＋６（ｘｆ—Ｘｉ－Ｉ）
／ｊ?缸，＋ｌ?Ａｙ＝：Ｉ口（ｘ；一ｘ；）＋６（ｘ６一Ｘｉ）
由于出现了新的未知量％，ｔｔ需奉ｔ．充一个方程

河海大学博士论文

０＝ａＸｂ＋ｂ

联立上述三个方程，求解可得

口：—２（ｆｎ．ａ—ｙ－ｂ）
Ｚｔ＋Ｘｉ一１

其中，。：Ｘｉ，，：垒』型
ｋｘｉ

ｂ＝２ｌ帆＋如）一

缈

缸？

ｃａｓｅ３

类似于ｃａｓｅ ２可列出方程组
Ｙｊ。ＩＴ－Ｘａ＋ｂ

厶缸：缈＝ｘ。一ｘＨ）Ａｙ＋罢ｂ； ｘ。２）＋６０。一工。）

厶止Ⅲ每＝昙ｃｃ磊，一ｘ；）＋ｂ（ｘ“
联立求解可得：

口：—２［Ａａ—ｙ－ｂ］
ｂ＝，＿【（２ｓ＋２ｔ一１）一２吮一２（ｓ—ｐ）Ａ

—ｚ石ｉ面刁‘面习Ｆ丽ｉ西可了了】．缈
其中５：二Ｌ

血．

”—■叠一
。一ｋｘｆ＋ｌ
ＬＵ．

ｘ，＋ｚ』＿１）

，：塾！！二兰二塾＝！
Ｘｆ＋Ｘｉ一１

Ｐ：

！型二！！＝１
２ｘＪ＋Ｉ—ｘｆ一３ｘｆ—

ｃａｓｅ４

Ｙｊ＝ａｘ日＋ｂ
Ｙ，一１＝ｃ“６＋ｂ

ｘｎ＋ｘ６＝２（ｆａ?血ｉ＋兀ｋｘｆ＋１）＋２ｘｊ

河海大学博士论文

＾缈．缸Ｍ＝昙僻一ｘ；）＋６（Ｘｂ－－Ｘｉ）
联立求解得：

口：坐
ｘｄ一工６

ｅ～川一ｉ△ｘｉ．ｘｂ
｛。ｂ ２——五『＿一
｝ Ｂ＋√Ｂ２—４ＡＣ

ｌＸａ＝２（／：－Ａｘｆ＋．，ｊ?血Ｊ＋１）＋２ｘｆ一１一ｘ６

其中：爿＝２ｙｊ－Ｉ＋等
Ｂ＝２（＿‘?△ｘｆ十ｆｂ?△ｘ“ｌ＋ｘｌ“＋ｘｔ）Ｙ』一１＋Ｇｆ＋２ｆｂ△ｘｆ＋１）?却

ｃ＝２眈?缸，＋＾－缸ｉ＋１＋Ｘｉ＋１）（ＸｉＹｊ－Ｉ＋厶血州?缈）＋等缈
容易验证，对于尺度为ｈ的方形网格，文献［６】中的结果即为本方法的特例。 显然在确定了ｎ，ｂ之后，利用Ｙ＝似＋ｂ即可完成自由液面的重构。

本文研究的运动界面是自由液面，自由液面单元的压强Ｐ。需要利用内插（或外

ｐ“２（１一＿）ｐ～＋秘ｓ
这里，７７＝ｄ。／ｄ，ｄ。是自由液面单元和邻近
自

的内部单元中心间的距离，ｄ是自由液面至内部
邻接单元中心的距离，用于内插的邻接单元的选 取标准是，联接自由液面单元中心和所选取内部 单元中心的联线最靠近自由液面的法线，该单元 又称内插邻接单元。

、平
／ｐ¨

—、～

王ｉ尘
ｐ＿

图５ ７自由面上的压强边界条件

河海大学博士论文

５．２．５运动界面的宽度与等值线的给定方法

由于运动界面的初始位置是由流体体积分数定义的㈣，随着流场的变化和时间的
推移，计算所得到的自由面轮廓是从流体体积分数为０的高度变到１的高度，即计算 流场中得到的自由面为Ｆ∈［Ｏ，１１的模糊界面，其宽度可以跨越几个阿格。如图５．８

图５．８运动界面的宽度

所谓等值线，就是根据不同问题的实际情况，流体体积分数取［ｏ，ｌ】区间内的某个 值，如令Ｆ＝Ｏ ５的等值线作为自由面的位置，当然该值不是固定的，有可能在其它 情况下取Ｆ＝０．６更合适。这是采用单一定值等值线确定自由面的方法。此外，还可 以利用Ｆ值的百分比宽度来给出边界面的形状，如９６％宽度的界面；显然，如果有 某种方法可以进一步将具有一定宽度的界面重构为单值轮廓线，则得到的自由西位置 将具有更高的精度。 ５．２．６数值稳定性条件 （１）对流稳定条件 为了正确的模拟自由液面或内界面的变化，我们只考虑相邻单元之间的通量，因 此，时间步长必须满足下面的不等式，

趴幽协尚｝
上式要求在计算网格确定之后，时Ｉ＇ｕｑ步长西的选取必须使流体在该时段内流过的 距离不能超也网格的尺寸，由于ｒａｉｎ是对所有的单元取的，故在实际计算中的毋一 般取上式右端的ｌ／４～１／３。

河海火学博士论文

（２）扩散稳定性 对于粘性流体，要求在每一个时间步内，动量的扩散同样不超过一个网格单元。
线性稳定分析给出如下的限制条件：

Ｖ藤２篇Ａｘ

２

．２＋△ｙ．

５．３运动边界模拟算例 算例１：Ｚａｌｅｓａｋ模型 为检验本文所建立的的自由面重构算法的可靠性，对经典的Ｚａｌｅｓａｋ界面模型在 平面旋转流场中的运动进行了计算模拟，计算状态为：计算区域【０，１］×【０，ｌ】，网格
节点数２００ｘ ２００．旋转速度场为：ｕ（ｘ，Ｙ）＝—石（ｙ—Ｏ．５）＇ｖ（ｘ，力＝ｘ（ｘ—Ｏ．５），初始位置

为一个下方有缺口的圆周；计算时间步长为母＝Ｏ．０００５。在图５．９中分别给出了经过 四个时间步长ｆ＝０．５ｓ，ｆ＝ｌｓ，ｆ＝１．５ｓ，ｆ＝２ｓ对应的二维界面图形，同时给出了文献［１００１
中采用Ｈｉｒｔ和Ｎｉｃｈｏｌｓ模型的计算结果。可以看出，Ｈｉｒｔ和Ｎｉｃｈｏｌｓ方法给出的界面

比较粗糙，无法准确地描述运动后的界面形状，而本文所给出的方法则能够获得非常
好的结果。

，＝０．５ｓ

ｆ＝１．０ｓ

ｆ＝１．５ｓ

ｔ＝２．０ｓ

（ａ）本方法

ｔ＝０．５ｓ

ｆ＝１．０ｓ

ｔ＝１．５Ｊ

ｆ＝２．０ｓ

（ｂ）Ｈｉｒｔ＆Ｎｉｃｈｏｌｓ方法 图５．９

Ｚａｌｅｓａｋ模型的计算结果

算例２：均匀流场中圆形界面模型
计算区域及网格设置与算例ｌ相同，圆的初始位置为：圆，１１，（０．３，０‘３），半径ｒ＝０．２

河海大学博士论文

均匀速度场为Ｕ＝０．２，ｖ＝０．２，经过四个时间步长圆的形状及位置如图５．１０。相对于 初始圆周，面积的变化量分别为ａｓ＝Ｏ．０００１５４，０．０００１５４，０．０００１５３，０．０００１５８，其中每一 步的面积计算是根据相应的面积分数相加而得到的。

ｆ＝０．５ｓ

ｔ＝１．０５

，＝１．５ｓ

ｆ＝２．Ｏｓ

图５．１０均匀流场中圆的移动

算例３：剪切流场中圆形界面的运动 平移和旋转流不会使运动界面的形状发生变化，而在真实的流体问题中会发生界

面不断变形的复杂物理现象，为了考察本文建立的算法对各种现实流场的模拟效果，
计算了剪切流场中圆形界面的变形运动，剪切速度场为
ｕ（ｘ，ｙ）＝７ｒｃｏｓ（ｘ（ｘ—Ｘｏ））ｓｉｎ（石（ｙ—Ｙｏ）） Ｖ（ｘ，ｙ）＝一，ｒｓｉｎ（ｘ（ｘ—Ｘｏ））ｃｏｓ（ｘ（ｙ一），０））

计算区域为【Ｏ，１］×［Ｏ，ｌ】，初始界面为圆心在（Ｏ．３，Ｏ．３）、半径为０．２的圆周，见 图５．１１（ａ），计算结果给出了ｒ＿Ｏ ５ｓ，忙ｌｓ，ｔ＝１．５５，，＝２ｓ圆形界面的变化情况（图 ５．１２），同时以ｆ＝２ｓ的结果作为计算初值，将速度反号进行反剪切２秒，结果在图５．１１
（ｂ）中给出。

ｌ

Ｉ

￡丫∞，。３一
、、＼ ｊ？ｆ

＼ｊ

（ａ）

（ｂ）

图５．１ｌ

界面初始位置（ａ）及反剪切后的恢复位置（ｂ）

河海大学博士论文

ｒ＝０．５ｓ

ｒ＝１．Ｏｓ

ｆ＝１．５ｓ

ｆ＝２．０ｓ

图５．１２剪切流场中界面的位置和形状

５．４小结

本章讨论了运动界面的模拟技术，重点研究了ＶＯＦ方法中运动界面的高精度重
构算法，以ＦＬＡＩＲ方法为基础，建立了可适应于非均匀网格的自由面重构模型，通 过对旋转流场中的马踢铁形界面运动、平移流场中圆形界面的运动、剪切流场中圆形 界面的运动进行模拟表明，本章建立的运动界面重构模型相对于其它方法能够更加准
确地确定运动界面的形状和位置。

河海大学博士论文

第六章溢流坝过坝水流二维流场数值模拟
６．１溢流坝过坝水流的研究方法
溢流坝的水力设计对于整个水力工程建设具有重要意义，其中溢流坝的流量系 数、自由水面线的形状和位置、坝面压力及坝面流速等都是水力设计中的重要参数， 也是选择合理溢流坝型的重要依据。这些数据多年来主要通过水力学模型试验或由经 验公式而获得，一方面要耗费大量的人力、物力和财力，同时受到测量条件的限制不

可能给出整个流场区域的全部信息。随着计算技术的迅速发展，利用数值模拟方法求
解溢流坝过坝水流场的各种计算模型相继产生，如有限差分方法、边界元方法［６０，１０７】、 有限元方法［１０３，１０４，１０５］、边界拟台曲线坐标变换方法【‘吲等。
有限差分法简易明了，易于理解，其主要不足是不能很好地处理不规则边界，因

此在过坝水流数值模拟中的应用受到了限制。

有限元法比有限差分法能更好地拟合复杂边界，因而在过坝水流的数值模拟中得
到广泛而成功的应用。其中最早将有限元方法应用到过坝水流求解的是日本的池川昌 弘ｌＩ“Ｊ，并对自由面的迭代提出了两种方法。该方法的基本思想是首先假设一个流量， 分别用缓流和急流的迭代方法进行自由面迭代，如果上游缓流边界能很好地与下游急 流边界交汇并光滑地过渡到下游，则该流量即为正确流量：如果不能很好地交汇，则

需要重新假设流量进行试算。然而通过计算，他认为两种流态边界不能很好地汇合，
因而未给出最终成果。文献［１０９］从可变区域的概念出发，把决定自由水面形状的坐 标和各节点的流函数作为未知量形成方程组同时求解。文献『１ ｌｏ］用有限元法计算了

过坝水流，提出了表面流带的概念，将一维水流的临界流量概念和水面线计算方法应
用于二维，找出了过坝流量和自由表面之间的联系，克服了寻找临界点的困难，文献

［１１１】分别给出了用有限元泛函极值加权余量法离散格式，并用此格式求解过坝水流问
题的方法，得到了与实验较一致的结果。文献ｆ１１２］提出了一种求解具有未知能头的
过坝水流迭代方法，根据自由表面能量极值条件，建立了确定未知水头的数学关系，

从而得以实现水头与自由面同步迭代求解。这一方法既适用于高坝，也适用于低堰。

由此可见，用有限元法对过坝水流进行数值计算已是硕果累累，占有很重要的位 置。用该方法对流量函数和自由表面位置的计算也是相当成功的，所得结果能够基本
满足水利工程的设计要求。 边界元法起步较晚，在近２０年来才得到了较快的发展，最早应用于势流计算， 同时也是应用最为成功的领域，但把边界元法直接应用于过坝水流计算目前还不多

见，文献［１１３］求解了溢流坝顶自由出流问题，文献［１１４］求解了水力学中的不定边界
流场问题。

河海大学博士论文

边界拟合曲线坐标变换法在近年来被更多地应用到过坝水流的数值模拟，该方法 是通过边界拟合坐标变换，将不规则的物理流动区域变换成规则的矩形计算区域，同 时将控制方程作相应变换，然后采用比较简单的有限差分法进行求解。文献［１０６１干０ 用该方法计算了过坝水流问题，计算结果与实验吻合良好，文献［２３］Ｎ用曲线坐标变 换对溢流坝反弧段水流进行了模拟。 需要注意的是，关于过坝水流的计算主要还是以势流理论为基础，本章利用前面 建立的非正交网格系统中的ＳＩＭＰＬＥ算法对溢流坝过坝水流二维粘流场进行计算，采 用分区网格技术对计算区域进行网格划分，湍流模型采用ＲＮＧｋ一占模型，自由水面 线的形状及位置采用ＶＯＦ方法确定，分别对高水头溢流坝的水面线形状、坝面压力、 流量系数及坝面流速分布等进行计算。 ６．２计算模型 目前世界上广泛使用的坝面溢流曲线是由美国陆军工程兵团水道实验站（ＷＥＳ） 提供的ＷＥＳ曲线‘“”，该型曲线具有工程量小、无有害的负压及泄流能力大的特点， 一般表达式为：

日Ｙ－ｃ（≯
其中：爿一不包括行近流速水头的溢流坝上的水头高度 ｃ，６一与上游坝面倾斜坡度有关的系数，
ｘ，Ｙ一以溢流坝最高点为原点的坐标值。 在坐标原点上游，坝面曲线一般由两个或三个圆弧线组成，其中三圆弧曲线的性能优 于二圆弧曲线。

本文选用的坝面曲线为三圆弧ＷＥＳ曲线Ｘ‘”＝２．０川”Ｙ，见图６．１，６．２，其中
上游坝面坐标在表６．１中给出，溢流坝高Ｈ。＝７６．２米，设计水头Ｈ。＝９．１４米，计算 区域与网格划分见图６．３。

６７

河海大学博士论文

图６．１

三圆弧溢流坝上游坝面

图６．２三圆弧溢流坝坝面曲线 表６．１
ｘ ｌ Ｈａ
．０．Ｏ

三圆弧ＷＥＳ坝面上游坐标（坐标原点取在坝顶中心线上）
．０．０５ ．０．１０ ．Ｏ．１５ ．Ｏ．１７５ ．０．２０ ．Ｏ．２２

Ｙ｜Ｈ ｄ

０．Ｏ

０．００２５

Ｏ．０１０ｌ

０．０２３

０．０３１６

０．０４３

０．０５５３

ｘ ｆ Ｈｄ

．Ｏ．２４

．Ｏ．２６

．０ ２７６

．０．２７８

．Ｏ．２８

．０．２８１８

Ｙ｜Ｈ ｄ

０ ０７１４

０ ０９２６

０．１１ ５３

０．１１９０

０．１２４ｌ

０．１３６０

河海大学博士论文

翻６．３计算区域与网格划分

从图６－３中可以看出，整个流场分为三个区域进行网格划分，网格类型均为结构 化网格，其中Ｉ区网格数为７５００，ＩＩ区网格数为４０００，为保证坝面压力的计算精度， 对ＩＩＩ区网格进行加密，网格数为８０００。 入１３边界速度按设定流量计算给出，出口边界条件取物理量沿水流方向梯度为 零：所有的气体边界都定义为压力边界，即压力为０，湍动动能ｋ和耗散率ｅ按经验 公式确定，固体边界速度按无滑移条件给定，并引入壁面函数技术，自由面的位置利 用ＶＯＦ方法确定：时间步长△ｆ＝Ｏ．００２，迭代次数２００００。 ６－３自由水面线计算 水面线位置的计算按表６．２进行，水面线形状在图６．４中给出，与试验结果的比 较在图６．ｊ中给出，可以看出，计算出的水面线位置与试验值非常一致。
表５．２水面线计算
Ｈ。／ＨＪ
（初值）

０．５０

１．Ｏ

１．３３

１．５

Ｑ｛Ｑｄ ＝（Ｈ。／ＨＪ）１
日
６

０．３２９８７７

１．Ｏ

１．５７８２０２

１．９１３１３７

（计算值）
Ｈ。／Ｈｄ
Ｌ计算值）

８０，８５

８５．３８

８８．２２

８９．６７

０．５０８７

１．００４４

１．３１３

１．４７８

河海大学博士论文

图６．４不同坝顶水头下的水面线（红色为水）

８０

７５

７０

６０ ｌＯ ０ １０ ２０ ３０

Ｘ（１１１） （ａ）Ｈｅ／ＨＦ０．５

河海大学博士论文

９０

８５

８０

５
｝

７５

７０

６５

６０
０ １０

２０

３０

Ｘ（ｍ）

（ｂ）Ｈｅ／Ｈａ＝ｉ．０
９０

８５

８０

占７５
｝

７０

６５

６０ ０ １０ ２０

Ｘ（ｍ）

图６．５

水面曲线的形状

６．４水位一流量关系
文献［９７】根据实验结果，利用图解法给出溢流坝的水位一流量关系曲线计算公式

若＝（鲁）ｌ

６，该实验曲线可以作为数值计算结果的验证比较。

首先计算设计流量Ｑ＝ＣＬＨ。３”米３，秒，其中Ｃ＝ｏ．５５２４瓦?ｍ为设计流量系数，
其值根据实验回归曲线确定。Ｌ为有效坝宽，二维情况下取值为１。 为了确定溢流坝的水位与流量的关系，给定８种状态下的不同流量，计算出对应

７

河海大学博士论文

的坝顶水头高度，绘制出水位一流量关系曲线，并与实验结果进行比较，分别列于表
６＾３和图６ ６中，可以看出，本文的计算方法能够比较准确地给出溢流坝的水位与流

量关系，从而一方面可以在溢流坝设计时准确计算出设计流量，另一方面，能够在全
部溢流水头范围内，对理论结果或模型实验值进行校核。
表６．３水位一流量关系 Ｑ｜Ｑｄ
Ｈ。／Ｈｄ
（计算）
Ｏ，０５ Ｏ．１ ０．２ ０．４ ０．６ Ｏ．８ １．Ｏ １．３３

ｏ １５３１７３

ｏ．２２９７５９

０．３６１０５

０ ５７９８６９

０．７１１１６

０．８５３３９２

１．００４３７６

１

１８１６１９

Ｈ：／Ｈ
Ｉ

ｄ

０ １５３７６５

０ ２３７１３７

０．３６５７１６

０ ５６４０１１

０．７２６６８２

０ ８６９８２４

１．０

１

１９５１０８

（实验） ｌ相对误差
０ ０００５９２ ０ ００７３７８ ０．００４６６６

０．０１５８６

０．０１５５２２

０ ０１５４３２

０．００４３８

０ ０１３４８９

ｌ

４

１

２

ｌ

０

８

ｌ＇＠ａ

０

６

０

４

Ｏ

２

Ｏ

图６．６水位．流量关系曲线

６．５流量系数
流量系数是决定溢流坝泄洪流量的重要参数，对于同一坝面曲线，不同的坝顶水 头将具有不同的流量系数，为了便于比较，本文在给定流量Ｑ／Ｑｄ＝Ｏ．３３，１．ｏ，１．２９， １．５８，１．９１时计算出对应的实际坝顶水头分别为

Ｈｅ／Ｈｄ＝０．５０８７，１．００４４，１．１６６，１．３１３，１．４７８，根据公式ｍ２了南即可得到对应的流量

河海大学博士论文

系数，计算和实验结果分别在表６．４和６．５中给出，同时将结果绘于图６．７中，二者
基本一致。
表６．４流量系数（计算）

Ｑ｝Ｑｄ
Ｈ：｝Ｈ
卅
ｄ

Ｏ．３３

１．Ｏ

１．２９

１．５８

１．９１

０．５０８７

１．００４４

１．１６６

１．３１３

１．４７８

０．４５５６４７

０．４９７９５８

Ｏ．５１１５８５

０．５２５８３２

０．５３３６０４

表６ ５流量系数（实验）

Ｑ旧ｄ
Ⅳ，／ＨＪ

０．３３

１．Ｏ

１．２９

１．５８

１．９１

Ｏ．５０

１．０

１．１７

１．３３

１．５

门，

０．４５５

０．５０ｌ

Ｏ．５１３

０．５２６

０．５３６

１．６

１

４

１

２

口

要
士

１

０．８

０．６

０．４

０．２
０．５ ｍ ０．５５

０．６

图６．７流量系数

需要说明的是，关于坝顶水头高度的计算是在水库上游面按流体体积分数Ｆ＝Ｏ．５ 时对应的水面线高度进行的。
６．６坝面压力 坝面压力的计算是溢流坝水力设计的重要内容，为此，本文分别计算了 日。／Ｈ。＝Ｏ．５，１．Ｏ，１．１７，１．３３，１．５共５种情况下的坝面压力，绘制成曲线并与文献【１７】中的

实验点进行比较，见图６．８，压力分布的基本趋势是一致的，其中工程中关心的负压

河海大学博士论文

极值，其出现的位置与大小均与实验符合良好。
０ ５
ｏ ４

ｏ ３

。弋
；＼ ， Ｉ
ｆ Ｌ ，ｊ

、

０ ２ ０ ｌ

＼． ｊ
√．
，
Ｉ

＂▲Ｕ
Ｊ
—Ｉ －－Ｊ
、一 ——－一

昭

０ ｌ

要

一

２
３ ４ ５ ６ Ｔ

Ｉ １
＼。

』
，，一

上
、

ｆ一
。、

—，一Ｊ

≯’

＼＿．
ｌ ｌ

＜?
＼

、 、 ／
＼／
－０
１

“’“‘叮。ｒ

叼∞加吣加加叼叼

８ －Ｏ．了
Ｕ

１

Ｕ．７３

Ｕ ｂ

０

Ｔ

０ ９

２７Ｈｄ

图６．８坝面压力分布曲线（圆点为实验值）

６．７坝面水流速度分布
为考察溢流坝过坝水流的流速变化规律，在图６．９中分别给出了断面 ｘ＝０．０，２．５７６，５，１０，２０，３０，坝顶水头分别为Ｈ。／Ｈｄ＝Ｏ．５，１．０，１．５时速度分布曲线图

９０
??＋??Ｘ＝２ ５７６

８５
、‘

———ｏ—一Ｘ＝５
８０
７５
———●●－—－Ｙ＝１ ｎ

０
’

盈

ｋ

，

—ｔｒ。—一Ｘ＝２０ 一
一

７０ 一 占６５
）、

一】
６０ ５５

５０
● ４５

４０
Ｏ Ｏ
Ｖ

册

巧

∞

惦弘∽

河海大学博士论文

９０

９０

８５

８５

８０

７５

。蔫＿｜｜｜＿ＩＩ．。一

８０

佰 阳 ∞ ∞ 踮

７０

６５

６０

］
５７６ １０ １５

５５

５０

４５

一
０ ５

印

Ｘｌ０２ ?－?－－ｘ＝

，

蛎 蚰 ％

－…一?ｘ＝２．５
—－——－ａ＿－－—－Ｘ＝５ ———啼÷——一Ｘ＝１０

７６

４０

———■ｒ—一ｘ。２０

３５ ２０ ２ ５ ３０

０

５

１０

ｌ ５

２０

２５

３０

Ｖ（ｍ／ｓ）

Ｖ（ｍ／ｓ）

（ｂ）

（ｃ）

图５．９坝面水流速度分布（ａ）Ｈ。／Ｈｄ＝Ｏ ５（ｂ）Ｈ。／Ｈｄ＝１．０（ｃ）Ｈ。／Ｈｄ＝１．５

为进一步考察水流沿坝面的速度变化，给出了个断面上最大速度的沿程变化曲线 见图６．１０。从图５．９，６．１０可以看出，水流速度沿坝面逐渐增加，且不同水头下的最
大速度的差值沿程减小。为考察流线的整体变化情况，在图６．１ｌ中给出了 Ｈ。／Ｈ。＝１．０时整个流场内的流线图。

图６’ＩＪ

６．８本章小结

镒流坝流线矢量图

河海大学博士论文

第七章带闸墩溢流坝三维过坝水流数值模拟
带闸墩溢流坝三维过坝水流数值模拟的工作开展不多，主要原因计算区域增大， 闸墩区域形状复杂，并且包含自由液面问题，这些都使得数值计算较为困难。因此过 坝水流的研究工作更多地集中在二维流场计算。尽管二维过坝流场计算具有计算工作 量小、计算成本低，并能够在一定程度上指导工程设计，比如文献『１１６１通过计算过 坝水流与闸墩绕流两个二维模型的方法近似模拟了三维带闸墩的溢流坝流场，可以作 为解决三维过坝水流问题的一条途径，但是将二维数值模型应用于实际的带闸墩的溢 流坝过坝水流计算还有一定的差距。主要表现在闸墩的存在对水面线位置、坝面压力、 溢流流量等均产生显著的影响。文献ｉ１１７］利用有限体积方法对溢流坝三维湍流场进 行数值模拟，但没有考虑闸墩的影响，本文直接建立溢流坝面与闸墩的三维整体模型， 对三维湍流场进行数值求解，其中湍流模型仍然采用ＲＮＧｋ—ｓ两方程模型，自由水 面线的位置采用ＶＯＦ方法确定，针对不同的墩型和布置形式进行计算，比较准确地
给出了过坝水流的水面线位置及坝面压力。 ７．１闸墩型式

溢流坝坝面曲线选用三圆弧ＷＥＳ曲线ｆ¨８】ｚ‘”＝２．０Ｈ≯Ｙ，见图６，ｌ，６．２。 溢流坝设计的重要参数是流量计算，对流量影响最大的是闸墩，因此必须对包含 闸墩的溢流坝三维流场进行数值模拟。闸墩的作用是将溢流坝前缘分隔为若干个孔
口，并支承闸门、启闭机和桥梁等传来的载荷，闸墩的断面形状应尽量减小孔口水流 的侧收缩，使水流能够平顺通过，因此关键在于墩头和墩尾形状的选择，头部主要影 响侧收缩，从而影响泄洪流量，尾部主要影响下游流态，如可能产生水冠或冲击波，
影响下游的消能。闸墩墩头的形状及闸墩的主要尺寸在图７．１中给出。
Ｋ
／

，
乇 州 ｏ

Ｊ

半径０ １３３Ｈｊ

、

Ｏ ２ｆ
● —、‘

７只、
ｒ

坝项轴线Ｌ
Ｉ

ｊ§ 习 磊 薄丽彰 悟
０．２ｔ

７Ｈａ．Ｊ
７Ｉ
—

ｆ．０２Ｉ 『片。。Ｉ坝顶轴线｝．０２ｆ
’ｉ

７Ｈ．，、Ｉ
’

—』乙?

Ｃ、／、

ｌ

ｌ

’＿ｒ

ｌ

ｉ¨

２型
圈７．１闸墩的型式

３Ａ型４型

河海大学博士论文

７．２计算区域与网格划分 含闸墩溢流坝三维流场计算区域在图７．２中给出，可以看出在闸墩区计算区域的 几何形状最为复杂，具有双向曲率变化，这也正是研究该类问题的难点所在。为了保 证计算精度，对该区域进行网格划分需要特别注意，本文采用分区网格系统，将整个 计算区域划分为３个子区，见图７ｔ３。由于闸墩区域边界变形较大，采用结构网格无 法生成可用于计算的网格系统，因此在该区内采用非结构网格，同时在壁面处网格进 行加密：在Ｉ区和ＩＩ区采用结构网格，由于ＩＩ区处于自由表面区，为提高自由面的模 拟精度，同样进行了网格加密。此外。为了保证Ｉ区、ＩＩ区与闸墩区网格的平顺过渡， 在闸墩区的上游和左下方增加一部分区域归并到闸墩区构成ＩＵ区进行网格划分，三个 区的网格数目约为２７０００，１４０００，７００００，共计约１ｌ万个网格节点。在图７．４中给出 了坝面及闸墩区的网格图。

．１００

图７．２含闸墩溢涟坝流场计算区域

河海大学博士论文

图７．３

含闸墩溢流坝流场网格的分区划分

图７．４坝面及闸墩区的网格划分

河海大学博士论文

７．３计算初始条件与边界条件

坝高Ｈ＝７６．２米，设计水头Ｈ。＝９．１４米，设计流量伤＝ＣＬ蟛３”，Ｃ为设计流
量系数，上为坝的净宽＝坝宽．闸墩厚度，分别对２型墩和３Ａ型墩进行计算，其中墩 长均为１．７７Ｈｄ，墩厚分别为０．２６７Ｈｄ和０．１７３Ｈｄ，墩厚与闸宽之比均为０．２。假定坝 项水头高度及流量，分别对三种不同水头高度Ⅳ。／Ｈｄ＝０．５，１．Ｏ，１－３３下的自由水面线、 坝面压力、流量系数进行计算，计算初值见下表
表７．１水头高度及流量初值 Ｈ！／Ｈｄ
（初值）

Ｏ．５０

１．Ｏ

１．３３

Ｑ『Ｑ。 ＝（Ｈ。／Ｈｄ）‘６

０．３２９８７７

１．Ｏ

１．５７８２０２

入口边界的速度己知，ｋ，ｓ值按公式々：０．００５“。２，￡：燮给出，出口处的速度
）Ｐ

及ｋ和￡值为未知．计算根据与出口断面相邻的上游两个断面的变量值，线性外延求 得出口断面之值，作为下一次迭代的边界条件．计算过程中，每次迭代都需根据上一 次迭代计算结果修正出口边界条件，固体壁面上流速按非滑移条件给定，固壁附近的 ｋ和占采用用壁函数技术；自由表面上的压力取为大气压，而ｋ和￡采用零梯度条件， 即令自由表面网格中的ｋ和ｓ的值等于相邻的流体内部网格的ｋ和ｓ的值。计算开始 流动为非定常流，经过一段时间的计算，流场各水力要素均趋不变，流动趋于定常， 计算即告结束。 ７．４自由水面线 按照上节的计算条件分别对２型墩和３Ａ型墩水面线形状及位置进行计算，
７，４．１

２型墩水面线计算 图７．５中给出了沿２型墩闸孔中心线水面线位置计算与实验结果的比较，三条曲

线分别对应日／Ｈ。＝０．５，１．Ｏ，１．３３三种情况。二者是比较一致的，在坝顶水头较低时， 误差较大，随着坝顶水头的增加，计算与实验结果的一致性更好。与沿闸孔中心线的 计算相比，沿闸墩表面水面线位置计算结果与实验值更加接近（见图７．６），说明本计 算模型是准确、可靠的。

８０

河海大学博士论文

９０

８５

８０

— ５

ｈ ７０

７５

６５

６０ １０ —５ ０ ５ １０ １５ ２０ ２５ ３０

ｘ（ｍ）

图７．５沿２型墩闸孔中心线水面线位置比较
（从上至下分别对应Ｈｅ／Ｈｄ＝ｌ
９０
３３．１ ０．０．５）

８５

８０

．、

县７５
＞、

７０

６５

６０

１０

—５

０

５

１０

１５

２０

２５

３０

ｘ（Ⅲ）

图７．６沿２型墩表面水面线位置比较
（从上至下分别对应Ｈｅ／Ｈｄ＝Ｉ
３３，１．０，０ ５）

为了比较闸墩对水面线位置的影响，在图７．７和图７．８中分别给出了无闸墩、沿 闸孔中心线和沿闸墩表面的自由水面线位置的计算曲线，显然闸墩的存在抬高了自由 水面线的位置，且随着水位的增加，抬高的幅度逐渐增加。 而沿闸墩表面的变化则

是闸墩头部水面抬升，在闸墩的平行段水面低于闸孔中心线处的水面。 在图７．９中比较直观地给出了三种水头情况下通过２型墩闸孔的过坝水流。图

河海大学博士论文

７．１０中给出了不同高度处闸墩附近的流线矢量图。

９０

８５

８０

５
ｈ

７５

６５

６０ ｌＯ 一５ ０ ５

ｘ（ｍ３０

１５

２０

２５

３０

图７．７无闸墩与２型墩沿闸孔中心线的自由水面线
（从上至下分别对应Ｈｅ／Ｈｄ＝ｌ
３３，１

０，０．５）

８５一

占７５，
＾

７０

６０。—一
一１０ ０ ５ １０ １５ ２０ ２５

ｘ（面 图７．８沿２型墩闸孔中心线与沿闸墩表面的自由水面线
（从上至下分别对应Ｈｅ／Ｈｄ＝Ｉ
３３，１

０，０．５）

河海大学博士论文

（ａ）

（ｂ） 图７．９ Ｈ／Ｈｄ＝０．５，１．Ｏ，１．３３时２型墩过坝水流 （红色为水，黄色为闸墩和坝面）

（ｃ）

图７．１０日／Ｈ。＝１．０时２型墩过坝水流不同高度处的流线矢量图 （ａ）ｈ＝７７ｍ，（ｂ）ｈ＝７８ｍ，（ｃ）ｈ＝７９ｍ，（ｄ）ｈ＝８０ｍ
７．４．２

３Ａ型墩水面线计算 图７．１ｌ，７．１２中给出了设计水头Ｈ／Ｈｄ＝Ｏ．５，１．０，１．３３三种情况下沿３Ａ型墩闸孔

中心线和沿闸墩的水面线位置计算结果，并与无闸墩情况进行了比较。与图７．７，７．８ 比较可以看出，两种墩型对孔中心处水面线位置影响不大，为了比较３Ａ型墩与２型 墩对水面线位置的影响，在图７．１３中给出了沿两种墩型表面的水面线位置计算结果， ２型墩使墩前水面产生更高的抬升。

河海大学博士论文

９０

８５

８０

３
ｈ

７５

７０

６５

６０ ５ ０ ５
ｘ

，１０

１

５

２０

２５

３０

ＬｍＪ

图７．１１无闸墩与３Ａ型墩沿闸孔中心线的自由水面线
（从上至下分别对应Ｈｅ／Ｈｄ＝１．３３，１
９０ ０，０ ５）

８５

８０

５
）、

７５

７０

６５

６０ １０ —５ ０ ５

ｘ（拶

１５

２０

２５

３０

图７．１２沿３Ａ型墩闸孔中心线与沿闸墩表面的自由水面线
（从上至下分别对应Ｈｅ／Ｈｄ＝ｌ
３３，１

０．０．５）

河海大学博士论文

９０

８５

８０

占７５
ｈ

７０

６５

６０ １０ —５ ０ ５ １０ １５ ２０ ２５ ３０

Ｘ（ｍ）

酗７．１３沿２型墩及３Ａ型墩表面的自由水面线
（从上至下分别对应Ｈｅ，／Ｈｄ＝Ｉ．３３，１ ０，０．５）

７．５坝面压力 坝面压力计算，特别是坝面最小压力的数值及出现位置是溢流坝设计的重要参 数，在图７．１４和图７．１５中给出了三种不同水头下２型闸墩沿闸孔中心线及沿闸墩的 坝面压力分布曲线，在图７．１６和图７．１７给出了三种不同水头下３Ａ型闸墩沿闸孔中

心线及沿闸墩的坝面压力分布曲线，并与实验结果进行对比，可以看出，沿闸孔中心

线坝面压力的计算结果与实验值非常一致，只是在以／ｚ‘＝１．０时在坝顶后段计算值
偏低，而坝面最小压力值的大小及出现位置均与实验吻合：对于沿闸墩面的压力，总

的趋势是计算结果略低于实验值，并且随着水头的增加其偏差量增大。在 Ｈ。／Ｈ。，＝１．Ｏ、１．３３两种情况下沿闸墩的坝面上的最小负压出现位置与实验一致，均在
坝顶前０．１Ｈ。处。

将２型墩与３Ａ型墩的坝面压力计算结果进行比较可以发现，在相同的墩厚闸宽

比情况下，２型墩将使坝面出现更低的负压，当水头Ｈ。／Ｈ。＝１．３３时负压极值为７．３ 米水柱，超过了工程中建议的不低于一６．１米水柱以保证坝面不出现空蚀现象的要求，

而３Ａ型墩在相同水头下对应的最低负压仅为一４．６水拄，不会出现坝面空蚀现象。

河海大学博士论文

０．６ ０，４ ０，２

≈

鼍０
￡
－０，２
—０，４

—０，６ —０，８ －ｉ －ｏ，３一Ｏ，１ Ｏ，ｌ

Ｏ，３

Ｏ，５

０，７

０，９

１，１

１，３

Ｘ／Ｈｄ
图７．１４沿２型墩闸孔中心线上的坝面压力分布
Ｏ

６

０．４

０

２

＼
．Ｃ

Ｚ

０
２

一０

—０ ４
—０ ６

—０ ８ １ —０．３—０．１ ０．１ ０．３ ｏ．５ ０．７ ０．９ Ｉ．１ Ｉ．３

Ｘ／Ｈｄ 图７．１５沿２型闸墩表面的坝面压力分布
０．６ Ｏ．４ ０．２ 弓

姜０
工 一０．２ —０．４ —０．６ —０．８
—１

降丧釜≤。。銎琴 ■＇：～：ｉ！ｊ｝－‘！ｉ亨
０．１
０．３

－０ ３—０．１

０．５

０．７

０．９

１．１

１．３

Ｘ，，Ｉ－Ｉｄ

图７．１６沿３Ａ型墩闸孔中心线上的坝面压力分布比较
（从上至下分别对应Ｈｅ．，／ｌ－Ｉｄ＝０ ５，Ｉ．０，ｌ ３３）

河海大学博士论文

图７．１７沿３Ａ型闸墩表面的坝面压力分布比较
（从上至下分别对应ＨｅＪＨｄ＝０．５，１．０，ｌ ３３）

７．６流量系数 流量系数关系到工程的泄洪规模和具体布置，但流量系数随工程的体型而变化， 由于水流均为三维，其流态和边界条件都比较复杂，故在一般的手册中难以找到准确 的数据。在图７．１８中比较了墩厚闸宽比ｔ／ｂ＝０．２时２型墩和３Ａ型墩在三种不同水头 下的流量系数，显然，３Ａ型墩具有更大的流量系数。
ｍ ６

５

呲ｌＩＪ
ｎ ５

¨５
ｍ ４ Ｏ，４ ０．６ ０．８
ｌ

１．２

１．４

Ｈｅ／ｌｉｄ

图７．１８有闸墩影响的流量系数（ｒ／ｂ＝０．２）

７．７墩厚闸宽比对坝面压力及流量系数的影响
除坝型和闸墩型式对泄流流量、坝面压力、水面线位置产生较大影响外，闸墩在 坝面上的布置同样至关重要，特别是对坝面压力影响更为显著，因此为了考察闸墩布

置对坝面压力等的影响，针对３Ａ型闸墩，给定初值ＨｅＰｒ－Ｉｄ＝１．０时情况，计算了墩厚

ｆ与闸宽ｂ之比ｔ／ｂ卸．１４，０．２，Ｏ．５三种情况下坝面压力分布，对应的实际水头高度

河海大学博士论文

分别为０．９６８２７１，１．００７６５９，１．０１７５０５，分别给出了沿ｉＹｑ孑ｔ，ｅｅ心线和沿闸墩的坝面压 力分布曲线（图７，１９，７．２０），并给出了对应的流量系数曲线（图７．２１）。
０，６
０，４

０，２
－。

暮０
上

一０，２
—０，４

－０，６
—０，８ —１

—０．３—０，１

０，ｌ

０．３

０，５

０，７

０，９

ｉ，１

１，３

Ｘ／Ｈｄ 图７．１９

Ｈｅ／Ｈｄ＝１．０时不同墩厚闸宽比沿闸孔中心线的压力分布

０，６

０，４ ０，２
口

罩０
上 一０，２

．拨
一¨＿Ｖ

ｒ．＿?！ｔ钿．．． 垃 －１ｌｔ－１一一：

—－?一ｔ／ｂ＝０ ５ －－－ＩＩ－－－ｔ／ｂ＝０．２ ——?广＿ｔ／ｂ＝０ １４
ｒ—一二
－

—０。４ —０，６

—０，８
。 ． ：

一１
一０．３—０．１

；

；

０，１

０，３

０，５

０，７

０．９

１，１

ｌ，３

Ｘ／Ｈｄ

图７．２０

Ｈｅ／Ｈｄ＝１．０时不同墩厚闸宽比沿闸墩面的压力分布

８８

河海大学博士论文

０，５５

０。５３

０，５１

ｇ
０ ４９

０，４７

０．４５ ０．４ ０，６ ０，８

ｔ／ｂ 图７．２１

Ｈｅ／Ｈｄ＝１．０时不同墩厚闸宽比对应的流量系数

７．８本章小结 本章在二维过坝水流流场计算的基础上，直接建立溢流坝和闸墩的三维整体模

型，对三维过坝水流粘流场进行数值计算，其中采用了分区网格系统来处理复杂的几
何边界，针对２型墩和３Ａ型墩计算了闸墩型式对水面线形状、坝面压力、流量系数

等的影响，通过与实验结果的对比表明，本文建立的模型能够比较准确地给出过坝水 流的水面线位置及坝面压力，计算指出，闸墩的存在抬高了水面线的位置，在闸墩头
部尤其明显，３Ａ型闸墩相对于２型墩具有更大的流量系数和更小的坝面压力。此外， 不同的墩厚闸宽比对泄流能力也将产生显著的影响，结果表明，随着墩厚闸宽比的增 加，坝面压力降低，而当ｔ／ｂ＝０．２时溢流坝具有更大的流量系数。

河海大学博士论文

第八章结论与展望
本文对粘流场数值模拟技术的研究现状和发展动态作了比较全面的回顾和展望， 重点研究了网格生成技术、数值求解技术、湍流模型技术和动边界模拟技术，在此基
础上建立了模拟自由面湍流场的数学模型，并成功地应用到二维溢流坝和带闸墩的三

维溢流坝过坝水流场数值计算。

８．１本文完成的主要工作与结论 １）对网格生成技术进行了较深入的研究，建立了一种简单的、能够独立于网格生成 技术之外的边界正交化方法，给出了详细的数学模型并编制了计算程序；通过网格生 成算例表明，该方法能够大大改善网格在边界处的正交性；同时针对分区结构网格系 统建立了分区交界面处的数据结构和计算模型： ２）利用有限体积方法在非正交同位网格系统中建立了ＳＩＭＰＬＥ求解算法，并对非正 交同位网格系统中的边界条件、延迟修正技术及计算节点的梯度计算等专题内容进行
了深入讨论。通过对倾斜方腔驱动流、后台阶绕流、非对称平板间的圆柱绕流等典型

流场的数值模拟，表明本文建立的非正交同位网格系统中的ＳＩＭＰＬＥ算法是合理的，
具有更加广泛的适用范围：

３）研究比较了两种湍流模型的特点及应用范围，以水翼粘流场为例完成了ＲＮＧｋ一￡ 湍流模型与标准ｋ一￡模型的计算比较，表明ＲＮＧｋ一占模型的模拟结果更优； ４）研究了运动边界模拟技术中的ＶＯＦ方法，详细建立了一种利用斜直线构造自由面
的高精度方法，通过对典型运动界面的数值模拟表明，该方法能够比其它自由面重构
方法更加准确地确定运动界面的形状和位置：

５）完成了对二维溢流坝过坝水流进行计算，分析了不同坝顶水头下的水面线位置、

坝面压力及流量、流速等水力特性，通过与典型实验资料的对比，计算的水面线高度
及坝面负压极值均与实验值具有非常好的一致性。

６）直接建立溢流坝和闸墩的三维整体模型，对三维过坝水流粘流场进行数值计算， 针对闸墩的型式及在坝面的布置计算了闸墩对水面线形状、坝面压力、流量系数等设 计参数的影响，并将不同墩型与布置形式的计算结果进行比较，结果表明，闸墩的存
在抬高了水面线的位置，其中在闸墩头部尤其明显；墩型对流量系数和坝面压力影响 较大，３Ａ型闸墩相对于２型墩具有更大的流量系数和更小的坝面压力；不同的墩厚 闸宽比对泄流能力也将产生显著的影响，随着墩厚闸宽比的增加，坝面压力降低，而

河海大学博士论文

当ｆ／ｂ＝０．２时溢流坝具有更大的流量系数。

７）本文建立的三维自由面湍流场模型能够比较准确地模拟带闸墩溢流坝三维过坝水 流湍流场，可以针对不同的坝型、墩型、坝面布置形式以及不同的坝顶水头完成系列
水力计算，为泄水工程建筑物的设计提供可靠的分析依据。

本论文的主要创新点： １）提出并建立了网格边界正交化算法模型：
２）建立了非正交同位网格系统中的ＳＩＭＰＬＥ算法，并以此为基础建立了非对接分区 结构网格数值求解模型； ３）提出了一种精度更高的自由面重构算法模型，通过典型运动界面的模拟验证了本 算法是精确的：

４）建立了带闸墩溢流坝三维整体计算模型，首次完成了带闸墩溢流坝三维过坝水流
粘流场数值模拟，得出了一些有意义的结论。

８．２展望 虽然本文对带自由面湍流场数值模拟技术开展了一定的研究工作，并成功地应 用到泄水工程中三维自由面湍流场的数值模拟，得到了一些有意义的结论。但是由 于三维自由面湍流是一个十分复杂的流体运动现象，其复杂的运动规律和内在机理
需要作进一步的分析和研究，因此，本论文只是在几个方面做了一些探讨性的研究， 今后仍有大量的工作要做，具体内容包括： １）进一步完善非正交同位网格系统中的ＳＩＭＰＬＥ算法，使之适用于网格强非正交的
情形；

２）对基于ＶＯＦ方法的运动界面追踪技术中的输运模型和重构方法仍需作更深入的研

究，进一步提高计算精度，以便更好地模拟水流的破碎、空化等现象； ３）开展带闸墩溢流坝上下游流场整体模型数值计算，探讨泄洪流量、空化空蚀等复 杂现象的内在机理。

９ｌ

河海大学博士论文

致

谢

作者在攻读博士论文的过程中得到了许多的指导、支持和帮助，在此表示衷心的
感谢。

首先感谢我的导师汪德燧教授，正是在导师一步一步的悉心指导和督促下，才得
以顺利完成毕业论文，导师深厚的学术造诣、严谨的治学态度、踏实的工作作风给我 留下了深刻的印象并将使我受益终生，在论文即将完成之际，特向导师表示最诚挚的 谢意和最美好的祝愿。

感谢华东船舶工业学院副校长王自力教授和科技处长蒋志勇教授，他们在我的
工作中给予了非常大的帮助，使我克服了工作中的许多困难，为顺利完成论文提供
了良好的工作和学习条件。

感谢杨松林教授、朱仁庆副教授、姚震球副教授在论文工作中提供的无私帮助
和建议。

感谢师弟赖锡军，在与他日常的讨论交流中作者获得了许多有益的帮助和启

迪，他的聪颖与豁达也给我留下了深刻的印象。感谢同学童朝峰、陈雄波，在与他 们的交往过程中获得了许多帮助。
在论文撰写过程中，得到了研究生周林慧的帮助，在此表示感谢。 还要感谢我的家人，他们帮助我分担了许多家务，提供了大力的支持，正是在

他们的帮助和鼓励下，我才得以将更多的时间、精力和饱满的热情投入到论文工作
中：感谢母亲对我的关心和牵挂，感谢并深深怀念我的父亲。 最后再一次向关心、支持和帮助过我的所有同仁表示最诚挚的谢意。

河海大学博士论文

参考文献
１．Ｔｈｏｍｐｓｏｎ Ｊ．Ｆ．，Ｗｅａｔｈｅｒｉｌｌ Ｎ．Ｐ，Ａｓｐｅｃｔｓ ｏｆ ｎｕｍｅｒｉｃａｌ ｇｒｉｄ ｇｅｎｅｒｔｉｏｎ：ｃｕｒｒｅｎｔ ａｎｄ ａｒｔ．

Ａ１ＡＡ（ｐａｐｅｒ），９３—３５３９
２．Ｓｈｉｈ Ｔ．Ｉ—Ｅａｎｄ Ｂａｉｌｅｙ Ｒ．Ｔ，Ａｌｇｅｂｒａｉｃ ｇｒｉｄ ｇｅｎｅｒａｔｉｏｎ ｆｏｒ ｃｏｍｐｌｅｘ ｇｅｏｍｅｔｒｉｅｓ，

Ｉｎｔ．Ｊ．Ｎｕｍ

  本文关键词：三维自由面湍流场数值模拟及其在水利工程中的应用，由笔耕文化传播整理发布。