• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    高凝相濃度噴管兩相流研究進展①

    2017-01-05 09:40:36劉平安
    固體火箭技術 2016年6期
    關鍵詞:粒徑方程發(fā)動機

    劉平安,王 良,王 璐,王 革,郜 冶

    (哈爾濱工程大學 航天與建筑工程學院,哈爾濱 150001)

    高凝相濃度噴管兩相流研究進展①

    劉平安,王 良,王 璐,王 革,郜 冶

    (哈爾濱工程大學 航天與建筑工程學院,哈爾濱 150001)

    為了選擇合適的方法研究高金屬含量發(fā)動機噴管內的高凝相濃度兩相流,從理論研究和實驗研究兩個方面,總結了噴管兩相流研究的進展,根據(jù)各種方法的研究結果,討論了各方法對高凝相濃度噴管兩相流研究的參考價值。發(fā)現(xiàn)實際流動介于兩相平衡流和兩相凍結流之間,在高凝相濃度噴管兩相流中,顆粒含量高,兩相間相互作用強,流動趨于等溫,兩相等溫流模型適用于簡化分析。兩相流數(shù)值模型包括顆粒軌道模型和顆粒擬流體模型,兩者均可用于高凝相濃度噴管兩相流計算,軌道模型忽略顆粒體積,但能描述顆粒的運動變化過程,擬流體模型能得到顆粒參數(shù)在空間的分布。兩相流實驗研究主要用于驗證數(shù)值計算的結果,并從實驗獲得顆粒粒徑、顆粒阻力系數(shù)、Nu數(shù)等參數(shù),以用于理論計算。采用理論和實驗相結合的方法,分析顆?;騼上嗔鲌龅哪承┚植刻卣?,最終可推廣到整個噴管內,從而大大減少實驗難度。通過分析噴管兩相流的研究進展,為高金屬含量發(fā)動機內高凝相濃度噴管兩相流的研究提供了參考。

    噴管兩相流;高凝相濃度;固體火箭發(fā)動機

    0 引言

    鋁、鎂等金屬很早就作為燃燒劑用于固體推進劑中,提高了發(fā)動機性能[1-2]。發(fā)動機性能提升程度受金屬使用量和粒徑的影響,隨著推進劑的發(fā)展,金屬含量逐漸增加,粒徑從微米級向納米級過渡。在納米鋁冰發(fā)動機和金屬燃料沖壓發(fā)動機的推進劑中,金屬含量達到40%~50%,其燃燒產物中凝相氧化物含量也高達94%,這類發(fā)動機被稱為高金屬含量發(fā)動機[3-5]。金屬燃料的使用,使發(fā)動機工作效率下降,且其含量越高,效率下降越多[6-7]。發(fā)動機工作效率是由熱力計算的理論性能與實測性能比較得到的,工作效率越高的發(fā)動機,理論性能預估越準確。如果用比沖效率來衡量,當發(fā)動機中金屬含量低于30%時,其比沖效率大于90%,性能預估較準確[8],而金屬含量為50%的鋁冰發(fā)動機,實測比沖效率僅為56%~60%,性能預估不準[9]。

    鋁冰發(fā)動機工作效率低,性能預估不準,導致預先研究時實驗屢屢失敗。因此,必須研究這類發(fā)動機性能損失的成因。在所有含金屬發(fā)動機中,由金屬燃料帶來的損失包括燃燒不完全損失和兩相流損失。一般而言,金屬的燃燒效率很高,燃燒不完全損失很小,如在復合推進劑發(fā)動機中,金屬燃燒效率在96%以上[10];在鋁冰發(fā)動機中,金屬燃燒效率也大于85%[11]。另一方面,兩相流損失是噴管膨脹過程的主要損失,對發(fā)動機性能影響很大,且兩相流凝相濃度越大,兩相流損失越大[8]。由此可見,兩相流研究是含金屬發(fā)動機性能損失研究的重點。在高金屬含量發(fā)動機中,推進劑燃燒產物中顆粒很多,兩相流凝相濃度很大。如鋁冰推進劑發(fā)動機的噴管中顆粒含量可達94%[12],在金屬鎂水沖壓發(fā)動機中,水燃比為0.5時,在噴管中顆粒含量達50%[13]。一般認為,凝相顆粒含量大于50%的噴管流動為高凝相濃度兩相流,這在高金屬含量發(fā)動機的噴管中很常見。高凝相濃度兩相流帶來較大的兩相流損失,這是高金屬含量發(fā)動機性能損失的主要原因。

    為了合理評估高凝相濃度噴管兩相流損失及其對發(fā)動機性能的影響,并最終找到準確預估發(fā)動機工作性能的方法,本文從理論和實驗兩方面,總結了發(fā)動機噴管兩相流研究的進展。鑒于專門針對高凝相濃度噴管兩相流的研究很少,在總結各個研究方法之后,還討論了各種方法對高凝相濃度噴管兩相流研究的參考價值。

    1 噴管兩相流理論研究的進展

    為了計算發(fā)動機的理論性能,必先得到噴管兩相流場參數(shù),而求解控制方程得到流場參數(shù)的過程,即為噴管兩相流理論研究。噴管兩相流控制方程描述了顆粒速度up、顆粒溫度Tp、氣相速度ug、氣相溫度Tg間關系及其隨時間空間的變化。

    在適當?shù)募僭O下,一維矢量形式的控制方程組如下[14]:

    (1)

    其中,U、F和S分別為待求未知量矩陣、通量矩陣和源相矩陣,形式如下:

    式中ms為相間質量源相;p為壓強;τv為顆粒的動量松弛時間;τv=mp/[3πμrp·f(Rep)];τT為顆粒的熱松弛時間,表示為τT=mpCp/[2πrpλ]。

    噴管兩相流求解的入口條件為Tg=Tp=T0,ug=up=u0,p=pc,出口條件為p=pe。其中,T0為發(fā)動機燃燒室溫度,u0為噴管頭部的平均流速,pc為燃燒室壓強,通常認為燃燒室到噴管頭部壓強變化很小,噴管頭部壓強用燃燒室壓強表示,pe為噴管出口壓強。確定求解條件之后,可用解析模型或數(shù)值模型,求解兩相流控制方程,得到噴管兩相流場參數(shù),用于發(fā)動機性能計算。1.1 噴管兩相流解析模型的研究

    在早期研究中,研究人員試圖通過簡化假設,得到兩相流控制方程的解析解。從顆粒與燃氣間速度、溫度的關系來看,顆粒相對燃氣可能處于4個狀態(tài),分別是平衡、凍結、滯后和非平衡。“平衡”表示兩相速度相等、溫度相同,即ug=up,Tg=Tp;“凍結”表示燃氣或顆粒的速度或溫度保持入口時的值(u0、T0)不變,其中顆粒速度凍結表示計算不考慮顆粒的影響;“滯后”是顆粒溫度比燃氣高,速度比燃氣低的現(xiàn)象,此處特指顆粒參數(shù)與燃氣參數(shù)呈一定已知函數(shù)關系的情況;“不平衡”表示顆粒雖然滯后于燃氣,但顆粒參數(shù)與燃氣參數(shù)的關系不能確定,顆粒參數(shù)與燃氣參數(shù)是兩組獨立的變量。當兩相間速度或溫度關系為平衡、凍結或滯后時,原本要求解的兩組參數(shù)(ug和up,Tg和Tp)變成了一組參數(shù)(u,T),求解過程將被大大簡化。采用這種方法,可得到兩相流控制方程的解析解。

    為了便于分析,定義顆粒的速度滯后數(shù)為(1-K),其中K為兩相速度系數(shù),K=up/ug;顆粒的溫度滯后數(shù)為(1-L),其中L為兩相溫度系數(shù),L=(T0-Tp)/(T0-Tg)。用K和L替換變量up和Tp,在適當?shù)募僭O下,把顆粒參數(shù)變換為氣相參數(shù)的等效形式,由此可通過理論推導得到噴管兩相流的解析解。已有研究成果可按不同K和L值列為表1中的各個解析模型。表1中,Dk為擴散系數(shù);C1、C2表示K或L為常數(shù);“×”表示系數(shù)未定;Tg或Tp等于T0表明,燃氣或顆粒溫度保持不變。

    表1 噴管兩相流的解析模型Table 1 Analytic models for two-phase nozzle flow

    上述解析模型主要可分為4類:無滯后/凍結模型、單種滯后模型、等溫流模型、全滯后模型。下面將分別進行闡述。

    (1)無滯后/凍結模型

    不考慮/過分考慮顆粒參數(shù)與燃氣參數(shù)不相等的情況,主要包括兩相凍結流模型和兩相平衡流模型[15]。從表1中可看出,在這些模型下,顆粒和氣相的參數(shù)變成了已知參數(shù)或一組未知參數(shù)(u、T),因此很容易求解析解。但在實際情況下,顆粒由于具有高密度和熱惰性,在隨燃氣通過噴管的過程中,其速度和溫度不可能與燃氣相同。因此,這類模型存在固有的局限性。“兩相平衡流模型”常用于噴管熱力計算[16],可得到發(fā)動機的最大理論性能,如發(fā)動機比沖。“兩相凍結流模型”假設燃氣流動不受顆粒影響,由于忽略了其貢獻的推力,計算得到的發(fā)動機理論性能很低。這兩種模型是噴管兩相流的極限情況,計算得到的發(fā)動機理論性能參數(shù)值分別為最大和最小,實際流動位于兩者之間。在高凝相濃度噴管兩相流中,顆粒對燃氣流動影響很大,且顆粒向外排出產生推力,是發(fā)動機推力的重要來源。因此,兩相凍結流模型并不適用。而當顆粒增多時,兩相間耦合作用增強,兩相接近平衡態(tài),平衡流模型可用于分析高凝相濃度噴管的最大理論性能。

    (2)單種滯后模型

    單種滯后模型只計算速度和溫度滯后中的一種,要么假定另一項平衡。相比無滯后/凍結模型,考慮了顆粒的滯后效應,主要包括單顆粒動力學模型、小滑移模型[17]、速度平衡溫度不平衡模型[18]、速度常滯后溫度平衡模型[19]。

    在凍結流模型的基礎上,“單顆粒動力學模型”用已經(jīng)求得的氣相流場計算顆粒沿軌道的參數(shù)?!靶』颇P汀奔僭O兩相間速度滯后等于顆粒作為混合物組元的擴散。這兩個模型均不考慮顆粒和燃氣的溫度,在噴管兩相流中使用很少。為此,研究人員假定兩相速度或溫度平衡,僅討論一種滯后對結果的影響。Soo[15]發(fā)現(xiàn),速度滯后對發(fā)動機性能影響更大,而線性加速假設能簡化非平衡情況的求解。Dillon[18]則在速度平衡假設下,對兩相凝相濃度為8的噴管進行計算,提出了“速度平衡溫度不平衡”模型,說明在高凝相濃度下,顆粒與燃氣溫度不平衡較小。于是,研究重點轉向了速度不平衡。Gilbert等[19]提出了“速度常滯后溫度平衡”模型,該模型采用了氣相線性加速的假設。計算結果發(fā)現(xiàn),顆粒粒徑越大,速度滯后越大,發(fā)動機性能下降越多。這說明速度滯后對燃氣在噴管中的膨脹過程影響很大。由于僅考慮了一種滯后,單種滯后模型目前已很少用,對高凝相濃度噴管兩相流計算的參考價值不大。

    (3)等溫流模型

    等溫流模型是單種滯后模型的特殊情況。當噴管兩相流中凝相的濃度很大時,由于顆粒多,兩相間熱交換量大,燃氣和顆粒的溫度可能保持平衡,且在流動過程中保持不變,即燃氣和顆粒的溫度均被“凍結”,這種現(xiàn)象稱為噴管的“等溫膨脹”。在溫度凍結假設下,研究模型包括速度常滯后且燃氣顆粒溫度均凍結模型(速度常滯后等溫流模型)、速度平衡且燃氣顆粒溫度均凍結模型(速度平衡等溫流模型)。

    其中,“速度常滯后等溫流模型”最常用于噴管中,由Rudinger[20-21]提出,遲偉等[12,14,22]用其解決了實際噴管流動問題。該模型求解方法簡單,凝相濃度較大(η=0.1~17)時,能得到準確的計算結果,適用于噴管兩相流的簡化計算。當凝相濃度更大時,需要考慮顆粒體積的影響。另外,研究發(fā)現(xiàn)凝相濃度越大,兩相速度越接近平衡,因此提出“速度平衡等溫流模型”。在高凝相濃度噴管兩相流計算中,應以等溫流為前提,顆粒速度滯后與速度平衡的等溫流模型均可采用。對于鋁冰發(fā)動機的噴管,由于兩相流凝相濃度高達17,應采用速度平衡的等溫流模型。

    (4)全滯后模型

    全滯后模型完整地考慮了顆粒速度和溫度滯后,包括兩相常滯后模型、兩相任意滯后模型和兩相全耦合模型。“兩相常滯后模型”假定兩相間溫度和速度滯后數(shù)為常數(shù),即K=C1,L=C2[23-27];“兩相任意滯后模型”假定滯后數(shù)為噴管軸向坐標x的特定函數(shù),即K=K(x),L=L(x),采用逆解法求解[28-29]。這2個模型均把兩相控制方程簡化為等效單一流體方程求解,但都需指定K和L,模型被過分簡化,計算結果與實際有較大偏差。因此,在實際計算中應用很少[30]?!皟上嗳詈稀蹦P筒恍柚付↘或L,而采用變量替換和迭代法直接求解控制方程,得到燃氣和顆粒各自參數(shù)[31-33]。因此,其與實際更接近,應用更多。計算發(fā)現(xiàn),顆粒尺寸對兩相滯后影響很大,在噴管喉部處,顆粒滯后最大。對高凝相濃度噴管兩相流,常滯后和任意滯后模型的結果與實際偏差較大,且很難通過模型改進來提高計算精度。因此,這兩個模型參考價值不大。但全耦合模型假設少、計算精度高,可用于高凝相濃度噴管兩相流計算。

    噴管兩相流本身很復雜,兩相流控制方程中包含的兩相相互作用、傳熱等,使方程求解很困難。因此,解析模型往往不能很好描述實際流動情況,在高精度求解時,需用數(shù)值模型代替。

    1.2 噴管兩相流數(shù)值模型的研究

    數(shù)值模型采用數(shù)值方法直接求解兩相流控制方程。因此,能完整考慮顆粒的速度和溫度滯后,計算結果更準確,但計算量更大。數(shù)值模型的使用主要分為三步:(1)在噴管物理模型上,建立適合求解的兩相流控制方程、(2)選用合適的數(shù)值計算方法、(3)設置合理的求解條件,以進行計算并輸出結果。其中,控制方程的建立和計算方法的選定是噴管兩相流計算的核心和難點。根據(jù)描述顆粒參數(shù)的控制方程和求解方法不同,數(shù)值模型主要可分為顆粒軌道模型(歐拉-拉格朗日模型)和顆粒擬流體模型(歐拉-歐拉模型)[34]。

    1.2.1 顆粒軌道模型

    顆粒軌道模型用拉格朗日表達方法描述噴管中的顆粒,并沿顆粒軌道求解其參數(shù),用歐拉方法求解噴管中燃氣的參數(shù)。

    (1)模型方程的處理方法

    早期的顆粒軌道模型未考慮顆粒質量、動量及能量的湍流輸運,并把模型方程寫為歐拉坐標系下的表達式[17]。穩(wěn)態(tài)模型控制方程組在噴管的跨音速段是奇異的,直接求解很困難。為此,研究者們按照方程維度不同進行了許多探索[21,35]。

    在一維模型方程下,消除方程奇異性的方法主要有3種:變量置換法、指定壓力法和時間相關法[36]。“變量置換法”通過變量替換消除方程的奇異性,“指定壓力法”通過選取合適的壓強,以減少未知量的個數(shù),從而消除方程的奇異性,并根據(jù)流道面積,判斷是否收斂[34-35]。一維控制方程求解時,需把噴管分為亞音速段、跨音速段和超音速段,上述兩種方法僅限于跨音速段的求解?!皶r間相關法”在穩(wěn)態(tài)控制方程中引入非穩(wěn)態(tài)項,以消除方程的奇異性,并用準穩(wěn)態(tài)的思想求得穩(wěn)態(tài)噴管兩相流場[37-38]。該方法不需對噴管分段,而是對整個噴管求解,相比變量置換法,方程更簡單;由于不需反復選取參數(shù),相比指定壓力法,計算更容易。因此,在在一維噴管兩相流求解中被廣泛應用,也被推廣到多維問題的求解上。

    二維軸對稱和三維模型更符合物理實際,應用更多,消除方程奇異性的方法主要有常滯后假設法、數(shù)值迭代松弛法、相容條件法和時間相關法?!俺蠹僭O法”使用常滯后假設,簡化求解噴管跨音速段的流場。因此,避免了方程的奇異問題。經(jīng)實驗證明,在該區(qū)域使用常滯后模型是合適的[21,39-42]。但由于計算精度低,該方法僅適用于小曲率噴管?!皵?shù)值迭代松弛法”首先用兩相平衡流模型,計算得到噴管初始流場;然后,假定燃氣參數(shù)不變,計算顆粒參數(shù),再用新的顆粒參數(shù)計算新的氣相流場,反復迭代直至收斂[43]。該方法也僅用于噴管跨音速段的求解,沒有考慮方程的混合型性質。因此,也沒有真正讓兩相全耦合,計算精度很低。“相容條件法”通過引入相容性條件,消除了方程的奇異性,在由流線及其法線所組成的坐標系中,求解控制方程組,假設燃氣的流線由初場確定,且不受顆粒的影響[44]。因此,改進了迭代松弛方法,提高了計算精度。上述數(shù)值方法由于在每次迭代時,僅計算燃氣或顆粒的參數(shù),都是兩相非耦合方法?!皶r間相關法”能消除方程的奇異性,并實現(xiàn)兩相全耦合,最早由Chang[37,45]應用在二維軸對稱噴管求解上。由于兩相耦合程度高,相比前幾種方法,其計算精度大大提高,因此應用更廣。該方法可對整個噴管和發(fā)動機內流場直接求解,且能用于邊界復雜的噴管。

    由上述討論可見,“時間相關法”適用于噴管穩(wěn)態(tài)兩相流場計算。為了提高計算準確度,還需考慮顆粒湍流對流場的影響。因此,在模型方程中,應包含湍流項。考慮湍流會增加計算的難度,但求解方法不需改變太大。根據(jù)方程形式及湍流處理方法不同,軌道模型可分為(有湍流影響的)確定軌道模型和隨機軌道模型[17]。

    “確定軌道模型”求解的是時均拉格朗日形式的顆??刂品匠?,并從擬流體假設中,引入顆粒的漂移速度和漂移力的概念,以考慮湍流擴散對顆粒軌道的影響?!半S機軌道模型”直接求解軌道拉格朗日形式的顆粒瞬時控制方程,沿顆粒軌道計算顆粒的湍流參數(shù)。

    (2)數(shù)值計算方法

    控制方程的形式確定后,需使用合適的方法求解。在一維模型方程下,時間相關法最終得到一階非線性常微分方程組,可采用Rugge-Kutta法求解。但二維軸對稱和三維模型方程很復雜,只能采用數(shù)值積分法求解。氣相的數(shù)值積分方法與單流體相同,而顆粒的數(shù)值方程是拉格朗日形式的控制方程,除顆粒參數(shù)外,還需求解顆粒軌道。因此,一般用Rugge-Kutta積分法和線性近似法計算,而在隨機軌道模型下,需使用Monte-Carlo法[46-54]。另外,邊界條件和數(shù)值積分格式對結果準確性的影響也很大。其中,數(shù)值格式直接影響計算的精度,研究人員進行了大量的研究[55-56]。

    顆粒軌道模型計算量小,能模擬有復雜經(jīng)歷的顆粒,沒有偽擴散,但不易考慮顆粒的湍流擴散過程,在復雜流場下,難以給出連續(xù)的顆粒參數(shù)分布。因此,很難和實驗對照。在噴管中,顆粒會經(jīng)歷相變、變形、成長、團聚、破碎、碰撞、沉積、輻射和化學反應等各種情況,用顆粒軌道模型計算時,能清楚反映顆粒變化過程[13,56,57-63]。對高凝相濃度噴管兩相流,噴管中顆粒很多,顆粒間會相互干擾,顆粒運動受湍流影響很大,應采用隨機軌道模型進行模擬。

    1.2.2 顆粒擬流體模型

    擬流體模型假設顆粒為擬流體(具有粘性、導熱和擴散等流體屬性,與氣相相互滲透,且占有一定體積分數(shù)[64]),模型方程是歐拉形式的顆粒方程和氣相方程,能完整考慮滑移、相間耦合和湍流擴散,最終得到連續(xù)的顆粒參數(shù)分布。

    (1)模型方程的處理方法

    顆粒擬流體模型求解湍流兩相流時均守恒方程組,相比顆粒軌道模型,更完整更嚴格地模擬氣相湍流和顆粒湍流,仍采用時間相關法分析和計算。顆粒湍流在擬流體模型研究中是一個重要概念,包括顆粒與流體的湍流相互作用、顆粒自身的湍流輸運(湍流粘性、擴散和導熱等)、顆粒間碰撞引起的湍流輸運、顆粒質量擴散引起的湍流輸運等[17]。擬流體模型中湍流的描述方法來源于單流體湍流模型,后者常用的有兩種:渦粘性模型中的雙方程模型(如k-ε模型)、雷諾應力模型中的代數(shù)應力模型。在擬流體模型研究過程中,以時間相關法為基礎,逐步發(fā)展了層流氣相-層流顆粒的、湍流氣相-層流顆粒的以及湍流氣相-湍流顆粒的擬流體模型,這些模型曾被用于噴管兩相流計算中,但隨著計算技術的進步,已逐漸被淘汰[65-73]。最新研究是完整考慮湍流影響的顆粒擬流體模型,最常見的是k-ε-kp-εp四方程湍流模型[73],其他模型也在研究中。這些模型方程應采用合適的數(shù)值計算方法求解。

    (2)數(shù)值計算方法

    顆粒擬流體模型方程采用數(shù)值積分法求解,可沿用并改進單相流的求解方法,如SIMPLE算法。方程離散方法主要有有限差分法(FDM)、有限體積法(FVM)和有限元法(FEM)。其中,有限體積法以有限差分法為基礎,目前大部分數(shù)值積分方法都是在有限體積法的基礎上發(fā)展出來的,這些方法[55]包括用黎曼精確解或近似解求解界面通量的Godunov方法、通過方程中通量的變形而衍生的通量求解方法(如通量分裂求解的Van leer方法和Steger-Warming方法),還有基于流動基本規(guī)律的一些積分方法(如總變差不增TVD格式)。在多維模型方程中,界面通量求解很困難,常用矢通量分裂法得到各方向上的一維黎曼算子,用Godunov法求解該方向上的黎曼近似解,最終組合為整個流場的解[56]。合適的數(shù)值積分格式能保證計算的完成,目前常用的數(shù)值格式包括一階/二階迎風格式、混合格式、乘方格式、HLL、AUSM、MUSCL、Lax-Wendroff、NND格式等,近年來有許多數(shù)值格式的研究,新的高精度格式不斷出現(xiàn)。數(shù)值格式在選取時,需保證穩(wěn)定性、守恒性等法則,應根據(jù)實際情況選取。擬流體模型方程的求解需使用TDMA算法、低松弛迭代、源相線性化等手段,與單相流類似。

    用時間相關法及合適的數(shù)值格式求解控制方程時,需要滿足CFL條件:Δt≤Δx/(ug+c),計算才能收斂。為提高計算速度,侯曉等[74]用“隱式近似因子分解法”,在差分格式上做了一些改變,減弱了該條件的影響,當噴管中顆粒較小(Dp<1 μm)時,該方法較適用。另外,噴管流場的邊界條件對數(shù)值計算結果影響也很大,為了確定合適的噴管入口邊界條件,方丁酉提出用虛擬噴管法確定兩相混合物的入口速度[36]。

    顆粒擬流體模型完整考慮了顆粒的湍流輸運,能得到顆粒在流場中連續(xù)的速度、溫度、體積分數(shù)和湍流參數(shù)的分布,還能考慮顆粒碰撞引起的顆粒參數(shù)的變化,結果易于與實驗對照。但該模型不能描述單個顆粒的運動和變化情況,存在偽擴散,當顆粒分組較多時計算量大。對于高凝相濃度噴管兩相流,由于顆粒含量高,碰撞可能性大,受湍流影響也很大,要重點考慮這些因素的影響時,用擬流體模型較合適。

    上述理論模型都是從不同角度對真實過程進行的簡化和近似,兩相平衡流模型是較大的簡化。實際中,常用的是顆粒軌道模型和顆粒擬流體模型。兩者側重點不同,顆粒軌道模型適用于描述顆粒的詳細運動和變化過程,顆粒擬流體模型適用于描述顆粒參數(shù)的空間分布。對于高凝相濃度噴管兩相流,應按照研究側重點不同,選擇合適的計算方法。

    2 噴管兩相流實驗研究的進展

    同一噴管流動過程在不同理論模型或計算方法下得到的結果可能不同,甚至互相矛盾[75]。兩相流實驗研究的主要目的是驗證理論計算結果,用以輔助工程設計。實驗研究的手段包括噴管兩相流場參數(shù)測量、兩相相互作用測量、顆粒粒徑測量和顆粒動力學特性測量等。隨著測量技術的進步,兩相流實驗研究發(fā)展很快。

    2.1 噴管兩相流場參數(shù)的測量

    兩相流場參數(shù)測量實驗包括冷流模擬實驗和發(fā)動機實驗兩種,前者由于和實際情況差別很大,難以模擬真實的流動狀態(tài),而后者對測量技術的要求很高,常規(guī)測量方法在高溫高速的氣固兩相噴管流動環(huán)境中往往失效。在發(fā)動機實驗(包括簡化或縮比的噴管流動特性實驗)中,噴管流場參數(shù)的測量方法主要是非接觸測量法[76]。該方法以攝影、射線、激光和圖像處理技術為基礎,包括X射線高速實時熒屏分析(RTR)、粒子圖像測速(PIV)、激光多普勒測速(LDV)、高速攝影和重曝攝影等。通??蓽y得顆粒的分布規(guī)律、顆粒速度和軌跡,這些測量結果與理論計算結果比較,可用于驗證模型的準確性。

    早期研究從兩相不平衡角度出發(fā),用消光測量法(Light Extinction Measurements)和鈉D-線反轉技術(Sodium D-line Reversal Technique),測量了顆粒的速度滯后和溫度滯后,證明了滯后效應的存在,并得到了與一維解析模型吻合的結果[77-78]。后來,Neilson等[79]為了驗證非耦合模型的計算結果,用重曝攝影法分析了顆粒運動速度,發(fā)現(xiàn)非耦合模型僅在顆粒凝相濃度很小時適用。20世紀70~80年代以來,噴管兩相流理論計算以顆粒軌道模型和顆粒擬流體模型為主,噴管兩相流實驗研究主要用于觀察擴張段是否存在無顆粒區(qū)、顆粒軌跡是否會越過噴管軸線、顆粒運動受重力影響的程度等問題[75-80]。在所有測試方法中,RTR技術不需使用可視噴管,應用最廣。對于高凝相濃度噴管兩相流,目前理論研究還很不完善,須用合適的實驗方法,驗證其計算結果。

    2.2 兩相相互作用實驗研究

    噴管中顆粒和燃氣間相互作用主要是作用力和傳熱,其中作用力包括顆粒的氣動阻力、升力、體積力等,一般情況以阻力為主;兩相傳熱則包括熱傳導、對流換熱和輻射,以對流和輻射為主。為便于計算,在控制方程中用阻力系數(shù)CD計算作用力,用Nu數(shù)計算對流傳熱。在均勻無粘低速兩相流中,可僅考慮氣動阻力和導熱,通過理論推導得到CD=24/Re,Nu=2,這2個值又被稱為斯托克斯阻力系數(shù)和傳熱系數(shù),用CDs和Nus表示[34]。實際情況比這一情況復雜得多,只能通過實驗測量得到CD和Nu數(shù)的經(jīng)驗計算式。

    (1)兩相阻力系數(shù)CD的實驗測量

    標準阻力曲線是測量不可壓無限大流場中剛性球體的受力得到的,可用式(2)擬合:

    (2)

    在實際情況下,阻力系數(shù)受許多因素的影響。Bailey等[81]在Ma=0.1~20、Re=0.2~106的參數(shù)范圍內,測量了顆粒阻力,并得到如下計算式:

    (3)

    式中 阻力系數(shù)的各個修正項依次為雷諾數(shù)函數(shù)、湍流效應、氣相稀薄效應、可壓縮性、顆粒非球形修正和溫差效應修正。

    式(3)考慮得比較全面,但公式形式復雜。目前,在噴管兩相流計算中,最常用的是標準阻力曲線,以及Carlson等[82]的經(jīng)驗公式:

    (4)

    另外,也常使用Hermsen[83]的經(jīng)驗公式:

    (5)

    其中,CD,inc為由標準阻力曲線擬合的阻力系數(shù),G(Re)和h(Ma)分別為Re數(shù)和Ma數(shù)的函數(shù)。

    除了上述公式外,還有許多用于噴管兩相流研究的阻力系數(shù)計算經(jīng)驗公式,遲鴻偉[14]比較了各個公式,結果發(fā)現(xiàn),Hermsen的公式是可靠的。

    對高凝相濃度噴管兩相流,顆粒間空隙小,大尺度湍流被分割,氣相稀薄效應影響大,阻力系數(shù)變大,CD計算式中應引入顆粒濃度影響[84]。另外,顆粒形狀不規(guī)則,其變化過程也會對阻力系數(shù)產生影響[62]。在實驗中,完全量化這些因素的影響很困難,Koike等[85-86]用PIV和LDV方法,實驗測量了顆粒的運動速度,并分析了顆粒受力,與理論結果比較,對阻力系數(shù)經(jīng)驗公式進行了修正。這一方法可推廣到高凝相濃度噴管兩相流中。

    (2)兩相傳熱Nu數(shù)的實驗測量

    在不可壓流中的球形顆粒,相間溫差較小時,兩相傳熱以對流換熱為主,實驗測量得到Nu數(shù)平均值的經(jīng)驗公式如下:

    Nu0=2+0.459Re0.55Pr1/3

    (6)

    其中,1

    為了考慮湍流或稀薄效應的影響,Kavanau等[87]通過實驗得到如下經(jīng)驗公式:

    (7)

    該公式在顆粒較多,克努森數(shù)Kn≤0.1時適用。

    另外,Koshmarov等[88]在整個噴管中得到了Nu數(shù)計算公式,Nelson等[89]在高溫噴管流動情況下,驗證了上述兩個經(jīng)驗公式的正確性。

    (3)兩相輻射的實驗測量

    在噴管兩相流中,輻射傳熱以顆粒輻射為主,包括顆粒對壁面、燃氣的輻射和顆粒間輻射。該問題的求解需要積分光譜輻射傳輸方程(spectral radiative transfer equation),常用的積分方法是有限體積法[90]。顆粒離散相的輻射特性參數(shù)的確定是兩相流輻射問題的難點,這些參數(shù)很難通過理論分析得到,往往需要通過實驗測量,并結合輻射反問題模型反演求解[91]。顆粒輻射特性參數(shù)一般是指顆粒的光學常數(shù)(復折射率),該參數(shù)與溫度相關,通過實驗和理論求得。實驗方法可分為測反射率法、測方向散射率法和測單色透射率光譜法,每種方法的實驗結果均可結合相應的理論模型反演計算,最終可求得顆粒復折射率。

    發(fā)動機噴管中顆粒的溫度很高,一般是顆粒對外輻射散熱。目前,沒有直接測量噴管內高溫顆粒輻射的研究,大部分是在爐膛內測得。最初的研究采用傅立葉紅外光譜儀(FTIR)和KBr吸收渣片等實驗手段,測量了煤粉或爐灰顆粒在高溫下的復折射率[92]。齊宏[91]設計并搭建了高溫粒子輻射特性測試綜合實驗平臺,通過測量自然狀態(tài)下Al2O3顆粒的光譜等效透射比,反演得到了顆粒的復折射率。Al2O3顆粒的復折射率與其溫度、粒徑、濃度等參數(shù)相關。在齊宏的實驗中,顆粒濃度為6.563×10-5,而在鋁冰發(fā)動機噴管中,顆粒凝相濃度為17時,濃度約為5.6×10-3??梢姡F(xiàn)有實驗結果與高凝相濃度噴管兩相流實際情況相差仍較大,在求解精度要求不高時,可用于近似計算。

    2.3 顆粒粒徑的實驗測定

    顆粒粒徑直接影響噴管兩相流動狀態(tài),并最終影響發(fā)動機性能,只能通過實驗測得。顆粒團聚和破碎會對粒徑產生影響,在不破壞顆粒原本狀態(tài)的情況下,實驗測得的粒徑包括單個顆粒的粒徑和顆粒團聚物的粒徑,統(tǒng)稱顆粒粒徑。

    在簡單噴管兩相流計算中,顆粒粒徑往往采用平均粒徑[93],一般用如下經(jīng)驗公式計算[8]:

    (8)

    式中p為噴管頭部壓強;ξ為顆粒含量;L*為發(fā)動機內流場特征長度;Dt為噴管喉徑。

    實際情況是粒徑呈一定規(guī)律分布,大小粒徑同時存在,僅用平均粒徑不足以描述兩相流動的實際情況。測量粒徑分布的實驗方法分為接觸測量法和非接觸測量法。接觸測量法先收集顆粒然后用粒度測定儀、顆粒計數(shù)器、電鏡觀察計數(shù)等方法得到粒徑分布[94]。顆粒收集方法包括:

    (1)把收斂噴管和裝有惰性冷卻液的容器相連,從而收集通過噴管喉部的顆粒[95]。

    (2)用冷卻容器或冷卻噴流法收集噴管出口羽流中的顆粒[94,96-97]。

    (3)發(fā)動機工作完成后,從噴管壁面殘留物中收集顆粒[8]。

    前2種方法應用較多,Sehgal[94]系統(tǒng)地研究了推進劑中金屬含量及粒徑、燃燒室溫度、壓強、噴管型面、顆粒駐留時間等對噴管出口處顆粒粒徑的影響,發(fā)現(xiàn)壓強是顆粒粒徑的主要影響因素。婁永春等[95,98-99]通過改進的顆粒收集裝置測量發(fā)現(xiàn),顆粒聚集程度越大,顆粒粒徑越大;壓強越大,顆粒粒徑越大;聚集情況下顆粒粒徑分布范圍較寬,往往呈雙峰或多峰分布,主峰所占比重很大。在高凝相濃度噴管兩相流中,顆粒含量高,在整個噴管流動范圍內都可稱為“顆粒聚集”的狀態(tài),上述顆粒粒徑分析結果可用于高凝相濃度噴管兩相流研究。

    非接觸測量法用RTR技術、高速攝影、激光衍射等方法,直接獲取噴管內或噴管出口的兩相流圖像/數(shù)據(jù),再用計數(shù)法、光譜分析等統(tǒng)計方法,分析圖像或光譜數(shù)據(jù),得到粒徑分布情況[100-102]。Koo等[100]總結了各種測量方法的優(yōu)缺點,把非接觸測量法分為無圖技術(nonimaging techniques)和有圖技術(imaging techniques),有圖技術包括攝影、全息攝影和自動圖像處理技術,無圖技術包括單粒子計數(shù)器和光衍射的集成粒徑測量技術(Ensemble Particle Sizing Technique)。后者精度更高,適合用于噴管尾流中顆粒粒徑的分析。張宏安等[101]用激光全息光學方法,測量了噴管出口顆粒粒徑分布曲線,發(fā)現(xiàn)噴管出口顆粒粒徑比原始金屬粒徑小,顆粒破碎速率大于團聚速率,顆粒粒徑隨鋁含量增高而增大。Youngborg等[102]用光衍射法和高速攝影法,測量了開窗的端燃發(fā)動機噴管入口和出口處顆粒粒徑,發(fā)現(xiàn)在噴管收斂段內,凝相粒徑是以15~19 μm為中心的單峰分布。在噴管出口,顆粒粒徑變小,呈多峰分布。非接觸測量法目前仍未用于噴管高凝相濃度兩相流中顆粒的測量,該方法不影響流場,采用RTR技術能得到噴管內任意位置顆粒粒徑,值得進一步研究。

    2.4 顆粒動力學特性的研究

    在帶潛入噴管的發(fā)動機中,顆粒會進入噴管背壁區(qū)沉積,并形成殘渣,極大程度地影響發(fā)動機的正常工作,這是兩相流研究的重點。用RTR技術可測量顆粒沉積的過程,但分辨率不高,而用冷流模擬可近似地描述顆粒的運動過程。韓新波[103]用冷流模擬和相位多普勒粒子分析儀(PDPA)測速技術,分析了潛入噴管和其背壁區(qū)內的兩相流動,驗證了雙流體數(shù)值模型。

    在有過載的發(fā)動機中,顆粒會定向運動,在發(fā)動機或噴管局部形成高凝相濃度兩相流,最終造成燃燒室或噴管材料的燒蝕。婁永春[95]用收斂噴管聚集形成高濃度兩相流,對多層復合材料進行了燒蝕實驗,實驗后發(fā)現(xiàn),高強度石墨噴管被燒蝕。陳劍等[104]研究了發(fā)動機正常工作和過載狀態(tài)下噴管的燒蝕,驗證了顆粒隨機軌道模型用于預測噴管燒蝕現(xiàn)象的準確性。在長尾噴管中,顆粒滯留時間長,也很容易造成噴管燒蝕[59]。另一方面,由于難燒蝕噴管喉襯材料的使用,兩相流更可能導致顆粒在噴管壁面和喉部的沉積,導致發(fā)動機工作性能改變。賈林祥等[105-106]通過總結噴管顆粒沉積層厚度的數(shù)據(jù),分析了沉積速率,討論了顆粒沉積對發(fā)動機性能的影響。

    在噴管中,顆粒的相變、碰撞、破碎、團聚等過程不僅會影響顆粒粒徑和運動狀態(tài),也會對氣相流動產生影響,最終影響發(fā)動機性能。劉叢林[62]研究了氧化鋁形成和非穩(wěn)態(tài)相變凝結過程,發(fā)現(xiàn)燃燒室內有20%的Al2O3會在鋁顆粒表面凝結形成大顆粒,剩余的80%凝結為小液滴(煙霧)。目前,尚未在實驗中觀察到氣態(tài)的Al2O3,而一般推進劑的燃燒溫度不高于Al2O3的沸點(3 800 K,常壓)。因此,可認為噴管內是由惰性燃氣和Al2O3顆粒、小液滴組成的兩相流。通過噴管時,小液滴由于布朗運動和慣性運動的作用,會向大顆粒表面凝結,這實際上是小顆粒與大顆粒的碰撞。夏勝勇[107]詳細研究了氧化鋁顆粒的碰撞過程,并指出氧化鋁液滴屬于高粘度、高表面張力系數(shù)的液滴,直接研究顆粒在噴管中的碰撞過程很困難。Averin等[108]用激光診斷方法,對發(fā)動機噴管內氧化鋁顆粒的碰撞團聚過程進行了測量,驗證了其顆粒碰撞模型的正確性。

    3 結束語

    本文綜述了噴管內兩相流理論和實驗研究的進展,討論了各種研究方法和研究結果對高凝相濃度噴管兩相流研究的參考價值。

    理論研究包括解析模型研究和數(shù)值模型研究,解析模型中,兩相平衡流模型能夠用于計算噴管的最大理論工作性能,等溫流模型是對高凝相濃度噴管兩相流的一種假設,可用于工程計算。數(shù)值模型包括顆粒軌道模型和擬流體模型,能完整地考慮兩相速度、溫度滯后,并計算湍流對兩相流動的影響,這2種模型都能用于計算高凝相濃度噴管兩相流,并得到較好結果。理論研究是從不同角度對真實兩相流動過程的近似和簡化,目前仍不完善,還需用實驗驗證。

    噴管兩相流的實驗研究主要用于理論模型的驗證,并可由實驗觀察得到的一些規(guī)律反演出新的理論模型。從目前研究結果看,實驗能很好地驗證數(shù)值計算結果。理論計算所需的顆粒參數(shù),如顆粒粒徑、阻力系數(shù)和傳熱等,都是從實驗測量得到。因此,實驗與理論密不可分。

    總之,高凝相濃度噴管兩相流研究還很不完善,還需要進一步研究。

    [1] 龐維強,樊學忠.金屬燃料在固體推進劑中的應用進展[J].化學推進劑與高分子材料,2009,7(2):1-5;14.

    [2] 龐愛民,黎小平.固體推進劑技術的創(chuàng)新與發(fā)展規(guī)律[J].含能材料,2015,23(1):3-6.

    [3] Risha G A,Connell Jr T L.Novel energetic materials for space propulsion[R].Nasa:DTIC Document:A546818,2011.

    [4] 劉巍.固體燃料沖壓發(fā)動機燃燒組織技術研究[D].長沙:國防科技大學,2010.

    [5] Hu J X,Han C.Experimental investigation on combustion of high-metal magnesium-based hydroreactive fuels[J].Journal of Propulsion and Power,2013,29(3):692-698.

    [6] Farrow D D.A theoretical and experimental comparison of aluminum as an energetic additive in solid rocket motors with thrust stand design[D].University of Tennessee,2011.

    [7] Geisler R L,Kinkead S A.The relationship between solid propellant formulation variables and motor performance[C〗//11th Propulsion Conference,AIAA 75-1199.

    [8] Coats D E,Nickerson G R.A computer program for the prediction of solid propellant rocket motor performance[R].CA:Edwards Air Force Base,1981.

    [9] Wood T D.Feasibility study and demonstration of an aluminum and ice solid propellant[D].Purdue University,2010:227.

    [10] 張勝敏,胡春波.固體火箭發(fā)動機燃燒室凝相顆粒燃燒特性分析[J].固體火箭技術,2010,33(3):256-259.

    [11] Risha G A,Sabourin J L.Combustion and conversion efficiency of nanoaluminum-water mixtures[J].Combustion Science and Technology,2008,180(12):2127-2142.

    [12] 楊寒.鋁冰固體火箭發(fā)動機氣固兩相流噴管參數(shù)研究[D].哈爾濱:哈爾濱工程大學,2013.

    [13] 胡凡.鎂基燃料水沖壓發(fā)動機理論分析與試驗研究[D].長沙:國防科技大學,2008.

    [14] 遲鴻偉.大負載比氣固兩相流噴管參數(shù)研究及通道形狀優(yōu)化[D].哈爾濱:哈爾濱工程大學,2011.

    [15] Soo S L.Gas dynamic processes involving suspended solids[J].AIChE Journal,1961,7(3):384-391.

    [16] Gordon S,Mcbride B J.Computer program for calculation of complex chemical equilibrium compositions and applications[R].Nasa,Cleveland,OH,United States:NASA Lewis Research Center,1994:58.

    [17] 周力行.湍流兩相流動與燃燒的數(shù)值模擬[M]. 北京:清華大學出版社,1991:81-164.

    [18] Dillon P L P.Heat transfer between solid particles and gas in a rocket nozzle[J].Journal of Jet Propulsion,1956,26(12):1091-1097.

    [19] Gilbert M,Davis L.Velocity lag of particles in linearly accelerated combustion gases[J].Journal of Jet Propulsion,1955,25(1):26-30.

    [20] Rudinger G.Gas-particle flow in convergent nozzles at high loading ratios[J].AIAA Journal,1970,8(7):1288-1294.

    [21] Rudinger G.Fundamentals of gas particle flow[M].Elsevier,2012.

    [22] 王天祥,何利民.氣液兩相噴嘴等溫流動模型[J].機械工程學報,2008,44(1):121-125.

    [23] Kliegel J R.One dimensional flow of a gas particle system[R].DTIC Document,1959.

    [24] Kliegel J R.Gas particle nozzle flows[C]//Symposium (International) on Combustion,1963:811-826,0782- 0784.

    [25] Bailey W S,Nilson E N.Gas particle flow in an axisymmetric nozzle[J].ARS Journal,1961,31(6):793-798.

    [26] 常顯奇.固體火箭發(fā)動機燃燒室一維兩相常滯后流動[J].宇航學報,1983,4(4):40-52.

    [27] 常顯奇.顆粒尺寸對顆粒速度滯后數(shù)的影響[J].推進技術,1985,6(2):1-6.

    [28] Hassan H A.Exact solutions of gas-particle nozzle flows[J].AIAA Journal,1964,2(2):395-396.

    [29] Dellinger T C,Hassan H A.Analysis of gas-particle flows in rocket nozzles[J].Journal of Spacecraft and Rockets,1966,3(4):601-603.

    [30] Ma Y C,Fendell F.Constant-fractional-lag model for axisymmetric two-phase flow[J].Journal of Propulsion and Power,1991,7(5):700-707.

    [31] Hultberg J A,Soo S L.Metallized propellant and gas-solid suspension flow through nozzle in high energy systems[J].Astronautica Acta,1965,11(3):207-216.

    [32] 王慧玉,張遠君.全耦合的一維兩相噴管流的數(shù)值解[J].航空動力學報,1986,1(1):53-58;93.

    [33] Rannie W D.Perturbation analysis of one-dimensional heterogeneous flow in rocket nozzles[M].Progress in Astronautics and Rocketry.1962:117-144.

    [34] 方丁酉.兩相流動力學[M]. 長沙:國防科技大學出版社,1988:187;545.

    [35] 左羅克,霍夫曼.氣體動力學 下冊[M]. 北京:國防工業(yè)出版社,1984.

    [36] 方丁酉.一維兩相噴管流動[J].宇航學報,1982,3(3):25-40.

    [37] Chang I S.One-and two-phase nozzle flows[J].AIAA Journal,1980,18(12):1455-1461.

    [38] Forde M.Quasi-one-dimensional gas/particle nozzle flows with shock[J].AIAA Journal,1986,24(7):1196-1199.

    [39] Hoffman J D,Conway C C.Analysis of the flow of gas-particle mixtures in two-dimensional and axisymmetric nozzles[R].Dtic:DTIC Document,1962.

    [40] Kliegel J R,Nickerson G R.Flow of gas-particle mixtures in axially symmetric nozzles[M].Detonation and Two-Phase Flow,1962:173-194.

    [41] Hoffman J D,Lorenc S A.A parametric study of gas-particle flows in conical nozzles[J].AIAA Journal,1965,3(1):103-106.

    [42] Kliegel J R,Nickerson G R.Axisymmetric two-phase perfect gas performance program[R].NASA:MSC- 11774,1968.

    [43] Regan J F,Thompson H D.Two-dimensional analysis of transonic gas-particle flows in axisymmetric nozzles[J].Journal of Spacecraft and Rockets,1971,8(4):346-351.

    [44] Jacques L J,Seguin J a M.Two dimensional transonic two-phase flow in axisymmetric nozzles[C]//10th Propulsion Conference,1974.

    [45] Sharma M P,Crowe C T.A novel physico-computational model for quasi one-dimensional gas-particle flows[J].Journal of Fluids Engineering,1978,100(3):343-349.

    [46] Crowe C T,Sharma M P.The particle-source-in cell (PSI-CELL) model for gas-droplet flows[J].Journal of Fluids Engineering,1977,99(2):325-332.

    [47] 方丁酉.兩相噴管流動研究進展[J].推進技術,1982,3(1):33-42.

    [48] 方丁酉.兩相跨音速噴管流動[J].宇航學報,1987,8(3):43-51.

    [49] Hayashi A K,Matsuda M.Numerical study on gas-solid two-phase nozzle and jet flows[J].Memoirs of the Faculty of Engineering,Nagoya University,1988,40(2):351-362.

    [50] 侯曉,何洪慶.固體火箭噴管顆粒尺寸分級兩相跨音速流場計算[J].推進技術,1991,12(1):1-8.

    [51] 侯曉,何洪慶.固體火箭噴管兩相粘性跨音速流場計算[J].推進技術,1991,12(2):9-15;25.

    [52] 陳林泉.固體火箭發(fā)動機噴管兩相流動計算[J].固體火箭技術,1994,17(1):27-33.

    [53] 何洪慶,侯曉.在兩相粘性跨音速噴管流動中簿層方程的一種隱式求解方法[J].應用數(shù)學和力學,1994,15(4):303-312.

    [54] Hwang C J,Chang G C.Numerical study of gas-particle flow in a solid rocket nozzle[J].AIAA Journal,1988,26(6):682-689.

    [55] 傅德薰,馬延文.計算流體力學[M]. 北京:高等教育出版社,2002.

    [56] 顧璇.噴管內氣粒兩相流場的數(shù)值模擬[D].哈爾濱:哈爾濱工程大學,2004.

    [57] Ciucci A,Iaccarino G.Numerical investigation of 3D two-phase turbulent flows in solid rocket motors[C]//34th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit,1998.

    [58] 昌澤舟,Berlemont A.考慮顆粒間碰撞的氣固兩相流拉格朗日模擬[J].計算力學學報,2001,18(4):388- 392.

    [59] 淡林鵬,張振鵬.長尾噴管中粒子運動軌跡的數(shù)值模擬[J].航空動力學報,2003,18(2):258-263.

    [60] 于勇,劉淑艷.固體火箭發(fā)動機噴管氣固兩相流動的數(shù)值模擬[J].航空動力學報,2009,24(4):931-937.

    [61] 楊丹.固體火箭發(fā)動機氣-固兩相流的數(shù)值模擬[D].哈爾濱:哈爾濱工程大學,2006.

    [62] 劉叢林.金屬顆粒在氣固多相熱流場的動力學特性研究[D].哈爾濱:哈爾濱工程大學,2012.

    [63] 韓超.高金屬含量鎂基燃料水沖壓發(fā)動機穩(wěn)態(tài)燃燒機理研究[D].長沙:國防科技大學,2011.

    [64] Rudinger G.Some effects of finite particle volume on the dynamics of gas-particle mixtures[J].AIAA Journal,1965,3(7):1217-1222.

    [65] 于勇,張夏.用雙流體模型模擬超聲速氣固兩相流動[J].航空動力學報,2010,25(4):800-807.

    [66] Di Giacinto M,Sabetta F.Two-way coupling effects in dilute gas-particle flows[C]//ASME Winter Annual Meeting 1982,1982:ASME,New York,NY,USA.

    [67] Coakley T J,Champney J M.Numerical simulation of compressible,turbulent,two-phase flow[C]//18th Fluid Dynamics and Plasmadynamics and Lasers Conference,1985.

    [68] Mehta R C,Jayachandran T.A fast algorithm to solve viscous two-phase flow in an axisymmetric rocket nozzle[J].International Journal for Numerical Methods in Fluids,1998,26(5):501-517.

    [69] Perrell E.Two-phase CFD calculations with continuous distributions of particle sizes[C]//34th AIAA/ASME/ SAE/ASEE Joint Propulsion Conference and Exhibit,1998.

    [70] 高波,葉定友.旋轉條件下固體火箭發(fā)動機燃燒室氣-固兩相湍流流動數(shù)值模擬[J].固體火箭技術,1999,22(3):6-10.

    [71] 曾卓雄,姜培正.可壓稀相兩相流場的數(shù)值模擬[J].推進技術,2002,23(2):154-157.

    [72] Mirzaei M,Shadaram A.Numerical simulation of supersonic gas-particle flow using eulerian-eulerian approach[C]//14th AIAA/AHI International Space Planes and Hypersonics Systems Technologies Conference,November 6,2006:382-389.

    [73] 王運良,徐忠.氣固兩相流k-ε-kp-εp雙流體模型及應用[J].航空動力學報,1994,9(3):82-84;111.

    [74] 侯曉,何洪慶.粘性跨音速噴管流場計算[J].推進技術,1990,11(5):11-16;67.

    [75] 李江,肖育民.固體火箭發(fā)動機燃燒室兩相流粒子運動軌跡的實驗研究[J].西北工業(yè)大學學報,1999,17(1):15-18.

    [76] 肖育民,何國強.用RTR技術研究固體發(fā)動機燃燒室中粒子運動軌跡[J].推進技術,1997,15(5):36-45.

    [77] Carlson D J.Experimental determination of thermal lag in gas-particle nozzle flow[J].AIAA Journal,1962,32(7):1107-1109.

    [78] Carlson D J.Experimental determination of velocity lag in gas-particle nozzle flows[J].AIAA Journal,1965,3(2):354-357.

    [79] Neilson J H,Gilchrist A.An analytical and experimental investigation of the velocities of particles entrained by the gas flow in nozzles[J].Journal of Fluid Mechanics,1968,33(1):131-149.

    [80] 孫敏,方丁酉.二維噴管兩相流動實驗理論研究[J].航空學報,1988,9(11):572-576.

    [81] Bailey A B,Hiatt J.Sphere drag coefficients for a broad range of Mach and Reynolds numbers[J].AIAA Journal,1972,10(11):1436-1440.

    [82] Carlson D J,Hoglund R F.Particle drag and heat transfer in rocket nozzles[J].AIAA Journal,1964,2(11):1980-1984.

    [83] Hermsen R W.Review of particle drag models[C]//JANAF Performance Standardization Subcommittee 12th Meeting Minutes,1979:113.

    [84] 周彥煌,孫興長.氣-固兩相流中顆粒群的稠密度及粒度對阻力系數(shù)影響的研究[J].兵工學報,1984,5(4):11-15.

    [85] Xiao Y M,Amano R S.Aluminized composite solid propellant particle path in the combustion chamber of a solid rocket motor[C]//Sixth International Conference on Advances in Fluid Mechanics,2006:153-164.

    [86] Koike S,Takahashi H.Correction method for particle velocimetry data based on the stokes drag law[J].AIAA Journal,2007,45(11):2770-2777.

    [87] Kavanau L L,Drake Jr R M.Heat transfer from spheres to a rarefied gas in subsonic flow[R].DTIC Report:AD-1911,1953.

    [88] Koshmarov Y A,Svirshevskii S B.Heat transfer from a sphere in the intermediate dynamics region of a rarefied gas[J].Fluid Dynamics,1972,7(2):343-346.

    [89] Nelson H F,Fields J C.Heat transfer in two-phase solid-rocket plumes[J].Journal of Spacecraft and Rockets,1996,33(4):494-500.

    [90] Jung J Y,Brewster M Q.Radiative heat transfer analysis with molten Al2O3dispersion in solid rocket motors[J].Journal of Spacecraft and Rockets,2008,45(5):1021-1030.

    [91] 齊宏.彌散顆粒輻射反問題的理論與實驗研究[D].哈爾濱:哈爾濱工業(yè)大學,2008.

    [92] 李劍云,柳朝暉.求解微粒吸收指數(shù)的改進算法[J].工程熱物理學報,1996,17(1):116-120.

    [93] Hermsen R W.Aluminum oxide particle size for solid rocket motor performance prediction[C]//19th Aerospace Sciences Meeting,1981.

    [94] Sehgal R.An experimental investigation of a gas-particle system[R].DTIC Document,1962.

    [95] 婁永春.高溫稠密兩相流沖刷條件下絕熱層燒蝕實驗研究[D].西安:西北工業(yè)大學,2005.

    [96] 張超才,孫緒全.含鋁固體復合推進劑燃氣中Al2O3顆粒收集與測量的實驗研究[J].推進技術,1986,7(4):56-60.

    [97] 張勝敏,胡春波.固體火箭發(fā)動機噴管喉部凝相顆粒粒度分布實驗[J].推進技術,2012,33(2):245-248.

    [98] 劉洋,何國強.聚集狀態(tài)下凝相顆粒的收集與測量[J].推進技術,2005,26(5):477-480.

    [99] 李江,婁永春.聚集狀態(tài)對固體火箭發(fā)動機顆粒粒度分布的影響[J].固體火箭技術,2005,28(4):265-267.

    [100] Koo J H,Hirleman E D.Review of principles of optical techniques for particle size measurements,in:Progress in Astronautics and Aeronautics[M].AIAA,1996.

    [101] 張宏安,葉定友.固體發(fā)動機噴管出口凝相微粒粒度分布研究[J].固體火箭技術,2001,24(3):14-18.

    [102] Youngborg E D,Pruitt T E.Light-diffraction particle size measurements in small solid-propellant rockets[J].Journal of Propulsion and Power,1990,6(3):243-249.

    [103] 韓新波.固體火箭發(fā)動機潛入噴管背壁區(qū)流動研究[D].西安:西北工業(yè)大學,2002.

    [104] 陳劍,魏祥庚.固體火箭發(fā)動機長尾噴管燒蝕實驗研究[J].固體火箭技術,2010,33(1):34-35;40.

    [105] 賈林祥.含鋁固體推進劑火箭發(fā)動機噴管沉積的實驗與傳熱分析[J].推進技術,1985,6(1):1-10.

    [106] 趙湘恒,方丁酉.三氧化二鋁顆粒在噴管內的沉積[J].固體火箭技術,1988,11(4):55-65.

    [107] 夏盛勇.三氧化二鋁液滴碰撞機理及模型研究[D].西安:西北工業(yè)大學,2015.

    [108] Averin V S,Arkhipov V A.Effect of a sudden change in cross-sectional area of the solid rocket motor duct on coagulation of condensed particles[J].Combustion,Explosion and Shock Waves,2003,39(3):316-322.

    (編輯:崔賢彬)

    Research on the high solid concentration two-phase nozzle flow

    LIU Ping-an,WANG Liang,WANG Lu,WANG Ge,GAO Ye

    (College of Aerospace and Civil Engineering,Harbin Engineering University,Ha'erbin 150001,China)

    In order to select proper methods to study the high solid concentration two phase nozzle flow in highly metalized Solid Rocket Motor(SRM),the research advance of two phase flow in rocket nozzles was summarized in two aspects:the theoretical research and the experimental research.Their results and reference value for the high solid concentration nozzle flow research was discussed.It is noticed that the actual two phase nozzle flow lies between the two phase equilibrium flow and the two phase frozen flow.The iso-thermal method is suitable for the primary analysis of the high concentration two phase nozzle flow,as high particle content strengthens the interactions between the particle and the exhaust gas.The numerical methods can be generally categorized in two models:the trajectory model and the pseudo-fluid model,they are both suitable for the simulation of the high solid concentration nozzle flow.The experimental research are mainly used to verify the numerical results,and acquire relevant data for theoretical calculation,such as the particle size distribution,the drag coefficient,and the Nusselt number.The experiment can be greatly simplified if we only study some local features of the particle or the two phase nozzle flow-field,then the results can be expanded to the whole nozzle.By analyzing the research advance of the two phase nozzle flow,the reference for the high solid concentration nozzle flow research was provided.

    nozzle two phase flow;high solid concentration;solid rocket motor

    2016-07-12;

    2016-10-09。

    中央高校基本科研基金(HEUCFD1404);(HEUCFD1502)。

    劉平安( 1980—) ,男,副教授,研究方向為金屬燃料發(fā)動機。E-mail:liupingan631@126.com

    V438

    A

    1006-2793(2016)06-0735-11

    10.7673/j.issn.1006-2793.2016.06.001

    猜你喜歡
    粒徑方程發(fā)動機
    方程的再認識
    方程(組)的由來
    木屑粒徑對黑木耳栽培的影響試驗*
    圓的方程
    發(fā)動機空中起動包線擴展試飛組織與實施
    基于近場散射的顆粒粒徑分布測量
    Oslo結晶器晶體粒徑分布特征的CFD模擬
    新一代MTU2000發(fā)動機系列
    SAPO-56分子篩的形貌和粒徑控制
    新型1.5L-Eco-Boost發(fā)動機
    国产探花在线观看一区二区| 狂野欧美白嫩少妇大欣赏| 男人狂女人下面高潮的视频| 国产精品不卡视频一区二区| 男人舔奶头视频| 亚洲熟妇中文字幕五十中出| av卡一久久| 国产单亲对白刺激| 中文字幕av在线有码专区| 国产一级毛片七仙女欲春2| 久久精品国产亚洲av香蕉五月| 噜噜噜噜噜久久久久久91| 国产片特级美女逼逼视频| 日本在线视频免费播放| 三级男女做爰猛烈吃奶摸视频| 18禁黄网站禁片免费观看直播| 九九热线精品视视频播放| 精品一区二区三区视频在线| 午夜视频国产福利| 精品人妻偷拍中文字幕| 97在线视频观看| 久久午夜亚洲精品久久| 国产又黄又爽又无遮挡在线| 亚洲欧美精品自产自拍| 精品久久国产蜜桃| 三级毛片av免费| 国产精品电影一区二区三区| 久久热精品热| 黄色视频,在线免费观看| 国产麻豆成人av免费视频| 成人欧美大片| 日本五十路高清| 亚洲最大成人av| 美女 人体艺术 gogo| 尤物成人国产欧美一区二区三区| 网址你懂的国产日韩在线| 最新中文字幕久久久久| 日韩av在线大香蕉| 久久久久精品国产欧美久久久| 国语自产精品视频在线第100页| 日韩中字成人| 天堂网av新在线| 久久精品91蜜桃| 日韩欧美三级三区| 乱系列少妇在线播放| 免费人成在线观看视频色| 我的老师免费观看完整版| 国产高清视频在线观看网站| 少妇人妻精品综合一区二区 | 可以在线观看毛片的网站| 免费观看在线日韩| 亚洲av免费在线观看| 一级毛片我不卡| 国产精品永久免费网站| 三级男女做爰猛烈吃奶摸视频| 日韩制服骚丝袜av| 97热精品久久久久久| 日本三级黄在线观看| 日韩在线高清观看一区二区三区| 亚洲成人av在线免费| 在线天堂最新版资源| 大香蕉久久网| 在线观看66精品国产| 熟妇人妻久久中文字幕3abv| 欧美日韩国产亚洲二区| 赤兔流量卡办理| 青春草视频在线免费观看| 精品久久久久久久久久久久久| 色视频www国产| 99热6这里只有精品| 亚洲四区av| 午夜老司机福利剧场| 高清毛片免费观看视频网站| av卡一久久| 久久久久九九精品影院| 1024手机看黄色片| 亚洲av成人精品一区久久| 亚洲经典国产精华液单| 欧美精品国产亚洲| 国产精华一区二区三区| 国产又黄又爽又无遮挡在线| 日韩av不卡免费在线播放| 亚洲一级一片aⅴ在线观看| av在线播放精品| 精品乱码久久久久久99久播| 午夜影院日韩av| 日本成人三级电影网站| 亚洲人成网站在线播放欧美日韩| 免费观看在线日韩| 亚洲电影在线观看av| 亚洲国产精品合色在线| 欧美日韩国产亚洲二区| 少妇熟女aⅴ在线视频| 欧美不卡视频在线免费观看| 在线观看免费视频日本深夜| 小说图片视频综合网站| 在线观看免费视频日本深夜| 午夜福利在线观看免费完整高清在 | 亚洲精品日韩av片在线观看| 国产 一区精品| 草草在线视频免费看| 亚洲精品色激情综合| 天堂动漫精品| 长腿黑丝高跟| 长腿黑丝高跟| 国产高清有码在线观看视频| 一边摸一边抽搐一进一小说| 亚洲精品影视一区二区三区av| 免费观看人在逋| 十八禁网站免费在线| 美女黄网站色视频| 五月玫瑰六月丁香| 国产国拍精品亚洲av在线观看| 黑人高潮一二区| 亚洲av美国av| 免费av毛片视频| 精品乱码久久久久久99久播| 精品久久久久久久久av| 国产精品女同一区二区软件| 一夜夜www| 最近手机中文字幕大全| 日本欧美国产在线视频| 国产中年淑女户外野战色| a级毛色黄片| 久久国内精品自在自线图片| 午夜激情欧美在线| 精品日产1卡2卡| 午夜久久久久精精品| 人人妻人人看人人澡| 国产精品亚洲美女久久久| 伦理电影大哥的女人| 婷婷精品国产亚洲av| 日本成人三级电影网站| 免费观看在线日韩| 国产视频一区二区在线看| 色综合色国产| 伦精品一区二区三区| 免费一级毛片在线播放高清视频| 免费一级毛片在线播放高清视频| 波多野结衣巨乳人妻| 中文字幕av在线有码专区| 成人一区二区视频在线观看| 欧美最黄视频在线播放免费| 欧美最黄视频在线播放免费| 亚洲,欧美,日韩| 国语自产精品视频在线第100页| 成人亚洲欧美一区二区av| 久久久久免费精品人妻一区二区| 国国产精品蜜臀av免费| 18禁在线播放成人免费| 国产片特级美女逼逼视频| 我的女老师完整版在线观看| 淫秽高清视频在线观看| 日韩欧美免费精品| 69人妻影院| 久久久久九九精品影院| 日本撒尿小便嘘嘘汇集6| 内地一区二区视频在线| 久久精品国产亚洲av天美| 中文字幕久久专区| 国产成人freesex在线 | 国产精品综合久久久久久久免费| 国产一区二区在线av高清观看| 欧美性猛交╳xxx乱大交人| 久久天躁狠狠躁夜夜2o2o| 中文字幕精品亚洲无线码一区| 级片在线观看| 97超碰精品成人国产| 久久精品国产亚洲av涩爱 | 99热精品在线国产| 国产老妇女一区| 大香蕉久久网| 性插视频无遮挡在线免费观看| 亚洲乱码一区二区免费版| 国产三级在线视频| 成人三级黄色视频| 国产高清激情床上av| АⅤ资源中文在线天堂| 69av精品久久久久久| 国产精华一区二区三区| 韩国av在线不卡| 亚洲精品影视一区二区三区av| 在线观看美女被高潮喷水网站| 亚洲aⅴ乱码一区二区在线播放| 国产精品人妻久久久影院| 国产一区亚洲一区在线观看| a级毛片免费高清观看在线播放| 日韩强制内射视频| 最近中文字幕高清免费大全6| 两个人视频免费观看高清| 久久久久九九精品影院| 人妻少妇偷人精品九色| 在线观看午夜福利视频| 欧美色视频一区免费| 免费在线观看成人毛片| 国产精品美女特级片免费视频播放器| 亚洲欧美日韩高清专用| 夜夜爽天天搞| 99久久精品国产国产毛片| 成人av在线播放网站| 亚洲人与动物交配视频| 精品日产1卡2卡| 天天一区二区日本电影三级| 日韩制服骚丝袜av| 日韩欧美免费精品| 日本三级黄在线观看| 免费av不卡在线播放| 亚洲综合色惰| 国产黄色视频一区二区在线观看 | 国产av不卡久久| 99在线视频只有这里精品首页| 中文字幕久久专区| 国产一区二区三区av在线 | 男插女下体视频免费在线播放| 国产乱人偷精品视频| 日韩在线高清观看一区二区三区| 在线播放无遮挡| 亚洲精华国产精华液的使用体验 | 欧美zozozo另类| 欧美xxxx黑人xx丫x性爽| 在现免费观看毛片| 免费看av在线观看网站| 一卡2卡三卡四卡精品乱码亚洲| 成年免费大片在线观看| 日本撒尿小便嘘嘘汇集6| 一本精品99久久精品77| 亚洲国产日韩欧美精品在线观看| 变态另类成人亚洲欧美熟女| 日本五十路高清| 亚洲自偷自拍三级| 狠狠狠狠99中文字幕| 在线免费十八禁| 精品日产1卡2卡| 美女被艹到高潮喷水动态| 日韩精品中文字幕看吧| 老熟妇仑乱视频hdxx| 搡女人真爽免费视频火全软件 | 欧美另类亚洲清纯唯美| 国产亚洲精品综合一区在线观看| 免费人成在线观看视频色| 亚洲中文日韩欧美视频| 两性午夜刺激爽爽歪歪视频在线观看| 国语自产精品视频在线第100页| 亚洲精品日韩av片在线观看| 久久久久久久久久成人| 日韩欧美 国产精品| 少妇人妻一区二区三区视频| 亚洲av.av天堂| 少妇的逼好多水| 国产精品国产三级国产av玫瑰| 欧美精品国产亚洲| 久久精品影院6| 一级毛片电影观看 | 黄色视频,在线免费观看| 九色成人免费人妻av| 听说在线观看完整版免费高清| 天天躁日日操中文字幕| 成人性生交大片免费视频hd| 麻豆乱淫一区二区| 欧美极品一区二区三区四区| 直男gayav资源| 亚洲不卡免费看| 国产 一区 欧美 日韩| 亚洲av中文字字幕乱码综合| 国产亚洲精品av在线| 少妇熟女aⅴ在线视频| 熟女电影av网| 在线观看66精品国产| 看免费成人av毛片| 18+在线观看网站| 天堂av国产一区二区熟女人妻| 免费看光身美女| 亚洲av二区三区四区| .国产精品久久| 三级经典国产精品| 婷婷亚洲欧美| 黄色配什么色好看| 欧美人与善性xxx| 国产精品不卡视频一区二区| 99热精品在线国产| 成人欧美大片| 亚洲在线观看片| 观看免费一级毛片| 午夜免费男女啪啪视频观看 | 久久热精品热| 一进一出抽搐gif免费好疼| 蜜臀久久99精品久久宅男| 亚洲国产精品sss在线观看| 久久99热这里只有精品18| 欧美一区二区国产精品久久精品| 18禁在线播放成人免费| 18禁裸乳无遮挡免费网站照片| 亚洲国产精品成人综合色| 久久久国产成人免费| 美女高潮的动态| 丝袜喷水一区| 日韩欧美在线乱码| 三级经典国产精品| 在线观看免费视频日本深夜| 免费在线观看影片大全网站| 亚洲自拍偷在线| 国产精品国产三级国产av玫瑰| 激情 狠狠 欧美| 岛国在线免费视频观看| 久久天躁狠狠躁夜夜2o2o| 精品欧美国产一区二区三| av女优亚洲男人天堂| 精品无人区乱码1区二区| 晚上一个人看的免费电影| 国产在视频线在精品| 欧美潮喷喷水| 男女那种视频在线观看| 免费av观看视频| av天堂中文字幕网| 91久久精品电影网| 99国产极品粉嫩在线观看| 精华霜和精华液先用哪个| 菩萨蛮人人尽说江南好唐韦庄 | 99精品在免费线老司机午夜| 中文字幕av在线有码专区| 久久久久国产精品人妻aⅴ院| 12—13女人毛片做爰片一| 22中文网久久字幕| 日韩av不卡免费在线播放| 婷婷精品国产亚洲av在线| 人妻夜夜爽99麻豆av| 99久久久亚洲精品蜜臀av| 欧美日韩在线观看h| 欧美日韩在线观看h| 亚洲欧美成人综合另类久久久 | 国产精品人妻久久久影院| 中文字幕熟女人妻在线| 久久久久九九精品影院| 国产中年淑女户外野战色| 色尼玛亚洲综合影院| 日日摸夜夜添夜夜添小说| 伦理电影大哥的女人| 欧美一区二区精品小视频在线| 最新中文字幕久久久久| 久久久久久伊人网av| 国产一区二区三区在线臀色熟女| 精品久久久久久久末码| 最近的中文字幕免费完整| 乱系列少妇在线播放| 日本免费a在线| 18禁在线无遮挡免费观看视频 | 成人精品一区二区免费| 中文资源天堂在线| 麻豆一二三区av精品| 色综合色国产| 亚洲成人久久性| 村上凉子中文字幕在线| 久久中文看片网| 成人漫画全彩无遮挡| 非洲黑人性xxxx精品又粗又长| 又爽又黄无遮挡网站| 男女那种视频在线观看| 亚洲专区国产一区二区| 久久久成人免费电影| 亚洲精品久久国产高清桃花| 国产av一区在线观看免费| 高清午夜精品一区二区三区 | 国产白丝娇喘喷水9色精品| 大型黄色视频在线免费观看| 波多野结衣高清作品| 午夜免费激情av| 老女人水多毛片| 亚洲成人av在线免费| 久久这里只有精品中国| 久久久久久久久久黄片| 中国国产av一级| 精品一区二区三区视频在线| 久久精品国产鲁丝片午夜精品| 熟妇人妻久久中文字幕3abv| 不卡一级毛片| 国产成人91sexporn| 天美传媒精品一区二区| 乱系列少妇在线播放| 特大巨黑吊av在线直播| 一本精品99久久精品77| 人妻制服诱惑在线中文字幕| 免费在线观看影片大全网站| 久久精品夜夜夜夜夜久久蜜豆| 欧美一区二区国产精品久久精品| 性插视频无遮挡在线免费观看| 午夜激情欧美在线| 啦啦啦观看免费观看视频高清| 一级黄色大片毛片| 成人三级黄色视频| 女生性感内裤真人,穿戴方法视频| 亚洲aⅴ乱码一区二区在线播放| 麻豆成人午夜福利视频| 一卡2卡三卡四卡精品乱码亚洲| 国产探花在线观看一区二区| 亚洲av中文字字幕乱码综合| 国产精品1区2区在线观看.| 国产成人精品久久久久久| 欧美日韩在线观看h| a级毛片免费高清观看在线播放| 99久久无色码亚洲精品果冻| 成年女人毛片免费观看观看9| 国产成人一区二区在线| 日日啪夜夜撸| 国产av不卡久久| 精品人妻熟女av久视频| 免费av观看视频| 亚洲欧美日韩无卡精品| 亚洲aⅴ乱码一区二区在线播放| 在线免费十八禁| 在线播放无遮挡| 2021天堂中文幕一二区在线观| 中文字幕熟女人妻在线| aaaaa片日本免费| 美女被艹到高潮喷水动态| 国产一级毛片七仙女欲春2| 久久久久久久久久黄片| 男女边吃奶边做爰视频| 一本一本综合久久| 亚洲av成人精品一区久久| 久久久欧美国产精品| www日本黄色视频网| 午夜精品在线福利| 欧美激情国产日韩精品一区| 日日撸夜夜添| 久久人人爽人人片av| 国产精品1区2区在线观看.| 精品久久久久久久久av| 精品午夜福利视频在线观看一区| 熟女电影av网| 欧美性猛交╳xxx乱大交人| 亚洲图色成人| 91午夜精品亚洲一区二区三区| 日本爱情动作片www.在线观看 | 99热精品在线国产| 精品不卡国产一区二区三区| 九色成人免费人妻av| 一区二区三区四区激情视频 | 国产视频一区二区在线看| 国产真实伦视频高清在线观看| 免费人成视频x8x8入口观看| 夜夜夜夜夜久久久久| 中文字幕av在线有码专区| 久久热精品热| 天堂动漫精品| 久久久久九九精品影院| 老师上课跳d突然被开到最大视频| 免费av不卡在线播放| 成人国产麻豆网| 色播亚洲综合网| 亚洲精品乱码久久久v下载方式| 亚洲成av人片在线播放无| 蜜臀久久99精品久久宅男| 国产亚洲精品av在线| 中文字幕av在线有码专区| 日韩一本色道免费dvd| 一个人免费在线观看电影| 日本一本二区三区精品| 精品一区二区免费观看| 国产av一区在线观看免费| 国语自产精品视频在线第100页| 日本精品一区二区三区蜜桃| 久久精品国产亚洲av天美| 日本五十路高清| 久久人人爽人人片av| 日产精品乱码卡一卡2卡三| 好男人在线观看高清免费视频| 两性午夜刺激爽爽歪歪视频在线观看| 天堂动漫精品| 99热这里只有是精品50| 亚洲av免费在线观看| 中文在线观看免费www的网站| 成人高潮视频无遮挡免费网站| 熟妇人妻久久中文字幕3abv| 亚洲国产高清在线一区二区三| 我的女老师完整版在线观看| 老师上课跳d突然被开到最大视频| 真人做人爱边吃奶动态| 欧美三级亚洲精品| 在线观看66精品国产| 不卡视频在线观看欧美| 亚洲成人久久性| 国产三级中文精品| 精品乱码久久久久久99久播| 91久久精品国产一区二区成人| 极品教师在线视频| 日本熟妇午夜| 国产免费一级a男人的天堂| 日本色播在线视频| 丰满乱子伦码专区| 一本一本综合久久| 亚洲天堂国产精品一区在线| 亚洲欧美成人精品一区二区| 中文字幕av在线有码专区| 国产爱豆传媒在线观看| 嫩草影院精品99| 美女 人体艺术 gogo| 欧美日韩国产亚洲二区| 床上黄色一级片| 老女人水多毛片| 亚洲国产精品sss在线观看| 51国产日韩欧美| 在线a可以看的网站| 日本免费a在线| 成人精品一区二区免费| 免费不卡的大黄色大毛片视频在线观看 | 国产探花在线观看一区二区| 少妇丰满av| 国产单亲对白刺激| 日日摸夜夜添夜夜添av毛片| 日本精品一区二区三区蜜桃| 国产精品永久免费网站| 高清午夜精品一区二区三区 | 在线观看66精品国产| 色播亚洲综合网| 午夜久久久久精精品| 欧美日韩精品成人综合77777| 美女大奶头视频| 国产麻豆成人av免费视频| 自拍偷自拍亚洲精品老妇| 欧洲精品卡2卡3卡4卡5卡区| 亚洲内射少妇av| 欧美绝顶高潮抽搐喷水| 色播亚洲综合网| 成人毛片a级毛片在线播放| 麻豆乱淫一区二区| 国产伦精品一区二区三区视频9| 人妻丰满熟妇av一区二区三区| or卡值多少钱| 精品久久久久久久久亚洲| 精品日产1卡2卡| 黑人高潮一二区| 搡老熟女国产l中国老女人| 久久天躁狠狠躁夜夜2o2o| 国产爱豆传媒在线观看| 内地一区二区视频在线| 精品久久久久久久久av| 国内精品一区二区在线观看| 在现免费观看毛片| 国产午夜精品久久久久久一区二区三区 | 免费观看精品视频网站| 欧美日韩乱码在线| 亚洲国产高清在线一区二区三| 99热全是精品| 国产精品久久久久久精品电影| 看片在线看免费视频| 日日摸夜夜添夜夜爱| 1000部很黄的大片| 欧美成人精品欧美一级黄| 亚洲精品色激情综合| 国产成人一区二区在线| 免费人成视频x8x8入口观看| 国产精品综合久久久久久久免费| 精品日产1卡2卡| 长腿黑丝高跟| 看免费成人av毛片| 一本一本综合久久| 日产精品乱码卡一卡2卡三| 欧美成人免费av一区二区三区| 国产精品人妻久久久影院| 久久亚洲国产成人精品v| 综合色丁香网| 精品一区二区三区人妻视频| 极品教师在线视频| 长腿黑丝高跟| 两性午夜刺激爽爽歪歪视频在线观看| 久久国内精品自在自线图片| 亚洲成人精品中文字幕电影| 亚洲第一电影网av| 激情 狠狠 欧美| 在线观看免费视频日本深夜| 日本撒尿小便嘘嘘汇集6| 日韩制服骚丝袜av| 国产精品人妻久久久久久| 午夜福利在线观看吧| 久久久午夜欧美精品| 尤物成人国产欧美一区二区三区| 露出奶头的视频| 日本成人三级电影网站| 国产日本99.免费观看| av免费在线看不卡| 不卡视频在线观看欧美| 能在线免费观看的黄片| 丰满的人妻完整版| 特级一级黄色大片| 好男人在线观看高清免费视频| 在线观看一区二区三区| 午夜视频国产福利| 中国美女看黄片| 国产精品无大码| 大又大粗又爽又黄少妇毛片口| 亚洲av熟女| 日韩av不卡免费在线播放| 国内久久婷婷六月综合欲色啪| 中文字幕熟女人妻在线| 在线天堂最新版资源| 成人鲁丝片一二三区免费| 天堂动漫精品| 国产毛片a区久久久久| 亚洲人与动物交配视频| 国产成人91sexporn| av在线亚洲专区| 一本一本综合久久| 亚洲美女视频黄频| 国产精品不卡视频一区二区| 美女cb高潮喷水在线观看| 97人妻精品一区二区三区麻豆| 99热精品在线国产| 一边摸一边抽搐一进一小说| 午夜影院日韩av| 黄色欧美视频在线观看| 69av精品久久久久久| 男人舔奶头视频| 国产色爽女视频免费观看| 两个人视频免费观看高清| 婷婷精品国产亚洲av在线|