999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

電弧加熱器湍流平板試驗流場計算分析

2017-05-24 14:46:22中國空氣動力研究與發展中心四川綿陽62000四川大學制造科學與工程學院成都60065
實驗流體力學 2017年2期
關鍵詞:模型

羅 躍, 周 瑋, 楊 鴻, 陳 衛(. 中國空氣動力研究與發展中心, 四川 綿陽 62000; 2. 四川大學 制造科學與工程學院, 成都 60065)

電弧加熱器湍流平板試驗流場計算分析

羅 躍1,2,*, 周 瑋1, 楊 鴻1, 陳 衛1
(1. 中國空氣動力研究與發展中心, 四川 綿陽 621000; 2. 四川大學 制造科學與工程學院, 成都 610065)

近年來新發展的變迎角方法增加了電弧加熱器湍流平板試驗的可控參數,大幅提升了試驗模擬能力,但同時也增加了狀態調試的復雜性。本文基于FLUENT軟件,采用CFD方法分析了不同試驗工況下的流場特性。通過與典型試驗狀態的測試數據比較,驗證了計算方法的有效性,定量地給出不同弧室壓力和模型迎角條件下校測平板表面主要氣動熱參數分布,結合流場激波反射、干擾和附面層分離等現象,總結分析了試驗控制參數的影響規律,并根據研究結果提出了試驗的基本原則。

電弧加熱器;平板試驗;流場特性;數值模擬;校測

0 引 言

電弧加熱器湍流平板試驗技術是針對高超聲速飛行器熱防護材料考核的一種地面試驗技術,主要用于表面為平板或以平板為基板的防熱材料試驗[1]。由電弧加熱器產生的高溫氣流,經轉接過渡段時形成層流邊界層并在噴管喉道附近轉捩為湍流邊界層,在超聲速矩形型面噴管出口處,與氣流有一定迎角地放置平板試件,兩者密接齊平無縫隙。試件上的邊界層是噴管型面壁上邊界層的自然延續,從而在試件上得到充分發展的湍流邊界層流動,而試件迎角造成的前緣斜激波,可以提高試件上模擬的參數[2]。這種試驗方法能夠模擬試件飛行條件下的表面壓力、熱流和剪切力等氣動參數,可以承擔大面積防熱材料燒蝕試驗研究、天線窗口周圍局部燒蝕試驗研究、控制翼局部燒蝕試驗研究以及槽孔和縫隙局部燒蝕試驗研究等,是電弧加熱設備上最常用的試驗方法之一。

在正式試驗之前,需要對電弧加熱器和試件當地流場進行調試和校測,以確認流場特性和當地參數分布,這也是確保試驗質量的關鍵環節。這個過程通常是由外形與試件相同的傳感器校測模型來完成,需要對各項控制參數進行反復調試和匹配。近年來隨著CFD技術的發展和其明顯的優勢,國外相關研究機構已越來越多地將CFD計算應用到電弧加熱器流場預測中,它能夠獲得校測試驗難以測量的局部區域,以及整個流場其他參數的分布,可以指導試驗調試的實施,減少流場校測環節的探索過程[3]。NASA的Ames研究中心基于N-S方程開發了CFD程序DPLR,應用于大量的電弧加熱設備試驗中[3-11],并將其作為地面試驗的固有流程和標準配置[4]。應用分析顯示,即便DPLR對試件表面熱流的預測不確定度較大,仍不失為一種有效的輔助分析工具。文獻[5]對改建的平板試驗設備TPTF進行不同迎角的流場模擬,與校測試驗數據誤差小于15%,獲得了多種工況下平板試驗的流場特性。文獻[6]利用ADSI程序對已有數據庫進行插值,以預測數據庫以外的流場狀態,分析結論指出基于CFD數據庫插值的結果與實際測量偏差小于10%,且優于以試驗數據庫進行插值的結果。國內CFD技術發展較快,但在電弧加熱試驗中并未大量廣泛應用,僅限于單一狀態或問題的模擬。張友華、涂建強等對湍流導管試驗流場進行了工程估算和數值模擬[12-13],以指導試驗方法。李媛、于勝春等利用商業軟件對超聲速噴管自由射流場進行了分析[14-15],楊鴻對湍流平板燒蝕試驗流場的工程估算、數值模擬和測試數據進行了比較,證明了CFD計算對試驗前的狀態估算是完全可行的[16-17],隆永勝利用不同氣體狀態方程計算了半橢圓噴管流場,得到了湍流和真實氣體效應對流場結構的影響[18]。

隨著高超聲速飛行器的發展,迫切需要進行更精確的防熱設計,防熱材料的地面考核試驗也要求盡可能模擬到更寬范圍的彈道參數變化。因此,近年來在電弧加熱器平板試驗上發展了變迎角試驗技術,即在同一車次中可連續變化試件迎角,以增大表面氣動參數的匹配范圍。這種方法提升了試驗的模擬能力,但同時也增加了狀態調試的復雜性,使得調試和校測環節占據整個試驗的大部分周期和費用。本文基于FLUENT軟件,采用CFD方法對電弧加熱器湍流平板試驗的不同壓力和迎角流場進行模擬和分析,研究試件表面氣動參數特性分布,以期為試驗的實施提供參考,幫助前期的試驗策劃。

1 計算方法

1.1 控制方程與離散方法

計算采用的控制方程為一般直角坐標系中無量綱可壓縮N-S方程[19]:

對控制方程采用二階迎風精度,基于微元中心有限體積法(控制容積法)進行空間離散求解,基于密度的隱式耦合求解方法,即同時求解連續方程、動量方程、能量方程及輸運方程的耦合方程組,然后逐一地求解湍流標量方程。

1.2 物理模型及邊界條件

電弧湍流平板試驗的原理如圖1所示,據此建立的計算域如圖2所示。噴管內流區和射流外流區組成整個流體計算域,噴管為出口尺寸100mm×32mm、名義馬赫數2.3的二元矩形超聲速型面噴管,平板模型長度為150mm。噴管左邊入口的高溫高壓空氣經內部型面加速后,作用到與噴管出口齊平安裝的平板模型上,其后半自由地射流到大氣環境中。為了減小噴管出口射流對邊界的影響,將與噴管出口臨近的邊界左移20倍特征長度的距離,并設置為壓力遠場條件,上下空間遠場邊界距離射流區60倍特征長度,噴管入口采用壓力入口條件,右邊出口采用壓力出口條件,模型表面采用固定無滑移等溫壁條件。

1.3 網格劃分

為了保證計算精度和收斂速度,整個計算域采用四邊形結構網格,保證網格物面正交并向物面和變化

梯度較大的區域按指數規律加密。由于涉及模型表面附面層加熱,近壁面的網格劃分尤為重要,除了確保足夠的網格數量和密度之外,需要根據采用的非平衡壁面處理方法,將壁面第一層無量綱網格距離Yplus控制在30~100之間[20]。Yplus按照如下公式估算:

式中:τw為壁面剪切力,Δy在非平衡壁面處理方法中為四分之一近壁面第1和第2個網格點之間距離。劃分的網格結構如圖3所示。

1.4 湍流模型

慢阻肺診斷標準:(1)有或無呼吸困難、慢性咳嗽、咳痰癥狀;(2)有吸煙、生物燃料、粉塵或其它風險因素暴露史;(3)肺功能提示存在明顯的氣流受限(舒張后FEV1/FVC<0.7)。

湍流模型采用兩方程的RNGk-ε湍流模型。RNGk-ε模型的輸運方程采用與標準k-ε模型相同的形式, 然而它的導出依據了更為嚴格的統計技術——重整化群理論, 顯著提高了解決快速變形流動問題的精度。在該模型中, 模式常數由重整化群理論給出:C1ε= 1.42,C2ε=1.68。

2 測試試驗及結果對比

為驗證計算方法,進行了有限車次的湍流平板流場校測試驗。試驗在CARDC(中國空氣動力研究與發展中心)的20MW管式電弧加熱器上進行。試驗方法及裝置布置如圖1及前文所述。

利用平板校測模型(見圖4)測量表面熱流和壓力分布。表面壓力是在需要的位置上,法向方向開Φ1mm的測壓孔,并連接微型壓力傳感器進行測量。測量熱流的量熱塞嵌入校測板中,并與表面平齊,采取措施隔離量熱塞與校測板之間的傳熱。在量熱塞瞬態感受流場熱流的過程中,可以忽略量熱塞側向的傳熱,認為表面熱流只向其法向方向進行一維熱傳導,根據一維熱傳導理論,表面熱流可以近似為:

因此,通過測量量熱塞的溫升速率即可獲得表面熱流。校測板上測點位置分布如圖5所示,通過調整校測板的位置,也可改變測點在流場中的相對位置。

電弧加熱器點火前,校測板處于流場外部,待電弧加熱器運行且狀態穩定后,送進機構快速將校測板送入指定的流場位置進行測量,校測板通常在流場中停留1~2s,獲得數據后再由送進機構快速送出流場完成測試,圖6為試驗過程中的流場照片。

驗證試驗選用了2個典型試驗狀態condA和condB進行測試并與數值模擬結果進行比較。condA的工況為:總溫1700K,總壓1.4MPa,迎角0°;condB的工況為:總溫3200K,總壓1.3MPa,迎角19°。圖7和8顯示了2個狀態下水冷校測板表面的熱流、壓力的計算結果和測試數據。從圖中可以看出,測試和計算結果在分布趨勢上有著較好的一致性,相對誤差小于20%,說明文中的計算方法能夠有效預測試驗條件下校測板的表面氣動熱參數。

Fig.7 Comparison of computed surface heat flux of the calibration plate with the test data

Fig.8 Comparison of computed surface pressure of the calibration plate with the test data

3 試驗流場模擬結果及分析

流場模擬分為2個部分:第1部分給定入口氣流狀態,改變模型迎角進行模擬,為了凸顯模型迎角的影響,設定狀態使噴管出口靜壓約為0.1MPa即環境大氣壓的氣流狀態;第2部分給定模型迎角,改變入口氣流狀態進行模擬。

第1部分(定入口氣流狀態)計算:給定入口總壓1.4MPa,總溫1500K的入口條件,分別對模型迎角0°、3°、6°、9°、12°、15°和18°的情況進行了穩態計算。圖9和10為不同模型迎角情況下噴管出口附近氣流馬赫數和靜壓云圖。

從計算結果可以看出模型無迎角時,出射氣流在模型上均勻發展,流場內部結構一致。在有迎角時,由于流動方向強迫改變,出射氣流在模型頭部形成斜激波,使得模型位置馬赫數減小,壓力增加,激波在氣流邊界上反射為膨脹波,觸及模型壁面后又反射為膨脹波。隨著模型迎角的增大,斜激波角度和強度隨之增大,反射波系更為明顯,并在模型后形成較強的激波串。可見計算方法較好地捕捉了超聲速氣流激波/激波干擾的現象,波系明顯,且波系結構符合基本的氣動力學規律和試驗結果。

利用修正牛頓理論計算不同迎角校測板前部10mm處表面壓力,并與數值計算對比(見圖11)。2種方法對壓力的計算一致性較好,相對誤差小于5%。

Fig.11 CFD computed and inferred surface pressure of the calibration plate with various angle

為了了解流場對模型表面參數的影響,將不同迎角情況下模型表面冷壁熱流和壓力分布列于圖12和13。從圖中可以看出,0°迎角時,模型表面熱流和壓力沿氣流方向緩慢下降,下降幅度較小,基本可認為整個模型表面參數是均勻的。有迎角時,由于前部斜激波的作用,模型表面熱流和壓力升高,且受到反射波的影響,尾部壁面參數出現驟降。隨著迎角的增大,模型表面參數隨之升高,然而激波反射位置逐漸前移,使得模型表面的均勻區逐漸減小。由此可以看出,通過增大模型迎角,能夠明顯提升表面參數,從0°~18°,均勻區的冷壁熱流大約增加至原來的260%。然而這是以犧牲均勻區范圍為代價的,18°迎角的均勻區只有約70mm。需要注意的是,此次計算選用的是較早期的小口徑噴管,對于后期配置的大口徑噴管,這種情況會改善很多,甚至不會出現,這取決于模型長度和噴管出口高度。

Fig.12 Computed surface heat flux of the calibration plate with various angle

Fig.13 Computed surface pressure of the calibration plate with various angle

值得注意的是,從15°迎角開始,模型前部出現了一段低熱流區域,在18°迎角時這種現象更為明顯,區域長度約為25mm。分析認為這是由于大迎角情況下,氣流在模型前部的拐角處出現了附面層分離,并且從圖14中的流線能夠看出噴管與模型相連的拐角位置的確沒有氣流貼附壁面。在更大迎角情況下,由于附面層分離再附而導致的低參數區域可能更長,鑒于此,在進行電弧加熱器湍流平板試驗時,模型迎角不宜超過20°,這也給本文的試驗方法提供了參考和指導。

第2部分(定迎角)計算:為了捕捉波系干擾,設定模型的迎角為9°,入口氣流總溫為1500K,分別對入口氣流總壓為0.8、1.0、1.4、1.6、1.8和2.0MPa的流場進行穩態數值計算。圖15和16為不同入口總壓噴管出口附近氣流馬赫數和靜壓云圖。從計算結果可以看出,入口壓力1.4MPa時,由于噴管出口靜壓約為0.1MPa,與環境壓力相近,噴管口只有由模型引起的斜激波,當入口壓力大于1.4MPa時,噴管口上緣出現膨脹波,使氣流加速、壓力減小,氣流邊界膨脹。膨脹波向下與激波交匯,由于此區域氣流受激波影響壓力較大,膨脹波經交匯點后依然是膨脹波,并在模型壁面反射為新的膨脹波,此后又與氣流邊界上反射的膨脹波交匯。波系交匯、反射以及干擾形成菱形區,隨著入口壓力的增大,這種現象更加明顯。當入口壓力小于1.4MPa時,噴管口上緣出現激波,使氣流壓力增加,氣流邊界收縮。上下兩條激波交匯,并分別在模型壁面和氣流邊界上反射為激波和膨脹波。由于氣流的上邊界收縮,這種情況的菱形區域更加復雜。可見計算方法有效捕捉了超聲速氣流波系交匯、干擾和反射現象。

Fig.15 Computed flowfield Mach number contours of various inlet total pressure

將不同入口壓力情況下模型表面冷壁熱流和壓力分布列于圖17和18中。從圖中可以看出,當入口壓力大于1.4MPa時,模型表面參數能夠得到提升,但由于膨脹波的反射干擾,反射點之后的表面參數驟降,均勻區范圍大幅減小至50mm,而且若以此方法來提升表面參數,將犧牲大量的加熱器功率消耗。隨著入口壓力的增加,均勻區范圍變化很小,但反射點之后的表面參數下降幅度逐漸增大。當入口壓力小于1.4MPa時,模型前部表面參數降低,但由于2道激波的反射干擾,表面參數在反射點之后大幅度跳躍變化,且均勻區范圍只有不到40mm,無法滿足試驗要求。鑒于此,在進行電弧加熱器湍流平板試驗時,要求噴管出口靜壓大于等于環境壓力。

Fig.17 Computed surface heat flux of the calibration plate with various total pressure

Fig.18 Computed surface pressure of the calibration plate with various total pressure

4 結 論

(1) 計算方法能夠正確捕捉流場激波反射、干擾和附面層分離等現象,典型試驗狀態的計算結果與試驗數據相符,可以有效預測試驗條件下的流場結構和參數。

(2) 湍流平板試驗中,改變入口氣流狀態和迎角均可以改變模型表面參數,2種方法各有其局限性。前者的局限性在于均勻區較小,且耗費大量的電弧功率,優勢在于提升空間大,且隨著參數的提高均勻區變化小;后者的局限性在于隨著迎角的增加,均勻區也會逐漸減小(比前者大),且過大的迎角將導致附面層分離,優勢在于變化效率高,控制靈活經濟,而且對于大尺寸噴管情況,模型表面參數可以避免波系干擾,均勻性較好。當然,對于不同的具體情況,需要結合2種方法以滿足試驗的需求。

(3) 根據試件表面參數分布特性均勻性要求,噴管出口靜壓應約等于或大于環境壓力,模型迎角大于15°時,應充分考慮附面層分離的影響。

(4) 本文僅針對電弧湍流平板試驗的中低溫度和二維流場條件,后續將進一步開展涉及真實氣體效應的高溫氣流三維流動研究。

[1]中國空氣動力研究與發展中心. GJB7050-2010導彈防熱材料電弧加熱器湍流平板試驗方法[S]. 北京:總裝備部軍標出版發行部, 2010: 1.China Aerodynamics Research & Development Center. GJB 7050-2010 Plat test method in turbulent flow for thermal protection materials of missile at arc-heater[S]. Beijing: Military Standard Publication and Distribution Section of General Armament Department, 2010: 1

[2]張志成. 高超聲速氣動熱和熱防護[M]. 北京: 國防工業出版社, 2003: 278-279.Zhang Z C. Hypersonic aerothermodynamic and thermal protection[M]. Beijing: National Defence Industry Press, 2003: 278-279

[3]Tahir G, Kristina S. Computational analysis of arc-jet wedge calibration tests in IHF 6-inch conical nozzle[R]. AIAA-2009-1348, 2009

[4]Tahir G, Antonella I A, Kristina A S. Computational simulations of panel test facility flow: compression-pad arc-jet tests[R]. AIAA-2011-3635, 2011

[5]John A B, Tahir G, Frank C, et al. Calibration of the truncated panel test arc-jet facility[R]. AIAA-2009-4090, 2009

[6]Grant P, Dinesh P, Imelda T S. Arc jet facility test condition predictions using the ADSI code[R]. AIAA-2015-2663, 2015

[7]Balboni J A, Dinesh P, Imelda T S. Consolidating NASA’s arc jets [R]. AIAA-2015-2667, 2015

[8]Balboni J A, Dinesh P, Imelda T S. On laminar to turbulent transition of arc-jet flow in the NASA Ames panel test facility[R]. AIAA-2012-3304, 2012

[9]Balboni J A, Dinesh P, Imelda T S. Computational simulation of high enthalpy arc heater flows[R]. AIAA-2006-1183, 2006

[10]Balboni J A, Dinesh P, Imelda T S. CFD analysis framework for arc-heated flowfields, II: shear testing in arc-jets at NASA ARC[R]. AIAA-2009-4081, 2009

[11]Balboni J A, Dinesh P, Imelda T S. Computational simulations of the 10-MW TP3 arc-jet facility flow[R]. AIAA-2015-3103, 2015

[12]張友華, 陳連忠. 超聲速湍流導管燒蝕流場穩定性研究[J]. 宇航材料工藝, 2010, (4): 64-67.Zhang Y H, Chen L Z. Stability of ablative flow field of supersonic turbulent duct[J]. Aerospace Materials & Technology, 2010, (4): 64-67

[13]涂建強, 陳連忠, 許考. 高溫含水氣流條件下燃燒室材料考核的電弧加熱試驗模擬方法[J]. 實驗流體力學, 2015, 29(4): 81-85.Tu J Q, Chen L Z, Xu K. Testing of combustor chamber material in arc jet flow mixing with transverse injected water[J]. Journal of Experiments in Fluid Mechanics, 2015, 29(4): 81-85

[14]李媛. 基于Fluent的傳熱風洞流場品質研究[J]. 機械制造與自動化, 2012, 41(5): 144-149.Li Y. Study of flow quality of heat transfer wind-tunnel base on fluent[J]. Machine Building Automation, 2012, 41(5): 144-149.

[15]于勝春, 湯龍生. 固體火箭發動機噴管及羽流流場的數值分析[J]. 固體火箭技術, 2004, 27(2): 95-97.Yu S C, Tang L S. Numerical analysis of nozzle and plume flow field of a solid rocket motor[J]. Journal of Solid Rocket Technology, 2004, 27(2): 95-97

[16]楊鴻. 超聲速湍流平板燒蝕試驗技術研究[D]. 長沙: 國防科技大學, 2005.Yang H. The research on the supersonic turbulent flat ablationexperiment technique[D]. Changsha: National University of Defense Technology, 2005

[17]楊鴻, 陳偉芳, 柳森. 電弧湍流平板燒蝕矩形噴管研制及應用[J]. 實驗流體力學, 2006, 20(1): 27-30.Yang H, Chen W F, Liu S. The development and application on rectangular nozzle of arc turbulent flat plate ablation test[J]. Journal of Experiments in Fluid Mechanics, 2006, 20(1): 27-30

[18]隆永勝, 胡振震, 袁竭, 等. 橢圓噴管設計與數值模擬[J]. 實驗流體力學, 2015, 29(3): 80-86.Long Y S, Hu Z Z, Yuan J, et al. Design and numerical simulation of an elliptical nozzle[J]. Journal of Experiments in Fluid Mechanics, 2015, 29(3): 80-86

[19]鄒寧. 超聲速噴管設計及其數值模擬和實驗研究[D]. 南京: 南京航天航空大學, 2009.Zou N. Supersonic nozzle design and its investigation with numerical and experimental methods[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2009

[20]王福軍. 計算流體動力學分析[M]. 北京: 清華大學出版社, 2004: 161-172.Wang F J. Computational fluid dynamic analysis[M]. Beijing: Tsinghua University Press, 2004: 161-172.

(編輯:楊 娟)

CFD analysis of the arc heater turbulent flow field of flat plate testing

Luo Yue1,2,*, Zhou Wei1, Yang Hong1, Chen Wei1

(1. China Aerodynamics Research and Development Center, Mianyang Sichuan 621000, China; 2. College of Manufacturing Science and Engineering, Sichuan University, Chengdu 610065, China)

The statement debugging is important but time-consuming work in the arc-heated turbulent flat plate testing. Specially, with the development of the alerting angle method in recent years, the debugging becomes more complex and further consuming in both time and cost. In order to guide the experiments, some CFD analyses of the flow characteristics under different testing conditions are completed with the FLUENT software. By comparing to the experimental measurements of a calibration-plate under some representative conditions, the numerical method is verified firstly. At the same time, the distributions of the primary aerothermal parameters over the plate surface are provided quantitatively under various attack angles and total pressures. Furthermore, by analyzing the shock wave reflection, interaction and boundary layer separation, the rules of the impact of the control parameters to the flow field and some basic principles in experiments are summarized.

arc heater;flat plate test;flow field characteristic;numerical simulation;calibration

2016-06-01;

2016-09-26

LuoY,ZhouW,YangH,etal.CFDanalysisofthearcheaterturbulentflowfieldofflatplatetesting.JournalofExperimentsinFluidMechanics, 2017, 31(2): 86-92. 羅 躍, 周 瑋, 楊 鴻, 等. 電弧加熱器湍流平板試驗流場計算分析. 實驗流體力學, 2017, 31(2): 86-92.

1672-9897(2017)02-0086-06

10.11729/syltlx20160088

V211.74+4

A

羅 躍(1980-),男,四川內江人,工程師。研究方向:高超聲速飛行器熱防護試驗及測試技術研究。通信地址:四川省綿陽市中國空氣動力研究與發展中心(621000)。E-mail:myheater@sohu.com

*通信作者 E-mail: myheater@sohu.com

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 亚洲Av综合日韩精品久久久| 一级毛片免费不卡在线| 精品无码人妻一区二区| 欧美日韩中文字幕在线| 日韩在线成年视频人网站观看| 亚洲精品爱草草视频在线| 99热在线只有精品| av手机版在线播放| 午夜在线不卡| 亚洲免费人成影院| 成人一区在线| 伊人成人在线| 天天摸夜夜操| 久久五月视频| 欧美午夜理伦三级在线观看| 性欧美在线| 国产人成乱码视频免费观看| 国产高清在线丝袜精品一区| 91在线视频福利| 欧美激情成人网| 久久综合伊人77777| 亚洲女同一区二区| 国产熟睡乱子伦视频网站| 亚洲国产在一区二区三区| 成年人国产网站| 国产精品深爱在线| 国产又黄又硬又粗| 污网站在线观看视频| 亚卅精品无码久久毛片乌克兰| 狼友视频国产精品首页| 一区二区午夜| 国禁国产you女视频网站| 亚洲视频影院| 国产综合色在线视频播放线视| 99无码熟妇丰满人妻啪啪 | 久久久亚洲色| 71pao成人国产永久免费视频| 亚洲福利视频网址| 91精品国产无线乱码在线| 免费一级毛片| 亚洲欧美一级一级a| 少妇被粗大的猛烈进出免费视频| 国产高清不卡| 热久久国产| 亚洲精品无码专区在线观看| 国产91小视频在线观看| 亚洲最新地址| 久久婷婷六月| 欧美a级完整在线观看| 亚洲人成人伊人成综合网无码| 亚洲色图综合在线| 亚洲妓女综合网995久久| 日韩毛片免费| 一本大道香蕉久中文在线播放 | 四虎永久免费网站| 久久国产精品影院| 欧美啪啪一区| 国产福利免费视频| 麻豆精品在线视频| 久久6免费视频| 国产va在线观看| 欧美亚洲国产精品第一页| 很黄的网站在线观看| AV无码一区二区三区四区| 国产亚洲一区二区三区在线| 亚洲AⅤ综合在线欧美一区| 国产玖玖玖精品视频| 精品一区二区三区四区五区| 免费中文字幕在在线不卡| 国产91视频免费| 色天天综合| 亚洲黄色成人| 国产区免费| 久久熟女AV| 欧美性猛交xxxx乱大交极品| 欧美性精品不卡在线观看| 波多野吉衣一区二区三区av| 免费激情网站| 欧美精品黑人粗大| 亚洲成在线观看 | 毛片在线播放a| 九色视频线上播放|