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

超臨界CO2熱泵蛇形管氣冷器傳熱特性研究

2023-11-09 01:37:44崔海亭
河北科技大學學報 2023年5期
關鍵詞:質(zhì)量

崔海亭,王 超,王 晨

(河北科技大學機械工程學院,河北石家莊 050018)

在“碳達峰、碳中和”的目標背景下,節(jié)能減排是一項重要舉措。在制冷行業(yè),傳統(tǒng)的CHCs和HCHCs類制冷劑仍會產(chǎn)生溫室氣體,使臭氧層遭受破壞。因此,尋找節(jié)能環(huán)保的制冷劑對目標的實現(xiàn)意義重大[1-2]。CO2作為一種自然工質(zhì),因其具有無毒、不燃、安全、成本低、環(huán)境友好等特點[3],引起了制冷、熱泵等領域的研究人員和工業(yè)界的關注。氣體冷卻器作為CO2熱泵系統(tǒng)中的關鍵部件,決定著整個系統(tǒng)的運行效率,超臨界CO2在氣冷器中放熱有著較大的溫度滑移,可以將冷卻水加熱至較高的溫度。最近有研究表明,CO2跨臨界系統(tǒng)有著較高的節(jié)能潛力[4-8]。然而,由于CO2在臨界點處急劇變化的物理性質(zhì),超臨界壓力下CO2換熱比跨臨界壓力下?lián)Q熱機理更為復雜,因此有必要深入研究超臨界CO2在換熱管中的換熱特性。

超臨界CO2在管內(nèi)的流動和傳熱過程受到廣大學者的廣泛研究[9-12],主要是超臨界CO2在直管和螺旋管流動的研究。由于在臨界點處CO2密度發(fā)生突變,產(chǎn)生浮升力和流動加速效應,從而對超臨界CO2的流動換熱產(chǎn)生影響。當CO2在管內(nèi)流動時,不同的操作參數(shù)會對管內(nèi)工質(zhì)的傳熱特性產(chǎn)生影響。顧騫[13]對不同操作參數(shù)下超臨界CO2在直管和螺旋管的對流傳熱過程進行了數(shù)值模擬,模擬結(jié)果表明在直管中,擬臨界點處浮升力對傳熱的影響不可忽略,增大入口雷諾數(shù)會強化換熱,增大管內(nèi)壓力會使傳熱系數(shù)峰值減小;在螺旋管中,浮升力的影響在流體達到臨界點前不可忽略。白萬金等[14]對不同操作參數(shù)下超臨界CO2在水平直管的流動傳熱進行實驗研究,結(jié)果表明質(zhì)量流速越大,傳熱系數(shù)越大;壓力越大,傳熱系數(shù)峰值點所處的溫度越高。崔海亭等[15]通過數(shù)值模擬對不同操作參數(shù)下超臨界CO2在直管套管內(nèi)的熵產(chǎn)進行了計算,結(jié)果表明壓力和熵產(chǎn)成正比,質(zhì)量流量和熵產(chǎn)成反比。

目前,對超臨界CO2在蛇形管[16-18]內(nèi)流動換熱的研究比較有限。羅峰等[19]對超臨界CO2在微細蛇形管內(nèi)層流對流換熱系數(shù)的影響進行了研究,對比分析了微細蛇形管豎直和水平對CO2流動換熱的影響,結(jié)果表明,當重力和入口流動方向相反時,截面處溫度和速度對稱分布,否則對稱分布效果減弱。黃騰等[20]對不同結(jié)構(gòu)蛇形管內(nèi)的超臨界CO2進行了數(shù)值模擬,結(jié)果表明曲率直徑或內(nèi)徑的增加均會降低傳熱系數(shù)。對于超臨界CO2對蛇形管傳熱的影響的研究多針對結(jié)構(gòu)參數(shù),而對不同操作參數(shù)下的蛇形管換熱研究較少。

本文針對套管式蛇形管氣體冷卻器建立物理模型,通過改變操作參數(shù)對管內(nèi)超臨界CO2的流動和換熱特性進行研究,具體分析不同壓力、不同質(zhì)量流量下超臨界CO2的傳熱性能。

1 模型與求解

1.1 物理模型

蛇形管結(jié)構(gòu)示意圖如圖1所示,該氣體冷卻器采用套管式。蛇形管管長L=1 200 mm,曲率直徑D=109.20 mm,外管內(nèi)徑d1=14 mm,材料為不銹鋼,內(nèi)管為管徑d2=3.80 mm、壁厚δ=1.1 mm的不銹鋼管。物理模型由3部分區(qū)域組成:CO2流動區(qū)域、內(nèi)管壁厚區(qū)域、冷卻水流動區(qū)域。氣冷器內(nèi)管工質(zhì)為CO2流體,內(nèi)管和外管的空腔為冷卻水,在本文中二者呈逆流換熱。為了簡化模型,作出如下假設:

圖1 蛇形管氣冷器物理模型Fig.1 Physical model of serpentine tube air cooler

1)假設模型為絕熱系統(tǒng),不與外界進行能量交換;

2)忽略套管外管壁厚的影響。

1.2 邊界條件與網(wǎng)格劃分

本文模擬采用壓力求解器,壓力-速度耦合采用SIMPLEC算法,壓力插值格式為PRESTIO!,動量、能量、湍動能等采用二階快速Q(mào)UICK格式,設置收斂殘差為10-6,當滿足殘差要求且進出口流量守恒時,計算完成。管內(nèi)介質(zhì)為超臨界CO2,物性參數(shù)通過REFPROP獲得,在Fluent軟件中通過piecewise-liner輸入變物性數(shù)據(jù),節(jié)點數(shù)取8個。采用RNGk-ε模型和增強壁面函數(shù),冷卻水和CO2入口均為質(zhì)量流量入口,出口為壓力出口,外管外壁面設為絕熱,內(nèi)管內(nèi)壁和外壁設為流固耦合壁面,邊界條件具體參數(shù)如表1所示。

表1 蛇形管氣冷器邊界條件設置Tab.1 Setting of boundary conditions of serpentine tube air cooler

蛇形管網(wǎng)格劃分使用Gambit軟件,網(wǎng)格類型為結(jié)構(gòu)化網(wǎng)格,近壁處由于流體受邊界層存在的影響,會對溫度等特性產(chǎn)生影響,故在內(nèi)管內(nèi)壁和外壁與流體交界處劃分邊界層,從而加強計算的準確性。本模型網(wǎng)格質(zhì)量評價指標Equisize Skew<0.4的網(wǎng)格數(shù)占總網(wǎng)格的99.67%,網(wǎng)格質(zhì)量良好。蛇形管網(wǎng)格劃分示意圖如圖2所示。

圖2 蛇形管網(wǎng)格劃分示意圖Fig.2 Schematic diagram of serpentine tube meshing

1.3 數(shù)學模型

本文采用RNGk-ε模型進行計算,數(shù)學模型包括連續(xù)性方程、動量方程、能量方程、湍動能方程(k方程)和耗散率方程(ε方程),具體公式如下。

連續(xù)性方程:

(1)

動量方程:

(2)

能量方程:

(3)

湍動能方程:

(4)

耗散率方程:

(5)

式中湍流黏度μt定義如下:

(6)

1.4 數(shù)據(jù)處理

為了更好地看出內(nèi)管CO2的流動狀態(tài),設置了重力g,其方向與Z軸正方向一致,且模擬過程所有的參數(shù)均選用國際單位。

蛇形管橫截面上的主流溫度Tf取質(zhì)量加權(quán)平均溫度,公式如下:

(7)

超臨界CO2在管內(nèi)流動的局部傳熱系數(shù)hi和平均傳熱系數(shù)h為

(8)

(9)

式中:u,ρ,Cp,qw,Tw,i和Tf,i分別表示速度、流體密度、定壓比熱容、熱流密度、局部壁溫、截面主流溫度。

2 模型驗證

網(wǎng)格無關性驗證數(shù)值模擬驗證結(jié)果如圖3所示。在CO2和冷卻水的質(zhì)量流量分別為0.004和0.030 kg/s,CO2壓力為8 MPa的工況下,對模型進行網(wǎng)格無關性驗證。圖3 a)展示了網(wǎng)格數(shù)分別在465 124,1 411 766和3 192 000下,沿管長方向CO2主流溫度的變化。由圖3 a)可以看出,CO2出口溫度分別為306.96,305.75和305.69 K,網(wǎng)格數(shù)為1 411 766和3 192 000時,CO2出口溫度基本不變,且曲線基本保持一致。為加快計算速度和保證計算精度,選取網(wǎng)格數(shù)為141 766進行數(shù)值模擬。

圖3 無關性驗證和數(shù)值模型驗證Fig.3 Verification of independence and validation of numerical models

為驗證數(shù)值模型和計算方法的準確性,采用文獻[21]中的實驗條件,用等比例方法建立物理模型。參考文獻中蛇形管內(nèi)徑為0.953 mm,外徑為2.1 mm,管長L為88 mm;文獻中對蛇形管進行加熱,故模擬時設置恒熱流壁面條件。以CO2壓力為7.65 MPa,CO2質(zhì)量流量為1 kg/h,壁面熱流為50 kW/m2工況為例進行數(shù)值模型驗證。圖3 b)為沿管長方向CO2主流溫度和壁面溫度變化情況,由圖3 b)可以看出,模擬值和實驗值變化趨勢基本一致,且CO2主流溫度模擬值較實驗值平均誤差為0.69%,壁面溫度模擬值較實驗值平均誤差為1.59%,誤差均在允許范圍內(nèi)。因此,本數(shù)值模型具有可靠性。

3 結(jié)果與分析

本文采用Fluent軟件對蛇形管進行仿真模擬,由于超臨界CO2物性參數(shù)在臨界點附近變化較大,故本文只對內(nèi)管即CO2流動區(qū)域進行分析討論。冷卻水入口溫度為290.15 K,CO2入口溫度為333.15 K。

3.1 蛇形管內(nèi)部流動特性分析

圖4和圖5分別為CO2入口質(zhì)量流量為0.004 kg/s,入口壓力為8 MPa,冷卻水入口質(zhì)量流量為0.03 kg/s時,在不同截面處的速度云圖和溫度云圖。由圖4、圖5可知,CO2沿著流動方向受到二次流的影響,不同截面處的速度云圖和溫度云圖發(fā)生了明顯的變化。最初CO2在云圖上呈現(xiàn)一對渦流,且在2L/7處之前云圖呈現(xiàn)中間向兩側(cè)擴散的趨勢,這表明在入口處由于主流體和壁面溫差較大,在重力和離心力的作用下,使得湍流程度較大;沿管長方向,渦流數(shù)量逐漸減小,湍流程度減弱,且云圖呈現(xiàn)由內(nèi)側(cè)向外側(cè)和由外側(cè)向內(nèi)側(cè)交替擴散的趨勢,從截面流動方向也能看出這種變化,這是由于蛇形管的結(jié)構(gòu)所致,蛇形管周期性的彎曲反向使得離心力也周期反向,在離心力和重力的作用下,CO2向密度較小的一側(cè)擴散。由速度云圖可知,沿著管長方向CO2流速逐漸降低,且流速在入口處變化較大,之后變化較小,在2L/7處較L/7處減小了0.24 m/s,在5L/7處較4L/7處減小了0.17 m/s。由溫度云圖可知,沿著管長方向CO2主流溫度逐漸降低,且溫度呈現(xiàn)先快速下降再緩慢下降的趨勢,入口處由于溫差較大,溫度下降較快,到x=L/7處溫度已經(jīng)從333.15 K降低到321.63 K,此時溫度已經(jīng)降低了3.46%;在x=3L/7之后由于溫差和湍流強度的減小,導致截面處二次流減弱,因此截面處溫度云圖變化不再顯著。

圖4 蛇形管不同截面速度云圖和流線圖Fig.4 Cloud and streamline diagrams of velocity on different sections of serpentine tube

圖5 蛇形管不同截面溫度云圖Fig.5 Cloud chart of temperature on different sections of serpentine tube

圖6為沿管長方向湍動能的變化。由圖6可知,湍動能在入口附近先增大再減小,這是因為在入口附近溫差較大,同時由于蛇形管連續(xù)彎曲的結(jié)構(gòu)導致離心力周期反向,從而湍流程度增大。由此可得,蛇形管周期性的彎曲反向使得離心力也周期反向,從而云圖呈現(xiàn)周期性的內(nèi)側(cè)和外側(cè)交互擴散的趨勢,強化了換熱效果。

圖6 沿管長方向湍動能的變化Fig.6 Changes in turbulent kinetic energy along a rectangular tube

3.2 CO2壓力對傳熱特性的影響

為考察CO2壓力變化對蛇形管內(nèi)超臨界CO2流動換熱的影響,設置CO2入口質(zhì)量流量為0.004 kg/s,冷卻水入口質(zhì)量流量為0.03 kg/s,分別改變CO2入口壓力為8,9和10 MPa,對蛇形管進行傳熱模擬。CO2在不同壓力下沿管長方向主流溫度變化和傳熱系數(shù)變化如圖7和圖8所示。由圖7和圖8可以看出,不同壓力下,沿管長方向CO2主流溫度均呈下降趨勢,傳熱系數(shù)均呈先上升再下降趨勢。

圖7 不同壓力下沿管長方向主流溫度變化Fig.7 Temperature variation of the mainstream along the length direction of the pipe under different pressures

圖8 不同壓力下沿管長方向傳熱系數(shù)變化Fig.8 Change of heat transfer coefficient along the tube length under different pressures

如圖7所示,沿管長方向CO2主流溫度先快速降低再緩慢降低,壓力越低,主流溫度快速下降得越快,緩慢下降得越慢。入口處由于較大的溫差,主流溫度快速下降,壓力越接近臨界壓力,湍流動能越大,因此,低壓下主流溫度快速下降的越快,壓力越低,主流溫度緩慢下降的越慢,這可以減小和壁面之間的溫差,使得CO2充分換熱。同時,CO2進出口溫差也不同,壓力越小,進出口溫差越小。當壓力為8 MPa時,出口溫度為305.76 K,進出口溫差為26.12 K;當壓力為9 MPa時,出口溫度為304.38 K,進出口溫差為27.64 K;當壓力為10 MPa時,出口溫度為301.69 K,進出口溫差為30.50 K。因此,減小壓力可以使CO2充分換熱。

如圖8所示,在不同壓力下,沿著管長方向傳熱系數(shù)先增大,到臨界點達到峰值,隨后再減小。壓力越大,峰值越小,這是因為在臨界點處,CO2比熱容的急劇變化導致?lián)Q熱系數(shù)的變化,同時壓力越高,擬臨界溫度就越高,導致到達峰值點的位置提前。不同壓力下超臨界CO2的擬臨界溫度如表2所示。

表2 不同壓力下的擬臨界溫度值Tab.2 Critical temperature values under different pressures

不同CO2溫度狀態(tài)下的傳熱系數(shù)明顯不同,在CO2溫度較高時,壓力越大,傳熱系數(shù)越高;相反,在CO2溫度較低時,壓力越大,傳熱系數(shù)越低。這是由CO2變物性所致,在臨界點之前,CO2溫度較高,隨著CO2壓力的升高,定壓比熱容和熱導率逐漸升高,在臨界點之后則相反。當P=8 MPa,時,峰值點傳熱系數(shù)為10 855.89 W/(m2·K),平均傳熱系數(shù)為6 135.89 W/(m2·K),較P=9 MPa時分別提高了63.59%和24.37%;較P=10 MPa時分別提高了98.76%和42.53%。因此,壓力越低,越靠近準臨界點,換熱系數(shù)越大,越有利于換熱。

3.3 CO2質(zhì)量流量對傳熱特性的影響

為了解CO2質(zhì)量流量變化對蛇形管內(nèi)超臨界CO2流動換熱的影響,設置CO2入口壓力為8 MPa,冷卻水入口質(zhì)量流量為0.03 kg/s,分別改變CO2入口質(zhì)量流量為0.003,0.004和0.005 kg/s對蛇形管進行傳熱模擬。不同CO2質(zhì)量流量下沿管長方向的溫度分布和傳熱系數(shù)變化分別如圖9和圖10所示。

圖9 不同CO2質(zhì)量流量下沿管長方向溫度分布圖Fig.9 Temperature distribution along the tube length under different CO2 mass flow rates

圖10 不同CO2質(zhì)量流量下沿管長方向傳熱系數(shù)變化Fig.10 Change of heat transfer coefficient along the tube length under different CO2 mass flow rates

由圖9可知,在不同質(zhì)量流量下主流溫度和壁面溫度隨著管長變化而變化的趨勢基本相同,都呈下降趨勢。但CO2的進出口溫差以及CO2和壁面的溫差卻有所不同,當CO2質(zhì)量流量為0.003 kg/s時,出口溫度為302.41 K,進出口溫差為29.18 K;當CO2質(zhì)量流量為0.004 kg/s時,出口溫度為305.76 K,進出口溫差為26.12 K;當CO2質(zhì)量流量為0.004 kg/s時,出口溫度為307.12 K,進出口溫差為24.97 K,質(zhì)量流量越大,進出口溫差越小。CO2和壁面的溫差變化也是如此,在流量分別為0.003,0.004,和0.005 kg/s時,最小溫差分別為5.16,4.28和3.81 K。由此可得,隨著CO2質(zhì)量流量的增加,超臨界CO2冷卻換熱具有更小的溫差,具有更小的溫度滑移,可以充分換熱。

由圖10中傳熱系數(shù)模擬值可知,不同CO2質(zhì)量流量下傳熱系數(shù)的變化趨勢一致,沿管長方向先增加,在臨界點達到達峰值后再下降。這是由于CO2的變物性所致,在臨界點處比熱和熱導率最大。質(zhì)量流量越低,到達峰值點的位置越靠前,這是因為低質(zhì)量流量下的CO2主流溫度會較早地降低到擬臨界溫度點。CO2質(zhì)量流量越大,峰值點越大,CO2質(zhì)量流量為0.005 kg/s的峰值點較0.003和0.004 kg/s分別提高了54.18%和16.88%,質(zhì)量流量為0.005 kg/s的平均傳熱系數(shù)較0.003和0.004 kg/s分別提高了57.92%和19.83%。由此可得,隨著CO2質(zhì)量流量的增加,湍流強度加強,邊界層變薄,導致傳熱得到改善,平均傳熱系數(shù)提高。

3.4 冷卻水質(zhì)量流量對傳熱特性的影響

為了解冷卻水質(zhì)量流量變化對蛇形管內(nèi)超臨界CO2流動換熱的影響,設置CO2入口壓力為9 MPa,質(zhì)量流量為0.004 kg/s,分別改變冷卻水質(zhì)量流量為0.02,0.03和0.04 kg/s對蛇形管的傳熱性進行模擬。不同冷卻水質(zhì)量流量下沿管長方向傳熱系數(shù)的變化,如圖11所示。從圖11可以看出,沿管長方向不同冷卻水質(zhì)量流量下的傳熱系數(shù)都是先增加后減小,在臨界點附近達到峰值,且峰值點傳熱系數(shù)相差不大。在CO2溫度較高時,冷卻水的質(zhì)量流量越大,傳熱系數(shù)越高;相反,傳熱系數(shù)越低。這是由于CO2流體的變物性所致,冷卻水質(zhì)量流量的變化引起CO2主流溫度的變化,從而影響CO2的定壓比熱容和熱導率。

圖11 不同冷卻水質(zhì)量流量下沿管長方向傳熱系數(shù)變化Fig.11 Change of heat transfer coefficient along tube length under different cooling water mass flow

冷卻水質(zhì)量流量越大,到達峰值點的位置越靠前,冷卻水質(zhì)量流量分別為0.02,0.03和0.04 kg/s時的平均傳熱系數(shù)分別為5 029.89,4 933.66 和4 821.84 W/(m2·K)。由此可知,冷卻水質(zhì)量流量的增加不會影響傳熱系數(shù)的峰值點,但會使峰值點所在的位置提前,質(zhì)量流量越大,平均傳熱系數(shù)越低。

4 結(jié) 論

1) 通過對蛇形管內(nèi)部CO2流動特性分析可知,由于離心力的周期性反向,在重力的作用下,導致溫度和速度梯度呈現(xiàn)出內(nèi)側(cè)和外側(cè)周期性交互擴散的變化趨勢,從而強化換熱。在入口處由于溫差較大,蛇形結(jié)構(gòu)使得溫度和速度云圖均產(chǎn)生多個渦,湍流強度增大。

2) 相同工況下,CO2壓力越接近臨界點,平均傳熱系數(shù)越大,壓力為8 MPa時的平均傳熱系數(shù)較為9 MPa和10 MPa分別提高了24.37%和42.53%。由于受超臨界CO2物性參數(shù)的影響,在CO2溫度較高時,壓力越大,傳熱系數(shù)越高;相反,在CO2溫度較低時,壓力越大,傳熱系數(shù)越低。

3) 相同工況下,CO2質(zhì)量流量越大,平均傳熱系數(shù)越高,這是由于隨著超臨界CO2質(zhì)量流量的增加,邊界層厚度不斷減薄,湍流更加劇烈,換熱更充分。

4) 相同工況下,冷卻水質(zhì)量流量越大,平均傳熱系數(shù)越低。冷卻水質(zhì)量流量的增加對峰值點的傳熱系數(shù)沒有影響,但會使峰值點出現(xiàn)的位置提前。

本文采用數(shù)值模擬的方法對超臨界CO2在蛇形管中的冷卻過程進行了分析,通過改變操作參數(shù),分析其換熱性能。未來需進一步研究結(jié)構(gòu)參數(shù)變化對超臨界CO2在蛇形管中換熱特性的影響,以綜合評價蛇形管的傳熱性能。

猜你喜歡
質(zhì)量
聚焦質(zhì)量守恒定律
“質(zhì)量”知識鞏固
“質(zhì)量”知識鞏固
質(zhì)量守恒定律考什么
做夢導致睡眠質(zhì)量差嗎
焊接質(zhì)量的控制
關于質(zhì)量的快速Q(mào)&A
初中『質(zhì)量』點擊
質(zhì)量投訴超六成
汽車觀察(2016年3期)2016-02-28 13:16:26
你睡得香嗎?
民生周刊(2014年7期)2014-03-28 01:30:54
主站蜘蛛池模板: 国产精品女人呻吟在线观看| 色综合久久无码网| 欧美日本激情| 国产高清免费午夜在线视频| 亚洲最黄视频| 中文字幕第4页| 婷婷综合缴情亚洲五月伊| 韩国福利一区| 免费激情网址| 思思热精品在线8| 国产尹人香蕉综合在线电影| 网友自拍视频精品区| 永久在线精品免费视频观看| 亚洲最新地址| 高潮毛片免费观看| 国产成人毛片| 波多野结衣在线一区二区| 久久综合色88| 青青操视频免费观看| 亚洲va在线观看| 久久公开视频| 国产视频你懂得| 婷婷色中文网| 国产爽爽视频| 亚洲国产天堂在线观看| 欧美日韩在线国产| 亚洲国产一区在线观看| 国内精品伊人久久久久7777人| 亚洲综合专区| a在线亚洲男人的天堂试看| 中文字幕 91| 久久精品视频亚洲| 亚洲成人www| 99久久精品国产综合婷婷| 国产亚洲欧美日本一二三本道| 日本成人精品视频| 精品亚洲国产成人AV| 欧美一级特黄aaaaaa在线看片| 色婷婷国产精品视频| 国产精品va| 国产精品亚洲五月天高清| 专干老肥熟女视频网站| 91精品啪在线观看国产| 亚洲大尺码专区影院| 广东一级毛片| 亚洲天堂网在线观看视频| 欧美精品啪啪| 欧美 亚洲 日韩 国产| 色欲不卡无码一区二区| 国产资源免费观看| 国产专区综合另类日韩一区| 日韩在线播放中文字幕| av性天堂网| 国产高潮视频在线观看| 亚洲一本大道在线| 国产精品女同一区三区五区| 91久久天天躁狠狠躁夜夜| 丁香五月婷婷激情基地| 国产日本视频91| av免费在线观看美女叉开腿| 白丝美女办公室高潮喷水视频| 真实国产乱子伦高清| 亚洲区欧美区| 久久无码av一区二区三区| 国产素人在线| 国产成人精品视频一区二区电影| 无码啪啪精品天堂浪潮av| 91蝌蚪视频在线观看| 亚洲欧美日韩精品专区| 1769国产精品视频免费观看| 九九热这里只有国产精品| 亚洲人妖在线| 91精品国产福利| 中文字幕在线观看日本| 国产人妖视频一区在线观看| 好吊妞欧美视频免费| 丁香婷婷激情网| 99久久这里只精品麻豆 | 亚洲精品无码AⅤ片青青在线观看| 亚洲欧美色中文字幕| 91蜜芽尤物福利在线观看| 亚洲天堂福利视频|