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

不同翼型厚度和來流馬赫數下的槳渦干擾噪聲分析

2023-09-14 11:09:38相倩余啟梁王鋼林劉勇葉海劍
航空科學技術 2023年6期

相倩 余啟梁 王鋼林 劉勇 葉海劍

摘 要:槳渦干擾噪聲是直升機氣動噪聲主要組成之一,為了正確預測和降低直升機噪聲,必須開展氣動噪聲相關物理參數研究。在對聲場進行計算流體力學(CFD)直接數值模擬的基礎上,分析了不同厚度和來流馬赫數下二維平行槳渦干擾噪聲傳播特性和聲源位置,分析了翼型厚度和來流馬赫數對槳渦干擾噪聲的影響,并得到了可壓縮情況下遠場聲壓預測公式。研究表明,低馬赫數下,翼型厚度對噪聲指向性影響不大,高馬赫數下,翼型厚度對噪聲指向性影響程度增大;噪聲強弱主要隨來流馬赫數變化,翼型厚度對其影響較小;翼型厚度和來流馬赫數變化不會改變聲源點位置。開展不同翼型厚度和來流馬赫數下的槳渦干擾噪聲分析可以為進一步了解并控制直升機槳渦干擾噪聲提供一定的參考。

關鍵詞:翼型厚度; 來流馬赫數; 二維平行BVI; 氣動噪聲; 直接法

中圖分類號:V211 文獻標識碼:A DOI:10.19452/j.issn1007-5453.2023.06.002

基金項目:國家自然科學基金(11962018,12262023)

在直升機低速下降、機動等飛行過程中,前行槳葉產生的槳尖渦脫落后會靠近或通過后行槳葉,形成槳渦干擾(BVI)現象。當槳尖渦軸線與槳葉平行時,可以簡化為二維平行BVI。BVI噪聲一旦出現,會顯著增大直升機的總體噪聲水平,帶來嚴重的噪聲污染[1]。因此,對影響BVI氣動噪聲的物理參數的研究,一直是預測和降低直升機噪聲的前提[2]。

先前,在理論和試驗方面研究者們已開展了許多與影響BVI氣動噪聲的物理參數相關的研究工作[3-7]。孟曉偉等[3]通過建立懸停狀態直升機傾斜式尾槳誘導速度的方法有效地對尾槳BVI進行了預測。趙俊等[4]基于計算流體力學(CFD)方法對多旋翼飛行器氣動力和噪聲特性進行了分析。H. Y. Yung等[5]通過理論計算、全機試驗和風洞試驗等對BVI噪聲聲壓級和飛行狀態影響進行了研究,但對于影響BVI的物理參數、噪聲波的形成和傳播過程等方面沒有深入研究。F. Caradonna等[6]在風洞中用旋轉機翼研究了BVI噪聲,分析了機翼與渦流相互作用,比較了Kirchhoff方法和FW-H方法等各種預測BVI噪聲的方法。遠場壓力的計算結果和測量結果對比表明,當渦特征(渦核位置、渦強度和渦核半徑)已知時,可以較準確地預測BVI噪聲。C. Kitaplioglu 等[7]為減少旋翼轉動時自身尾流的影響,利用外部產生的渦流與表面裝有壓力傳感器的靜止旋翼相互作用,并通過固定的近場麥克風和可移動的遠場麥克風陣列收集聲學數據。試驗研究了翼型與渦干擾垂直間距、渦核強度以及葉尖馬赫數等幾個影響BVI的參數,該試驗方法雖準確可靠,但試驗周期長且成本高。

旋翼BVI現象復雜且隨機性強,在旋翼各個角度和展向位置上可能同時發生BVI現象,通過試驗方式測量噪聲存在局限。隨著計算機的發展,CFD方法逐漸成為直升機流場分析最有效的方法之一[8]。J. M. Chen等[9]通過在渦流發生器的下游放置NACA0012翼型,對二維BVI進行了試驗研究,研究結果表明,當旋渦與非升力翼型相互作用時,它在前緣附近產生相當劇烈的局部表面壓力變化。對于迎角參數的影響,該研究者先后試驗了迎角為5°和10°的工況,試驗結果顯示,相比于5°迎角,10°迎角下翼型上的瞬態升力幅度更大,瞬態升力的幅值也隨著渦流發生器和機翼之間垂直間隔的減小而增大。I. Marcel等[10-12]使用大渦模擬(LES)方法對不同迎角、不同渦特征以及翼型與渦不同干擾垂直間距下的BVI流場進行了數值模擬,減小了數值耗散,并捕捉到了干擾時的流場信息,研究結果表明,氣動系數的大小隨著干擾垂直間距的增加和渦核尺寸的減小而減小,LES方法為獲得非定常壁面壓力場提供了一種很有前景的工具。

在國內,史勇杰等[13]通過Euler方程分析BVI流場的特征,將自由渦在流場的誘導速度等效成網格畸變運動,并入真實網格運動中保持渦核結構不變。描述了干擾流場的升力變化特性,然后分析了渦核強度、翼型與渦干擾垂直間距和干擾角度等參數對BVI流場特性的影響。喬渭陽等[14]采用有限翼展槳葉模型計算葉片表面的非定常力,對渦強度、渦核半徑、翼型與渦干擾的垂直間距和夾角、馬赫數等影響參數進行實例計算。聲場的計算使用了聲比擬法。結果表明,該方法得出的BVI噪聲聲壓數據與試驗結果吻合良好。

目前,對于影響BVI氣動噪聲的物理參數研究主要圍繞流場分析,而物理參數對BVI氣動噪聲產生和傳播方面的影響還需要進一步計算和分析。

BVI氣動噪聲生成過程是一個強烈的非線性過程[15],對于氣動噪聲的數值計算,大多采用聲比擬方法預測遠場噪聲;而LBM-LES相結合的直接計算氣動噪聲方法,不僅能準確預測遠場噪聲,還能捕捉氣動噪聲的生成、傳播和衰減實時變化過程,有利于認識不同翼型厚度和來流馬赫數對BVI脈沖氣動噪聲的影響問題,所以本文采用LBM-LES氣動噪聲直接計算方法,對不同翼型厚度和來流馬赫數下,二維平行BVI產生的脈沖非線性氣動噪聲進行直接計算。

1 數值方法及邊界條件

1.1 直接法

直接法計算氣動噪聲具有同時得到精確流場和聲場的優點,可用于對氣動噪聲進行傳播過程、產生機理等研究,但需要采用高精度的時空離散格式和無反射邊界條件兩個關鍵技術來實施。

1.2 LBM-LES方法

格子玻耳茲曼(LBM)方法可高效求解非定常低速不可壓的湍流問題[16],多松弛時間(MRT)-LBM方法具有四階空間精度和二階時間精度,時間離散采用顯式格式,具有數值穩定性高、低耗散、精度高等特性[17]。在相同的耗散情況下,該方法計算效率高,因此逐漸發展出了基于LBM直接法計算氣動噪聲。

LES方法是通過亞格子模型求解小尺度渦,對含能較高的大尺度渦則直接求解。BVI噪聲屬于中低頻窄帶噪聲,具有含能高且尺度大的特點,和高雷諾數低馬赫數湍流噪聲相比,對計算精度、時間尺度、空間尺度、計算機及耗散性等的要求可以降低,因而比較容易實施,這使采用LES方法直接計算BVI近場和遠場噪聲成為可能。

1.3 聲場計算方法

因為流場與聲場的基本方程是相同的,所以可通過求解LBM方程直接得出流場和聲場的統一解。

1.4 網格分布及邊界條件

由于BVI干擾主要作用于翼型附近區域,所以采用多尺度模式來劃分計算網格。將計算區域劃分成壁面附近、近場、遠場和聲波吸收4個區域。圖1所示為區域劃分及網格分布示意圖,圖1(a)翼型附近區域放大為圖1(b)。計算域大小為80C×80C,其中C為翼型弦長。

上述4個邊界均采用求解局部一維無黏方程[20-21]實現無反射邊界條件,避免了邊界處波向聲場內反射,遠場外有網格尺度較大的聲波吸收區域,LES方法可過濾掉來自邊界的小尺度偽反射波。

2 計算結果與分析

2.1 網格無關性驗證

首先進行網格無關性分析,以保證計算結果的可靠性?;贐VI引起的脈沖聲壓峰值建立網格數目影響分析。表1為網格數目和尺寸,N表示網格數目。

圖2表示三種不同網格數目在監測點(0,20C)處的聲壓時間歷程。由圖2可以發現,N=107萬和N=220萬之間的脈沖聲壓曲線更為接近,而N=70萬的脈沖聲壓曲線與另兩條脈沖聲壓曲線存在差異;隨著網格數目的增加,N=220萬和N=107萬的聲壓峰值差異更小??紤]到節省計算資源以及計算精度,本文采用N=107萬網格進行數值計算。

2.2 噪聲傳播

通過直接法計算氣動噪聲可以直接獲得聲場信息,一般用脹量Θ=??V顯示聲波[23]:脹量值大小可表示聲波強弱,脹量值最大值連線可表示輻射方向。圖3所示脹量云圖表示不同翼型厚度和來流馬赫數下槳渦干擾脈沖噪聲聲波向遠場傳播的過程,其中相同來流馬赫數下脈沖聲波截取時刻相同,即圖3(a)~圖3(c)為相同時刻脹量云圖,圖3(d)~圖3(f)為相同時刻脹量云圖,圖3(g)~圖3(i)為相同時刻脹量云圖。以圖3(a)為例,t1~t4代表脈沖噪聲聲波傳播的不同時刻;角度表示聲波上下輻射角度值(輻射角以翼型下游對應的方向為0°,逆時針為正);正負數字表示t1時刻聲波沿輻射角方向脹量值,即t1時刻脹量最大值。

由圖3可知,不存在因為網格分區和遠場邊界產生與脈沖噪聲同量級強度偽反射波;從上下遠場脹量的變化規律可知,上下遠場噪聲幅值大小相近,但相位相反,脹量圖呈8字形,說明聲場的輻射特性符合偶極子聲源的輻射特性;上遠場輻射角絕對值(簡稱上輻射角)和下遠場輻射角絕對值(簡稱下輻射角)表示脈沖噪聲聲源指向性;上遠場脹量絕對值(簡稱上脹量值)均略小于下遠場脹量絕對值(簡稱下脹量值),說明來流渦為順時針旋轉時,上源場噪聲比下源場噪聲稍弱,是旋渦的下洗導致;相同馬赫數下比較,隨著厚度的增加,相同時刻相同位置處脹量絕對值減小,說明聲波隨厚度增加而變弱,但變化幅度不大。

如圖3(a)~圖3(b)所示,在來流馬赫數Ma為0.14時,t1時刻脹量絕對值大小在0.3~0.43范圍內變化;上輻射角在82°附近變化,下輻射角在89°附近變化,上下輻射角有明顯差別;上下輻射角隨著厚度增加而增大,但變化不明顯。如圖3(d)~圖3(f)所示,來流馬赫數增加為0.4時,t1時刻脹量絕對值大小在1.84~2.26范圍內變化;上、下輻射角均在77°附近,隨厚度增加變化不明顯。如圖3(g)~圖3(i)所示,來流馬赫數Ma=0.6時,t1時刻脹量絕對值大小在2.59~3.03范圍內變化;同Ma情況一樣,Ma上下輻射角大小相近,說明隨著來流馬赫數增加,渦的旋轉方向對上、下輻射角的影響變弱;隨著厚度增加,上、下輻射角由63°增至70°,而在Ma=0.14和Ma=0.4情況下,隨著厚度增加,上、下輻射角則基本不變。這表明來流馬赫數為0.6時,翼型厚度對輻射角的影響增強,聲波指向性與厚度變化情況和來流速度有關。

對相同翼型厚度、不同來流馬赫數脹量云圖進行比較可知,隨著馬赫數增加,相同位置處脹量值絕對值增加明顯,說明來流馬赫數對聲波強弱影響較大;隨著來流馬赫數增大,聲波輻射角趨近于90°,且低馬赫數下,上、下輻射角差別在7°左右,高馬赫數下,上、下輻射角差別在1°左右,說明隨著來流馬赫數增加,渦的旋轉方向對輻射角的影響變弱。

綜上可知,來流馬赫數和翼型厚度可以影響槳渦干擾脈沖噪聲強弱和方向。脈沖噪聲強弱隨來流馬赫數增加而增強,且變化幅度較大;隨翼型厚度增加而減弱,且變化幅度較小。脈沖噪聲指向隨來流馬赫數增加趨向于來流方向,隨厚度增加趨向于來流反方向;翼型上下聲源強弱相近,且來流渦為順時針旋轉時,上源場噪聲比下源場噪聲弱;來流馬赫數增加,渦的旋轉方向對噪聲輻射角的影響變弱。

2.3 聲源聲壓大小和位置分析

翼型前緣附近壓強發生強度不同的兩次突變,對應峰值時刻的渦量—壓強差云圖隨翼型厚度和來流馬赫數的變化如圖4、圖5所示。圖4為第一次壓強變化時刻渦量和壓力差云圖,由圖 4可知,相同翼型厚度下,聲源附近壓力差隨著來流馬赫數增加而增加;相同來流馬赫數下,聲源附近壓力差隨翼型厚度增加而減小,變化不大;來流馬赫數改變,聲源附近壓力差隨厚度變化規律相同;翼型厚度改變,聲源附近壓力差隨來流馬赫數變化規律相同;翼型厚度和來流馬赫數對聲源位置均沒有影響,在翼型弦長的10%附近。由圖5可知,聲源附近壓力差隨來流馬赫數和翼型厚度變化規律和聲壓正峰值時刻相同,聲源位置依然位于翼型弦長的10%附近,且不隨翼型厚度和來流馬赫數變化。

對比圖4和圖5可知,不同工況下,聲壓負峰值時刻聲源附近壓力差均更大,噪聲更強,此時旋渦位于翼型最大厚度點附近。

2.4 BVI遠場聲壓預測公式

在此基礎上考慮壓縮性和翼型厚度影響,由于壓縮性和翼型厚度直接影響輻射角方向,可以用輻射角方向表示壓縮性和厚度對遠場聲壓峰值的影響。代入計算得到的不同翼型厚度、不同來流馬赫數情況下升力系數波動幅值和輻射角方向遠場聲壓峰值進行線性擬合。如圖6所示。

3 結論

基于LBM-LES方法,對二維平行槳葉的槳渦干擾噪聲的聲場進行了直接計算,得到了精確的槳渦干擾脈沖氣動噪聲聲場數據及非定常流場數據。在此基礎上,分析了翼型厚度和來流馬赫數對槳渦干擾噪聲的影響,得出以下結論:

(1)隨著來流馬赫數的變化,噪聲強度和翼型厚度的變化對槳渦干擾噪聲的影響不同。噪聲強度隨來流馬赫數的變化程度遠大于隨翼型厚度的變化程度,這可能是由于來流馬赫數對氣動聲學特性的影響更為顯著。此外,隨著翼型厚度的增加,槳渦干擾噪聲會相應減小。

(2)噪聲傳播方向的變化規律是本文研究的另一個重要結果。在低馬赫數下,厚度變化對槳渦干擾噪聲的指向性影響不大,說明厚度對氣動聲學特性的影響較小。但是,在高馬赫數下,厚度對指向性影響程度增加,指向靠近來流的法線方向,這可能是由于隨著來流馬赫數的增加,氣動聲學特性的變化更加顯著,導致噪聲傳播方向的變化也更加明顯。

(3)進一步研究發現,噪聲傳播方向隨來流馬赫數增加而靠近來流方向,而隨著厚度的增加,隨來流馬赫數增加的程度降低。這表明隨著來流馬赫數的增加,氣動聲學特性對噪聲傳播方向的影響逐漸增強,而翼型厚度的變化對其影響程度降低。

本文研究得到不同來流馬赫數下的遠場聲壓峰值預測公式,這為槳渦干擾噪聲控制提供了實用的工程手段。需要指出的是,由于受試驗條件的限制和氣動聲學特性的復雜性,預測公式的適用范圍有限,需要在實際應用中加以驗證。

參考文獻

[1]馮劍波,陸洋. 直升機旋翼槳渦干擾噪聲主動控制技術綜述[J]. 噪聲與振動控制,2018,38(3):1-9. Feng Jianbo, Lu Yang. Review of active control techniques of blade vortex interaction noise for helicopter rotor[J]. Noise and Vibration Control, 2018, 38(3): 1-9. (in Chinese)

[2]鄧景輝. 直升機技術發展與展望[J]. 航空科學技術,2021,32(1): 10-16. Deng Jinghui. Development and prospect of helicopter technology[J]. Aeronautical Science & Technology, 2021, 32(1): 10-16. (in Chinese)

[3]孟曉偉,張宏林,楊文鳳. 直升機傾斜式尾槳渦環預測與試飛研究[J]. 航空科學技術,2018,29(1): 58-62. Meng Xiaowei, Zhang Honglin, Yang Wenfeng. Theoretic calculation and flight test for the helicopter tilting tail rotor vortex ring[J]. Aeronautical Science & Technology, 2018, 29(1): 58-62. (in Chinese)

[4]趙俊,李志彬.變距四旋翼飛行器氣動力及噪聲特性計算研究[J]. 航空科學技術,2022,33(4): 57-66. Zhao Jun, Li Zhibin. Computational research on aerodynamic and noise characteristics of variable pitch quadrotor aircraft[J]. Aero‐nautical Science & Technology, 2022, 33(4): 57-66. (in Chinese)

[5]Yung H Y. Rotor blade-vortex interaction noise[J]. Progress in Aerospace Sciences, 2000, 36(2): 97-115.

[6]Caradonna F, Kitaplioglu C, McCluer M, et al. Methods for the prediction of blade-vortex interaction noise[J]. Journal of the American Helicopter Society, 2000, 45(4): 303-317.

[7]Kitaplioglu C, Caradonna F X, Burley C L. Parallel blade‐vortex interactions: An experimental study and comparison with computations[J]. Journal of the American Helicopter Society, 1997, 42(3): 272-281.

[8]Lele S. Computational aeroacoustics-a review[C]//35th Aero‐space Sciences Meeting and Exhibit, 1997: 18.

[9]Chen J M, Chang D M. Unsteady pressure measurements for parallel vortex-airfoil interaction at low speed[J]. Journal of Aircraft, 1997, 34(3): 330-336.

[10]Marcel I. Numerical study of helicopter blade-vortex mechanism of interaction using large-eddy simulation[J]. Computers & Structures, 2009, 87(11-12): 758-768.

[11]Marcel I, Nitzsche F, Matida E. Influence of vortex characteris‐tics on the blade-vortex mechanism of interaction using largeeddy simulation[C]. 47th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, 2009: 1282.

[12]Marcel I, Nitzsche F, Matida E. Influence of angle of attack on the helicopter blade-vortex mechanism of interaction using large eddy simulation[C]. 46th AIAA Aerospace Sciences Meeting and Exhibit, 2008: 551.

[13]史勇杰,招啟軍,徐國華. 旋翼槳—渦干擾氣動特性計算及參數影響研究[J]. 航空學報,2010,31(6): 1106-1114. Shi Yongjie, Zhao Qijun, Xu Guohua. Numerical calculation and paramertric study of aerodynamics of rotor blade-vortex interaction[J]. Acta Aeronautica et Astronautica Sinica, 2010, 31(6): 1106-1114. (in Chinese)

[14]喬渭陽,唐狄毅,李文蘭. 旋翼BVI噪聲的理論模擬與分析[J]. 航空學報,1994,15(6): 725-730. Qiao Weiyang, Tang Diyi, Li Wenlan. Theoretical study on the blade-vortex interaction noise of helicopter rotor[J]. Acta Aeronautica et Astronautica Sinica, 1994, 15(6): 725-730. (in Chinese)

[15]尹堅平,胡章偉. 用測量的葉片表面非定常壓力預估直升機旋翼槳渦干擾(BVI)噪聲[J]. 航空學報, 1996,17(5): 87-89. Yin Jianping, Hu Zhangwei. Prediction of rotor blade vortex interaction noise from measured unsteady pressure of blade[J]. Acta Aeronautica et Astronautica Sinica, 1996, 17(5): 87-89.(in Chinese)

[16]Brionnaud R, Chávez M, Trapani G, et al. Direct noise compu‐tation with a Lattice Boltzmann method and application to in‐dustrial test cases[C]. 22nd AIAA/CEAS Aeroacoustics Confer‐ence, 2016: 2969.

[17]Marié S, Ricot D, Sagaut P. Comparison between Lattice Boltzmann method and Navier-Stokes high order schemes for computational aeroacoustics[J]. Journal of Computational Physics, 2009, 228(4): 1056-1070.

[18]Chen S, Doolen G D. Lattice Boltzmann method for fluid flows[J]. Annual Review of Fluid Mechanics, 1998, 30(1): 329-364.

[19]馬致遙,單鋒,章東. 格子玻爾茲曼方法對不同張角聚焦聲束的建模[J]. 聲學學報,2018,43(2): 217-223. Ma Zhiyao, Shan Feng, Zhang Dong. Modeling of focused acoustic beams with varied aperture angles using Lattice Boltzmann method[J]. Acta Acustica, 2018, 43(2): 217-223. (in Chinese)

[20]Heubes D, Bartel A, Ehrhardt M. Characteristic boundary conditions in the Lattice Boltzmann method for fluid and gas dynamics[J]. Journal of Computational and Applied Mathematics, 2014, 262: 51-61.

[21]Stoll S J B. Lattice Boltzmann simulation of acoustic fields, with special attention to non-reflecting boundary conditions[Z]. Institutt for Elektronikk Telekommunikasjon, 2014.

[22]Loiodice S, Drikakis D, Kokkalis A. Influence of vortex models on the prediction of 2D airfoil vortex interaction[C]. 34th European Rotorcraft Forum 2008,2008: 360-369.

[23]Colonius T, Lele S K, Moin P. Sound generation in a mixing layer[J]. Journal of Fluid Mechanics, 1997, 330: 375-409.

[24]鄒森,劉勇,王琦. 二維平行槳渦干擾氣動噪聲特性分析[J].聲學學報,2020,45(4): 587-593. Zou Sen, Liu Yong, Wang Qi. Analysis on aeroacoustic characteristic of the two-dimensional parallel blade-vortex interaction[J]. Acta Acustica, 2020, 45(4): 587-593.(in Chinese)

Analysis on Rotor Vortex Interference Noise Under Different Airfoil Thickness and Free-stream Mach Number

Xiang Qian1, Yu Qiliang2, Wang Ganglin1, Liu Yong2, Ye Haijian2

1. Chinese Aeronautical Establishment, Beijing 100029, China

2. Nanchang Hangkong University, Nanchang 330063, China

Abstract: Blade-vortex interaction noise is one of the main components of helicopter aerodynamic noise. In order to correctly predict and reduce helicopter noise, it is necessary to conduct research on physical parameters related to aerodynamic noise. Based on the direct numerical simulation of the acoustic field using Computational Fluid Dynamics (CFD), the propagation characteristics and sound source location of two-dimensional parallel propeller blade-vortex interference noise under different thicknesses and free-stream Mach numbers are analyzed. The effects of airfoil thickness and free-stream Mach number on propeller blade-vortex interference noise are analyzed, and the prediction formula for far field sound pressure in compressible cases is obtained.The research shows that at low Mach number, the influence of airfoil thickness on noise directivity is not significant, while at high Mach number, the influence of airfoil thickness on noise directivity increases;The noise intensity mainly varies with the free-stream Mach number, and the influence of airfoil thickness is small; Variations of airfoil thickness and free-stream Mach number do not change the location of the sound source point.This research provides a certain reference value for further study on controlling the blade-vortex interaction noise of helicopters.

Key Words: airfoil thicknesses; free-stream Mach number; two-dimensional parallel BVI; aeroacoustics; direct simulation method

主站蜘蛛池模板: 中国国产A一级毛片| 91 九色视频丝袜| 一本久道热中字伊人| 91福利一区二区三区| 国产精品流白浆在线观看| 欧美中文一区| 色综合激情网| 久久大香伊蕉在人线观看热2| 人妻无码一区二区视频| 97在线视频免费观看| 老色鬼欧美精品| 亚洲中文字幕23页在线| 免费在线成人网| 国产一二三区在线| 国产AV无码专区亚洲精品网站| 麻豆AV网站免费进入| 亚洲中文字幕精品| 2021国产v亚洲v天堂无码| 亚洲无码高清一区二区| 在线另类稀缺国产呦| 亚洲国产欧美自拍| 992Tv视频国产精品| 亚洲制服中文字幕一区二区 | 亚洲欧美不卡| 91麻豆久久久| 免费一级成人毛片| 91精品aⅴ无码中文字字幕蜜桃| 欧美午夜小视频| 国产精品欧美亚洲韩国日本不卡| 黄色网在线| 3p叠罗汉国产精品久久| 好紧好深好大乳无码中文字幕| www.youjizz.com久久| 熟妇丰满人妻av无码区| 国产精品一老牛影视频| 免费激情网址| 亚洲国产av无码综合原创国产| 欧美精品不卡| 亚洲成人黄色在线| 中文毛片无遮挡播放免费| 欧美成人免费一区在线播放| 国产精品爽爽va在线无码观看| 国产香蕉在线视频| 日本久久久久久免费网络| 尤物国产在线| 第九色区aⅴ天堂久久香| 专干老肥熟女视频网站| 香蕉国产精品视频| 中文天堂在线视频| 久久综合色88| 亚洲欧美另类久久久精品播放的| 国产成人精品亚洲77美色| 精品成人一区二区三区电影| 欧美97色| 日韩高清无码免费| 国产三区二区| 波多野结衣亚洲一区| 久久国产精品夜色| 在线va视频| 国产97视频在线观看| 无码中文字幕乱码免费2| 一级毛片不卡片免费观看| 99久久精彩视频| 一区二区午夜| 91人妻日韩人妻无码专区精品| 伊人久久大线影院首页| 伊人色综合久久天天| av色爱 天堂网| 热久久这里是精品6免费观看| 97视频在线观看免费视频| 怡红院美国分院一区二区| 国产麻豆福利av在线播放| 欧美亚洲激情| 国产精品99久久久久久董美香| 男女男精品视频| 人妻少妇乱子伦精品无码专区毛片| 99久久精品国产自免费| 91在线高清视频| 精品久久久久无码| 国产视频资源在线观看| 青青国产成人免费精品视频| 国产精品露脸视频|