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

管道內(nèi)冰漿流動壓降特性模擬和實驗研究

2018-03-27 09:17:15
制冷學(xué)報 2018年2期
關(guān)鍵詞:實驗模型

(1 天津商業(yè)大學(xué) 天津市制冷技術(shù)重點實驗室 天津 300134; 2 北京理工大學(xué) 北京 100081)

冰漿是一種水、冰晶粒子以及凝固抑制劑(比如氯化鈉、乙醇等)組成絮狀混合物,其中冰粒直徑≤1×10-3m,一般冰粒子濃度大于5%時認(rèn)為是非牛頓流體。冰漿作蓄冷介質(zhì)具有良好的熱物性、傳輸性及相變特性[1],因此在建筑制冷、礦井降溫、食品冷藏、醫(yī)療衛(wèi)生和消防運輸?shù)阮I(lǐng)域都有應(yīng)用[2-6]。近年來國內(nèi)外許多研究者對冰漿的流動特性進行了研究,但主要是結(jié)合冰漿流動實驗現(xiàn)象與非牛頓流體流變特性建立冰漿阻力特性預(yù)測模型。常用的非牛頓流體流變模型[7]有Bingham、Herschel-Bulkley、Powder Law和Casson,一些研究者針對這4種流變模型研究冰漿水平管道流動壓降模型,M. Grozdek等[8]實驗分析了在水平直管段內(nèi)含冰率0%~30%條件下的壓降特性。F. Illán等[9]以9%的NaCl溶液作為載流流體,分析了直管內(nèi)的傳熱和壓降特性。V. Ayel等[10]分析了冰漿的流動和傳熱特性,其中流動特性主要研究了直管內(nèi)的壓降特性。A. Kitanovski等[11]同樣分析了冰漿的流動特性。然而,冰漿流動壓降與管徑、管型、流速、冰粒子濃度等有關(guān),這些模型不能兼顧所有因素對流動壓降的影響,而且由于模型采用的勻質(zhì)流假設(shè),只能用于低速分層流動,通用性和精度不高。

隨著計算機技術(shù)和計算流體力學(xué)的發(fā)展,固液兩相流的數(shù)值模擬也隨之有了很大發(fā)展。常用的固液兩相模型有以離散單元法(DEM)為基礎(chǔ)的Euler-Langrian模型[12]和以顆粒動力學(xué)為基礎(chǔ)的Euler-Euler模型[13]。Euler-Euler模型將顆粒相作為連續(xù)相,相對于Euler-Langrian模型計算量大為減小。將冰漿看作固液兩相流體,利用Ansys Fluent研究冰漿的流動壓降特性,但此方面研究尚不多見[14]。B. Niezgoda-Zelasko等[15]提出利用兩相雙流模型研究冰漿阻力特性,但冰粒子黏度模型仍將冰漿等效為均質(zhì)流體,未考慮冰粒子的非均勻分布對黏度的影響,因此數(shù)學(xué)描述與實際流動存在一定的偏差。本文采用以顆粒相動力學(xué)為基礎(chǔ)的Euler-Euler模型模擬冰漿在水平直管、90°彎管和T型管中的無相變流動壓降,并用實驗進行模型驗證。

1 數(shù)學(xué)模型

描述冰漿流體等溫流動的雙流體計算模型[13],數(shù)學(xué)描述可表示為:

1)連續(xù)性方程

(1)

由質(zhì)量守恒方程,得:

(2)

2)動量守恒方程

將冰漿看做液固混合的顆粒相流體,用液相動力參數(shù)描述動量方程[16]:

αqρqg+αq(Fq+Flift,q+FVm,q)+

(3)

固體顆粒相動量方程為:

αsρsg+αs(Fs+Flift,s+FVm,s)+

(4)

液體相與顆粒相之間的作用力可用以下Gidaspow模型來[13]確定。

當(dāng)αl>0.8時,液體-固體交換系數(shù)Ksl有如下形式:

(5)

其中:

(6)

當(dāng)αl≤0.8時,Ksl有如下形式:

(7)

式中:Re為相對雷諾數(shù)。

主相q和第二相p的相對雷諾數(shù):

(8)

3)湍流方程

本文采用k-ε湍流模型,該模型以N-S方程組為基礎(chǔ),同時考慮顆粒相和液相的相互作用[16]。

(9)

(10)

其中:

(11)

Gk,m=μt,m(um+

(12)

μt,m=ρmCμk2/ε

(13)

(14)

式中:C1ε=1.44;C2ε=1.92;σε=1.3;σk=1;Cμ=0.09;下標(biāo)m為混合相;k為湍流脈動動能;μt,m為混合相湍流動力黏度, Pa·s;σk為湍流脈動動能的普朗特數(shù);ε為耗散項。

4)黏度模型

冰漿是一種固液兩相混合物,物理密度、動力黏度以及結(jié)冰點不僅受添加物濃度影響,還與含冰率有關(guān),D.G.Thomas[17]通過反復(fù)實驗和理論推算,得出黏度關(guān)于含冰率的公式:

μ=μl(1+2.5IPF+10.05IPF2+

0.002 73(e16.6IPF-1))

(15)

式中:IPF為含冰率。

2 數(shù)值模擬結(jié)果分析

借助模擬軟件FLUENT進行流體計算,設(shè)定含冰率為10%和20%,冰晶粒子直徑約為1×10-4m,管徑為0.02 m;雙精度并考慮重力影響,流速為0.2~10 m/s;邊界條件在入口處設(shè)為速度入口,出口設(shè)為出流,壁面設(shè)為無滑移;求解器采用一階迎風(fēng)格式的SIMPLE算法;不考慮相變模型但開啟能量源項,收斂精度設(shè)置為10-4,冰漿的基本參數(shù)見表1。

表1 冰漿的基本參數(shù)

1)水平直管

冰漿通過水平直管時,管道縱截面平均壓力如圖1所示,可以看出冰漿在水平直管中流動時壓力一直下降,且入口的壓降明顯比出口大,造成這種現(xiàn)象的原因是入口處擾動較大,能量損失嚴(yán)重。

圖1 水平直管縱截面平均壓力分布云圖Fig.1 The contour of pressure on the longitudinal section of horizontal tube

2)90°彎管

圖2所示為90°彎管縱截面平均壓力的分布云圖。可以看出,拐角處的壓力變化劇烈,壓力損失嚴(yán)重,這是由于在拐彎處存在摩擦和二次流現(xiàn)象[18],引起能量的損失。由于入口處擾動的影響,使入口壓降大于出口壓降。

圖2 90°彎管縱截面平均壓力分布云圖Fig.2 The contour of pressure on the longitudinal section of the elbow tube

圖3所示為彎管彎曲處冰粒子的速度矢量圖,由此看出在彎管的彎曲處存在明顯的二次流現(xiàn)象。當(dāng)流體運動時遇到變向阻礙,彎道外側(cè)近壁面邊界層速度減小,離心力減小,而中心層速度較高,提供向心力的壓力也較高,因此造成中心層向外層流動的趨勢,形成附加流動,與此同時,還形成外側(cè)向內(nèi)側(cè)的流動以保持質(zhì)量守恒。二次流是主流與螺旋流動的疊加,因此能量損失極大。

圖3 90°彎管的拐彎處管道橫截面速度矢量云圖Fig.3 The contour of velocity on the cross section of elbow tube

圖4所示為90°彎管的管道縱截面冰漿速度分布云圖,可以看出在流體轉(zhuǎn)彎之前,固體冰粒速度逐漸增大,且有明顯的邊界層。在彎管彎曲處,有一個明顯的高速漩渦出現(xiàn),這是因為在離心力作用下流體向管道外側(cè)運動,并逐漸外側(cè)堆積,進而出現(xiàn)外側(cè)向內(nèi)側(cè)的流動。流體離開轉(zhuǎn)彎后,彎管外側(cè)的速度較大,而內(nèi)側(cè)速度較小。

3)T型管

T型管形狀較上兩種都復(fù)雜,圖5所示為T型管縱截面壓力分布云圖,左側(cè)為冰漿進口,流體在三叉處分流至兩個出口,其中上側(cè)為出口1,右側(cè)為出口2。

圖4 90°彎管的管道縱截面冰漿速度分布云圖Fig.4 The contour of velocity for the ice crystals on the longitudinal section of the tubes

圖5 T型管縱截面壓力分布云圖Fig.5 The contour of pressure on the longitudinal section of T-type tube

由圖5可以看出,沿著冰漿流動方向壓力變化趨勢分明,從流動進口到分流之前,壓力逐漸下降。在拐角處分流至兩個方向,一個方向是垂直于來流方向,另一個方向為沿著原速度方向向前。拐角分流處液體壓力劇烈變化,可以看出冰漿流經(jīng)此處時有劇烈的湍流現(xiàn)象和能量交換,經(jīng)過動量交換,冰漿在沿著來流方向的分流方向上壓力劇增,然后又逐漸下降;而在垂直于來流方向的分流方向上,壓力基本保持不變。值得注意的是垂直分流方向上的兩個拐角出現(xiàn)一個壓力最高點和最低點,最低點的出現(xiàn)是由于邊界層脫落分離,造成壓力突然降低;最高點的出現(xiàn)與90°彎管外側(cè)的情況相似,都是由于能量交換使速度減小而壓力增大。

圖6所示為T型管縱截面冰漿平均速度分布云圖,冰漿的平均速度除了在入口處有明顯的變化外,在分流前基本保持不變,在分流處急劇下降至一定值,兩個分流管內(nèi)的流體平均速度基本相等。

圖6 T型管縱截面冰漿平均速度分布云圖Fig.6 The contour of average velocity for ice slurry

3 實驗方法及過程

參考H. Kumano等[19-20]的實驗方案設(shè)計如圖7所示,實驗系統(tǒng)包括4部分:冰漿存儲裝置、冰漿輸送系統(tǒng)、流動實驗測試段、數(shù)據(jù)采集系統(tǒng)。其中,冰漿存儲裝置是兩個尺寸為0.395 m×0.395 m×0.025 m的不銹鋼箱,外部包有10 mm厚的保溫材料。冰漿輸送系統(tǒng)由冰漿輸送泵、過濾器、流量調(diào)節(jié)裝置和輸送管構(gòu)成,冰漿輸送泵采用可用于漿體的磁力驅(qū)動循環(huán)泵,過濾裝置用于篩除粒徑大于1×10-3m的冰粒子。流動實驗段需要測量溫度、壓降和流量等物理量,實驗測量裝置的參數(shù)見表2。數(shù)據(jù)采集主要由GP10便攜式無紙記錄儀記錄。測量段的3種管型采用的是管徑為2×10-2m的黃銅管,管總長為1 m,且彎管和T型管均在中間位置有改變,彎管曲率半徑為0.1 m,測試管水平放置,外包裹厚度約為1×10-2m的聚氨酯材料進行保溫。

圖7 流動實驗段原理Fig.7 Principle of the flow experiments

表2 實驗測量裝置參數(shù)

實驗前,先對儀表進行標(biāo)定,溫度校準(zhǔn)設(shè)備為OPTI-CAL溫濕度校驗儀,精度±0.1 ℃;壓力變送器經(jīng)過Beamex多功能校驗儀標(biāo)定,INT160壓力模塊的量程為0~16 MPa,分辨率為0.000 1,精度為0.005%FS+0.0125%RDG。

實驗時,儲存箱1中的冰漿溫度在-2 ℃左右,冰粒直徑約為4×10-4~8×10-4m。IPF為10%~20%,流速在0.3~1.5 m/s選定5組。

冰漿的含冰率IPF用密度測量法計算為:

(16)

式中:ρslurry為溶液平均密度(質(zhì)量與體積的比值,體積用量筒測量),kg/m3;ρice、ρwater分別為冰和水的密度,kg/m3。

4 實驗結(jié)果分析

1)管型對壓降的影響

管道形狀分為直管、90°彎管和T型管,3者之間對比分析的意義在于研究阻力部位的形狀和阻力程度的關(guān)系。

圖8 不同含冰率下3種管型內(nèi)的壓降隨流速變化Fig.8 Variation trend of pressure drop in three piping shapes

由圖8可以看出,壓降受管型影響較大。含冰率為10%和20%,壓降值從大到小依次為:T型管>彎管>直管。流速較低時,3種管道的壓降值較接近,相差不大。流速越高,由于管型阻力不同,造成同流速下T型管與90°彎管、90°彎管與直管間壓降差值變大。直管中沒有阻力部件,因此壓降主要是沿程損失,且直管部分的曲線可看出速度越高,壓力損失越大。彎管中有一個彎曲部位,流動方向改變時,發(fā)生二次流,冰漿會撞擊外側(cè)壁面造成動量的損失,在豎直部分流動時,外側(cè)流速大,質(zhì)量多,故摩擦劇烈,壓力損失較大。T型管分流后一支路沿原速度方向,另一支路轉(zhuǎn)向90°,冰漿流經(jīng)分流處遇到岔口,會形成漩渦,擾動更大。因此,在冰漿輸送中,應(yīng)盡量使用直管和彎管,T型管壓降過大,盡量減少使用。

圖9 直管壓力損失系數(shù)與流速的關(guān)系Fig.9 Pressure lose coefficients in horizontal piping

流體流動時會產(chǎn)生沿程阻力,為了評價速度對壓力損失的影響,參考流體力學(xué)中沿程損失和文獻[21],研究壓力損失系數(shù):

(17)

圖9所示為直管壓力損失系數(shù)與流速的關(guān)系,可以看出,水平直管中流速越大,壓力損失系數(shù)越小,即對壓降的影響越小。而流速與雷諾數(shù)成正比,且密度和黏度變化不大時,壓力損失系數(shù)與雷諾數(shù)的關(guān)系與圖9相同。

2)含冰率IPF對壓降的影響

如圖8所示,3種管型,IPF為20%的壓降比IPF為10%的大。這是由于IPF越大,冰漿的非牛頓流體特征越明顯,黏度也隨固體顆粒的增多而增大。3種管型在高流速時的壓降差高于低流速。這是因為速度小,冰水分層現(xiàn)象較為明顯;隨著速度增大,冰水混合逐漸均勻,顆粒效應(yīng)使得黏度變大,因此壓降差增大。當(dāng)流速為0.9 m/s時,IPF從10%增大到20%時,直管壓降從1.1 kPa/m增大到1.5 kPa/m; 90°彎管壓降從1.9 kPa/m增大到2.1 kPa/m; T型管壓降從2.27 kPa/m增大到2.73 kPa/m。由圖8還可以看出,當(dāng)速度增大到一定值時,彎管和T型管壓降差值趨于恒定,原因是阻力部件存在,使得冰水摻混劇烈,易達到充分混合狀態(tài),此時,顆粒效應(yīng)穩(wěn)定,黏度差恒定,壓降差基本不變。

3)實驗與模擬結(jié)果對比

圖10 壓降模擬值與實驗值對比Fig.10 Comparison of simulation results and experimental data for the pressure drop

圖10所示為水平直管中實驗數(shù)據(jù)與模擬數(shù)據(jù)的對比。由圖可知數(shù)值模擬和實驗的壓降隨流速的變化趨勢一致,但在流速較大時,模擬值要大于實驗值且誤差達到20%,這是由于流速越大,冰漿流動的湍流程度就越大,冰顆粒與水的相互碰撞越劇烈,二者之間能量交換越多,相變影響越大,而數(shù)值模擬不考慮相變,所以模擬值與實驗值偏差較大。

5 結(jié)論

本文從模擬和實驗兩個角度對冰漿在帶有阻力部件的管道中流動特性進行了研究。管道的形狀有直管、90°彎管和T型管3種,其中直管不含局部阻力部件,彎管帶有一個彎曲處,T型管的連接部位是三通,帶有兩個彎曲處。主要研究了壓降受含冰率、流速、管型的影響,結(jié)論如下:

1)從數(shù)值模擬的結(jié)果看,不考慮相變的影響時,3種管型在入口處都有明顯的壓降。對于彎管,在拐彎處出現(xiàn)二次流,壓力變化較大;而T型管在三通處壓力變化明顯,甚至在垂直分流方向上的兩個拐角分別出現(xiàn)壓力最高點和最低點。

2)當(dāng)實驗條件和含冰率相同時,直管的壓降最小,90°彎管次之,而T型管的壓降最為嚴(yán)重。因此,在實際應(yīng)用過程中應(yīng)盡量避免T型管的使用。

3)冰漿各管型中水平方向流動時,壓降都隨著含冰率和流速的增大而增大。

4)水平直管中計算的沿程壓力損失系數(shù)隨速度的增大而降低。

5)直管的流動模擬結(jié)果與實驗值誤差在20%以內(nèi),在低流速時吻合較好。

本文受天津市應(yīng)用基礎(chǔ)與前沿技術(shù)研究計劃項目(15JCYBJC21600)和天津市高等學(xué)校創(chuàng)新團隊項目(TD12-5048)資助。(The project was supported by the Tianjin Research Program of Application Foundation and Advanced Technology (No.15JCYBJC21600) and the Science Research Innovation Team Project in Tianjin (No.TD12-5048).)

[1] EGOLF P W, KAUFFELD M. From physical properties of ice slurries to industrial ice slurry applications[J]. International Journal of Refrigeration, 2005, 28(1):4-12.

[2] DAVIES T W. Slurry ice as a heat transfer fluid with a large number of application domains[J]. International Journal of Refrigeration, 2005, 28(1):108-114.

[3] KAUFFELD M, WANG M J, GOLDSTEIN V, et al. Ice slurry applications[J]. International Journal of Refrigeration, 2010, 33(8):1491.

[4] BELLAS I, TASSOU S A. Present and future applications of ice slurries[J]. International Journal of Refrigeration, 2005, 28(1):115-121.

[5] 杜衛(wèi)新, 鄭光相. 冰漿制冷降溫系統(tǒng)在平煤六礦的應(yīng)用[J]. 煤炭工程, 2009, 365(4):63-65. (DU Weixin, ZHENG Guangxiang. The application of the ice slurry system for six coal mine in Pingdingshan Coal Group Co. Ltd.[J]. Coal Engineering, 2009, 365(4):63-65.)

[6] OTAKE H, SHITE J, PAREDES O L, et al. Catheter-based transcoronary myocardial hypothermia attenuates arrhythmia and myocardial necrosis in pigs with acute myocardial infarction[J]. Journal of the American College of Cardiology, 2007, 49(2):250-260.

[7] 何國庚, 王忠衡. 冰漿流體流動與換熱研究綜述[J]. 制冷學(xué)報, 2005, 26(4):1-5. (HE Guogeng, WANG Zhongheng. Review of study on flow and heat transfer of ice slurry[J]. Journal of Refrigeration, 2005, 26(4):1-5.)

[8] GROZDEK M, KHODABANDEH R, LUNDQVIST P. Experimental investigation of ice slurry flow pressure drop in horizontal tubes[J]. International Journal of Refrigeration, 2010, 33(2):357-370.

[10] AYEL V, LOTTIN O, PEERHOSSAINI H. Rheology, flow behaviour and heat transfer of ice slurries: a review of the state of the art[J]. International Journal of Refrigeration, 2003, 26(1):95-107.

[11] KITANOVSKI A, VUARNOZ D, ATA-CAESAR D, et al. The fluid dynamics of ice slurry[J]. International Journal of Refrigeration, 2005, 28(1):37-50.

[12] KUO H P, KNIGHT P C, PARKER D J, et al. The influence of DEM simulation parameters on the particle behaviour in a V-mixer[J]. Chemical Engineering Science, 2002, 57(17):3621-3638.

[13] GIDASPOW D. Multiphase flow and fluidization[J]. Multiphase Flow & Fluidization, 1994(95):1-29.

[14] MONTEIRO A C S, BANSAL P K. Pressure drop characteristics and rheological modeling of ice slurry flow in pipes[J]. International Journal of Refrigeration, 2010, 33(8):1523-1532.

[15] NIEZGODA-ZELASKO B, ZALEWSKI W. Momentum transfer of ice slurry flows in tubes, experimental investigations[J]. International Journal of Refrigeration, 2006, 29(3):418-428.

[16] 王繼紅, 王樹剛, 張騰飛, 等. 水平管道內(nèi)冰漿流體阻力特性CFD模擬[J]. 大連理工大學(xué)學(xué)報, 2012, 52(6):845-849. (WANG Jihong, WANG Shugang, ZHANG Tengfei, et al. Simulation of pressure drop for ice slurry flow in horizontal pipes by CFD[J]. Journal of Dalian University of Technology, 2012, 52(6):845-849.)

[17] THOMAS D G. Transport characteristics of suspension: VIII. A note on the viscosity of Newtonian suspensions of uniform spherical particles[J]. Journal of Colloid Science, 1965, 20(3):267-277.

[18] 王繼紅, 王樹剛, 張騰飛, 等. 水平90°彎管內(nèi)冰漿流體流動特性的數(shù)值模擬[J]. 高校化學(xué)工程學(xué)報, 2012, 26(4):581-586. (WANG Jihong, WANG Shugang, ZHANG Tengfei, et al. Numerical simulation of ice slurry flow in a horizontal 90° elbow pipe[J]. Journal of Chemical Engineering of Chinese Universities, 2012, 26(4):581-586.)

[19] KUMANO H, HIRATA T, SHIRAKAWA M, et al. Flow characteristics of ice slurry in narrow tubes[J]. International Journal of Refrigeration, 2010, 33(8): 1513-1522.

[20] MIKA. Ice slurry flow in a poppet-type flow control valve[J]. Experimental Thermal & Fluid Science, 2013, 45:128-135.

[21] MIKA. Pressure loss coefficients of ice slurry in horizontally installed flow dividers[J]. Experimental Thermal and Fluid Science, 2013, 45(2): 249-258.

猜你喜歡
實驗模型
一半模型
記一次有趣的實驗
微型實驗里看“燃燒”
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
做個怪怪長實驗
3D打印中的模型分割與打包
NO與NO2相互轉(zhuǎn)化實驗的改進
實踐十號上的19項實驗
太空探索(2016年5期)2016-07-12 15:17:55
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
主站蜘蛛池模板: 国产极品美女在线播放| 精品无码国产一区二区三区AV| 五月婷婷亚洲综合| 思思热精品在线8| 亚洲综合香蕉| 国产精品福利在线观看无码卡| 国产网站在线看| 456亚洲人成高清在线| 国产免费a级片| 一本大道AV人久久综合| 狠狠色狠狠色综合久久第一次| 丁香综合在线| 波多野结衣二区| 欧美日韩在线成人| 国产成人精品男人的天堂| 色综合天天娱乐综合网| 国产主播在线一区| 亚洲国产91人成在线| 在线精品亚洲一区二区古装| 丁香五月婷婷激情基地| 国产成人无码Av在线播放无广告| 国产丰满成熟女性性满足视频| 91青青视频| 91福利免费| 国内精自视频品线一二区| 制服丝袜无码每日更新| 玖玖免费视频在线观看| 亚洲第一区在线| 亚洲欧美综合精品久久成人网| 色香蕉影院| 国产亚洲第一页| 成人91在线| 一本综合久久| 91精品视频网站| 日韩欧美高清视频| 天堂久久久久久中文字幕| 伊人五月丁香综合AⅤ| 久久综合色天堂av| 国产成人精品一区二区秒拍1o| 国产成人亚洲精品蜜芽影院| 久久不卡精品| 99热在线只有精品| 国产真实乱了在线播放| 亚洲无线一二三四区男男| 国产亚洲美日韩AV中文字幕无码成人 | 成人国产免费| 99久久免费精品特色大片| 国产成人h在线观看网站站| 亚洲一区二区黄色| 亚洲欧美日韩成人高清在线一区| 麻豆国产原创视频在线播放| 欧美色综合久久| 国产欧美在线视频免费| 日本精品视频一区二区| 久久精品国产国语对白| 国产黄色片在线看| 制服丝袜亚洲| 国产日本一线在线观看免费| 精品欧美一区二区三区久久久| 欧美综合中文字幕久久| 美女无遮挡免费视频网站| 亚洲婷婷六月| 久久国产精品波多野结衣| 久久精品视频亚洲| 中文字幕在线欧美| 免费一级全黄少妇性色生活片| 久久人人妻人人爽人人卡片av| 亚洲欧美人成电影在线观看| 久久久久免费看成人影片| 欧美激情视频一区| 一级毛片免费不卡在线视频| 国产第一页屁屁影院| 国产精品蜜芽在线观看| 青青操视频免费观看| 亚洲精品少妇熟女| 国产精品美女网站| 久久综合九色综合97婷婷| 亚洲bt欧美bt精品| 亚洲成人黄色在线观看| 国产精品内射视频| 国产91色| 一本久道久久综合多人|