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

空化對軸流式水輪機尾水管壓力脈動和轉輪振動的影響

2021-09-04 12:01:24朱國俊馮建軍羅興锜
農業工程學報 2021年11期
關鍵詞:振動信號

朱國俊,李 康,馮建軍,2※,羅興锜,2

(1.西安理工大學水利水電學院,西安 710048;2.浙江富安水力機械研究所,杭州 311121)

0 引 言

目前,經濟指標優良、可開發性好的中高水頭水能資源已基本開發完畢,低水頭水能資源是未來水電開發的重點方向之一。軸流式水輪機是進行低水頭水能資源開發的兩大主力機型之一,其性能的優劣直接影響著水能資源的轉換效率[1-2]。軸流式水輪機的主要應用限制在于空化造成的轉輪空蝕破壞及水輪機穩定性劣變。振動和壓力脈動是反映水輪機穩定性的常用指標,空化通過惡化壓力脈動和振動來影響軸流式水輪機的穩定性,所以開展空化對軸流式水輪機壓力脈動和轉輪振動的影響規律研究對軸流式水輪機穩定性的優化有重要意義。

水輪機壓力脈動及振動是一種非穩態過程,其產生原因十分復雜,主要與水力激勵以及系統響應有關[3-4]。目前,關于常規、無空化條件下的水輪機壓力脈動主要通過數值模擬[5-9]和模型試驗方法[10-15]來進行研究,這些相關研究的重點均著眼于分析水輪機尾水管的壓力脈動特性及揭示其產生機理。近年來,隨著水輪機的安全穩定性要求的不斷提升,空化對水輪機壓力脈動或者振動特性的研究也逐漸涌現。Rus等[16-17]對軸流式水輪機進行了模型試驗,采集了不同空化工況下的振動信號數據,結果發現振動幅值隨空化系數的降低而上升。唐巍等[18]采用數值模擬方法對中低比轉速混流式水輪機空化工況下的水力振動進行了分析,發現機組的振動頻率為轉頻的整數倍。李琪飛等[19]結合試驗與數值模擬方法研究了空化條件下水泵水輪機的壓力脈動特性,發現水泵水輪機無葉區的壓力脈動主頻為轉輪轉頻,且幅值隨空化系數的降低而增加。徐洪泉等[20]采取將理論分析與模型試驗相結合的方法,發現尾水管高倍轉速頻率壓力脈動由轉輪葉片或導葉端面間隙空化所引起。但大部分相關研究都單獨分析空化對壓力脈動特性的影響或空化對水輪機尾水管部位振動特性的影響,同時探究空化對水輪機壓力脈動及轉輪振動特性影響規律,并分析兩者關聯關系的研究很少。轉輪是軸流式水輪機最重要的旋轉部件,既是引發軸流式水輪機內部渦流、壓力脈動的主因,也是流道內水力激振力在轉子系統上的唯一作用部位,同時還是軸流式水輪機內空化作用首先出現并影響最嚴重的部件,因此,有必要開展研究揭示空化對轉輪振動和水輪機尾水管壓力脈動的影響規律,以便為軸流式水輪機穩定性的改善提供技術支撐。

本文以軸流式模型水輪機為研究對象,通過構建包含高速攝影系統、激光測振儀(Laser Doppler Vibrometer,LDV)和高頻壓力脈動傳感器的同步測試系統采集了變空化系數下的轉輪徑向振動信號、尾水管壓力脈動信號和轉輪空化圖像資料,綜合變分模態分解算法(Variational Mode Decomposition,VMD)和去趨勢互相關分析法對采集到的振動和壓力脈動數據進行了主成分提取,結合轉輪空化圖像定量分析了空化對軸流式水輪機尾水管壓力脈動和轉輪振動的影響規律,以期為水電站中軸流式水輪機設備運行穩定性的改善提供理論依據。

1 材料與方法

1.1 試驗裝置和測試方法

構建包含高速攝像組件、LDV及高頻壓力脈動傳感器的同步測試系統,并采用該同步測試系統對某 4葉片軸流式水輪機模型進行測試。該軸流式水輪機模型對應的真機額定水頭為14.83 m,參數如表1所示。

表1 水輪機模型參數Table 1 Parameters of turbine model

本文中的軸流式水輪機模型試驗于浙江富安水力機械研究所滿足IEC60193測試要求的高精度水力機械試驗臺上完成。根據水輪機試驗臺流量計的誤差為±0.150%FS,差壓變送器測量水頭的誤差為±0.080%FS,測功機測力矩的誤差為±0.075%FS,轉速傳感器的誤差為±0.003‰FS,可得效率測試的系統誤差為±0.191%。由系統誤差與效率測試的隨機誤差±0.035%可得水輪機模型效率最終測試的綜合誤差為±0.194%。水輪機模型試驗系統圖如圖1所示。

本文構建的同步測試系統包含高速攝像組件,該組件由頻閃儀、高速電荷耦合元件(Charge Coupled Device,CCD)攝像機和相應的采集軟件組成,可以記錄不同空化數下轉輪內部的空化圖像。同步測試系統還包含LDV、高頻壓力脈動傳感器和數據采集卡。本系統中采用的LDV為德國Polytec公司的VGO-200高精度數字便攜式LDV,根據檢定證書可知其 95%置信區間內測量不確定度為1%,測量分辨率可達0.02μm/(s?Hz0.5)。因為軸流式水輪機模型的轉輪室采用高透明有機玻璃制造,所以LDV的激光束可穿透有機玻璃聚焦在轉輪體上采集轉輪的徑向振動速度。測試系統中高頻壓力脈動傳感器(M112A22,美國PCB公司)的檢定證書表明其95%置信區間內測量不確定度為 1%。同步測試系統采集的轉輪徑向振動信號和壓力脈動信號都通過屏蔽電纜接入美國NI公司的便攜式數據采集模塊中。同步測試系統儀器的現場布置如圖2所示。

尾水管錐管處是水輪機壓力脈動測試時的重點部位。因此,本文的壓力脈動測點布置在軸流式水輪機尾水管錐管段,位于軸流式轉輪旋轉中心以下 0.2D1(D1為轉輪標稱直徑,cm)處。壓力脈動測點數目為2個,二者之間的周向間隔為180°,測點位置如圖3所示。

水輪機額定工況是水輪機性能的重點考核工況,也是軸流式水輪機必須實現無空化運行的工況,為研究空化對軸流式水輪機壓力脈動及轉輪振動的影響,試驗過程中選擇水輪機額定工況作為基礎工況。通過改變試驗裝置的空化系數σ實現了軸流式水輪機即轉輪由無空化到完全空化的變化過程,即轉輪由無空化到完全空化的過程。該過程中各個工況下σ的數值如表2所示,其數值可根據同參數轉輪的試驗經驗來調整下游低壓水箱內的絕對壓力進行設置,確保能準確捕捉空化初生。在表2中的各個工況下均采集振動和壓力脈動數據,為了準確分析空化對測試數據的影響,傳感器的采樣頻率設置為25.6 kHz。

表2 空化試驗過程中空化系數σ的數值Table 2 Value of cavitation coefficient σ in the process of cavitation experiment

通過前述的試驗臺和同步測試系統獲得轉輪的徑向振動和壓力脈動數據以后,對采集到的信號數據進行降噪和主分量提取,以便進行進一步的分析。

1.2 變分模態分解算法理論

獲取非平穩時序信號不同成分的分解方法較多,如小波分解、經驗模態分解(Empirical Mode Decomposition,EMD)以及EMD的改進型算法等。其中,EMD及其改進型算法是使用較多的分解算法,但它們仍在一定程度上存在模態混疊、邊界效應等問題。為解決上述問題,本文采用變分模態分解算法(VMD)對采集到的信號數據進行分解。VMD算法將信號數據的分解約束到變分框架內,通過構造并求解約束變分問題實現原始信號的分解[21-22]。這種方法的優勢在于其采用了完全非遞歸的處理策略,相比于EMD類算法的遞歸模式分解策略,能有效抑制或完全避免模態混疊、邊界效應等問題。

VMD算法的目的是將復雜的信號分解為K個模態分量,并且要求每個模態分量都有各自的中心頻率和帶寬,且使得各個模態分量帶寬之和最小,其約束表達式如下:

式中K為信號分解得到的模態分量數量;uk(t)為分解信號所得的第k個模態分量;f(t)為原始信號;j為旋轉算子;δ(t)為狄拉克函數;ωk為第k個模態分量的中心頻率,Hz;t為時間,s。

引入二次懲罰因子α和Lagrange乘數算子λ,將約束變分問題轉化為非約束變分問題,則可得增廣 Lagrange函數L(uk,ωk,λ)如下:

對轉化后的無約束變分方程采用交替方向乘子法[23]來求解,通過不斷迭代求取到擴展Lagrange方程的鞍點,收斂條件如式(3)所示。

式中n代表迭代序號,表示第n次迭代;分別為的傅里葉變換;ω為頻率變量,Hz;ε為一極小值,取為10-7可滿足本文計算精度要求。

如果式(3)滿足,則輸出K個分解所得的模態分量。如果式(3)不滿足,則利用傅里葉變換,將迭代后的結果按式(4)~式(6)轉換到頻域上進行更新,然后再重新根據式(3)進行判斷,直到式(3)的條件被滿足,輸出K個分解所得的模態分量。

1.3 去趨勢互相關分析

通過VMD算法將信號分解為K個分量以后,選取與原始信號相關性最高的分量作為原始信號的主成分。在信號分析領域,互相關[24-25]通常用來度量兩組信號之間的相似特征,分析信號間的相關程度。對于任意2組非平穩時間序列信號{xt}和{yt},則去趨勢互相關系數[26-28]ρDCCA表示為{xt}和{yt}的去趨勢協方差與去趨勢方差間的比值,公式如下:

式中F2DCCA表示{xt}和{yt}的去趨勢協方差;FDFA{xt}和FDFA{yt}分別表示{xt}和{yt}的去趨勢方差。若ρDCCA=0,表示2個信號序列間無相關性;若0<ρDCCA<1,表示2個信號序列呈正相關,且ρDCCA=1時,說明2個信號序列呈高度嚴格正相關;若?1<ρDCCA<0,表示 2個信號序列呈負相關,且ρDCCA=?1時,說明2個信號序列呈高度嚴格負相關。具體的算法流程見圖4。

1.4 壓力脈動無量綱幅值計算

在空化試驗過程中,通過壓力脈動傳感器采集了表2中不同空化系數下的水輪機尾水管壓力脈動信號。為了便于數據分析,本文對壓力脈動的幅值進行了無量綱化處理,具體計算表達式[29]如下:

式中Cp表示無量綱壓力脈動幅值;Δp表示實際壓力脈動值,Pa;ρ表示流體密度,kg/m3;H表示模型試驗中水輪機進、出口的水頭差,m。

2 結果與分析

2.1 壓力脈動與徑向振動的時域特性

2.1.1 壓力脈動數據的時域分析

通過高速攝影圖像分析,轉輪的空化過程可分為無空化、空化初生、空化發展和完全空化 4個階段。限于篇幅,在圖5中針對轉輪空化的4個階段分別給出對應的空化圖像,圖6中則給出了代表性σ下的壓力脈動時域信號波形圖(單個傳感器)。

從圖5中可以看出,在空化初生階段(σ=1.50),轉輪葉片上稀疏的空泡最先出現在轉輪葉片輪緣翼型的中后部。隨著空化系數σ的降低,在空化發展階段(σ=0.93),該位置處的空泡逐漸變多,并發展成包含間隙空化和云空化的混合空化現象,同時輪轂處也出現了明顯的云狀空化現象。當空化系數繼續降低、空化繼續發展到了完全空化階段(σ=0.63),轉輪葉片輪緣翼型的下游已由大片云狀空化和大尺度空泡混合而成的復雜空化現象占據,此外,輪緣翼型在距離其前緣約 20%翼型弦長位置處也出現了明顯的由間隙泄漏渦導致的渦空化現象。轉輪空化程度和形態的急劇變化對壓力脈動幅值的影響也由圖6反映出來。根據圖6可知,隨著空化程度的不斷加劇,壓力脈動的幅值明顯增大,特別是完全空化階段(σ=0.63),此時壓力脈動的幅值已遠高于空化發展的其他階段。

2.1.2 轉輪徑向振動數據的時域分析

不同空化階段下軸流式轉輪的徑向振動速度信號如圖7所示。從圖中可以看出,隨著空化系數的降低,振動速度信號的幅值隨之增大,當空化系數σ達到0.63時,振動速度幅值的均方根值達到最大值1.36 mm/s。將圖6與圖7進行比較可以發現,壓力脈動與振動速度信號的幅值變化趨勢基本一致。將壓力脈動峰峰值ΔCp與轉輪徑向振動速度的峰峰值隨σ的變化曲線進行比較,如圖8所示。

從圖8中可以看出,壓力脈動峰峰值與轉輪徑向振動速度峰峰值隨σ的變化曲線非常相似。在軸流式轉輪由無空化到空化初生的階段,壓力脈動峰峰值與轉輪徑向振動速度峰峰值的增加幅度相對較小。空化初生時的壓力脈動峰峰值是無空化時的1.49倍,轉輪徑向振動速度峰峰值是無空化時的2.32倍。當轉輪空化進一步發展達到空化發展階段,壓力脈動峰峰值和轉輪徑向振動速度峰峰值均出現了較為明顯的增加,分別提高至無空化時的2.58和3.36倍。當轉輪空化進入到充分發展階段直至完全空化時,壓力脈動峰峰值和轉輪徑向振動速度峰峰值陡增,分別達到了無空化時的9.16和10.12倍。為了探討空化后壓力脈動與轉輪徑向振動速度峰峰值的增長率隨σ的變化規律,以無空化時的壓力脈動與轉輪徑向振動速度峰峰值為基準,在圖9中給出了不同σ工況下兩類峰峰值的增長率。

從圖9中可以看出,壓力脈動峰峰值和轉輪徑向振動速度峰峰值的增長率隨σ的變化規律均呈現明顯的非線性。2類峰峰值的增長率均在空化發展前呈平緩增長趨勢,然后在空化發展至完全空化階段急劇增加直到最大值。此外,通過圖9可知,從無空化到空化初生時,轉輪徑向振動速度和壓力脈動峰峰值均出現了明顯的增加,但轉輪徑向振動速度峰峰值的增長率達到了132.2%,遠高于壓力脈動峰峰值的增長率(僅為49.2%),表明空化對轉輪振動的影響程度遠高于壓力脈動。因為轉輪是空化產生和直接作用的部件,因此空泡聚合、潰滅所產生的沖擊波首先作用于轉輪,然后再通過水體傳播到水輪機其他部位。而由于沖擊波能量在水體中的不可逆耗散,當其傳播到距空化區域一定距離的壓力脈動測點時,已無法引發與空化核心區等效的壓力波動[30],所以壓力脈動的峰峰值增長率小于轉輪徑向振動速度峰峰值的增長率。

2.2 壓力脈動與徑向振動的頻域特性

基于VMD算法和去趨勢互相關分析技術,本文將測得的壓力脈動信號和轉輪徑向振動速度信號進行降噪并提取主成分,然后再進行頻譜分析。VMD算法中的K值經對比篩選后取 7。限于篇幅,下面只以完全空化工況(σ=0.63)下的壓力脈動信號為例展示主成分提取過程的中間結果。圖10為分解該工況下壓力脈動信號得到的固有模態函數(Intrinsic Mode Function,IMF)分量圖。由圖10可知,信號分解后得到的各固有模態分量的振幅差異表明各固有模態分量的能量也存在區別。為了獲取信號中的主要成分,采用本文第1.3節中的去趨勢互相關分析技術計算各IMF分量與原始信號的互相關系數ρDCCA,ρDCCA值最高的 IMF分量即被選擇作為主成分。在圖10中,經計算得到ρDCCA值最大的分量為IMF1,所以選擇IMF1作為該信號的主分量。采用同樣方法獲得圖6中不同σ下的壓力脈動信號主成分時域波形圖,并分別進行快速傅里葉變換,最終得到不同σ下的壓力脈動信號主成分頻域圖如圖11所示。

從圖11中可知,在無空化階段,壓力脈動的主頻是明顯的低頻,為0.2fn(fn為轉輪轉頻,Hz)。根據國際電工委員會的水輪機試驗驗收標準IEC60193對壓力脈動誘因的總結,導致尾水管錐管上該低頻壓力脈動現象的原因可能為轉輪出口的渦旋流動。根據圖5的圖像資料可知,當空化初生時,空化流動為零散、稀疏的小尺度空泡組成的空泡流,并且只集中在葉片輪緣中后部的極小區域,因此無法對主流產生影響,所以壓力脈動信號的低頻區域相比無空化時基本沒有變化。結合陳廣豪[31]研究中得到的空泡非均勻潰滅會產生高頻壓力脈沖向四周傳播的結論,在空化初生階段,壓力脈動信號頻域圖中頻率為24.1fn的主頻成分為葉片輪緣翼型中后部附著型空泡群非均勻潰滅所導致。當空化進一步發展,進入空化發展階段時,顯著的間隙渦空化現象已對下游流動產生了明顯擾動,壓力脈動信號的主頻等于轉輪葉片通頻,且空化導致的壓力脈動高頻部分幅值也明顯增加。隨著σ進一步降低達到完全空化階段時,通過高速攝影圖像可知大尺度的空化區域已接近壓力脈動的測點位置。因此,已完全空化的葉片周期性掃過測點所在的周向位置時引發了強烈的壓力脈動,頻率為葉片通頻,幅值則急劇增加至空化發展階段(σ=0.93)主頻幅值的14.3倍。上述現象表明,當軸流式轉輪進入空化發展階段以后,會影響下游壓力脈動的頻域。

采用本節方法對轉輪的徑向振動速度信號進行主成分提取,并對主成分信號進行快速傅里葉變換,得到圖12所示的頻域圖。

從圖12中可以發現,除了σ=0.93的工況外,轉輪的徑向振動的主頻基本為轉輪葉片通頻(f=4.0fn)。在σ為2.08、1.50及0.93工況下,σ的下降主要引發的是轉輪徑向振動信號高頻段能量的變化,對葉片通頻成分的幅值影響很小。

在σ=0.93工況,轉輪葉片的徑向振動信號中頻率為36倍轉頻的高頻成分幅值增加到最大值0.180 mm/s,并超過了提升幅度很小的葉片通頻對應的幅值0.164 mm/s,所以此刻振動信號主頻為 36fn的高頻成分。原因是此時轉輪內部的空化處于發展的臨界狀態,由圖5可知該工況下非定常云狀空化區域已形成一定的規模,云狀空化區域密集的小尺度空泡潰滅給葉片造成了高頻沖擊,所以此時轉輪徑向振動高頻成分的幅值出現了增長。而葉片通頻成分的幅值基本不變,此消彼長下,振動信號的高頻成分成為主頻。即便如此,該工況下轉輪葉片通頻對應的幅值0.164 mm/s也僅略低于主頻的0.180 mm/s。

誘發轉輪葉片通頻頻率振動的原因為轉輪體安裝 4個轉輪葉片后無法達到完全的周向質量均勻,即周向質量不均衡所導致。這種周向質量不均衡在轉輪完全空化后愈發劇烈,所以σ=0.63工況下葉片通頻對應的幅值遠高于其他工況。此外,圖12表明空化導致的中高頻幅值提升基本只局限于12.0fn~200.0fn的范圍內。

2.3 壓力脈動與徑向振動的能量分布

為了詳細分析原始壓力脈動和轉輪徑向振動速度信號在整個頻段的能量,通過巴塞伐爾(Parseval)定理計算獲得了2種信號在不同σ下的能量譜。由于不同頻率成分的信號能量差異大,所以縱坐標采用對數坐標。此外,為了便于表示壓力脈動信號的能量,計算能量譜時壓力脈動的幅值不做無量綱化處理,采用實測脈動壓力幅值進行計算。壓力脈動和轉輪徑向振動速度測試數據的能量譜如圖13所示。

綜觀圖13可知,不管是壓力脈動還是轉輪的徑向振動,隨著σ的下降,高頻段(f/fn>50)能量提高。此外,在壓力脈動信號能量譜圖(圖13a)中可以發現,隨著σ的降低,高頻段的局部能量極值發生了朝低頻區域遷移的現象,從而導致空化以后的壓力脈動低頻區域能量明顯提高。聯合圖5中不同σ下的轉輪空化形態圖像以及陳廣豪[31]揭示的空化形態演變與壁面壓力脈動頻幅特性間關聯關系可知,σ的下降使得轉輪上的空化形態由小尺度稀疏泡狀空化(空化初生)向云狀空化與大尺度空泡混合的復雜形態(完全空化)進行演變,空化形態的演變不斷導致壁面壓力脈動的低頻能量成分增強,進而使不同σ下的壓力脈動信號在能量譜圖上表現出高頻段的局部能量極值向低頻區域遷移的現象。大型水輪發電機組的共振頻率都是低頻,尾水管壓力脈動的這種高頻能量局部極值遷移現象提高了誘發機組共振的可能性,增加了機組運行的不確定性。而根據圖13b可知,轉輪的徑向振動速度信號則沒有這種現象。

由圖13b可以發現,空化程度的增加致使振動速度信號整個頻段的能量都出現了提升,并且高頻段(f/fn>50)的能量分布形態基本類似。因為空化直接作用于轉輪葉片及轉輪體,因此空化程度的增加直接導致了空化附著物高頻振動能量的增加,與尾水管壓力脈動這種受空化間接影響的現象明顯不同。

目前,上述研究中發現的規律只局限于軸流式轉輪葉片及轉輪體。因為空化引發的振動特性取決于空化發展過程中空化的形態,而水流在繞流不同對象時產生的空化形態各異,引發的振動特性變化規律也各有不同。當前無法確保本研究中發現的空化誘導振動能量變化的規律可普適于所有研究對象,還有待其他學者共同開展進一步的研究。

3 結 論

本文采集了不同空化工況下的軸流式水輪機尾水管壓力脈動和轉輪徑向振動數據,并對測試數據進行了分析,主要結論如下:

1)轉輪空化發生后,隨著空化系數的下降,尾水管錐管上的壓力脈動與轉輪徑向振動速度峰峰值的增長率呈明顯的非線性變化規律。空化充分發展以后,錐管上的壓力脈動與轉輪徑向振動速度峰峰值會出現陡增,分別達到了無空化時的9.16和10.12倍。

2)轉輪空化程度的增加使得尾水管錐管上的壓力脈動主頻發展為葉片通頻,但對轉輪徑向振動速度的主頻影響較小。空化程度的增加會致使轉輪徑向振動信號中12~200倍轉頻范圍內的振動幅值明顯增加。

3)空化系數的降低會導致尾水管錐管壓力脈動產生高頻能量局部極值遷移現象,從而提升壓力脈動低頻區域的能量,增加引發機組共振的可能性,而轉輪的徑向振動速度則無此特性。

猜你喜歡
振動信號
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
噴水推進高速艇尾部振動響應分析
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
This “Singing Highway”plays music
孩子停止長個的信號
振動攪拌 震動創新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動性
基于LabVIEW的力加載信號采集與PID控制
一種基于極大似然估計的信號盲抽取算法
主站蜘蛛池模板: 欧美日韩国产在线播放| 制服丝袜在线视频香蕉| 欧美黄色a| 久精品色妇丰满人妻| 好紧太爽了视频免费无码| 69av在线| 精品在线免费播放| 在线综合亚洲欧美网站| 91小视频版在线观看www| 国产美女主播一级成人毛片| 精品国产一区二区三区在线观看 | 国产成人无码综合亚洲日韩不卡| 日韩无码黄色| 强奷白丝美女在线观看| 97亚洲色综久久精品| 美女视频黄频a免费高清不卡| 国产SUV精品一区二区| 午夜啪啪福利| 高清色本在线www| 亚洲一区二区三区在线视频| 亚洲欧洲日本在线| 欧美 国产 人人视频| 免费国产小视频在线观看| 国产激爽爽爽大片在线观看| 伊人五月丁香综合AⅤ| 国产精品永久不卡免费视频| 亚洲国产91人成在线| 成色7777精品在线| 国产91无码福利在线| 亚洲二区视频| 中文一区二区视频| 国产成年无码AⅤ片在线| 久久这里只精品热免费99| 亚洲码一区二区三区| 波多野结衣国产精品| 国产女同自拍视频| 国产精品美女自慰喷水| 精品中文字幕一区在线| a级免费视频| 黄网站欧美内射| 噜噜噜综合亚洲| 女人18毛片一级毛片在线| 久久精品嫩草研究院| 国产香蕉97碰碰视频VA碰碰看| 亚洲成a人在线观看| 精品久久久久久成人AV| 91青青在线视频| 亚洲系列中文字幕一区二区| 1级黄色毛片| 伊伊人成亚洲综合人网7777| 日韩无码白| 免费看的一级毛片| 国产1区2区在线观看| 天天躁夜夜躁狠狠躁躁88| 一级毛片在线直接观看| 麻豆精品在线播放| 国产91透明丝袜美腿在线| 99视频在线看| 91国内外精品自在线播放| 国产av无码日韩av无码网站| 亚洲国产精品国自产拍A| 在线无码九区| 欧美一级大片在线观看| 亚洲欧美另类专区| 香蕉精品在线| 国产女人18毛片水真多1| 亚洲精品无码AV电影在线播放| 在线观看无码a∨| 国产裸舞福利在线视频合集| 国产一级裸网站| 蜜桃视频一区二区三区| 久久国产乱子| 国产精品人成在线播放| 2021国产精品自产拍在线观看| 亚洲妓女综合网995久久| 亚洲水蜜桃久久综合网站| 亚卅精品无码久久毛片乌克兰| 青青操视频免费观看| 成人福利视频网| 97av视频在线观看| 国产在线第二页| 亚洲视频色图|