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

基于大渦模擬的風琴管噴嘴結構參數優化及流動特性分析

2024-07-02 02:48:36楊新霞閆月娟徐艷李森劉海水劉幸倩
化工機械 2024年3期

楊新霞 閆月娟 徐艷 李森 劉海水 劉幸倩

DOI:10.20031/j.cnki.0254?6094.202403009

摘 要 基于ZGB空化模型和WALE亞格子模型,采用大渦模擬(LES)方法對風琴管噴嘴空化射流流場進行數值模擬研究,在實驗與數值模擬空化云演化規律基本吻合的基礎上開展風琴管噴嘴結構參數優化設計。基于三因素四水平的正交試驗設計,分析了圓柱段及擴散段對空化射流流場特性的影響規律。結果表明,擴散角角度為空化效果的主要影響因素,圓柱段長度為射流效果的主要影響因素,擴散段長度對空化射流影響相當,通過極差分析獲得了風琴管噴嘴最優結構參數組合。對比分析風琴管噴嘴結構參數優化前后的空化射流流場發現,優化后的風琴管噴嘴空化射流效果顯著提升。

關鍵詞 風琴管噴嘴 大渦模擬 正交試驗 結構參數優化 空化射流

中圖分類號 TQ055?? 文獻標志碼 A?? 文章編號 0254?6094(2024)03?0382?10

Large Eddy Simulation?based Optimization of

Organ Nozzle Structural Parameters and Flow Characteristics Analysis

YANG Xin?xia, YAN Yue?juan, XU Yan, LI Sen, LIU Hai?shui, LIU Xing?qian

(School of Mechanical Science and Engineering, Northeast Petroleum University)

Abstract?? Based on ZGB cavitation model and WALE subgrid model, having large eddy simulation (LES) method adopted to simulate cavitation jet flow field of the organ nozzle was implemented, including basing on the general coincidence between experimental and numerical simulation of the evolution law of cavitation clouds to optimize the organ nozzles structure parameters. In addition, based on orthogonal experimental design of three factors and four levels, the influence of both cylindrical and diffusion sections on the cavitation jet flowfield characteristics was analyzed to show that, the diffusion angle mainly influences cavitation effect, and the cylindrical sections length mainly influence jet effect, and the diffusion sections length has a similar effect on the cavitation jet. Through analysis of range, the optimal structural parameters combination of the organ nozzle was obtained. Comparing and analyzing the cavitation jet flowfield simulated before and after optimizing structural parameters of the organ nozzle indicates that, the optimized organ nozzle has significantly improved the cavitation jet effect.

Key words?? organ nozzle, large eddy simulation, orthogonal experiment, structural parameter optimization,cavitation jet

(Continued from Page 359)

基金項目:黑龍江省自然科學基金(批準號:LH2022E016,LH2022A004)資助的課題。

作者簡介:楊新霞(1997-),碩士研究生,從事機器學習與空化射流數值模擬研究。

通訊作者:閆月娟(1971-),教授,從事流體機械及內流特性的研究,1072254366@qq.com。

引用本文:楊新霞,閆月娟,徐艷,等.基于大渦模擬的風琴管噴嘴結構參數優化及流動特性分析[J].化工機械,2024,

51(3):382-391.

高壓水射流清洗領域中[1~4],風琴管空化射流噴嘴利用空泡潰滅產生的高溫高壓可以大幅提高射流的沖擊效果,因此近年來得到了廣泛應用。在綠色低碳發展的大背景下,一些石油化工設備由于所處理介質的復雜性,導致產生的污油或垢蝕難以清洗,因此優化設計一種更加高效節能的風琴管噴嘴結構尤為重要。

風琴管噴嘴空化射流具有強烈的多相湍流特征[5]。目前,大多數學者采用RANS方法對風琴管噴嘴進行數值模擬研究,然而RANS模型的渦黏假設僅能對時均化流動進行模擬,缺少瞬時流場信息對空化射流流動的影響。大渦模擬(LES)在求解過程中對大尺度的渦直接求解,對小尺度的渦建模,因此在理論上模擬精度優于RANS模型。王連安等基于大渦模擬方法對風琴管空化噴嘴流場進行數值模擬,發現大渦模擬可以更好地捕捉圓柱段以及擴散段的渦環分布[6]。李貴東等采用大渦模擬方法對射流式離心泵射流器進行研究,并與SST k?ω湍流模型的模擬結果進行對比,結果表明,基于時間平均的SST k?ω湍流模型無法很好地模擬泵喉管內部的壓力脈動情況,大渦模擬結果則能得到壓力脈動隨時間的變化,且與試驗結果一致[7]。暴春航等利用大渦模擬對推進裝置的強剪切空化射流特性進行了數值模擬研究,結果表明,大渦模擬可以捕捉到外流場渦分布以及空化云演變規律。因此,大渦模擬可以更準確地模擬空化射流流場[8]。

空化噴嘴作為典型的空化發生器,其流道結構參數對射流空化程度、空化效果影響巨大。大渦模擬準確度高但計算成本也高,因此針對空化噴嘴結構參數優化問題,大多數學者采用k?ε湍流模型進行數值模擬研究。王萍輝和方湄驗證了風琴管自激振蕩空化噴嘴的清洗效果與噴嘴諧振腔的直徑、長度以及靶距都有關,同時也驗證了用理論公式設計風琴管噴嘴結構是可靠的[9]。郭世建根據理論公式確定了圓柱段直徑,同時入口段與諧振腔的尺寸可根據最佳長徑比確定,采用RNG k?ε模型對風琴管噴嘴圓柱段長度、擴散段長度、擴散角和靶距進行了結構參數優化,得到了最佳參數組合[10]。楊涵根據自振空化射流噴嘴的基本結構,采用RNG k?ε模型對其主要結構參數(如噴嘴擴散角、擴散段長度及圓柱段長度等)對噴嘴性能的影響進行了分析,得到了最佳參數組合[11]。CAI T F等對風琴管噴嘴結構參數進行了實驗研究,結果表明,當擴散角為10~25°時空化效果最好[12]。

綜上所述,LES模型可以更好地捕捉渦環分布和空化云演變規律,采用LES方法可以更精確地優化風琴管噴嘴結構參數,從而提高清洗效率。當壓力不同時,優選出的最佳結構參數可能不同,為此,筆者采用數值模擬研究方法在12 MPa低壓工況下進行風琴管噴嘴結構參數優化。根據理論公式確定圓柱段直徑,借鑒已有的研究成果,采用大渦模擬方法對風琴管噴嘴圓柱段長度、擴散段長度和擴散角角度進行參數優選。對比分析優化前后的風琴管噴嘴空化射流流場特性,從而驗證優化后的風琴管噴嘴結構空化效果,獲得低壓高效節能的風琴管噴嘴結構參數。

1 物理模型

1.1 計算模型

根據蔡騰飛等對風琴管噴嘴結構的優化結果[13],結合實際清洗應用情況,確定了噴嘴初始結構(圖1),主要包括入口腔、諧振腔、圓柱段和擴散段,尺寸參數如下:

入口直徑D 6.4 mm

諧振腔直徑D 3.2 mm

圓柱段直徑d 1 mm

入口長度L 5.2 mm

諧振腔長度L 5 mm

圓柱段長度L 1 mm

擴散段長度L 2 mm

擴散角角度θ 15°

1.2 網格劃分

為保證噴嘴入口來流的充分發展,將入口長度向上游延長至5倍入口直徑,即32 mm。為使流場充分發展,所將噴嘴外部設置為出口段,出口段屬于外流場研究區域。為了避免出口處產生回流從而影響射流區域,將出口長度設為80 mm,直徑為20 mm。風琴管噴嘴空化射流流體計算域如圖2所示。

用ICEM軟件建模,對流體域進行結構性網格劃分,網格單元為六面體。由于圓柱段、擴散段部位屬于發生空化射流的核心區域,所以進行加密處理。網格無關性驗證如圖3所示,網格節點數量分別為72萬、188萬、307萬和420萬。通過對比模擬后的最大氣相體積分數可以發現,當網格數量增加到307萬時,最大氣相體積分數變化不超過2%。綜合考慮計算精度與時間成本,最終選擇307萬網格數作為后續計算模型,其網格劃分如圖4所示。

2 控制方程與模擬參數設置

2.1 大渦模擬控制方程

大渦模擬通過濾波函數將湍流的瞬時運動分解成大尺度渦與小尺度渦,將連續性方程與動量方程濾波后得:

+=0(1)

+=-+(μ-τ)(2)

τ=-2μS+τδ(3)

其中,ρ為混合相密度,[u][~]、[u][~]為i、j方向濾波后的速度分量,[p][~]為濾波后的壓力,μ為混合相動力黏度系數,τ為亞格子應力項,S為平均應變率張量,δ為Kronecker delta(克羅內克)函數,τ為亞網格尺度應力的各向同性部分。

2.2 空化模型

選擇合適的空化模型對準確描述空化噴嘴內部流場具有重要影響。依據文獻[14]的相關研究,Zwart?Gerber?Belamri(ZGB)模型計算精度高,適用于描述復雜流場,且具有較好的收斂性,其蒸發率R與凝結率R表達式如下:

R=F(4)

R=F (5)

其中,F為蒸發率系數,α為氣泡體積分數,α為氣相體積分數,ρ為氣相密度,R為氣泡半徑,p為飽和蒸氣壓,p為局部壓力,ρ為液相密度,F為凝結率系數。

2.3 模擬參數設置

空化射流屬于氣液兩相流動,氣相為水蒸氣,液相為25 ℃的水,液相密度1 000 kg/m3,飽和蒸汽壓3 169 Pa。使用雙精度壓力基求解器進行數值模擬,多相流采用Mixture模型,壁面無滑移速度。空化模型選ZGB模型,亞格子模型選用壁面適應的局部渦黏(WALE)模型。邊界條件分別為壓力入口和壓力出口,入口壓力12 MPa,出口壓力為大氣壓,連續性曲線收斂精度10-5。選用Coupled算法進行速度-壓力耦合,瞬態項方程采用Bounded Second Order Implicit,壓力項為PRESTO!方法,初始化方式選用全局初始化,時間步長設置為10-5 s。

3 數值模擬與實驗空化云對比

游離型空化泡潰滅時,近壁處產生微射流,從而產生速度快、打擊力強的高壓脈沖以實現對污垢的清洗。因此,通過分析空化云演變規律可以預測噴嘴清洗效率。在進行數值模擬時,時間步長設為10-5 s,每10步自動保存一次,即每0.1 ms記錄一次空化云狀態。空化射流實驗臺主要由儲水箱、柱塞泵、壓力表、風琴管噴嘴、實驗水箱及高速攝像機等組成,根據已有實驗研究設置實驗具體參數[6]。為保證空化云發展周期與模擬時間保持一致,高速攝像機曝光時間設為100 μs,拍攝頻率為10 000 fps(1 fps=0.304 m/s)。數值模擬結果與實驗研究空化云演變規律如圖5、6所示。

空化云的發展具有連續的周期性,具有初生期、發展期、脫落期和潰滅期4個階段,從圖5、6可以看出,空化云開始脫落的時間為0.0、0.8、

1.6 ms,脫落周期在0.8 ms左右。0.0~1.5 ms是一個完整的空化云演變周期,空化泡在0 ms時產生,在0.1~0.7 ms時間段內屬于發展期,在0.8 ms開始脫落,脫落后的空化云不斷潰滅直至1.5 ms后消失。將數值模擬結果與實驗研究對比發現,LES模型空化云分布情況與實驗結果基本吻合,驗證了筆者所建立的LES空化射流數值計算模型可以準確地模擬風琴管噴嘴空化射流規律。下面采用LES模型進行正交試驗優選噴嘴結構參數。

4 正交試驗

4.1 正交試驗因素水平及指標確定

正交試驗設計利用正交性挑選出具有代表性的點進行試驗,是一種高效率的試驗設計方法,其利用較少的試驗次數即可得到最優的結構參數組合[15]。筆者采用正交試驗方法進行數值模擬分析,從而更準確地找到風琴管噴嘴結構參數最優組合。

風琴管噴嘴的圓柱段直徑是確定其他參數的依據[16]。圓柱段直徑的最佳理論公式為:

d=()×1.05-b(6)

其中,q為泵的額定流量,p′為泵的額定壓力,n為噴嘴個數,修正系數b的范圍為0.1~0.3。代入相關數據計算可得d=1 mm。

圓柱段直徑與壓力值確定后,研究圓柱段長度、擴散段長度和擴散角角度對空化射流的影響,并對這3個參數進行模擬優選[11]。從空化效果與射流速度兩方面評價噴嘴清洗性能。圓柱段和擴散段為空化產生的核心區域,該位置的空泡含量影響著空化作用的強弱,因此選取核心區域位置的最大氣相體積分數作為評價結構性能的指標。同時,空泡潰滅產生的射流速度是打擊力的關鍵參數,從而影響清洗效果,因此將選取目標靶距的平均射流速度作為評價結構性能的指標。

設計采用三因素四水平(表1)的正交試驗法來確定影響清洗效果的主次因素。

4.2 試驗結果分析

根據正交試驗設計方案,采用大渦數值模擬方法對16組試驗進行數值模擬研究,通過試驗指標的極差分析得到最優參數組合。

4.2.1 最大氣相體積分數極差分析

最大氣相體積分數極差分析結果見表2。

三因素各水平下的最大氣相體積分數如圖7所示。可以看出,擴散角角度(因素C)為最大氣相體積分數的主要影響因素,其值呈現先增大后減小的趨勢,在水平3時,氣相體積分數達到最大值,為97.10%;圓柱段長度(因素A)為最大氣相體積分數的次要影響因素,其值呈現先增大后減小的趨勢,在水平3時,氣相體積分數達到最大值,為96.57%;擴散段長度(因素B)為最弱影響因素,最大氣相體積分數呈現先減小后增大的趨勢,在水平4處達到最大值,為96.72%。

因此,當最大氣相體積分數為性能檢測指標時,最優參數結構的組合為A3B4C3,即最優方案。

4.2.2 目標靶距平均速度極差分析

針對內徑62 mm的管道清洗時,目標靶距選取12 mm,射流打擊力大小由目標靶距位置射流速度決定,因此分析出口位置12 mm處平均速度極差情況,結果見表3。

三因素各水平下的平均速度如圖8所示。可以看出,隨著圓柱段長度的增加,目標靶距位置處平均速度也逐漸增大,在水平4處平均速度達到了最大值,為43.10 m/s;隨著擴散段長度的增加,目標靶距位置處的平均速度呈出出先增大再減小的趨勢,在水平2處平均速度達到了最大值,為43.49 m/s;隨著擴散角角度的增加,目標靶距位置處平均速度先減小后增大再減小,在水平3處平均速度達到最大值,為43.90 m/s。其優選出的組合分布與極差分析一致,因此,當平均速度為性能檢測指標時,最優參數結構組合為A4B2C3,即最優方案。

4.2.3 風琴管噴嘴最優結構參數

圓柱段長度(因素A)對最大氣相體積分數與平均速度的影響相當,但最大氣相體積分數是產生空泡的主要原因,而空泡含量決定了空化作用的強弱,從而對清洗效率產生重要影響,因此選用A3作為最優結果。擴散段長度(因素B)為平均速度影響的主要因素,是最大氣相體積分數影響的次要因素,平均速度影響著打擊力的大小,因此選用B2作為最優結果。擴散段角度(因素C)是最大氣相體積分數影響的主要因素,是平均速度影響的次要因素,由于最大氣相體積分數間接性地決定了空化的強弱,因此選取C3作為最優結果。

綜上,風琴管噴嘴最優結構組合為A3B2C3,具體參數如下:

圓柱段長度 3 mm

擴散段長度 3 mm

擴散角角度 20°

5 噴嘴結構優化前后流場對比

5.1 壓力場

5.1.1 壓力云圖

當射流過程中局部氣壓低于相對溫度下水的飽和蒸氣壓時,部分液體由液態轉變為氣態產生“空泡”,壓力增大到一定程度時,空泡集中破裂產生“空化”現象。因此噴嘴內外流場的壓力分布影響著噴嘴的空化效果。結構參數優化前后模擬得到的負壓云圖如圖9所示。圖9顯示擴散段至出口位置均有負壓產生,其中,中心區域為負壓最大區域,壓力達到-97.8 kPa。優化后的結構負壓區主要在噴嘴出口軸向位置66.8 mm內,負壓范圍為-97.8~-69.8 kPa,壓力波動比較明顯,負壓大且波及范圍廣。初始結構負壓區主要在出口軸向位置59.3 mm內產生,除中心區域外,負壓范圍主要集中在-69.8~-41.9 kPa,負壓值與負壓范圍均較小。負壓大有利于空化產生,負壓范圍大有利于空化泡發展,因此優化后的結構空泡更多,發展更充分。

5.1.2 目標靶距位置壓力分布

在目標靶距位置徑向上創建一條過圓心的線段,近壁面總壓為0,因此截取-5~5 mm區域觀測壓力分布情況,優化前后風琴管噴嘴模擬的壓力分布曲線如圖10所示。可以看出,優化前后壓力分布呈尖塔形狀,分布比較集中,均處于中心位置處。優化前的噴嘴最大壓力為10.8 MPa,優化后最大壓力為12.0 MPa。因此,優化后的風琴管噴嘴打擊力更大,清洗效率更高。

5.2 速度場

5.2.1 速度云圖與矢量圖

速度云圖及矢量圖分布如圖11所示。可以看出,噴嘴優化前后的核心位置處速度均達到

160 m/s并產生射流。隨著射流卷吸、附面層分離和空泡的產生、長大與破裂,射流影響區域不斷增大。壓力與速度發生振蕩在矢量圖中體現為渦流。初始結構模擬的速度矢量主要在徑向9.8 mm范圍內發生振蕩,但振蕩不明顯。優化后的噴嘴結構模擬得到的流場徑向16.6 mm范圍內有明顯振蕩。因此,優化后的噴嘴空化射流流場產生渦量更多。

5.2.2 軸線速度分布

優化前后風琴管噴嘴空化射流流場的軸向速度分布曲線如圖12所示。可以看出,兩種噴嘴結構在入口段軸向速度均無明顯變化,諧振腔處直徑由6.4 mm縮減至3.2 mm,因此速度均輕微增大。圓柱段直徑縮小至1 mm后,速度均急劇增大到最大值160 m/s,產生射流現象,并在噴嘴出口附近出現了明顯的等速段。優化后的噴嘴模擬得到的等速段長度約為初始結構的1.27倍,等速段越長,射流速度越高。在噴嘴出口位置附近初始

結構速度最小值為20 m/s,優化后的結構速度最小值為47 m/s。在噴嘴出口位置附近,優化后的噴嘴模擬得到的軸向速度值振蕩更明顯,渦量更多。因此,優化后的噴嘴在射流過程中,由于速度較大從而產生了更大的打擊力,提高了清洗效率。

5.2.3 外流場不同位置速度分布

優化前后風琴管噴嘴空化射流流場在x=

50 mm與x=80 mm位置徑向速度變化曲線如圖13所示。可以看出,x=50 mm位置處,由于流體高速流出噴嘴且形成等速核,因此振蕩不明顯。初始結構速度達到158 m/s,但射流寬度僅為2.5 mm。優化后的結構速度為155 m/s,在峰值處有輕微振蕩,射流寬度約6 mm,極大地提高了清洗效率。在x=80 mm位置處,由于流體在出口位置充分發展,因此會向邊壁發散,從而引起速度振蕩。初始結構模擬的速度為45 m/s,振蕩范圍較小,僅在徑向-2.5~2.5 mm之間有輕微振蕩。優化后的噴嘴速度為92 m/s,在徑向-5~5 mm范圍內有明顯振蕩。因此,優化后的結構有更多渦產生,空化作用更強。

5.3 空化云對比分析

優化前后同一時刻的氣相分布云圖如圖14所示。可以看出,優化后的噴嘴結構模擬得到的空化云長度約89.8 mm,寬度為12.3 mm。初始結構的噴嘴模擬得到的空化云長度約80.2 mm,寬度為9.2 mm。空化云長度與寬度決定了噴嘴清洗范圍與有效清洗面積。因此,優化后的風琴管噴嘴空化效果顯著提高。

6 結論

6.1 采用LES數值模擬方法,基于三因素四水平正交試驗設計,通過極差分析得到最優結構是圓柱段長度為3 mm,擴散段長度為3 mm,擴散角角度為20°。

6.2 通過對比分析風琴管噴嘴結構參數優化前后的壓力場可知,優化后的結構負壓區集中在出口段66.8 mm內,初始結構負壓區集中在出口段59.3 mm內。優化后的噴嘴空化云發展更充分。

6.3? 通過對比分析風琴管噴嘴結構參數優化前后的速度場可知,優化后的結構等速核長度約為初始結構的1.27倍,等速段越長射流速度越高。優化后的噴嘴射流速度更大,對管道清洗更加有利。

6.4 通過對比分析風琴管噴嘴結構參數優化前后的氣相分布云圖可知,優化后的結構較初始結構的空化云長度增加了9.2 mm,寬度增加了

3.1 mm。優化后的噴嘴空化云長度、寬度都有所增加,因此空化效果更好,有效清洗面積與清洗效率更高。

參 考 文 獻

[1] 李根生,沈忠厚.常壓下淹沒自振空化射流沖蝕巖石效果的試驗研究[J].華東石油學院學報(自然科學版),1987(3):12-22.

[2] 吳同鋒,蔡曉君,劉湘晨,等.基于高壓水射流的管道清洗方案設計[J].化工機械,2017,44(1):43-46.

[3] 劉峰,朱南文,王亞林,等.射流空化技術處理乳化含油廢水的研究[J].石油與天然氣化工,2005,34(5):416-419.

[4] 龐雷,王永強,魯飛,等.鉆井泥漿罐自動化清洗與回收成套設備研制[J].流體機械,2020,48(7):12-15;21.

[5] DESANTES J M,PAYRI R,SALVADOR F J,et al.Influence of cavitation phenomenon on primary break?up and spray behavior at stationary conditions[J].Fuel,2010,89(10):3033-3041.

[6] 王連安,徐艷,王尊策,等.風琴管噴嘴空化水射流流場的大渦模擬[J/OL].機械科學與技術:1-10[2024-05?03].https://doi.org/10.13433/j.cnki.1003?8728.2023

0072.

[7] 李貴東,王洋,楊學明,等.基于大渦模擬的射流式離心泵射流器內部的流動特性[J].排灌機械工程學報,2017,35(5):369-374.

[8] 暴春航,龍新平,梁蘊致,等.高壓水射流噴水推進裝置空化流動及推力特性仿真分析[J].液壓與氣動,2023,47(1):32-43.

[9] 王萍輝,方湄.空化水射流清洗的實驗研究[J].煤礦機械,2004(1):37-40.

[10] 郭世建.基于空化射流的油管內壁清洗技術研究[D].沈陽:東北石油大學,2018.

[11] 楊涵.空化射流噴嘴流場數值模擬研究[D].成都:西華大學,2015.

[12]?? CAI T F,PAN Y,MA F.Effects of nozzle lip geometry on the cavitation erosion characteristics of self?excited cavitating waterjet[J].Experimental Thermal and Fluid Science,2020,117:110137.

[13] 蔡騰飛,潘巖,馬飛,等.噴嘴出口結構參數對風琴管射流空化作用的影響[J].機械工程學報,2019,55(18):150-156.

[14] 曹東剛,何國強,潘宏亮,等.三種空穴模型在可調汽蝕文氏管數值模擬中的對比研究[J].西北工業大學學報,2013,31(4):596-601.

[15] 李云雁,胡傳榮.試驗設計與數據處理[M].第2版.北京:化學工業出版社,2008.

[16] 李根生,沈忠厚,周長山,等.自振空化射流研究與應用進展[J].中國工程科學,2005,7(1):27-32.

(收稿日期:2023-04-21,修回日期:2024-05-07)

主站蜘蛛池模板: 又爽又大又黄a级毛片在线视频| 国产成人综合亚洲网址| 波多野结衣一区二区三区88| 99re视频在线| 久久免费精品琪琪| 爽爽影院十八禁在线观看| 精品在线免费播放| 日韩欧美国产中文| 国产剧情一区二区| 国产尤物视频在线| 欧美一道本| 免费看的一级毛片| 亚洲午夜综合网| 综合人妻久久一区二区精品| 亚洲AV无码不卡无码| AV不卡无码免费一区二区三区| 四虎成人免费毛片| 日本一区二区三区精品国产| 国产精品视频导航| 国产欧美在线观看视频| 国产不卡网| 一区二区午夜| 亚洲午夜福利精品无码不卡| 国产微拍一区| 精品欧美一区二区三区久久久| 天天综合网亚洲网站| 免费国产黄线在线观看| 国产爽歪歪免费视频在线观看| 97se亚洲综合| 成人av专区精品无码国产| 真实国产精品vr专区| 国产成人区在线观看视频| 天天色天天操综合网| 久草青青在线视频| 欧美啪啪视频免码| 精品人妻无码区在线视频| 大陆精大陆国产国语精品1024| 国产特级毛片| 亚洲色欲色欲www网| 亚洲日本精品一区二区| 国产91精品调教在线播放| 99中文字幕亚洲一区二区| 久久精品无码国产一区二区三区 | 永久免费无码成人网站| 国产毛片片精品天天看视频| 国产成年无码AⅤ片在线| 久草视频精品| 四虎永久在线精品影院| 久久6免费视频| 免费看a级毛片| 国产激情无码一区二区APP | 97久久精品人人| 5555国产在线观看| 最新精品久久精品| 亚洲精品成人片在线播放| 国产精品尤物铁牛tv | 欧美精品xx| 成人在线综合| 欧美综合区自拍亚洲综合天堂| 欧美在线免费| 亚洲 成人国产| 九九精品在线观看| 日韩精品久久无码中文字幕色欲| 亚洲综合18p| 免费国产小视频在线观看| 亚洲AV人人澡人人双人| 综合亚洲网| 经典三级久久| 青青青伊人色综合久久| 2048国产精品原创综合在线| 99精品热视频这里只有精品7| 国产不卡网| 日韩无码视频网站| 亚洲国产看片基地久久1024| 91精品国产一区自在线拍| 国产在线八区| 欧美a在线看| 97在线碰| 欧美在线精品怡红院| 国产精品视频观看裸模| 欧美综合成人| 日本少妇又色又爽又高潮|