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

單級軸流風扇單音噪聲源及噪聲輻射特性的數值研究

2015-08-03 07:28:14張建華楚武利
動力工程學報 2015年7期

張建華,楚武利,2,劉 煒

(1.西北工業大學動力與能源學院,西安710072;2.先進航空發動機協同創新中心,北京100191;3.華中科技大學機械科學與工程學院,武漢430074)

由于環境保護的需要,風機噪聲指標越來越引起人們的重視.近10年來旋轉機械噪聲成了一個熱門研究課題.目前,對軸流風機噪聲產生機理和預測方法的研究主要有3個方面:單轉子寬頻噪聲、級轉靜干涉單音噪聲和葉尖渦噪聲.Wright等[1-3]研究了軸流葉輪機械轉子氣動力與噪聲的關系,得到了單轉子軸流風機的噪聲產生機理:轉子葉片是主要的噪聲源,噪聲以寬頻為主,端壁邊界層和轉子葉尖的干涉形成的湍流尾跡脫落是噪聲產生的主要原因,氣動力作用于葉尖泄漏流推動了轉子尾跡流動和通道內的二次流流動.然而對于亞音速(Ma<0.3)單級風扇而言,級轉靜干涉噪聲是此類葉輪機械主要的離散單音噪聲分量,靜葉表面的非定常壓力脈動是此類噪聲的主要噪聲源,其來源于轉子尾跡對下游靜葉的周期性沖擊[4-5].這種非定常壓力脈動通常可以通過高效和高精度的商用CFD 軟件獲得.目前,對于軸流風扇[6-9]及軸流泵[10]噪聲的預測多采用混合氣動聲學方法(HCAA),該方法的本質是將復雜的噪聲計算分成2個獨立的過程,首先基于三維黏性非定常流場(URANS),采用商用CFD軟件計算獲得噪聲源,再采用Lighthill聲類比理論的FW-H 方程耦合URANS聲源信息獲得遠場噪聲輻射,但這種方法并沒有考慮實際存在的外殼體對聲波的作用,誤差較大.為了獲得更準確的預測結果,Lee 等[11-12]在預測中引入kirchhoff-helmholtz BEM 技術,這種技術在計算中考慮了聲源外部固體邊界對聲傳播的散射和反射作用.Cai等[13]和Hu等[14]采用多區域邊界元方法(MDBEM)分別預測了工業用離心風機和單級軸流風機噪聲輻射,取得了較好的效果.然而上述聲學邊界元方法的本質是利用自由空間的格林函數積分離散Lighthill方程,目前的格林函數積分解法只能針對空間簡單幾何邊界求解,而對于空間復雜幾何邊界難以得到準確解,對于這類復雜邊界問題只能將其在自由空間內簡化求解,因此忽略了聲波在復雜固壁上的反射和散射等作用.為了克服這個缺陷,Koopmann 等[15]基于Lighthill方程的變分形式,采用有限元方法求解,這種方法可以將傳播中的聲源在聲場空間離散以考慮結構和聲的相互作用,適用于任何復雜的幾何結構.

對于離散單音噪聲而言,壁面非定常流場的壓力脈動是最為直接的噪聲源,為了快速獲取風扇的噪聲特性,Liu等[16]根據Lighthill聲類比理論中聲源信息和噪聲的關系,通過分析聲源區頻譜定性地預測了工業用離心風機的噪聲輻射特性.Tan等[17]應用此方法定性分析了離心壓縮機噪聲源的強度和位置.

筆者以單級軸流風扇為研究對象,采用數值模擬方法研究了該軸流風扇的非定常氣動力和氣動噪聲特性,主要分為3個部分:(1)噪聲源流場特性的分析;(2)基于噪聲源表面壓力脈動頻譜分析定性確定噪聲的強度和位置;(3)采用基于Lighthill聲類比理論的FW-H 方程耦合URANS方法計算風扇遠場噪聲輻射.為了考慮復雜壁面對聲波的散射和反射影響,聲場計算采用聲學有限元方法,即將FW-H 方程轉變為頻域Helmholtz弱積分形式并采用Galerkin方法離散.

1 風扇計算模型

1.1 風扇模型參數

以某單級軸流風扇為研究對象,具體參數見表1.

表1 風扇轉子和靜子參數Tab.1 Rotor/Stator parameters of the fan

1.2 流場計算網格及邊界條件

1.2.1 網格特點

基于三維設計軟件UG 建立風扇各結構的實體模型,并在Blade-Gen中分別生成風扇的轉子和靜子葉片通道,在Ansys Turbo-Grid中劃分通道的高質量六面體結構化網格,網格拓撲采用正交性較好的H-J-C-L型網格,葉片表面周圍采用O 型網格拓撲,近壁面網格節點加密處理,通道進口到葉片前緣采用J型網格拓撲,其余部分采用H 型網格拓撲,這樣既可以使網格正交性達到最佳又可以避免全部采用J型網格拓撲時前緣和尾緣網格長寬比過大.轉子葉頂間隙為2.5mm,采用H 型網格.轉子通道最小網格面角大于25°,網格長寬比小于100,由于風扇為軸對稱結構,為了節省計算資源,只計算一個級通道(轉子/靜子).級通道計算域分為3部分:進口流域及轉子通道計算域、靜子通道及擴壓段計算域和出口延伸段計算域.進口流域及轉子通道網格數為349 208,靜子通道及擴壓段網格數為290 904,出口延伸段網格數為68 080.圖1給出了轉子、靜子和輪轂的表面網格分布.為驗證網格無關性,將各計算域網格分別加密一倍,圖2為不同數量網格下風機的效率和全壓升對比圖.由圖2可知,網格增大為原來的2倍后全壓升和效率變化不大,說明此套網格已達到網格無關性要求,現用網格數量已滿足計算精度需要.

圖1 葉片及輪轂表面網格Fig.1 Surface mesh of the blade and hub

圖2 網格無關性驗證Fig.2 Grid independence validation

1.2.2 流場計算模型

采用商業CFD 軟件Ansys CFX 求解此單級風扇的流場,基于有限元體積法離散Navier-Stokes方程和RNGk-ε湍流方程[18].RNGk-ε湍流方程能準確預測旋轉和曲率對流動的影響,滿足所研究軸流風扇流場計算需要.空間離散采用高精度的高階格式差分,非定常計算的時間推進格式為雙時間步長全隱式格式,時間項使用二階后向歐拉差分離散,使用壓力耦合求解算法來求解連續性方程和動量方程.由于樣機葉片葉尖馬赫數為0.18(<0.3),認為氣流不可壓縮,流動被認為是絕熱的且不需要考慮進出口溫差,所以計算中舍棄了能量方程.

1.2.3 邊界條件

計算域介質為20 ℃、101 325Pa的空氣.邊界條件為:給定進口質量流量,出口靜壓邊界條件,出口壓力值為大氣壓101 325Pa;所有壁面采用無滑移邊界條件,對近壁面流動區域采用可伸縮壁面函數法處理,基于RNGk-ε湍流模型,近壁面邊界層網格節點數約為10,用以保證Y+在30~300之間.數值計算中引入了多重參考坐標系(MRF),葉輪流域設定為轉動域,其他流域設定為靜止域.交界面處理:定常計算轉子和靜子交界面采用混合面法[11],交界面上游(轉子)流動參數周向平均后傳遞給下游(靜子)交界面,相鄰的2個通道之間采用節點一一對應的周期性交界面.非定常計算采用滑移網格技術的轉/靜交界面法處理,由于計算模型動、靜葉片數量分別為15和16,為了節約計算資源采用了單通道處理,且通道比近似接近1∶1,采用此方法會使計算出的動靜干涉非定常脈動力略高于實際結果.葉輪轉動一周分成600個時間步,因此一個柵距通道內設定40個時間步,每個時間步長為3.424 6×10-5s,這個時間步長滿足動態壓力信號的采集需求.連續性方程、3個方向速度方程、湍動能方程和湍流耗散率方程的標準化均方根殘差均設置為10-5,且計算域進出口質量流量差值小于0.5%,作為流場計算的收斂標準.非定常計算以定常計算結果為初場,每個時間步下所有求解方程殘差降到10-5以下,葉輪轉動3 000 個時間步,即轉動5 周后,設定的監測點達到穩定的周期性波動,此時判定計算收斂.

1.3 噪聲計算模型及網格

基于URANS/FW-H 方程的混合計算聲學方法求解低速軸流風扇的遠場噪聲輻射.對于低馬赫數的單級風扇,級轉靜干涉單音噪聲是主要的單音噪聲[19].因此,聲源信息中與內部湍流流動相關聯的四極子聲源項和與葉片厚度相關的單級子聲源項被忽略,只保留由于周期性非穩定壓力脈動引起的偶極子聲源項.相對于時域信息而言,更關注頻域信息,FW-H 方程經過Fourier變換后得到頻域形式的Helmholtz波動方程,并采用聲學有限元方法離散該Helmholtz波動方程,聲學網格建模過程中考慮了葉片和機匣壁面對聲波的散射和反射作用,具體網格劃分見圖3.聲學網格采用非結構化的四面體網格,最大網格尺寸為18 mm,總的網格數為1 600 094,為保證在最大頻率下的空間求解精度,有限元網格單元長度L必須滿足:L≤c/(6fmax),c為聲速,fmax為最高計算頻率(3倍頻=2 190Hz).

圖3 軸流風扇聲學有限元網格Fig.3 Acoustic finite element mesh of the axial flow fan

2 結果分析與討論

2.1 轉靜干涉流場

任何葉輪機械的非定常流動行為是以非定常流動特征的出現為依據的,這種非定常流動特征來源于葉片表面的相對運動,而在葉片表面出現的邊界層產生的尾跡結構與主流摻混,并會隨著主流的傳播產生新的非定常流動現象.這些非定常流動特征會在流動中進行能量交換,此外,周期性非定常流動對下游葉片的干涉是此類葉輪機械最主要的單音噪聲源.因此,精準地描述非定常流動行為對于提高葉輪機械性能和降噪非常重要.

圖4給出了轉子出口50%葉展截面上各節點的速度隨周向位置的變化,圖中T表示轉子轉過一個通道所需的時間.轉子尾跡向下游傳播過程中有非常明顯的尾跡虧損,在此截面的軸向速度云圖(圖5)中發現一個狹長的貫穿整個通道的尾跡區,此尾跡與葉頂的間隙泄漏渦摻混,加劇了葉頂區域的流動損失.同時,這種尾跡結構會與主流摻混進入下游的靜葉通道,打斷靜葉前緣的一致性流動.

圖4 轉子出口50%葉展截面軸向速度隨周向位置的變化Fig.4 Circumferential distribution of axial velocity on 50% span section at rotor outlet

圖5 0.5T 時刻轉子出口截面軸向速度分布Fig.5 Axial velocity distribution on section of rotor outlet at 0.5T

圖6給出了葉片通道內50%葉展截面的熵增云圖.從圖6可以明顯地觀察到轉子出口的尾跡對靜葉前緣的周期性沖擊作用,這種周期性沖擊導致的靜葉前緣的非定常壓力波動正是此類葉輪機械最主要的單音噪聲來源.

圖7給出了靜葉通道50%葉展截面上熵增和軸向渦量矢量圖.圖7清晰顯示了轉子尾跡在靜葉通道中的演變過程:轉子尾跡周期性掃過靜葉,在靜葉前緣被靜葉切割,并在靜葉通道內拉伸變形以及與主流摻混等.在尾跡內部由于強烈的剪切作用產生嚴重的氣動損失,圖中尾跡區的高熵分布充分說明了這一點.從圖7還可以看出,尾跡在靜葉通道傳播過程中逐漸變寬,這是因為尾跡與靜葉通道中的主流摻混降低了尾跡區的速度虧損,增加了尾跡寬度;尾跡向下游輸送過程中通道內逆壓梯度增大,增加了尾跡寬度.

從圖7(a)可以看出,由于負射流作用[19],在轉子尾跡前后誘導出一對旋轉方向相反的漩渦結構.這些漩渦與尾跡區邊界層相互作用,改變了靜葉表面邊界層流動的發展,造成靜葉表面產生更劇烈的壓力波動.此外,由于周向逆壓梯度作用,尾跡在向下游輸送過程中使得漩渦向靜葉吸力面一側聚集,并在轉子尾跡掃過靜葉前,首先與靜葉前緣接觸而被靜葉切割,并與靜葉表面邊界層相互作用,這些漩渦使得邊界層內速度加快,加速了邊界層內流動的動量和能量交換,使得通道內的尾跡虧損降低,尾跡寬度增加.

為了比較動、靜葉表面壓力的大小,獲得準確的噪聲源強度和位置,圖8給出了動、靜葉表面均方根壓力分布.圖中均方根壓力的定義如下:

式中:為時均壓力;p為瞬時壓力;T為計算周期.

圖8 動、靜葉表面均方根壓力分布Fig.8 RMS pressure distribution on rotor and stator surface

為了表達清晰,圖8中給出了均方根壓力的對數表示形式.從圖8可以看出,靜葉前緣吸力面一側是壓力脈動擾動最大的區域,這里受到上游轉子尾跡和尾跡誘導漩渦的強烈沖擊作用,而此處也是最主要的噪聲源區域.同時,在靜葉吸力面還發現大范圍的壓力波動區域,這是因為轉子尾跡誘導出的漩渦由于周向逆壓梯度的作用向吸力面移動,此漩渦與靜葉吸力面表面的邊界層相互作用,破壞了邊界層內流動的一致性,增強了邊界層內的流動,從而導致吸力面邊界層較大范圍的波動.

2.2 偶極子源辨析

圖9 給出了靜葉在基頻(BPF)及其二倍頻(2BPF)下壓力波動的均方根值,左側為靜葉壓力面載荷分布,右側為吸力面載荷分布.從圖9 可以看出,最大載荷分布在靜葉前緣靠近吸力面區域,這是因為此處受上游轉子尾跡以及尾跡誘導漩渦的周期性強烈干涉作用,使得此區域靜葉表面的邊界層流動劇烈.此外,從圖9還可以看出,基頻的幅值遠高于二倍頻.綜上所述,基頻下靜葉前緣靠近吸力面區域是最主要的噪聲源區域.

圖9 靜葉非定常載荷分布Fig.9 Distribution of unsteady aerodynamic load on stator

2.3 噪聲預測結果

以葉片表面的氣動載荷作為聲源信息,并經過頻域的FFT 變換后加載給Ffowcs Williams-Hawking方程的頻域形式(即頻域聲學Helmholtz波動方程),來求解風扇的遠場噪聲輻射.計算遠場噪聲時,依據標準ISO 3744—1994《聲學聲壓法測定噪聲源聲功率級和聲音能量級 反射面上方近似自由場的工程法》建立包圍聲源區的場點網格.依據聲學測量要求在場點網格不同位置建立9個測點,各測點距離樣機對應實體表面1m,附近1m 內無發射面,依據式(3)得到測量表面平均聲壓級:

式中:N為測點數;Lpi為第i個測點的聲壓級;pe,i表示第i個測點的有效聲壓;基準聲壓取pref=20×10-6Pa.

分別計算動、靜葉表面偶極子激發的外場噪聲輻射,圖10分別給出了動、靜葉表面偶極子源噪聲輻射的平均聲壓級在不同頻率(BPF、2BPF 和3BPF)下的計算結果.由圖10可知,靜葉表面偶極子源激發的噪聲(95dB)遠大于動葉(77.5dB),且以基頻噪聲為主,這說明由靜葉表面非定常力激發的轉靜干涉噪聲是此類葉輪機械的主要單音噪聲.

圖11給出了基頻及其諧波軸截面聲壓級云圖,圖中上方為管道出口位置,下方為管道進口位置.從圖11可以看出,噪聲在管道壁面上以交替周向模態向管道進出口兩側傳播,且以管道出口為主.此外基頻的輻射聲壓級最大,最大聲壓級達到120dB,隨著頻率值衰減到2BPF,聲壓級幅值急劇減小,最大值減小到105dB,噪聲的管道周向傳播模態也有較大程度的衰減.當頻率值衰減到3BPF時,其聲壓級值相對于基頻聲壓幾乎可忽略不計,其管道周向傳播模態進一步衰減,但是變得更加密集且沒有發現明顯的聲場指向性.

圖10 基頻及其諧波的噪聲模態Fig.10 SPL acoustic mode at BPF and its harmonics

圖11 基頻及其諧波軸截面聲壓級云圖Fig.11 Sound pressure profile on axial section at BPF and its harmonics

以葉輪軸心位置為圓點,半徑1.5m 的軸截面上間隔10°均勻布置36個監測點,監測點以管道出口位置為起始零點.圖12給出了監測點的聲壓級分布.從圖12可以看出,聲壓級在軸截面上呈對稱分布,管道出口處的聲壓級輻射高,聲壓級隨著頻率的衰減而降低,即基頻下聲壓級最高,這與前面的聲壓級云圖分析一致.

3 結 論

(1)上游的轉子尾跡及葉頂間隙泄漏渦對下游靜葉前緣的周期性干涉破壞了靜葉前緣的一致性流動.這種周期性非定常干涉改變了靜葉表面前緣的壓力波動.尾跡與葉頂的間隙泄漏渦摻混加劇了葉頂的流動損失,使靜葉前緣葉頂區域產生較大的壓力波動.

圖12 聲場指向性分布Fig.12 Directivity distribution of far-field SPLs

(2)由于負射流效應,在轉子尾跡前后誘導出一對旋轉方向相反的漩渦,這種漩渦結構在和尾跡一起向下游輸送過程中首先與靜葉前緣接觸,并由于靜葉壓差作用向靜葉吸力面聚集,漩渦與靜葉吸力面表面邊界層相互作用,加速了邊界層內的動量和能量交換,造成了靜葉吸力面表面更大的壓力波動.

(3)噪聲源的頻譜分析表明,在基頻下,靜葉前緣靠近吸力面區域是最主要的噪聲源區域.

(4)平均聲壓級對比分析表明,靜葉表面非定常力(偶極子源)激發的轉靜干涉噪聲是此類葉輪機械的主要單音噪聲.噪聲主要以基頻為主,聲壓級隨著頻率的衰減而降低.此類轉靜干涉噪聲主要從風扇管道出口向外輻射傳播.

[1]WRIGHT S E.The acoustic spectrum of axial flow machines[J].Journal of Sound and Vibration,1976,45(2):165-223.

[2]CUMPSTY N A.A critical review of turbomachinery noise[J].Journal of Fluids Engineering,1977,99(2):278-293.

[3]BLAKE W K.Mechanics of flow-induced sound and vibration[M].Orlando,Fla,USA:Academic Press,1986.

[4]FLEETER S.Discrete frequency noise reduction modeling for application to fanjet engines[J].The Journal of the Acoustical Society of America,1980,68(3):957-965.

[5]WOODWARD R P,ELLIOTT D M,HUGHES CE,etal.Benefits of swept and leaned stators for fan noise reduction[J].Journal of Aircraft,2001,38(6):1130-1138.

[6]YOUNSI M,BAKIR F,KOUIDRI S,etal.Numerical and experimental study of unsteady flow in a centrifugal fan[J].Proceeding of the Institution of Mechanical Engineers,Part A:Journal of Power and Energy,2007,221(7):1025-1036.

[7]KHELLADI S,KOUIDRI S,BAKIR F,etal.Predicting tonal noise from a high rotational speed centrifugal fan[J].Journal of Sound and Vibration,2008,313(1/2):113-133.

[8]DíAZ A K M,FERNáNDEZ O J M,MARIGOTA E B,etal.Numerical prediction of tonal noise generation in an inlet vaned low-speed axial fan using a hybrid aeroacoustic approach[J].Proc IMechE,Part C:J Mechanical Engineering Science,2009,223(9):2081-2098.

[9]萬劍鋒,楊愛玲,戴韌.風機葉片表面分離渦與寬頻噪聲輻射特性的分析[J].動力工程學報,2014,34(9):736-741.

WAN Jianfeng,YANG Ailing,DAI Ren.Analysis of separation vortex and broadband noise radiation of fan blade surfaces[J].Journal of Chinese Society of Power Engineering,2014,34(9):736-741.

[10]王風華,楊愛玲,戴韌,等.前彎葉片對軸流泵噪聲輻射性能的影響[J].動力工程學報,2011,31(6):440-444.

WANG Fenghua,YANG Ailing,DAI Ren,etal.Effect of blade's forward-skewed angles on flow sound radiation of axial pumps[J].Journal of Chinese Society of Power Engineering,2011,31(6):440-444.

[11]LEE D J,JEON W H,CHUNG K H.Development and application of fan noise prediction method to axial and centrifugal fan[C]//Proceeding of ASME FEDSM'02 Fluids Engineering Division Summer Meeting,Montreal.Quebec,Canada:ASME,2002:987-992.

[12]POLACSEK C,DESBOIS-LAVERGNE F.Fan interaction noise reduction using a wake generator:experiments and computational aeroacoustics[J].Journal of Sound and Vibration,2003,265(4):725-743.

[13]CAI Jiancheng,QI Datong,LU Fuan,etal.Study of the tonal noise of a centrifugal fan at the blade passing frequency.PartⅠ:aeroacoustics[J].Journal of Low Frequency Noise,Vibration and Active Control,2010,29(4):253-266.

[14]HU Binbin,OUYANG Hua,WU Yadong,etal.Numerical prediction of the interaction noise radiated from an axial fan[J].Applied Acoustic,2013,74(4):544-552.

[15]KOOPMANN G H,NEISE W,CUNEFARE K A.Fan casing noise radiation[J].Journal of Vibration and Acoustics,1991,113(1):37-42.

[16]LIU Q,QI D,MAO Y.Numerical simulation of centrifugal fan noise[J].Proceedings of Institute Mechanical Engineers,Part C:Journal of Mechanical Engineering Science,2006,220(8):1167-1178.

[17]TAN Hongtao,QI Datong,LIU Fuan.Identification of acoustic dipole source in a propylene centrifugal compressor stage [C]//Proceeding of ASME Turbo Expo 2010:Power for Land,Sea and Air.Glasgow,UK:ASME,2010:14-18.

[18]YAKHOT V,ORSZAG S A,THANGAM S,etal.Development of turbulence models for shear flows by a double expansion technique[J].Physics of Fluids A,1992,4(7):1510-1520.

[19]WECKMULLER C,GUERIN S,ASHCROFT G.CFD-CAA coupling applied to DLR UHBR-fan:comparison to experimental data[C]//15thAIAA/CEAS Aeroacoustics Conference (30th AIAA Aeroacoustics Conference).Miami,Florida:AIAA,2009:1-11.

主站蜘蛛池模板: 亚洲视频四区| 成年午夜精品久久精品| 激情乱人伦| 一级不卡毛片| 亚洲国产精品日韩专区AV| 久久99国产精品成人欧美| 夜夜操天天摸| 一级毛片在线播放| 男人天堂伊人网| 国产成人精品一区二区秒拍1o| 亚洲高清国产拍精品26u| 日本免费a视频| 国产爽妇精品| 亚洲码一区二区三区| 国产视频入口| 国产一区二区三区夜色| 国产h视频免费观看| 成人小视频网| 超清无码一区二区三区| 欧美国产菊爆免费观看| 毛片基地视频| 亚洲国产综合精品一区| 亚洲欧洲天堂色AV| 亚洲人成人伊人成综合网无码| 日韩中文无码av超清| 91精品人妻互换| 日韩精品一区二区三区swag| 国产99精品久久| 狠狠做深爱婷婷久久一区| 亚洲中文制服丝袜欧美精品| 中日韩欧亚无码视频| 这里只有精品免费视频| 国产精品.com| 欧美区一区二区三| 欧美一级大片在线观看| 亚洲国产一区在线观看| 无码内射中文字幕岛国片| 青青热久免费精品视频6| 国产日韩av在线播放| 在线观看免费国产| 成人在线天堂| 国产三级精品三级在线观看| 狠狠干综合| 久久精品66| 久久国产精品麻豆系列| 亚洲欧美一区二区三区图片| 欧美日韩v| 久久久波多野结衣av一区二区| 成人在线综合| 亚洲欧洲日产国码无码av喷潮| 亚洲高清在线播放| 国产自产视频一区二区三区| 韩日无码在线不卡| 国产欧美性爱网| 亚洲三级a| 欧美日韩综合网| 人妻91无码色偷偷色噜噜噜| 在线观看亚洲成人| 一级毛片在线播放| 国产精品久久久久久搜索| 免费无码又爽又黄又刺激网站| 国产本道久久一区二区三区| 国产 在线视频无码| 青青草欧美| 欧美一级黄片一区2区| 久久这里只有精品2| 免费一级无码在线网站| 四虎影视8848永久精品| аⅴ资源中文在线天堂| 国产成a人片在线播放| 综合天天色| 久久熟女AV| 福利在线不卡| 欧美伦理一区| 国内精品免费| 国产极品美女在线观看| 2020最新国产精品视频| 狠狠色狠狠色综合久久第一次| 夜夜拍夜夜爽| 亚洲a级在线观看| 日韩第一页在线| 日韩精品亚洲一区中文字幕|