李昆鵬，魏成柱，梁曉鋒

(1. 陸軍研究院作戰(zhàn)保障研究所，江蘇 無錫 214035；2. 中國艦船研究設(shè)計(jì)中心上海分部，上海 201108；3. 上海交通大學(xué) 海洋智能裝備與系統(tǒng)教育部重點(diǎn)實(shí)驗(yàn)室，海洋工程國家重點(diǎn)實(shí)驗(yàn)室，上海 200240)

0 引 言

滑行艇在軍事和民用領(lǐng)域的應(yīng)用非常廣泛，其相關(guān)研究工作也在大量進(jìn)行。傳統(tǒng)上，滑行艇的設(shè)計(jì)和優(yōu)化工作基于經(jīng)驗(yàn)公式和模型試驗(yàn)。隨著CFD 技術(shù)的發(fā)展，模型試驗(yàn)結(jié)合CFD 成為滑行艇設(shè)計(jì)優(yōu)化和性能分析的主流手段之一。CFD 技術(shù)在滑行艇水動(dòng)力性能分析中扮演了愈來愈重要的角色。相比潛艇和排水型船，高速滑行艇的數(shù)值仿真難度更大。Brizzolara等[1]和王碩等[2]分別進(jìn)行了滑行艇數(shù)值仿真計(jì)算精度的研究，均使用了NACA 公布的滑行楔形體模型。Brizzolara 等在數(shù)值仿真計(jì)算中使用了六面體網(wǎng)格，王碩等在數(shù)值仿真計(jì)算中使用了六面體網(wǎng)格和切割體網(wǎng)格。他們的研究工作表明基于CFD 的滑行艇數(shù)值仿真精度很高，誤差在實(shí)驗(yàn)測(cè)量誤差之內(nèi)，表明基于CFD 的滑行艇數(shù)值仿真計(jì)算可以應(yīng)用于實(shí)際工程中。

網(wǎng)格是CFD 分析的基礎(chǔ)?，F(xiàn)實(shí)中的滑行艇通常具有復(fù)雜幾何形狀的船體表面，經(jīng)典的六面體網(wǎng)格很難對(duì)其進(jìn)行網(wǎng)格劃分。因此，采用非結(jié)構(gòu)性網(wǎng)格成為解決此類網(wǎng)格貼面問題的重要手段。傳統(tǒng)上，四面體網(wǎng)格被用來適配一些滑行艇的復(fù)雜幾何表面，但近年來切割體網(wǎng)格（trimmed mesh 或cut cell）開始得到比較廣泛的應(yīng)用。除了非結(jié)構(gòu)性的四面體網(wǎng)格和切割體網(wǎng)格外，多面體網(wǎng)格也逐漸得到廣泛應(yīng)用。多面體網(wǎng)格已經(jīng)在空氣動(dòng)力學(xué)相關(guān)領(lǐng)域有了廣泛的應(yīng)用并贏得良好的口碑[3]。但是多面體網(wǎng)格在海洋工程領(lǐng)域中的應(yīng)用十分有限。多面體網(wǎng)格是一種非結(jié)構(gòu)網(wǎng)格。二維多面體網(wǎng)格主要由六邊形網(wǎng)格構(gòu)成，其幾何形式廣泛存在于自然界中；三維多面體網(wǎng)格單元為任意多面體形狀。多面體網(wǎng)格在引入可靠的基于面的有限體積法和相關(guān)的多面體網(wǎng)格生成器后才逐漸可用，并在近幾年得到了商業(yè)解算器的推廣。相較于同為非結(jié)構(gòu)網(wǎng)格的四面體網(wǎng)格，多面體網(wǎng)格有更多的臨近單元，能夠更加精確地計(jì)算梯度及局部的流動(dòng)狀況，并且對(duì)幾何變形沒有四面體網(wǎng)格敏感。對(duì)于一種新興的并具有良好特性的網(wǎng)格，有必要驗(yàn)證其在海洋工程領(lǐng)域中的適用性。

研究人員在對(duì)滑行艇進(jìn)行數(shù)值計(jì)算時(shí)經(jīng)常會(huì)遇到偽擴(kuò)散的問題。 一層薄空氣層會(huì)附在船體，與實(shí)際船體表面水氣分布不符合。若使用VOF 法[4]捕捉自由面，偽擴(kuò)散意味著船體浸沒表面部分或全部的水組分不等于1。這個(gè)空氣層會(huì)直接影響船體表面的切應(yīng)力。偽擴(kuò)散解方程歐拉型模式所特有的，其大小與所用的有限差分格式有關(guān)。為了克服偽擴(kuò)散，須采取特殊的技術(shù)措施和各種不同的差分格式。在網(wǎng)格設(shè)置方面，克服偽擴(kuò)散的方法有增加網(wǎng)格密度，沿流線方向布置網(wǎng)格。王志剛等[5]指出，網(wǎng)格加密10 倍偽擴(kuò)散效應(yīng)就減小1 個(gè)量級(jí)。每個(gè)方向按10 倍進(jìn)行加密，那么3 個(gè)方向的網(wǎng)格數(shù)將是原網(wǎng)格數(shù)的1 000 倍，需要付出極大的計(jì)算代價(jià)。故采用網(wǎng)格加密不能有效減小偽擴(kuò)散效應(yīng)，而沿流線方向布置網(wǎng)格總體上僅對(duì)靜態(tài)仿真有效，當(dāng)網(wǎng)格發(fā)生改變時(shí)原有優(yōu)化布置的網(wǎng)格很可能不再與流線方向相符合。此外，F(xiàn)RISK[6]在研究中通過添加源匯，強(qiáng)制性將船體濕表面處水組分小于0.7 的地方設(shè)置為1。Spiegel 等[7]對(duì)比分析了四面體網(wǎng)格和多面體網(wǎng)格在腦血管動(dòng)力學(xué)仿真計(jì)算中的應(yīng)用效果，認(rèn)為多面體網(wǎng)格收斂性更好，且對(duì)WSS（wall shear stress）的計(jì)算更準(zhǔn)確，這對(duì)計(jì)算滑行艇的摩擦阻力是有益的。參考文獻(xiàn)[7]并考慮到多面體網(wǎng)格有更多的臨近單元，能夠更加精確地計(jì)算梯度及局部的流動(dòng)狀況，本文嘗試通過使用多面體網(wǎng)格處理滑行艇數(shù)值仿真中遇到的偽發(fā)散問題。

研究證實(shí)，基于多面體網(wǎng)格滑行艇的水動(dòng)力預(yù)報(bào)具有相當(dāng)高的精度，完全滿足工程需要，多面體對(duì)控制滑行艇仿真計(jì)算中出現(xiàn)的偽擴(kuò)散問題有很大改善。

1 多面體網(wǎng)格計(jì)算精度分析

使用NACA 公開的滑行楔形體幾何模型及相關(guān)的實(shí)驗(yàn)設(shè)置和實(shí)驗(yàn)數(shù)據(jù)[8]?；行ㄐ误w的幾何尺寸參見文獻(xiàn)[8]，幾何模型如圖1 所示。選擇底部斜升角為20°的模型。仿真計(jì)算采用固定模式，即模型初始狀態(tài)設(shè)置為水池實(shí)驗(yàn)中測(cè)量到的狀態(tài)。選擇其中幾個(gè)典型的工況進(jìn)行仿真計(jì)算，如表1 所示。實(shí)驗(yàn)中測(cè)得的升力平均相對(duì)誤差為3%，最大相對(duì)誤差為8%；測(cè)得的阻力平均相對(duì)誤差為20%，最大相對(duì)誤差為50%；速度平均相對(duì)誤差為0.70%，最大相對(duì)誤差為1.30%；吃水平均相對(duì)誤差為5%，最大相對(duì)誤差為10%；縱傾角平均相對(duì)誤差為2%，最大相對(duì)誤差為5%。

圖 1 NACA 滑行楔形體幾何輪廓Fig. 1 NACA planing wedge

表 1 基于NACA 滑行楔形體的數(shù)值計(jì)算工況統(tǒng)計(jì)Tab. 1 Statistics on simulations of NACA planing wedge

數(shù)值仿真計(jì)算使用重疊網(wǎng)格，背景網(wǎng)格使用切割體網(wǎng)格控制網(wǎng)格密度，減小網(wǎng)格數(shù)量和提高網(wǎng)格的分層性和正交性。隨體網(wǎng)格通過合理布置獲得良好分層性和正交性的多面體網(wǎng)格。合理布置網(wǎng)格密度及改善網(wǎng)格的正交性和分層性對(duì)自由面鋸齒現(xiàn)象的消除及興波與飛濺的捕捉非常關(guān)鍵。計(jì)算域和網(wǎng)格設(shè)置如圖2所示。由于使用重疊網(wǎng)格，網(wǎng)格數(shù)量在260 萬左右。

計(jì)算選用STAR CCM 作為求解器來求解非穩(wěn)態(tài)RANS方程。數(shù)值計(jì)算選用SST k-ω[9]湍流模型。為了對(duì)自由面的捕捉，使用VOF（volume of fluid）模型。為了減小計(jì)算代價(jià)并參考文獻(xiàn)[10]，使用高yplus（又稱y+）壁面函數(shù)和較為粗糙的近壁面網(wǎng)格在一定程度上減少網(wǎng)格數(shù)量來提高計(jì)算速度。

圖 2 基于NACA 滑行楔形體的計(jì)算域設(shè)置及網(wǎng)格Fig. 2 Computational domain settings and mesh

數(shù)值仿真計(jì)算給出了光順的興波圖，在傳統(tǒng)非結(jié)構(gòu)網(wǎng)格在波形描繪中出現(xiàn)的等高線鋸齒現(xiàn)象有很大改善，如圖3 所示。將本文的計(jì)算結(jié)果同水池實(shí)驗(yàn)結(jié)果及文獻(xiàn)[1] 和文獻(xiàn)[2] 中的阻力和升力結(jié)果進(jìn)行匯總（見表2 和表3）。表4 和表5 進(jìn)一步匯總了本文及文獻(xiàn)[1]和文獻(xiàn)[2]中的計(jì)算誤差。文獻(xiàn)[1]基于六面體網(wǎng)格的計(jì)算結(jié)果標(biāo)記為HEX-1；基于六面體網(wǎng)格的計(jì)算結(jié)果標(biāo)記為HEX-2，文獻(xiàn)[2]基于切割體網(wǎng)格的計(jì)算結(jié)果標(biāo)記為TRIMMED；本文的計(jì)算結(jié)果標(biāo)記為POLY。對(duì)比文獻(xiàn)[8]中NACA 給出的水池一般性固有誤差可見，基于多面體網(wǎng)格的數(shù)值仿真計(jì)算的升力誤差與其3%固有平均誤差相比稍微偏大，但誤差均小于最大相對(duì)誤差，而阻力計(jì)算誤差遠(yuǎn)小于其20%的參考值。同其他計(jì)算結(jié)果對(duì)比表明，基于多面體網(wǎng)格的阻力預(yù)報(bào)的精度稍高，而升力的預(yù)報(bào)的誤差稍微偏大。對(duì)比每一個(gè)工況下不同網(wǎng)格的阻力和升力預(yù)報(bào)可知，基于多面體網(wǎng)格的預(yù)報(bào)結(jié)果誤差值和分布具有一定的相似性。

通過以上對(duì)比可見，基于多面體網(wǎng)格的滑行艇水動(dòng)力預(yù)報(bào)具有相當(dāng)高的精度，完全滿足工程需要。

2 多面體網(wǎng)格對(duì)偽擴(kuò)散的改善

圖 3 NACA 滑行楔形體興波（case2）Fig. 3 Wave making of NACA planing wedge (case2)

表 2 不同網(wǎng)格的NACA 滑行楔形體阻力計(jì)算結(jié)果Tab. 2 Numerical results of drag of NACA planing wedge based on different meshes

表 3 不同網(wǎng)格的NACA 滑行楔形體升力計(jì)算結(jié)果Tab. 3 Numerical results of lift of NACA planing wedge based on different meshes

表 4 不同網(wǎng)格的NACA 滑行楔形體阻力計(jì)算結(jié)果誤差Tab. 4 Drag error of NACA planing wedge based on different meshes

在基于切割體網(wǎng)格、嘗試使用STAR CCM 作為求解器對(duì)某一細(xì)長尖肶滑行艇進(jìn)行數(shù)值仿真計(jì)算時(shí)遇到了偽擴(kuò)散問題，表現(xiàn)為船體濕表面處的水組分值應(yīng)為1 的地方卻小于1。但在使用多面體網(wǎng)格（填充隨體區(qū)域）和六面體網(wǎng)格（填充靜止區(qū)域）構(gòu)成的混合網(wǎng)格并利用FLuent 求解時(shí)未遇到偽擴(kuò)散問題，船底濕表面處水組分值為1，如圖4 所示。在研究過程中發(fā)現(xiàn)，通過合理設(shè)置第一層網(wǎng)格的高度和選擇合適的壁面函數(shù)可以改善基于切割體網(wǎng)格的偽擴(kuò)散問題，但不能完全解決。為了驗(yàn)證多面體網(wǎng)格改善偽擴(kuò)散的能力，對(duì)比不同yplus 下，F(xiàn)n=1.63 時(shí)分別基于切割體網(wǎng)格和多面體網(wǎng)格的船體表面水組分分布及阻力預(yù)報(bào)結(jié)果。網(wǎng)格采用上一節(jié)中的設(shè)置，區(qū)別在于隨體區(qū)域分別由切割體網(wǎng)格填充和由多面體網(wǎng)格填充。船體固定，姿態(tài)設(shè)置為基于四面體網(wǎng)格給出的預(yù)報(bào)值。計(jì)算基于求解非穩(wěn)態(tài)RANS 方程。

表 5 不同網(wǎng)格的NACA 滑行楔形體升力計(jì)算結(jié)果誤差Tab. 5 Lift error of NACA planing wedge based on different meshes

圖 4 基于多面體網(wǎng)格的船體表面水組分分布Fig. 4 Volume fraction of water on the hull bottom predicted by polyhedral mesh

由圖5 可知，不同yplus 下切割體網(wǎng)格對(duì)偽擴(kuò)散的敏感性要明顯高于多面體網(wǎng)格。對(duì)比工況a 和工況b，可見在壁面函數(shù)選擇和yplus 值分布區(qū)間不合理時(shí)，基于多面體網(wǎng)格的數(shù)值仿真計(jì)算仍然給出了合理的船底水組分分布，而基于切割體網(wǎng)格的數(shù)值仿真計(jì)算則出現(xiàn)了嚴(yán)重的偽擴(kuò)散。水組分分布直接影響了摩擦阻力的預(yù)報(bào)。由表6 可知，工況a 所得的摩擦阻力系數(shù)遠(yuǎn)小于其他工況所得的摩擦阻力系數(shù)。當(dāng)選合理選擇壁面函數(shù)和yplus 值分布區(qū)間時(shí)，基于切割體網(wǎng)格的數(shù)值仿真計(jì)算遇到的偽擴(kuò)散問題有很大改善，但是沒有徹底消除。而基于多面體網(wǎng)格的數(shù)值仿真計(jì)算未遇到偽擴(kuò)散的問題，對(duì)偽擴(kuò)散不敏感。因此在不能很好確定流場(chǎng)狀況及第一層網(wǎng)格高度的情況下，使用多面體網(wǎng)格可以減小預(yù)報(bào)誤差。此外，本文研究結(jié)果也同文獻(xiàn)[7]中提到的多面體網(wǎng)格對(duì)WSS（wall shear stress）的計(jì)算更準(zhǔn)確形成呼應(yīng)。

圖 5 不同yplus 下的基于多面體網(wǎng)格和切割體網(wǎng)格船底面水組分分布Fig. 5 Volume fraction of water on the hull bottom predicted by polyhedral and trimmed meshes under different yplus

表 6 不同工況下網(wǎng)格參數(shù)及阻力系數(shù)Tab. 6 Mesh information and drag coefficient of different cases

3 結(jié) 語

本文研究得到以下結(jié)論：

1）基于多面體網(wǎng)格的滑行艇數(shù)值仿真計(jì)算具有相當(dāng)高的精度，完全滿足工程需要；

2）多面體可以改善滑行艇仿真計(jì)算中遇到的偽擴(kuò)散問題。通過對(duì)比切割體網(wǎng)格可知，基于多面網(wǎng)格的滑行艇的數(shù)值計(jì)算對(duì)偽擴(kuò)散問題不敏感。對(duì)于一些具有復(fù)雜外形設(shè)計(jì)、難以使用六面體網(wǎng)格計(jì)算域填充的滑行艇仿真案例，若使用切割網(wǎng)格遇到偽擴(kuò)散問題時(shí)則可以考慮使用多面體網(wǎng)格，或者直接使用多面體網(wǎng)格填充隨體區(qū)域來嘗試解決遇到的偽擴(kuò)散問題。