鄔 明，孫善春

（中國(guó)船舶重工集團(tuán)公司第七一○研究所，湖北 宜昌 443003）

0 引 言

超空泡技術(shù)是一種可以使水下航行體高效減阻技術(shù)，其原理是在水與航行體表面形成穩(wěn)定的氣層，達(dá)到減小水的粘性阻力，從而提高水下航行體的運(yùn)行速度。水下超空泡航行體所采用的操縱面，要兼顧全沾濕狀態(tài)和超空泡狀態(tài)時(shí)的流體動(dòng)力布局要求：航行體處于全沾濕狀態(tài)時(shí)，航行器流體動(dòng)力布局與常規(guī)水下航行體相似，需依靠操縱面提供穩(wěn)定力矩和操縱力矩。而在超空泡狀態(tài)時(shí)要盡量減小其所受阻力，同時(shí)操縱面必須穿透空泡壁進(jìn)入水中才能提供所需力矩，該狀態(tài)下的操縱面應(yīng)具有較高的升阻比，因此對(duì)于操縱面的研究具有重要意義。

隨著計(jì)算性能的快速提高和計(jì)算數(shù)學(xué)理論的不斷發(fā)展完善，數(shù)值方法已經(jīng)成為研究問題和解決問題的重要手段之一，它不僅在理論研究領(lǐng)域得到普遍應(yīng)用，而且在工程實(shí)際中也被廣泛使用。求解雷諾平均N－S方程式是當(dāng)前數(shù)值計(jì)算的主要方法，需補(bǔ)充湍流模型對(duì)方程進(jìn)行封閉。湍流模型對(duì)于數(shù)值模擬空化流場(chǎng)的精度具有決定性影響。當(dāng)前應(yīng)用最普遍的湍流模型是k－ε 湍流模型［1－2］，然而該模型有以下不足之處：在模型中湍流尺度未知；僅限于湍流邊界層壓力相對(duì)穩(wěn)定的情況；其壁面函數(shù)在邊界層的修正中難以彌補(bǔ)計(jì)算模型與實(shí)際現(xiàn)象之間的差距。相對(duì)于k－ε湍流模型的以上不足之處，SSTk－ω湍流模型［3］卻擁有以下優(yōu)點(diǎn)：該模型能夠適應(yīng)逆壓梯度變化的各種物理現(xiàn)象；可應(yīng)用于粘性內(nèi)層，通過對(duì)壁面函數(shù)的應(yīng)用，能夠精確模擬邊界層的現(xiàn)象，無需使用較容易失真的粘性衰減函數(shù)。本文將SSTk－ω湍流模型應(yīng)用于操縱面空化流場(chǎng)的仿真分析中，取得了較好的結(jié)果。

SSTk－ω湍流模型以下簡(jiǎn)稱SST湍流模型。

1 控制方程與算法

1．1 連續(xù)方程

式中：ρ表示密度；U表示速度；t為時(shí)間。

其中

μt和αl分別表示流體動(dòng)力粘度和水相的體積分?jǐn)?shù)。

1．2 湍流模型

若分別以φ1、φ2、φ3表示k－ω 模型、k－ε 模型和SST湍流模型中的函數(shù)關(guān)系，則SST湍流模型可表示為：

SST湍流模型考慮到湍流剪切應(yīng)力的輸運(yùn)，不但能夠?qū)Ω鞣N來流進(jìn)行準(zhǔn)確的預(yù)測(cè)，還能在各種壓力梯度下精確模擬分離現(xiàn)象，其綜合了近壁面k－ω模型的穩(wěn)定性及邊界層外部k－ε模型獨(dú)立性的優(yōu)點(diǎn)，它的計(jì)算模擬性能優(yōu)于后兩者。各系數(shù)取值為β′＝0．09，α1＝5／9，β1＝0．075，σk1＝2，σω1＝2，α2＝0．44，β2＝0．0828，σk2＝1，ωω2＝1．168，各個(gè)數(shù)據(jù)的取值取自參考文獻(xiàn)［4］。

2 空化模型

選擇Singhal空化模型［5］，其中 Rayleigh Plesset方程是控制汽化產(chǎn)生和空泡凝結(jié)的理論基礎(chǔ)，該方程描述了液體中空泡的變化過程：

式中r為空泡半徑，Pc為汽化壓強(qiáng)，σvl是汽液兩相間的表面張力系數(shù)，質(zhì)量傳輸率為：

式中，ρg表示飽和蒸汽密度；sgn表示符號(hào)函數(shù)；NB代表單位體積的汽泡數(shù)；F代表經(jīng)驗(yàn)因子，汽化過程和凝結(jié)過程取值不同，F(xiàn)vap＝50，F(xiàn)cond＝0．01。

3 計(jì)算模型

3．1 計(jì)算區(qū)域和網(wǎng)格劃分

針對(duì)水下高速航行體在不同運(yùn)動(dòng)階段的流體動(dòng)力特性，對(duì)構(gòu)成操縱面的二維基本翼型進(jìn)行了相應(yīng)研究——選用典型的小頂角楔形，經(jīng)驗(yàn)表明以上區(qū)域邊界對(duì)流場(chǎng)的計(jì)算結(jié)果影響很小。

計(jì)算網(wǎng)格的好壞直接影響到數(shù)值計(jì)算的可行性、收斂性以及計(jì)算精度。前處理軟件ICEM CFD是一款成熟的網(wǎng)格劃分軟件，它向用戶提供業(yè)界領(lǐng)先的高質(zhì)量網(wǎng)格技術(shù)，其強(qiáng)大的網(wǎng)格劃分功能可以滿足CFD仿真計(jì)算的嚴(yán)格要求。

對(duì)于本文中的問題若采用非結(jié)構(gòu)網(wǎng)格，生成的空泡邊界會(huì)出現(xiàn)凹凸不平的毛刺，影響對(duì)仿真結(jié)果的分析。若想消除毛刺現(xiàn)象，就要對(duì)網(wǎng)格進(jìn)行細(xì)分。因此采用基于六面體結(jié)構(gòu)化網(wǎng)格劃分方法，對(duì)計(jì)算區(qū)域進(jìn)行全六面體結(jié)構(gòu)網(wǎng)格劃分，并利用與湍流強(qiáng)度相關(guān)的Yplus對(duì)所建模型進(jìn)行考核。楔形操縱面流場(chǎng)網(wǎng)格圖及局部放大圖如圖1和圖2所示。

圖1 計(jì)算區(qū)域與網(wǎng)格劃分Fig．1 The computational domain and meshing

圖2 操縱面局部放大圖Fig．2 Zooming in images locally in control surface

3．2 邊界條件設(shè)定

邊界條件設(shè)定：上游邊界給定來流速度和各流體體積分?jǐn)?shù)；下游邊界給定壓強(qiáng)；流域外邊界為光滑物面邊界。初始條件：設(shè)定初始來流速度和各流體體積分?jǐn)?shù)。

3．3 數(shù)值方法

采用有限元的有限體積法（FVM），N－S方程中的對(duì)流項(xiàng)采用CFX高精度格式，它利用調(diào)整混合因子提高計(jì)算精度且平均收斂精度。以一階迎風(fēng)為例，對(duì)流項(xiàng)：

式中：φi表示變量φ第i個(gè)節(jié)點(diǎn)值；φu為變量φ迎風(fēng)節(jié)點(diǎn)值；ξ為混合因子；▽?duì)諡橛L(fēng)節(jié)點(diǎn)的節(jié)點(diǎn)梯度；從迎風(fēng)節(jié)點(diǎn)到第i個(gè)節(jié)點(diǎn)的矢量。

針對(duì)不同條件下，ξ會(huì)作出相應(yīng)的調(diào)整來滿足計(jì)算要求獲得精確的解，例如在第i節(jié)點(diǎn)臨近迎風(fēng)節(jié)點(diǎn)的節(jié)點(diǎn)梯度變化劇烈，則混合因子趨近于0，反之趨近于1。

4 楔形操縱面的典型空泡

根據(jù)相關(guān)實(shí)驗(yàn)文獻(xiàn)［6］，楔形操縱面的典型空泡形態(tài)分為三種，如圖3。

圖3 楔形操縱面的典型空化形態(tài)（全空化、部分空化和底部空化）Fig．3 Typical characteristics of cavities of control surface（entirely、substantially and bottom enveloped by cavities）

從圖3中不難發(fā)現(xiàn)，當(dāng)來流攻角小于楔形半頂角時(shí)，僅在翼型底部產(chǎn)生空化現(xiàn)象，而在上下翼面完全處于全沾濕狀態(tài)。此時(shí)所產(chǎn)生流體動(dòng)力和力矩與攻角近似成線性關(guān)系，且較為穩(wěn)定。當(dāng)來流攻角大于楔形半頂角時(shí)，上翼面處于超空泡狀態(tài)，而下翼面為全沾濕狀態(tài)，此時(shí)，同樣可產(chǎn)生穩(wěn)定的流體動(dòng)力和力矩，且與攻角成線性關(guān)系。當(dāng)來流攻角與楔形半頂角相當(dāng)時(shí)，除在翼型底部產(chǎn)生空化外，上翼面還將出現(xiàn)間歇性片狀空泡，此時(shí)翼型所產(chǎn)生的流體動(dòng)力及力矩呈非定常脈動(dòng)分布。

5 仿真計(jì)算結(jié)果及分析

空化數(shù)是衡量流體空化程度的物理量，計(jì)算公式為：

其中P∞、PV、ρ、V∞分別為某基準(zhǔn)位置上的絕對(duì)靜壓力、空泡內(nèi)的絕對(duì)靜壓力、流體密度、相對(duì)于物體的未擾動(dòng)流體的速度。研究中設(shè)定初始條件P∞、25℃時(shí)水的飽和蒸汽壓、V∞和ρ。

L為翼型弦長(zhǎng)，升、阻力系數(shù)基于速度坐標(biāo)系。

表1給出了繞楔形操縱面的超空化仿真與試驗(yàn)結(jié)果的對(duì)比。

表1 繞楔形操縱面的超空化仿真與試驗(yàn)結(jié)果Table 1 Results of simulation and test of control surface

圖4給出了計(jì)算得到的阻力系數(shù)和升力系數(shù)隨攻角的變化仿真曲線及實(shí)驗(yàn)曲線，兩者基本一致，從中可以看出當(dāng)保持自然空泡數(shù)不變，僅改變來流攻角時(shí)，操縱面的阻力系數(shù)呈逐漸增大趨勢(shì)，而升力系數(shù)在攻角增至一定數(shù)值后，其增大幅度減緩。隨著攻角的增大，空泡不斷增大，并形成超空泡，使得操縱面全部處于空泡的包裹之中，此時(shí)，操縱面的控制效率也將相應(yīng)地降低。圖5中的升阻比曲線分布表明，當(dāng)攻角增至一定數(shù)值后，升阻比將有所減小，這與臨界攻角和局部空泡的形成有關(guān)。

圖4 阻力系數(shù)和升力系數(shù)隨攻角的變化曲線Fig．4 The change of drag characteristics and lift characteristics by angles of attack

圖5 升阻比隨攻角的變化曲線Fig．5 The change of ratio of lift over drag by angles of attack

CFX后處理可以更加直觀地觀察空泡形態(tài)，如圖6所示，不同攻角下，操縱面的密度云圖。從圖中可以看出，當(dāng)來流攻角很小時(shí)（α＜3°），楔形操縱面底部產(chǎn)生空泡；當(dāng)攻角較小時(shí)（3°≤α≤3．5°），在操縱面上端面產(chǎn)生局部空泡；而在攻角繼續(xù)增大時(shí)（α≥3．5°），操縱面完全被空泡包裹，并且呈增大趨勢(shì)。這與理論結(jié)果可較好地吻合。

圖6 不同攻角下的密度分布圖Fig．6 The contour of the charge density in different angles of attack

圖7 給出的是空化數(shù)為0．323時(shí)，操縱面的壓力系數(shù)分布，可看出，其下翼面上壓力系數(shù)分布呈單調(diào)遞減趨勢(shì)，當(dāng)攻角a＝3°時(shí)，上翼面出現(xiàn)局部空泡的壓力系數(shù)分布形態(tài)，增大來流攻角后，上翼面的壓力系數(shù)將降至一恒定值，即上翼面已處于超空化狀態(tài)。

圖7 不同攻角，操縱面表面的壓力分布曲線Fig．7 The curve of pressure distribution for control surface in different angles of attack

圖8 給出了不同空化數(shù)下的空泡輪廓圖。從圖中可以發(fā)現(xiàn)隨著空化數(shù)的降低空泡尾部的回射流強(qiáng)度在減弱，這是因?yàn)榭张蓍L(zhǎng)細(xì)比增大，閉合時(shí)空泡邊界入射角減小的緣故。這和文獻(xiàn)［7］中的理論分析取得了很好地一致性。

圖8 不同空化數(shù)下的空泡輪廓Fig．8 The profile of cavities in different citation numbers

6 結(jié) 論

采用SST湍流模型，文章對(duì)典型二維操縱面在不同來流空化數(shù)及攻角條件下的升阻力特性、空泡形態(tài)、表面壓力分布規(guī)律進(jìn)行了詳細(xì)的研究，并與相關(guān)實(shí)驗(yàn)做了分析比較，結(jié)果表明：

（1）減小空化數(shù)后，相同攻角條件下的操縱面的阻力系數(shù)和升力系數(shù)都有所下降。隨著攻角的增大，由于上翼面處于超空化狀態(tài)，其升阻比曲線將有所減小；

（2）運(yùn)用N－S方程及SST湍流模型對(duì)均質(zhì)流場(chǎng)求解，并利用Rayleigh－Plesset空泡模型對(duì)二維操縱面的空化流動(dòng)進(jìn)行數(shù)值模擬，與相關(guān)實(shí)驗(yàn)相對(duì)比，取得了較好的一致性，證明了所選湍流模型的可行性，為以后的超空泡實(shí)驗(yàn)研究提供參考依據(jù)；

（3）SST湍流模型能夠?qū)Σ倏v面的空化流場(chǎng)進(jìn)行較為準(zhǔn)確的模擬，尤其在空泡流場(chǎng)、空泡外形以及尾部回射流方面，可以準(zhǔn)確地預(yù)測(cè)空化流場(chǎng)，從而指導(dǎo)超空泡航行器的外形設(shè)計(jì)，顯示出了優(yōu)越性，具有廣闊的應(yīng)用前景。

［1］ KUNZ R，BOGER D，STINEBRING D，et al．A preconditioned implicit method for two－phase flows with application to cavitation prediction［J］．Comput．Fluids，2000，29（8）：849－875．

［2］ LEROUX J B，COUTIER－DELGOSHA O，ASTO－LFI J A．A joint experiment and numerical study of mechanisms associated to instability of partial cavitationon two－dimensionalhydrofoil［J］．Phys．Fluids，2005，17（5）：052－101．

［3］ COUTIER－DELGOSHA O．Numerical prediction of cavitation flow on a two－dimensional symmetrical hydrofoil and comparison to experiments［J］．Journal of Fluids Engineering，2007，129（3）：279－291．

［4］ KUBOTA A，KATO H．A new modeling of cavitating flows：a numerical study of unsteady cavitation on a hydrofoil section［J］．FluidMech．，1992，240（1）：59－96．

［5］ SINGHAL A K，ATHAVALE M M，et al．Mathematical basis and validation of the full cavitation model［J］．ASME Journal of Fluids Engineering，2002，124：617－624．

［6］ GARABEDIAN P R．Calculation of axially symmetric cavities and jets［J］．Pac．J．Math．，1956，6（4）：611－684．

［7］ SEMENENKOV N．Artificial supercavitation．physics and calculation［A］．RTO AVT lecture series on supercavitating flows［C］．Von Karman Institute Brussels，Belgium：2001，206237．

［8］ 張宇文．空化理論與應(yīng)用［M］．西安：西北工業(yè)大學(xué)出版社，2007．