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

偏心旋轉閥氣動噪聲產生及控制機理研究

2022-12-20 15:44:32郝嬌山劉柏圻楊恒虎王偉波
噪聲與振動控制 2022年6期
關鍵詞:閥門模型

廖 靜,郝嬌山,劉柏圻,楊恒虎,王偉波

(1.重慶川儀自動化股份有限公司 技術中心調節閥研究所,重慶 400707;2.重慶川儀調節閥有限公司,重慶 400707)

偏心旋轉閥作為煤氣化裝置中的關鍵調節閥,在氣化爐升壓、泄壓、保壓運行、停車放空及氮氣置換的工藝操作環境中,承擔著背壓調節、控制系統壓力、安全放空等重要作用[1]。但其在功能實現過程中可能產生的噪聲不僅對環境產生污染,還會影響閥門本身性能,同時,伴隨噪聲產生的結構振動將造成相鄰管路和設備疲勞損傷,甚至可能導致安全事故,該閥門的高噪聲問題在煤氣化領域日益受到關注[2]。因此,研究偏心旋轉閥氣動噪聲的產生機理,是實現噪聲控制的重要基礎。

近年來,隨著計算機能力的大幅提高以及計算流體力學和工程聲學軟件開發日益成熟,國內外學者在閥門氣動噪聲問題上進行了大量數值模擬研究。氣動噪聲的數值模擬方法主要有氣動聲學方法(簡稱CAA)、萊特希爾聲類比方法(Lighthill 聲類比)和混合計算方法(Hybrid Method)[3]。Wei 等[4]采用基于寬帶噪聲源模型的直接模擬方法對高參數減壓閥的內部流場和聲場進行了數值分析,確定了射流產生的渦流和激波的交替形成和脫落是氣動噪聲產生的主要原因。孟令雅等[5]采用直接模擬方法和混合模擬方法計算氣流流過球閥產生的氣動噪聲,結果表明閥門的噪聲源為偶極子聲源和四極子聲源。徐號鐘等[6]采用基于Lighthill 聲類比模型的CFD/CAA 多步驟混合計算方法對截止閥腔內氣動聲學共振特性進行數值計算,得到了閥門腔內聲共振特性對噪聲的影響,以及小孔消聲裝置對閥門內流動聲源的控制效果。孫長周等[7]通過Lighthill 聲類比方法計算得到某型號流量調節閥內壁面處的聲壓,通過聲振模型獲得閥門外部聲場以及外部監測點的聲壓級頻譜曲線。綜上所述,閥門氣動噪聲的數值模擬大多采用基于Lighthill聲類比模型的混合計算方法,該模型對于低馬赫數氣動噪聲的計算效果較好,但對于偏心旋轉閥的高馬赫數氣動噪聲的數值模擬需采用考慮了聲源區域流體的可壓縮性和聲輻射區域流體的對流效應的M?hring 聲類比模型。該模型已應用于在發動機和壓氣機的噪聲預測[8-9],但有關其在閥門領域中的應用文獻較少。

鑒于此,本文選取煤氣化裝置中用于合成氣洗滌塔高壓差壓力調節的偏心旋轉閥為研究對象,按照實際工況參數建立偏心旋轉閥及附屬管道的有限元分析模型,應用ANSYS Fluent 2020 R2 流體仿真軟件和MSC Actran 15.0 聲學仿真軟件開展聲流固耦合數值模擬計算。通過Conval 11.0閥門選型軟件求解氣動噪聲標準值,驗證數值模擬方法的可靠性。采用加裝降噪孔板、減小流量和增大進口壓力的降噪措施研究噪聲控制機理,并為偏心旋轉閥的高噪聲問題提供合理的降噪設計方案。

1 多步驟混合計算方法

本文采用多步驟混合計算方法開展偏心旋轉閥氣動噪聲的數值模擬計算,應用ANSYS Fluent計算流體力學軟件中的SST-DDES模型進行瞬態流場計算,獲得壓力、速度、密度等流場信息,然后通過MSC Actran聲學軟件的M?hring聲類比模型將流場信息轉換為氣動噪聲源,通過網格映射將平均流場信息、氣動聲源和湍流壁面壓力映射到聲場流體域網格上,并采用傅里葉變換進行時域和頻率間變換,如圖1所示。

圖1 多步驟混合計算方法流程圖

1.1 控制方程

SSTk-ω湍流模型方程[10]如式(1)所示:

式中:WT為守恒變量;Fc,T為對流項;Fv,T為黏性項;QT為湍流運輸方程的源項。

基于SSTk-ω兩方程湍流模型定義DES 尺度,替換k方程中的耗散項實現SST-DDES模型[11],如式(2)所示:

式中:Δ為3 個方向上的最大網格步長;f1由湍流模型給出

聲類比方法最初由Lighthill[12-13]提出,通過對N-S方程進行變換重組,將其寫成波動方程格式,忽略了流動對聲波的影響。M?hring[14]使用類似方法對此進行了擴展,提出了M?hring聲類比模型,該模型考慮了聲源區域介質的可壓縮性和聲輻射區域介質的對流效應,更適用于高馬赫數(Ma>0.3)氣動噪聲計算[8]。M?hring聲類比方程如式(3)所示:

式中:ρ為密度;ρT為總密度;vi為速度;b為標量焓值;c2=為聲速為旋渦。

式(3)中等式左邊(LHS)為聲傳播項,右邊(RHS)為聲源項,對RHS 進行傅里葉變換,得到M?hring 聲類比頻率方程,然后應用有限元法、部分積分法及高斯定理得到有限元方程,如式(4)所示:

1.2 幾何模型和工況參數

偏心旋轉閥的幾何模型由閥體、閥座、閥桿、V球、進口管道和出口管道組成,材料簡化為結構鋼,密度為7 850 kg/m3,泊松比為0.3,楊氏模量2.1×1011Pa。閥門公稱通徑為DN150,額定流量系數為380,流量特性曲線為近似等百分比,壓差比系數為0.4,閥門類型修正系數為0.25,壓力恢復系數為0.9,介質由閥座方向流入。建立幾何模型時忽略了介質非流動區域結構,并簡化填料部件、底軸、底蓋等結構??紤]氣動噪聲監測點位置要求[15]和工業管道對閥門的影響,進口管道長度為2倍管道公稱通徑[16],出口管道長度為1.3 米[17],總長共計1 828.5 mm,如圖2所示。

圖2 偏心旋轉閥幾何模型

合成氣洗滌塔的高壓差壓力調節用偏心旋轉閥以合成氣為介質[18],介質參數和實際工況參數均按某煤氣化裝置實際參數選定,摩爾質量為24.344 kg/kmol,進口壓力為2 MPa(G),出口壓力為0.5 MPa(G),流量為10 kg/s,溫度為180℃,相應的等熵指數為1.38,操作密度為13.5 kg/m3,黏性系數為1.54×10-5kg/(m·s),聲速為469.23 m/s。

1.3 網格模型

網格模型包括瞬態流場模擬所需的流體域網格和聲場模擬所需的聲流體域網格、結構域網格和聲輻射域網格,如圖3所示。

圖3 聲場網格模型

流體域網格采用poly-hexcore 網格劃分技術進行劃分,在近壁面的流動區域劃分邊界層,采用最小面網格尺寸Smin細化處理流動變化劇烈區域,采用最大面網格尺寸Smax控制網格數量。通常網格數量越多計算精度越高,但計算時間越長,網格單元達到一定數量后計算結果基本不變,這被稱為網格無關化[19]。本文以流量值為網格無關化的衡量指標。如表1 所示,采用3 組網格尺寸建立40°開度的流場有限元模型a、b、c,基于上述介質和工況參數,選用SSTk-ω湍流模型開展穩態流場數值模擬,最大流量誤差為0.05%。因此,考慮模擬計算精度、時間成本和工作量,本文認為模型b已達到網格無關化。

表1 流場流體域網格無關化模擬結果/(kg?s-1)

聲流體域網格和聲輻射域網格采用MSC Actran軟件的波長計算工具和網格劃分工具進行處理,得到最高頻率為5 000 Hz 的聲場流體域網格單元大小為10 mm,單元數為399 130,節點數為89 976,聲場輻射域網格單元大小為20 mm,單元數為514 739,節點數為119 417,滿足了最高頻率對應的聲波波長內劃分6~8 個單元的要求??紤]在氣動噪聲計算中結構域網格與聲流體域的流固耦合,以及與聲輻射域的聲振耦合作用,結構域網格密度需滿足振動波長內8~10 個網格單元的要求,本文采用ANSYS Meshing軟件劃分結構域網格模型,設置網格單元大小為5.5 mm,單元數為676 218,節點數為160 829。

1.4 數值模擬監測點

為了監測瞬態流場中壓力、流量、馬赫數等參數隨采樣時間的變化趨勢,同時獲得偏心旋轉閥氣動噪聲模擬值和聲指向性曲線,本文以偏心旋轉閥閥體中心為原點,以介質流動方向為X軸正方向,以豎直向上的安裝方向為Y軸建立坐標系。流場對稱面為XZ平面,流場監測點“B”坐標為(1114.5,0,0)。氣動噪聲模擬值監測點“A”坐標為(1 114.5,1 084.15,0),符合文獻[11]中的標準規定;以監測點“B”坐標為圓心,以1 084.15 mm 為半徑建立360 個監測點,用以觀察氣動噪聲在遠離閥門出口1米處的聲指向性,如圖4所示。

圖4 數值模擬監測點示意圖

2 數值模擬結果與分析

2.1 流場數值模擬結果與分析

采用多步驟混合計算方法建立40°開度偏心旋轉閥及附屬管道的有限元模型(簡稱為N 模型),以合成氣為介質模擬穩態流場數值,以實際進出口壓力為邊界條件,開啟能量方程,選用SSTk-ω湍流方程,將監測“B”坐標點的壓力脈動曲線作為流場收斂性判斷依據。將收斂后的穩態流場結果作為瞬態流場計算的初始值,選用SST-DDES模型,設置瞬態流場采樣步數為800 步,采樣步長為0.0001 s,根據Nyquist 采樣定理可知,聲場計算的最大頻率為5 000 Hz,最小頻率為12.5 Hz。

提取0.04 s 時XZ平面的馬赫數,如圖5 所示。合成氣以2 MPa(G)進口壓力流入進口管道,此時馬赫數約為0.1,合成氣處于低馬赫數流動狀態,因管道內面積不變,該區域內的馬赫數基本保持恒定;隨后,合成氣進入閥體內腔區域,馬赫數迅速上升至約1.7,到達出口管道后馬赫數逐漸趨于穩定,馬赫數約為0.35。由文獻[12]可知,高馬赫數氣體所在的區域極容易產生沖擊波、噴射流、旋渦流等凌亂流體,這種流體產生的氣動噪聲沿著出口管道傳遞到各處,嚴重時將因振動過大而破壞管道系統。

圖5 0.04s時XZ平面的馬赫數散點圖(N模型)

N 模型在0.04 s 時的流線云圖和渦量圖如圖6所示。流線云圖采用150 條流線繪制,渦量圖采用Q-Criterion方法進行繪制。流線以恒定流速進入閥座區域,V球的阻擋使得流動方向發生改變,因閥體內腔受限,在擠壓效應的作用下,狹窄區域的流速迅速增加至533.6 m/s,V 球兩側位置均出現了沖擊射流和回流現象,圍繞射流的流線分布較混亂。大量的渦旋結構分布于閥體內腔和出口位置并呈脫落現象。這些渦旋對偏心旋轉閥內部結構產生的不平衡力容易引起振動,渦旋中心容易產生氣穴;渦旋脫離的頻率與閥門及附屬管道結構的固有頻率接近或相同時容易產生共振現象,這些都是氣動噪聲源[20]。

圖6 0.04 s時的流線云圖和渦量圖(N模型)

2.2 聲場數值模擬結果與分析

提取2.1 小節瞬態流場的速度、壓力和密度信息,采用M?hring 聲類比方法將流場信息轉換為聲源。采用直接頻率分析方法,分別設置聲場流體域網格和聲輻射域網格為有限域組件,結構域網格為固體域組件,聲輻射域外壁面為無限域組件,流體域進出口面為自由模態管組件,在結構域進出口面設置位移約束,在聲場流體域施加聲源,以此實現閥門氣動噪聲的聲流固耦合計算。

以2×10-5Pa 作為參考聲壓,繪制N 模型“A”監測點在頻率內的A加權聲壓級曲線和1/3倍頻程圖,如圖7 所示。聲壓級曲線呈上下波動狀態,但總趨勢是隨著頻率的增加而下降;聲壓級主要分布在500 Hz~3 000 Hz區間,而又以500 Hz~1 500 Hz為主,“A”監測點的最大聲壓級99.3 dB(A)對應于頻率500 Hz,其次是99.3 dB(A)對應1 313 Hz,處于中低頻區域;聲壓級在高頻區呈明顯的寬頻特性。對頻率內的A 加權聲壓級進行總聲壓級計算,獲得距偏心旋轉閥出口和管壁1 米處(“A”監測點)的氣動噪聲數值模擬值為114.6 dB(A)。

圖7 N模型“A”監測點聲壓級和1/3倍頻程聲壓級曲線

提取最大聲壓級對應頻率500 Hz和1 313 Hz下距閥門出口1米處的360個監測點的聲壓級,并轉換為極坐標的形式,其中Y軸正方向對應極坐標的0°,Z軸正方向對應極坐標的90°,如圖8所示。

圖8 N模型的最大聲壓級對應頻率下的聲指向性曲線

聲壓級曲線呈中心對稱分布,500 Hz 的聲壓級曲線表現出顯著的偶極子聲源指向特性,最大聲壓級位于偏離Z軸45°方向;最小聲壓級位于偏離Y軸50°方向。1 313 Hz 的聲壓級曲線表現出顯著的四極子聲源指向特性,最大聲壓級分別位于偏離Y軸60°和150°對稱軸;最小聲壓級分別位于偏離Y軸20°和110°對稱軸。

提取聲場流體域和聲輻射域的聲壓和結構域的位移,采用Actran 軟件后處理模塊獲得500 Hz 和1 313 Hz頻率下XZ平面流體域聲壓云圖、聲輻射域聲壓云圖和結構域位移云圖,如圖9 所示。聲輻射壓力云圖和結構域位移云圖基本對稱;聲輻射域的最大聲壓主要分布于出口管道附近,說明了閥體內腔和出口位置渦旋結構的形成和脫落引起的噪聲沿出口管道傳播到各處,500 Hz 的聲輻射域的聲壓主要由偶極子產生,1 313 Hz 的聲輻射域的聲壓主要由四極子產生。流體域聲壓明顯高于聲輻射域聲壓,這是由于噪聲在穿過閥門和管道結構時產生了傳遞損失;結構域振動位移變化源自偶極子聲源和四極子聲源的共同作用。

圖9 N模型的聲壓云圖和振動位移云圖

綜上所述,介質為合成氣,進口壓力為2 MPa(G),出口壓力為0.5 MPa(G),開度為40°時偏心旋轉閥在“A”監測點的氣動噪聲模擬值為114.6 dB(A),聲壓級主要分布在500 Hz~3 000 Hz 區間,并以中低頻500 Hz~15 00 Hz 為主,最大聲壓級對應頻率為500 Hz,其次為1 313 Hz。聲輻射域的最大聲壓位于出口管道附近;此數值模型的氣動噪聲聲源包含了介質與閥門和管道內壁的相互作用形成的偶極子聲源,以及閥門和管道內渦旋結構形成的四極子聲源。

2.3 數值模擬方法驗證

Conval 軟件是一款專業的閥門選型軟件,其氣動噪聲標準值的算法來源于標準《IEC60534-8-3:2010》[21]。根據第1 節所述偏心旋轉閥的結構參數、介質參數和工況參數進行氣動噪聲標準值計算,結果表明:偏心旋轉閥產生的氣動噪聲處于Ⅳ狀態,此時,馬赫面形成、分子碰撞減少、激波紊流作用占主要因素;閥門計算流量系數為173.4,開度為40°,出口管道位置馬赫數為0.32,節流位置最大馬赫數為1.72,管壁外1 米處的A 加權聲壓級為117 dB(A)。與數值模擬結果比較可知:標準噪聲值與氣動噪聲模擬值誤差小于3 dB(A);閥門開度基本一致;出口管道位置處馬赫數誤差約為8%;節流位置的最大馬赫數誤差約為2%?;谝陨辖Y論判定該數值模擬方法具有可靠性。

3 氣動噪聲控制機理研究

3.1 閥后加裝降噪孔板

在偏心旋轉閥與出口管道間位置安裝降噪孔板是一種常見的噪聲控制方法,該方法基于聲路處理法原理,在保持偏心旋轉閥原有結構的基礎上,降低從聲源到監測點的傳播聲路的聲壓級[20]。本文對第2節所述N模型的氣動噪聲進行降噪處理,基于數值模擬閥門開度和邊界條件一致性,在閥門出口位置加裝厚度為10 mm、線性分布有345 個?5 mm 小孔的降噪板,建立閥后加裝降噪孔板后的有限元模型(簡稱“K 模型”),并按1.3 節網格繪制方法劃分網格,在降噪孔板部分進行邊界層劃分及局部加密處理,流體域網格如圖10所示。

圖10 K模型的流體域網格模型

繪制K 模型在0.04 s 時的流線云圖和渦量圖如圖11 所示。比較圖6 可知,閥后加裝降噪孔板的K模型有效抑制了閥體內腔和出口位置的沖擊射流和回流現象,流場的最高流速由526.2 m/s 減小到433 m/s,且最大流速位置轉移至降噪孔板出口位置,有效保護了閥體內腔的零部件。K模型的渦旋結構主要分布于閥體內腔和降噪孔板的出口位置。渦旋結構的不斷形成和脫落導致的壓力脈動作用于閥體內腔和管道內壁,因此產生振動進而導致噪聲的產生。

圖11 0.04 s時的流線云圖和渦量圖(K模型)

提取K 模型在采樣時間內瞬態流場監測點“B”的壓力信息,以及采樣頻率內監測點“A”的聲壓信息,對比N模型繪制壓力脈動曲線和聲壓級曲線,如圖12所示。采樣時間內N模型在監測點“B”位置的最大壓差為6.9 kPa(G),K模型的最大壓差為3.5 kPa(G)。因此,閥后加裝降噪孔板可以通過抑制渦旋結構的形成和脫落,減小壓力脈動而控制噪聲的產生。采樣頻率內的聲壓級曲線均隨頻率呈整體下降的趨勢,最大聲壓級對應頻率由500 Hz轉移到1 800 Hz,仍處于中低頻區域;K 模型的聲壓級明顯低于N 模型;對頻率內的聲壓級進行總聲壓級計算,獲得“A”監測點的氣動噪聲模擬值為106 dB(A)。因此,閥后加裝降噪孔板通過抑制閥體內腔和出口管道內的渦旋結構的生成和脫落,降低壓力脈動強度,從而實現降噪。

圖12 K模型和N模型的壓力脈動和聲壓級對比曲線

3.2 減小流量

以合成氣為介質,進口壓力為2 MPa(G),出口壓力為0.5 MPa(G),以流量為10 kg/s、6 kg/s、2 kg/s 3種工況參數為邊界條件,采用多步驟混合數值方法分別建立有限元模型(N 模型、O 模型、P 模型)并開展聲流固耦合數值模擬。

提取O模型和P模型在采樣時間內瞬態流場監測點“B”的壓力信息,以及采樣頻率內監測點“A”的聲壓信息,對比N 模型繪制壓力脈動曲線和聲壓級曲線,如圖13 所示。N 模型的壓力脈動曲線的最大壓差為6.9 kPa(G),O 模型為4.8 kPa(G),P 模型為3.28 kPa(G),壓力脈動隨著流量的減小而降低。3種模型的聲壓級曲線趨勢基本一致,最大聲壓級對應頻率分別為500 Hz、1 325 Hz、1 938 Hz,均處于500 Hz~3 000 Hz 頻率的中低頻區域,且最大聲壓級隨流量減小而降低。對頻率內的聲壓級進行總聲壓級計算,獲得“A”監測點的氣動噪聲模擬值分別為102.6 dB(A)和88.1 dB(A)。因此,流量減小后,壓力脈動強度變弱,噪聲也隨之變小。

圖13 不同流量下的壓力脈動和聲壓級對比曲線

3.3 增大進口壓力

以合成氣為介質,出口壓力為0.5 MPa(G),流量為10 kg/s,以進口壓力分別為5 MPa(G)、3 MPa(G)、2 MPa(G)3 種工況參數為邊界條件,采用多步驟混合數值方法分別建立有限元模型(L模型、M模型、N模型)并開展聲流固耦合數值模擬。

提取L模型和M模型在采樣時間內瞬態流場監測點“B”的壓力信息,以及采樣頻率內監測點“A”的聲壓信息,對比N模型繪制壓力脈動曲線和1/3倍頻程聲壓級曲線,如圖14 所示。N 模型的壓力脈動曲線的最大壓差為6.9 kPa(G),M 模型為15.6 kPa(G),L 模型為23.2 kPa(G),壓力脈動隨著進口壓力的增大變強。3種模型的聲壓級在0~1 500 Hz頻率區間無趨勢變化,1500 Hz以上的聲壓級隨著進口壓力的增大而變大,趨勢較明顯。對頻率內的聲壓級進行總聲壓級計算,獲得“A”監測點的氣動噪聲模擬值分別為117.6 dB(A)和119.5 dB(A)。因此,增大進口壓力,壓力脈動強度變大,噪聲也隨之變大。

圖14 不同進口壓力下壓力脈動和1/3倍頻程聲壓級曲線

4 結語

本文以煤氣化裝置中合成氣洗滌塔的高壓差壓力調節用偏心旋轉閥為研究對象,根據實際工況參數,采用多步驟混合計算方法開展氣動噪聲產生及控制機理研究,結論如下:

(1)多步驟混合計算方法基于計算流體動力學CFD 和計算氣動聲學CAA 提出,應用ANSYS Fluent 計算流體力學軟件中的SST-DDES 模型實現瞬態流場的計算,采用Actran 聲學軟件的M?hring聲類比模型、網格映射、傅里葉變換方法實現流場與聲場信息轉換和噪聲求解,獲得了管壁外1 米處的氣動噪聲模擬值,以及瞬態流場和聲場結果信息。

(2)介質為合成氣,進口壓力為2 MPa(G),出口壓力為0.5 MPa(G),流量為10 kg/s 時偏心旋轉閥在管壁外1米處“A”監測點的氣動噪聲模擬值為114.6 dB(A),流場和聲場結果表明:偏心旋轉閥的氣動噪聲來源于介質與閥門和管道內壁的相互作用形成的偶極子聲源,以及閥門和管道內渦旋結構形成的四極子聲源共同作用的結果,其聲壓級頻率主要分布在500 Hz~3 kHz 區間。采用Conval 閥門選型軟件計算得到的噪聲標準值為117 dB(A),兩者計算結果誤差在±3 dB(A)內,證明了數值模擬方法具有可靠性。

(3)基于閥后加裝降噪孔板、減小流量、增大進口壓力3 種噪聲控制方法進行噪聲控制機理研究。研究結果表明:閥后加裝降噪孔板削弱了閥體內腔和出口管道沖擊波和回流的形成,以及出口管道內渦旋結果的形成和脫落,降低了壓力脈動強度,“A”監測點氣動噪聲降低了9 dB(A);流量減小,降低了“B”監測點壓力脈動強度,氣動噪聲隨之減??;增大進口壓力加強“B”監測點壓力脈動,氣動噪聲隨之增加。

猜你喜歡
閥門模型
一半模型
美嘉諾閥門(大連)有限公司
流程工業(2022年3期)2022-06-23 09:41:08
VANESSA始終引領三偏心閥門的未來發展
裝配式玻璃鋼閥門井的研發及應用
煤氣與熱力(2021年3期)2021-06-09 06:16:18
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
核電站閥門緊急采購系統的構建
智富時代(2018年5期)2018-07-18 17:52:04
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
省力閥門瓶蓋
中學科技(2014年11期)2014-12-25 07:38:53
主站蜘蛛池模板: 亚洲午夜福利精品无码不卡| 欧美一级在线| 久久精品一卡日本电影 | 欧美国产日产一区二区| 国产亚卅精品无码| 亚洲视频三级| 四虎影视8848永久精品| 国产第二十一页| 国产精品久线在线观看| 色欲色欲久久综合网| 午夜免费小视频| 欧美特级AAAAAA视频免费观看| 91在线精品免费免费播放| 五月综合色婷婷| 青青青国产视频手机| 激情综合图区| 成人精品午夜福利在线播放| 亚洲女同欧美在线| 国产乱人乱偷精品视频a人人澡 | 无码一区二区三区视频在线播放| 国产精品永久不卡免费视频| 国产成人av一区二区三区| 丰满人妻一区二区三区视频| 日韩福利在线观看| 国产成人综合日韩精品无码首页| 国产午夜一级淫片| 97成人在线观看| 国产精品hd在线播放| 精品人妻一区二区三区蜜桃AⅤ| 天天操精品| 欧美va亚洲va香蕉在线| 久久婷婷人人澡人人爱91| 国产91丝袜在线播放动漫| 波多野结衣一区二区三视频| AV天堂资源福利在线观看| 重口调教一区二区视频| 日韩国产高清无码| 久久久久夜色精品波多野结衣| 国产欧美日韩在线在线不卡视频| 国产自产视频一区二区三区| 美女无遮挡免费网站| 日本欧美成人免费| 熟女日韩精品2区| 国产精品漂亮美女在线观看| 色综合久久久久8天国| 国产无套粉嫩白浆| lhav亚洲精品| 亚洲日本韩在线观看| 久久6免费视频| 国产精品一区不卡| 波多野一区| 米奇精品一区二区三区| 国产午夜人做人免费视频| 91精品日韩人妻无码久久| 国产欧美日韩另类精彩视频| 欧美在线中文字幕| 欧亚日韩Av| 欧美日韩精品在线播放| 国产97公开成人免费视频| 午夜久久影院| 天天综合天天综合| 巨熟乳波霸若妻中文观看免费| 一区二区影院| 一级高清毛片免费a级高清毛片| 亚洲a级在线观看| 免费在线色| 国内精品久久久久鸭| 亚洲va在线观看| 中文无码精品a∨在线观看| 99精品视频在线观看免费播放| 亚洲毛片在线看| 热这里只有精品国产热门精品| 国产精品亚欧美一区二区| 伊人国产无码高清视频| 91欧美在线| 国产精品久久自在自线观看| 91久久精品国产| 欧美黑人欧美精品刺激| 青青青视频91在线 | 国产精品七七在线播放| 久久国产精品娇妻素人| 2019年国产精品自拍不卡|