本文關(guān)鍵詞：三維自由面湍流場數(shù)值模擬及其在水利工程中的應用，由筆耕文化傳播整理發(fā)布。

河海大學 博士學位論文 三維自由面湍流場數(shù)值模擬及其在水利工程中的應用 姓名：王志東 申請學位級別：博士 專業(yè)：水力學及河流動力學 指導教師：汪德爟 20040322

河海大學博士論文

摘要
自然界中的流體流動問題一般都具有三維性和湍動性，對于最常見的流體一水則
同時還擁有自由表面，水利工程中繞閘墩的泄洪水流、船

舶工程中的船體繞流等均具 有上述流動的三種特性，這三種特性是計算流體力學（ＣＦＤ）領域非常重要的研究內(nèi) 容和方向。
．

本文針對粘流場數(shù)值模擬技術(shù)的研究現(xiàn)狀和發(fā)展動態(tài)作了比較全面的回顧與展 望，重點研究了網(wǎng)格生成技術(shù)、數(shù)值求解技術(shù)、湍流模型技術(shù)和動邊界模擬技術(shù)，在 此基礎上建立了模擬自由面湍流場的數(shù)學模型，并成功地應用Ｎ－維溢流壩和帶閘墩
的三維溢流壩過壩水流的數(shù)值計算。

本文完成的具體工作和主要結(jié)論如下： １）以代數(shù)網(wǎng)格生成方法為基礎提出了一種簡單的、可獨立于網(wǎng)格生成方法之外的邊 界正交化技術(shù)；針對分區(qū)結(jié)構(gòu)網(wǎng)格系統(tǒng)建立了分區(qū)交界面處的數(shù)據(jù)結(jié)構(gòu)與計算模型： ２）利用有限體積方法在非正交同位網(wǎng)格系統(tǒng)中建立了ＳＩＭＰＬＥ求解算法，對非正交 網(wǎng)格系統(tǒng)中的邊界條件、延遲修正技術(shù)及計算節(jié)點的梯度計算等專題進行了深入討 論。通過對傾斜方腔驅(qū)動流、后臺階繞流、圓柱繞流等典型流場的數(shù)值模擬，表明本 文建立的非正交同位網(wǎng)格系統(tǒng)中的ＳＩＭＰＬＥ算法是合理的； ３）研究比較了各種湍流模型的特點及應用范圍，以水翼繞流場為例完成了ＲＮＧｋ一￡ 模型與標準ｋ—ｓ模型的計算比較，結(jié)果表明ＲＮＧｋ—ｓ模型能夠更好地模擬流場的湍
流特性：

４）研究了運動邊界模擬技術(shù)中的ＶＯＦ方法，詳細建立了一種精度更高的自由面重構(gòu) 模型，通過對典型運動界面的數(shù)值模擬表明，該方法能夠更加準確地確定運動界面的
形狀和位置：

５）利用本文建立的非正交網(wǎng)格系統(tǒng)中自由面湍流場數(shù)學模型，對二維溢流壩過壩水 流進行了數(shù)值計算，分析了不同壩頂水頭下的水面線位置、壩面壓力及流量、流速等 水流特性，通過與典型實驗資料的對比，水面線高度及壩面負壓極值均與實驗值具有
非常好的一致性。 ６）通過建立溢流壩和閘墩的三維整體模型，完成了對三維過壩水流粘流場的數(shù)值計

算，針對閘墩的型式及在壩面的布置計算了閘墩對水面線形狀、壩面壓力、流量系數(shù) 等設計參數(shù)的影響，并將不同墩型與布置形式的計算結(jié)果進行比較。結(jié)果表明，閘墩

ＩＩ

河海大學博士論文

的存在抬高了水面線的位置，其中在閘墩頭部尤其明顯：墩型對流量系數(shù)和壩面壓力 影響較大．３Ａ型閘墩相對于２型閘墩具有更大的流量系數(shù)和更小的壩面壓力：不同
的墩厚閘寬比ｆ／６對泄流能力也將產(chǎn)生顯著的影響．隨著墩厚閘寬比的增加，壩面壓

力降低，而當ｔｌｂ＝０．２時溢流壩具有更大的流量系數(shù)。 ７）本文建立的三維自由面湍流場模型能夠比較準確地模擬帶閘墩溢流壩三維過壩水 流湍流場，可以針對不同的壩型、墩型、壩面布置形式以及不同的壩頂水頭完成系列
水力計算，為泄水工程建筑物的設計提供可靠的分析依據(jù)。

關(guān)鍵詞：數(shù)值模擬，網(wǎng)格，湍流模型，流體體積法，溢流壩，閘墩

ＩＩ

河海大學博士論文

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

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

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

ｎａｔｕｒｅ

ａｎｄ

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

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

ｐｒｏｊｅｃｔ，

ａｒｏｕｎｄ ｈｕｌｌ

ｏｎ，ａｌｌ

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

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

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

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

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

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

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

ｔｈｅ

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

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

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

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

ｏｎ

２Ｄ ｓｐｉｌｌｗａｙ

ａｎｄ

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

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

ｊｏｂ ａｎｄ ｍａｉｎ

ｃｏｎｃｌｕｓｉｏｎ

ａｒｅ

ａｓ ｆｏｌｌｏｗｓ：

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

ｔｈｅ ｂａｓｅ ｏｆ

ａｌｇｅｂｒａｉｃ鰣ｄ

ｇｅｎｅｒａｔｉｏｎ

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

ｓｔｒｕｃｔｕｒｅ

ａｎｄ

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

ｏｎ

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

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

（２）ＳＩＭＰＬＥ

ｏｆ
ａｓ

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

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

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

ａｎｄ

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

ｎｏｄｅｓ，ａｎｄ

ＳＯ

ｏｎ．Ｔｈｅ

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

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

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

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

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

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

ａｎｄ

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

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

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

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

ａ

河海大學博士論文

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

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

ｔｈｒｏｕｇｈ

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

ｏｖｅｒ

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

ＢＳｅ

ｏｆ

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

ｆｒｅｅ

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

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

ｄａｍ

ａｎｄ

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

ｏｎ

ｔｈｅ ｓｕｃｈ

ｆｌｕｘ、 ｔｈｅ

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

ａｒｅ

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

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

ｄａｍ

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

（６）３Ｄ

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

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

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

ｓｐｉｌｌｗａｙ ｉｓ

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

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

ａＳ

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

ｂｒａｋｅ ａｎｄ

ｄｉｓｐｏｓａｌ

ｄａｍ

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

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

ｈｉｇｈｅｒ

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

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

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

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

ａｎｄ

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

ａｎｄ

ｐｒｅｓｓｕｒｅ

ｏｎ

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

ａｎｄ

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

（７）Ｔｈｅ
ｂｒａｋｅ

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

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

ａｎｄ

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

ｏｖｅｒ

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

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

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

ａｎａｌｙｓｉｓ

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

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

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

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

ｏｆＦｌｕｉｄｓ，

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

Ｖ

河海大學博士論文

．１上?

聃

——１一

日

本論文以水利工程中泄水工程的水力設計為背景，研究了自由面湍流場數(shù)值模擬
技術(shù)，并將之應用于帶閘墩溢流壩三維過壩水流場的數(shù)值計算。

隨著計算機科學和數(shù)值計算方法的發(fā)展，數(shù)值模擬在科學技術(shù)研究中具有更大的 優(yōu)勢和發(fā)展空間，一方面可以確定物理問題整體舶和局部的細致行為，另一方面在解 決實際工程問題的過程中便于進行方案優(yōu)選，從而降低費用，縮短周期，此外還可以 替代一些昂貴的、危險的甚至難以實現(xiàn)的實驗，因此數(shù)值模擬在整個科學技術(shù)進步中
起到非常重要的作用。
｝

然而，利用數(shù)值模擬技術(shù)解決實際工程問題必須要求數(shù)值模擬過程中的每一個環(huán) 節(jié)．如：數(shù)學模型、控制條件、數(shù)值算法等都是高效的ｊ先進的，并經(jīng)過充分驗證和 確認的，它們決定了數(shù)值模擬成果最終的可靠性和實用性。 在前人研究工作的基礎上，本論文對數(shù)值模擬的關(guān)鍵技術(shù)：網(wǎng)格生成技術(shù)、數(shù)值 求解技術(shù)、湍流模型和動邊界模擬等進行了較深入的研究，建立了非正交同位網(wǎng)格系 統(tǒng)中自由面湍流場數(shù)學模型，并成功地應用于二維溢流壩和帶閘墩的三維溢流壩過壩
水流數(shù)值模擬。
．

通過本論文的工作，取得了以下研究成果： １）提出了一種簡單的、可獨立于網(wǎng)格生成方法之外的網(wǎng)格邊界正交化方法； ２）利用分區(qū)結(jié)構(gòu)網(wǎng)格技術(shù)和非正交同位網(wǎng)格系統(tǒng)中ＳＩＭＰＬＥ算法建立了求解復雜流 場的數(shù)學模型，提高了計算效率和計算精度； ３）詳細建立了一種精度更高的自由面重構(gòu)方法，能夠更加準確地確定運動界面的形
狀和位置：

４）利用以上模型首次計算了帶閘墩溢流壩三維過壩水流粘流場。針對不同的壩頂水 頭、閘墩的型式以及在壩面上的布置完成了系列水力計算，為泄水工程建筑物的水力 設計提供可靠的分析依據(jù)。

學位論文獨創(chuàng)性聲明：

本人所呈交的學位論文是我個人在導師指導下進行的研究工作 及取得的研究成果。盡我所知，除了文中特別加以標注和致謝的地方
外，論文中不包含其他人已經(jīng)發(fā)表或撰寫過的研究成果。與我一同工

作的同事對本研究所做的任何貢獻均己在論文中做了明確的說明并 表示了謝意。如不實，本人負全部責任。

論文作者（簽名）：
（注：手寫親筆簽名）

￡色：壘
墜＆：堅

旦壘年

３

月

衛(wèi)２日

學位論文使用授權(quán)說明： 河海大學、中國科學技術(shù)信息研究所、國家圖書館、中國學術(shù)
期刊（光盤版）電子雜志社有權(quán)保留本人所送交學位論文的復印件或

電子文檔，可以采用影印、縮印或其他復制手段保存論文。本人電子
文檔的內(nèi)容和紙質(zhì)論文的內(nèi)容相一致。除在保密期內(nèi)的保密論文外， 允許論文被查閱和借閱。論文全部或部分內(nèi)容的公布（包括刊登）授

權(quán)河海大學研究生院辦理。

論文作者（簽名）：
（注：手寫親筆簽名）

互盤壘

竺年

３月２工日

河海大學博士論文

符號表
松弛因子，控制網(wǎng)格邊界正交化的程度
Ｐ，Ｑ
ｐ

阻尼參數(shù)，控制網(wǎng)格正交性的衰減程度
流體密度 流體運動的速度分量

Ｕ，（ｆ＝１，２，３）

Ｐ，
－“

粘性流體中的法向應力
粘性流體中的切向應力
流體的動力粘度

“

Ｒｅ＝ＵＬ／ｙ Ｕ
上

雷諾數(shù)

特征速度
特征長度 流體的運動粘度

＾
膏
占

渦粘系數(shù) 湍動能 湍動能耗散率
壁面切應力

Ｆ Ｈｓ

流體體積分數(shù)
溢流壩高 溢流壩頂設計水頭 溢流壩頂水頭
溢流壩流量

Ｈ１

日。
０
，竹

溢流壩流量系數(shù)
溢流壩面壓力水頭
閘墩厚度
閘孔寬度
．

＾
，

ｂ

ＩＸ

河海大學博士論文

第一章緒論
１．１數(shù)值模擬技術(shù)的工程研究意義
數(shù)值模擬、理論分析和模型試驗是研究流體力學（水力學）的三種主要手段， ～般來說，這三種研究手段和方法必須相互配合，相互補充，相互促進，才能共同 推進流體力學學科的發(fā)展并解決各種工程實際問題。

理論分析方法是在研究流體運動規(guī)律的基礎上提出各種簡化流動模型，建立各
種控制方程，在一定的假使和條件下，經(jīng)過一系列的解析推導和運算得到問題的解 析解，理論分析的許多方法仍是目前解決實際問題，主要是在初步設計階段中常常 采用的方法，但理論分析法常常無法用于研究復雜的、以非線性為主的流動現(xiàn)象。

模型試驗方法是研究流動機理、分析流動現(xiàn)象、探討并獲得流動新概念，推動
流體力學發(fā)展的主要手段，其不足是完成一個完整的試驗過程需要解決一系列復雜 的技術(shù)問題，所需周期長、費用高。 數(shù)值模擬方法的特點是在用比模型試驗所需花費少得多的情況下給出流場內(nèi)部 細節(jié)的詳細描述。并且只要數(shù)值模擬的數(shù)學提法（控制方程、邊界條件、數(shù)值算法 等）是正確的，則可以在較廣泛的流動參數(shù)（如雷諾數(shù)，空化數(shù)等）和物理設計參

數(shù)范圍內(nèi)較快地給出流場的定量結(jié)果，不受試驗中固有約束條件的影響。因此隨著 計算流體力學和計算機科學的迅速發(fā)展，數(shù)值模擬在科學技術(shù)研究中將具有更大的 優(yōu)勢和發(fā)展空間，主要體現(xiàn)在：①數(shù)值模擬在某種意義下比理論和實驗研究對流體
運動過程的認識更加深刻、細致，不僅可以得到運動的結(jié)果，而且可以了解整體的

和局部的細致行為；②在解決實際工程問題的過程中，利用數(shù)值模擬可以選擇最優(yōu) 的設計方案或?qū)嶒灧桨�，從而減少實驗次數(shù)，周期及費用：③數(shù)值模擬可以從理論 上探索流體運動新的現(xiàn)象和規(guī)律．而且可以替代一些昂貴的、危險的甚至難以實現(xiàn)
的實驗，比如水壩坍塌造成的災害預報、水下爆炸與核爆炸的過程和效應、大氣運 動現(xiàn)象等。因此數(shù)值模擬在流體力學研究乃至整個科學技術(shù)進步中起到非常重要的
作用。

雖然數(shù)值模擬在流體力學研究中占有如此重要的地位，然而要使其一方面能夠為
工程設計提供高質(zhì)量、短周期、可靠的分析設計依據(jù)，另一方面為探索流動新概念、 新機理提供新途徑，則要求數(shù)值模擬過程中的每～個環(huán)節(jié)，如：數(shù)學模型、控制條件、 數(shù)值算法等都是高效的、先進的，并經(jīng)過充分驗證和確認的，它們決定了數(shù)值模擬成 果最終的可靠性和實用性。

河海大學博士論文

１．２粘流場數(shù)值模擬關(guān)鍵技術(shù)研究現(xiàn)狀 １．２．１網(wǎng)格生成技術(shù)
網(wǎng)格生成技術(shù)是計算流體力學（ＣＦＤ）發(fā)展的一個重要分支，在ＣＦＤ高度發(fā)展的

美國，網(wǎng)格生成所需的人力時間占計算任務全部人力時間的６０％，因此將ＣＦＤ技術(shù)進
一步推向工程應用，首先必須完善網(wǎng)格生成技術(shù)“。。對于復雜程度一般的計算區(qū)域邊

界，采用結(jié)構(gòu)網(wǎng)格系統(tǒng)完全能夠滿足工程要求的精度，由于算法相對簡單，因而目前 仍被廣泛采用。常用的三維結(jié)構(gòu)網(wǎng)格生成方法大致可分為代數(shù)生成法“。、橢圓微分方 程生成法”’”和雙曲微分方程生成法。”等三類。 代數(shù)生成法是利用已知的邊界值進行中間插值來生成網(wǎng)格，其核心是確定插值函
數(shù)，晟常用的是“超限插值”（Ｔｒａｎｓｆｉｎｉｔｅ ｉｎｔｅｒｐｏｌａｔｉｏｎ），內(nèi)部網(wǎng)格節(jié)點的分布 可以由邊界上的伸縮函數(shù)來控制，通過一些專門的、不是十分復雜的技術(shù)能夠?qū)崿F(xiàn)網(wǎng) 格在邊界的正交或近似正交“““，因此，對于幾何形狀不是十分復雜的流場，代數(shù)生

成方法不失為一種簡單有效的方法。 橢圓型網(wǎng)格生成法是通過求解橢圓型微分方程來實現(xiàn)網(wǎng)格的生成，該方法的最
早代表為著名的ＴＴＭ方法”’，即通過求解Ｌａｐｌａｃｅ方程來生成區(qū)域網(wǎng)格，該方法的 優(yōu)點是所生成的網(wǎng)格是光滑的，能夠處理復雜的幾何邊界，同時若在方程右端增加

源項交成泊松方程則可以控制網(wǎng)格分布的疏密及傾斜度，其缺點是這樣的源項很難 確定，且較難實現(xiàn)內(nèi)部節(jié)點的控制。目前最常用的源項控制方法是根據(jù)邊界正交性 和網(wǎng)格間距的要求壹接推導出右端項的表達式”’…，橢圓型網(wǎng)格生成方法是最常用
的一種單域貼體網(wǎng)格生成方法。 雙曲型網(wǎng)格生成法是通過求解雙益型微分方程來實現(xiàn)的，該方法的優(yōu)點是能夠保 證網(wǎng)格在邊界處的正交性，容易控制兩格的尺寸，且計算時間比橢圓方法少卜２個量

級，缺點是邊界上的任何不連續(xù)都會傳播到內(nèi)部網(wǎng)格節(jié)點，外邊界位置不能事先給定。 但隨著對一些敏感問題，如邊界條件、物面不連續(xù)和角點的處理等所開展的討論和改
進。提高了該方法的適應能力，在近年來得到了比較廣泛的應用”’”‘“。。

隨著邊界復雜程度的提高，生成單域貼體計算網(wǎng)格更加困難，同時網(wǎng)格質(zhì)量不 能得到保證。影響最終的數(shù)值計算結(jié)果，為此人ｆ｜、］提出了分區(qū)結(jié)構(gòu)網(wǎng)格、非結(jié)構(gòu)網(wǎng) 格和結(jié)構(gòu)／非結(jié)構(gòu)混合網(wǎng)格技術(shù)。 分區(qū)結(jié)構(gòu)網(wǎng)格系統(tǒng)是目前比較成熟的、并經(jīng)常采用的構(gòu)造復雜區(qū)域網(wǎng)格的方 法，分區(qū)結(jié)構(gòu)網(wǎng)格方法的基本思想是根據(jù)邊界的特點將整體流場劃分為若干個子
區(qū)，對每一個子區(qū)分別建立網(wǎng)格，并在其中對流動控制方程求解。但在區(qū)域劃分對

需遵循以下原則：盡量使每個子域的邊界簡單以方便結(jié)構(gòu)網(wǎng)格生成；各子區(qū)的大小
也盡量相近以實現(xiàn)計算負荷的平衡，這一點對并行計算尤為重要。常用的分區(qū)網(wǎng)格

２

河海大學博士論文

方法有：對接網(wǎng)格方法和重疊網(wǎng)格方法。
國外分區(qū)網(wǎng)格生成技術(shù)于８０年代興起并逐漸形成實用軟件，具有代表性的有

ＥＡＧＬＥ…。和ＧＲＩＤＧＥＮ““。前者是以橢圓型網(wǎng)格生成理論為基礎，兼有代數(shù)生成方法 （超限插值）和表面網(wǎng)格生成技術(shù)：后者是以代數(shù)網(wǎng)格生成方法為基礎，功能不如 前者完善。到了９０年代，二者又補充了許多新的技術(shù)，如基于弧長的超限插值函 數(shù)，提高正交性的Ｈｅｒｍｉｔ三次型函數(shù)，良好的用戶界面（６ＵＩ）和加強的與ＣＡＤ接 口能力等，從而發(fā)展成為功能很強的網(wǎng)格生成系統(tǒng)。其基本思想是依靠超限插值方 法生成表面和空間網(wǎng)格，然后由橢圓方法光順，通過調(diào)整代數(shù)方法中的Ｈｅｒｍｉｔｅ超

限插值或橢圓方法中的源項控制函數(shù)來使網(wǎng)格線在交界面處的斜率具有連續(xù)性，以
保證交界面兩側(cè)網(wǎng)格的正交性。 非結(jié)構(gòu)網(wǎng)格技術(shù)于８０年代末９０年代初得到迅速發(fā)展“““。，其基本思想是：由 于四面體是三維空間最簡單的形狀，任何空間區(qū)域都可以被四面體單元所填滿，即

任何空間區(qū)域都可以被以四面體為單元的網(wǎng)格所劃分。顯然，非結(jié)構(gòu)網(wǎng)格方法舍去
了網(wǎng)格節(jié)點的結(jié)構(gòu)性限制，節(jié)點和單元的分布是任意的，能夠較好地處理復雜邊界。 非結(jié)構(gòu)網(wǎng)格方法在其生成過程中都采用一定準則進行優(yōu)化判定，因而能生成高質(zhì)量 的網(wǎng)格，且很容易控制網(wǎng)格的大小和節(jié)點的疏密度，一旦在邊界上指定網(wǎng)格分布， 在兩個邊界之間即可自動生成網(wǎng)格，因此非結(jié)構(gòu)網(wǎng)格系統(tǒng)在近年來受到了很大的重

視，有了較大的發(fā)展，如文獻［１６］針對二維粘流場構(gòu)造出適應于非結(jié)構(gòu)網(wǎng)格的有限 體積時間推進格式；文獻［１７］利用點自動生成的Ｄｅｌａｕｎａｙ三角化方法構(gòu)造非結(jié)構(gòu)
網(wǎng)格，同時處理成自適應網(wǎng)格形式，能夠計算帶移動邊界的瞬時流動問題等。非結(jié) 構(gòu)網(wǎng)格的主要缺點是所需計算內(nèi)存較大，計算效率較低。 結(jié)構(gòu)／非結(jié)構(gòu)混合網(wǎng)格系統(tǒng)則兼顧了二者的優(yōu)點，對于幾何形狀復雜的三維問 題可以由近邊界處的非結(jié)構(gòu)網(wǎng)格和遠離邊界的結(jié)構(gòu)網(wǎng)格構(gòu)成，具有靈活模擬復雜邊
界和計算效率高的特點。

目前分區(qū)結(jié)構(gòu)網(wǎng)格的生成方法和并行計算方法、結(jié)構(gòu)／非結(jié)構(gòu)混合網(wǎng)格系統(tǒng)是計 算具有復雜幾何形狀流場的主要方法，其進一步的發(fā)展方向是減少工作量、實現(xiàn)自動

生成和自適應加密，具有良好的人機界面及可視化，強調(diào)更有效的數(shù)據(jù)結(jié)構(gòu)等，最終
的目的是能夠快速生成適用于粘流場計算的高質(zhì)量計算網(wǎng)格。

今后，網(wǎng)格生成技術(shù)的研究重點應當是網(wǎng)格的專項技術(shù)研究，如：邊界正交化技
術(shù)、分區(qū)結(jié)構(gòu)網(wǎng)格及并行計算技術(shù)、網(wǎng)格的自適應加密技術(shù)以及動態(tài)網(wǎng)格技術(shù)等。

河海大學博士論文

１．２．２數(shù)值求解技術(shù) 數(shù)值求解技術(shù)實際上是一種離散近似的計算方法，各種各樣的流體力學問題必須
從給定的微分方程或基本定律出發(fā)，建立物理上合理、數(shù)學上適定、能夠在計算機上 進行計算的離散化數(shù)學模型。大多數(shù)數(shù)值求解技術(shù)的基本思想可歸納為：把原來在時 間和空間上連續(xù)的物理量的場（如速度場、溫度場等），用有限個離散點上的值來代 替，然后按一定的方式建立關(guān)于這些值的代數(shù)方程組并求解，從而獲得物理量場的近 似解。目前流場計算的主要數(shù)值求解方法有：有限差分法（ＦＤＭ）、有限元法（ＦＥＭ）、

有限分析法（刪肘）、邊界元法（ＢＥＭ）、有限體積法（ＦＶＭ）等，其中前四種數(shù)
值計算方法的特點及適應范圍在眾多的文獻中均可查閱［１８，１９，２０，２ｌ，２２，２３］，下
面重點介紹有限體積法。

有限體積法由ＭｃＤｏｎａｌｄ在１９７１年首次用于求解二維歐拉方程，１９７２年被 Ｐａｔａｎｋｅｒ等人用于ＳＩＭＰＬＥ算法計算恒定不可壓流。有限體積法的基本思想是：將計 算區(qū)域劃分成一系列規(guī)則或不規(guī)則的單元或控制體積，并使每一個控制體積包含一個 網(wǎng)格節(jié)點，將待解的微分方程對每一個控制體積分得到一組離散方程，其中的未知數(shù) 是網(wǎng)格節(jié)點上因變量的值，首先計算出通過每個控制體邊界沿法向輸入（輸出）流量 的動量通量，然后對每一個控制體分別進行流量及動量平衡計算，便可得到計算時刻

各個控制體內(nèi)流動參數(shù)的平均值（平均流速、平均壓強等）。因此，有限體積法實際 上是對推導原始微分方程所用控制體方法的回歸，與ＦＤＭ和ＦＥＭ的數(shù)值逼近相比 其物理意義更加直接明確，因而在流場計算中獲得了廣泛的應用“。。’“�！�。。
若邊界通量的計算只使用前一個時刻的數(shù)值，為顯式ＦＶＭ，反之，當涉及到下 一個時刻的數(shù)值時，則為隱式ＦＶＭ。此外還有半隱式ＦＶＭ“““。。

１．２．２．１數(shù)值方法對離散格式的要求
數(shù)值求解首先要對微分方程進行離散，然而不同的離散格式將帶來完全不同的 結(jié)果，因為同一微分方程不同的離散格式并不等價，也不一定能保持原來微分方程 的性質(zhì)，對漸變光滑的情況其結(jié)果可能相差不大，對迅變以致間斷的情況其解則可

能相差很大。因此利用ＦＶＭ在建立和使用離散格式時需要滿足以下條件才能保證 數(shù)值求解的精度和可靠性。①相容性：即當時空步長趨于零時截斷誤差（離散方程 與原始方程之差）為零；②守恒性：是指在每一個空間有限的控制體積內(nèi)保持質(zhì)量 和動量平衡，具有守恒性的數(shù)值格式在計算過程中不會產(chǎn)生守恒誤差，并且計算精
度和穩(wěn)定性較好，同時只有守恒格式才能處理計算空間中可能產(chǎn)生的間斷面，即流

場中可能產(chǎn)生的激波：③逆風性：即空間導數(shù)離散時根據(jù)流速方向選擇差分方向，
從而使粘性的耗散是合理的：④無虛假振蕩：是指在計算間斷（或解梯度大）時要

河海大學博士論文

求格式具有足夠耗散，不出現(xiàn)虛假振蕩，當然，無假振要求與計算精度的要求是相

互矛盾的，目前大多采用冊格式計算有大梯度和間斷的解；⑤高分辨率地捕捉間
斷：是指要求人工粘性和格式粘性的強度與解的光滑性相適應，在光滑區(qū)粘性自動

減少以提高精度，在間斷點附近或梯度大的區(qū)域又自動加強，以防假振并達高分辨 率，目前一些優(yōu)良格式的間斷分辯率可達二到三個網(wǎng)格：⑥計算穩(wěn)定性：是指當數(shù)
值解在受到微小的擾動后不發(fā)生突變現(xiàn)象，目前很難通過數(shù)學推導給出實用的穩(wěn)定 性條件，應用中幾乎總是通過數(shù)值試驗來確定時間步長，對時間步長起決定作用的

因素有時間離散格式、網(wǎng)格的規(guī)則性及邊界條件的處理等；⑦解的收斂性：數(shù)值格
式的準確解與原微分方程的真解之差成為離散化誤差，當時空步長趨于零時，離散 誤差亦趨于零，則稱格式具有收斂性，顯然相容性只保證數(shù)值格式收斂于原微分問

題，更重要的是數(shù)值解收斂于真解；⑨精度：數(shù)值解的精度取決于算法、網(wǎng)格、邊
界處理等諸多因素，除網(wǎng)格形狀要求規(guī)則外，還要求網(wǎng)格光滑，步長漸變，否則要

按照步長比進行校正：⑨健全性及通用性：健全性是指能保證數(shù)值解不園所使用的 網(wǎng)格而異，例如可選擇一個簡單的明渠均勻流問題，所用格式在不同網(wǎng)格上的計算
結(jié)果必須一致，此外格式中不含專門為計算而引入的、需要人工加以調(diào)整的參數(shù)（如 人工粘性系數(shù)），目的是保證結(jié)果的客觀性，有文獻預言，隨著格式其它性能的完 善，未來十年內(nèi)對格式的主要要求為健全性；通用性則要求格式適用于各種流態(tài)（恒 定流／非恒定流，低速流／高速流，連續(xù)／間斷）、任意幾何形狀，多種網(wǎng)格布置、多 種邊界條件等盡可能多的情況。 能夠滿足上述大多數(shù)要求的格式可以稱作為高性能計算格式”…。 為了達到上述要求，首先必須了解各種離散格式的優(yōu)缺點及其改進方法。下面 對微分方程目前主要的空間和時間離散方法進行綜述。

（１）空間離散化方法： ①中心格式：中心格式的原理是泰勒級數(shù)展開式，其最大優(yōu)點是計算簡單、精度
高，有算例表明，在加粘后其計算量仍只是二階ＴＦＤ格式的］／５，因而成為計算流 體動力學的代表方法之一，常用來計算光滑流動；缺點是計算解在計算空間中會產(chǎn)

生假性振蕩，同時穩(wěn)定性較差，目前實際計算中最有代表性的中心型數(shù)值方法是
Ｊａｍｅｓｏｎ““等首先提出的。 ②逆風格式：逆風格式是根據(jù)流速方向來選擇差分方向，優(yōu)點是穩(wěn)定性較好，如 二階單側(cè)逆風格式允許的ｆ為同樣點數(shù)的對稱格式的兩倍，對急流和出流邊界必須 或經(jīng)常使用逆風格式，入流邊界使用單側(cè)逆風格式是不可取的，但是，一階逆風的 耗散往往太強，精度低，二階線性逆風格式在計算間斷解時常產(chǎn)生和中心格式相似

的振蕩，這就要求在經(jīng)典逆風格式的基礎上研究非線性的單調(diào)性保持或ＴＦＤ格式。

河海大學博士論文

文獻［３２，３３］以中心格式和逆風格式為基礎建立了通量修正格式，該格式既防止假
性振蕩又可達到高階精度，且格式較為健全。

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

Ｄｉｍｉｎｉｓｈｉｎｇ）的含義是總變差減小，對于微分方程

的數(shù)值解“（一）來說，其總變差對應網(wǎng)格點的變差為：ＴＶ（ｕ（?））＝∑ｌｕ（ｘ，）－－Ｕ（Ｘ。）Ｉ
ｉ＝ｌ

如果對于不同時間步長的近似解有：ＴＶ（ｕ”、）≤ＴＶ（ｕ“） 則稱總變差是減小的，即具有ＴＶＤ性質(zhì)。

ＴＶＤ格式的建立是以方程罷＋罷：０的柯西問題為基本對象，在間斷形成前， 拼盤
解的全變差Ｚｙ保持常數(shù)，以后ｒ礦下降，同時在解的遞推過程中不產(chǎn)生新的極值，
原有極值也不增加，從而初值的單調(diào)性得以保持，據(jù)此可建立嚴格意義下的ＴＶＤ格
式。

建立二階ＴＶＤ格式的途徑有：修改通量（在粘性過大的～階逆風格式中加入

受限制的反擴散通量）：通量限制（通過引入限制因子對反擴散項加以限制以提高
精度的階數(shù)）：變量插值（在變量插值時對坡度進行限制）［１８，３０，３４Ｊ。

④ＮＮＤ格式 ＮＮＤ格式即為無波動、無自由參數(shù)的高分辨率格式，由張涵信所構(gòu)造““，其 基本出發(fā)點是在Ｅｕｌｅｒ方程和Ｎ．ｓ方程的右端附加一個三階導數(shù)項，并使該項的系
數(shù)在激波上游為正，在激波下游為負，以此來抑制差分解在激波上下游的非物理波 動，格式中不包含任何需要人工調(diào)整的參數(shù)，并具有二階精度的ＴＩｅＤ性，對應不同 的時間推進形式，可以得到不同的顯式全離散格式ＮＮＤ—ｌ到ＮＮＤ一５，并可構(gòu)造

出隱式時間推進無條件穩(wěn)定的ＮＮＤ格式，其中ＮＮＤ一２到ＮＮＤ一５格式對時間均 為二階精度。因此ＮＮＤ格式以其準確可靠，程序簡單，計算量小等特點在科學研 究及工程計算領域得到了廣泛的應用并取得了良好的效果。 需要說明的是中心格式和逆風格式可以相互轉(zhuǎn)化，中心格式添加人工粘性后可 解釋為逆風格式，中心格式也可改寫為逆風格式加上反擴散項，對后者加以限制即 構(gòu)成ＴＶＤ格式。反之，逆風格式也寫作中心格式加耗散項，對后者加以限制也構(gòu) 成ＴＶＤ格式。因此兩類格式的主要區(qū)別在于如何處理耗散，中心格式的人工粘性
和逆風或ＴＶＤ格式的粘性是相通的，故可以相互借鑒。

（１）時間離散化格式
時間坐標軸與空間坐標軸不同，它是單向的，只能是過去影響將來，時閶積分

河海大學博士論文

必須是采用相對ｒ。的單側(cè)格式，且往往與空間離散化分別對待。時間積分主要分為 顯式和隱式兩類。顯式是指可以由已知的ｆ。時的初值直接計算出ｔ。時刻的未知解，

不必求解方程組。其優(yōu)點是每個時間步的計算工作量和存儲量均較小，且程序簡單， 缺點是時間步長受到穩(wěn)定性要求的限制，當局部網(wǎng)格很密時，時間步長必須非常小
而使整個計算次數(shù)很大。隱式是指可同時由，。和，。時的解逼近空間導數(shù)和源項，

故通常需要求解非線性微分方程組。其優(yōu)點是在線性分析時往往是無條件穩(wěn)定的，
計算時間步長可以取較大的值從而提高總體效率。缺點是需要對每一個時間步長求 解方程組，因而計算工作量和內(nèi)存占用量均很大，程序設計也比較復雜，主要用于 恒定流或解的時間變率小、流動時間歷程長、高雷諾數(shù)流等情況。

應當注意，時間離散格式的選擇對址的確定和計算量至關(guān)重要，址應在保證 穩(wěn)定性的前提下由精度、單調(diào)性等綜合考慮確定，通常需要通過數(shù)值試驗來最終確
定。

１．２．２．２線性方程組的解法 線性方程組的解法有兩類：直接法和迭代法。直接法可通過有限步得到準確解， 求解精度高，但占有內(nèi)存大，僅適用于小型的方程組，而迭代法是從給定初值出發(fā)，
通過某種算法反復加以校正，直至收斂。迭代法主要用于求解大型多維問題的不對

稱稀疏線性方程組，缺點是不一定能保證迭代收斂，且終止迭代的準則不明確。
常用的迭代法是Ｊａｃｏｂｉ法和Ｇａｕｓｓ．Ｓｅｉｄｅｌ（ＧＳ）法，在ＧＳ方法中為加速收斂一般

均引入松弛因子田，稱為ＳＯＲ法（逐次超松弛法），若ｍ選取適當，可大大加速收 斂，∞的取值范圍為［０，２］，若∞＞ｌ，則稱超松弛法，ＳＯＲ法通常性能更好，因而
在計算流體力學中應用廣泛。ＳＯＲ法的一種變型為ＳＳＯＲ（對稱逐次超松弛法），

其收斂速度常遠遠超過ＳＯＲ，ＳＳＯＲ的每步迭代由兩小步組成，一是向前ＳＯＲ，一 是向后ＳＯＲ，實際計算中需通過數(shù)值試驗來選取甜。 １．２１３湍流模型技術(shù)
工程實際中存在的流動大部分都是湍流流動，湍流問題一直是計算流體力學中

的一個重要難題，湍流的特點是具有很強的隨機性、非恒定性，并且無一例外地都 是三維問題。由于湍流現(xiàn)象的復雜性，湍流運動以及與之相聯(lián)系的熱和物質(zhì)的輸運

現(xiàn)象都很難描述，從而也極難進行理論預測。如果要通過計算來解決這一問題，就
不可避免地對湍流輸運過程提出各種假設，采用一些經(jīng)驗性的結(jié)果和假定，把湍流 輸運過程中的各種物理量與時均流場聯(lián)系起來，這就是湍流模型研究的內(nèi)容。所謂

湍流模型就是以雷諾時均運動方程與脈動運動方程為基礎。根據(jù)現(xiàn)有的理論和經(jīng)驗 引進一系列假設而建立起來的一組描述湍流平均量的封閉方程組。湍流模型的選取

河海大學博士論文

不僅要考慮其模擬復雜流動的可靠度和精確度，而且還需要確定計算所需的費用和 處理問題所需的時間，后者對于工程應用尤為重要。通常評價湍流模型優(yōu)劣的標準 是：能夠適用于較多種類型的水流現(xiàn)象；具有足夠的模擬精度；人力和計算機費用 適度；復雜程度適當。顯然，那些適用范圍大，精度高的湍流模型必然復雜，消耗

費用較大，而好的湍流模型正是在這些相互矛盾的要求之間達到某種理想的平衡。 目前的湍流模型主要有雷諾時均模型、大渦模擬模型（Ｌａｒｇｅ Ｓｉｍｕｌａｔｉｏｎ－ＬＥＳｌ和直接數(shù)值模擬湍流模型（Ｄｉｒｅｃｔ
Ｔｕｒｂｕｌｅｎｃｅ—ＤＮＳｎ。
Ｎｕｍｅｒｉｃａｌ
Ｅｄｄｙ Ｓｉｍｕｌａｔｉｏｎ

雷諾時均模型：該模型在總體上可劃分為基于Ｂｏｕｓｓｉｎｅｓｑ假定的紊動粘性模型

和基于雷諾應力輸運方程的二階矩封閉模型，前者又可分為混合長度模型、一方程
模型和二方程（如ｋ—ｅ，ｋ—ｕ，ｋ－ｒ等）模型，后者則有代數(shù)應力模型（ＡＳＭ）和輸
運方程模型（ＤＳＭ）兩大類。

Ｂｏｕｓｓｉｎｅｓｑ假定的基本思想是：認為雷諾應力各向異性張量與平均場的應變率

張量成正比，即湍流應力的作用類似于層流應力，但其比例系數(shù)為一流動變量一渦
粘性系數(shù)，正是對渦粘性系數(shù)不同的計算方法導出了零方程、一方程、二方程等湍 流模型。不難想象，在渦粘性模型中，雷諾剪切應力的方向與平均場的應變率相同 這一特點可以說是渦粘性模型的主要缺陷之一，即雷諾應力與平均場應變率之間的

關(guān)系為線性關(guān)系，使得ｋ－ｅ模型無法再現(xiàn)流動中（如在方管中）出現(xiàn)的二次流“…， 這是由于復雜流動中，湍流場的變化總是滯后于平均場的變化。 在二階矩封閉模型中，雷諾應力與平均速度之間的關(guān)系隱含在一組聯(lián)立的偏微
分方程中，它們之間的相互關(guān)系包括一些復雜的物理機制，同時在雷諾應力方程中

的壓力一應變率相關(guān)項為雷諾應力輸運過程中提供了快慢兩種不同的變化機理，在 一定程度上反映了湍流的非線性作用，然而通常形式的二階矩封閉模型對雷諾應力 來說基本上是線性的，所包含的非線性作用非常小，還不能令人滿意地用來解決復
雜的工程問題。

上述各類雷諾時均湍流模型的特點及應用可參見文獻［１８，３３，３７，３８，３９，４０］，這 里就不再詳細展開。

大渦模擬模型（ＬＥＳ）：通過不斷研究人們認識到，在湍流流動中除存在許多隨機
性很強的小尺度渦運動外，還有一些組織很好的大尺度渦結(jié)構(gòu)，小尺度渦運動受流
動邊界條件的影響很小，且近似是各相同性的，流動中大部分質(zhì)量、動量或能量的

輸運主要來自于大渦運動。因此，如果將比網(wǎng)格尺度大的大渦運動通過數(shù)值求解 Ｎ．Ｓ方程來計算，而將比網(wǎng)格尺度小的小渦運動利用湍流模型來模擬，使得通過模 型提供的計算在整個求解中只占很小的比重，降低最終結(jié)果對模型的依賴性，從而 達到提高計算精度的目的，這就是大渦模擬的基本思想。目前大渦模擬方面的工作

河海大學博士論文

開展不多，主要集中在氣象學的流動問題或基礎性研究及簡單的流動問題，如文獻 ［４ｌ，４２】利用大渦模型分別對方柱繞流和明渠流動進行了數(shù)值模擬。關(guān)于ＬＥＳ更詳細
的介紹可參見文獻［３８］

直接數(shù)值模擬模型（ＤＮｓＴ）：模型化總是存在缺陷驅(qū)使人們尋求更好的解決
端流問題的途徑，由于Ｎ．Ｓ方程本身就是封閉的，不需要建立模型，由此提出了不 引入任何湍流模型，直接求解完整的三維非定常Ｎ．Ｓ方程，對湍流的瞬時運動進行 直接數(shù)值模擬，各種感興趣的物理量都可通過對瞬時量的平均運算來取得，這就是 直接數(shù)值模擬湍流。該方法顯然有許多優(yōu)點，如精確、可提供流場的全部信息，能

實現(xiàn)實驗室中很難做到的流動條件控制且比實驗更為經(jīng)濟，但是直接數(shù)值模擬的計 算工作量卻非常巨大，由于湍流運動所包含的單元（即小渦）比運動區(qū)域的尺度要
小得多，典型的數(shù)量級是流動區(qū)域尺度的ｌＯ。倍，而數(shù)值計算所需要的網(wǎng)格必須比

湍動單元的尺度更小，在三維情況下，為了覆蓋流動區(qū)域至少需要ｌＯ。個網(wǎng)格點，

如此多的網(wǎng)格點和相應的計算次數(shù)是現(xiàn)代計算機的儲存能力和運算速度所不能達 到的。目前國際上進行的湍流直接數(shù)值模擬還僅限于較低雷諾數(shù)（Ｒｅ＝２００）和簡單
幾何邊界條件的問題，主要用于湍流的基礎研究．如發(fā)現(xiàn)新結(jié)構(gòu)，揭示新機理，提 供新概念，檢驗和改進湍流模型等。 顯然，無論是大渦模擬還是直接數(shù)值模擬，主要的困難不僅在于計算枧的限制，

還有方法本身的問題，限于目前計算機的發(fā)展水平，短期內(nèi)還不能夠利用上述兩種方 法對工程問題中的復雜流動進行數(shù)值模擬，因此對現(xiàn)有的湍流模型加以改進，提出更
加合理的新的湍流模型依然是解決工程實踐中湍流問題的主要途徑。 重整化群（Ｒｅｎｏｒｍａｌｉｚａｔｉｏｎ Ｇｒｏｕｐ。ＲＮＧ）是一種用于構(gòu)筑許多物理現(xiàn)象模型的 通用方法，該方法起源于量子力學和高能物理中對基本粒子場的研究，自７０年代后

期開始把ＲＮＧ方法引入到湍流研究領域，Ｙａｋｈｏｔ和０ｒｓｚａｇ…。于１９８６年應用ＲＮＧ 方法建立了第一個湍流模型，簡稱為ＲＮＧｋ－ｅ。通過數(shù)值模擬表明它較之傳統(tǒng)湍流 建模方法具有顯著的優(yōu)越性和發(fā)展?jié)摿�。自從基于ＲＮＧ方法的湍流模型提出以后�?br />許多學者應用該模型進行了數(shù)值模擬，Ｓｐｅｚｉａｌｅ Ｇ．Ｇ．”“和Ｙａｋｈｏｔ““利用該模型求

解了后臺階繞流中的湍流分離現(xiàn)象，蔣莉““將ＲＮＧ ｋ—ｅ湍流模型與壁面率相結(jié)合， 應用于９０。彎曲槽道湍流流動的數(shù)值模擬，結(jié)果表明，該湍流模式可以有效地模擬有 曲率影響的湍流流動，李玲”“利用該模型對鈍體繞流的尾流場進行了數(shù)值模擬，結(jié) 果優(yōu)于標準的ｋ—ｅ湍流模型。由于利用ＲＮＧ方法建立的湍流模型中不包含任何經(jīng)驗 常數(shù)和可調(diào)節(jié)的參數(shù)，其模型參數(shù)是利用ＲＮＧ理論精確推導出來的，因而是通用的：
該模型適用于各種雷諾數(shù)范圍，包括層流、過渡流以及充分發(fā)展的湍流。此外，ＲＮＧ

模型能夠較好地反映流動的各向異性，對于與時間相關(guān)的大尺度運動也能給出真實的
模擬。因此ＲＮＧ湍流模型在湍流計算中具有進一步的推廣和應用的價值。

河海大學博士論文

１．２．４自由面模擬技術(shù) 自由液面問題是計算流體力學（ＣＦＤ）中非常重要的研究內(nèi)容，自由液面模擬的
合理性直接關(guān)系到含自由液面問題數(shù)值模擬成果的可靠性及精度。由于自由液面形狀 的復雜性和整體位置的不斷變化，特別是當自由液面發(fā)生破碎麗使得該闖題成為強非 線性問題，所以精確計算和模擬自由液面的難度非常大，可以說，自由液面問題已經(jīng)

成為計算流體力學領域的重要分支和前沿課題。目前在該領域的研究處于領先地位的 是日本的東京大學和日本船研所（ＳＲＩ），此外還有美國的ＬＡＷＯ大學和韓國的ＳＲＩ，，
我國的７０２研究所、上海交通大學、大連理工大學等單位對自由面模擬技術(shù)開展了一 些卓有成效的研究工作，但由于問題的復雜性，上述結(jié)果還遠遠沒有達到令人滿意的 程度ｃ４８－４９’５引。

下面對處理自由液面的一般方法及其應用作簡單描述，重點介紹ＶＯＦ方法。 剛蓋假定方法是將自由面看作是一個可移動的固體壁面，直接采用固體壁面的 “不可穿入”條件，僅適用大水體的宏觀運動； 靜壓假定法”“則認為流體中的壓強沿水深的分布服從靜壓規(guī)律，目前幾乎所有
的大體積水流數(shù)值模擬都以此假定為基礎；

高度函數(shù)法”“是將自由面處的水深ｈ定義為坐標（ｘ，ｙ，ｚ）和時間ｆ的函數(shù)，通過求 解高度函數(shù)ｈ滿足的微分方程來確定自由液面的位置，該方法適用于非恒定自由表面
問題．但無法處理標高函數(shù)為多值函數(shù)的情況，如射流、波浪破碎等問題。

流線迭代法”“是直接從Ｎａｖｉｅｒ－Ｓｔｏｋｅｓ方程出發(fā)，迭代求解壓力分布及流線位置。 文獻［２３１采用此方法計算了溢流壩反弧曲線段的自由水面。 除上述幾種方法外，目前處理自由液面問題的典型方法是歐拉模型中的ＭＡＣ法
和ＶＯＦ法。ＭＡＣ（Ｔｈｅ Ｍａｒｋｅｒ ａｎｄ Ｃｅｌｌ Ｍｅｔｈｏｄ）是設想在整個流場中均勻分布沒有

體積和質(zhì)量的微小顆粒一標記點，這些標記點以當?shù)亓黧w質(zhì)點的速度隨流體一起運 動，顯然只要計算出標記點在空間的位置即可確定自由表面的形狀和位置”‘…。該方 法的主要缺點是需要很大的計算機內(nèi)存和計算工作量，并且對于非均勻流場在網(wǎng)格中 會出現(xiàn)虛假的密度很高或很低的標記點，造成自由液面形狀的失真，而ＶＯＦ方法則 克服了上述不足，是目前研究自由液面問題的理想方法。
ＶＯＦ（Ｖｏｌｕｍｅ

ｏｆＦｌｕｉｄ）”“方法是定義這樣一個體積分數(shù)Ｆ，在有流體的點取值

為ｌ，在無流體的點取值為０，這樣，在一個網(wǎng)格單元內(nèi)，Ｆ的平均值則代表單元內(nèi) 流體所占的份額，若Ｆ的均值為１．則網(wǎng)格充滿流體，若Ｆ的均值為０，則為空網(wǎng)格，

若Ｆ的均值在Ｏ和ｌ之閻，則網(wǎng)格包含自由液面，為自由面網(wǎng)格。由于ＶＯＦ法追蹤
的是網(wǎng)格中的流體體積，而不是追蹤流體質(zhì)點的運動，因而具有容易實現(xiàn)、計算量小

和精度高等優(yōu)點，并且可以處理自由面折疊、自由面入水等強非線性悶題，成為目前
Ｏ

河海大學博士論文

自由面運動模擬的主流方法。

在ＶＯＦ方法中，關(guān)鍵技術(shù)是建立運動界面的輸運與重構(gòu)模型，如文獻［５７１研究了 運動界面追蹤技術(shù)中的界面重構(gòu)算法，并討論了由此引起的數(shù)值誤差問題，文獻［５８１ 以ＶＯＦ理論為基礎建立了三維空間自由面輸運及重構(gòu)的一階模型，文獻［５９１將有限體 積法和動網(wǎng)格技術(shù)應用到運動界面追蹤，自由面的位置和形狀利用運動邊界條件和空
間守恒率確定等，總之７０Ｆ方法中有許多問題有待進一步深入研究。

Ｉ＇３本論文的工作意義 本論文的主要工作是研究流場數(shù)值模擬中的關(guān)鍵技術(shù)，并將之應用于水利工程中 帶閘墩溢流壩過壩水流三維流場的水力計算。 ２１世紀人類面臨的主要問題是生態(tài)保護型可再生能源的安全供電、用于抵抗饑 荒、貧窮和疾病的優(yōu)質(zhì)而充足的供水以及可滿足人類基本需求的能源和食物，這些問 題的解決在很大程度上依靠設計開發(fā)新的水電工程來實現(xiàn)，水電工程對世界經(jīng)濟繁榮 和可持續(xù)發(fā)展起著至關(guān)重要的作用。在我國隨著水電技術(shù)的迅速發(fā)展，高壩的建設日 益增多，其特點是水流落差大，流量大，功率大，多位于河谷狹窄地帶，帶來的一個 突出問題是泄洪消能、防沖及安全．可以說我國在該領域遇到問題的技術(shù)難度居世界
第一，因此需要更加重視大壩安全問題。

在泄水建筑物一溢流壩的水力設計中需要解決的問題不僅包括壩面下游反弧段 上的摻氣、振動、空蝕防蝕及霧化等問題，同時需要解決壩面上游壓力分布、水面線 的位置及形狀、泄洪流量、閘墩對流量及水面線形狀的影響等問題。目前主要的研究 手段是物理模擬（水工模型試驗）和數(shù)值模擬．物理模擬有著公認的實際意義和科學 價值．它可以預演工程中非常復雜的三維水流現(xiàn)象，檢驗設計方案的可靠性和準確性， 但是它具有費用高、周期長、無法準確測取工程的某些湍流特征量和不便于進行多方 案優(yōu)選分析等弱點，因而現(xiàn)在物理模擬中已有相當部分的工作為數(shù)值模擬所替代。數(shù) 值模擬便于進行方案比較和優(yōu)選，同時研究高速水流、空化和空蝕現(xiàn)象已在相當程度 上可借助于湍流的數(shù)值模擬來實現(xiàn)。 在水利工程中，基于水深方向的二維模型在一定程度上能反映出壩面的基本流動 規(guī)律并解決了水利工程中的一些實際問題，如溢流壩反弧曲線的設計，因此在８０—９０ 年代得到了廣泛的研究和應用。許協(xié)慶利用有限元方法計算了泄洪洞進口水流和溢流 壩反弧水流…。；顧兆勛等利用Ｎ—ｓ方程及相應的湍流模型模擬了二維非定常帶自由 水面的湍流流動�！�；王奇峰采用極坐標系下的水流控制方程和ｋ－ｉ湍流模型模擬了 反弧水流的水力特性”１；魏文禮利用正交曲線坐標系和曲率修正的ｋ—ｅ模型模擬了 二維溢流壩流動，并提出了一種抗蝕性能較優(yōu)的反弧曲線形式“…。但是這些成果還 不能滿足水利工程實踐的需要，因為泄水工程中的水流具有明顯的空間三維特性，二

河海大學博士論文

維模型無法揭示復雜流場內(nèi)部的空間特性，必須借助于三維數(shù)學模型。目前國外在水
利工程中河口航道的三維水流數(shù)值模擬方面開展了較多的研究工作，提出了許多新模

型和新算法．國內(nèi)學者對河口三維水動力學的研究也給予了很大的關(guān)注……“，開發(fā) 出一些有特色的三維潮流場計算軟件。與此相比，在泄水工程的水動力學研究方面， 利用三維湍流模型全面揭示整個流場水力特性的研究工作開展不多，其中馬福喜利用 修正的ｋ－ｅ湍流模型和流體體積分數(shù)法模擬自由表面計算了重力壩及其下游消力池 中的水流流場和一個三維水躍流場…。，金忠青利用ｋ—ｅ模型數(shù)值模擬了帶有自由水 面的三維強紊動流場“…。對溢流壩壩型和平面布置（如：頂部溢流曲線、閘墩的形
狀、位置及尺寸等）的綜合三維數(shù)值模擬研究工作開展不足，而這些因素將直接關(guān)系

到工程的泄流能力、抗蝕體型的優(yōu)化乃至整個大壩的安全。 正是由于泄水工程的水力設計對水利大壩工程的泄洪效率及安全營運起著至關(guān)
重要的作用，因此必須進一步加強對泄水工程復雜三維流場特性的研究，其中包括建 立具有足夠精度的水動力學數(shù)學模型和提供高效、可靠的數(shù)值計算方法，在此基礎上 則可以在較多的設計方案和大壩體型范圍內(nèi)快速而準確地給出整個流場的定量結(jié)果，

從而把大型水利工程的建設建立在更扎實的科學研究的基礎之上。 １．４本論文的研究內(nèi)容 １）對網(wǎng)格生成技術(shù)進行深入研究，詳細給出了網(wǎng)格的代數(shù)生成方法：研究了邊界正 交化技術(shù)，給出了一種簡單的網(wǎng)格邊界正交化方法；建立了非對接分區(qū)結(jié)構(gòu)網(wǎng)格生
成方法，包括分區(qū)交界面上的數(shù)據(jù)結(jié)構(gòu)形式和數(shù)值計算方法：

２）利用有限體積方法建立了非正交同位網(wǎng)格系統(tǒng)中的ＳＩＭＰＬＥ算法，詳細給出了運 動方程的離散形式，數(shù)值計算采用二階精度的通量修正格式，并完成了傾斜方腔驅(qū) 動流、后臺階繞流、圓柱繞流等流場的數(shù)值模擬： ３）研究比較了各種湍流模型的特點及應用范圍，建立了ＲＮＧｋ一占湍流模型，并以水 翼粘流場為例完成了與標準ｋ一占模型的計算比較：

４）對運動界面模擬技術(shù)進行了比較深入的研究，針對ＶＯＦ方法中關(guān)鍵技術(shù)一自由面
重構(gòu)方法詳細給出了利用傾斜直線段模擬運動界面的方法，并對旋轉(zhuǎn)流場中的馬蹄
鐵形界面運動、平移流場中圓形界面的運動、剪切流場中圓形界面的運動進行模擬：

５）利用建立的數(shù)學模型對ＷＥＳ型溢流壩過壩水流二維流場進行數(shù)值計算，給出了三 種不同壩頂永頭下的自由水面線位置以及五種壩頂水頭下的壩面壓力分布曲線，并 與實驗結(jié)果進行了比較。同時計算了流量系數(shù)及壩面流速分布。
６）研究了帶閘墩溢流壩三維過壩水流粘流場，針對閘墩的型式及壩面的布置計算了

閘墩對水面線形狀、壩面壓力、流量系數(shù)等設計參數(shù)的影響，并將不同墩型與布置
形式的計算結(jié)果進行比較。

２

河海大學博士論文

第二章網(wǎng)格生成技術(shù)研究
在對流動問題進行數(shù)值計算之前，首先要解決如何將流動區(qū)域離散化的問題，目 前已經(jīng)發(fā)展出許多種將物理區(qū)域進行離散以生成計算網(wǎng)格的方法，這些方法統(tǒng)稱為網(wǎng)
格生成技術(shù)（ｇｒｉｄ
ｇｅｎｅｒａｔｉｏｎ

ｔｅｃｈｎｉｑｕｅｓ）。必須注意，網(wǎng)格生成技術(shù)在流場數(shù)值模擬過

程中占有非常重要的地位，網(wǎng)格生成質(zhì)量的好壞直接關(guān)系到計算結(jié)果的精度，因此需
要引起足夠的重視。

由于眾多的網(wǎng)格生成方法各有特點，在第一章中己經(jīng)介紹，關(guān)于一般的網(wǎng)格生成
方法可參閱文獻［６７，６８，６９，７０，７１，７２１，因此本章主要針對網(wǎng)格生成方法中的專項技術(shù)進

行討論，這些討論有利于在數(shù)值計算之前幫助建立合理的計算網(wǎng)格系統(tǒng)。 ２．１有限體積方法中處理不規(guī)則區(qū)域的網(wǎng)格生成方法 工程中大多數(shù)的流動問題都具有一個不規(guī)則的復雜幾何邊界，如水利工程中的溢
流壩、泄洪道、水輪機蝸殼等過水建筑物，船舶工程中的船體及附體，天然河道中的 水流等問題都包含了非常復雜的幾何邊界，為此需要研究針對不規(guī)則區(qū)域的網(wǎng)格生成
方法，常用的方法有：

２．１．１采用階梯形邊界逼近真實邊界 當流場邊界為曲線或斜直線時，一種最簡單的邊界逼 近方法是采用階梯形直線段逼近，如圖２．１，這種方法的
特點是同一網(wǎng)格線上的網(wǎng)格點（或控制體）的個數(shù)是不等
＼
，＿

真實邊界

ｌ

的，同時計算邊界為相互垂直的鋸齒狀粗糙表面，當網(wǎng)格
較稀疏時會帶來較大的計算誤差，隨著網(wǎng)格在曲線邊界處
的不斷細化，這種誤差會逐漸減小。
／ ／

≯
，

由于這種網(wǎng)格構(gòu)造簡單，可以適應任何形狀的邊界， 若同時配合網(wǎng)格的局部加密技術(shù)，該方法仍然不失為一種
處理不規(guī)則邊界的方法。

，

ｔ ｔ

計算邊界
圖２．１階梯形邊界逼近曲線邊界

需要注意的是，當采用這種方法逼近真實邊界時，如果把求解區(qū)域限制在計算邊
界之內(nèi)，則對每一個具體問題都需要單獨開發(fā)計算程序，如果把計算區(qū)域擴充為一個

規(guī)則區(qū)域，則通過特殊的處理可以針對不同的問題使用同一個規(guī)則區(qū)域內(nèi)的計算程 序，這是處理形狀不是特別復雜的計算區(qū)域的有效方法，稱之為區(qū)域擴充法（Ｄｏｍａｉｎ
ｅｘｔｅｎｓｉｏｎ

ｍｅｔｈｏｄ），具體細節(jié)可參見文獻［７３］。

河海大學博士論文

２．１．２采用曲線坐標變換 在水利工程中經(jīng)常采用的一種處理不規(guī)則邊界的方法是引入廣義曲線坐標進行 曲線坐標變換‘６８ｔ２３，７４，＇７５，７６】，即將原來的物理空間（ｘ，Ｙ，ｚ）變換到廣義曲線坐標空間 （孝，刁，ｆ），同時保證物理空間的幾何邊界與曲線坐標系中的坐標軸重合或平行，從而 將（ｘ，Ｙ，ｚ）空間中不規(guī)則的物理區(qū)域變換到（｛，７７，ｆ）空間中形狀規(guī)則的計算區(qū)域，如圖 ２．２所示，坐標變換公式為
孝＝毒（ｘ，Ｙ，ｚ）’
ｒ／＝ｒ／（ｘ，Ｙ，：），

ｆ＝ｆ（ｘ，Ｙ，２）

因此如果在變換后的計算平面上進行數(shù)值計算，則網(wǎng)格的生成就十分簡便。

物理空間

計算空間

圖２．２曲線坐標變換

但是，由于數(shù)值計算是在變換后的廣義曲線坐標系中進行的，因此描述流體運動 的控制方程也需要進行相應的變換，從而在方程中會出現(xiàn)一些附加項，使原本十分復
雜的控制方程變得更加復雜，同時這些附加項的物理意義不明確，而且在對附加項的 導數(shù)進行離散時增加了離散誤差，這是曲線坐標變換方法的主要不足之處。

２．１．３邊界貼體非正交網(wǎng)格 在計算復雜幾何區(qū)域的流體運動時，經(jīng)常使用邊界貼體非正交網(wǎng)格系統(tǒng)【３”，其中
網(wǎng)格形式可以是結(jié)構(gòu)網(wǎng)格、分區(qū)結(jié)構(gòu)網(wǎng)格，也可以是非結(jié)構(gòu)網(wǎng)格，這類網(wǎng)格的優(yōu)點是

可以適應任意的幾何形狀，同時比正交曲線坐標變換更容易實現(xiàn)。由于非正交網(wǎng)格的 網(wǎng)格線是沿邊界布置的，因而相對于階梯形逼近曲線邊界的情形更容易應用邊界條 件，特別是當使用分區(qū)結(jié)構(gòu)網(wǎng)格或非結(jié)構(gòu)網(wǎng)格時，還可以沿流動方向布置網(wǎng)格以及在 物理量變化劇烈的區(qū)域進行網(wǎng)格加密，因此目前在眾多的商用軟件中均采用邊界貼體
非正交網(wǎng)格系統(tǒng)。

當然，非正交網(wǎng)格也存在一些不足，如方程中由于網(wǎng)格非正交而產(chǎn)生的修正項增 加了方程求解的難度；非正交網(wǎng)格可能導致出現(xiàn)非物理解；此外變量在非正交網(wǎng)格系 統(tǒng)中的布置方式對求解精度和算法的有效性都將產(chǎn)生影響。

１４

河海大學博士論文

２．１．４分區(qū)結(jié)構(gòu)網(wǎng)格
結(jié)構(gòu)化網(wǎng)格是指網(wǎng)格系統(tǒng)中節(jié)點排列有序，相鄰節(jié)點的關(guān)系明確；分區(qū)結(jié)構(gòu)網(wǎng)格 是指把一個復雜的計算區(qū)域劃分成若干個塊，每～塊內(nèi)均采用結(jié)構(gòu)化網(wǎng)格，但不同塊

內(nèi)的網(wǎng)格系統(tǒng)可以不同。分區(qū)結(jié)構(gòu)網(wǎng)格的優(yōu)點是能夠用于幾何形狀復雜的流場計算， 同時可以使用結(jié)構(gòu)網(wǎng)格各個系統(tǒng)中成熟的求解技術(shù)，適合于并行計算。便于提高計算
效率［７７，７８ ７９１。

２．１．５非結(jié)構(gòu)兩格 與結(jié)構(gòu)化網(wǎng)格不同，在非結(jié)構(gòu)化網(wǎng)格中節(jié)點的位置無法用一個固定的法則進行有 序地定義，但是它可以適應任何形狀的幾何邊界，若在非結(jié)構(gòu)網(wǎng)格系統(tǒng)中應用有限差 分方法和有限體積方法，則可以使這兩種數(shù)值計算方法對不規(guī)則邊界的適應性增強到 與有限元方法等同的程度，這是非結(jié)構(gòu)網(wǎng)格的主要優(yōu)點１１７１。非結(jié)構(gòu)網(wǎng)格的缺點是網(wǎng)格
生成的工作量大，離散方程的求解速度較慢。
２．２

網(wǎng)格生成的代數(shù)方法一雙邊界方法 在眾多的網(wǎng)格生成方法中，代數(shù)方法以其簡單、通用性強而被廣泛使用，若同時

配合網(wǎng)格的區(qū)域加密技術(shù)、邊界正交化技術(shù)和分區(qū)網(wǎng)格技術(shù)，其適用的范圍和精度將

被大大提高，因此首先介紹一種具有通用意義的代數(shù)網(wǎng)格生成方法一雙邊界方法Ｉｔ，９，舯１。
設在物理平面上有一不規(guī)則區(qū)域ＡＢＣＢ，其中ＡＢ，ＣＤ為兩條不直接相聯(lián)的邊界， 首先選定這兩條邊界上的刁值，分別為ｒｉｂ，礬，于是這兩條邊界上的ｘ，Ｙ僅隨ｆ變化，
即

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

其中：６，ｔ分別表示底邊ＡＢ和頂邊ＣＤ。 雙邊界方法提供了～種在兩個邊界之問內(nèi)插的方法，從而確定內(nèi)部節(jié)點的位置， 具體公式如下：

ｘ（善，Ｊ７）＝工ｒ（ｆ）＾（７７）＋ｘｚ（亭）也（Ｊ７）＋�。�！墮ｊ；；三塑２ｈ，（ｒｉ）＋‘��！叢ｊ尹＾４（叩）
彤，野）＝Ｙｌ（嬲（椰＋Ｙｚｇ）＾：（７）＋魚生ｉ導盟縞（叩）＋魚塵粵業(yè)‰（彳）
ｄｒ／

ｄ１７

其中：

河海大學博士論文

（１）ＸＩ

Ｙ－，ｘ２，Ｙ２為邊界１和邊界２上的離散點坐標：

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

若將＾（ｆ＝１，２，３，４）表達式中的叩（２了乞擊）用一維伸縮函數(shù)

．

剛叩＋（Ｉ－Ｐ�。币�、ｔｈ（Ｑ（１ｒ－ｒ／））］

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

來代替，則可實現(xiàn)網(wǎng)格密度的控制，上式中的Ｐ，Ｑ即為控制節(jié)點分布的參數(shù)。

ｄＣ

ｄ，７

ｄ‘

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

Ｏ‘

系數(shù)Ｋｌ，Ｋ２用于控制正交網(wǎng)格深入?yún)^(qū)域內(nèi)部的程度。
說明：

（１）插值公式中的后兩項導數(shù)項是為了保證區(qū)域邊界處的正交性。

（２）當ｈ．＝２叩３—３叩２＋１，ｈ２＝－２ｒ／３＋３７／２時，可同時實現(xiàn)網(wǎng)格在兩個邊界處的加密
而當ｈ。＝１一，７，ｈ：＝ｒ／（線性情況）時，生成的網(wǎng)格比較均勻。 （３）當Ｐ＜１時，網(wǎng)格線向Ｊ＝ｌ的邊界處加密 當Ｐ＞Ｉ時，網(wǎng)格線向Ｊ＝ＪＭ的邊界處加密 （４）Ｑ用于控制善與刁的關(guān)系偏離線性的程度，對網(wǎng)格影響不大。

（５）若將紅表達式中的叩（２孟矗）用下式替代，也能夠?qū)崿F(xiàn)網(wǎng)格的加密

其中：成為大于ｌ的常數(shù)，若玎接近于１，則可使網(wǎng)格在叩＝０和玎＝ｌ（即Ｊ＝ｌ
和Ｊ＝ＪＭ）處加密。

２．３網(wǎng)格邊界的正交化技術(shù) 網(wǎng)格的正交性，特別是在物理區(qū)域邊界附近的正交性，對于整個流場的計算精度 是至關(guān)重要的。實現(xiàn)網(wǎng)格正交化的方法有許多種，比如：在微分方程網(wǎng)格生成方法中，
改善正交性的途徑是選取合適的控制函數(shù)，根據(jù)網(wǎng)格的正交條件ｇ．，＝０給出了控制

函數(shù)及控制方程，并實現(xiàn)了網(wǎng)格在邊界上的近似正交”４”。：在網(wǎng)格生成過程中引入
Ｈｅｒｍｉｔｅ插值，通過求解一個泊松方程和一個拉氏方程來實現(xiàn)網(wǎng)格的正交化”。。 在上述方法中，網(wǎng)格的正交性算法隱含在整個網(wǎng)格生成過程當中，并且與網(wǎng)格的

河海大學博士論文

生成是同時實現(xiàn)的，此外由于控制函數(shù)是按正交性的要求而設定的，故無法同時實現(xiàn)
網(wǎng)格的加密。還有一種正交化的方法是通過移動邊界點來實現(xiàn)的，但是當邊界曲線的 曲率較大時容易引起網(wǎng)格的變形“…。 文獻【６］介紹了一種網(wǎng)格邊界正交化的代數(shù)方法，本節(jié)以該方法為基礎，對其進

行了改進和完善，給出了更具有一般性，容易實現(xiàn)且獨立于網(wǎng)格生成過程之外的代數(shù) 網(wǎng)格邊界正交化方法，通過控制參數(shù)的設定很容易控制網(wǎng)格正交化的程度和范圍。 ２．３．１正交化算法
．、ｊ

本算法的基本思想是給定邊界點的分
一——

布，將區(qū)域內(nèi)部網(wǎng)格節(jié)點向?qū)吔琰c的法 矢量方向移動，從而實現(xiàn)網(wǎng)格在邊界區(qū)域附
近的正交化，下面以Ｊ＝１邊界為例說明該正 交性算法。如圖２．３，利用超限插值方法生
ｖ、ｍ

、、、ｊ．＾
。，．

一一。Ｊ－
（卜Ｉ，１）

，

、．蕊、～
，ｉ

成的區(qū)域內(nèi)初始網(wǎng)格點為ｘ：，Ｙ：，首先求
出邊界曲線Ｊ：ｌ在（ｉ，１）點處的法向矢量 月，則可以得到初始網(wǎng)格點在ｎ上的投形點
圖２．３網(wǎng)格在邊界處的正交化

《，Ｙ：。將初始網(wǎng)格點向投影點移動，則
可以實現(xiàn)網(wǎng)格在邊界處的正交性，移動的程度可通過一個權(quán)重數(shù)∞來描述。

ｘ。＝（１一％）ｘ：＋％ｘ； ＹⅣ：（１一曲Ⅳ）ｙ：＋Ｄ口Ｙ；
顯然，權(quán)重出，，應滿足０≤∞，，≤ｌ 甜。＝０表示網(wǎng)格點沒有移動 ∞。＝１表示網(wǎng)格點無阻尼地移動到對應邊界點的法線上。

在算法的具體實現(xiàn)中需要解決兩個問題： 第一個問題是網(wǎng)格點在對應邊界點法線上投影點的計算。通過解析幾何的方法推得投 影點的計算分式為：

ｘ：＝隊，砘）＋吉ｈ＋吲他＋％，） ｙ：＝一÷ｘ：＋ｙ。＋古ｘ。
其中Ｉ．為邊界線在點（ｉ，１）處的切線斜率

第二個問題是權(quán)重甜。的給定，本文采用文獻［６］推薦的公式，并進行改進，給出了更

河海丈學博士論文

一般形式的權(quán)重計算公式。

ｏａ，，＝ｅｘｐ｛枷卜老）－２＿ｌ】卜，烈‰酬＿。一ｐ】。
其中：ｄ，ｊ＝√（ｘｇ—ｚ叫）２十（ｊ，Ⅱ一ｙ，｛）２
Ｐ，Ｑ稱為阻尼參數(shù)，Ｐ值控制網(wǎng)格正交性在Ｊ方向上衰減程度，Ｐ值越大，則沿
Ｊ方向向區(qū)域內(nèi)部傳播的正交性衰減越快，Ｑ值控制網(wǎng)格正交性沿Ｉ方向的衰減程度， Ｑ值越大，則表明只有在邊界的中心區(qū)域網(wǎng)格是正交的。 在解決了上述兩個問題之后，即可實現(xiàn)網(wǎng)格的正交化。

２．３．２正交網(wǎng)格生成算例
為檢驗正交算法的可靠性，給出了正交網(wǎng)格生成的幾個算例…，為便于比較，在

圖２．４中給出了利用雙邊界法，但未采用正交化技術(shù)生成的初始網(wǎng)格。從圖２．４、２。５ 中可以看出，使用正交化技術(shù)能夠?qū)崿F(xiàn)網(wǎng)格在邊界處的完全正交，區(qū)域內(nèi)部正交化的 程度可根據(jù)情況進行調(diào)整，周時網(wǎng)搐的加密與正交化可同時進行。

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

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

審器至毒害耋薹輦ｊ尊毫享三量 �！鞭睫啡浪颍橹`：摹霉 蔓芝三三主三三妄；｝■÷如ｉ、‘‘：一．ｊ

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

～——

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

Ｏ瑚０．８０

４∞０．４０

～．一一
Ｏ∞

．

０∞

１∞

（１）

（２）

圖２．４網(wǎng)格生成算例１一未使用正交化技術(shù)，２－－ｆｌ�。瑁瘢榻换夹g(shù)

《三～
圖２．５正交網(wǎng)格生成算例
培

Ｉ∞

河海大學博士論文

２．４分區(qū)結(jié)構(gòu)網(wǎng)格技術(shù) ２．４．１分區(qū)結(jié)構(gòu)網(wǎng)格的一般特點
實際流動問題往往具有非常復雜的幾何物理區(qū)域，傳統(tǒng)的處理方法是通過曲線坐 標變換將復雜的物理區(qū)域轉(zhuǎn)化為規(guī)則的計算區(qū)域，由此導致變換后的方程組更加復

雜，方程中各項的物理含義也模糊不清，計算速度和精度大大降低，目前流行的處理 方法是采用非結(jié)構(gòu)網(wǎng)格技術(shù)和分區(qū)結(jié)構(gòu)網(wǎng)格技術(shù)，由于非結(jié)構(gòu)網(wǎng)搐的節(jié)點和單元分布 是任意的，因此具有非常優(yōu)越的幾何靈活性，可應用于任意形狀流場的模擬，其主要 缺點是需要更多的計算機內(nèi)存和ＣＰＵ時間，結(jié)構(gòu)網(wǎng)格中已經(jīng)成熟的算法在非結(jié)構(gòu)網(wǎng) 格中無法應用，同時在復雜粘性流場中仍存在一些尚需解決的問題。分區(qū)結(jié)構(gòu)網(wǎng)格的 基本思想則是把幾何外形復雜的計算區(qū)域劃分為若干個盡可能規(guī)則的子區(qū)域，在每個 子區(qū)域中獨立生成網(wǎng)格并進行流場計算，其最大的優(yōu)點是可以使用結(jié)構(gòu)網(wǎng)格中非常有 效、可靠的求解技術(shù)，適合于并行計算，計算效率高，因此分區(qū)網(wǎng)格技術(shù)越來越受到
人們的重視。 分區(qū)網(wǎng)格系統(tǒng)一般分為相鄰子區(qū)有重疊部分的重疊網(wǎng)格（ｏｖｅｒｌａｐｐｉｎｇ ｇｒｉｄｓ）１８
２＇馴

和相鄰子區(qū)無重疊網(wǎng)格的拼接網(wǎng)格【３２】（ｐａｔｃｈｅｄ ｇｒｉｄｓ），分區(qū)拼按網(wǎng)格系統(tǒng)又分為相鄰 予區(qū)網(wǎng)格線在交界面上對接和非對接兩種情況，非對接情況比對接情況要復雜，但其 應用的空間和優(yōu)勢也更大，因為可以實行逐塊加密而實現(xiàn)每一區(qū)的網(wǎng)格最優(yōu)，同時允 許沿交界面移動網(wǎng)格，因此交界面處網(wǎng)格線非對接的分區(qū)結(jié)構(gòu)網(wǎng)格擁有非結(jié)構(gòu)網(wǎng)格的 幾何靈活性和局部加密能力。同時可以利用可靠的結(jié)構(gòu)網(wǎng)格求解技術(shù)，使計算過程更
加簡單、高效．因而對于大多數(shù)流動問題都是非常有吸引力的。

２．４。２分區(qū)結(jié)構(gòu)網(wǎng)格及其數(shù)據(jù)結(jié)構(gòu) 分區(qū)結(jié)構(gòu)網(wǎng)格的優(yōu)點從圖２．６中可以看出，如果采用單區(qū)網(wǎng)格技術(shù)，則不論是利 用代數(shù)生成方法還是解橢圓方程生成方法，都將會產(chǎn)生網(wǎng)格線在物面處非正交、網(wǎng)格 扭曲、尖點多、網(wǎng)格走向與流動方向不一致等缺點，從而嚴重影響計算精度及效果，
若將流場按照其幾何外形和流動特點劃分為圖中所示的Ｉ，ｌＩ。ＩＩｌ區(qū)，則問題迎刃而
解。

為了保證在對分區(qū)結(jié)構(gòu)網(wǎng)格數(shù)值計算時能夠使用 單區(qū)結(jié)構(gòu)網(wǎng)格求解器，在分區(qū)結(jié)構(gòu)網(wǎng)格中需要采用兩 種數(shù)據(jù)結(jié)構(gòu)，即在每一區(qū)中的規(guī)則數(shù)據(jù)結(jié)構(gòu)（ｆ√）和
分區(qū)交界面的數(shù)據(jù)結(jié)構(gòu)，典型的分區(qū)交界面如圖２．７。
一

在規(guī)則數(shù)據(jù)結(jié)構(gòu)中，整個流場的網(wǎng)格節(jié)點和計算節(jié)點
按逐塊和逐列的方法利用一維數(shù)組來標記和存放。

ｊ
９

河海大學博士論文

，＝Ⅳ。ｋ�！ǎ妫保�；＋Ｊ

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

其中：ＮＢ～總的分區(qū)數(shù)；Ⅳ？，川ｋ一第ｋ?yún)^(qū)網(wǎng)格在Ｉ?Ｊ方向的網(wǎng)格節(jié)點數(shù)
Ⅳ：＝Ｎ？×Ｎ：一第ｋ?yún)^(qū)網(wǎng)格的節(jié)點總數(shù)。
此外，在交界面數(shù)據(jù)結(jié)構(gòu)中，需要對交界 面左右兩側(cè)控制體中心以及交界面上所有網(wǎng)格 節(jié)點進行序號標記。由于交界面上兩側(cè)子區(qū)的 網(wǎng)格線是非對稱的，故同一個控制體積在交界
面上可能擁有幾個控制面，如圖２．７中典型的
ｒ，ｔｒ冉 工 月 月 衛(wèi)

舡卉Ｊ

分區(qū)交界面，控制體積Ｌ在交界面上有３個控 制面．因此需要確認交界面上每一個網(wǎng)格面所 對應的左右兩側(cè)的控制體。關(guān)于分區(qū)結(jié)構(gòu)網(wǎng)格

的數(shù)值計算方法在第三章中單獨給出。
２．５小結(jié)

圖２．７分區(qū)網(wǎng)格交界面

本章以通用形式的代數(shù)網(wǎng)格生成方法為基礎研究了網(wǎng)格的邊界正交化技術(shù)，提出 了一種簡單的、能夠獨立于網(wǎng)格生成技術(shù)之外的邊界正交化方法，給出了詳細的數(shù)學

模型并編制了計算程序：通過網(wǎng)格生成算例表明，該方法能夠大大改善網(wǎng)格在邊界處 的正交性。此外對分區(qū)結(jié)構(gòu)網(wǎng)格技術(shù)進行深入的討論，并建立了非對接分區(qū)結(jié)構(gòu)網(wǎng)格
的數(shù)據(jù)結(jié)構(gòu)形式。

河海大學博士論文

第三章非正交同位網(wǎng)格系統(tǒng)中的ＳＩＭＰＬＥ算法
３．１交錯網(wǎng)格系統(tǒng)與同位網(wǎng)格系統(tǒng)
數(shù)值計算的基本思想就是將求解區(qū)域劃分成許多小的區(qū)域，即網(wǎng)格，在每一個網(wǎng) 格中定義求解變量并進行計算，網(wǎng)格中計算變量所在的點稱之為計算點，計算點的布 置方式有多種，可以在網(wǎng)格節(jié)點，可以在網(wǎng)格中心點，也可以在網(wǎng)格線中心點，在圖 ３．１中給出了三種情況計算節(jié)點對應的控制體。

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

ｊ－：Ｅ， ７∥，■

＿■一ｊ—ｊ．

１『一再，７

～Ｍ

７，≯謄ｎ，≯ 。；：Ａ ｊ夠

幽３．Ｉ計算節(jié)點及對應的控制體

（ａ）交錯網(wǎng)格系統(tǒng)

（ｂ）同位網(wǎng)格系統(tǒng)

圖３．２非正交網(wǎng)格系統(tǒng)中變量的布置形式

如果求解變量定義在不同的計算節(jié)點且擁有不同的控制體，則稱為交錯布置，對 應的網(wǎng)格系統(tǒng)稱為交錯網(wǎng)格（ｓｔａｇｇｅｒｅｄ ｇｒｉｄ），采用交錯網(wǎng)格系統(tǒng)的主要目的是解 決Ｎ—ｓ方程中壓力梯度項及連續(xù)方程的離散困難，同時避免出現(xiàn)波動的壓力場和速度 場１４”，但是在進行編程計算時，需進行各物理變量的定位及復雜的插值運算。 如果求解域中所有的物理變量均定義在同一個計算點上且擁有相同的控制體，則 稱為同位布置，對應的網(wǎng)格系統(tǒng)稱為同位網(wǎng)格（ｃｏｌｏｃａｔｅｄ ｇｒｉｄ）。典型的交錯網(wǎng)格
系統(tǒng)與同位網(wǎng)格系統(tǒng)見圖３．２。

顯然，只要能解決數(shù)值解的波動問題，同位網(wǎng)格在編程的簡易性和通用性方面更 具有優(yōu)勢，因為所有的變量擁有同一個控制體，它們對應離散方程中的許多項是～樣 的，計算量及存儲量均減少，但是由于網(wǎng)格面上沒有布置速度變量，在計算通過控制 面的通量時其速度需要由相鄰兩個計算節(jié)點上的速度插值得到，從而在程序中增加了

河海大學博士論文

更多的插值運算，盡管如此，同位網(wǎng)格系統(tǒng)在應用于復雜計算區(qū)域的非ｊ＿匕父嗍格和非

結(jié)構(gòu)網(wǎng)格中更具有交錯網(wǎng)格無法比擬的優(yōu)判２ ５１。
３．２

Ｎａｖｉｅｒ－Ｓｔｏｋｅｓ方程及其在非正交同位網(wǎng)格系統(tǒng)中的離散方法

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

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

…蘇．

（３．１）

ａｘ．

為便于理解，寫出分量形式

掣ｏｔ塒?嘞＝六一警＋Ｃ等＋等］ 掣塒?刪＝六一警＋（魯＋魯］
辦
ｌ
ｄ．
’

＠ｚ，

ｄ，Ｊ

＠ｓ，

掣塒?咖＝正一警＋［等＋等］
其中：Ｐ，＝一Ｐ＋２∥警一姜腳ｖ礦為粘性流體中的法向應力； 戚
ｊ

慨ａ，

。＝∥（等＋斟“一燁一Ⅲ粘性流體中的切向魷溈切向應
力所在面的法線方向，Ｊ表示應力本身的方向。

代入（３．１）式，并注意到研究對象是不可壓縮的牛頓流體：

掣＋掣＝∥一若＋毒卜考］
西＝一ｇｒａｄ（ｐｇｈ）。
其中：ｈ＝ｚ一‰，ｚｏ為參考點的垂向坐標。

ｃ。．ｓ，

其中：，為重力加速度ｇ在三個坐標軸方向上的分量。顯然，若ｐ＝ｃｏｎｓｔａｎｔ，則

則Ｐ＋觸＝蘆
從而可以將質(zhì)量力項合并到壓力項中處理，為方便起見，在以下各式中，將聲仍 記作Ｐ，從而可以得到微分形式的Ｎ?Ｓ方程：

河海大學博士論文

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

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

（３．６）

蘇。

將上述方程對控制體Ｑ積分，并利用Ｇａｕｓｓ定理即可得到積分形式的Ｎ．Ｓ方程：

曇￡艘地＋肛，礦．蕊＝一ｆＰｉ，?贏＋ｐ｝?蕊

（。。，）

其中：ｉ（ｆ－１，２，３）為三個坐標軸方向上的單位矢量，等式左端第一項為隨時間變化的
非定常項，第二項為對流項，等式右端第一項為包含質(zhì)量力的壓力項，第二項為粘性 力項，又稱擴散項。 積分形式的連續(xù)方程直接寫出：

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

（３．８）

需要說明的是，關(guān)于有限體積方法的特點在文獻ｆ３０，３３，５５】中有非常詳細的介紹，在 這里就不詳細展開，直接利用有限體積方法推導非正交同位網(wǎng)格系統(tǒng)中的離散模型。
３．２．２

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

３．２．２．１網(wǎng)格點的布置
本文采用非『Ｅ交同位網(wǎng)格系統(tǒng)，計算點布置在網(wǎng)格中心，對應的控制體即為網(wǎng)格
本身。

３．２．２．２動量方程的離散 以２Ｄ情況為例，推導非正交同位網(wǎng)格系統(tǒng)中動量方程離散形式，離散過程中涉
及到的計算節(jié)點、控制體、控制面等信息參見圖３．３
ｚ．
～

～～

ｆ～
釅
●Ｅ

Ｃ

，ｉ

．

一～～—鑫＆蘭
?

；
ｆ

ｓ

圖３．３計算節(jié)點與控制體

河海大學博士論文

１）對流通量Ｉｐｕ，Ｖ?ｒｉｄｓ
僅以控制面中的ｅ面為例說明，其它控制面的形式完全相同，速度分量“，１１－－Ｉ殳 符號≯替代。同時為保證離散形式具有二階精度，并避免產(chǎn)生波動解，采用中心格式
與上風格式聯(lián)合的方法，具體方法在３．４中介紹。

①顯式對流通量（中心差分格式）

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

（３?９）

卅，＝￡ｐ礦－贏＝ｂ礦?二ｋ ｓ。＝ｐｐ。?以＋ｓ７丸）為通過ｅ面的質(zhì)量通量；
ｓ。，Ｓｙ為控制面積Ｓ。沿Ｘ，Ｙ方向的投影面積： 以，≯，為速度在ｘ，Ｙ方向的分量：

以為控制面ｅ中心處ｅ點的速度，由線性插值確定，如圖３．４所示。

吮噸＋（覿ｏ。－－。）＋㈣（ｙ。Ｍ）

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

（３．１０）

其中，相鄰計算節(jié)點連線與ｅ面的交點ｅ處的坐標、遞廈及遠廈梯度ｉ電遼線住艷

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

丸咆＋罷∽－－。ｐ等帆誡．）
砘”辦（１－－五ｅ）＋罷（Ｘｅ－－Ｘｅ，）＋詈（壙虬）

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

ｃ。?Ⅲ

河海大學博士論文

其中：

丸＝Ｉ÷蘭ｉ為線性插值矢量；計算節(jié)點上的速度梯度計算在后面單獨給出。

圖３．４控制面上的變量插值

②隱式對流通量（上風差分格式）
ｆ

ｋ。尸＝上廁蕊：卜辦…
Ｉｍ。‘屯…

礦（ｙ?扦）。＞０

礦艫．’ｎｌ＜０
（３．１２）

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

則總的對流通量采用下式計算

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

（３．１３）

其中【】中的項稱為延遲修正項（ｄｅｆｆｅｒｅｄｃｏｒｒｅｃｔｉｏｎ），放在方程右端的源項之中， 即為離散方程中的Ｑ‘項，并利用上次的迭代值計算。采用延遲修正的原因是一方面 可使迭代數(shù)值解具有二階精度，另一方面可以避免因單純上風格式造成系數(shù)矩陣的非 對稱性，從而導致迭代算法的發(fā)散。

２）擴散通量［∥挈ｉ淼
�。掀�

仍以控制面ｅ為例說明，動力粘性系數(shù)∥及速度”，分別用一股變量ｒ和礦替代，
于是擴散通量的一般形式為：

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

①顯式擴散通量

河海大學博士論文

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

＝Ｌ伊＋－�。妫妫蟆�

（３ １４）

上式中需要計算物理量在網(wǎng)格面上的梯度值（—蘭）。，一般近似采用中心差分格式，即 佩
直接由控制體中心處的導數(shù)插值，只是需要對網(wǎng)格的非正交性影響而引入修正。

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

其中當網(wǎng)格在ｅ面處正交時，上式第二項等于０。

②隱式擴散通量

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

２等¨¨鄧。（刪∥（禹以） 一Ｃ¨ｓｏ（ａ彘＃ａｘｍ＋霧知”卜蹦刪刪“ａ隔ｒＥ－ｒｐ噸）（３
其中：ｆ，２ｒａｄｄ，）“有網(wǎng)格中一Ｉ）處的梯度線性插值，

１５）

禹噸。與掣
則擴散通量

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

扛乏面＝

Ｆ：

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

（３．１６）

上式［】中的項稱為延遲修正（ｄｅｆｆｅｒｅｄ

放在離散方程的源項之中，即

為∥項，并利用上一次的迭代值計算。
３）源項計算
首先給出壓力源項的計算

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

（３．１７）

河海大學博士論文

對于“勰酣＝一￡駟＝一（新△Ｑ
具甲：Ａｆ２為控制體體積，任意形狀控制體體積的計算方法見３．４．４ 對于ｖ方程，情形類似。
非定常項的計算方法如下：

（曇脅堿＝害＠？＋１－－４Ｕ～一）
＝４；“ｉ；１一Ｑ＝．

煦恥等卅，＝害㈣？叫）
此外，在對流項與擴散項的計算中，分別出現(xiàn)了對流項修正和擴散項修正，也作
為離散方程的源項進行處理，即

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

ｇ“＝ｋ哼…一ｋ“嚴
以“方程中的ｅ面為例給出具體的表達式：

餅＝川�！就锜o＋辦（１一ｔ）＋Ｘ。－－Ｘｅ，ＯＸ 一ｍｉｎ（ｍ。，Ｏ）屯＋ｍａｘ（ｍ。，Ｏ）紼

）＋警（咒一以）】 刪
。

Ｑ？２

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

最后給出由于網(wǎng)格非正交而引起的修正項

Ｑｇ＝呱刪Ｈ囂高詞】［１蠅】
３．２．３離散形式的Ｎａｖｉｅｒ－Ｓｔｏｋｅｓ方程
將（ｇ．１３）、（３．１６）式代入（３．７）式，得動量方程的離散形式：

４辦＋∑鉑＝ｇ
ｆ

（３．１８）

其中：Ｉ＝Ｅ，Ｗ，Ｎ，Ｓ為與計算點Ｐ相鄰的計算點

河海大學博士論文

。

喊

咿

鵬一ｋ

一 鵬一‰ 一 叢‰
娜
ｍ

一

卟

髓一ｋ

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

Ｑ一：鐘＋Ｑ；＋彰＋Ｑ；＋Ｑ８．該式右端分別為壓力項、對流項修正項、擴散
項修正項、非定常項以及非正交網(wǎng)格修正項對源項的貢獻。 ３．３非正交同位網(wǎng)格系統(tǒng)中的壓力校正方程及其離散
由于隱式離散方法對時間步長沒有嚴格的限制。故正如上節(jié)中介紹的那樣，動量

方程采用隱式方法進行離散，為推導壓力校正方程，我們將動量方程中的壓力項從源
項中分離出來，以ｕ方程為例進行說明。

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

（３．１９）

若動量方程的第ｍ次迭代值為“，，則

彳；“；’＋∑鐘“？＋＝（Ｑ。一餅）＋鐘
，

（３．２０）

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

群一

（３．２１）
‘

’

記右端第一項為印同時咖一１１駟－－（孰垃
蠆亥ＪＰ
△Ｑ，印、
（３．２２）

河海大學博士論文

顯然，速度“：’滿足動量方程，但不滿足連續(xù)方程，為此需要對該速度進行修正，
以滿足連續(xù)方程（３．８），在ＳＩＭＰＬＥ算法中，是通過修ｉＥ壓，Ｊ７場來實現(xiàn)的。

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

（３ ２３）

“；＝”；＋＋“；為滿足連續(xù)方程的速度，“ｊ為校正速度。
Ｐ”＝Ｐ＋Ｐ’為相應的壓力場，Ｐ’為校正壓強。

若將（３．２３）式代入連續(xù)方程塑竺衛(wèi)：ｏ，可得壓力泊松方程

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

（３．２ａ）

通過直接求解壓力泊松方程（３．２４），即可得到滿足連續(xù)方程的速度場Ⅳ？，但此 時的速度場及壓力場并不滿足動量方程，因此需要將迭代進行下去，上述的迭代過程，

即通過壓力場校正來滿足速度場，稱為外部迭代，而內(nèi)部迭代是指線性代數(shù)方程組的
內(nèi)部求解迭代，這種直接求解壓力泊松方程來確定壓強的方法稱為Ｐｒｏｊｅｃｔｉｏｎ
ｍｅｔｈｏｄ。

在實際計算中最常用的一種方法是壓力校正法，（３．２３）式減去（３．２２）式可得
到校正壓力Ｐ’與校正速度”’之間的關(guān)系

啦即籌（參
值之后控制面ｅ上的修正速度與修正壓強梯度依然保持上述關(guān)系，即

（３．２５）

為了應用連續(xù)方程（３．８），需要知道控制面上的速度修正和壓強修正，可以通過 插值的方法來確定。由于在修正速度的插值計算中包含了修正壓強梯度的插值，故插

咖∥一等《）。
∑Ａ，ｕｌ
中，該項忽略不計。
即

（３ ２６）

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

河海大學博士論文

同理

吩一等（答

（３ ２８）

將（３．２７）、（３．２８）式代入離散型連續(xù)方程ｍ。＋ｍ。＋ｍ，＋ｍ。＝０，并注意到由動量方 程第ｍ次迭代解出的速度場并不一定滿足連續(xù)方程，因此

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

其中△祝為迭代速度場對應的質(zhì)量通量誤差。

即（ｐｓｕ’）。一（ｐｓｕ’）。＋（ｐｓｕ’）。－（ｐｓｕ’），＋Ａｔｈ：＝０
在非正交網(wǎng)格系統(tǒng)中，修正質(zhì)量通量計算采用下面的公式：

（３．２９）

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

＝一器咖喊�。ǖ梗饛摹В� 鼽（拳ｚ南№ｈ卜（ｇｒａｄｐ’Ｈ瑤刊】
（ｇｒａｄｐｌ）。＝（ｇｒａｄｐ‘）Ｅ丑＋（ｇｒａｄｐ’）Ｐ（１～２０）

（３ ３０）

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

“（嬰）。丑。＋（孕），（１一五。）慨一ｙｐ） 們 卻
應當注意，式（３．３０）的第二項是由于網(wǎng)格的非正交性而引入的，當網(wǎng)格非正交 性較弱時，該項影響不大。但當網(wǎng)格是強非正交的，則必須計入該項的影響，將另外

三個控制面上的質(zhì)量通量代入（３．２９）式，并只進行一次修正，可得非正交網(wǎng)格中的 壓力校正方程為：

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

（３．３１）

舯班一（等）。而１
３０

河海大學博士論文

儼一（等）。瓦１麗
爿：＝＿（、ｐＡ雒ＤＳ）。百殺

班一ｃ等，，赤
爿ｆ＝一∑鐘 ａ；＝∑（ｇｒａｄｐ。）；（不一‘），ｋ＝ｅ，ｗ，ｓ，ｎ，Ｋ＝Ｅ，∥，Ｓ，Ｎ
△�。海揭唬ㄋ唬海埃海�，ｆｌ：＋嘲：）”１為前一次迭代速度場對應的控制體質(zhì)量通量。
至此，給出了非正交同位網(wǎng)格系統(tǒng)中ＳＩＭＰＬＥ算法中的動量方程和壓力校正方程 的離散形式，下面將ＳＩＭＰＬＥ算法的計算步驟簡述如下：

（１）將ｎ時刻的計算值“？，ｐ”作為初始值計算ｎ＋ｌ時刻的∥、ｐ∥；
（２）求解動量方程得到“？‘（該速度場不滿足連續(xù)方程，要通過修正壓力場來修正速
度場）

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

（４）修正速度及壓力，得到滿足連續(xù)方程的速度場和壓力場 “●＝“？＋＋“’，Ｐ”＝Ｐ＋Ｐ。 （５）重復（２）．利用妒和Ｐ”作為對∥”，Ｐ”’的修正，計算直到所有的修正量足夠
小。

（６）推進到下一個時間步長。
（２）～（５）稱為外部迭代。

３．４非正交同位網(wǎng)格系統(tǒng)中的專題討論 ３．４．１邊界條件 ３．４．１．１壁面邊界條件 壁面邊界條件通常采用非滑移邊界條件－（Ｄｉｒｉｃｈｌｅｔ條件），即流體速度等于壁面速 度，在有限體積方法中一般利用法向應力等于零的Ｎｅｕｍａｎｎ邊界條件，以二維情況
的南邊界為例進行說明，見圖３．５

河海大學博士論文

Ｎ

—７一—７７＿，一，—！，——≥—＿———７—７—７——７———，——．ｒ ／／／／／７／，／一。，／，／／７／ ／／／／／
圖３．５壁面邊界條件

“＝２∥（罷）。＝０ 砂
擴散通量的貢獻為

（３．３２）

由于壁面邊界是不可穿透的，故對對流通量項沒有貢獻，只有對擴散通量的貢獻。對

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

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

（３．３３）

其中Ａ，６，Ｓ。６分別為邊界條件對離散方程中系數(shù)及源項的貢獻。 對于３Ｄ情況，與２Ｄ情況相類似，下面以底面（Ｂｏｔｔｏｍ）為例．

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

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

＝∥�！辏üP＋警產(chǎn)。
Ｍｊ，摯，

Ｍ＆警

河海大學博士論文

２警曠－以Ｌ‘－－“ｓ
對于ｖ方程，類似可得：

（３ ３４）

ｆ．。＝警Ｖ，一警Ｖｓ
其中，底面對Ｗ方程中的擴散通量沒有貢獻

（３．３５）

３．４．１．２動量方程中的入１３及出１３邊界條件
由于入口邊界處的速度是已知的，故對流通量可直接計算

ｍ，＝ｐＶｎｓ，＝ｐ吣，＋＂ｖｓ，Ｊ
如果采用隱式離散格式，則入口邊界對對流通量的貢獻為

（３．３６）

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

（３．３７）

其中：第一項為入口邊界對ｕ方程源項的貢獻，ｍａｘ■，，ｏ）為入１３邊界對離散方程中系
數(shù)一．的貢獻。入口邊界對擴散通量的貢獻可類似于一般情形擴散通量的計算，只是 速度梯度采用單向插值計算：
Ｆｄ

ｆＦｇｒａｄｕ‘贏＝ｒ（黔，＝號¨訓

㈣。ｓ，

出口邊界條件對對流及擴散通量的影響，其處理方法類似于入口邊界，只是出口

邊界要放在距離入口較遠的下游，同時出口邊界上的速度采用一階上風近似，即 廬ｆ＝辦，≯＝“，Ｖ，‘ｗ ３．４．１．３壓力校正方程中的入口及出口邊界條件
壓力校正方程中的邊界條件同樣需要特別注意。若通過某邊界的質(zhì)量通量是給定 的，則壓力校正方程中的質(zhì)量通量修正等于０，即采用零梯度的Ｎｅｕｍａｎｎ條件．

在出口，若入口質(zhì)量通量是給定的，對于定常運動，則出口邊界上的速度可采用 外推法而獲得（即“，＝＂，），但是需要對由外推法得到的出口速度進行修正，以保證
各個流場內(nèi)的質(zhì)量守恒（即出口質(zhì)量通量與入口質(zhì)量通量相等），出口速度修正方法

見圖３．６。利用修正后的速度進行下一次的外部迭代，同時出口邊界上的質(zhì)量通量修
正值為零。

為保證解的唯一性，需要選取一個固定點的壓強作為參考壓強，所有其它點的壓 強修正計算都要減去同一個參考點的壓強，但是需要注意，在分區(qū)結(jié)構(gòu)網(wǎng)格中，參考

點只能取一個，不能每區(qū)各取一個。

河海大學博士論文

圖３．６出口邊界的速度修正

３．４．２延遲修正技術(shù)（ｄｅｆｆｅｒｅｄ ｃｏｒｒｅｃｔｉｏｎ）
在有限體積方法中．網(wǎng)格控制面上的變量值都需要通過插值來確定，若利用高階 插值的方法會導致計算量迅速增大，如對二維情況，利用四階辛普生（Ｓｉｍｐｓｏｎ）公式， 則每個通量需要計算１５個節(jié)點上的變量值，并且每個控制體對應的代數(shù)方程包含２５ 個節(jié)點上的變量值，對于非線性方程組來說，其求解工作量是非常大的。

我們知道，利用前一步的迭代值顯式計算高階通量是簡單的，而隱式的低階通量 僅需要與計算節(jié)點最鄰近的節(jié)點上的變量值，計算工作量小。若能將二者結(jié)合起來，
則將同時兼顧計算的精度與簡單性，這種思想由Ｋｈｏｓｌａ
出【８卯，具體形式如下：
ａｎｄ

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

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

ｌｏｌｄ

ｔ３．３９）

其中：Ｆ州為某種低階隱式離散格式； ｒ州為某種高階顯式離散格式； 式中第二項由前一步的迭代值計算，與第一項的隱式部分相比是很小的，因而顯 式處理不會影響計算的收斂性。這里將高階與低階離散格式聯(lián)合的目、的有兩個：一是
避免出現(xiàn)物理上不合理的波動解：二是改善某一迭代求解算法的收斂性［３２，３３１。

河海大學博士論文

３．４．３控制體中心的梯度計算
在對流通量，擴散通量及壓力源項計算中 計算方法如下，以擴散通量為例：
Ｆｅｄ

都需要知道控制體中心Ｐ處的梯度

ｔ莩（黔。Ｊ

（３．４０）

其中交界面ｅ處的梯度可以通過控制體積中心點上的梯度插值得到

（考］。＝［考］，屯＋［等］，ｏ以，
控制體積中心處梯度可利用Ｇａｕｓｓ定理來確定

ｃ。．ｔ?，

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

ｚ莩箬２面１融＋丸Ｓ。＇－＃ｉｗ吖一∥）
３．４．４三維任意形狀網(wǎng)格控制體體積計算

（３．ｔ２）

任意形狀的控制體積計算的原理是Ｇａｕｓｓ定理， 即將體積分轉(zhuǎn)換成面積分
由于

擊ｖ（ｘｉ）＝１，故控制體體積

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

（３．４３）

其中ｃ為控制體擁有的控制面，疋。為控制體的面積矢量在ｘ方向上的分量，即

控制面積Ｓ。在ｙｚ坐標面上的投影面積

ｓ。＝ｓ。。ｉ＋ｓ。ｙｊ＋墨。ｉ
下面以空間六面體為例，ｃ＝ｌ、２、３、４、５、６，工�？扇。憧刂�

面上四個頂點Ｘ坐標的平均值計算，如圖３．７，Ｓ。１為第ｃ個控
制面在ｙｚ坐標面上的投影面積
圖３．７控制面投影

ｓ。Ｊ＝圭Ｉ舅一ｉ）×Ｅ—ｉＩ＝吉［（ｙ．一ｙ，）（Ｚ２－－Ｚ４）一（２＇Ｉ－－Ｚ３）（ｙ：一ｙ。）】
應為負號。

（３．４４）

將每個控制面上的Ｘｃ和Ｓ�！嬎愠鰜碇�，利用（３．４３）式即可計算出該控制體體積。 注意：若控制面外法線方向與Ｘ坐標軸方向的夾角大于９０。，則（３．４３）式中該項前

河海大學博士論文

３．５分區(qū)結(jié)構(gòu)網(wǎng)格中ＳＩＭＰＬＥ算法的實現(xiàn)

動量方程的離散形式為：Ａｐ辦＋∑Ａｋ≯ｋ＝Ｑ尸，女＝Ｅ，∥，Ｎ，Ｓ，在單區(qū)網(wǎng)格內(nèi)部，
ｌ

每一個控制體擁有４個控制面，而交界面處的控制體則可能擁有多于４個的控制面， 故實際計算時每個網(wǎng)格面上的對流通量和擴散通量需單獨處理，然后補充到交界面兩 側(cè)控制體的離散方程系數(shù)中去，如計算節(jié)點上對應的ｅ面為分界面，則�。�，＝０，在 離散方程的系數(shù)Ａ，中補充交界面的貢獻一Ｒ，Ａ。，，Ａ馬，同理，若控制體Ｒ：的∥面

為分界面，�。粒剑�，在對應的系數(shù)Ａ，中補充交界面的貢獻Ａ。。 在計算對流通量時，需要知道物理量在控制面上的值，對于正交網(wǎng)格可直接有網(wǎng)
格面兩側(cè)的計算點插值得到。但在一般情況下，交界面上的網(wǎng)格線是非正交的，如圖 ３．８所示，即，與，’不重合，因此交界面上每一個網(wǎng)格面的變量不能由兩側(cè)計算節(jié)點直 接插值得到，而需要進行修正，即

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

?Ｎ

／
，

?

其中：

，

一一’一７Ｌ——’－

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

“？≮～
圖３．８

一…ｉ。————；ｉ

學），＝（缸”嘗Ｍ，１）
以上兩式均由網(wǎng)格面兩側(cè)計算節(jié)點線 性插值確定。
非對接分區(qū)交界面

Ａ，＝；｝—魯為線性插值矢量，
‘

ｈ一■ｌ

計算節(jié)點的梯度值采用下式計算

。耠挲
靠
ｙ

其中Ｓ：包括該控制體在分區(qū)交界面上的所有網(wǎng)格面。

上述格式可以保證交界面上的通量計算具有二階精度。 在計算交界面上擴散通量的貢獻時，需要知道交界面上的梯度矢量，若采用對流 通量的處理格式，則將出現(xiàn)高階偏導數(shù)從而增加了計算工作量。常用的方法是直接采

河海大學博士論文

用中心差分格式，即由控制體中心的導數(shù)線性插值，同時對交界面處網(wǎng)格的非正交性 進行修正。

ｃ等�。憧迹唬察o焉一ｃ等，產(chǎn)ｃ音�！�，，
其中當網(wǎng)格在交界面處正交時右端的第二項等于０。

非正交網(wǎng)格中的壓力校正方程的離散形式為：爿；ｐ：＋∑彳，ｐ，，ｉ＝一△�。唬r
其中需要注意的是對交界面上各個網(wǎng)格面上的質(zhì)量通量進行修正，以，控制面為 例：

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

一器（扣小礎如ｄｐ＇Ｍ¨）］
其中：（！蘭），ｚ＝＿�。剑椋郏ǎ穑阂唬穑海┮唬ǎ纾颍幔洌穑�，（乓一五）】 靠（ｒＲ一丘妒
（ｇｒａｄｐ’），＝（ｇｒａｄｐ’）Ｒ＾＋（ｇｒａｃｌｐ’）￡（１—４）

（ｇｒａｄｐ‘），．（不一五）＝【冬）。丑＋（冬）。（１一刪‰一吒） 盤 劣 ＋【冬）。＾＋（嬰）。（１－丑）艦一ｙＬ） 卯
ｃｒｙ

應當注意，上式的第二項是由于網(wǎng)格的非正交性而引入的，當網(wǎng)格非正交性較弱時， 該項影響不大，但當網(wǎng)格是強非正交的，則必須計入該項的影響。 ３．６非正交同位網(wǎng)格系統(tǒng)中ＳＩＭＰＬＥ算法求解算例 ３．６．１傾斜方腔頂板驅(qū)動流 頂板移動的方腔內(nèi)流場計算常被用來檢驗數(shù)值算法的穩(wěn)定性和精度，為了便于比 較，本文首先計算了正交網(wǎng)格下方腔的內(nèi)流場，網(wǎng)格數(shù)為８０×８０，迭代過程中速度 及壓強的欠松弛矢量為０，８，０．８，０．３，給出了速度分布曲線，并與文獻［３３】的結(jié)果進 行比較，非常一致，見圖３．９。為了檢驗非正交網(wǎng)格系統(tǒng)下的ＳＩＭＰＬＥ算法，對方腔 壁傾斜角為４５。時的傾斜方腔內(nèi)流場進行了數(shù)值模擬，網(wǎng)格數(shù)仍為８０×８０，見圖

河海大學博士論文

３．１０，但迭代過程中速度及壓強的欠松弛矢量為０．７。０．７，０．２。在圖３．１ｌ中給出了雷

諾數(shù)Ｒｅ＝１００時的斜方腔中縱剖面及中橫剖面上的速度分布曲線，其中Ｒｅ：墮，“為
，，

頂板移動速度，，為方腔長度，ｙ為流體的運動粘性系數(shù)，圖３．１２，３．１３中分別給出 了Ｒｅ＝１００和Ｒｅ＝１０００時對應的流線矢量圖和壓力等值線分布圖，可以看出，在斜方

腔右下方產(chǎn)生的次渦中心隨著雷諾數(shù)的增加向上移動。
’∞］
：—４

ｙ州
ｏ∞一

，．一一一
彳

…ｊｔ忡ｍＩ

…
？

五
『『－
／．

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

¨、ｊ

７

。。ｊ／
Ｄ∞ ａ∞
ａ∞

了
１∞

＼
０∞ 口棚ｏ∞
ｔＭ

圖３．９

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

１∞

８－－９０。，Ｒｅ＝１０００，方腔中縱斷面及中橫斷面上的速度分布曲線： 圖中離散點為文獻【２５】中的計算結(jié)果

圖３．１０傾斜方腔中的非正交網(wǎng)格

河海大學博士論文

一一㈠
圖３
，

一

ｌ １

Ｂ＝４５０，Ｒｅ＝１００，方腔中縱斷面及中橫斷面上的速度分布曲線

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

２０

●∞

¨ ∞ ¨ ∞

ｏ ０ １ １ ５ ２

（ａ）

Ｒｅ＝１００

，

¨ ％ ¨ 艙

ｏ ０ ２

（ｂ）

Ｒｅ＝１０００

圖３．１２ １３＝４５。傾斜方腔內(nèi)的流線矢量曲線

河海大學博士論文

圖３ １３

Ｂ－＝４５。傾斜方腔內(nèi)的流線矢量曲線

３．６．２后臺階繞流
后臺階繞流是典型的流動問題，如文獻［８６，８７］分別對三維后臺階的層流流動及后

臺階繞流的出流邊界條件進行討論。為檢驗分區(qū)結(jié)構(gòu)網(wǎng)格算法的可靠性，對計算流體 力學中的典型算例一后臺階繞流進行分區(qū)計算，臺階的尺寸為入口高度５．２ｍｍ，臺階 高度４．９ｒａｍ，臺階上游長度２００ｒａｍ，臺階下游長度５００ｒａｍ，計算結(jié)果給出了不同Ｒｅ下

４０

河海大學博士論文

的回流區(qū)長度的變化曲線，其中Ｒｅ中的特征速度和特征長度分別為臺階上游的平均 速度及臺階高度，同時給出了Ａｒｍａｌｙ
ｅｔ

ａ１．的試驗曲線【７８】和文獻【７９】的計算結(jié)果見圖

３．１４，可以看出，在Ｒｅ＜４００時，計算結(jié)果的一致性很好，而當Ｒｅ＞４００時與試驗曲 線偏離較大，其主要原因是試驗數(shù)據(jù)是在三維流動的情況下獲得的。

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

２

３

４

Ｒｅ

５

６

７

８

圖３．１４二二維后臺階回流區(qū)長度計算結(jié)果比較

３．６．３非對稱平板間圓柱繞流 當研究位于兩個平板間的圓柱繞流時【Ｅ”，如圖３．１５，重點關(guān)注的是柱面壓力與尾 渦，因而需要對柱面附近和尾部區(qū)域使用較密的網(wǎng)格，此時使用單區(qū)網(wǎng)格是不合適的，

要么整體利用非結(jié)構(gòu)網(wǎng)格，或者使用分區(qū)結(jié)構(gòu)網(wǎng)格．本文采用六分區(qū)網(wǎng)格對其進行了
數(shù)值計算，計算條件為：圓柱直徑為Ｄ，圓柱中心至上下平板的距離分別為２．０８Ｄ和 ２，ＯＤ，圓柱上游距離為８Ｄ，下游距離為２０Ｄ，來流速度為ｕ，雷諾數(shù)Ｒｅ＝ＵＤ／，，， 這里僅給出Ｒｅ＝１００和Ｒｅ＝１０００對應的壓力等值線，速度等值線和流線圖，見圖３．１６。 從圖中可以清楚地看到圓柱后尾渦的形成和脫落過程，隨著Ｒｅ的增加，流線的波動 逐漸加劇。

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

一

一

圖３ １５平板間圓柱繞流

４Ｉ

河海大學博士論文

圈３ １６

Ｒｅ＝１００（ａ）和Ｒｅ＝１０００（ｂ）對應的壓力，速度等值線與流線圖

３．７小結(jié) 本章在非正交同位網(wǎng)格系統(tǒng)中推導出Ｎａｖｉｅｒ－Ｓｔｏｋｅｓ方程的離散模型、壓力校正方 程及其離散模型，對非正交同位網(wǎng)格系統(tǒng)中的邊界條件、延遲修正技術(shù)、控制體中心 的梯度以及任意形狀控制體體積計算等進行專項討論，建立了非正交同位網(wǎng)格系統(tǒng)中 的ＳＩＭＰＬＥ算法，并以此為基礎建立了非對接分區(qū)結(jié)構(gòu)網(wǎng)格數(shù)值求解模型。利用上述 算法編制計算程序，對傾斜方腔驅(qū)動流、后臺階繞流、非對稱平板間圓柱繞流進行計 算，通過與實驗結(jié)果的比較表明，本文建立的數(shù)學模型是可靠的和正確的。

河海大學博士論文

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

湍流模型

大多數(shù)工程中的流體力學問題都是湍流問題，應當用湍流理論進行處理，但是由

于實際工程中的流動甚為復雜，目前的湍流理論和計算技術(shù)都未達到可以大量解決實
際工程中流體力學問題的水平，人們?yōu)榱藵M足工程應用與設計方面的要求，往往顧及

不上湍流場的分辨率，而僅對湍流影響作簡單化處理。有些涉及到流體內(nèi)部結(jié)構(gòu)的工 程流體力學問題，不用湍流理論不可能較好地得到解決，如環(huán)境工程、水利工程中的
高速水流和挾沙水流、空氣動力學及船舶動力學中的繞流問題等都是如此，因而近

２０年來，國內(nèi)外都在研究如何用湍流理論解決工程中的湍流問題。其主要內(nèi)容包括： ①利用目前比較成熟的湍流理論解釋一些特定的流動現(xiàn)象， ②應用不同的湍流模型預測工程中的流體力學問題。， ③以湍流理論為依據(jù)，提出一些解決工程問題的新理論和新方法。 上述研究內(nèi)容緊密聯(lián)系，相互促進，一方面湍流理論已被廣泛地用于解決工程問
題，同時在工程應用中又促進了湍流理論本身的發(fā)展，如直接服務于工程的湍流模式
理論等。

根據(jù)目前湍流研究的趨勢，除對湍流基本結(jié)構(gòu)大力進行研究外，結(jié)合工程中的問 題進行工程湍流研究，己成為國際湍流研究的重要方面，如湍流模式、湍流控制、湍
流測量、湍流與傳質(zhì)傳熱等。

因此，以工程為背景開展湍流研究，對于國民經(jīng)濟建設和科學技術(shù)的發(fā)展都具有 十分重要的意義，為此，必須以發(fā)展和利用已有的湍流知識解決實際技術(shù)問題為目標， 在工程中廣泛應用湍流理論，為解決工程湍流問題提供實踐基礎，同時積極促進湍流
理論及其相關(guān)技術(shù)的發(fā)展，正如周培元教授指出的那樣需要“積極開展流體湍流運動 的實驗、理論和應用研究”１３ ７１
ｏ

４．２湍流模式理論 由于湍流的復雜性，目前要想獲得復雜的計算結(jié)果幾乎是不太可能的，需要引入 各種假定，這些假定是對實際湍流流動的某種近似模擬。所謂湍流模式就是指將真實 的湍流運動�；癁槿藶樵O計的模型，這種模型的平均行為應與實際湍流統(tǒng)計平均行為
基本一致。
。

雷諾時均模型在總體上可劃分為基于Ｂｏｕｓｓｉｎｃｓｑ假定的湍動粘性模型和基于 雷諾應力輸運方程的二階矩封閉模型，前者又可分為混合長度模型、一方程模型和二 方程（如ｋ．￡，ｋ．ｍ，ｋ－ｒ等）模型，后者則有代數(shù)應力模型（ＡＳＭ）和輸運方程模型（ＤＳＭ）
兩大類。

河海大學博士論文

Ｂｏｔｔｓｓｉｎｅｓｑ假設的基本思想在于，認為湍流應力的作用類似于層流應力

ｆ＝∥罷，但其比例系數(shù)為～流動變量——渦粘系數(shù)以，即ｆ；∥．—ｄ＿ｕ。其中＾是
ｕｙ

８ｙ

與，完全不同的量，它不是流體的一種性質(zhì)，而是決定于湍流的流動狀態(tài)，因而幾在 不同流動狀態(tài)，甚至在同一流場中的不同點，其值都很不相同。根據(jù)對渦粘系數(shù)一的 近似計算方法的不同，即可構(gòu)造出不同的湍流模型，如零方程，一方程、二方程等湍 流模式。需要注意，在渦粘性模式中，雷諾剪切應力的方向與平均場的應變率相同， 即渦粘系數(shù)是各向同性的標量，這一特點是Ｂｏｕｓｓｉｎｅｓｑ假設模式的主要缺陷之一。 與Ｂｏｕｓｓｉｎｅｓｑ假設模式不同，在二階矩封閉模式中，雷諾應力與平均速度之間的 關(guān)系隱含在一組聯(lián)立的偏微分方程中，在一定程度上反映了湍流的非線性作用。然而 通常形式的二階矩封閉模式對雷諾應力來說基本上是線性的，所包含的非線性作用非 常小，還不能令人滿意地用來解決復雜工程問題，可以說非線性湍流模式在理論和應
用上目前都已成為工程湍流計算研究中的核，Ｉｉ，問題之一【８９】。

４．３湍流運動的基本方程 在工程計算中，往往是流動參數(shù)的平均量對工程設計更具有重要作用。因此，流 體的運動方程都建立在雷諾平均的概念下，即對Ｎ．Ｓ方程取某種形式的平均，如時間
平均，空間平均、系統(tǒng)綜合平均（指對相同初始條件和邊界條件下大量實驗結(jié)果求平
均）。

對不可壓流體來說，平均場的動量方程（即雷諾方程ＲＡＮＳ＝Ｒｅｙｎｏｌｄｓ
Ｎａｖｉｅｒ－Ｓｔｏｋｅｓ

Ａｖｅｒａｇｅｄ

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

魯五籌一羆＋毒（，謄ｄｘ一瓦）壇 ａ。教ｋ礎：融｜、．�！� —二＋“Ｉ一一十一ｌｙ一一“．”｝＋Ｐ．

其中：“∽為雷諾應力，是方程中的未知量，Ｕ．和ＪＰ分別為平均速度分量和壓力， 島為體積力，，，為流體的動力粘度。當平均方式為時間平均時，－拋＇－；７。：０。如果流動 是非定常的，則不能取是平均。
不可壓縮流體運動的連續(xù)方程為：

墮：０
敏．

（４．２）

由于在方程中出現(xiàn)了雷諾應力，因而方程（４．Ｉ）、（４．２）是不封閉的，即未知量

河海大學博士論文

的個數(shù)大于方程的個數(shù)，因此求解控制方程（４．１）、（４．２）的首要問題是確定雷諾應

力（一“，“，），在湍流模式理論中是通過引進低階的關(guān)系或平均量來近似表示這些湍動 量的，由這些微分形式或代數(shù)關(guān)系表示的湍動量模擬關(guān)系式與控制方程（４．１）、（４．２）
共同組成封閉方程組。

需要注意的是，雷諾應力與粘性應力是不同的，粘性應力是流體分子運動與分子 相互作用的效應，而雷諾應力則是流體宏觀的湍動效應，不是流體的物質(zhì)本性，是流

體運動到一定階段后的產(chǎn)物，從力的觀點看，雷諾應力是遷移慣性力（非線性項）因 取時均而出現(xiàn)的，因此在實質(zhì)上是屬于慣性力的范疇，只有在統(tǒng)計意義下雷諾應力才
會出現(xiàn)。 在雷諾平均框架下的流動控制方程中，湍流的全部信息可以被認為是包含在二階 速度相關(guān)項“。“，之中，若“，“，可以被忽略不計時，則方程（１）回歸到層流的控制方 程。由此可見，層流可被看作是湍流的一個極限特例，在大部分工程湍流問題中，甜．“， 不僅不能忽視，而且它對流動混合，擴散等起決定性作用。 在建立“弘的封閉模式之前，有必要討論雷諾應力輸運的控制方程，該方程的具 體推導過程可參見文獻［９０１，這里僅給出表達式：

警：導（，學一而一叢小絲西）
一—－＿ 一—－＿

＿—＿．

＿２ｖ當當 面 生％ 麗 墮奶 卜 ‰ｐ 贏一％ 巫魄
＋ ＋ ０冀．Ｏ％．

（４．３）

該方程等式左邊為應力的對流輸運項（ｑ），右邊的第一個括號中的項（吒）代表
湍流的擴散特性，控制湍流的空間分布，由三個不同的物理機理組成：粘性擴散、三 階速度相關(guān)和壓力與速度的相互作用。第二個括號表示湍流的生成機制（只），即湍
流在通常情況下由雷諾應力和速度梯度的相互作用所產(chǎn)生。第三個括號為湍動壓力與

應變率的相關(guān)項，最后一項為應力的耗散率（毛）。為了方便起見，（４．３）可寫作

ｑ＝吒＋弓＋ｍ，一毛
通過對方程中的毛，中。和６１ｉ作適當�；纯傻玫讲煌耐牧髂Ｊ健�
４．４

（４?４）

ｔ一￡兩方程模型。
對方程（４．４）中的ｆ√進行縮并，根據(jù)縮并原理及連續(xù)條件有中。＝０，可以得到

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

河海人學博士論文

嬖：§ｙ婺一一ｋ＇ｕｔ一土麗卜瓦婺一ｙ挈晏 面２面咿面一 一石掣ｐ“１，嘶葛吖菇瓦
項為湍動擴散，右端第二項為湍動能發(fā)生項，最后一項為湍動能耗散項。

（４．５）

其中：左端為動能變化率，右端第一項為擴散項，【】中內(nèi)第一項為分子擴散，后兩

在ｋ一￡模型中，為了封閉湍流運動方程，還需要增加耗散率占的控制方程：

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

：隔２
（４．６）

式中宰為耗散能變化率，右端第一項為擴散（包括紊流擴散和分子擴散）項，第二
及第四項為產(chǎn)生能項，第三及第五項為破毀能項

方程（４．１）一（４．６）共１２個方程，若取Ｂｉ，Ｍ：“：，ｋ，ｓ共１２個量作為未知量?
則構(gòu)成一個近似封閉的方程組，為了能夠求解，還需要將方程組中其余的未知量與 １２個選定的量建立聯(lián)系，因此還需要對方程（４．４）、（４．５）、（４．６）進行摸化，給出

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

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

一“

可

面ｉ

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

韓只＝＿（而蕁＋麗習～再等。

河海大學博士論文

Ｃｋ＝０．２２５，Ｃｌ＝２．２，Ｃ２＝ｏ．４，ｅ＝ｏ．１３，巳＝１．４５，ｅ，＝１．９２。
這些常數(shù)都是由實驗或特例的計算給定。 方程（４．１）、（４．２）、（４．７）構(gòu)成了可以求解的封閉方程組，需要指出的是，方

程（４．７）中出現(xiàn)的六個常數(shù)Ｇ，ｃｌ，ｃ：，ｃ。，ｔ．，ｃ。是在�；^程中采用近似的結(jié)果， 顯然在�；胁捎貌煌募俣〞玫讲煌某�(shù)。
在標準ｋ一占模型中，雷諾應力的計算并不采用它的摸化方程，而是直接采用 Ｂｏｕｓｓｉｎｅｓｑ假設，即雷諾應力與平均場應變率應有如下關(guān)系：

一瓦＝ｕ譬＋等）一；毛七
其中：渦粘系數(shù)＿在ｋ—ｓ模式中取

（４．ｓ）

ｖ，＝Ｃ。生

而湍流動能（膏＝寺蠔“，）和湍能的耗散率ｓ采用�；匠蹋�

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

（４．９）

ｃａ．?。，

方程（４．９）（４．１０）與方程（４．１）（４．２）即構(gòu)成了標準ｋ—ｓ湍流模型，在該模型中 的常數(shù)取值如下：

Ｉ
｜

ｃ， 。．。９

Ｃ引 １．４４

Ｃ；２ ｌ，９２

盯★

盯‘

ｌ，Ｏ

１．３

ｋ一￡兩方程模型是目前應用比較廣泛的一種湍流模型。在推演過程中，采 用了以下幾項基本處理：（１）用湍能七反映特征速度；（２）用湍能耗散率占反映特征 長度尺度；（３）引進了■＝Ｃ。ｋ２／Ｆ的關(guān)系；（４）利用了Ｂｏｕｓｓｉｎｅｓｑ假定進行簡化。 因此，ｋ一占模型具有以下優(yōu)點：（１）通過求偏微分方程考慮湍流物理量的輸運過程， 即通過求解偏微分方程確定湍動特征速度與平均場速度梯度的關(guān)系，而不是直接將兩 者聯(lián)系起來：（２）特征長度不是由經(jīng)驗確定，而是以耗散尺度作為特征長度，并由求 解相應的偏微分方程得到。由于脈動特征速度和特征長度是通過解相應的微分方程求 得，因而ｋ一占模型在一定程度上考慮了流動場中各點的湍能傳遞和流動的歷史作用。
４７

河海大學博士論文

計算結(jié)果表明，它能比較好地用于某些復雜的流動，例如環(huán)流、渠道流、邊壁射流和 自由湍流，甚至某些復雜的三維流等。但是它也有局限性，主要有以下幾點：（１）采 用Ｂｏｕｓｓｉｎｅｓｑ假定，即采用了梯度型和ｕ各向同性的概念，因而使ｋ一占模型難以準

確地模擬剪切層中平均流動方向的改變對湍流場的影響；（２）引入了一系列的經(jīng)驗系 數(shù)，而這些系數(shù)都是在一定試驗條件下得來的，因而也限制了模型的使用范圍。 ４．５ＲＮＧ女一占模型 ＲＮＧｋ—ｓ模型由Ｙａｋｈｏｔ和Ｏｒｓｚａｇ于１９８６年應用重整化群（Ｒｅｎｏｒｍａｌｉｚａｔｉｏｎ
Ｇｒｏｕｐ）的方法導出，對于高雷諾數(shù)湍流，ＲＮＧｋ一占模型具有和標準ｋ一￡模型相同

的形式，只是在ｓ方程中出現(xiàn)了一個附加生成項，當流動快速畸變時，這一項將顯著 增大。ＲＮＧｋ一占模型與標準ｋ一占的主要不同之處在于它們的模型系數(shù)不同，標準 ｋ一￡模型系數(shù)由經(jīng)驗給定，而ＲＮＧ＾一￡模型系數(shù)均由理論精確計算得出。此外
Ｙａｋｈｏｔ和Ｏｒｓｚａｇ認為ＲＮＧ ｋ一占模型不需要采用壁面率，可直接積分到壁面，而文

獻【４６】作者在應用中發(fā)現(xiàn)，若采用ＲＮＧｋ—Ｅ模型直接積分到壁面，必須在壁面附近 劃分足夠細的網(wǎng)格，這必然增加計算工作量，尤其在三維工程湍流計算中巨大的網(wǎng)格 數(shù)量往往是難以接受的，因此類似于文獻［４６］的作法，本論文采用考慮壁面率的
ＲＮＧ女一占模型，并將ＲＮＧｋ—ｓ湍流模型推廣應用于三維湍流模擬。

ＲＮＧｋ—ｓ模型與標準的ｋ—ｓ模型具有相似的形式，該模型中的兩個標量方程
ｋ方程與ｓ方程形式如下：

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

（４．１１）

?０動６＇＋ｕ，毒＝云（（，＋與Ｏ＂。罷Ｏ＇Ｘｉ）一矗＋｛（ｅ．ｎ—ｅ：占）
程中的Ｒ項為ＲＮＧｋ一￡模型與標準ｋ～￡的主要區(qū)別，Ｒ＝

（４．１２）

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

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

平均應變鄹㈣珊７＝厝孝船均湖擱瞄湍流時間尺趔Ｅ，

河海大學博士論文

湍流粘性系數(shù)∥。＝ｐ?Ｇ。七二／５，叩。＝４．３８是印在剪切流中的典型值， 口＝０．０１５。此外．與標準女一ｓ模型的不同之處還有方程中的常數(shù)不是用經(jīng)驗方法 確定的，而是采用ＲＮＧ理論推導出的精確值，具體數(shù)值如下：
２

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

１?６８

與標準女一ｓ模型相比，ＲＮＧｋ一￡模型具有以下特點：

①ＲＮ６ｋ一￡模型的占方程中有一附加項Ｒ，該項代表平均應變率對￡的影響，可以 改善迅交流中的模擬精度： ②在ＲＮＧ≈一ｓ模型中包含了旋渦對湍流的影響，提高了對旋渦流的模擬精度； ③方程中的常數(shù)不是用經(jīng)驗方法確定的，而是利用ＲＮＧ理論推導出的精確值； ④ｋ一￡模型為高雷諾數(shù)模型，而ＲＮＧｋ一￡模型則同時適應于低雷諾數(shù)的情形，
此時需要在近壁面處作特別處理

綜上所述，ＲＮＧ七一ｆ模型比標準七一占模型具有更高的精度和更廣的應用領域。
４．６

ｋ—ｓ湍流模型中的邊界條件
對ｔ一￡方程求解需要給出合理的邊界條件，這些邊界條件的應用對所有的方程

基本上是相類似的，只是對于固體壁面。問題可能變得比較困難。在壁面上取七＝０，
ｊ６

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

對于高見流動問題，由于粘性底層很薄，無法在壁面處采用非常細密的網(wǎng)格模 擬，解決的方法是采用壁面函數(shù)（ｗａｌＩ Ｆｕｎｃｔｉｏｎｓ），該函數(shù)與速度分布的對數(shù)有關(guān)。湍
流邊界層的速度分布在圖４．１中給出，從圖中可以看出，在對數(shù)區(qū)域的速度分布為：
，，

１

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

（４。１３）

Ｆ

∞礦俗
ｍ

，

ｏ ２

５

１０

２０

５０

１００礦

圖４．１湍流邊界層速度分布

河海大學博士論文

斯∞與壁面平行嗍溉“。黼應力戤或稱為摩擦斌鏟層，
ｆ。為壁面切應力，Ｋ稱為卡門常數(shù)（Ｋａｍａｎ ｃｏｎｓｔａｎｔ）（ｒ：ＯＡｌ），Ｂ表示粘性底層厚
度的實驗常數(shù)（平板邊界層中Ｂ 化。
ｚ

５．２），”為網(wǎng)格點到壁面的距離，ｎ＋為ｎ的無因次

ｎ＋：型
“

（４．１４）

由于通常假定流體是局部平衡的，即湍動能的產(chǎn)生與耗散是近似相等的，因此，

“，＝ｃ：４√．ｉ｝

（４．１５）

由（４。１４）（４．１５）兩式可以推導出壁面上第一個網(wǎng)格點的速度與壁面切應力的關(guān)系

鏟倒；＝斛４ｒ壓赤
建立更加精確的方程。

（４．１６）

其中：Ｅ＝Ｐｏ，由于緊靠壁面的控制體的一個控制面在壁面上，因此在動量方程中

對于平行于壁面的控制體，其壁面上的切應力值是需要的，因此利用該邊界條件可以

當采用壁面率邊界條件時，女在壁面上的擴散通量取為ｏ，即婺：ｏ，關(guān)于耗
散率占的邊界條件可以根據(jù)局部平衡的假定，即壁面處產(chǎn)生與耗散近似相等推導出 來。在壁面區(qū)，湍動能的產(chǎn)生項可由下式計算：

只一。魯
ａＨ

（４．１７）

同時該區(qū)內(nèi)切應力近似為常數(shù)，為了計算出控制體中心處的耗散率，通過對（１）式
的對數(shù)速度分布求導，可得：

舀Ｏ，：÷：掣
、

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

該式與（４．１６）式聯(lián)立即可計算出耗散率。
當采用上述邊界條件時，近壁面處的控制體中是不需要求解占方程的，而是由下

河海大學博士論文

式直接計算

。Ｐ

，：竺豎
砌口 必須注意，上述邊界條件是使用只有當靠近壁面的第一個網(wǎng)格點處于對數(shù)區(qū)，即

月：＞３０時才有效。 此外，在入口邊界處，ｋ，６通常未知，一般給膏賦一個小量值，如１０一云２，而￡的

值一般由式占。譬確定。
４．７算例：水翼流場模擬 利用非正交網(wǎng)格同位網(wǎng)格系統(tǒng)中的ＳＩＭＰＬＥ算法對ＮＡＣＡ００１５型二維水翼粘流 場進行數(shù)值計算，湍流模型分別采用標準ｋ一￡模型和ＲＮＧｋ—ｅ，雷諾數(shù)Ｒｅ＝１．９７ｅ～， 在圖４．２中給出了不同攻角下升力系數(shù)曲線，并與文獻［９１］的實驗結(jié)果進行比較，在 小攻角下兩種湍流模型的計算結(jié)果與實驗值均比較一致，ＲＮＧｋ－￡模型略優(yōu)于標準
ｋ一￡模型。
１．４ １．２ 甚１ ０．８ ０．６
０．４

０．２ ０
Ｏ

２

４

６

８

１，１４ ｇ（ｄｅａ）

１（１

１６

１８

２０

２２

圖４．２
４．８小結(jié)

ＮＡＣＡ００１５水翼升力系數(shù)一攻角變化曲線

本章討論了湍流模式理論研究的工程背景和意義，詳細給出了標準ｋ一占模型和 ＲａＮＧｋ一￡模型，比較了兩種模型的特點及應用范圍，為更好地處理固壁邊界，采用 與壁面率相結(jié)合的ＲＮＧｋ一占模型模擬湍流流動，通過對ＮＡＣＡ００１５型水翼湍流場的 模擬表明，ＲＮＧｋ—ｓ模型能夠更好地模擬出流場的湍流特性。

河海大學博士論文

第五章運動界面數(shù)值模擬技術(shù)
５．１運動邊界模擬技術(shù)簡介 自由面模擬技術(shù)在水利工程、船舶工程、海洋工程、機械工程等領域有著深刻的 現(xiàn)實意義和應用價值，隨著Ｈｉｒｔ等人提出用ＶＯＦ方法模擬自由面追蹤技術(shù)，目前已 經(jīng)發(fā)展和形成了許多實用軟件，如ＶＯＦ．３Ｄ，ＮＡＳＡ．ＶＯＦ，ＳＯＬＡ．ＶＯＦ等，同時推廣 到一般介質(zhì)面的運動界面模擬，該理論已經(jīng)成功應用于半導體工業(yè)、噴墨技術(shù)等現(xiàn)代
高新技術(shù)產(chǎn)業(yè)。據(jù)不完全統(tǒng)計，１９９９年度的Ｊ．Ｃｏｍｐｕｔ．Ｐｈｙｓ．雜志涉及運動界面追蹤問

題的文章不少于３０篇，可見運動邊界數(shù)值模擬的研究和分析方法受到越來越多的關(guān) 注，但是由于運動界面的形狀和整體位置不斷變化，所以精確模擬的難度非常大。 運動界面的模擬模型一般分為兩大類：歐拉模型和拉格朗日模型。歐拉模型的特 征是坐標系或者固定或者按給定的方式移動以適應不斷變化的計算域形狀，其主要優(yōu) 點是可以適應內(nèi)部界面的大變形和大幅度移動而能保持一定的計算精度。利用歐拉模 型需要解決的主要問題是用什么形式的控制方程最方便，如何選取合適的數(shù)值離散方
法和最佳的內(nèi)界面跟蹤技術(shù)。拉格朗日模型的特征是坐標隨流體雨運動，但每個網(wǎng)格 單元內(nèi)包含的流體質(zhì)點保持不變。該模型的優(yōu)點是可以精確跟蹤內(nèi)界面并準確描繪出 運動界面的形狀，便于使用界面的運動學和動力學邊界條件；拉格朗日模型還可以在

運動界面附近進行網(wǎng)格加密，以精確模擬運動界面附近的流場內(nèi)部結(jié)構(gòu)，同時精細網(wǎng) 格區(qū)可隨界面移動，大大提高了運動界面的模擬精度。拉格朗日模型的主要不足是隨 著計算的進行，網(wǎng)格會出現(xiàn)重疊、纏繞等現(xiàn)象，隨著網(wǎng)格越來越不規(guī)則，計算精度也
越來越低。 在歐拉模型中，模擬運動邊界的典型方法有三種：即ＭＡＣ法，ＶＯＦ法，Ｌｅｖｅｌ
Ｓｅｔ

方法。ＭＡＣ法的基本原理是在計算區(qū)域中包含流體的網(wǎng)格內(nèi)布設“標記點”（Ｍａｌ（ｅｒ），

標記點沒有質(zhì)量，不參與流場計算，只是用來標記流體質(zhì)點的運動軌跡，隨著流場計 算的進行，可獲得每個標記點的速度，并由速度確定出標記點的最新位置，根據(jù)標記 點的最新位置即可確定哪些網(wǎng)格單元被流體充滿，哪些網(wǎng)格單元被流空，從而確定運 動界面或自由液面在下一時刻的位置；ＶＯＦ法是Ｈｉｒｔ和Ｎｉｃｈｏｌｓ［９２１于１９８１年首先提

出，通過對主場的計算采用了Ｄｏｎｏｒ－－Ａｃｃｅｐｔｏｒ（施主一受主）的概念、Ｕｐｗｉｎｄ特征
思想和補償效應，成功地模擬出潰壩和涌浪（ｂｒｏｋｅｎｄａｍ，ｂｒｅａｋｉｎｇｂｏｒｅ）自由面，避 開了對于自由面問題通常采用的高度函數(shù)法和計算量十分浩大的Ｍａｒｋｅｒ點的方案，
引入了一個流體體積函數(shù)方程，與流體運動基本方程聯(lián)合求解，開辟了自由面計算的

新途徑。盡管當初的格式、方法是比較粗糙的，但其基本思想至今一直受到重視，特 別是對運動界面的重構(gòu)方法開創(chuàng)了一個新的研究方向：１９８８年，Ｏｓｈｅｒ［９３１等人提出了

河海大學博士論文

一種零等值面方法（Ｌｅｖｅｌ Ｓｅｔ方法），在許多復雜的界面追蹤問題中獲得成功，該方

法的主要思想是將運動界面定義為一個函數(shù)的零等值面（線），即妒（ｘ，，）＝０，然后始 終保持它是零等值面，并且在界面附近保持單調(diào)。
５．２

ＶＯＦ方法 ＶＯＦ方法的基本思想是定義一個函數(shù)Ｆ，在有流體的點取值為１，在無流體的點

取值為０，這樣在一個單元內(nèi)，Ｆ的平均值就可以代表單元內(nèi)流體所占的份額，如果 Ｆ的均值為１，則單元內(nèi)充滿流體，如果單元內(nèi)無流體，則Ｆ的值為零。而Ｆ值在０

和１之間則表示該單元必然包含自由液面，所以，ＶＯＦ法提供了與ＭＡＣ法類似的有 關(guān)內(nèi)部界面的信息。由于ＶＯＦ方法只需對每一個單元提供一個內(nèi)存就可以了，不像 ＭＡＣ法那樣需要對標志點提供內(nèi)存，從而節(jié)省的計算空間，提高了計算效率。 ５．２．１控制方程
流體體積分數(shù)Ｆ隨時間的變化過程由下面的方程控制：

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

缸

（５．１）

咖

對于不可壓縮流體，控制方程的守恒形式為：

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

孤

（５．２）

磚

自由面模擬問題的控制方程一般由上面的流體體積函數(shù)方程與流體運動基本方
程構(gòu)成。

５．２．２流體體積分數(shù)Ｆ的輸運模型 以Ｈｉｒｔ和Ｎｉｃｈｏｌｓ建立的施主一受主模型（ｄｏｎｏｒ—ａｃｃｅｐｔｏｒ）為例說明Ｆ輸運的基本思
想。

考慮網(wǎng)格的右側(cè)邊界，國時段內(nèi)，通過該網(wǎng)格面的流體和空氣的體積通量 Ｖ＝“田，其中”為網(wǎng)格面上的法向速度，其符號決定了網(wǎng)格是“施主網(wǎng)格”還是“受 主網(wǎng)格”，若“的符號為正，則網(wǎng)格面上游是施主網(wǎng)格，下游是受主網(wǎng)格。在一個時 間步長內(nèi)通過網(wǎng)格面的Ｆ的通量等于６Ｆ與網(wǎng)格面積的乘積。８Ｆ由下式計算：

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

（５．３）

其中第一項表示輸運的流體的量，第二項表示施主網(wǎng)格所擁有的流體的量。

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

（５．４）

河海大學博士論文

式中，下標Ａ代表受主網(wǎng)格，下標Ｄ代表施主網(wǎng)格，下標ＡＤ則表示可以是施

主網(wǎng)格，也可以是受主網(wǎng)格，根據(jù)邊界面上流動的方向確定，比如，同一個網(wǎng)格相對
于不同網(wǎng)格邊界的情況。（５．３）式中的Ｍｉｌｌ是為了防止通過界面的通量超過施主網(wǎng)格

所能提供的Ｆ的量，（５．４）式中的ｍａｘ是為了保證當空氣的通量超過可能提供的量時，
需要補充流體通量，這個補充的量就是ＣＦ。

在圖５．１中給出了施主一受主模型通量輸運的基本形式，當流動通過一個垂直的單 元邊界面時，在（ａ）中給出了施主網(wǎng)格與受主網(wǎng)格的定義：當下游網(wǎng)格為受主網(wǎng)格

即ＡＤ＝Ｄ時，則凹＝０，Ｆ的通量為Ｆ：Ｆｏｌｖｌ，其中Ｖ＝�。叮簦浦涤糜诖_定網(wǎng)格
邊界面通過流體的那部分斷面面積，當然要求川＜Ａｘ以保證施主網(wǎng)格內(nèi)的流體不會
流空，見（ｂ）：當ＡＤ＝Ａ時，則施主網(wǎng)格內(nèi)的Ｆ值用來確定網(wǎng)格邊界面上通過流體的
斷面面積（ｃ），此時施主網(wǎng)格內(nèi)的所有流體都將流失，即所有存儲在虛線和邊界面之 間的流體和空氣均進入受主網(wǎng)格；在（ｄ１中表明，當ＡＤ＝Ａ，但通過網(wǎng)格面的流體通

量超過ＹｏＩＶＩ時，則需要補充額外的流體以滿足需要的通量，其中虛線和邊界面之間
的額外流體就等于ＣＦ。

（ｂ）

（ｃ） 圖５．１
（ａ）

（ｄ）

施主．受主模型中Ｆ通量的計算方法舉例

施主?受主模型網(wǎng)格的配置，（ｂ）、（ｃ）、（ｄ）中雙陰影線區(qū)域表示實際傳輸?shù)模频牧?br />
事實上，使用施主網(wǎng)格中的Ｆ值還是受主網(wǎng)格中的Ｆ值來確定邊界面上通過流體 的斷面面積完全取決于自由液面的平均定向，這可分為兩種情況①當自由液西基本上

河海大學博士論文

沿其法向移動時，ＡＤ＝Ａ，②當自由液面基本上沿其切向移動時．ＡＤ＝Ｄ。此外如
果受主網(wǎng)格是空的，或者施主網(wǎng)格的上游是空的。則一律有ＡＤ＝Ａ，也就是說，施 主網(wǎng)格必須首先補充流體，然后才能向下游的空單元輸運流體。 用通量值５Ｆ乘以流體通過邊界的斷面面積就可以得到輸運流體的量，從施主網(wǎng)

格中減去該通量值添加到受主網(wǎng)格，然后對所有的網(wǎng)格邊界重復上述步驟即可完成Ｆ
的輸運計算。 在上述計算過程中可能會出現(xiàn)Ｆ值小于零和大于ｌ的情況，因此必須對Ｆ的 計算值歸整到０和１之間，同時將由此引起的誤差進行累計，以保證總的流體體積守

恒。此外由于在數(shù)值計算中產(chǎn)生的舍入誤差使得Ｆ的值不可能正好是１或０，為此通 常設定一個很小的數(shù)如占，＝１０一，當Ｆ小于卸時，網(wǎng)格是空的，當Ｆ大于１一￡時，
網(wǎng)格是滿的。

５。２。３運動界面的重構(gòu)方法 Ｈｉｒｔ和Ｎｉｃｈｏｌｓ的貢獻主要體現(xiàn)在：～是作者采用的施主一受主（Ｄｏｎｏｒ－Ａｃｃｅｐｔｏｒ） 方法給出了界面通量輸運計算的初始格式，同時該方法不僅考慮本網(wǎng)格單元上的流體
體積分數(shù)，而且用到相鄰網(wǎng)格上的流體體積分數(shù)，從而提高了計算精度：作者的另一 個重大貢獻是提出了一種自由面的重構(gòu)方法，因此掌握運動界面的重構(gòu)方法必須從

Ｈ—Ｎ模型出發(fā)。

一赫
．÷一１Ｉ歷黝嬲勉
圖５．２施主．受主模型中的自由面重構(gòu)方法

該模型重構(gòu)運動界面的基本思想是：將流體自由面看作是局部的單值函數(shù) ｒ（ｘ）或Ｘ（ｙ），計算中需要用到９個網(wǎng)格的信息，如圖５．２。首先計算ｆ一１，ｆ，ｆ＋ｌ網(wǎng)格

列的ｒ值和Ｊ一１，，，Ｊ＋１網(wǎng)格列的ｚ，值，估算出每個網(wǎng)格上自由面的斜率值譬，譬， 靠聊
然后根據(jù)流體體積分數(shù)和斜率的大小確定網(wǎng)格（ｆ，Ｊ）上的自由面位置和方向，即

河海大學博士論文

Ｋ＝∑Ｆ，ｋＳｙＩ，

？＿ｉ一１，ｉ，ｉ＋１

（５ ５）

（氅，：ｏ暨斗
、ｄｒ“

（５ ６） 、…７

蘇…＋２融，＋蘇Ｈ

ｘ，＝∑ｒ，ｋ４ｖ＾，

，＝產(chǎn)１，Ｊ，ｊ＋ｌ

（５ ７）

（箏，：善婆蘭業(yè) Ｉ一｝
、咖。 ￥Ｈ＋２�。妫А海�

理論上講，以上兩個導數(shù)應該互為倒數(shù)，但數(shù)值計算結(jié)果則不然，作者的意圖是比較
它們絕對值的大小，來決定自由面在該網(wǎng)格內(nèi)是“水平”的還是“垂直”的。即如果

＝一

（５。，

ｌ：，．６， 、 。

引爿
則定義自由面是水平的；否則自由面是垂直的，如圖５．２所示。

（５．９）

盡管Ｈ．Ｎ運動界面重構(gòu)方法比較簡單，但它卻起到了拋磚引玉的開創(chuàng)性作用， 在這一思想的啟迪下，自由面的重構(gòu)技術(shù)得到了迅速的發(fā)展，比如：Ｃｈｏｆｉｎｌ９４】引入

ｌ，４和３／４單元以連接水平和垂直的自由面單元；Ｐｕｃｋｅｔｔ［９５１用離散的線段來表示自由
面，通過相鄰網(wǎng)格的體積分數(shù)預測自由面的法向；Ｈｏｎｇｆａｎｇｗｅｎｌ９６ＪＮＮ斜線段來代替 自由面，其中斜線段方程中的系數(shù)由自由面網(wǎng)格周圍８個網(wǎng)格的體積分數(shù)確定：

Ａｓｈｇｒｉｚ和Ｐｏｏ的ＦＬＡＩＲ方法ｐ�！縿t利用每兩個水平相鄰網(wǎng)格的體積分數(shù)求出斜線段來 模擬自由液面，但只適用于均勻網(wǎng)格情況。Ｒｕｄｍａｒｌｌ９”對上述兩類自由面重構(gòu)算法進
行了比較計算指出，利用斜線段來逼近自由液面比其它方法精度更高，這些方法的基

本思想都是首先確定網(wǎng)格內(nèi)運動界面的具體位置及其運動方向，然后確定運動界面在 下一個時刻的新的位置，同時更加注意運動界面與網(wǎng)格線的交點的位置，因此得到的 運動界面也更加精確，普遍采用的格式是用網(wǎng)格內(nèi)具有法向量的斜直線段Ｙ＝ａｘ＋ｂ 來逼近，也有用二次曲線Ｙ＝Ｏ．Ｘ２＋ｂｘ＋ｃ來構(gòu)造運動界面１９８１，只是考慮得更加精細而
已。

本論文利用ＦＬＡＩＲ方法重構(gòu)自由面，并對該方法進行了推廣。
ＦＬＡＩＲ自由面重構(gòu)方法的基本思想是在兩個相鄰的網(wǎng)格內(nèi)構(gòu)造一條傾斜的直線

段來近似該網(wǎng)格邊界的界面，然后計算透過該網(wǎng)格邊界的流體逶量，作為修正流體體
積分數(shù)的數(shù)值通量，由于界面構(gòu)造僅涉及到兩個網(wǎng)格，故需要根據(jù)網(wǎng)格內(nèi)的流體體積

分數(shù)劃分為許多種情況分別計算，各個網(wǎng)格內(nèi)的流體體積分數(shù)可能有以下四種情況，
見圖５．３。

河海大學博士論文

（Ｉ）施主網(wǎng)格是滿網(wǎng)格，即瓦＝１ （２）施主網(wǎng)格是空網(wǎng)格，即屹＝０ （３）施主網(wǎng)格是半網(wǎng)格，受主網(wǎng)格也是半網(wǎng)格，即０＜ＦＤ＜１，０＜只＜１ （４）施主網(wǎng)格是半網(wǎng)格，受主網(wǎng)格滿網(wǎng)格或半網(wǎng)格，即０＜％＜１，只＝０，只＝０

（ｃ）

ｉ－，？Ｔ Ｔ二

二

ｊｊｉ彩
（ｅ） （ｆ）

∥％／

圖５．３相鄰網(wǎng)格間的關(guān)系

情況（１）、（２）比較簡單，情況（３）可以分為１６種，經(jīng)過轉(zhuǎn)換可以歸結(jié)為４種，

見圖５．４，即￡≥厶，且目標流體在運動界面以下。此時只要確定了直線的斜率及位
置即可計算出通過網(wǎng)格邊界的流體輸運量，如對情況（ａ），設左邊是施主網(wǎng)格，則舀 時間內(nèi)輸運到相鄰網(wǎng)格的體積量Ｖ＝”．西，即圖中的斜線陰影區(qū)域的面積；如果右邊 是施主網(wǎng)格，則應取右邊網(wǎng)格內(nèi)靠近網(wǎng)格邊界對應區(qū)域面積，施主網(wǎng)格和受主網(wǎng)格由 網(wǎng)格邊界上的速度方向確定。

河海大學博士論文

笏

］赫勿
（ｃ， （ｄ）

圖５．４

ＦＬＡＩＲ模型運動界面重構(gòu)方法

情況（４）最為復雜，以圖５．５（ａ）為例，施主網(wǎng)格（ｆ√）是半網(wǎng)格，而受主網(wǎng)格（ｉ＋ｌ，，）

是空網(wǎng)格，可以用網(wǎng)格（ｆ，『）垂直方向的兩個相鄰網(wǎng)格分別按照第三種情況構(gòu)造界面，
然后求通過網(wǎng)格邊界流入受主網(wǎng)格的流體體積輸運量，同樣可以歸結(jié)為圖５．５（ｂ）所示 的四種情況：若受主網(wǎng)格是空網(wǎng)格，則取靠右邊的寬度為“．西的區(qū)域面積作為輸運

量，如果受主網(wǎng)格是滿網(wǎng)格，則取左邊的寬度為Ｉｔ．毋的區(qū)域面積作為輸運量，即保
證流體分布靠近滿網(wǎng)格一邊。

Ｅ瓢 簍貔 遴～絲擻
（ａ）

（ｂ）

圖５．５運動界面的四種重構(gòu)方法

下面以情況（３）為例給出ＦＬＡＩＲ模型構(gòu)造運動界面的方法，見圖５．４，并將該方法

推廣到非均勻網(wǎng)格，丘，兀分別表示左、右網(wǎng)格內(nèi)的體積分數(shù)，見圖５．６。

河海大學博士論文

Ｃａｓｅ １

Ｃａｓｅ ３

Ｃａｓｅ ４

圖５．６
Ｃａｓｅｌ

運動界面結(jié)構(gòu)形式

不妨假設兩個相鄰網(wǎng)格中的體積分數(shù)滿足：Ｌ≥以，自由面認為是與網(wǎng)格邊界相交 的直線段，故自由面可由Ｙ＝刪＋６來描述，其中系數(shù)口，ｂ由已知的體積分數(shù)Ｚ，五確
定。因此可得到下面的方程組：

厶?血，緲＝裊，（ａｘ＋６協(xié)＝ｉ１口（ｘ｝ｘ２Ｉ）＋ｂ（ｘｉ－ｘ，＿Ｉ） 兀－Ａｘｉ＋ｌ緲＝璧＋１（ａｘ＋ｂ）ｄｘ＝去ｎ（毒ｌ一＃）＋６（鼉＋Ｊ—ｘｉ）
聯(lián)立求解，得到

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

２

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

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

河海大學博士論文

０＝ａＸｂ＋ｂ

聯(lián)立上述三個方程，求解可得

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

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

ｂ＝２ｌ帆＋如）一

緲

缸？

ｃａｓｅ３

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

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

厶止Ⅲ每＝曇ｃｃ磊，一ｘ；）＋ｂ（ｘ“
聯(lián)立求解可得：

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

—ｚ石ｉ面刁‘面習Ｆ麗ｉ西可了了】．緲
其中５：二Ｌ

血．

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

ｘ，＋ｚ』＿１）

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

Ｐ：

！型二�。。剑�
２ｘＪ＋Ｉ—ｘｆ一３ｘｆ—

ｃａｓｅ４

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

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

河海大學博士論文

＾緲．缸Ｍ＝曇僻一ｘ；）＋６（Ｘｂ－－Ｘｉ）
聯(lián)立求解得：

口：坐
ｘｄ一工６

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

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

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

ｃ＝２眈?缸，＋＾－缸ｉ＋１＋Ｘｉ＋１）（ＸｉＹｊ－Ｉ＋厶血州?緲）＋等緲
容易驗證，對于尺度為ｈ的方形網(wǎng)格，文獻［６】中的結(jié)果即為本方法的特例。 顯然在確定了ｎ，ｂ之后，利用Ｙ＝似＋ｂ即可完成自由液面的重構(gòu)。

本文研究的運動界面是自由液面，自由液面單元的壓強Ｐ。需要利用內(nèi)插（或外

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

的內(nèi)部單元中心間的距離，ｄ是自由液面至內(nèi)部
鄰接單元中心的距離，用于內(nèi)插的鄰接單元的選 取標準是，聯(lián)接自由液面單元中心和所選取內(nèi)部 單元中心的聯(lián)線最靠近自由液面的法線，該單元 又稱內(nèi)插鄰接單元。

、平
／ｐ¨

—、～

王ｉ塵
ｐ＿

圖５ ７自由面上的壓強邊界條件

河海大學博士論文

５．２．５運動界面的寬度與等值線的給定方法

由于運動界面的初始位置是由流體體積分數(shù)定義的㈣，隨著流場的變化和時間的
推移，計算所得到的自由面輪廓是從流體體積分數(shù)為０的高度變到１的高度，即計算 流場中得到的自由面為Ｆ∈［Ｏ，１１的模糊界面，其寬度可以跨越幾個阿格。如圖５．８

圖５．８運動界面的寬度

所謂等值線，就是根據(jù)不同問題的實際情況，流體體積分數(shù)�。郏铮臁繀^(qū)間內(nèi)的某個 值，如令Ｆ＝Ｏ ５的等值線作為自由面的位置，當然該值不是固定的，有可能在其它 情況下取Ｆ＝０．６更合適。這是采用單一定值等值線確定自由面的方法。此外，還可 以利用Ｆ值的百分比寬度來給出邊界面的形狀，如９６％寬度的界面；顯然，如果有 某種方法可以進一步將具有一定寬度的界面重構(gòu)為單值輪廓線，則得到的自由西位置 將具有更高的精度。 ５．２．６數(shù)值穩(wěn)定性條件 （１）對流穩(wěn)定條件 為了正確的模擬自由液面或內(nèi)界面的變化，我們只考慮相鄰單元之間的通量，因 此，時間步長必須滿足下面的不等式，

趴幽協(xié)尚｝
上式要求在計算網(wǎng)格確定之后，時Ｉ＇ｕｑ步長西的選取必須使流體在該時段內(nèi)流過的 距離不能超也網(wǎng)格的尺寸，由于ｒａｉｎ是對所有的單元取的，故在實際計算中的毋一 般取上式右端的ｌ／４～１／３。

河�；饘W博士論文

（２）擴散穩(wěn)定性 對于粘性流體，要求在每一個時間步內(nèi)，動量的擴散同樣不超過一個網(wǎng)格單元。
線性穩(wěn)定分析給出如下的限制條件：

Ｖ藤２篇Ａｘ

２

．２＋△ｙ．

５．３運動邊界模擬算例 算例１：Ｚａｌｅｓａｋ模型 為檢驗本文所建立的的自由面重構(gòu)算法的可靠性，對經(jīng)典的Ｚａｌｅｓａｋ界面模型在 平面旋轉(zhuǎn)流場中的運動進行了計算模擬，計算狀態(tài)為：計算區(qū)域【０，１］×【０，ｌ】，網(wǎng)格
節(jié)點數(shù)２００ｘ ２００．旋轉(zhuǎn)速度場為：ｕ（ｘ，Ｙ）＝—石（ｙ—Ｏ．５）＇ｖ（ｘ，力＝ｘ（ｘ—Ｏ．５），初始位置

為一個下方有缺口的圓周；計算時間步長為母＝Ｏ．０００５。在圖５．９中分別給出了經(jīng)過 四個時間步長ｆ＝０．５ｓ，ｆ＝ｌｓ，ｆ＝１．５ｓ，ｆ＝２ｓ對應的二維界面圖形，同時給出了文獻［１００１
中采用Ｈｉｒｔ和Ｎｉｃｈｏｌｓ模型的計算結(jié)果。可以看出，Ｈｉｒｔ和Ｎｉｃｈｏｌｓ方法給出的界面

比較粗糙，無法準確地描述運動后的界面形狀，而本文所給出的方法則能夠獲得非常
好的結(jié)果。

，＝０．５ｓ

ｆ＝１．０ｓ

ｆ＝１．５ｓ

ｔ＝２．０ｓ

（ａ）本方法

ｔ＝０．５ｓ

ｆ＝１．０ｓ

ｔ＝１．５Ｊ

ｆ＝２．０ｓ

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

Ｚａｌｅｓａｋ模型的計算結(jié)果

算例２：均勻流場中圓形界面模型
計算區(qū)域及網(wǎng)格設置與算例ｌ相同，圓的初始位置為：圓，１１，（０．３，０‘３），半徑ｒ＝０．２

河海大學博士論文

均勻速度場為Ｕ＝０．２，ｖ＝０．２，經(jīng)過四個時間步長圓的形狀及位置如圖５．１０。相對于 初始圓周，面積的變化量分別為ａｓ＝Ｏ．０００１５４，０．０００１５４，０．０００１５３，０．０００１５８，其中每一 步的面積計算是根據(jù)相應的面積分數(shù)相加而得到的。

ｆ＝０．５ｓ

ｔ＝１．０５

，＝１．５ｓ

ｆ＝２．Ｏｓ

圖５．１０均勻流場中圓的移動

算例３：剪切流場中圓形界面的運動 平移和旋轉(zhuǎn)流不會使運動界面的形狀發(fā)生變化，而在真實的流體問題中會發(fā)生界

面不斷變形的復雜物理現(xiàn)象，為了考察本文建立的算法對各種現(xiàn)實流場的模擬效果，
計算了剪切流場中圓形界面的變形運動，剪切速度場為
ｕ（ｘ，ｙ）＝７ｒｃｏｓ（ｘ（ｘ—Ｘｏ））ｓｉｎ（石（ｙ—Ｙｏ）） Ｖ（ｘ，ｙ）＝一，ｒｓｉｎ（ｘ（ｘ—Ｘｏ））ｃｏｓ（ｘ（ｙ一），０））

計算區(qū)域為【Ｏ，１］×［Ｏ，ｌ】，初始界面為圓心在（Ｏ．３，Ｏ．３）、半徑為０．２的圓周，見 圖５．１１（ａ），計算結(jié)果給出了ｒ＿Ｏ ５ｓ，忙ｌｓ，ｔ＝１．５５，，＝２ｓ圓形界面的變化情況（圖 ５．１２），同時以ｆ＝２ｓ的結(jié)果作為計算初值，將速度反號進行反剪切２秒，結(jié)果在圖５．１１
（ｂ）中給出。

ｌ

Ｉ

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

＼ｊ

（ａ）

（ｂ）

圖５．１ｌ

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

河海大學博士論文

ｒ＝０．５ｓ

ｒ＝１．Ｏｓ

ｆ＝１．５ｓ

ｆ＝２．０ｓ

圖５．１２剪切流場中界面的位置和形狀

５．４小結(jié)

本章討論了運動界面的模擬技術(shù)，重點研究了ＶＯＦ方法中運動界面的高精度重
構(gòu)算法，以ＦＬＡＩＲ方法為基礎，建立了可適應于非均勻網(wǎng)格的自由面重構(gòu)模型，通 過對旋轉(zhuǎn)流場中的馬踢鐵形界面運動、平移流場中圓形界面的運動、剪切流場中圓形 界面的運動進行模擬表明，本章建立的運動界面重構(gòu)模型相對于其它方法能夠更加準
確地確定運動界面的形狀和位置。

河海大學博士論文

第六章溢流壩過壩水流二維流場數(shù)值模擬
６．１溢流壩過壩水流的研究方法
溢流壩的水力設計對于整個水力工程建設具有重要意義，其中溢流壩的流量系 數(shù)、自由水面線的形狀和位置、壩面壓力及壩面流速等都是水力設計中的重要參數(shù)， 也是選擇合理溢流壩型的重要依據(jù)。這些數(shù)據(jù)多年來主要通過水力學模型試驗或由經(jīng) 驗公式而獲得，一方面要耗費大量的人力、物力和財力，同時受到測量條件的限制不

可能給出整個流場區(qū)域的全部信息。隨著計算技術(shù)的迅速發(fā)展，利用數(shù)值模擬方法求
解溢流壩過壩水流場的各種計算模型相繼產(chǎn)生，如有限差分方法、邊界元方法［６０，１０７】、 有限元方法［１０３，１０４，１０５］、邊界擬臺曲線坐標變換方法【‘吲等。
有限差分法簡易明了，易于理解，其主要不足是不能很好地處理不規(guī)則邊界，因

此在過壩水流數(shù)值模擬中的應用受到了限制。

有限元法比有限差分法能更好地擬合復雜邊界，因而在過壩水流的數(shù)值模擬中得
到廣泛而成功的應用。其中最早將有限元方法應用到過壩水流求解的是日本的池川昌 弘ｌＩ“Ｊ，并對自由面的迭代提出了兩種方法。該方法的基本思想是首先假設一個流量， 分別用緩流和急流的迭代方法進行自由面迭代，如果上游緩流邊界能很好地與下游急 流邊界交匯并光滑地過渡到下游，則該流量即為正確流量：如果不能很好地交匯，則

需要重新假設流量進行試算。然而通過計算，他認為兩種流態(tài)邊界不能很好地匯合，
因而未給出最終成果。文獻［１０９］從可變區(qū)域的概念出發(fā)，把決定自由水面形狀的坐 標和各節(jié)點的流函數(shù)作為未知量形成方程組同時求解。文獻『１ ｌｏ］用有限元法計算了

過壩水流，提出了表面流帶的概念，將一維水流的臨界流量概念和水面線計算方法應
用于二維，找出了過壩流量和自由表面之間的聯(lián)系，克服了尋找臨界點的困難，文獻

［１１１】分別給出了用有限元泛函極值加權(quán)余量法離散格式，并用此格式求解過壩水流問
題的方法，得到了與實驗較一致的結(jié)果。文獻ｆ１１２］提出了一種求解具有未知能頭的
過壩水流迭代方法，根據(jù)自由表面能量極值條件，建立了確定未知水頭的數(shù)學關(guān)系，

從而得以實現(xiàn)水頭與自由面同步迭代求解。這一方法既適用于高壩，也適用于低堰。

由此可見，用有限元法對過壩水流進行數(shù)值計算已是碩果累累，占有很重要的位 置。用該方法對流量函數(shù)和自由表面位置的計算也是相當成功的，所得結(jié)果能夠基本
滿足水利工程的設計要求。 邊界元法起步較晚，在近２０年來才得到了較快的發(fā)展，最早應用于勢流計算， 同時也是應用最為成功的領域，但把邊界元法直接應用于過壩水流計算目前還不多

見，文獻［１１３］求解了溢流壩頂自由出流問題，文獻［１１４］求解了水力學中的不定邊界
流場問題。

河海大學博士論文

邊界擬合曲線坐標變換法在近年來被更多地應用到過壩水流的數(shù)值模擬，該方法 是通過邊界擬合坐標變換，將不規(guī)則的物理流動區(qū)域變換成規(guī)則的矩形計算區(qū)域，同 時將控制方程作相應變換，然后采用比較簡單的有限差分法進行求解。文獻［１０６１干０ 用該方法計算了過壩水流問題，計算結(jié)果與實驗吻合良好，文獻［２３］Ｎ用曲線坐標變 換對溢流壩反弧段水流進行了模擬。 需要注意的是，關(guān)于過壩水流的計算主要還是以勢流理論為基礎，本章利用前面 建立的非正交網(wǎng)格系統(tǒng)中的ＳＩＭＰＬＥ算法對溢流壩過壩水流二維粘流場進行計算，采 用分區(qū)網(wǎng)格技術(shù)對計算區(qū)域進行網(wǎng)格劃分，湍流模型采用ＲＮＧｋ一占模型，自由水面 線的形狀及位置采用ＶＯＦ方法確定，分別對高水頭溢流壩的水面線形狀、壩面壓力、 流量系數(shù)及壩面流速分布等進行計算。 ６．２計算模型 目前世界上廣泛使用的壩面溢流曲線是由美國陸軍工程兵團水道實驗站（ＷＥＳ） 提供的ＷＥＳ曲線‘“”，該型曲線具有工程量小、無有害的負壓及泄流能力大的特點， 一般表達式為：

日Ｙ－ｃ（≯
其中：爿一不包括行近流速水頭的溢流壩上的水頭高度 ｃ，６一與上游壩面傾斜坡度有關(guān)的系數(shù)，
ｘ，Ｙ一以溢流壩最高點為原點的坐標值。 在坐標原點上游，壩面曲線一般由兩個或三個圓弧線組成，其中三圓弧曲線的性能優(yōu) 于二圓弧曲線。

本文選用的壩面曲線為三圓�。祝牛忧�Ｘ‘”＝２．０川”Ｙ，見圖６．１，６．２，其中
上游壩面坐標在表６．１中給出，溢流壩高Ｈ。＝７６．２米，設計水頭Ｈ。＝９．１４米，計算 區(qū)域與網(wǎng)格劃分見圖６．３。

６７

河海大學博士論文

圖６．１

三圓弧溢流壩上游壩面

圖６．２三圓弧溢流壩壩面曲線 表６．１
ｘ ｌ Ｈａ
．０．Ｏ

三圓�。祝牛訅蚊嫔嫌巫鴺耍ㄗ鴺嗽c取在壩頂中心線上）
．０．０５ ．０．１０ ．Ｏ．１５ ．Ｏ．１７５ ．０．２０ ．Ｏ．２２

Ｙ｜Ｈ ｄ

０．Ｏ

０．００２５

Ｏ．０１０ｌ

０．０２３

０．０３１６

０．０４３

０．０５５３

ｘ ｆ Ｈｄ

．Ｏ．２４

．Ｏ．２６

．０ ２７６

．０．２７８

．Ｏ．２８

．０．２８１８

Ｙ｜Ｈ ｄ

０ ０７１４

０ ０９２６

０．１１ ５３

０．１１９０

０．１２４ｌ

０．１３６０

河海大學博士論文

翻６．３計算區(qū)域與網(wǎng)格劃分

從圖６－３中可以看出，整個流場分為三個區(qū)域進行網(wǎng)格劃分，網(wǎng)格類型均為結(jié)構(gòu) 化網(wǎng)格，其中Ｉ區(qū)網(wǎng)格數(shù)為７５００，ＩＩ區(qū)網(wǎng)格數(shù)為４０００，為保證壩面壓力的計算精度， 對ＩＩＩ區(qū)網(wǎng)格進行加密，網(wǎng)格數(shù)為８０００。 入１３邊界速度按設定流量計算給出，出口邊界條件取物理量沿水流方向梯度為 零：所有的氣體邊界都定義為壓力邊界，即壓力為０，湍動動能ｋ和耗散率ｅ按經(jīng)驗 公式確定，固體邊界速度按無滑移條件給定，并引入壁面函數(shù)技術(shù)，自由面的位置利 用ＶＯＦ方法確定：時間步長△ｆ＝Ｏ．００２，迭代次數(shù)２００００。 ６－３自由水面線計算 水面線位置的計算按表６．２進行，水面線形狀在圖６．４中給出，與試驗結(jié)果的比 較在圖６．ｊ中給出，可以看出，計算出的水面線位置與試驗值非常一致。
表５．２水面線計算
Ｈ。／ＨＪ
（初值）

０．５０

１．Ｏ

１．３３

１．５

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

０．３２９８７７

１．Ｏ

１．５７８２０２

１．９１３１３７

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

８０，８５

８５．３８

８８．２２

８９．６７

０．５０８７

１．００４４

１．３１３

１．４７８

河海大學博士論文

圖６．４不同壩頂水頭下的水面線（紅色為水）

８０

７５

７０

６０ ｌＯ ０ １０ ２０ ３０

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

河海大學博士論文

９０

８５

８０

５
｝

７５

７０

６５

６０
０ １０

２０

３０

Ｘ（ｍ）

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

８５

８０

占７５
｝

７０

６５

６０ ０ １０ ２０

Ｘ（ｍ）

圖６．５

水面曲線的形狀

６．４水位一流量關(guān)系
文獻［９７】根據(jù)實驗結(jié)果，利用圖解法給出溢流壩的水位一流量關(guān)系曲線計算公式

若＝（魯）ｌ

６，該實驗曲線可以作為數(shù)值計算結(jié)果的驗證比較。

首先計算設計流量Ｑ＝ＣＬＨ。３”米３，秒，其中Ｃ＝ｏ．５５２４瓦?ｍ為設計流量系數(shù)，
其值根據(jù)實驗回歸曲線確定。Ｌ為有效壩寬，二維情況下取值為１。 為了確定溢流壩的水位與流量的關(guān)系，給定８種狀態(tài)下的不同流量，計算出對應

７

河海大學博士論文

的壩頂水頭高度，繪制出水位一流量關(guān)系曲線，并與實驗結(jié)果進行比較，分別列于表
６＾３和圖６ ６中，可以看出，本文的計算方法能夠比較準確地給出溢流壩的水位與流

量關(guān)系，從而一方面可以在溢流壩設計時準確計算出設計流量，另一方面，能夠在全
部溢流水頭范圍內(nèi)，對理論結(jié)果或模型實驗值進行校核。
表６．３水位一流量關(guān)系 Ｑ｜Ｑｄ
Ｈ。／Ｈｄ
（計算）
Ｏ，０５ Ｏ．１ ０．２ ０．４ ０．６ Ｏ．８ １．Ｏ １．３３

ｏ １５３１７３

ｏ．２２９７５９

０．３６１０５

０ ５７９８６９

０．７１１１６

０．８５３３９２

１．００４３７６

１

１８１６１９

Ｈ：／Ｈ
Ｉ

ｄ

０ １５３７６５

０ ２３７１３７

０．３６５７１６

０ ５６４０１１

０．７２６６８２

０ ８６９８２４

１．０

１

１９５１０８

（實驗） ｌ相對誤差
０ ０００５９２ ０ ００７３７８ ０．００４６６６

０．０１５８６

０．０１５５２２

０ ０１５４３２

０．００４３８

０ ０１３４８９

ｌ

４

１

２

ｌ

０

８

ｌ＇＠ａ

０

６

０

４

Ｏ

２

Ｏ

圖６．６水位．流量關(guān)系曲線

６．５流量系數(shù)
流量系數(shù)是決定溢流壩泄洪流量的重要參數(shù)，對于同一壩面曲線，不同的壩頂水 頭將具有不同的流量系數(shù)，為了便于比較，本文在給定流量Ｑ／Ｑｄ＝Ｏ．３３，１．ｏ，１．２９， １．５８，１．９１時計算出對應的實際壩頂水頭分別為

Ｈｅ／Ｈｄ＝０．５０８７，１．００４４，１．１６６，１．３１３，１．４７８，根據(jù)公式ｍ２了南即可得到對應的流量

河海大學博士論文

系數(shù)，計算和實驗結(jié)果分別在表６．４和６．５中給出，同時將結(jié)果繪于圖６．７中，二者
基本一致。
表６．４流量系數(shù)（計算）

Ｑ｝Ｑｄ
Ｈ：｝Ｈ
卅
ｄ

Ｏ．３３

１．Ｏ

１．２９

１．５８

１．９１

０．５０８７

１．００４４

１．１６６

１．３１３

１．４７８

０．４５５６４７

０．４９７９５８

Ｏ．５１１５８５

０．５２５８３２

０．５３３６０４

表６ ５流量系數(shù)（實驗）

Ｑ舊ｄ
Ⅳ，／ＨＪ

０．３３

１．Ｏ

１．２９

１．５８

１．９１

Ｏ．５０

１．０

１．１７

１．３３

１．５

門，

０．４５５

０．５０ｌ

Ｏ．５１３

０．５２６

０．５３６

１．６

１

４

１

２

口

要
士

１

０．８

０．６

０．４

０．２
０．５ ｍ ０．５５

０．６

圖６．７流量系數(shù)

需要說明的是，關(guān)于壩頂水頭高度的計算是在水庫上游面按流體體積分數(shù)Ｆ＝Ｏ．５ 時對應的水面線高度進行的。
６．６壩面壓力 壩面壓力的計算是溢流壩水力設計的重要內(nèi)容，為此，本文分別計算了 日。／Ｈ。＝Ｏ．５，１．Ｏ，１．１７，１．３３，１．５共５種情況下的壩面壓力，繪制成曲線并與文獻【１７】中的

實驗點進行比較，見圖６．８，壓力分布的基本趨勢是一致的，其中工程中關(guān)心的負壓

河海大學博士論文

極值，其出現(xiàn)的位置與大小均與實驗符合良好。
０ ５
ｏ ４

ｏ ３

。弋
；＼ ， Ｉ
ｆ Ｌ ，ｊ

、

０ ２ ０ ｌ

＼． ｊ
√．
，
Ｉ

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

昭

０ ｌ

要

一

２
３ ４ ５ ６ Ｔ

Ｉ １
＼。

』
，，一

上
、

ｆ一
。、

—，一Ｊ

≯’

＼＿．
ｌ ｌ

＜?
＼

、 、 ／
＼／
－０
１

“’“‘叮。ｒ

叼∞加吣加加叼叼

８ －Ｏ．了
Ｕ

１

Ｕ．７３

Ｕ ｂ

０

Ｔ

０ ９

２７Ｈｄ

圖６．８壩面壓力分布曲線（圓點為實驗值）

６．７壩面水流速度分布
為考察溢流壩過壩水流的流速變化規(guī)律，在圖６．９中分別給出了斷面 ｘ＝０．０，２．５７６，５，１０，２０，３０，壩頂水頭分別為Ｈ。／Ｈｄ＝Ｏ．５，１．０，１．５時速度分布曲線圖

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

８５
、‘

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

０
’

盈

ｋ

，

—ｔｒ�！唬兀剑玻� 一
一

７０ 一 占６５
）、

一】
６０ ５５

５０
● ４５

４０
Ｏ Ｏ
Ｖ

冊

巧

∞

惦弘∽

河海大學博士論文

９０

９０

８５

８５

８０

７５

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

８０

佰 陽 ∞ ∞ 踮

７０

６５

６０

］
５７６ １０ １５

５５

５０

４５

一
０ ５

印

Ｘｌ０２ ?－?－－ｘ＝

，

蠣 蚰 ％

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

７６

４０

———■ｒ—一ｘ。２０

３５ ２０ ２ ５ ３０

０

５

１０

ｌ ５

２０

２５

３０

Ｖ（ｍ／ｓ）

Ｖ（ｍ／ｓ）

（ｂ）

（ｃ）

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

為進一步考察水流沿壩面的速度變化，給出了個斷面上最大速度的沿程變化曲線 見圖６．１０。從圖５．９，６．１０可以看出，水流速度沿壩面逐漸增加，且不同水頭下的最
大速度的差值沿程減小。為考察流線的整體變化情況，在圖６．１ｌ中給出了 Ｈ。／Ｈ。＝１．０時整個流場內(nèi)的流線圖。

圖６’ＩＪ

６．８本章小結(jié)

鎰流壩流線矢量圖

河海大學博士論文

第七章帶閘墩溢流壩三維過壩水流數(shù)值模擬
帶閘墩溢流壩三維過壩水流數(shù)值模擬的工作開展不多，主要原因計算區(qū)域增大， 閘墩區(qū)域形狀復雜，并且包含自由液面問題，這些都使得數(shù)值計算較為困難。因此過 壩水流的研究工作更多地集中在二維流場計算。盡管二維過壩流場計算具有計算工作 量小、計算成本低，并能夠在一定程度上指導工程設計，比如文獻『１１６１通過計算過 壩水流與閘墩繞流兩個二維模型的方法近似模擬了三維帶閘墩的溢流壩流場，可以作 為解決三維過壩水流問題的一條途徑，但是將二維數(shù)值模型應用于實際的帶閘墩的溢 流壩過壩水流計算還有一定的差距。主要表現(xiàn)在閘墩的存在對水面線位置、壩面壓力、 溢流流量等均產(chǎn)生顯著的影響。文獻ｉ１１７］利用有限體積方法對溢流壩三維湍流場進 行數(shù)值模擬，但沒有考慮閘墩的影響，本文直接建立溢流壩面與閘墩的三維整體模型， 對三維湍流場進行數(shù)值求解，其中湍流模型仍然采用ＲＮＧｋ—ｓ兩方程模型，自由水 面線的位置采用ＶＯＦ方法確定，針對不同的墩型和布置形式進行計算，比較準確地
給出了過壩水流的水面線位置及壩面壓力。 ７．１閘墩型式

溢流壩壩面曲線選用三圓�。祝牛忧�ｆ¨８】ｚ‘”＝２．０Ｈ≯Ｙ，見圖６，ｌ，６．２。 溢流壩設計的重要參數(shù)是流量計算，對流量影響最大的是閘墩，因此必須對包含 閘墩的溢流壩三維流場進行數(shù)值模擬。閘墩的作用是將溢流壩前緣分隔為若干個孔
口，并支承閘門、啟閉機和橋梁等傳來的載荷，閘墩的斷面形狀應盡量減小孔口水流 的側(cè)收縮，使水流能夠平順通過，因此關(guān)鍵在于墩頭和墩尾形狀的選擇，頭部主要影 響側(cè)收縮，從而影響泄洪流量，尾部主要影響下游流態(tài)，如可能產(chǎn)生水冠或沖擊波，
影響下游的消能。閘墩墩頭的形狀及閘墩的主要尺寸在圖７．１中給出。
Ｋ
／

，
乇 州 ｏ

Ｊ

半徑０ １３３Ｈｊ

、

Ｏ ２ｆ
● —、‘

７只、
ｒ

壩項軸線Ｌ
Ｉ

ｊ§ 習 磊 薄麗彰 悟
０．２ｔ

７Ｈａ．Ｊ
７Ｉ
—

ｆ．０２Ｉ 『片。。Ｉ壩頂軸線｝．０２ｆ
’ｉ

７Ｈ．，、Ｉ
’

—』乙?

Ｃ、／、

ｌ

ｌ

’＿ｒ

ｌ

ｉ¨

２型
圈７．１閘墩的型式

３Ａ型４型

河海大學博士論文

７．２計算區(qū)域與網(wǎng)格劃分 含閘墩溢流壩三維流場計算區(qū)域在圖７．２中給出，可以看出在閘墩區(qū)計算區(qū)域的 幾何形狀最為復雜，具有雙向曲率變化，這也正是研究該類問題的難點所在。為了保 證計算精度，對該區(qū)域進行網(wǎng)格劃分需要特別注意，本文采用分區(qū)網(wǎng)格系統(tǒng)，將整個 計算區(qū)域劃分為３個子區(qū)，見圖７ｔ３。由于閘墩區(qū)域邊界變形較大，采用結(jié)構(gòu)網(wǎng)格無 法生成可用于計算的網(wǎng)格系統(tǒng)，因此在該區(qū)內(nèi)采用非結(jié)構(gòu)網(wǎng)格，同時在壁面處網(wǎng)格進 行加密：在Ｉ區(qū)和ＩＩ區(qū)采用結(jié)構(gòu)網(wǎng)格，由于ＩＩ區(qū)處于自由表面區(qū)，為提高自由面的模 擬精度，同樣進行了網(wǎng)格加密。此外。為了保證Ｉ區(qū)、ＩＩ區(qū)與閘墩區(qū)網(wǎng)格的平順過渡， 在閘墩區(qū)的上游和左下方增加一部分區(qū)域歸并到閘墩區(qū)構(gòu)成ＩＵ區(qū)進行網(wǎng)格劃分，三個 區(qū)的網(wǎng)格數(shù)目約為２７０００，１４０００，７００００，共計約１ｌ萬個網(wǎng)格節(jié)點。在圖７．４中給出 了壩面及閘墩區(qū)的網(wǎng)格圖。

．１００

圖７．２含閘墩溢漣壩流場計算區(qū)域

河海大學博士論文

圖７．３

含閘墩溢流壩流場網(wǎng)格的分區(qū)劃分

圖７．４壩面及閘墩區(qū)的網(wǎng)格劃分

河海大學博士論文

７．３計算初始條件與邊界條件

壩高Ｈ＝７６．２米，設計水頭Ｈ。＝９．１４米，設計流量傷＝ＣＬ蟛３”，Ｃ為設計流
量系數(shù)，上為壩的凈寬＝壩寬．閘墩厚度，分別對２型墩和３Ａ型墩進行計算，其中墩 長均為１．７７Ｈｄ，墩厚分別為０．２６７Ｈｄ和０．１７３Ｈｄ，墩厚與閘寬之比均為０．２。假定壩 項水頭高度及流量，分別對三種不同水頭高度Ⅳ。／Ｈｄ＝０．５，１．Ｏ，１－３３下的自由水面線、 壩面壓力、流量系數(shù)進行計算，計算初值見下表
表７．１水頭高度及流量初值 Ｈ�。龋�
（初值）

Ｏ．５０

１．Ｏ

１．３３

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

０．３２９８７７

１．Ｏ

１．５７８２０２

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

及ｋ和￡值為未知．計算根據(jù)與出口斷面相鄰的上游兩個斷面的變量值，線性外延求 得出口斷面之值，作為下一次迭代的邊界條件．計算過程中，每次迭代都需根據(jù)上一 次迭代計算結(jié)果修正出口邊界條件，固體壁面上流速按非滑移條件給定，固壁附近的 ｋ和占采用用壁函數(shù)技術(shù)；自由表面上的壓力取為大氣壓，而ｋ和￡采用零梯度條件， 即令自由表面網(wǎng)格中的ｋ和ｓ的值等于相鄰的流體內(nèi)部網(wǎng)格的ｋ和ｓ的值。計算開始 流動為非定常流，經(jīng)過一段時間的計算，流場各水力要素均趨不變，流動趨于定常， 計算即告結(jié)束。 ７．４自由水面線 按照上節(jié)的計算條件分別對２型墩和３Ａ型墩水面線形狀及位置進行計算，
７，４．１

２型墩水面線計算 圖７．５中給出了沿２型墩閘孔中心線水面線位置計算與實驗結(jié)果的比較，三條曲

線分別對應日／Ｈ。＝０．５，１．Ｏ，１．３３三種情況。二者是比較一致的，在壩頂水頭較低時， 誤差較大，隨著壩頂水頭的增加，計算與實驗結(jié)果的一致性更好。與沿閘孔中心線的 計算相比，沿閘墩表面水面線位置計算結(jié)果與實驗值更加接近（見圖７．６），說明本計 算模型是準確、可靠的。

８０

河海大學博士論文

９０

８５

８０

— ５

ｈ ７０

７５

６５

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

ｘ（ｍ）

圖７．５沿２型墩閘孔中心線水面線位置比較
（從上至下分別對應Ｈｅ／Ｈｄ＝ｌ
９０
３３．１ ０．０．５）

８５

８０

．、

縣７５
＞、

７０

６５

６０

１０

—５

０

５

１０

１５

２０

２５

３０

ｘ（Ⅲ）

圖７．６沿２型墩表面水面線位置比較
（從上至下分別對應Ｈｅ／Ｈｄ＝Ｉ
３３，１．０，０ ５）

為了比較閘墩對水面線位置的影響，在圖７．７和圖７．８中分別給出了無閘墩、沿 閘孔中心線和沿閘墩表面的自由水面線位置的計算曲線，顯然閘墩的存在抬高了自由 水面線的位置，且隨著水位的增加，抬高的幅度逐漸增加。 而沿閘墩表面的變化則

是閘墩頭部水面抬升，在閘墩的平行段水面低于閘孔中心線處的水面。 在圖７．９中比較直觀地給出了三種水頭情況下通過２型墩閘孔的過壩水流。圖

河海大學博士論文

７．１０中給出了不同高度處閘墩附近的流線矢量圖。

９０

８５

８０

５
ｈ

７５

６５

６０ ｌＯ 一５ ０ ５

ｘ（ｍ３０

１５

２０

２５

３０

圖７．７無閘墩與２型墩沿閘孔中心線的自由水面線
（從上至下分別對應Ｈｅ／Ｈｄ＝ｌ
３３，１

０，０．５）

８５一

占７５，
＾

７０

６０�！�
一１０ ０ ５ １０ １５ ２０ ２５

ｘ（面 圖７．８沿２型墩閘孔中心線與沿閘墩表面的自由水面線
（從上至下分別對應Ｈｅ／Ｈｄ＝Ｉ
３３，１

０，０．５）

河海大學博士論文

（ａ）

（ｂ） 圖７．９ Ｈ／Ｈｄ＝０．５，１．Ｏ，１．３３時２型墩過壩水流 （紅色為水，黃色為閘墩和壩面）

（ｃ）

圖７．１０日／Ｈ。＝１．０時２型墩過壩水流不同高度處的流線矢量圖 （ａ）ｈ＝７７ｍ，（ｂ）ｈ＝７８ｍ，（ｃ）ｈ＝７９ｍ，（ｄ）ｈ＝８０ｍ
７．４．２

３Ａ型墩水面線計算 圖７．１ｌ，７．１２中給出了設計水頭Ｈ／Ｈｄ＝Ｏ．５，１．０，１．３３三種情況下沿３Ａ型墩閘孔

中心線和沿閘墩的水面線位置計算結(jié)果，并與無閘墩情況進行了比較。與圖７．７，７．８ 比較可以看出，兩種墩型對孔中心處水面線位置影響不大，為了比較３Ａ型墩與２型 墩對水面線位置的影響，在圖７．１３中給出了沿兩種墩型表面的水面線位置計算結(jié)果， ２型墩使墩前水面產(chǎn)生更高的抬升。

河海大學博士論文

９０

８５

８０

３
ｈ

７５

７０

６５

６０ ５ ０ ５
ｘ

，１０

１

５

２０

２５

３０

ＬｍＪ

圖７．１１無閘墩與３Ａ型墩沿閘孔中心線的自由水面線
（從上至下分別對應Ｈｅ／Ｈｄ＝１．３３，１
９０ ０，０ ５）

８５

８０

５
）、

７５

７０

６５

６０ １０ —５ ０ ５

ｘ（拶

１５

２０

２５

３０

圖７．１２沿３Ａ型墩閘孔中心線與沿閘墩表面的自由水面線
（從上至下分別對應Ｈｅ／Ｈｄ＝ｌ
３３，１

０．０．５）

河海大學博士論文

９０

８５

８０

占７５
ｈ

７０

６５

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

Ｘ（ｍ）

酗７．１３沿２型墩及３Ａ型墩表面的自由水面線
（從上至下分別對應Ｈｅ，／Ｈｄ＝Ｉ．３３，１ ０，０．５）

７．５壩面壓力 壩面壓力計算，特別是壩面最小壓力的數(shù)值及出現(xiàn)位置是溢流壩設計的重要參 數(shù)，在圖７．１４和圖７．１５中給出了三種不同水頭下２型閘墩沿閘孔中心線及沿閘墩的 壩面壓力分布曲線，在圖７．１６和圖７．１７給出了三種不同水頭下３Ａ型閘墩沿閘孔中

心線及沿閘墩的壩面壓力分布曲線，并與實驗結(jié)果進行對比，可以看出，沿閘孔中心

線壩面壓力的計算結(jié)果與實驗值非常一致，只是在以／ｚ‘＝１．０時在壩頂后段計算值
偏低，而壩面最小壓力值的大小及出現(xiàn)位置均與實驗吻合：對于沿閘墩面的壓力，總

的趨勢是計算結(jié)果略低于實驗值，并且隨著水頭的增加其偏差量增大。在 Ｈ。／Ｈ。，＝１．Ｏ、１．３３兩種情況下沿閘墩的壩面上的最小負壓出現(xiàn)位置與實驗一致，均在
壩頂前０．１Ｈ。處。

將２型墩與３Ａ型墩的壩面壓力計算結(jié)果進行比較可以發(fā)現(xiàn)，在相同的墩厚閘寬

比情況下，２型墩將使壩面出現(xiàn)更低的負壓，當水頭Ｈ。／Ｈ。＝１．３３時負壓極值為７．３ 米水柱，超過了工程中建議的不低于一６．１米水柱以保證壩面不出現(xiàn)空蝕現(xiàn)象的要求，

而３Ａ型墩在相同水頭下對應的最低負壓僅為一４．６水拄，不會出現(xiàn)壩面空蝕現(xiàn)象。

河海大學博士論文

０．６ ０，４ ０，２

≈

鼉０
￡
－０，２
—０，４

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

Ｏ，３

Ｏ，５

０，７

０，９

１，１

１，３

Ｘ／Ｈｄ
圖７．１４沿２型墩閘孔中心線上的壩面壓力分布
Ｏ

６

０．４

０

２

＼
．Ｃ

Ｚ

０
２

一０

—０ ４
—０ ６

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

Ｘ／Ｈｄ 圖７．１５沿２型閘墩表面的壩面壓力分布
０．６ Ｏ．４ ０．２ 弓

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

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

－０ ３—０．１

０．５

０．７

０．９

１．１

１．３

Ｘ，，Ｉ－Ｉｄ

圖７．１６沿３Ａ型墩閘孔中心線上的壩面壓力分布比較
（從上至下分別對應Ｈｅ．，／ｌ－Ｉｄ＝０ ５，Ｉ．０，ｌ ３３）

河海大學博士論文

圖７．１７沿３Ａ型閘墩表面的壩面壓力分布比較
（從上至下分別對應ＨｅＪＨｄ＝０．５，１．０，ｌ ３３）

７．６流量系數(shù) 流量系數(shù)關(guān)系到工程的泄洪規(guī)模和具體布置，但流量系數(shù)隨工程的體型而變化， 由于水流均為三維，其流態(tài)和邊界條件都比較復雜，故在一般的手冊中難以找到準確 的數(shù)據(jù)。在圖７．１８中比較了墩厚閘寬比ｔ／ｂ＝０．２時２型墩和３Ａ型墩在三種不同水頭 下的流量系數(shù)，顯然，３Ａ型墩具有更大的流量系數(shù)。
ｍ ６

５

呲ｌＩＪ
ｎ ５

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

１．２

１．４

Ｈｅ／ｌｉｄ

圖７．１８有閘墩影響的流量系數(shù)（ｒ／ｂ＝０．２）

７．７墩厚閘寬比對壩面壓力及流量系數(shù)的影響
除壩型和閘墩型式對泄流流量、壩面壓力、水面線位置產(chǎn)生較大影響外，閘墩在 壩面上的布置同樣至關(guān)重要，特別是對壩面壓力影響更為顯著，因此為了考察閘墩布

置對壩面壓力等的影響，針對３Ａ型閘墩，給定初值ＨｅＰｒ－Ｉｄ＝１．０時情況，計算了墩厚

ｆ與閘寬ｂ之比ｔ／ｂ卸．１４，０．２，Ｏ．５三種情況下壩面壓力分布，對應的實際水頭高度

河海大學博士論文

分別為０．９６８２７１，１．００７６５９，１．０１７５０５，分別給出了沿ｉＹｑ孑ｔ，ｅｅ心線和沿閘墩的壩面壓 力分布曲線（圖７，１９，７．２０），并給出了對應的流量系數(shù)曲線（圖７．２１）。
０，６
０，４

０，２
－。

暮０
上

一０，２
—０，４

－０，６
—０，８ —１

—０．３—０，１

０，ｌ

０．３

０，５

０，７

０，９

ｉ，１

１，３

Ｘ／Ｈｄ 圖７．１９

Ｈｅ／Ｈｄ＝１．０時不同墩厚閘寬比沿閘孔中心線的壓力分布

０，６

０，４ ０，２
口

罩０
上 一０，２

．撥
一¨＿Ｖ

ｒ．＿?�。翕殻� 垃 －１ｌｔ－１一一：

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

—０。４ —０，６

—０，８
。 ． ：

一１
一０．３—０．１

；

；

０，１

０，３

０，５

０，７

０．９

１，１

ｌ，３

Ｘ／Ｈｄ

圖７．２０

Ｈｅ／Ｈｄ＝１．０時不同墩厚閘寬比沿閘墩面的壓力分布

８８

河海大學博士論文

０，５５

０。５３

０，５１

ｇ
０ ４９

０，４７

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

ｔ／ｂ 圖７．２１

Ｈｅ／Ｈｄ＝１．０時不同墩厚閘寬比對應的流量系數(shù)

７．８本章小結(jié) 本章在二維過壩水流流場計算的基礎上，直接建立溢流壩和閘墩的三維整體模

型，對三維過壩水流粘流場進行數(shù)值計算，其中采用了分區(qū)網(wǎng)格系統(tǒng)來處理復雜的幾
何邊界，針對２型墩和３Ａ型墩計算了閘墩型式對水面線形狀、壩面壓力、流量系數(shù)

等的影響，通過與實驗結(jié)果的對比表明，本文建立的模型能夠比較準確地給出過壩水 流的水面線位置及壩面壓力，計算指出，閘墩的存在抬高了水面線的位置，在閘墩頭
部尤其明顯，３Ａ型閘墩相對于２型墩具有更大的流量系數(shù)和更小的壩面壓力。此外， 不同的墩厚閘寬比對泄流能力也將產(chǎn)生顯著的影響，結(jié)果表明，隨著墩厚閘寬比的增 加，壩面壓力降低，而當ｔ／ｂ＝０．２時溢流壩具有更大的流量系數(shù)。

河海大學博士論文

第八章結(jié)論與展望
本文對粘流場數(shù)值模擬技術(shù)的研究現(xiàn)狀和發(fā)展動態(tài)作了比較全面的回顧和展望， 重點研究了網(wǎng)格生成技術(shù)、數(shù)值求解技術(shù)、湍流模型技術(shù)和動邊界模擬技術(shù)，在此基
礎上建立了模擬自由面湍流場的數(shù)學模型，并成功地應用到二維溢流壩和帶閘墩的三

維溢流壩過壩水流場數(shù)值計算。

８．１本文完成的主要工作與結(jié)論 １）對網(wǎng)格生成技術(shù)進行了較深入的研究，建立了一種簡單的、能夠獨立于網(wǎng)格生成 技術(shù)之外的邊界正交化方法，給出了詳細的數(shù)學模型并編制了計算程序；通過網(wǎng)格生 成算例表明，該方法能夠大大改善網(wǎng)格在邊界處的正交性；同時針對分區(qū)結(jié)構(gòu)網(wǎng)格系 統(tǒng)建立了分區(qū)交界面處的數(shù)據(jù)結(jié)構(gòu)和計算模型： ２）利用有限體積方法在非正交同位網(wǎng)格系統(tǒng)中建立了ＳＩＭＰＬＥ求解算法，并對非正 交同位網(wǎng)格系統(tǒng)中的邊界條件、延遲修正技術(shù)及計算節(jié)點的梯度計算等專題內(nèi)容進行
了深入討論。通過對傾斜方腔驅(qū)動流、后臺階繞流、非對稱平板間的圓柱繞流等典型

流場的數(shù)值模擬，表明本文建立的非正交同位網(wǎng)格系統(tǒng)中的ＳＩＭＰＬＥ算法是合理的，
具有更加廣泛的適用范圍：

３）研究比較了兩種湍流模型的特點及應用范圍，以水翼粘流場為例完成了ＲＮＧｋ一￡ 湍流模型與標準ｋ一￡模型的計算比較，表明ＲＮＧｋ一占模型的模擬結(jié)果更優(yōu)； ４）研究了運動邊界模擬技術(shù)中的ＶＯＦ方法，詳細建立了一種利用斜直線構(gòu)造自由面
的高精度方法，通過對典型運動界面的數(shù)值模擬表明，該方法能夠比其它自由面重構(gòu)
方法更加準確地確定運動界面的形狀和位置：

５）完成了對二維溢流壩過壩水流進行計算，分析了不同壩頂水頭下的水面線位置、

壩面壓力及流量、流速等水力特性，通過與典型實驗資料的對比，計算的水面線高度
及壩面負壓極值均與實驗值具有非常好的一致性。

６）直接建立溢流壩和閘墩的三維整體模型，對三維過壩水流粘流場進行數(shù)值計算， 針對閘墩的型式及在壩面的布置計算了閘墩對水面線形狀、壩面壓力、流量系數(shù)等設 計參數(shù)的影響，并將不同墩型與布置形式的計算結(jié)果進行比較，結(jié)果表明，閘墩的存
在抬高了水面線的位置，其中在閘墩頭部尤其明顯；墩型對流量系數(shù)和壩面壓力影響 較大，３Ａ型閘墩相對于２型墩具有更大的流量系數(shù)和更小的壩面壓力；不同的墩厚 閘寬比對泄流能力也將產(chǎn)生顯著的影響，隨著墩厚閘寬比的增加，壩面壓力降低，而

河海大學博士論文

當ｆ／ｂ＝０．２時溢流壩具有更大的流量系數(shù)。

７）本文建立的三維自由面湍流場模型能夠比較準確地模擬帶閘墩溢流壩三維過壩水 流湍流場，可以針對不同的壩型、墩型、壩面布置形式以及不同的壩頂水頭完成系列
水力計算，為泄水工程建筑物的設計提供可靠的分析依據(jù)。

本論文的主要創(chuàng)新點： １）提出并建立了網(wǎng)格邊界正交化算法模型：
２）建立了非正交同位網(wǎng)格系統(tǒng)中的ＳＩＭＰＬＥ算法，并以此為基礎建立了非對接分區(qū) 結(jié)構(gòu)網(wǎng)格數(shù)值求解模型； ３）提出了一種精度更高的自由面重構(gòu)算法模型，通過典型運動界面的模擬驗證了本 算法是精確的：

４）建立了帶閘墩溢流壩三維整體計算模型，首次完成了帶閘墩溢流壩三維過壩水流
粘流場數(shù)值模擬，得出了一些有意義的結(jié)論。

８．２展望 雖然本文對帶自由面湍流場數(shù)值模擬技術(shù)開展了一定的研究工作，并成功地應 用到泄水工程中三維自由面湍流場的數(shù)值模擬，得到了一些有意義的結(jié)論。但是由 于三維自由面湍流是一個十分復雜的流體運動現(xiàn)象，其復雜的運動規(guī)律和內(nèi)在機理
需要作進一步的分析和研究，因此，本論文只是在幾個方面做了一些探討性的研究， 今后仍有大量的工作要做，具體內(nèi)容包括： １）進一步完善非正交同位網(wǎng)格系統(tǒng)中的ＳＩＭＰＬＥ算法，使之適用于網(wǎng)格強非正交的
情形；

２）對基于ＶＯＦ方法的運動界面追蹤技術(shù)中的輸運模型和重構(gòu)方法仍需作更深入的研

究，進一步提高計算精度，以便更好地模擬水流的破碎、空化等現(xiàn)象； ３）開展帶閘墩溢流壩上下游流場整體模型數(shù)值計算，探討泄洪流量、空化空蝕等復 雜現(xiàn)象的內(nèi)在機理。

９ｌ

河海大學博士論文

致

謝

作者在攻讀博士論文的過程中得到了許多的指導、支持和幫助，在此表示衷心的
感謝。

首先感謝我的導師汪德燧教授，正是在導師一步一步的悉心指導和督促下，才得
以順利完成畢業(yè)論文，導師深厚的學術(shù)造詣、嚴謹?shù)闹螌W態(tài)度、踏實的工作作風給我 留下了深刻的印象并將使我受益終生，在論文即將完成之際，特向?qū)煴硎咀钫\摯的 謝意和最美好的祝愿。

感謝華東船舶工業(yè)學院副校長王自力教授和科技處長蔣志勇教授，他們在我的
工作中給予了非常大的幫助，使我克服了工作中的許多困難，為順利完成論文提供
了良好的工作和學習條件。

感謝楊松林教授、朱仁慶副教授、姚震球副教授在論文工作中提供的無私幫助
和建議。

感謝師弟賴錫軍，在與他日常的討論交流中作者獲得了許多有益的幫助和啟

迪，他的聰穎與豁達也給我留下了深刻的印象。感謝同學童朝峰、陳雄波，在與他 們的交往過程中獲得了許多幫助。
在論文撰寫過程中，得到了研究生周林慧的幫助，在此表示感謝。 還要感謝我的家人，他們幫助我分擔了許多家務，提供了大力的支持，正是在

他們的幫助和鼓勵下，我才得以將更多的時間、精力和飽滿的熱情投入到論文工作
中：感謝母親對我的關(guān)心和牽掛，感謝并深深懷念我的父親。 最后再一次向關(guān)心、支持和幫助過我的所有同仁表示最誠摯的謝意。

河海大學博士論文

參考文獻
１．Ｔｈｏｍｐｓｏｎ Ｊ．Ｆ．，Ｗｅａｔｈｅｒｉｌｌ Ｎ．Ｐ，Ａｓｐｅｃｔｓ ｏｆ ｎｕｍｅｒｉｃａｌ ｇｒｉｄ ｇｅｎｅｒｔｉｏｎ：ｃｕｒｒｅｎｔ ａｎｄ ａｒｔ．

Ａ１ＡＡ（ｐａｐｅｒ），９３—３５３９
２．Ｓｈｉｈ Ｔ．Ｉ—Ｅａｎｄ Ｂａｉｌｅｙ Ｒ．Ｔ，Ａｌｇｅｂｒａｉｃ ｇｒｉｄ ｇｅｎｅｒａｔｉｏｎ ｆｏｒ ｃｏｍｐｌｅｘ ｇｅｏｍｅｔｒｉｅｓ，

Ｉｎｔ．Ｊ．Ｎｕｍ

  本文關(guān)鍵詞：三維自由面湍流場數(shù)值模擬及其在水利工程中的應用，由筆耕文化傳播整理發(fā)布。