亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        單回路液氦脈動(dòng)熱管氣液兩相流數(shù)值模擬

        2023-03-20 05:46:18呂秉坤施雅然信紀(jì)軍方志春李來(lái)風(fēng)
        低溫工程 2023年1期
        關(guān)鍵詞:液氦工質(zhì)氣液

        呂秉坤 徐 冬 王 維 施雅然 信紀(jì)軍 方志春 李來(lái)風(fēng)

        (1 松山湖材料實(shí)驗(yàn)室 東莞 523808)

        (2 中國(guó)科學(xué)院低溫工程學(xué)重點(diǎn)實(shí)驗(yàn)室(理化技術(shù)研究所)北京 100190)

        1 引言

        脈動(dòng)熱管(Pulsating Heat Pipe,PHP)是20 世紀(jì)90 年代出現(xiàn)的一種新型熱管,由日本的Akachi 最早提出。它是將一根內(nèi)徑較小(一般為0.5—3 mm)的毛細(xì)管在蒸發(fā)段和冷凝段之間反復(fù)彎折并將首尾連接形成的蛇形結(jié)構(gòu),其余部分為絕熱段。其工作原理是將管內(nèi)抽成真空后充入工作介質(zhì),由于管徑足夠小,在表面張力的作用下,接近飽和態(tài)的兩相流體工質(zhì)在毛細(xì)管內(nèi)以氣塞和液塞的形式隨機(jī)交替分布。當(dāng)管的一端(蒸發(fā)段)被加熱時(shí),工質(zhì)吸熱導(dǎo)致氣泡生長(zhǎng)和蒸氣壓力的增大;與此同時(shí),氣泡在另一端(冷凝段)凝結(jié),氣泡會(huì)收縮破裂,壓力下降。氣泡的生長(zhǎng)和破裂導(dǎo)致蒸發(fā)段和冷凝段存在壓差以及相鄰管子之間存在壓力不平衡,從而推動(dòng)工質(zhì)在管內(nèi)脈動(dòng)運(yùn)動(dòng)或單向循環(huán)流動(dòng),通過(guò)氣液相變的潛熱和液塞流動(dòng)時(shí)的顯熱傳遞熱量。

        脈動(dòng)熱管由于其結(jié)構(gòu)緊湊、性能良好、可靠、靈活、無(wú)噪聲且耐用,經(jīng)過(guò)優(yōu)化還可實(shí)現(xiàn)無(wú)重力運(yùn)行,已經(jīng)成為最具創(chuàng)新性、有效性的兩相換熱系統(tǒng)之一,在微電子領(lǐng)域和工業(yè)領(lǐng)域得到了廣泛應(yīng)用,例如太陽(yáng)能利用、電子設(shè)備冷卻、空間熱管理和余熱回收等。在低溫領(lǐng)域,脈動(dòng)熱管在超導(dǎo)磁體高效冷卻應(yīng)用中也有很大的潛力。

        雖然脈動(dòng)熱管的結(jié)構(gòu)非常簡(jiǎn)單,但由于脈動(dòng)熱管傳熱傳質(zhì)過(guò)程中流體力學(xué)和熱力學(xué)的耦合作用,使得脈動(dòng)熱管的運(yùn)行機(jī)制非常復(fù)雜。對(duì)于低溫脈動(dòng)熱管,循環(huán)流動(dòng)是其高效傳熱的基礎(chǔ)。然而低溫脈動(dòng)熱管的理論研究處于起步階段,低溫脈動(dòng)熱管的可視化到目前為止還沒(méi)有被實(shí)現(xiàn),因此通常采用數(shù)學(xué)建模和數(shù)值模擬兩種方式來(lái)研究低溫脈動(dòng)熱管的運(yùn)行機(jī)理。但是,一維的數(shù)學(xué)模型由于過(guò)多的簡(jiǎn)化與假設(shè),模擬結(jié)果與實(shí)際的兩相流流動(dòng)仍有很大的差別。隨著計(jì)算機(jī)計(jì)算能力的發(fā)展,利用計(jì)算機(jī)進(jìn)行二維或三維的數(shù)值模擬成為主要的研究手段。印度薩爾達(dá)爾瓦拉布巴伊國(guó)家技術(shù)學(xué)院的Sagar 等人[1]采用多相流VOF(Volume of Fluid,VOF)模型對(duì)近臨界直徑二維低溫脈動(dòng)熱管進(jìn)行了數(shù)值模擬,發(fā)現(xiàn)與地面條件相比,在低重力條件下,低溫脈動(dòng)熱管的換熱性能得到了改善。隨后,他們對(duì)另一個(gè)蒸發(fā)段、冷凝段溫度均設(shè)為恒定的二維低溫脈動(dòng)熱管進(jìn)行了數(shù)值模擬[2],研究了不同運(yùn)行條件下液氮的體積分?jǐn)?shù)結(jié)果,對(duì)工質(zhì)溫度、絕熱壁面溫度和流動(dòng)循環(huán)速度進(jìn)行了分析。上海理工大學(xué)的唐愷等人[3-5]利用Fluent 軟件中的多相流VOF 模型,分別研究了氮和氦工質(zhì)低溫脈動(dòng)熱管的流型和充液率、加熱功率、內(nèi)徑和傾斜角等參數(shù)對(duì)性能的影響。哈爾濱理工大學(xué)的徐金柱等[6]對(duì)單環(huán)路液氫溫區(qū)脈動(dòng)熱管進(jìn)行了CFD(Computational Fluid Dynamics)數(shù)值模擬,模擬的傳熱熱阻與實(shí)驗(yàn)值的誤差不超過(guò)15%。

        上述模擬研究中,由于模型的簡(jiǎn)化與假設(shè),模擬結(jié)果與脈動(dòng)熱管內(nèi)部的兩相流仍存在一定程度上的差異,且液氦脈動(dòng)熱管的模擬數(shù)據(jù)十分有限。為了準(zhǔn)確描述液氦脈動(dòng)熱管的氣液兩相流,本研究建立了單回路液氦脈動(dòng)熱管CFD 模型,并對(duì)其初始狀態(tài)和運(yùn)行狀態(tài)的流體流動(dòng)與傳熱進(jìn)行了分析。

        2 液氦脈動(dòng)熱管模型建立

        2.1 幾何模型

        本研究應(yīng)用ANSYS 2020R2 建立液氦脈動(dòng)熱管的模型并求解。采用建模軟件ANSYS DesignModeler軟件繪制二維單回路液氦脈動(dòng)熱管的幾何模型,如圖1所示。液氦脈動(dòng)熱管模型的冷凝段、絕熱段和蒸發(fā)段長(zhǎng)度分別為5 mm、5 mm、5 mm,彎頭直徑為3 mm,彎頭部分也視為絕熱段,內(nèi)徑為0.5 mm,采用豎直放置底部加熱的方式。

        圖1 二維單回路液氦脈動(dòng)熱管幾何模型和網(wǎng)格劃分Fig.1 Geometric model and mesh generation of two-dimensional single loop helium pulsating heat pipe

        采用ANSYS Meshing 軟件對(duì)該脈動(dòng)熱管模型進(jìn)行結(jié)構(gòu)化網(wǎng)格劃分。太少的網(wǎng)格計(jì)算時(shí)間雖然短,但氣液界面不夠清晰,影響模擬結(jié)果,太多的網(wǎng)格必然耗費(fèi)更多時(shí)間計(jì)算。在完成網(wǎng)格獨(dú)立性檢驗(yàn)后,確定節(jié)點(diǎn)數(shù)為32 844,網(wǎng)格數(shù)為31 280。

        2.2 控制方程

        低溫脈動(dòng)熱管內(nèi)部的流動(dòng)為氣液兩相流動(dòng)。應(yīng)用VOF 方法來(lái)跟蹤脈動(dòng)熱管每個(gè)單元中基于氣相和液相體積分?jǐn)?shù)αv和αl的表面界面,其中下標(biāo)v 和l 分別代表氣相和液相。在控制體積中

        在這個(gè)模型中,氣相被認(rèn)為是可壓縮的理想氣體,液相不可壓縮,給出了兩相體積分?jǐn)?shù)的連續(xù)性方程

        式中:ρv和ρl和分別為蒸氣和液體的密度,kg/m3;為速度矢量,m/s;Sl和Sv分別為液相和氣相的源項(xiàng),kg/(m3·s);和分別為從蒸氣到液體的質(zhì)量轉(zhuǎn)移以及從液體到蒸氣的質(zhì)量轉(zhuǎn)移,kg/(m3·s)。

        其中蒸發(fā)冷凝過(guò)程可以用Lee 模型表示。如果Tl>Tsat,發(fā)生蒸發(fā)過(guò)程

        如果Tv

        式中:Tl和Tv分別為液體和蒸氣的溫度,K;Tsat為飽和溫度,K;rlv和rvl為傳質(zhì)時(shí)間的松弛因子,這些系數(shù)影響數(shù)值收斂的精度和界面溫度。根據(jù)Lee 等人[7]的研究結(jié)果,通常設(shè)置為默認(rèn)值0.1,使界面溫度在數(shù)值上保持接近工質(zhì)的飽和溫度。

        動(dòng)量方程在式(6)中求解,它依賴(lài)于所有相的體積分?jǐn)?shù)。

        式中:ρ為控制體積內(nèi)的平均密度,kg/m3;為速度矢量,m/s;p為壓力,Pa;μ為黏度,Pa·s;為重力加速度,m/s2;為由彎曲界面引起的拉普拉斯方程計(jì)算的體積力,N/m3。

        表面張力是由于流體中分子間的內(nèi)聚力而產(chǎn)生的,它產(chǎn)生的表面力在小通道中占主導(dǎo)地位。對(duì)于表面張力的建模,采用ANSYS Fluent 中的連續(xù)表面力(Continuum Surface Force,CSF)模型[8]。在VOF 模型中,將表面張力作為動(dòng)量方程的源項(xiàng)加入到動(dòng)量方程中,如式(7)所示:

        本模型的能量方程為:

        式中:E為熱力學(xué)能,J/kg;λeff為有效導(dǎo)熱系數(shù),W/(m·K);Sh為相變引起的能量源項(xiàng),W/m3,它是由傳質(zhì)速率乘以潛熱Sh=-hLHlv=hLHvl得到的。

        2.3 數(shù)值模擬方法

        本研究應(yīng)用ANSYS Fluent 對(duì)建立的模型進(jìn)行求解。氦氣的物性隨溫度的變化而變化,為了使模擬更接近真實(shí)情況,采用美國(guó)國(guó)家標(biāo)準(zhǔn)技術(shù)研究所(NIST)開(kāi)發(fā)的Refprop 軟件獲得不同溫度下氦氣的物性并擬合成多項(xiàng)式的形式,其中密度、定壓比熱、熱導(dǎo)率、運(yùn)動(dòng)黏度系數(shù)的多項(xiàng)式擬合公式分別如式(9)、(10)、(11)、(12)所示。

        液氦被視為不可壓縮流體,因此其物性隨溫度變化可以忽略不計(jì),選用4.215 K 溫度下的物性參數(shù)進(jìn)行計(jì)算。

        兩相流模型VOF 方程采用顯式形式,無(wú)量綱庫(kù)朗數(shù)設(shè)置為0.25,選用隱式體積力方程。為使模型更快收斂,飽和溫度設(shè)置為4.215 K。其它關(guān)于模型的求解設(shè)置列于表1 中。

        表1 二維液氦脈動(dòng)熱管數(shù)值模擬Fluent 求解設(shè)置Table 1 Solution setting of two-dimensional helium pulsating heat pipe by Fluent numerical simulation

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

        3.1 液氦脈動(dòng)熱管初始狀態(tài)

        為了模擬液氦脈動(dòng)熱管的初始狀態(tài),將初始狀態(tài)每個(gè)網(wǎng)格的溫度設(shè)置為4.215 K,充液率為50%。邊界條件設(shè)置為恒溫4.215 K。采用ANSYS Fluent 求解計(jì)算,計(jì)算穩(wěn)定后得到脈動(dòng)熱管的初始?xì)庖悍植紶顟B(tài),如圖2 所示。圖2a 展示了二維單回路液氦脈動(dòng)熱管初始狀態(tài)的氣液分布云圖。可以看出在表面張力和重力的共同作用下,管內(nèi)工質(zhì)形成氣塞和液塞交替分布的靜止?fàn)顟B(tài),氣塞和液塞的大小和長(zhǎng)度是隨機(jī)的,且氣液分界面明顯。初始狀態(tài)下液氦脈動(dòng)熱管內(nèi)的溫度均保持在4.215 K 左右,蒸發(fā)段和冷凝段沒(méi)有溫差,如圖2b 所示。同時(shí),如圖2c 所示,初始狀態(tài)下液氦脈動(dòng)熱管處于飽和狀態(tài),不存在質(zhì)量傳遞。當(dāng)工質(zhì)充注到液氦脈動(dòng)熱管后,由于毛細(xì)作用,管內(nèi)的工質(zhì)會(huì)形成隨機(jī)分布的氣塞與液塞,初始狀態(tài)的氣液分布為液氦脈動(dòng)熱管的啟動(dòng)提供了基礎(chǔ)。

        圖2 二維單回路液氦脈動(dòng)熱管初始狀態(tài)云圖Fig 2 Initial state contour diagram of two-dimensional single loop nitrogen pulsating heat pipe

        3.2 液氦脈動(dòng)熱管運(yùn)行狀態(tài)的流動(dòng)分析

        改變邊界條件為蒸發(fā)段和冷凝段施加恒熱通量邊界,加熱和冷卻功率均為0.4 W,絕熱段設(shè)定為絕熱邊界條件,二維單回路液氦脈動(dòng)熱管進(jìn)入運(yùn)行狀態(tài)。圖3 展示了二維單回路液氦脈動(dòng)熱管運(yùn)行狀態(tài)的氣液分布云圖。從圖中可以看出,當(dāng)蒸發(fā)段和冷凝段的邊界改變后,液氦脈動(dòng)熱管管內(nèi)工質(zhì)不再保持靜止?fàn)顟B(tài)。在加入熱負(fù)荷的前0.4 s 內(nèi),由于蒸發(fā)段壁面被加熱,液膜吸收熱量并蒸發(fā),氣塞不斷產(chǎn)生并膨脹,同時(shí)冷凝段內(nèi)氣體溫度低于飽和溫度被冷凝為液體,氣塞收縮并破碎,因此蒸發(fā)段和冷凝段之間形成壓差,壓差驅(qū)動(dòng)力推動(dòng)兩根通道內(nèi)的氣塞會(huì)由于膨脹作用同時(shí)向冷凝段流動(dòng)。此外,在50%充液率下,管內(nèi)工質(zhì)的流型為塞狀流。在0.4 s 至1.4 s 時(shí)間段內(nèi),隨著蒸發(fā)段液膜的不斷蒸發(fā),氣塞逐漸變長(zhǎng),脈動(dòng)熱管內(nèi)的壓力差也不斷增加。此時(shí)左側(cè)通道內(nèi)的壓力低于右側(cè)通道內(nèi)的壓力,因此左側(cè)管道內(nèi)的工質(zhì)改變流動(dòng)方向向下流動(dòng),右側(cè)工質(zhì)繼續(xù)向上流動(dòng)。這是由于蒸發(fā)段氣化核心產(chǎn)生的位置具有一定的隨機(jī)性,蒸發(fā)段內(nèi)氣塞的位置和大小因此會(huì)產(chǎn)生差異。當(dāng)壓差足夠克服重力以及毛細(xì)力時(shí),管內(nèi)工質(zhì)向同一方向運(yùn)動(dòng)。此階段流動(dòng)方式為單方向循環(huán)流動(dòng)。在1.4 s至1.6 s 內(nèi),由于蒸發(fā)段的持續(xù)加熱,蒸發(fā)段內(nèi)溫度持續(xù)升高,核態(tài)沸騰發(fā)生,蒸發(fā)段液膜吸熱蒸發(fā)更加劇烈,這時(shí)左側(cè)通道由于強(qiáng)烈的氣液相變,產(chǎn)生了更大的推動(dòng)力,工質(zhì)改變流動(dòng)方向,左側(cè)通道為上升管,右側(cè)通道為下降管,工質(zhì)流動(dòng)方式為逆向流動(dòng)。因此,在運(yùn)行階段,該單回路液氦脈動(dòng)熱管管內(nèi)工質(zhì)經(jīng)歷了以下3 種流動(dòng)狀態(tài):兩側(cè)通道內(nèi)同時(shí)向冷凝段運(yùn)動(dòng)、單方向循環(huán)流動(dòng)、逆向循環(huán)流動(dòng)。

        圖3 二維單回路液氦脈動(dòng)熱管運(yùn)行狀態(tài)氣液分布云圖Fig.3 Gas-liquid distribution contour diagram of working state of two-dimensional single loop helium pulsating heat pipe

        3.3 液氦脈動(dòng)熱管運(yùn)行狀態(tài)的傳熱特性

        圖4 展示了二維單回路液氦脈動(dòng)熱管的溫度分布云圖。該模擬工況下,隨著蒸發(fā)段和冷凝段壁面以恒定功率傳熱,脈動(dòng)熱管內(nèi)蒸發(fā)段的溫度逐漸升高,反之,冷凝段的溫度逐漸降低,蒸發(fā)段的溫度在1.2 s時(shí)升至最高5.6 K,冷凝段的溫度則降至3.4 K。在0.6 s 至1.2 s 期間,由于工質(zhì)攜帶熱負(fù)荷往逆時(shí)針?lè)较蜻\(yùn)動(dòng),右側(cè)通道為上升管,左側(cè)通道為下降管,因此,右側(cè)通道內(nèi)蒸發(fā)段的熱流也隨著工質(zhì)流動(dòng)方向朝冷凝段運(yùn)動(dòng)。同時(shí),左側(cè)通道內(nèi)的熱流朝蒸發(fā)段運(yùn)動(dòng)。在1.2 s 以后,由于脈動(dòng)熱管管內(nèi)工質(zhì)的運(yùn)動(dòng)方向改變?yōu)轫槙r(shí)針循環(huán)流動(dòng),右側(cè)通道內(nèi)熱流轉(zhuǎn)變?yōu)橄蛘舭l(fā)段運(yùn)動(dòng),左側(cè)通道內(nèi)的熱流轉(zhuǎn)變?yōu)橄蚶淠芜\(yùn)動(dòng)。由此可見(jiàn),單回路液氦脈動(dòng)熱管運(yùn)行過(guò)程中蒸發(fā)段溫度升高,冷凝段溫度降低,其傳熱方向與工質(zhì)流動(dòng)方向相同。

        圖4 二維單回路液氦脈動(dòng)熱管運(yùn)行狀態(tài)溫度分布云圖Fig.4 Temperature contour diagram of working state of two-dimensional single loop helium pulsating heat pipe

        圖5 展示了二維單回路液氦脈動(dòng)熱管運(yùn)行狀態(tài)的質(zhì)量傳遞云圖。從圖中可以看到,單回路液氦脈動(dòng)熱管運(yùn)行時(shí)相變和質(zhì)量傳遞主要在蒸發(fā)段和冷凝段發(fā)生。由于蒸發(fā)段有熱量輸入,通道內(nèi)的液氦溫度超過(guò)飽和溫度,產(chǎn)生核態(tài)沸騰,液膜蒸發(fā),產(chǎn)生氣泡,且處于蒸發(fā)段的氣塞也會(huì)隨著液膜蒸發(fā)而膨脹。反之,冷凝段內(nèi)有熱量輸出,工質(zhì)溫度低于飽和溫度,處于過(guò)冷狀態(tài),氦氣冷凝為液相,氣塞逐漸收縮。蒸發(fā)段和冷凝段由于相變傳質(zhì)而具有壓力梯度,驅(qū)動(dòng)管內(nèi)的工質(zhì)在蒸發(fā)段和冷凝段之間運(yùn)動(dòng),從而實(shí)現(xiàn)熱量的傳遞。由此可知,液氦脈動(dòng)熱管運(yùn)行過(guò)程中由于冷熱端的溫度差產(chǎn)生的相變是造成壓力差和工質(zhì)流動(dòng)的原因,而工質(zhì)流動(dòng)實(shí)現(xiàn)了熱量的傳遞。由模擬結(jié)果可知,液氦脈動(dòng)熱管自身具有自激振蕩的特性,因此不需要增加額外的驅(qū)動(dòng)功率就能傳遞熱負(fù)荷。

        圖5 二維單回路液氦脈動(dòng)熱管運(yùn)行狀態(tài)質(zhì)量傳遞云圖Fig.5 Mass transfer contour diagram of working state of two-dimensional single loop helium pulsating heat pipe

        4 結(jié)論

        為了探究液氦脈動(dòng)熱管的運(yùn)行機(jī)理,使用ANSYS Fluent 中的多相流VOF 模型、表面張力模型(CSF 模型)以及傳熱傳質(zhì)模型(Lee 模型),模擬了二維單回路液氦脈動(dòng)熱管的初始狀態(tài)氣液分布以及運(yùn)行狀態(tài)的流動(dòng)與傳熱。其中初始狀態(tài)的氣液交替分布為液氦脈動(dòng)熱管的啟動(dòng)提供了基礎(chǔ)。該單回路液氦脈動(dòng)熱管運(yùn)行時(shí)管內(nèi)工質(zhì)經(jīng)歷了3 種流動(dòng)狀態(tài):兩側(cè)通道內(nèi)同時(shí)向冷凝段運(yùn)動(dòng)、單方向循環(huán)流動(dòng)、逆向循環(huán)流動(dòng)。在50% 充液率下,該脈動(dòng)熱管的流型為塞狀流。通過(guò)數(shù)值模擬發(fā)現(xiàn),在薄膜蒸發(fā)/冷凝和氣泡壓縮/膨脹的共同作用下,蒸發(fā)段內(nèi)的氣塞溫度大于相同壓力下的飽和溫度,處于過(guò)熱狀態(tài);冷凝段內(nèi)的液塞溫度小于相同壓力下的飽和溫度,處于過(guò)冷狀態(tài)。由模擬結(jié)果可知,液氦脈動(dòng)熱管依靠?jī)?nèi)部的工質(zhì)流動(dòng)及蒸發(fā)冷凝來(lái)實(shí)現(xiàn)熱量的傳遞,運(yùn)行的動(dòng)力來(lái)源于管內(nèi)的壓力差。模擬結(jié)果對(duì)理解液氦脈動(dòng)熱管的運(yùn)行和傳熱機(jī)制有很大幫助。

        猜你喜歡
        液氦工質(zhì)氣液
        海洋溫差能發(fā)電熱力循環(huán)系統(tǒng)的工質(zhì)優(yōu)選
        我國(guó)自主研發(fā)液氦罐箱首次進(jìn)入批量生產(chǎn)
        液氦儲(chǔ)罐發(fā)展現(xiàn)狀及關(guān)鍵技術(shù)
        低溫與特氣(2021年5期)2021-11-12 07:31:00
        微重力下兩相控溫型儲(chǔ)液器內(nèi)氣液界面仿真分析
        采用R1234ze(E)/R245fa的非共沸混合工質(zhì)有機(jī)朗肯循環(huán)系統(tǒng)實(shí)驗(yàn)研究
        氣液分離罐液位計(jì)接管泄漏分析
        歷史的變革!MRI即將進(jìn)入“100%無(wú)液氦時(shí)代”
        采用二元非共沸工質(zhì)的有機(jī)朗肯循環(huán)熱力學(xué)分析
        CO2 驅(qū)低液量高氣液比井下氣錨模擬與優(yōu)化
        若干低GWP 純工質(zhì)在空調(diào)系統(tǒng)上的應(yīng)用分析
        久久久国产精品五月天伊人| 国产成人无码免费看片软件| 国产精品第一国产精品| 国产av人人夜夜澡人人爽麻豆| 国产真实乱人偷精品人妻| 婷婷一区二区三区在线| 最新国产精品精品视频| 亚洲精品色播一区二区| 日本在线观看一二三区| 粗大的内捧猛烈进出少妇| 粗大猛烈进出白浆视频| 亚洲AV无码一区二区三区日日强| 亚洲丁香五月激情综合| 亚洲av自偷自拍亚洲一区| 国产精品国产三级国产av品爱| 无码国产精品久久一区免费| 双腿张开被9个男人调教| 欧美极品少妇性运交| 亚洲AV小说在线观看| 最新亚洲视频一区二区| 亚洲youwu永久无码精品| 极品少妇小泬50pthepon| 人妻av中文字幕无码专区| 午夜视频网址| 亚洲码专区亚洲码专区| 日韩不卡的av二三四区| 国产又大又黑又粗免费视频| 日韩一卡2卡3卡4卡新区亚洲| 国产精品无码Av在线播放小说| 国产经典免费视频在线观看| 蜜桃视频永久免费在线观看| 久久一二区女厕偷拍图| 日韩亚洲欧美中文在线| 大地资源中文第三页| 人妻少妇精品系列一区二区| 日韩女优视频网站一区二区三区 | 97午夜理论片影院在线播放| 久久婷婷香蕉热狠狠综合| 国产精品久久久久久久久免费观看 | 区二区欧美性插b在线视频网站| 亚洲色欲色欲www成人网|