劉祥，袁壽其，朱涵，司喬瑞，邱寧

(江蘇大學國家水泵及系統(tǒng)工程技術(shù)研究中心，江蘇 鎮(zhèn)江212013)

空化是流體機械中一種常見的物理現(xiàn)象，普遍認為是在高流速流場中由于液體內(nèi)部局部壓力低于當?shù)仫柡驼羝麎毫Χ鴮е乱合嘟橘|(zhì)發(fā)生相變[1].空化現(xiàn)象廣泛存在于各種流體機械中，例如泵、噴嘴、推進器、水翼等，對流體機械運行性能有顯著的影響[2].空化的發(fā)生常伴有水體阻力增大、劇烈相間變化、振動噪聲及動力損失等現(xiàn)象，空化嚴重時可能引起劇烈壓力脈動和材料表面空蝕[3].近年來，由于流體機械不斷向高轉(zhuǎn)速設(shè)計，空化的發(fā)生不可避免[4].因此，研究空化非定常的流動特征，對流體機械的設(shè)計與改進有著至關(guān)重要的作用.

近年來，隨著計算流體動力學發(fā)展和計算機性能的提高，數(shù)值模擬成為人們認識和探索復(fù)雜空化現(xiàn)象與空化流動機理的基本方法.針對水翼繞流問題，早期的研究主要分析空化的發(fā)生區(qū)域與形態(tài)[5]，對空化的形成機理以及相應(yīng)響應(yīng)特征沒有統(tǒng)一的定論.云空化或片空化向云空化的轉(zhuǎn)捩過程常被認為是回射流引起的，即由液體組成的回射流形成于空泡的閉合區(qū)域，并且沿著空腔下方的壁面向上運動，導致空泡脫落形成周期性云空化[6].GANESH等[7]發(fā)現(xiàn)由于片狀空化向云狀空化轉(zhuǎn)捩過程的不穩(wěn)定性，會在壁面形成非定常壓力脈動.GNANASKANDAN等[8]通過數(shù)值模擬對非定?？栈@流進行了研究，結(jié)果表明名義進口空化數(shù)(normalized inlet cavitation number)σ/2α較高時，回射流對空化流動影響占主導地位；σ/2α較低時，空化云在翼型尾部脫落、潰滅后產(chǎn)生的高幅值壓力波是片狀空化從翼型表面脫落的主要原因.SAITO等[9]通過有限體積法模擬了固定側(cè)壁之間的NACA0015水翼周圍三維非定?？栈魈匦裕l(fā)現(xiàn)兩側(cè)側(cè)壁上的邊界層由層流向湍流的不斷發(fā)展與回射流不斷向水翼前緣靠近是形成U型云空化的主要原因.而在以Navier-Stokes方程為基礎(chǔ)的空化全流場數(shù)值模擬中，湍流模型的選擇對高雷諾數(shù)空化可壓縮性流動的精確求解有著重要的影響.為了捕捉大尺度的周期性云空化脫落與小尺度下的流動特征，JI等[10]采用大渦模擬(large eddy simulation, LES)方法對NACA0015水翼片空化與云空化進行了精確求解，但該方法要求對精細的網(wǎng)格劃分以及對時間步長的限制使其應(yīng)用范圍有相當?shù)木窒扌裕y以應(yīng)用在工程實際.傳統(tǒng)的湍流模型無法合理預(yù)測可壓縮性空化流動，需要對湍流模型進行修正[11].研究表明，湍流渦黏性系數(shù)對于空化流動中空泡的形成與壓力脈動的求解起決定性因素[12]，為解決標準雷諾時均化湍流模型(Reynolds-averaged Navier-Stokes, RANS)對湍流黏度過度預(yù)測，JOHANSEN等[13]提出的濾波模型(filter-based model, FBM)與COUTIER等[14]提出的密度修正模型(density corrected model, DCM)能在一定程度上解決此問題.HUANG等[15]采用密度函數(shù)將FBM與DCM的優(yōu)點相結(jié)合，提出的一種基于密度分域的湍流模型(FBDCM)，在使用較少的計算資源情況下，既能充分模擬近翼型近壁面空化核心區(qū)域的氣液混合相可壓縮性，又能充分模擬遠壁面大尺度云空化的非定常演變，使得此方法能夠精確求解空化流動中兩相之間的動態(tài)特征與空泡脫落頻率.

文中針對空化的非定常特性及其誘導產(chǎn)生的壓力脈動，采用ANSYS商用軟件CFX對NACA0015水翼空化繞流進行數(shù)值模擬.為精確求解空化可壓縮兩相流動，采用密度分域模型(FBDCM)對湍流黏性系數(shù)進行修正.此次研究主要分析云狀空化的脫落機理，及其脫落過程中水翼表面不同位置的壓力脈動特性.

1 控制方程與計算模型

1.1 控制方程

空化流動屬于氣-液兩相范疇，在流動過程中包含氣液相間介質(zhì)轉(zhuǎn)換，文中將兩相視為均質(zhì)氣泡混合物，連續(xù)性方程與動量方程表達式為

(1)

(2)

式中：下標i和j為坐標方向；u為速度；p為壓力；δij為雷諾切應(yīng)力；ρm為氣液混合相密度，定義為ρm=ρlαl+ρvαv；αv與αl為氣相體積分數(shù)與液相體積分數(shù)；μt為湍流黏度.

1.2 空化模型

空化模型用于求解氣液兩相之間質(zhì)量交換速率.為準確模擬空化非定常過程中流體密度改變與相變過程，均相流模型被普遍采用.根據(jù)對混合密度不同定義，均相流可分為狀態(tài)類模型與運輸方程類模型.文中采用的Zwart空化模型為通過推導R-P方程而建立氣液兩相之間運輸關(guān)系模型，其蒸發(fā)與凝結(jié)速率定義為

(3)

(4)

式中：αnuc為氣核體積分數(shù)，αnuc=1.0×10-4；RB為氣泡半徑，RB=1.0×10-6m；pv為飽和蒸汽壓力；Cdest與Cprod分別為空泡的蒸發(fā)與凝結(jié)速率,Cdest=50,Cprod=0.01[16].

許多研究表明，湍流對空化的流動特性有顯著影響，因此將飽和蒸汽壓力定義[17]為

pv=psat+0.5ptur,

(5)

ptur=0.39ρmk,

(6)

式中:psat為蒸汽壓力;k為湍動能.

1.3 湍流模型與可壓縮性修正

采用ANSYS商用軟件CFX中RNGk-ε湍流模型求解控制方程組.但CFX用于求解湍流的RANS方法是基于單相且不可壓縮流動而建立的，在模擬空化流動時會過度預(yù)測流動過程中湍流黏度，導致無法準確模擬云狀空化的脫落過程，也無法精確捕捉空泡脫落過程中的動力特征.混合密度分域湍流模型修正方法對遠離壁面區(qū)域使用FBM模型對空化旋渦脈動進行濾波修正，在近壁面區(qū)域使用DCM模型充分考慮氣液兩相流可壓縮性.湍流黏性系數(shù)修正表達式[18]為

(7)

fhybrid=ξ(ρm/ρl)fFBM+[1-ξ(ρm/ρl)]fDCM,

(8)

2tanhC1,

(9)

(10)

(11)

式中：C1，C2,C3和Cμ為計算模型的調(diào)節(jié)參數(shù)；k為湍動能；ε為湍流耗散率；λ為網(wǎng)格濾波尺寸；fDCM為密度修正函數(shù)表達式；fFBM為網(wǎng)格濾波函數(shù)表達式.

DCM模型和FBM模型通過修正湍流模型橋接函數(shù)ξ(ρm/ρl)聯(lián)立用于限制空化流動中湍流黏性系數(shù)μt-FBDCM.研究使用的量綱一化參數(shù)分別定義為：空化數(shù)σ=(p∞-pv)/(0.5ρlU2);壓力系數(shù)Cp=(p-p∞)/(0.5ρlU2);雷諾數(shù)Re=ρmU∞c/μ,其中U∞，p∞，ρl和p分別為計算流體域中進口斷面平均流速、出口斷面平均壓力、液體密度與監(jiān)測位置處壓力值；pv為飽和蒸汽壓力.

1.4 算例設(shè)置

為研究水翼繞流非定?？栈捌涿撀溥^程中的脈動機理，采用ANSYS進行模擬，通過CFX表達式語言輸入FBDCM湍流模型.為保證模擬準確性與穩(wěn)定性，計算模型的調(diào)節(jié)參數(shù)C1，C2和Cμ取值分別為4.00,0.20和0.09[15].試驗在浙江大學航空航天學院水洞中進行，試驗段長1 020 mm，截面尺寸200 mm×200 mm[19].模擬采用與試驗相同的NACA0015水翼空化流動條件，雷諾數(shù)Re=1.3×106,所選取翼型弦長為c=100 mm,展長為2c，布置在計算域中部，攻角為α=8°.圖1—3分別為計算域邊界條件、翼型表面Y+值分布和網(wǎng)格劃分示意圖.

圖1 計算區(qū)域及邊界條件

圖2 翼型表面Y+值分布

圖3 翼型周圍網(wǎng)格分布

計算域采用速度進口，取值U∞=10 m/s；壓力出口，其值根據(jù)空化數(shù)σ=1.25設(shè)為54 088.2 Pa；上下邊及翼型表面均為無滑移固壁，兩側(cè)為對稱壁面.為限制翼型表面附近對湍流渦黏性系數(shù)過度預(yù)測，網(wǎng)格濾波尺寸設(shè)為近壁面附近加密區(qū)域網(wǎng)格單元最大值[18]，即λ≥Δgird=(Δx·Δy·Δz)1/3，取值λ=1.3×10-5m.在非定常模擬中，時間步長設(shè)為Δt=5.0×10-5s，為使得平均的庫朗數(shù)CFL處在一個合理的區(qū)間，CFL=U∞Δt/Δx≤5，其中Δx定義為網(wǎng)格濾波尺寸特征長度.

2 結(jié)果分析

2.1 云空化非定常演變過程

數(shù)值模擬的監(jiān)測點布置與試驗一致，即監(jiān)測點布置在對稱軸線上，分別位于x/c=0.18，0.32，0.46，0.60弦長位置,其中x為翼型軸面弦向方向位置坐標.文中分析空泡產(chǎn)生與潰滅的周期性演變過程，以及對應(yīng)時刻壓力脈動的變化，數(shù)值模擬結(jié)果雖與試驗結(jié)果存在一定差異，但從宏觀空泡演變非定常過程仍然能較好地模擬空泡脫落過程的特征.

較薄附著型空泡首先在翼型前緣形成，發(fā)展過程中僅小尺度空泡脫落.當空泡接近翼型中部，新空泡不斷生成，生成的空泡與未潰滅的附著型空泡融合，空泡的體積分數(shù)不斷增加.試驗結(jié)果指出非定常云空化雖然具有高湍流度特征與非定常演變特性，但存在逆壓梯度，回射流附著空泡的剪切力導致空泡從翼型表面脫落，形成大尺度脫落型云狀空泡.當云狀空泡群完全脫離翼型尾緣后在下游高壓區(qū)域發(fā)生潰滅，如此周而復(fù)始形成完整周期.

對比數(shù)值模擬與試驗結(jié)果，云空化演化過程具有明顯周期性.圖4為一個典型云空化演化周期內(nèi)空泡變化過程，數(shù)值計算結(jié)果由空泡體積分布與空泡體積分數(shù)等值面分布共同呈現(xiàn).定義數(shù)值計算收斂曲線穩(wěn)定后，連續(xù)變化10個周期演變時間的平均值為Tc，圖4從上至下依次為t1,t2,…,t5時刻.

圖4 云空化演化過程空泡體積分數(shù)變化與等值面結(jié)構(gòu)變化

通過對比分析，將云空化演化過程分為3個階段：第一階段(t1=Tc+0.250Tc)—附著型空泡形成與生長；第二階段(t3=Tc+0.500Tc)—附著型空泡脫落與云空化形成；第三階段(t5=Tc+0.875Tc)—云狀空泡發(fā)展與潰滅.在第一階段，在翼型尾緣處上個周期的云空泡完全脫落并潰滅.在翼型前緣，由于局部壓力低，附著型空泡形成漸向翼型中部發(fā)展.t1時刻，附著型空泡已發(fā)展至水翼中部位置.在第二階段，附著型空泡開始脫落，水翼外側(cè)空泡移動速度高于近壁面空泡移動速度，形成順時針旋轉(zhuǎn)的云狀空泡群，且在翼型表面不斷生成的空泡與脫離的附著型空泡混合后形成空泡群.t4時刻，回射流在逆壓梯度的作用下由下游發(fā)展至翼型中部，與翼型中部的空泡相遇后致使空泡群發(fā)生劇烈變化.空泡群在回射流的作用下潰滅，并在翼型表面發(fā)生脫落.在第三階段，云狀空泡在回射流的作用下脫離翼型表面，隨著主流移動到下游高壓區(qū)域，發(fā)生潰滅.此次數(shù)值模擬的回射流機制下云狀空泡的脫落過程與之前的研究成果[20-21]一致.

2.2 云空化的周期特性與壓力脈動的關(guān)系

圖5為數(shù)值模擬4個監(jiān)測點在6個周期內(nèi)的壓力變化情況，可以看出，4個監(jiān)測點壓力系數(shù)的幅值變化均為-1.0～1.0，周期性變化的時間跨度均為毫秒量級.由于監(jiān)測點處于翼型表面不同位置，壓力系數(shù)在單個周期內(nèi)變化情況有所不同，總體周期性改變趨勢相同.根據(jù)空泡演變的整個非定常周期，空泡在翼型前緣生成，隨主流移動到下游后潰滅，而空泡在移動過程中會引起翼型表面壓力分布變化.因此，各個監(jiān)測點壓力信號有一定的相位差，且下游的壓力信號相對于上游的壓力信號具有一定滯后性.

圖5 數(shù)值模擬展向監(jiān)測點壓力系數(shù)變化過程

根據(jù)整個非定常云空化演變過程，附著型空泡在水翼前緣產(chǎn)生，發(fā)展至翼型中部后在水翼表面脫落，在下游高壓區(qū)域發(fā)生潰滅.在弦向方向，壓力系數(shù)的改變與云空化演化具有明顯時間傳遞特性，在一個新的周期內(nèi)，翼型前緣附著型空泡生成的初始時刻，下游整體的壓力水平高于翼型前緣附近區(qū)域，因此在逆壓梯度的作用下會形成回射流.在空泡發(fā)展過程中，空泡覆蓋翼型表面區(qū)域后，該區(qū)域壓力系數(shù)會迅速降低至較低水平，當空泡脫離該區(qū)域后，當?shù)貕毫ο禂?shù)水平又逐漸回升.云空化周期演變的反復(fù)過程導致了各個監(jiān)測點壓力系數(shù)也呈周期性的變化，且監(jiān)測點壓力系數(shù)變化過程能夠量化云空化的演化過程.

圖6為翼型對稱軸線上弦向方向壓力系數(shù)分布，可以看出，空泡形成初期，翼型中部的壓力系數(shù)高于翼型前緣區(qū)域，脫落并潰滅的空泡對壁面產(chǎn)生局部高壓.后處理提取翼型表面壓力脈動系數(shù)分布、空泡體積等值面分布，如圖7所示，圖中LE為水翼前緣(leading edge)，TE為水翼尾緣(trailing edge).結(jié)果顯示，空泡發(fā)展過程中，空泡覆蓋翼型表面區(qū)域時，該區(qū)域壓力系數(shù)迅速降低，空泡脫離該區(qū)域后，壓力系數(shù)水平又迅速增大.在t1～t2時刻空泡迅速向下游發(fā)展，但此時翼型尾部壓力系數(shù)水平較高，形成了較大逆壓梯度，導致空泡發(fā)展至中下游時從翼型表面脫落，這種現(xiàn)象與當前研究認為逆壓梯度產(chǎn)生的回射流致使空泡脫落的結(jié)論[22]一致.在附著型空泡形成初期t1時刻，前緣區(qū)域處于低壓狀態(tài)，而在翼型尾部上個周期的脫落型云空泡的大尺度潰滅，該區(qū)域整體處于高壓狀態(tài)，形成了較大的逆壓梯度.

圖6 弦向壓力系數(shù)瞬態(tài)值分布

圖7 翼型表面壓力系數(shù)分布俯視圖和等值面分布

結(jié)合圖5吸力面壓力系數(shù)云圖可知，在翼型壁面尾緣附近區(qū)域壓力水平明顯高于兩側(cè)靠近壁面區(qū)域，由此產(chǎn)生的近壁面回射流速度中部小于兩端，所以在整個周期內(nèi)兩側(cè)空泡脫落時間優(yōu)先于中間區(qū)域.因此，在展向方向上空泡尾部閉合區(qū)域呈曲線型，如圖中t2～t5時刻.而當回射流與主流在前緣相遇時，由于回射流主流方向相反，將導致壁面附近區(qū)域水的流速為0，當?shù)貏幽軒缀跬耆D(zhuǎn)變?yōu)閴耗?，致使該區(qū)域內(nèi)局部壓力系數(shù)水平升高，這一現(xiàn)象在文獻[23]中也有類似報道.

2.3 空泡體積分數(shù)與壓力脈動的關(guān)系

通過數(shù)值模擬得到的計算域空泡總體積分數(shù)變化規(guī)律如圖8所示，清晰地反映了云狀空泡體積分數(shù)周期變化過程中逆壓梯度隨時間改變過程以及空泡發(fā)展的3個階段.

圖8 空泡總體積分數(shù)變化規(guī)律

在整個變化周期過程中，發(fā)生空化的核心區(qū)域為低壓區(qū)域，而非空化區(qū)域為高壓區(qū)域.在空化演變的第一階段時刻A之后，附著型空泡在翼型前緣生成，非空化區(qū)域壓力遠高于附近發(fā)生空化的核心區(qū)域，同時逆壓梯度的產(chǎn)生回射流且不斷向上游發(fā)展.在第二階段臨近結(jié)束時，如時刻B所示，此周期內(nèi)計算域內(nèi)空泡體積分數(shù)發(fā)展至最大，且此時空泡幾乎覆蓋整個翼型表面，致使翼型表面壓力系數(shù)處于較低水平.在時刻B之后，空泡隨著主流移動向下游高壓區(qū)域，空泡體積分數(shù)逐漸減小，在空泡尾部閉合區(qū)域產(chǎn)生的回射流導致空泡脫落并且伴隨部分空泡潰滅.在時刻C之后，空泡從翼型表面脫落之后，翼型表面的壓力系數(shù)又逐漸恢復(fù)到較高水平.在一個周期完全結(jié)束之后，下一周期的附著型空泡在水翼前緣逐漸生成，而此周期內(nèi)的空泡完全脫離翼型表面并在下游高壓區(qū)發(fā)生部分潰滅，空泡體積分數(shù)又開始逐漸上升.因此，從翼型表面完整的非定?？栈葑冎芷诳矗芷谛援a(chǎn)生脫落與潰滅的空泡伴隨翼型表面壓力非定常脈動.

3 結(jié) 論

通過數(shù)值模擬與試驗相結(jié)合的方法，對繞水翼空泡脫落過程及其空化誘導產(chǎn)生的壓力脈動特性進行研究，得到如下結(jié)論：

1) FBDCM模型結(jié)合了DCM模型與FBM模型的優(yōu)點，能夠有效捕捉繞NACA0015水翼云空化演化中空泡的初生、生長、脫落、潰滅的整個過程以及演變過程中誘發(fā)的壓力脈動.數(shù)值模擬預(yù)測的演化過程與試驗具有較高的吻合度.

2) 通過數(shù)值模擬結(jié)果與試驗結(jié)果的對比分析，云空化演化可分為3個階段：①附著型空泡形成與生長;②附著型空泡脫落與云空化形成；③云狀空泡發(fā)展與潰滅.在第二階段，脫離翼型表面的空泡會引起大幅值壓力系數(shù)波動，由逆壓梯度導致的回射流是空泡從翼型表面脫落的主要原因.同時，附著型空泡通過充分發(fā)展，流域內(nèi)空泡總體積分數(shù)會增至整個周期內(nèi)最大，使得翼型表面壓力系數(shù)處于較低水平，在回射流與來流空泡相遇后，附著型空泡從表面脫落.在第三階段，空泡在下游高壓區(qū)域潰滅，空泡總體積分數(shù)降至最小，翼型表面壓力系數(shù)隨著空泡的脫落有所回升.

3) 翼型表面脫落型的云狀空泡尾部區(qū)域壓力系數(shù)的回升與兩側(cè)壁面效應(yīng)共同導致了展向方向上空泡非均勻性發(fā)展與較大幅值的壓力脈動，并且攜帶能量的回射流會使空泡從翼型表面脫落，脫落的空泡呈連續(xù)性狀態(tài).各個監(jiān)測點壓力信號有一定的相位差，且下游壓力信號相對于上游壓力信號具有一定滯后性.