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

基于聯(lián)合射流的翼型動態(tài)失速流動控制研究

2018-11-28 11:49:32楊慧強許和勇葉正寅
航空工程進(jìn)展 2018年4期

楊慧強,許和勇,葉正寅

(1.西北工業(yè)大學(xué) 航空學(xué)院,西安 710072)(2.西北工業(yè)大學(xué) 深圳研究院,深圳 518057)

0 引 言

動態(tài)失速現(xiàn)象廣泛存在于直升機旋翼、風(fēng)力機葉片、機動飛行器、低雷諾數(shù)撲翼和鳥類昆蟲的繞流中。顧名思義,動態(tài)失速是發(fā)生在攻角隨時間動態(tài)變化的升力部件上的失速現(xiàn)象。由于動態(tài)效應(yīng),動態(tài)失速的失速攻角一般大于相應(yīng)的靜態(tài)失速攻角,具有增大最大升力的作用。但是,隨之而來的是升力部件吸力面動態(tài)失速渦的產(chǎn)生、移動和脫落,形成非常復(fù)雜的流動現(xiàn)象,產(chǎn)生很大的阻力和俯仰力矩峰值,給結(jié)構(gòu)帶來嚴(yán)重的振動載荷和穩(wěn)定性問題。對于飛行器來說,動態(tài)失速會極大地限制機翼或槳葉的氣動性能及飛行包線。因此,采用流動控制手段對動態(tài)失速進(jìn)行控制一直是空氣動力學(xué)領(lǐng)域的研究熱點之一。

按照是否需要外部能量輸入來分類,動態(tài)失速控制分為被動控制和主動控制。被動控制方法包括渦流發(fā)生器[1-2]、格尼襟翼[3-4]、固定前緣下垂[3]、固定前緣縫翼[5]、仿生波狀前緣[6]等。其中,組合使用格尼襟翼和固定前緣下垂方法可以起到較佳的控制效果,既可以減小遲滯環(huán)面積,又可以增加升力并減小力矩系數(shù)負(fù)峰值。主動控制方法由于其有效性和靈活性一直備受關(guān)注,形成了各種各樣的主動控制策略,且均可以達(dá)到不同程度的控制效果,例如后緣偏轉(zhuǎn)舵面[7-11]、動態(tài)前緣下垂[12-16]、前緣或表面變形[17-21]、射流型渦流發(fā)生器[22]、表面吹氣[23-24]、合成射流[25-29]、等離子體激勵[30-32]等。G.C.Zha等[33]于2004年提出了一種新穎的零質(zhì)量組合吹吸氣控制方法——Co-Flow Jet(簡稱CFJ)。該方法在翼型前緣設(shè)計噴氣口進(jìn)行吹氣,同時在后緣設(shè)計吸氣口進(jìn)行等量吸氣,從而達(dá)到零質(zhì)量射流的效果。翼型內(nèi)部的氣泵裝置連接前緣噴氣口和后緣吸氣口,實現(xiàn)氣源的循環(huán)供給,因而無需從發(fā)動機引氣。不同于傳統(tǒng)的基于薄膜振動的單孔合成射流方法,該方法形成的射流屬于連續(xù)射流,可以靈活控制射流強度,能夠向外流場中注入更多的能量。此外,前緣噴氣口處于吸力峰附近,壓強比后緣吸氣口低,氣泵所消耗的能量也更小。靜態(tài)翼型的風(fēng)洞實驗[34]和數(shù)值模擬[35]的結(jié)果表明,CFJ可以大幅增升減阻并增大失速攻角,使得CFJ翼型具有非常優(yōu)異的氣動性能。A.Lefebvre等[36]通過數(shù)值模擬初步研究了CFJ對旋翼翼型動態(tài)失速的抑制效果,表明CFJ對不同程度的動態(tài)失速均有顯著的抑制效果,但未進(jìn)行更詳細(xì)的參數(shù)影響研究和流動分析。國內(nèi),劉沛清等[37]最早對CFJ(論文稱CFJ為聯(lián)合射流)的靜態(tài)增升效果進(jìn)行了跟蹤研究;隨后,朱敏等[38]將CFJ(論文稱CFJ為協(xié)同射流)應(yīng)用至螺旋槳的增效研究中;Xu H Y等[39]將CFJ應(yīng)用至風(fēng)力機翼型的流動控制中。本文沿用劉沛清等[37]最早對CFJ的譯法,下文統(tǒng)一稱CFJ為聯(lián)合射流。鑒于聯(lián)合射流方法的新穎性和其在流動控制中的應(yīng)用前景,有必要進(jìn)一步針對其動態(tài)失速控制開展更為深入和詳細(xì)的研究。

本文采用計算流體力學(xué)方法,對翼型動態(tài)失速的聯(lián)合射流控制進(jìn)行研究,主要對比不同射流動量系數(shù)下的動態(tài)失速控制效果,并分析在射流關(guān)閉情況下聯(lián)合射流在上翼面的氣流通道對動態(tài)失速特性的影響。

1 計算方法

1.1 計算模型和網(wǎng)格

本文以NACA0012翼型為基準(zhǔn)翼型,通過對基準(zhǔn)翼型的修形得到相應(yīng)的聯(lián)合射流翼型NACA0012-CFJ,二者輪廓線的對比如圖1所示。CFJ翼型的前緣噴氣口設(shè)置在距離前緣7%弦長處,高度為0.45%弦長,后緣吸氣口設(shè)置在距離前緣85%弦長處,高度為0.90%弦長。吸氣口高度為噴氣口高度的2倍,是為了避免等質(zhì)量吸氣時在吸氣口發(fā)生壅塞。

圖1 NACA0012與NACA0012-CFJ翼型對比

為了使數(shù)值模擬的噴氣口和吸氣口的氣流更加貼近實際情況,計算中保留高壓氣室和低壓氣室。高壓氣室的右端邊界為射流的入流邊界,低壓氣室的左端邊界為射流的出流邊界。此外,為了盡可能保持射流的附壁效應(yīng),在上翼面設(shè)計射流通道,通過將原翼型表面進(jìn)行向下平移和微幅旋轉(zhuǎn)得到。

計算網(wǎng)格采用結(jié)構(gòu)化網(wǎng)格,如圖2所示,在NACA0012基準(zhǔn)翼型輪廓線之外采用常規(guī)的O網(wǎng)格拓?fù)洌诟邏簹馐摇⑸淞魍ǖ馈⒌蛪簹馐业膮^(qū)域內(nèi)采用H形網(wǎng)格。為了使噴氣口和吸氣口附近表面網(wǎng)格點過渡更加光滑,進(jìn)行局部加密。

圖2 NACA0012-CFJ翼型的計算網(wǎng)格圖

本文主要研究翼型俯仰運動所引起的動態(tài)失速控制。翼型俯仰運動規(guī)律可描述為

α(t)=α0+αm·sin(ωt)

(1)

式中:α0為平均攻角;αm為攻角振幅;ω為振動角頻率。

1.2 數(shù)值模擬方法

流場模擬采用的控制方程為非定常雷諾平均Navier-Stokes方程,湍流模型為一方程Spalart-Allmaras模型[40]。方程采用中心有限體積法離散[41],其中無黏通量項采用Roe格式[42]離散和三階MUSCL格式插值[43],黏性通量項采用二階中心差分格式離散。物理時間推進(jìn)采用雙時間方法[44],偽時間推進(jìn)采用隱式LU-SGS迭代[45]。為了進(jìn)一步加快計算收斂速度,采用偽時間步的當(dāng)?shù)貢r間步長方法以及OpenMP并行計算技術(shù)[46]。壁面為無滑移邊界條件,遠(yuǎn)場邊界條件為基于黎曼不變量分析的無反射條件[47]。高壓氣室的進(jìn)氣邊界采用入流邊界條件,總壓、總溫、入流角給定,入流速度由內(nèi)場插值得到;低壓氣室的出氣邊界采用出流邊界條件,靜壓給定,密度、水平速度和垂直速度均由內(nèi)場插值得到。對于聯(lián)合射流方法,通常采用射流動量系數(shù)來反映射流強度的大小,其定義為

(2)

在非定常計算的偽時間迭代中,每隔一定步數(shù)修正入流邊界的總壓值和出流邊界的靜壓值,以達(dá)到設(shè)定的射流動量系數(shù)值并保證入流質(zhì)量流與出流質(zhì)量流相等。

2 結(jié)果與討論

2.1 算例驗證

為了驗證本文計算方法的可靠性,首先對實驗?zāi)P瓦M(jìn)行計算對比,并進(jìn)行網(wǎng)格無關(guān)性和時間步長無關(guān)性驗證。實驗?zāi)P秃蜖顟B(tài)取自文獻(xiàn)[48],Ma=0.283,Re=3.45×106,α0=15°,αm= 10°,κ=0.151。計算域取半徑為50倍弦長的圓形區(qū)域,生成三套單元數(shù)量不同的網(wǎng)格:稀網(wǎng)格200×100,中等網(wǎng)格300×150,密網(wǎng)格400×200。第一層網(wǎng)格高度均為1×10-5c,附面層內(nèi)網(wǎng)格的增長比分別為1.20,1.15,1.10。無量綱時間步長為dt=0.01(以c/V∞無量綱化)。采用此三套網(wǎng)格計算得到的結(jié)果與實驗值[48]的對比情況如圖3所示,可以看出:在翼型上仰過程的附著流階段,升阻力和力矩的計算值與實驗值均吻合得很好,但是在失速發(fā)生之后,計算值和實驗值的吻合度有所下降,特別是接近最大攻角時的峰值存在顯著差別;動態(tài)失速發(fā)生以后,流場以大范圍分離渦結(jié)構(gòu)為主導(dǎo),本文采用的RANS方法對于分離流動模擬存在固有的缺陷,可能是出現(xiàn)上述差異的主要原因;但是,從總體上看,計算值和實驗值吻合較好。通過對比三套網(wǎng)格的計算值可以看出,三套網(wǎng)格的結(jié)果差異不大,特別是中等網(wǎng)格和密網(wǎng)格之間差別很小,可以認(rèn)為中等網(wǎng)格已經(jīng)基本滿足了網(wǎng)格無關(guān)性的要求。此外,通過對比dt分別為0.005、0.010、0.015三個物理時間步長的計算結(jié)果,驗證了dt=0.010基本滿足物理時間步長的無關(guān)性要求。因此,下文統(tǒng)一采用中等網(wǎng)格規(guī)模進(jìn)行計算研究,相應(yīng)的NACA0012-CFJ翼型網(wǎng)格在此NACA0012中等網(wǎng)格基礎(chǔ)上通過補充射流通道和氣室內(nèi)部的網(wǎng)格而得到,物理時間步長均采用dt=0.010。

(a) 升力系數(shù)

(b) 阻力系數(shù)

(c) 力矩系數(shù)

2.2 射流通道對動態(tài)失速特性的影響

對于一種主動流動控制方法,其對原始外形的改動所造成的氣動特性影響是評估該方法可靠性的一個重要方面。本文的聯(lián)合射流方法在實施過程中對上翼面的外形做了一定改動(如圖1所示),相比于原始翼型,其上翼面流動區(qū)域多出了射流通道部分。在聯(lián)合射流系統(tǒng)出現(xiàn)故障或者因其他原因停止工作時,原射流通道將會對繞流產(chǎn)生影響。射流關(guān)閉時,噴氣口和吸氣口均無射流通過,因此可以認(rèn)為噴氣口和吸氣口均為固壁邊界,計算中不考慮兩個氣室。翼型的俯仰運動狀態(tài)與2.1節(jié)相同。基準(zhǔn)NACA0012翼型和射流關(guān)閉后的NACA0012-CFJ_jet-off翼型的計算結(jié)果對比如圖4所示,可以看出:射流通道的影響主要體現(xiàn)在翼型上仰時的附著流階段,升力系數(shù)和力矩系數(shù)的絕對值明顯下降,阻力系數(shù)略有增加;在動態(tài)失速發(fā)生之后的分離流階段,射流通道對氣動特性的影響很小。

(a) 升力系數(shù)

(b) 阻力系數(shù)

(c) 力矩系數(shù)

上仰過程中瞬時攻角α=16.71°時基準(zhǔn)翼型和射流關(guān)閉翼型的流場對比如圖5所示。

(a) 基準(zhǔn)翼型

(b) 射流關(guān)閉翼型

從上翼面前緣附近流場的放大圖可以看出:基準(zhǔn)翼型前緣附近流動附著壁面,流線光滑,而CFJ翼型在噴氣口處的臺階使得流動產(chǎn)生后臺階流動,顯著降低了臺階后上翼面的平均流速,因此CFJ翼型上翼面的低壓水平下降,導(dǎo)致升力系數(shù)下降。這解釋了圖4(a)中CFJ翼型上仰階段升力曲線有明顯下移量的原因。另外,該臺階的存在也會使翼型的分離時刻有所推遲。

下俯過程中α=24.14°時基準(zhǔn)翼型和射流關(guān)閉翼型的流場對比如圖6所示,可以看出:二者流場差別不大,射流通道的存在并沒有顯著影響流動結(jié)構(gòu)。實際上,從空氣動力學(xué)原理的角度來看,固壁形狀的微小差異能夠顯著影響附著流狀態(tài)的氣動性能,但是對大范圍分離流動的影響很小。對于發(fā)生了動態(tài)失速的下俯階段,翼型背風(fēng)面流動從前緣開始即已處于大范圍分離狀態(tài),而射流通道完全處于分離區(qū),故二者氣動特性差異很小。

(a) 基準(zhǔn)翼型

(b) 射流關(guān)閉翼型

2.3 射流動量系數(shù)對動態(tài)失速控制的影響

基準(zhǔn)NACA0012翼型和射流打開后的NACA0012-CFJ翼型的計算結(jié)果對比如圖7所示,可以看出:聯(lián)合射流極大地改善了翼型的動態(tài)失速特性,翼型的失速程度和失速后氣流再附著的時間大幅降低,升力系數(shù)極大提升,阻力系數(shù)也顯著下降,在低攻角階段由于較大的射流反推力甚至出現(xiàn)負(fù)阻力,同時阻力系數(shù)和力矩系數(shù)的峰值顯著降低,且射流動量系數(shù)越大,上述改善效果越明顯。

(a) 升力系數(shù)

(b) 阻力系數(shù)

(c) 力矩系數(shù)

在此計算狀態(tài)下,翼型從前緣開始出現(xiàn)分離,由于聯(lián)合射流從噴氣口高速噴出的氣流加速了翼型上表面的流動,極大地提高了吸力峰處的流速,而翼型前緣過高的氣流速度在繞過相同曲率的前緣時貼附翼型表面顯得更加困難,因此聯(lián)合射流使此狀態(tài)下的分離提前。且Cμ越大,提前得越多,各翼型上表面前緣出現(xiàn)分離的時刻以及此時前緣吸力峰處最大馬赫數(shù)如表1所示,可以看出:在α=16.37°時,原始翼型的前緣最大馬赫數(shù)為1.22,此時Cμ=0.10的CFJ翼型的前緣最大馬赫數(shù)高達(dá)1.55,因此它最早出現(xiàn)了分離。類似地,每種翼型在出現(xiàn)分離時的前緣最大馬赫數(shù)相比其他翼型是最高的,過了該時刻前緣分離泡逐漸發(fā)展成分離渦,前緣最大馬赫數(shù)均有所下降。

表1 出現(xiàn)分離時的攻角以及此時前緣吸力峰處最大馬赫數(shù)對比

為了深入地研究聯(lián)合射流在此深失速狀態(tài)下的控制機理,各翼型在上仰時攻角分別為23.37°、25.00°和下俯時攻角分別為24.14°、21.71°、20.02°時的流場圖如圖8所示。雖然RANS模擬出的流場渦只是平均大尺度結(jié)構(gòu),但對其非定常脫落過程進(jìn)行分析仍有助于理解失速特性。

從圖8可以看出:

隨著上仰攻角的增大,各翼型的前緣分離泡逐漸發(fā)展擴大成分離渦①,如圖8(a)所示,帶標(biāo)號的渦均是指時間渦。此時α=23.37°,可以看出:此時的Cμ越大,射流對分離渦的控制力度越大,使其向后緣發(fā)展的程度也更明顯。

在上仰至最大攻角α=25.00°時,如圖8(b)所示,翼型均失速,只是失速程度有所不同。此時CFJ翼型的流場也比原始翼型的稍復(fù)雜,相比原始翼型正在脫落的順時針渦①,Cμ=0.06的CFJ翼型的順時針渦①已經(jīng)完全脫落,Cμ=0.08的CFJ翼型的順時針渦①不僅完全脫落,且已在逆時針渦②的前面形成了新的順時針渦③,Cμ=0.10的CFJ翼型的順時針渦①也已經(jīng)完全脫落,且由于較強的射流將新形成的順時針渦③擠壓在逆時針渦②前面發(fā)展演變,渦的整體位置也更加靠后。

下俯至α=24.14°時,如圖8(c)所示,原始翼型的升力系數(shù)跌至最低點附近,此時CFJ翼型的升力系數(shù)已到達(dá)最低點后的峰值點附近,且射流動量系數(shù)越大,該峰值點越大,升力系數(shù)從最低點陡增到該峰值點的時間也越短。此時原始翼型原來的順時針渦①已經(jīng)脫落,原來的逆時針渦②也即將脫落,新形成的順時針渦③向后緣運動的同時逐漸與前緣一系列小渦合并。與原始翼型相比,CFJ翼型的順時針渦③都已運動到尾緣,Cμ=0.10的CFJ翼型的順時針渦③更是已經(jīng)脫落,并且所有CFJ翼型的尾緣處又出現(xiàn)新的逆時針渦④,射流強度越大,該渦的發(fā)展越快,同時噴氣口附近的一系列小渦擴散得也更快。

下俯至α=21.71°時,如圖8(d)所示,Cμ=0.10的CFJ翼型上表面的所有渦均已脫落完畢,流動基本變?yōu)楦街鳎ο禂?shù)從陡增峰值點降到新的極值點;Cμ=0.06和Cμ=0.08的CFJ翼型的升力系數(shù)還未降到新極值點,原來的順時針渦③進(jìn)一步發(fā)展擴大,且原來的逆時針渦④已經(jīng)脫落形成新的逆時針渦⑤;而此時原始翼型原來的順時針渦③才開始脫落,并開始形成新的逆時針渦④。

下俯至α=20.02°時,如圖8(e)所示,Cμ=0.08的CFJ翼型上表面的所有渦基本快要脫落完畢,只剩下最后一個剛開始脫落的逆時針渦⑤,流動開始逐漸轉(zhuǎn)為附著流,升力系數(shù)也從陡增峰值點降到了新的極值點,而此時Cμ=0.06的CFJ翼型雖然也完成了順時針渦③的脫落并剩下一個尚未脫落的逆時針渦⑤,但它的升力系數(shù)還未降到新的極值點,還會繼續(xù)形成新的順時針渦。此時原始翼型的流場較上一個狀態(tài)無太大變化,只是渦的尺度隨著運動略微變大。

此后Cμ=0.06的CFJ翼型上表面的逆時針渦⑤脫落以后形成新的順時針渦⑥,待該渦脫落以后流動逐漸轉(zhuǎn)為附著流,當(dāng)所有CFJ翼型上表面的渦都脫落完畢時,原始翼型的順時針渦③和逆時針渦④還一直保持緩慢脫落,待順時針渦③先脫落后還會再形成最后一個新的順時針渦⑤,待該渦脫落以后流動約在α=9°左右也轉(zhuǎn)為附著流。

綜合整個俯仰過程,可以看出,聯(lián)合射流在深失速狀態(tài)下對翼型動態(tài)失速有著很好的控制效果。雖然它提前誘導(dǎo)了前緣分離泡的出現(xiàn),但也加速了翼型上表面分離渦的運動和脫落,對比不同時刻的流場圖可以明顯看到CFJ翼型上表面失速渦的脫落頻率要比原始翼型快得多,使得CFJ翼型在失速后能更快地達(dá)到陡增峰值點,縮短整個失速周期,更快地恢復(fù)到附著流狀態(tài),且射流強度越大,該效果越明顯。

(a1) 基準(zhǔn)翼型 (a2)Cμ=0.06 (a3)Cμ=0.08 (a4)Cμ=0.10

(a)α=23.37°,上仰

(b1) 基準(zhǔn)翼型 (b2)Cμ=0.06 (b3)Cμ=0.08 (b4)Cμ=0.10

(b)α=25°,最大攻角

(c1) 基準(zhǔn)翼型 (c2)Cμ=0.06 (c3)Cμ=0.08 (c4)Cμ=0.10

(c)α=24.14°,下俯

(d1) 基準(zhǔn)翼型 (d2)Cμ=0.06 (d3)Cμ=0.08 (d4)Cμ=0.10

(d)α=21.71°,下俯

(e1) 基準(zhǔn)翼型 (e2)Cμ=0.06 (e3)Cμ=0.08 (e4)Cμ=0.10

(e)α=20.02°,下俯

圖8 俯仰過程中不同時刻的流場對比

Fig.8 Comparison of flow fields at different AoAs during pitching process

兩個時刻各翼型的壓強系數(shù)分布曲線如圖9所示。

(a) α=23.37°,上仰

(b) α=21.71°,下俯

從圖9可以看出:由于射流的作用翼型前緣吸力峰附近有更強的低壓區(qū),在上仰α=23.37°時,壓強系數(shù)波動較大的區(qū)間即為分離渦所在的位置,顯然較大Cμ的分離渦要更加靠近后緣;在下俯階段,隨著攻角的減小與失速渦的不斷脫落,翼型壓強系數(shù)分布的區(qū)別也越來越小。

力矩回線之間的面積是一個氣動阻尼指標(biāo),回線的走向確定系統(tǒng)是正阻尼(逆時針)還是負(fù)阻尼(順時針)。而阻尼指標(biāo)直接關(guān)系到系統(tǒng)的穩(wěn)定性,當(dāng)負(fù)的氣動阻尼(即力矩系數(shù)遲滯閉環(huán)中順時針部分所包含的區(qū)域)在整個振蕩周期中超過正的氣動阻尼時,翼型會失去穩(wěn)定性而引起失速顫振,從而導(dǎo)致結(jié)構(gòu)的疲勞損傷和系統(tǒng)失穩(wěn)。因此動態(tài)失速控制的主要目標(biāo)是:既要保持動態(tài)失速能使最大升力增加的優(yōu)勢,同時又要大幅降低阻力、負(fù)的俯仰力矩峰值以及負(fù)的氣動阻尼[49-50]。一個俯仰周期內(nèi)升力系數(shù)峰值的變化情況如表2所示,可以看出:在此深失速狀態(tài)下,施加聯(lián)合射流控制以后,最大升力系數(shù)有不同程度的提高,最大阻力系數(shù)和力矩負(fù)峰得到不同程度的降低,且射流動量系數(shù)越大,降低的程度也越大。同時,從圖7(c)可以看出:原始翼型在上仰和下俯時的力矩曲線發(fā)生交疊,出現(xiàn)了一定區(qū)域的負(fù)阻尼,而CFJ翼型則均未出現(xiàn)負(fù)阻尼,整個回線均為正的氣動阻尼。

表2 一個振蕩周期內(nèi)氣動力系數(shù)峰值的變化

3 結(jié) 論

(1) 在聯(lián)合射流關(guān)閉的情況下,射流通道對動態(tài)失速特性有一定影響,主要體現(xiàn)在翼型上仰時的附著流階段,而對處于失速分離階段的氣動特性影響很小。

(2) 在打開聯(lián)合射流的情況下,動態(tài)失速特性得到了極大的改善,升阻力系數(shù)遲滯環(huán)和力矩系數(shù)遲滯環(huán)的面積均顯著減小,升力系數(shù)大幅提高,阻力系數(shù)顯著減小,且阻力系數(shù)和力矩系數(shù)曲線的峰值顯著減小。此外,聯(lián)合射流可以完全消除基準(zhǔn)翼型力矩系數(shù)曲線所反映出的負(fù)阻尼區(qū)域,使得整個俯仰周期內(nèi)力矩系數(shù)曲線均表現(xiàn)為正阻尼。

(3) 聯(lián)合射流方法對于翼型動態(tài)失速問題具有很好的控制效果和應(yīng)用前景。

主站蜘蛛池模板: 国产午夜看片| 国产高清精品在线91| 午夜精品影院| 四虎永久免费地址在线网站 | 国产成人AV综合久久| 九色免费视频| 一级一毛片a级毛片| jizz国产在线| 国产精品久久自在自线观看| 丰满人妻久久中文字幕| 伊人久久大香线蕉影院| 亚洲人成人无码www| 中文字幕调教一区二区视频| 超级碰免费视频91| 国产精品无码AⅤ在线观看播放| 白丝美女办公室高潮喷水视频| 欧美日韩国产高清一区二区三区| 99久久国产综合精品2020| 美女免费黄网站| 一级毛片在线免费看| 国产精品美女免费视频大全| 国产成人综合亚洲网址| 波多野结衣国产精品| 亚洲日本一本dvd高清| 国产亚洲精品自在久久不卡| 99久久精品久久久久久婷婷| 国产成人综合久久精品下载| 亚洲人成网线在线播放va| 久久黄色小视频| 亚洲综合中文字幕国产精品欧美| 国产95在线 | 国产欧美高清| 亚洲精品无码不卡在线播放| 成人国产小视频| 久久96热在精品国产高清| 久操线在视频在线观看| 亚洲a级毛片| 亚洲日韩欧美在线观看| 国产综合另类小说色区色噜噜| 亚洲精品少妇熟女| 亚洲午夜片| 国产99免费视频| 久久a毛片| 国产在线精品香蕉麻豆| 日韩黄色精品| 在线国产91| 国产 日韩 欧美 第二页| 成人年鲁鲁在线观看视频| 久久久精品无码一区二区三区| 全部毛片免费看| 精品国产毛片| 拍国产真实乱人偷精品| 亚洲一区二区三区国产精华液| 第一区免费在线观看| 日本日韩欧美| 日韩精品久久久久久久电影蜜臀| 国产午夜一级淫片| 久久夜色精品| 欧美色图第一页| 欧美成人午夜在线全部免费| 久久午夜夜伦鲁鲁片无码免费| 亚洲资源在线视频| 欧美亚洲国产一区| 日韩欧美中文| 91在线丝袜| 亚洲日本www| www.99精品视频在线播放| 手机在线国产精品| 国产福利微拍精品一区二区| 国产欧美成人不卡视频| 日本亚洲欧美在线| 不卡国产视频第一页| 久久精品娱乐亚洲领先| 亚洲自拍另类| 久久婷婷色综合老司机| 大陆精大陆国产国语精品1024| 亚洲女同一区二区| 久久免费视频6| 国产精品va免费视频| 一级香蕉视频在线观看| 欧美亚洲第一页| 亚洲天堂网在线观看视频|