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

R32毛細管流量特性模型計算和諾模圖建立

2023-12-12 03:27:46胡晉珽杜仲星張智鋌吳鵬展
制冷學報 2023年6期
關鍵詞:實驗模型

谷 波 胡晉珽 杜仲星 張智鋌 吳鵬展

(上海交通大學機械與動力工程學院 上海 200240)

毛細管是一種節流降壓元件,具有結構簡單、價格低廉、能降低壓縮機啟動力矩等優點,常應用于小型空調系統中。毛細管與壓縮機、換熱器部件之間的流量平衡會直接影響系統蒸發溫度,使系統輸入功率、制冷量和性能系數均隨之改變。因此,研究毛細管的流量特性非常關鍵,雖然毛細管結構簡單,但內部流動機理復雜,目前已有較多研究學者通過實驗測試和數學建模的方式進行相關研究。在實驗方面,何欽波等[1]實驗研究了R410A在變徑毛細管內的流量特性,并基于實驗數據提出經驗模型,確定冷暖空調上使用的毛細管結構參數。劉業鳳等[2]搭建了CO2跨臨界流在毛細管內流動特性的測試實驗臺,同時使用無量綱分析法提出了毛細管流量計算和長度計算的經驗關聯式,并繪制毛細管選型圖。Y. J. Youn等[3]對毛細管內出現的壅塞流進行深入研究,實驗測試了絕熱毛細管內振蕩壅塞流的液膜厚度,發現加減速效應會對液膜厚度產生影響,并進一步修正了所提出的液膜厚度預測經驗關聯式[4],整體精度誤差小于15%。

在建模仿真方面,一些學者通過建立數學模型,對毛細管內的質量流量進行迭代求解,M. M. Bolstad等[5]提出了基于熱力學平衡狀態下的均相流模型,并在均相流動基礎上、摩擦因子為常數的假設條件下求解絕熱毛細管的解析解,通過對質量、能量、動量守恒方程用較為簡單的方法處理,得到毛細管流動方程。S. M. Liang等[6]將漂移流模型應用于兩相區計算,發現與實際工況吻合程度較高。劉海峰等[7]使用非線性法對冰箱毛細管建立計算模型,得到冰箱實際運行工況下毛細管內R12及R134a的流動特性,同時提出一種具有回熱的毛細管優化方法。尚峰等[8]采用遺傳算法對毛細管流量特性神經網絡關聯式進行修正,提高了R407C毛細管流量特性計算模型的精度和計算速度。王江翠等[9]基于文獻中的實驗數據,對比分析了11種絕熱毛細管無量綱流量關聯式,覆蓋工質包括R218、R134a、R290、R407C,結果表明,基于機理模型的神經網絡關聯式效果最佳。王棟等[10]基于均相流模型建立了毛細管流通特性數學模型,并搭建了CO2系統測試實驗臺以驗證模擬計算精度,結果表明,流量變化趨勢吻合,數據精度在高排氣壓力下較低。谷波等[11-12]就過冷區和兩相區建立了R407C 絕熱毛細管的能量、動量、質量方程和壓降計算公式,同時還分析了絕熱毛細管中壅塞流對制冷劑流量的影響及在計算過程中的處理方法。邵雙全等[13]在對空調系統性能影響進行數值仿真分析的基礎上,建立了空調器制冷劑充注量與毛細管長度的優化匹配模型。

此外,也有一些學者基于大量的實驗數據建立了毛細管內流量的顯性計算公式。一種是π關聯式,文獻[14-17]均基于該方法建立了不同結構的毛細管質量流量經驗關聯式。此類計算方法的優點是使用簡單、精度較高,缺點是推廣性較差,難以適用于不同的制冷劑。另一種是基于ASHRAE諾模圖的快速計算方法,翟千鈞等[18]采用此方法,針對4張毛細管諾模圖中的標準流量、流量系數、臨界壓力系數提出了顯性計算公式,該方法相比于π關聯式還可以判定毛細管出口是否為臨界流。

目前針對制冷劑在毛細管內的流動特性及數學模型已有較多研究,但在實際的工程應用中,相比于各種計算模型,從業人員一般更傾向于直接查找諾模圖以獲得制冷劑在毛細管內的流量。R32制冷劑是一種環境友好型制冷劑,被認為是R410A眾多替代品中最有前景的一種,目前在家用空調中已有較多應用,但目前仍無R32毛細管諾模圖,這對從業人員帶來諸多不便。因此,本文基于毛細管均相流假設及能量守恒、動量守恒、質量守恒方程,計算了不同工況下R32在毛細管內的流量數據,同時搭建了R32毛細管流量測試平臺,最后,在模型仿真數據及實驗數據的基礎上,給出R32在毛細管內流量計算的諾模圖。

1 R32毛細管快速計算模型

為建立毛細管穩態數學模型,主要采用如下假設:1)流動為一維絕熱流動,忽略徑向參數的變化;2)忽略重力對制冷劑流動的影響;3)毛細管內的兩相流動為均相流,各流動界面均處于熱力學平衡狀態。

1.1 基本模型

質量守恒:

(1)

動量守恒:

(2)

能量守恒:

(3)

式中:ρ為密度,kg/m3;u為流速,m/s;x為長度,m;p為壓力,kPa;h為焓值,J/kg;f為摩擦系數;di為管內徑,m。

若將模型離散化后,則每個微元內相應的差分方程為:

質量守恒:

ρ1u1=ρ2u2

(4)

動量守恒:

(5)

能量守恒:

ρ1u1(h2-h1)+ρ1u1(u22-u12)/2=0

(6)

式中:ΔL為微元長度,m;下標1、2分別表示微元進、出口;fm為平均摩擦系數。

1.2 求解方法

本文建立的毛細管仿真模型,通過在各結構參數確定的情況下假設毛細管流量,并以此計算毛細管管長,與實際管長進行對比,當兩者誤差較小時,即認為此流量為實際流量。以下分為過冷區和兩相區敘述毛細管管長的求取。

1.2.1 過冷區解法

毛細管進口截面上的壓力p1受到流通面積減小和局部阻力的影響小于進口前的壓力p0,計算如下:

p0υ0-p1υ1=(1+ξ01)u2/2

(7)

由于過冷區內制冷劑密度隨溫度的變化較小,且制冷劑流量保持不變,因此在過冷區內制冷劑流速基本保持穩定,將過冷區看作等溫、等焓、等速流動。則毛細管在過冷區的長度計算如下:

L過冷=2(p1-psat)di/(fρ1u12)

(8)

式中:ξ01為進口截面的局部阻力系數;υ為比容,m3/kg;下標sat表示飽和狀態。

1.2.2 兩相區解法

毛細管在進入兩相區后,為得到更準確的出口參數需采用分布參數模型。在選擇網格劃分依據時,考慮毛細管流動過程中可能存在臨界流現象。當發生臨界流時,毛細管流量即出口參數不會改變,此時毛細管出口會發生自由膨脹,出口的臨界壓力自由膨脹至背壓。當接近臨界點時,即使毛細管的長度發生微小變化,也會引起壓力參數的顯著變化,因此為保證精度,選擇以微小壓力降為單位劃分網格。

微元出口兩相區比容,計算如下:

υ2=xυ2g+(1-x)υ2f

(9)

通過流量不變方程求取出口流速:

u2=Gυ2/A

(10)

式中:G為質量流量,kg/h;A為毛細管截面積,m2。

出口焓值根據能量守恒方程計算如下:

h2=h1+(u12-u22)/2

(11)

兩相區微元長度可通過微元段的進出口參數計算:

(12)

1.3 壅塞流

當制冷裝置處于正常運行工況時,毛細管出口制冷劑為兩相狀態,當流速達到當地音速時會發生流動壅塞。發生流動壅塞后,毛細管背壓的變化將不再影響制冷劑在毛細管內的流動。因此,需在模型中判斷是否發生壅塞現象,以避免計算失真。由于節流過程為不可逆過程,必須滿足熵增原理,因此,當式(13)等號成立時,則認為流動達到壅塞。

ds≥0或ds/dp≤0

(13)

為使模型更方便計算,將式(13)化為式(21)等價形式。根據當地音速的定義:

(14)

(15)

引入兩相區的比容公式得:

(16)

(17)

將式(16)帶入式(17),整理得:

(18)

(19)

(20)

綜上,式(13)的等價形式為:

G≤Gc

(21)

式中:s為比熵,kJ/(kg·K);c為音速,m/s;下標g、f、s分別為氣體、液體、等熵狀態。

由式(18)~(20)可知,臨界質量流量僅取決于制冷劑熱力性質和出口干度。利用臨界流量判據(21),僅需計算毛細管出口處的比熵,即可得到毛細管臨界流量。此時可直接將對應背壓的臨界流量值作為流量的初始值,提高了算法的效率。

1.4 摩擦系數

在理論守恒模型中,摩擦系數的選擇是影響流量計算的重要因素,直接關系到模型求解的準確性和可靠性。本文采用S. W. Churchill[19]關聯式,該關聯式覆蓋范圍廣,且考慮了內管壁粗糙度的影響。

過冷區:

f=8[(8/Re)12+1/(Y+Z)3/2]1/12

(22)

(23)

Z=(37 530/Re)16

(24)

兩相區:阻力系數的計算式仍為式(22),但雷諾數、動力粘度計算如下:

Ret=Gd/μt

(25)

μt=xμg+(1-x)μf

(26)

式中:Re為雷諾數;Y、Z為過程計算參數;ε為實際管內徑,m;ε/di為相對粗糙度;μ為動力粘度,Pa·s;下標t為兩相狀態。

1.5 R32熱力性質

R32物性參數數據來源于軟件REFPROP 9.0,本文基于邱琳禎[20]提出的關聯式基礎形式進行擬合。

飽和氣體蒸氣壓力:

(27)

式中:T為溫度,K;a1~f1為飽和氣體蒸氣壓力系數,取值如表1所示。

表1 飽和氣體蒸氣壓力系數

飽和液體密度:

ρsat=a2+b2(1-T/Tc)1/3+c2(1-T/Tc)2/3+

e2(1-T/Tc)4/3+f2(1-T/Tc)1/2+

g2(1-T/Tc)2

(28)

式中:Tc為臨界溫度,K;a2~g2為飽和液體密度系數,取值如表2所示。

表2 飽和液體密度系數

采用馬丁候狀態方程式描述R32在飽和氣體線即過熱區溫度-壓力-比容間關系:

(29)

式中:V為氣體體積,m3;k、b0均為特性常數,取值分別為5.475、3.481×10-4;R為摩爾氣體常數,J/(mol·K);Ai~Ci為狀態方程系數,取值如表3所示。

表3 狀態方程系數

根據Clapeyron-Clausius方程式和式(26)可推導出蒸發潛熱計算式:

(30)

式中:r為蒸發潛熱,kJ/kg;J為單位轉換系數,取值為100;υ′為飽和液體比容,m3/kg;υ″為飽和氣體比容,m3/kg;b1~f1取值見表1。

氣體定壓比熱容:

cpg=

(31)

式中:cpg為氣體比定壓熱容,J/(kg·K);Ag1~Ag12為氣體比定壓熱容系數,取值如表4所示;ΔT為溫差,℃。

表4 氣體比定壓熱容系數

飽和液體焓:

hf=Al0+Al1(1-T/Tc)+Al2(1-T/Tc)2+

Al3(1-T/Tc)3+Al4(1-T/Tc)4+

Al5(1-T/Tc)5

(32)

式中:Al0~Al5為飽和液體焓系數,取值如表5所示。

表5 飽和液體焓系數

飽和氣體及過熱區氣體焓:

(33)

式中:Ah1~Ah11為飽和液體焓系數,取值如表6所示。

表6 飽和氣體及過熱區氣體焓系數

飽和氣體及過熱區氣體熵:

(34)

式中:As1~As11為飽和氣體及過熱區氣體熵系數,取值如表7所示。

表7 飽和氣體及過熱區氣體熵系數

1.6 計算方法

首先假定臨界流量,調用毛細管長度計算模塊將所得的毛細管長度與給定長度對比,當兩者相等時,得到毛細管臨界流量,同時獲得臨界壓力。再對比臨界壓力與給定背壓的大小,判斷毛細管是否處于臨界流狀態。若處于臨界流,則實際流量為臨界流量,反之以背壓作為出口壓力計算流量。

2 R32毛細管流量測試實驗

2.1 實驗裝置

本實驗裝置由R32住宅空調改造而成,如圖1所示。R32制冷劑通過滾動轉子壓縮機在系統內循環,在壓縮機出口設置高效油分離器以防止潤滑油進入毛細管,影響實驗結果。將額定功率為3 kW的套管式換熱器作為過冷器,以保證制冷劑為過冷液體,并安裝干燥過濾器去除系統中的外來雜質。在過冷器出口設置科氏流量計測量制冷劑質量流量。在毛細管前放置電動膨脹閥調節制冷劑進入毛細管的流量大小,如圖2所示。實驗中通過冷卻水流量、電子膨脹閥等調節毛細管入口參數。實驗樣品是兩根長為0.5、1.0 m,內徑均為1 mm的銅管。室內側實驗系統如圖3所示。

圖1 實驗系統原理

圖2 毛細管段

圖3 室內側實驗系統

本實驗參考熱泵的實際工作工況,在冷凝壓力為2.4~3.2 MPa、過冷度為5~20 ℃工況范圍內進行毛細管流量測試實驗。用鉑電阻和絕對壓力傳感器測量制冷劑溫度和壓力。所有鉑電阻均由標準的汞溫度計提前校準(精度:±0.05 ℃)。測量儀表參數如表8所示。所有數據信號均由數據采集系統(DAS)進行采集和轉換。

表8 測量儀表參數

2.2 實驗與仿真結果對比

本實驗在不同工況和結構下共進行了36組實驗,上述仿真程序計算結果與實驗的對比如圖4所示,其平均相對誤差為5.1%,最大相對誤差為17.1%,所有工況點的相對誤差范圍為-16.7%~17.1%。

圖4 實驗值與仿真值對比

3 R32毛細管諾模圖的建立

毛細管諾模圖通常指研究者對常用制冷劑進行毛細管流量實驗并根據實驗數據繪制的圖表。毛細管諾模圖中,在一個寬泛的運行范圍內,隨著毛細管長度、內徑、制冷劑進口過冷度(或干度)以及進口壓力的改變,均能找到相應的制冷劑流量。因此,諾模圖是一種非常適合工程實際應用快速、準確的毛細管設計和分析數據源。

在諾模圖快速計算模型中,毛細管流量由標準流量(Gbz)、流量系數(Ce)及臨界壓力修正系數(Kc)計算得到。其中,標準流量是指在標準尺寸毛細管(內徑di=1.625 mm、長度L=2.03 m)所計算的流量,毛細管質量流量計算如下:

G=GbzCeKc

(35)

臨界壓力系數是與壓差比(C0)相關的函數,壓差比計算如下:

(36)

式中:pcri為出口臨界壓力,kPa。

R32是一種新型環保制冷劑,其毛細管流量一般需通過仿真模型計算而得,該模型雖然通用性強,但公式為隱式形式,需不斷迭代求解,對程序的編寫要求較高,一般工程人員不易掌握,且存在迭代時計算速度較慢等問題。基于毛細管仿真模型及實驗數據所作的諾模圖可提供快速而準確的R32毛細管設計分析數據。

本文基于仿真結果制作的R32制冷劑諾模圖如圖5~圖8所示,tsc為過冷度(℃)。由R32毛細管諾模圖和實驗結果可知,毛細管內的質量流量主要取決于毛細管結構參數和運行參數,即毛細管長度、內徑和進口壓力、出口壓力、入口過冷度/干度。具體表現如下:1)內徑越大,制冷劑在管內所受流動阻力更小,質量流量會變大;2)R32在長度越大的毛細管內,總沿程阻力更大,則導致質量流量變小;3)R32入口過冷度的變化影響毛細管內過冷段和兩相段的長度,R32在兩相段內的壓降比過冷段更大,所以當過冷度越小或干度越大時,管內兩相段會增加,導致質量流量不斷減小;4)毛細管進出口壓差越大,則單位長度上的壓降越大,這導致了質量流量的增加。

圖5 標準毛細管流量圖(L=2.03 m,di=1.625 mm)

圖6 流量系數

圖7 出口臨界壓力

圖8 臨界壓力修正系數

4 結論

本文基于毛細管內均相流假設及能量守恒、動量守恒和質量守恒方程,建立了R32毛細管絕熱計算模型,得到結論如下:

1)仿真模型中采用臨界流作為壅塞流是否出現的判斷依據,通過仿真模型得到不同工況下R32在毛細管內的流量數據,為諾模圖的建立提供了大量數據基礎。

2)搭建了使用R32的毛細管測試平臺,實驗測試了長度為0.5、1 m,內徑為1 mm的毛細管在冷凝壓力為2.4~3.2 MPa、過冷度為5~20 ℃工況下的質量流量,共測試36組工況點。

3)仿真數據與實驗數據的平均相對誤差為5.1%,最大相對誤差為17.1%。在趨勢性上較為一致:毛細管內徑的增加,長度的減小,過冷度的增加或干度的減小,進出口壓差的增加,均會使R32在毛細管內的質量流量增加。

4)基于理論模型計算所得的數據和實驗數據繪制了針對R32的毛細管諾模圖,避免了傳統隱式計算對編程要求高且計算速度受影響的問題,為制冷空調領域的研究者和工程師提供了快速而準確的R32毛細管設計分析數據。

[1] 何欽波, 閆格尼, 徐言生. 變徑毛細管在R410A冷暖空調器中的流量特性[J]. 制冷學報, 2018, 39(3): 39-43. (HE Qinbo, YAN Geni, XU Yansheng. Flow characteristic of variable diameter capillary tube for heat pump type air conditioner with R410A[J]. Journal of Refrigeration, 2018, 39(3): 39-43.)

[2] 劉業鳳, 賈世偉, 李續, 等. CO2跨臨界循環中絕熱毛細管數值模擬計算與分析[J]. 流體機械, 2016, 44(11): 78-83. (LIU Yefeng, JIA Shiwei, LI Xu, et al. Numerical simulation and analysis on adiabatic capillary tube with CO2transcritical cycle[J]. Fluid Machinery, 2016, 44(11): 78-83.)

[3] YOUN Y J, HAN Y, SHIKAZONO N. Liquid film thicknesses of oscillating slug flows in a capillary tube[J]. International Journal of Heat and Mass Transfer, 2018, 124: 543-551.

[4] YOUN Y J, MURAMATSU K, HAN Y, et al. The effect of bubble deceleration on the liquid film thickness in microtubes [J]. International Journal of Heat and Fluid Flow, 2016, 58: 84-92.

[5] BOLSTAD M M, JORDAN R C. Theory and use of the capillary tube expansion device [J]. Refrigerating Engineering, 1948, 56: 519-523.

[6] LIANG S M, WONG T N. Numerical modeling of two-phase refrigerant flow through adiabatic capillary tubes[J]. Applied Thermal Engineering, 2001, 21(10): 1035-1048.

[7] 劉海峰,何茂剛,陰建民. 替代工質應用于冰箱毛細管的特性研究及優化方法[J]. 西安交通大學學報,1997,41(2): 46-78. (LIU Haifeng, HE Maogang, YIN Jianmin. Study on the characteristics and optimization method of substitute working medium applied to refrigerator capillary tubes[J]. Journal of Xi′an Jiaotong University, 1997, 41(2): 46-78.)

[8] 尚峰, 吳鋼, 李想, 等. 基于遺傳算法-神經網絡的絕熱毛細管流量特性關聯[J]. 低溫與超導, 2014, 42(5): 51-54. (SHANG Feng, WU Gang, LI Xiang, et al. The optimization design of adiabatic capillary flow characteristics based on genetic algorithm-neural networks[J]. Cryogenics &Superconductivity, 2014, 42(5): 51-54.)

[9] 王江翠,金曉辰,邵亮亮,等. 絕熱毛細管無量綱流量關聯式評估[J]. 制冷技術, 2011, 31(3): 30-33. (WANG Jiangcui, JIN Xiaochen, SHAO Liangliang, et al. Correlation evaluation of dimensionless flow in adiabatic capillary tubes[J]. Chinese Journal of Refrigeration Technology. 2011, 31(3): 30-33.)

[10] 王棟, 李蒙, 戚利利, 等. 二氧化碳制冷系統毛細管的設計及實驗研究[J]. 化工學報, 2011, 62(10): 2753-2758. (WANG Dong, LI Meng, QI Lili, et al. Design and experimental study of capillary tube in transcritical carbon dioxide refrigeration system[J]. CIESC Journal, 2011, 62(10): 2753-2758.)

[11] 谷波, 魏躍文, 韓榮梅. 絕熱毛細管的臨界流動分析[J]. 上海交通大學學報, 2001, 35(8): 1187-1190. (GU Bo, WEI Yuewen, HAN Rongmei, et al. Analysis of the critical flow in adiabatic capillary tube[J]. Journal of Shanghai Jiao Tong University, 2001, 35(8): 1187-1190.)

[12] 谷波, 黎遠光, 裴勇華, 等. 絕熱流動的R407C毛細管計算模型[J]. 上海交通大學學報, 2003, 37(7): 1040-1044. (GU Bo, LI Yuanguang, PEI Yonghua, et al. Analysis on the adiabatic flow of R407C in capillary tube[J]. Journal of Shanghai Jiao Tong University, 2003, 37(7): 1040-1044.)

[13] 邵雙全,石文星,李先庭,等. 制冷劑充灌量和毛細管長度對空調系統性能的影響[J]. 低溫工程, 2002(2): 48-53, 61. (SHAO Shuangquan, SHI Wenxing, LI Xianting, et al. Effect of refrigerant filling amount and capillary length on the performance of air conditioning system[J]. Cryogenics, 2002(2): 48-53, 61.)

[14] RASTI M, JEONG J H. A generalized continuous empirical correlation for the refrigerant mass flow rate through adiabatic straight and helically coiled capillary tubes[J]. Applied Thermal Engineering, 2018, 143: 450-460.

[15] PUNIA S S, SINGH J. An experimental study of the flow of LPG as refrigerant inside an adiabatic helical coiled capillary tube in vapour compression refrigeration system[J]. Heat and Mass Transfer, 2015, 51(11): 1571-1577.

[16] DEODHAR S D, KOTHADIA H B, IYER K N, et al. Experimental and numerical studies of choked flow through adiabatic and diabatic capillary tubes[J]. Applied Thermal Engineering, 2015, 90: 879-894.

[17] MITTAL M K, KUMAR R, GUPTA A. An experimental study of the flow of R-407C in an adiabatic helical capillary tube[J]. International Journal of Refrigeration, 2010, 33(4): 840-847.

[18] 翟千鈞, 谷波. 絕熱毛細管流量快速計算模型[J]. 上海交通大學學報, 2010, 44(8): 1140-1144. (ZHAI Qianjun, GU Bo. Fast calculation model of adiabatic capillary tube flow[J]. Journal of Shanghai Jiao Tong University, 2010, 44(8): 1140-1144.)

[19] CHURCHILL S W. Friction-factor equation spans all fluid-flow regimes[J]. Chemical Engineering, 1977, 84(24): 91-92.

[20] 邱琳禎. R32/潤滑油混合物性能及沸騰換熱特性研究[D]. 上海: 上海交通大學, 2020. (QIU Linzhen. Study on properties of R32/lubricating oil mixture and boiling heat transfer characteristics[D]. Shanghai: Shanghai Jiao Tong University, 2020.)

猜你喜歡
實驗模型
一半模型
記一次有趣的實驗
微型實驗里看“燃燒”
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
做個怪怪長實驗
3D打印中的模型分割與打包
NO與NO2相互轉化實驗的改進
實踐十號上的19項實驗
太空探索(2016年5期)2016-07-12 15:17:55
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 中日韩欧亚无码视频| 成人免费午间影院在线观看| 麻豆精品国产自产在线| 91麻豆久久久| 国产黄色片在线看| 国产拍在线| 亚洲中文字幕日产无码2021| 国产精品欧美亚洲韩国日本不卡| 中文字幕在线播放不卡| 成人欧美日韩| 色老头综合网| h网站在线播放| 精品福利视频网| 欧美一级特黄aaaaaa在线看片| 日本一本正道综合久久dvd| 999精品视频在线| 免费午夜无码18禁无码影院| 无码啪啪精品天堂浪潮av| 91精品国产福利| 欧美啪啪视频免码| 国产一级小视频| 91无码人妻精品一区二区蜜桃| 精品夜恋影院亚洲欧洲| 日本三级精品| 国产v精品成人免费视频71pao | 中文字幕欧美日韩| 欧美午夜久久| 99人妻碰碰碰久久久久禁片| 成人免费网站久久久| 中文字幕在线观| 日韩毛片在线播放| 91国内外精品自在线播放| 欧美区一区| 伊人色在线视频| 狠狠五月天中文字幕| 四虎国产精品永久一区| 亚洲精品无码在线播放网站| 久久久久亚洲av成人网人人软件 | 99在线国产| 嫩草国产在线| 亚洲一区二区约美女探花| 国产女同自拍视频| 久草网视频在线| 亚洲美女一区二区三区| 久久精品日日躁夜夜躁欧美| 丁香五月亚洲综合在线| 无码区日韩专区免费系列| 久久这里只有精品66| 99久久免费精品特色大片| 欧美中文字幕无线码视频| 久久精品aⅴ无码中文字幕| 亚洲天堂在线免费| 最新国语自产精品视频在| 免费在线色| 国产91高跟丝袜| 成年午夜精品久久精品| 在线播放真实国产乱子伦| 亚洲清纯自偷自拍另类专区| 九九热在线视频| 久久精品这里只有国产中文精品| 国产成人艳妇AA视频在线| 在线观看国产精品一区| 99在线观看视频免费| 久久精品国产精品青草app| 国产凹凸一区在线观看视频| 国产不卡网| 成人年鲁鲁在线观看视频| 伊人色天堂| 亚洲免费黄色网| 亚洲成人在线网| 欧美a在线| 亚洲天堂网在线播放| 亚洲午夜国产片在线观看| 成年人福利视频| 99久久精品美女高潮喷水| 黄片一区二区三区| 国产一区二区三区在线精品专区 | 日韩欧美国产三级| 成年av福利永久免费观看| 久久不卡国产精品无码| 亚洲成人动漫在线观看| h视频在线观看网站|