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

雷諾數(shù)變化對翼型邊界層發(fā)展及失速特性的影響

2019-07-09 01:45:52張彥軍趙軻張同鑫陳利麗
航空工程進(jìn)展 2019年3期
關(guān)鍵詞:模型

張彥軍,趙軻,張同鑫,陳利麗

(第一飛機(jī)設(shè)計(jì)研究院 總體氣動(dòng)設(shè)計(jì)研究所,西安 710089)

0 引 言

在航空航天領(lǐng)域,失速問題對飛行器設(shè)計(jì)和使用造成惡劣影響,甚至導(dǎo)致機(jī)毀人亡。為了滿足飛行安全和改善起降性能,失速問題一直是氣動(dòng)設(shè)計(jì)師必須面對和解決的基礎(chǔ)科學(xué)問題。探索飛行器的失速問題隱含的空氣動(dòng)力學(xué)機(jī)理是明確失速機(jī)制及建立控制措施的必由之路,但是復(fù)雜的部件干擾和雷諾數(shù)變化增加了研究失速機(jī)理的難度。幸運(yùn)的是,雖然難以給出飛行器失速問題精確的解析關(guān)系式,但組成這種強(qiáng)非線性現(xiàn)象的核心機(jī)制仍然是翼型的失速問題,因此研究翼型的失速問題是建立飛行器失速問題理論體系的基石。

從20世紀(jì)上半葉開始,學(xué)者們采用試驗(yàn)方法對翼型失速的表面氣動(dòng)力和邊界層變化進(jìn)行了諸多分析,E.N.Jacobs[1]第一次推測:升力線的拐點(diǎn)意味著翼型表面出現(xiàn)了分離現(xiàn)象,即提出失速問題;隨后,E.N.Jacobs等[2]第一次提出雷諾數(shù)與翼型分離之間存在關(guān)系,重點(diǎn)關(guān)注層流轉(zhuǎn)捩現(xiàn)象對翼型前緣和后緣分離的影響;Lissaman[3]闡述低速翼型與雷諾數(shù)之間的影響作用,并分析雷諾數(shù)影響下分離泡的延伸與破裂是導(dǎo)致翼型失速的主要原因;V.O.Michael等[4]基于PIV技術(shù)研究低雷諾數(shù)翼型分離泡的相關(guān)特性;Hu H等[5]使用PIV檢測GA(W)-1翼型流場和表面壓力分布,分析層流分離泡與壓力梯度的相互影響;H.Yamato等[6]使用動(dòng)態(tài)控制裝置影響前緣分離泡破碎,從而抑制翼型分離。國內(nèi), Wang N等[7]通過試驗(yàn)研究雷諾數(shù)對RAE2822翼型氣動(dòng)特性的影響;孫智偉等[8]在某NF-3和NF-6風(fēng)洞中進(jìn)行超臨界翼型的設(shè)計(jì)和實(shí)驗(yàn)研究,給出超臨界翼型在高、低速時(shí)的分離演化特性。

在數(shù)值模擬方面,A.C.Aranake等[9]使用SA湍流模型耦合γ-Reθt轉(zhuǎn)捩模型對S827和S809翼型分別進(jìn)行失速特性數(shù)值模擬,相比于全湍流計(jì)算結(jié)果,無論是升力特性還是阻力特性均與試驗(yàn)結(jié)果更加接近,且考慮邊界層轉(zhuǎn)捩時(shí)計(jì)算所得氣動(dòng)力系數(shù)能夠反應(yīng)失速狀態(tài)下的輕失速和深失速特性;Ke J等[10]采用耦合轉(zhuǎn)捩模型的RANS方法和混合方法對A-Airfoil及NACA0012翼型靜態(tài)失速特性進(jìn)行模擬,分析前緣層流分離泡及后緣分離發(fā)展,其結(jié)果與試驗(yàn)結(jié)果吻合良好。國內(nèi),馮濤等[11]針對典型氣動(dòng)負(fù)載渦輪葉片在不同雷諾數(shù)下邊界層特性變化進(jìn)行敏感性研究;靳允立等[12]采用不同湍流模型研究翼型失速及風(fēng)速變化引起的雷諾數(shù)改變對翼型和風(fēng)力機(jī)氣動(dòng)性能的影響;王強(qiáng)等[13]采用基于SST模型的γ-Reθt轉(zhuǎn)捩模型研究S809翼型動(dòng)態(tài)失速過程中分離泡的作用;袁尚科[14]通過理論分析、數(shù)值計(jì)算及外場實(shí)驗(yàn)相結(jié)合的方法,分別對風(fēng)力機(jī)常用翼型及葉片的靜態(tài)失速、動(dòng)態(tài)失速及三維旋轉(zhuǎn)效應(yīng)等進(jìn)行研究。

綜上所述,國內(nèi)外諸多學(xué)者將翼型失速的研究對象重點(diǎn)集中于失速點(diǎn)的精確預(yù)測,對于失速時(shí)邊界層的變化,層流轉(zhuǎn)捩現(xiàn)象的影響作用和雷諾數(shù)對邊界層演化的完整影響過程并未有過多涉足。

本文使用Menter SST湍流模型和擾動(dòng)放大因子輸運(yùn)模型對轉(zhuǎn)捩翼型邊界層進(jìn)行數(shù)值模擬,分析雷諾數(shù)變化時(shí)邊界層的發(fā)展,總結(jié)雷諾數(shù)變化對翼型失速的影響規(guī)律。

1 數(shù)值模擬方法

流場模擬采用Navier-Stokes方程組進(jìn)行計(jì)算,其守恒積分形式為

(1)

式中:Q為守恒向量;S為控制體V所有表面S的集合;F為通過控制體的S面上的通量;n為S面的單位外法向矢量。無粘項(xiàng)、粘性項(xiàng)的空間離散方法分別為二階迎風(fēng)Roe格式及二階中心差分方法。時(shí)間離散則采用近似因子分解方法。

Menterk-ωSST湍流模型[15]是一個(gè)兩方程k-ω渦粘性模型,在工程界廣泛應(yīng)用。SST模型是一個(gè)耦合模型,在近壁面保持Wilcoxk-ω模型[16]的特性,而在湍流邊界層外及自由來流中保持Launder-Sharmak-ε模型[17]的特性。通過混合函數(shù)保留兩個(gè)模型的優(yōu)點(diǎn)并克服模型各自的不足。其方程形式為

(2)

(3)

渦粘性由下式計(jì)算得到

(4)

公式中模型參數(shù)(α,β,σ,F(xiàn)1,F(xiàn)2)的數(shù)值及處理方法可以參照文獻(xiàn)[15]。

在模擬層流-湍流轉(zhuǎn)捩邊界層時(shí),采用Coder等[18]提出的擾動(dòng)放大因子輸運(yùn)模型。該輸運(yùn)方程模型形式如式(5)所示。

(5)

式中的Fcrit,F(xiàn)growth,dn/dRδ2均與邊界層形狀因子相關(guān),關(guān)系式的具體表示可參考文獻(xiàn)[18]。因此模型的核心點(diǎn)為使用層流邊界層的當(dāng)?shù)貕毫μ荻葏?shù)HL計(jì)算非當(dāng)?shù)貐?shù)形狀因子H12,后者常用于描述速度型的特征。當(dāng)?shù)貕毫μ荻葏?shù)定義為

(6)

式中:d為壁面法向距離;S為當(dāng)?shù)丶羟袘?yīng)變率;Ue為邊界層邊緣速度。

由Falkner-Skan二維邊界層相似性解可以得到不同壓力梯度下HL最大值與形狀因子的函數(shù)關(guān)系,即

(7)

從而實(shí)現(xiàn)在流場中形狀因子H12的求解,邊界層中max(H12)為精確值。

將擾動(dòng)放大因子輸運(yùn)方程與SST模型耦合,其湍動(dòng)能輸運(yùn)方程為

(8)

(9)

ft2=[1-exp(2(n-Ncrit))]exp(ct4χ2)

(10)

式中:模型參數(shù)F1,ct4,ct5,χ,Ncrit等的設(shè)置與文獻(xiàn)[19]一致。

2 翼型轉(zhuǎn)捩邊界層發(fā)展及失速影響分析

2.1 翼型和網(wǎng)格

使用NLF(1)-0416自然層流翼型繞流[20]研究雷諾數(shù)的影響。該翼型有較長的層流區(qū)且具備高升力特性,為典型后緣失速翼型,一般應(yīng)用于通用航空飛機(jī)翼型和風(fēng)力機(jī)葉型等領(lǐng)域。對其進(jìn)行轉(zhuǎn)捩數(shù)值模擬,計(jì)算網(wǎng)格分布如圖 1所示,計(jì)算域半徑為100倍弦長,翼型表面布置481個(gè)網(wǎng)格單元,法向布置181個(gè)網(wǎng)格單元,網(wǎng)格量約為18萬。壁面法向第一層網(wǎng)格高度約為1.0×10-6m,使y+小于1。風(fēng)洞試驗(yàn)雷諾數(shù)工況分別為2×106,4×106,9×106。風(fēng)洞自由來流馬赫數(shù)為0.1,湍流度設(shè)置為0.03%。

圖1 翼型繞流數(shù)值模擬網(wǎng)格分布

2.2 數(shù)值模擬結(jié)果

為了與試驗(yàn)狀態(tài)保持一致,便于對比驗(yàn)證和分析,本文選擇了與試驗(yàn)雷諾數(shù)相同的三個(gè)雷諾數(shù)進(jìn)行數(shù)值模擬研究,鑒于雷諾數(shù)對失速特性的影響一直是一個(gè)復(fù)雜的問題,對于更大雷諾數(shù),我們將在后續(xù)工作中開展研究。使用擾動(dòng)放大因子輸運(yùn)模型對三個(gè)雷諾數(shù)下翼型繞流進(jìn)行模擬,得到升力系數(shù)-阻力系數(shù)極曲線并與全湍數(shù)據(jù)和試驗(yàn)數(shù)據(jù)進(jìn)行對比,結(jié)果如圖2所示。

(a) Re=2×106

(b) Re=4×106

(c) Re=9×106

從圖2可以看出:相比全湍模型的數(shù)值模擬結(jié)果,轉(zhuǎn)捩模型的預(yù)測結(jié)果與試驗(yàn)數(shù)據(jù)更加貼合,與真實(shí)流動(dòng)更加吻合。

雷諾數(shù)為4×106,0°迎角工況的翼型表面壓力系數(shù)分布如圖3所示。

圖3 Re=4.0×106,0°迎角翼型表面的壓力系數(shù)分布

從圖3可以看出:在該工況下,翼型繞流轉(zhuǎn)捩類型為層流分離泡轉(zhuǎn)捩,且轉(zhuǎn)捩模式對該分離再附著過程描述較好。模型對轉(zhuǎn)捩過程的準(zhǔn)確預(yù)測為分析邊界層發(fā)展奠定基礎(chǔ)。

2.3 雷諾數(shù)影響分析

三個(gè)雷諾數(shù)下翼型升力系數(shù)隨迎角變化曲線如圖4所示,可以看出:雷諾數(shù)對翼型的氣動(dòng)特性有較大的影響,當(dāng)雷諾數(shù)增大時(shí),翼型的失速迎角變大,同時(shí)最大升力系數(shù)增大。產(chǎn)生此現(xiàn)象的主要原因是雷諾數(shù)變化對翼型邊界層轉(zhuǎn)捩及發(fā)展有重要影響,因此通過對比不同雷諾數(shù)下的轉(zhuǎn)捩過程及邊界層發(fā)展情況,對雷諾數(shù)的影響規(guī)律進(jìn)行分析。

圖4 不同雷諾數(shù)下翼型升力系數(shù)隨迎角變化曲線

不同雷諾數(shù)下翼型上下表面的轉(zhuǎn)捩位置隨升力系數(shù)變化曲線如圖5所示。

圖5 不同雷諾數(shù)翼型轉(zhuǎn)捩位置隨升力系數(shù)的變化

從圖5可以看出:隨著迎角增大,升力系數(shù)增大,翼型上表面的轉(zhuǎn)捩位置逐漸前移,而下表面變化趨勢相反。相同升力系數(shù)下,當(dāng)雷諾數(shù)增大時(shí),上、下表面的轉(zhuǎn)捩位置均逐漸前移。

對0°迎角工況下翼型邊界層參數(shù)隨雷諾數(shù)的變化進(jìn)行分析。0°迎角時(shí)翼型的表面壓力分布圖如圖6所示,可以看出:雷諾數(shù)變化對壓力系數(shù)的整體影響較小,其主要影響體現(xiàn)在壓力系數(shù)突變位置的變化,即出現(xiàn)層流分離泡、發(fā)生轉(zhuǎn)捩的位置。

圖6 0°迎角時(shí)翼型表面壓力系數(shù)分布對比

不同雷諾數(shù)在0°迎角時(shí)翼型的表面摩擦力系數(shù)分布如圖7所示,可以看出:在雷諾數(shù)為2×106和4×106的工況下,翼型上下表面的摩擦力系數(shù)均出現(xiàn)負(fù)值,表面附面層中出現(xiàn)層流分離泡,并伴隨轉(zhuǎn)捩的發(fā)生;而雷諾數(shù)為9×106時(shí),翼型繞流發(fā)生層流到湍流的轉(zhuǎn)捩過程,但是并沒有層流分離泡出現(xiàn)。初步的結(jié)論與傳統(tǒng)的研究結(jié)論一致,低雷諾數(shù)環(huán)境雖然轉(zhuǎn)捩位置比較靠后,但是層流非常容易被逆壓梯度誘導(dǎo)產(chǎn)生層流分離泡,觸發(fā)轉(zhuǎn)捩。該現(xiàn)象在高空長航時(shí)無人機(jī)表面比較常見,因?yàn)閿?shù)萬米高空中空氣相對稀薄,雷諾數(shù)較小,容易產(chǎn)生層流分離泡。

(a) Re=2×106

(b) Re=4×106

(c) Re=9×106

翼型上、下表面層流分離泡形態(tài)如圖8所示,可以看出:當(dāng)雷諾數(shù)增大時(shí),翼型上下表面的層流分離泡均逐漸減小,在雷諾數(shù)為9×106時(shí),層流分離泡已經(jīng)消失。

(a1) 上表面分離泡

(a2) 下表面分離泡

(b1) 上表面分離泡

(b2) 下表面分離泡

(c1) 上表面分離泡

(c2) 下表面分離泡

對比不同雷諾數(shù)下翼型上表面邊界層的邊界層動(dòng)量損失厚度δ2、形狀因子H12、渦量雷諾數(shù)Reν等參數(shù)分布,如圖9~圖11所示。

圖9 0°迎角上表面邊界層動(dòng)量損失厚度分布

圖10 0°迎角翼型上表面邊界層形狀因子分布

圖11 翼型上表面邊界層渦量雷諾數(shù)分布

從圖9可以看出:當(dāng)雷諾數(shù)減小時(shí),動(dòng)量損失厚度增大,這是由于雷諾數(shù)減小使流動(dòng)粘性作用相對增強(qiáng),有利于邊界層的發(fā)展。

邊界層形狀因子的定義為

(11)

式中:δ1為邊界層位移厚度。

從圖10可以看出:當(dāng)層流邊界層通過轉(zhuǎn)捩點(diǎn)變?yōu)橥牧鬟吔鐚樱吔鐚雍穸葧?huì)突然增厚。湍流邊界層由于流速分布更趨均勻化而使邊界層位移厚度減小,同時(shí)由于阻力增加而使邊界層動(dòng)量損失厚度增大,于是邊界層形狀因子在轉(zhuǎn)捩后將減小。翼型上表面在三個(gè)雷諾數(shù)下的轉(zhuǎn)捩位置均在40%弦長附近,因此上表面邊界層形狀因子在40%弦長處出現(xiàn)峰值,而且由于雷諾數(shù)為2×106和4×106時(shí)翼型出現(xiàn)層流分離泡,所以這兩個(gè)雷諾數(shù)對應(yīng)的形狀因子峰值更高;而雷諾數(shù)為9×106時(shí)邊界層未出現(xiàn)分離泡,且轉(zhuǎn)捩位置前移,40%弦長處流動(dòng)已基本發(fā)展為湍流,對應(yīng)的形狀因子值較低。在轉(zhuǎn)捩發(fā)生后,形狀因子的發(fā)展情況與前文的分析對應(yīng):三組形狀因子均迅速下降,雷諾數(shù)越大,最終邊界層形狀因子的值越小。

從圖11可以看出:相同弦長站位處,雷諾數(shù)越大,邊界層中的渦量雷諾數(shù)越高。因?yàn)闇u量雷諾數(shù)與常用轉(zhuǎn)捩判定參數(shù)動(dòng)量厚度雷諾數(shù)幾乎是正比關(guān)系,因此可以作為邊界層發(fā)展的標(biāo)志性參數(shù)進(jìn)行分析。邊界層轉(zhuǎn)捩觸發(fā)位置隨雷諾數(shù)升高而向前緣移動(dòng),在轉(zhuǎn)捩后,渦量雷諾數(shù)迅速增大,使邊界層保持在湍流狀態(tài)。

0°迎角時(shí)翼型上表面不同站位流向速度型分布對比如圖12所示。

(a) 0.2c

(b) 0.3c

(c) 0.4c

(d) 0.5c

(e) 0.6c

(f) 0.7c

(g) 0.8c

(h) 0.9c

從圖12可以看出:在翼型前緣,轉(zhuǎn)捩前較高雷諾數(shù)對應(yīng)邊界層的速度型更加均勻、飽滿,邊界層更薄。在40%弦長處,由于翼型邊界層在較低雷諾數(shù)時(shí)出現(xiàn)層流分離泡,速度型分布在壁面附近呈凹曲線,速度型出現(xiàn)拐點(diǎn)。轉(zhuǎn)捩后,三個(gè)雷諾數(shù)對應(yīng)的速度型更加相近,高雷諾數(shù)的速度型分布更加飽滿。

為了分析雷諾數(shù)對失速迎角附近邊界層發(fā)展的影響,取12°迎角工況下翼型邊界層進(jìn)行觀察。12°迎角時(shí)翼型的表面摩擦力系數(shù)分布如圖13所示,隨雷諾數(shù)增高翼型轉(zhuǎn)捩點(diǎn)前移,轉(zhuǎn)捩處層流分離泡逐漸縮小消失,僅雷諾數(shù)為2×106時(shí)存在比較明顯的層流分離泡。

圖13 12°迎角翼型表面摩擦力系數(shù)分布

12°迎角時(shí)翼型上表面不同站位流向速度型分布對比如圖14所示。由于雷諾數(shù)為2×106時(shí)在約7%弦長處產(chǎn)生層流分離泡,其對應(yīng)的速度分布出現(xiàn)明顯的拐折。雷諾數(shù)為9×106時(shí)的速度型更早趨于均勻,因其轉(zhuǎn)捩點(diǎn)更加靠近前緣。轉(zhuǎn)捩后,隨雷諾數(shù)的增高,流速分布變得更加均勻而壁面附近流速梯度增大。

翼型上表面邊界層內(nèi)各站位處的湍動(dòng)能分布如圖15所示。從翼型表面摩擦力系數(shù)分布可知,三個(gè)雷諾數(shù)下翼型上表面均已在10%弦長前完成轉(zhuǎn)捩,因此下圖中的結(jié)果均為轉(zhuǎn)捩后邊界層湍動(dòng)能的分布。

(a) 0.005c(b) 0.010c(c) 0.030c

(d) 0.060c(e) 0.100c(f) 0.200c

(g) 0.400c(h) 0.600c(i) 0.900c

圖14 12°迎角翼型上表面不同站位速度型分布

Fig.14 The velocity distribution of different stations on the upper surface of airfoil at 12° angle of attack

(a) 0.2c

(b) 0.4c

(c) 0.6c

(d) 0.9c

從圖15可以看出:轉(zhuǎn)捩后相同站位低雷諾數(shù)對應(yīng)的湍動(dòng)能峰值更高。這是由于翼型前緣發(fā)生轉(zhuǎn)捩,在較低雷諾數(shù)時(shí)伴隨有層流分離泡形成,這一過程使邊界層流動(dòng)的脈動(dòng)特性得到加強(qiáng),其對應(yīng)的湍動(dòng)能、渦粘性增大。但是層流分離泡會(huì)消耗邊界層的動(dòng)能,使流動(dòng)自持能力即抗分離能力減弱,從而產(chǎn)生分離。當(dāng)雷諾數(shù)增高,分離泡減小甚至消失,流動(dòng)能量消耗減小,動(dòng)能更充沛,壓力分布可以維持較長距離的梯度,后緣分離滯后發(fā)生,從而增大翼型升力系數(shù)及失速迎角。

最后對比相同迎角不同雷諾數(shù)下翼型流場的后緣分離情況,結(jié)果如圖16~圖17所示。

(a) Re=2×106

(b) Re=4×106

(c) Re=9×106

(a) Re=2×106

(b) Re=4×106

(c) Re=9×106

從圖16~圖17可以看出:在相同迎角下,隨著雷諾數(shù)增大,流動(dòng)穩(wěn)定性增強(qiáng),抵抗逆壓梯度干擾的能力提高,流動(dòng)分離趨勢減弱。后緣分離逐漸減小,這與本節(jié)中對邊界層發(fā)展及其參數(shù)變化趨勢的觀察結(jié)果相一致。

3 結(jié) 論

(1) 對于層流-湍流轉(zhuǎn)捩翼型邊界層,雷諾數(shù)增大,使流動(dòng)粘性作用相對減弱,會(huì)增大邊界層當(dāng)?shù)販u量雷諾數(shù),從而使翼型轉(zhuǎn)捩位置前移,減小層流分離泡尺寸。

(2) 分離泡減小使邊界層能量損耗減小,同時(shí)轉(zhuǎn)捩位置前移使湍流區(qū)增長,這些使得流動(dòng)動(dòng)能更加充沛,流動(dòng)穩(wěn)定性得到加強(qiáng),抵抗逆壓梯度干擾的能力增強(qiáng),流動(dòng)分離趨勢減弱。

(3) 當(dāng)雷諾數(shù)增大時(shí),翼型后緣分離將延遲,失速迎角及最大升力系數(shù)提高。

猜你喜歡
模型
一半模型
一種去中心化的域名服務(wù)本地化模型
適用于BDS-3 PPP的隨機(jī)模型
提煉模型 突破難點(diǎn)
函數(shù)模型及應(yīng)用
p150Glued在帕金森病模型中的表達(dá)及分布
函數(shù)模型及應(yīng)用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 欧美无遮挡国产欧美另类| 在线中文字幕网| 熟女视频91| 青青青国产在线播放| 亚洲Va中文字幕久久一区| 国产欧美在线视频免费| 日韩欧美在线观看| 97免费在线观看视频| 精品伊人久久久香线蕉 | 国产永久免费视频m3u8| 亚洲成A人V欧美综合| 欧美日韩北条麻妃一区二区| 亚洲精品制服丝袜二区| 国产在线观看高清不卡| 激情午夜婷婷| 99国产精品免费观看视频| 欧美一区日韩一区中文字幕页| a级高清毛片| 波多野结衣无码中文字幕在线观看一区二区| 国产JIZzJIzz视频全部免费| 好吊色妇女免费视频免费| 中文国产成人久久精品小说| 99久久精品国产精品亚洲 | 综合网久久| 怡红院美国分院一区二区| 天堂亚洲网| 免费看av在线网站网址| 久久九九热视频| 欧美成人手机在线视频| 欧美精品成人| 成人小视频网| 国产三级韩国三级理| 亚洲欧美色中文字幕| 欧美日韩精品一区二区在线线| 久久久久中文字幕精品视频| 欧美在线视频不卡第一页| 国产永久无码观看在线| 亚洲另类色| 亚洲有无码中文网| 日韩东京热无码人妻| 久久久亚洲色| 国产精品久久自在自线观看| 国内精品视频区在线2021| 中文字幕欧美日韩高清| 一区二区自拍| 亚洲Aⅴ无码专区在线观看q| 5555国产在线观看| 亚洲成人www| 99精品视频九九精品| 大乳丰满人妻中文字幕日本| 久久人妻xunleige无码| 亚洲大尺度在线| 中文字幕永久在线看| 激情视频综合网| 毛片最新网址| 欧美日韩中文国产| 黄色网站不卡无码| 国产精品第三页在线看| 91原创视频在线| 国产不卡在线看| 国产成人精品在线1区| 美女国内精品自产拍在线播放 | 国产日本欧美在线观看| 久一在线视频| 国产日产欧美精品| 伊人久久大香线蕉影院| 亚洲国产91人成在线| 992tv国产人成在线观看| 国产三级毛片| 亚洲成人在线网| 亚洲狠狠婷婷综合久久久久| 国产在线八区| 亚洲成A人V欧美综合| 国产精品妖精视频| 大乳丰满人妻中文字幕日本| 午夜福利在线观看入口| 免费毛片网站在线观看| 亚洲欧美一区二区三区麻豆| 日韩欧美国产中文| 在线观看精品国产入口| 欧美色综合网站| 国产毛片不卡|