超空化流動數值模擬

0空化器特性仿真研究超空氣泡法是一種高效、低可見水支持技術的方法。其原理是在水與航行體表面之間形成穩定的氣層,達到減小水的粘性阻力,從而提高水下航行體的運行速度。為了能夠生成穩定的、高效減阻的、形態可控的超空泡,位于航行體頭部的空化器起著決定性作用。空化器形狀直接影響超空泡形態及水下航行體減阻效果,因此優化空化器形狀是獲得穩定高效減阻超空泡的前提,其中圓盤空化器具有控制超空泡形態的特點。因此國內外的學者針對空化器作了大量實驗研究即理論分析工作Savchenko利用控制空化器參數來獲得穩定的超空泡,并總結出經大量實驗驗證過的經驗公式。Stinebring等對頭部為圓盤空化器進行了研究,并給出了空泡長度及厚度與空化器參數關系的公式。國內何友聲等對系列頭體的空化器作了實驗研究,在速度比較低即空泡數比較大的情況下,給出了一系列的經驗公式。賈力平、王海斌對空化器與阻力關系進行了研究。袁緒龍、張宇文應用商業軟件FLUENT建立多相流CFD模型對典型空化器進行了數值模擬。以往對空化器進行的數值模擬都是針對二維軸對稱體空化器進行研究。但是,當考慮重力影響以及帶攻角空化器生成的超空泡問題時,二維仿真不能解決此問題,因此本文針對三維圓盤空化器超空化流動進行數值模擬。首先,對在不同空化數下無攻角圓盤空化器生成超空泡進行仿真,將仿真結果同Savchenko經驗公式進行對比,總結其變化規律;其次,針對帶攻角圓盤空化器進行數值模擬,將模擬結果與無攻角圓盤空化器所得結果進行對比,總結其變化規律。文中主要是針對帶攻角圓盤空化器進行研究。為了解決三維數值模擬中龐大計算量問題和獲得更好的仿真結果,本文采用構架于深騰1800高性能服務器上的商業軟件CFX進行仿真。選擇深騰1800高性能服務器是利用其多個高速CPU以及大容量內存的優點解決龐大計算量問題,而擇商業軟件CFX是因為其高效的并行計算方法,以及可以獲得較精確仿真結果。1控制方程和算法運用N-S方程及k-ε湍流模型對均質流場求解,引入Rayleigh-Plesset方程模擬超空泡問題。1.1連續方程?ρ?t+??(ρU)=0,(1)式中:ρ表示密度;U表示速度;t表示時間。1.2流體動力sm廣義源項{??t(ρU)+??(ρU?U-μ(?U+(?U)Τ))=SΜ；ρ=ΝΡ∑a=1γaρa；μ=ΝΡ∑a=1γaμa?(2)式中:a是相的序號且1≤a≤Np(Np是相數);μ為流體動力粘度;SM為廣義源項;γa為a相體積分數。1.3t+u的選取{?(ρk)?t+?(ρUk)=?[(μ+μtσk)?k]+pk-ρε；?(ρε)?t+?(ρUε)=?[(μ+μtσε)?ε]+εk(Cε1pk-Cε2ρε).(3)式中:pk=μt?U(?U+?UΤ)-23?U(3μ1?U+ρk);μt=Cμρk2ε;Cε1=1.44;Cε2=1.92;Cμ=0.09;σk=1.0;σε=1.3.1.4密度和流場密度自然空化狀態是相變過程,為描述均質多相流中此相變過程定義傳輸方程:??t(γaρa)+?(γaρaUa)=Sa+Γa,(4)式中:Sa為廣義源項;Γa為由于相間質量傳遞產生的進入a相的單位體積質量源。其中Γa=Νp∑b=1Γab(Γab:由b相到a相單位質量流率,本文中a相為水,b相為水蒸汽),Γab=˙mabAab,(5)式中:˙mab是b相到a相單位界面密度的單位質量流率;Aab為相間界面密度。根據空泡動力學知,Rayleigh-Plesset方程描述空泡與水之間的相變過程。方程如下:RBd2RBdt2+32dRBdt2+2τρLRB=pc-pρL,(6)式中:RB為氣泡半徑;νL為運動粘度;ρL為液體密度;τ為表面張力;p為流場參考壓力;pc為水的蒸汽壓。在忽略表面張力以及高次項后化簡公式dRBdt=±√23pc-pρL,令蒸汽體積分數γb=VBΝB=43πRB3ΝB,其中NB、VB分別是單位體積氣泡數和氣泡體積。單位體積總質量傳遞率m˙lg=ΝBdmBdt=F3γbρgRB23|pc-p|ρLsgn(pc-p),(7)式中:F為通過實驗獲得的經驗系數,Fvap=50,Fcond=0.01.令代γnuc(1-γb)替γb,則m˙lg=F3γnuc(1-γb)ρgRB23|pc-p|ρLsgn(pc-p),(8)式中:γnuc為成核位置的體積率。通常在數值模擬中計算空泡模型時采用人為給定單位體積氣泡數NB=9×106(單位:個),而對于RB<1mm的氣泡無法得到好的仿真結果。采用γnuc(1-γb)代替γb方法,可以計算小尺寸氣泡半徑如RB=1μm,即使氣泡半徑再縮小2個數量級均可以仿真。這是因為RB=1μm,其體積分數γb<10-12,根據方程(8)可知,生成的空泡由成核位置的體積率γnuc決定,這樣避免由于計算時產生舍入誤差影響仿真結果。1.5流域外邊界的確定基于有限元的有限體積法,N-S方程中的對流項采用CFX高精度格式,它利用調整混合因子提高計算精度且平衡收斂精度。以一階迎風為例,對流項?i=?u+ξ??Δζ˙?(9)式中:?i為變量?第i個節點值;?u為變量?迎風節點值;ζ為混合因子;??為迎風節點的節點梯度;ζ˙為從迎風節點到i節點的矢量。針對不同條件下,ζ會作出相應的調整滿足計算要求獲得精確的解,例如在i節點臨近迎風節點的節點梯度變化劇烈,則混合因子趨近于0,反之趨近于1.邊界條件為:上游邊界給定來流速度和各流體體積率;下游邊界給定壓強;流域外邊界采用速度為零物面邊界,或者根據文獻采用上游邊界,但是根據模擬研究發現流域外邊界采用上游邊界時,仿真的收斂效果不如采用物面邊界,即RMS殘差曲線振蕩較劇烈;另一方面,當v>12m/s,阻塞比(即,模型截面面積與流域截面面積之比)小于0.01時,流域外邊界采用速度為零的物面邊界和采用上游邊界模擬結果無差別;當阻塞比大于0.01時,流域截面直徑影響空泡的形態。因此,本研究工作均采用小阻塞比的物面邊界作為流域的外邊界;初始條件:設定初始來流速度和各流體體積分數以及流域初始壓力。三維物體數值仿真具有計算量龐大,精度高等特點,這就使得進行高精度數值仿真時普通計算機無法完成任務。因此本仿真實驗使用了哈爾濱工業大學能源學院的深騰1800高性能服務器。2下覆蓋物的數值模擬利用商業CFD軟件CFX,針對三維圓盤空化器模型進行數值模擬,所得結果將為水下超空泡水下高速航行體減阻機理研究提供重要的理論基礎和依據。2.1空化器網格劃分運動體幾何形狀如圖1所示,(a)圖表示無攻角圓盤空化器,圓盤空化器尺寸為R=6cm,L=1cm;(b)圖表示帶α攻角圓盤空化器正視圖。尺寸與無攻角圓盤空化器一致。在仿真時二者表面均認為是速度為零的光滑物面。模型所在坐標系為三維直角坐標系,XY平面為正視圖,XZ平面為俯視圖,YZ平面為側視圖。對于本文問題,若采用非結構化網格,生成的空泡外邊界會出現凸凹不平的毛刺,影響對于仿真結果的判斷。若想消除毛刺現象,就要對網格進行細分,但會帶來網格質量差和計算量呈指數倍增長的缺點。因此采用基于六面體結構化網格劃分方法。無攻角圓盤空化器網格分布圖和流場域網格分布圖如圖2所示。總共六面體網格數為2165124,節點數為2204421,其中圓盤空化器表面網格數為5378.流域總長為14m,空化器前面流域長度為1m,為滿足阻塞比要求,流域截面積為圓形,其半徑為1.5m.其中流域設定14m是因為在小空化數下空泡的長度可達十幾米長。如表1,在速度為100m/s時,空泡長度達到11.04m.此外,根據軟件ICEMCFD評定網格質量方法可知,結構化網格的網格質量大于等于0.2足以滿足仿真要求。此空化器模型的網格質量為0.23,完全滿足仿真要求。帶攻角圓盤空化器所在流域的網格數量與其相差無幾,只是空化器帶有攻角,網格質量為0.2.因此,根據本文研究內容,采用基于六面體結構化網格劃分方法,滿足數值模擬的必要條件。3初始條件及研究目標的確定空化數是衡量流體空化程度的一個物理量,計算公式σ=(p∞-pc)/0.5ρv2,其中p∞,pc,ρ,v分別為某基準位置上的絕對靜壓力、空泡內絕對壓力、液體密度、相對于物體的未擾動液體的流速。研究中設定初始條件p∞、25℃時的水的pc和ρ.雷諾數反應流體粘性的物理量,其公式Re=v∞Dn/ν,在本文研究中其變化范圍為3.56×106～1.19×107.由于空化數較大情況下無法生成超空泡,反之,則由于情況比較復雜超出軟件的計算能力,因此設定空化數0.019～0.23為研究目標。由于計算軟件定義的上游邊界條件和初始條件是以速度或者壓力為參數,所以在仿真的時候需要將空化數轉換成速度參數,即取速度范圍在30～100m/s范圍內為研究目標,并以每次仿真速度遞增10m/s進行采樣。3.1仿真結果與實驗結果對比直徑Dn=120mm無攻角圓盤空化器,速度范圍為30～100m/s,仿真結果取空泡形態參數,即空泡的長度Lc和中截面最大直徑Dc,見表1,空泡長度與厚度隨著速度地增加而增加,即隨空化數的減小而變大。表1是以初始速度為30m/s并以10m/s為步長增加至速度為100m/s的仿真結果.烏克蘭專家Savchenko依據實驗結果和空泡截面擴大獨立性原理,利用自由流線擴展的漸近律,得到并用實驗驗證了:在自然空化情形下,空泡的基本幾何形狀是一個近似的橢球體。基本的超空泡尺寸有中截面直徑Dc和長度Lc,在2%～4%的精度范圍內可以由以下漸近關系式決定:DcDn=Cx0(1+σ)σ,LcDn=1σCx0(1+σ)ln1σ,(10)式中:Dc為中截面直徑即超空泡最大直徑;Lc為超空泡最大長度;Dn為空化器直徑;σ為空化數,Cx0空化器阻力系數。由于本文采用的是圓盤空化器,從文獻中可以得知Cx0=0.82.為進一步驗證仿真結果的可靠性,將仿真結果與經驗公式(1)進行對比,如圖3、4所示,由公式所得超空泡長度和中截面最大直徑變化曲線可知,其長度和中截面最大直徑是隨著空化數的增加而減少,仿真結果的數據點大部分落在該曲線上。即使在較大空化數下,略有不同,但是均在2%～4%精度范圍內,因此二者具有很好的一致性。其中D′c=Dc/ζ和L′c/ζ分別為經過無量綱化的空泡中截面最大直徑與長度,ζ為特征長度。此外,對比結果也說明了使用商業軟件CFX可以進行三維自然超空化的數值模擬研究。它的計算結果可以為以后的超空泡實驗研究提供參考依據。以速度為70m/s的仿真結果為例,如圖5所示,空泡最大中截面直徑Dc=0.55m和長度Lc=4.94m,圖5為三維視圖中的正視圖(XY平面)。選取方法是在25℃水蒸氣的等高圖中選取最小含水蒸汽率為50%處定為輪廓邊界,利用軟件自身的測量工具測定Dc和Lc.其中設定最小含水蒸汽率為50%是以在上述邊界條件下,速度范圍在30～100m/s中進行仿真結果與(1)式相比時最匹配為標準而設定的。另一個方面,根據(1)式計算Lc=4.94458m和Dc=0.55286m,兩者比較之后,誤差在2%～4%精度范圍之內,說明二者結果相吻合,其中自然超空泡所取正視圖的水蒸汽體積率為50%～100%.圖5中,在超空泡尾部呈現雙渦管閉合方式,這與文獻中關于超空泡尾部閉合方式相一致。如圖6所示,是其三維視圖中的三維形狀的超空泡,取含水蒸氣50%的等值曲面圖形。速度70m/s時3-D超空泡形狀是個橢球體,從其側面圖圖3看到的雙渦管形狀是超空泡尾部形成高逆壓梯度區域,迫使部分水倒流,形成回注射流。對于回注射流現象在后面會進一步證明。3.2帶攻角圓空化器的模擬采用帶攻角圓盤空化器可以為航行體提供必要升力,對于水下航行體航行起控制作用,因此具有很重要的研究意義。3.2.1無攻角部分空化器回注射流仿真條件帶攻角圓盤空化器生成超空泡形態可由圖7(a)～(b)所示,以60°攻角為例,Dn=12cm圓盤空化器,σ=0.219,航行速度v=30m/s,水蒸汽含汽率為50%～100%.從俯視圖和側視圖對比可知兩種角度觀測到的空泡厚度是不一致的,而長度一致;空泡尾部均有凹入部分,其含汽率為0,表明此處的流體為水,因此可以證明產生了十分明顯的回注射流,且形狀不同。無攻角圓盤空化器生成超空泡形態可由圖7(c)～(d)所示。仿真條件除無攻角外均與上面一致。與帶攻角圓盤空化器生成超空泡形態比較,從側視圖和俯視圖觀察其空泡厚度相等,長度也一致,尾部產生的回注射流形狀相同。如圖8所示,帶攻角圓盤空化器生成的超空泡正面投影為橢圓形,中間區域是回注射流所占位置,產生的位置不在橢圓形中心,且上下不對稱,左右對稱。而無攻角圓盤空化器生成的超空泡正面投影較均勻,呈正圓形,回注射流的位置在圓心處。回注射流現象的產生是根據汽水百分比含量推斷出來的。為了進一步驗證該區域是回注射流我們引入矢量圖9.從圖9中的水流速度矢量方向可以知道,水流繞過空化器到達空泡尾部低壓區水流速度方向發生變化,水向相反方向流動,因此進一步證明回注射流現象的發生。如圖10所示,3D帶攻角圓盤空化器生成帶回注射流的空泡,圖11所示無攻角圓盤空化器生成帶回注射流的空泡,此時二者的空化數均為0.219,所取含水蒸汽率為50%的等值面。由圖7～圖11可知,帶攻角圓盤空化器生成的超空泡長度比無攻角圓盤空化器生成的超空泡長度短。產生回注射流現象與空化器是否帶攻角無關。4數值模擬結果對比運用N-S方程及k-ε湍流模型對均質流場求解,并且利用Rayleigh-P