?

有序渦旋對三角槽道脈動流強化傳熱的影響

2016-10-14 00:03黃其王勛廷楊志超鐘英杰
化工學報 2016年9期
關鍵詞:溫度梯度壓力梯度脈動

黃其,王勛廷,楊志超,鐘英杰

?

有序渦旋對三角槽道脈動流強化傳熱的影響

黃其,王勛廷,楊志超,鐘英杰

(浙江工業大學能源與動力工程研究所脈動技術工程研究中心,浙江 杭州 310014)

以水為工質對三角槽道內單相液體充分發展層流脈動傳熱特性進行了研究。應用粒子圖像測速技術(PIV)測得流場內渦的變化規律,從“渦及渦運動”的角度揭示了“有序”的渦生長及遷移過程對脈動流強化傳熱的影響。此外,應用“場協同”理論,通過數值模擬深入分析了流場特性與傳熱之間的關系。研究發現,“有序”的渦生長及遷移過程,破壞了流體邊界層,促進了近壁區熱流場與速度場的協同,同時,強化了三角槽道內流體與主流區流體的摻混,熱量輸運能力提升;存在最佳的Strouhal數(),使得渦旋既能充分發展又能在較短時間脫落進入主流,實現最大效率的壁面換熱;有序渦旋對速度場、溫度梯度場以及壓力梯度場三者協同性的改善是換熱性能提升的關鍵。

傳熱;脈動流;層流;渦;協同性;實驗測量;數值分析

引 言

熱量的傳遞過程廣泛存在于動力、核能、制冷、化工、石油、航空等行業中。換熱器在上述各行業中不僅是保證工程設備正常運轉不可或缺的部件,而且在金屬消耗、動力消耗和投資方面占有重要份額。例如在制冷機中,換熱器的動力消耗占總能耗的20%~30%[1];化工企業里,換熱器的能耗占到企業總能耗的40%以上。因此,高效節能的傳熱傳質設備的研發與應用,一直是流程工業所需解決的關鍵技術之一。脈動流傳熱技術在新一代的換熱技術領域表現出了巨大的潛力,脈動流傳熱的機理研究對節能高效新型換熱器的開發以及我國的節能戰略意義重大。

對于脈動流強化傳熱的機理,目前尚無統一的結論,學者根據自己研究的工況提出了很多脈動流傳熱機理假設,但并不全面。如“減薄邊界層”無法解釋脈動流對圓管層流流動不能強化換熱的現象;“回流”機理很難解釋頻率與傳熱間存在最佳值的現象;“流體共振”分析了頻率對于脈動傳熱的影響,卻無法應用于大幅度脈動換熱問題[2-4]。

值得注意的是,在對不同脈動流強化傳熱機理的詮釋過程中,均提及“流體的摻混現象”。近10年來,通過對不同周期性幾何結構流道中脈動流強化傳熱的模擬研究,諸多學者相繼提出:流體的摻混是脈動流強化傳熱的主要原因之一[5-11]。不過,對于摻混現象的研究依然存在如下問題:第一,通過可視化實驗,以“渦及渦運動”為視角[12-13],研究具有特殊結構流道的液態脈動流強化傳熱過程的報道相對較少;第二,有關“摻混現象對于強化傳熱的貢獻”的解釋還有待進一步深入。

由此,本文借助脈動流傳熱實驗臺和PIV測試系統,將周期性三角槽道作為研究對象,從實驗角度觀察周期性幾何結構中脈動流傳熱的近壁流層熱流場與速度場的協同關系以及與主流區的摻混過程,解析流動與換熱之間的關系。同時,基于場協同理論[14-15],通過數值模擬進一步闡述“渦及渦運動”對脈動流場不同區域的速度場、溫度梯度場、壓力梯度場所造成的影響。

1 實驗部分

1.1 脈動傳熱實驗臺

脈動傳熱實驗臺可分為:實驗段、電加熱系統、水循環系統、數據采集系統、可視化實驗系統。其中,實驗段選用三角槽道,尺寸的確定以有利于細致觀察渦的生長、遷移過程為原則而定。

圖1為自行搭建的脈動傳熱實驗臺。增壓水泵作為動力裝置,為系統提供水循環動力;隔膜計量泵輸出脈動壓力,在混合室中將恒定流轉變為有平均流速的脈動流;變壓器用以調節粘貼在槽道金屬壁面上的加熱片供熱功率,實現等熱流換熱的實驗條件;NI-DAQ多功能數據采集儀用于實時獲取溫度、壓力等實驗數據,并與計算機進行連接;外接PIV測試系統,對實驗流體進行圖像分析。

1.2 流場可視化(PIV)實驗系統

PIV系統主要由高速CCD相機、雙腔激光器及同步器、后處理計算機三大部分組成,旨在獲取流道中的流場結構。

實驗中,示蹤粒子選用平均直徑為25 μm的聚丙乙烯。利用高速CCD相機,拍攝相應流場的示蹤粒子圖片,通過自適應互相關算法處理得到流場結構圖。所有的流場結構數據均取自流動狀態已經充分發展的三角槽道。

速度矢量由式(1)計算得到

式中,表示迭代區域內粒子的平均位移;為粒子的瞬時速度。

圖像后處理中,采用中心差分格式將式(1)轉變為代數方程再進行求解[16]。

1.3 實驗參數設定

脈動流傳熱實驗工況為:300≤≤900,脈動頻率范圍為0~50 Hz,脈動振幅為平均流速的20%~100%,并保證實驗段壁溫與換熱工質溫差在5~10℃。

2 模擬部分

2.1 物理模型

如圖2所示,三角槽深為9 mm,流道高度為15 mm,三角槽長度為18 mm,流道寬度為187 mm。換熱段總長為590 mm,遠大于寬度和高度。因此,簡化為二維模型進行數值計算。選用水作為換熱工質。

2.2 計算方法及邊界條件

以文獻[17-18]對于周期性脈動流動的假設為參考,本文選擇非穩態層流模型,主體部分四邊形網格大小為0.4 mm×0.4 mm,三角槽和主流區上部區域進行網格加密,網格總數為135533個。使用SIMPLE算法構造壓力修正方程,設置連續性方程和動量方程的松弛因子分別為0.3和0.7。將所有方程的殘差收斂標準均設為1×10?6,迭代時間步長設置為0.01 s。

邊界條件的設置如下。

(1)入口 穩態計算時設為恒定速度邊界條件,脈動流計算時用UDF自定義函數

=s+ssin() (2)

式中,s為平均流速;為速度脈動振幅。

(2)出口 設定為pressure-out。

(3)上壁面和下壁面均為wall邊界條件,上壁面設置恒定熱流。

模擬參數:平均流速為0.02 m·s-1(=300),脈動頻率為0.5 Hz,脈動振幅為平均流速的50%。

3 實驗結果與討論

3.1 傳熱性能實驗結果分析

定義為

式中,為流道截面的當量直徑;為流體的熱導率;為換熱量;Δ為壁面與流體之間的溫差。

表征量綱1化的脈動頻率,即

式中,為脈動頻率;為平均流速。

用強化傳熱因子反映相同槽道條件下,脈動流與穩定流傳熱的差異,計算式為

式中,p為脈動流在一個脈動周期內的平均值;s為恒定入口流速的流體在相同時間段的平均值。

圖3為不同條件下,強化傳熱因子隨的變化。結果表明,在不同條件下,隨的變化趨勢一致,均為先增后減,存在最佳使達到極值點;在和一定的條件下,越小,越大。

圖4為不同下,強化傳熱因子隨的變化。結果表明,在不同條件下,隨的變化趨勢一致,均為先增后減,且達到最佳傳熱效率時的相近;在和一定的條件下,的大小決定了脈動流傳熱效果的優劣。

綜合考量,在實驗工況內,均大于1,相較于穩定流,脈動流具有更好的換熱效果;越小,脈動流傳熱強化效果越明顯;存在最佳脈動頻率使得脈動流傳熱效果最優;高脈動振幅可以促進脈動流傳熱強化。因此,在層流脈動流換熱情況下,三角槽流道中采用低、高、合適,可以大幅提升換熱性能。

3.2 可視化實驗結果分析

(1)渦的形成狀態分析

圖5為一個脈動周期內,三角槽道中渦變化的一系列“有序”過程。從=0到=/6為渦的發展階段,渦形成后,渦中心逐漸由上壁面向槽道中心移動,并且渦徑不斷增大,表現出近壁區良好的熱流-速度協同現象;從=2/6到=3/6為充分發展階段,渦在三角槽中心處不斷擴張,從主流區不斷吸入流體,表現出劇烈的槽道流體區與主流區的流體摻混現象;從=4/6到=為渦消逝階段,渦中心開始向主流區域移動,并且渦流逐漸融入主流,渦尺寸逐漸減少直至消失。

從渦的變化可知,三角槽道內層流脈動流可以強化傳熱的原因與渦的形成和消散具有密切的聯系。渦的出現,有效破壞了流體邊界層,促進了近壁區熱流場與速度場的協同,同時,不僅強化了槽道內流體的摻混,而且使得主流區與三角槽道內流體的交換能力加強,主流區的流體不斷進入槽道,經過換熱的流體也被主流快速帶走,熱量輸運的總體能力提高。

(2)不同對渦的影響分析

圖6為在不同下,三角槽道中渦的形成狀況。當=0.15時,渦也經歷了形成、發展、消逝的過程,但是渦在槽道中心處并未充分發展擴張就趨向于消逝;當=0.6時,三角槽道內出現了上下兩個渦,未經充分發展,兩渦就迅速融合后消散;當=0.36時,相較于前兩種,渦的發展較充分,換熱效果明顯提高??梢?,對三角槽道中渦的形成發展產生重大影響。

圖7為在不同條件下,三角槽道中渦中心的運動軌跡。當=0.15時,渦中心在=3/6后長時間貼近上壁面;當=0.36時,渦中心較長時間處于三角槽道中央,具有充分發展的空間;當=0.6時,一號渦中心不進入三角槽道中心,與二號渦發生碰撞融合,此后迅速脫離三角槽道區域并進入主流??梢?,合適的能夠優化渦的運動軌跡。

圖8為不同條件下,渦徑隨著時間的變化。結果表明,相較于其他兩種工況,同時刻=0.36渦徑明顯較大,渦的發展程度較好,受到渦影響的換熱區域也最大。

綜上所述,從渦在不同工況下的變化規律可知,當過低時,渦具有充足的發展時間,但受限于流速較低的影響,使渦開始擴張時過于靠近壁面而無法充分發展。此外,由于周期時間較長,造成一定時間內進入主流的渦的數量相對較少;當過高時,渦在三角槽道內未充分發展就迅速脫落進入主流,未實現流體充分摻混。因此,存在合適的,使得渦既能充分發展又能在較短時間脫落進入主流,以最大效率實現壁面換熱。

4 模擬結果與討論

4.1 模擬可靠性驗證

從圖9可以發現,模擬與實驗的結果變化趨勢相似,最大偏差控制在10%以內,平均偏差為5%左右。出現誤差的原因存在兩方面:第一,實驗本身存在的測量誤差;第二,由于實驗系統的加熱片直接貼于換熱面,可能出現貼合不緊密以及絕熱保溫不充分的現象,但模擬計算是按絕熱面處理的,這也會產生一定的誤差。因此,綜合考慮后,本模擬結果相對準確,并具有一定的數值實驗價值。

4.2 場協同分析

為了表征速度場和溫度梯度場的協同性,定義場協同數為[19]

式中,為量綱1化的速度矢量;?為量綱1化的溫度梯度矢量。

全場的速度和溫度梯度的平均協同數為

式中,dS為第個微元面。

由圖6(b)和圖10對比分析可知,=0時,渦在三角槽道上游形成,主流區上游與初始渦交匯處出現一條主高帶,但槽道中流體與主流交匯界面處依然存在較大范圍的低區;=/6時,渦開始變大,并向三角槽道中心運動,渦與主流交匯界面增大,主高帶也隨之移動和增大,低區減小,槽道平均變大;=2/6至=3/6時,渦充分發展,并在三角槽中心處不斷擴張,槽道內流體與主流區進行了充分的摻混,交匯界面處的主高帶影響區域擴大,前兩個時刻產生的低區域幾乎消失,有效遏制了主低帶的發展,槽道達到最大值;=4/6時,渦的中心越過三角槽道中心線,渦開始變小,交匯界面縮小并逐漸向下游移動,主高帶也隨之向下游移動,主低帶在槽道上游開始擴張并逐漸開始影響整個槽道。由于依然受到主高帶的較大影響,槽道平均雖有所減少但依然處于較高值;=5/6時,渦逐漸消散,主高帶消失,三角槽道內出現較大部分的低區域。

通過上文對云圖的分析可以發現,槽道內脈動流可以起到強化傳熱作用的關鍵在于“渦及渦運動”。渦旋的出現,使得槽道內渦旋與主流區交匯處出現高帶,此處將槽道流體與主流進行了充分的摻混和能量交換,主流輸運能力得到大幅提升,并改善了整個流道的速度場與溫度梯度場的協同性;槽道與渦徑變化規律一致,渦的尺寸對流場熱交換能力呈現正相關的關系。

速度場與溫度梯度場的協同性能改善是提高換熱的有效手段;同時,流動阻力也會影響能耗的產生量。因此,對層流脈動流流動阻力做進一步分析。

由動能方程得到壓力梯度的做功功率為[20]

式中,為速度矢量;?為壓力梯度;為速度矢量與壓力梯度之間的夾角。夾角越大,壓力做功能力越強,產生的壓降越小。

由式(7)推導可得,速度矢量與壓力梯度的協同角為

全場的周期平均協同角為

由上文表述可知,速度場與溫度梯度場的協同性較優時,壓力梯度場與速度場的協同程度也較好。從三角槽道流態的角度做進一步分析,和一定的條件下,處于高的脈動流態,渦旋表現出良好的形成狀態,對邊界層的破壞相對較大,能夠遏制一定范圍內的黏滯力增長,壁面剪切力作用有所下降,故而呈現壓損下降的現象。因此,在工況一定的條件下,選擇合適的,使得速度場、溫度梯度場以及壓力梯度場的三場協同性較優,不但可以促進強化傳熱,而且也可以實現壓損增量的減小。

綜上所述,由脈動流與三角槽道結構共同作用而引起的渦的“有序”生長及遷移過程,對流場的換熱性能及流動阻力產生重大影響。當處于合適范圍,有序渦的發展程度較好,速度場、溫度梯度場以及壓力梯度場均可達到較優的協同程度,使得流體摻混充分、熱量輸運迅速、壓損增量減小,從而實現高效低阻的換熱目標。

5 結 論

(1)在300<<900實驗工況內,當脈動流處于低、高、合適時,可以大幅提升換熱性能,比相同結構槽道的穩定流換熱效率提高近2倍。

(2)由層流脈動流與三角槽道結構共同作用而引起的渦的“有序”生長及遷移,有效破壞了流體邊界層,強化了槽道中流體的摻混,經過換熱的流體被主流快速帶走,熱量輸運能力提高。

(3)對槽道中渦的形成狀況產生重大影響,直接影響近壁區流體與主流區流體的熱量輸運。當過低或者過高時,渦均無法進行較好的發展,無法實現流體的充分摻混。存在一個合適范圍的使得渦能夠良好發展而提高換熱效率。

(4)從脈動流傳熱三場協同角度分析,三角槽道內層流脈動流能夠提升換熱性能的原因為:“有序”的渦生長及遷移過程,改善了速度場、溫度梯度場以及壓力梯度場三者的協同性,實現了提高換熱性能并盡量減少流動阻力的目的。

符 號 說 明

A——脈動振幅 d——直徑,mm D——平均位移 E——強化傳熱因子 f——脈動頻率,s?1 Fc——場協同數 k——熱導率,W·m?1·K?1 Np——做功功率 Nu——Nusselt數 Pr——Prandtl數 ?p——壓差,Pa q——熱通量,W·m?2 Re——Reynolds數 Si——微元面 ΔT——溫差,K U——速度,m·s?1 θ——速度矢量與壓力梯度的夾角

References

[1] 林宗虎, 汪軍, 李瑞陽, 等. 強化傳熱技術[M]. 北京: 化學工業出版社, 2007.
LIN Z H, WANG J, LI R Y,. The Technology of Heat Transfer Augmentation [M]. Beijing: Chemical Industry Press, 2007.

[2] TAO H J, SEO Y K, JAE M H. Experiments on heat transfer enhancement from a heated square cylinder in a pulsating channel flow [J].International Journal of Heat and Mass Transfer, 2008, 51(5/6): 1130-1138.

[3] MOHAMMAD J, MOUSA F, KUROSH S. Pulsating flow effects on convection heat transfer in a corrugated channel: a LBM approach [J]. International Communications in Heat and Mass Transfer, 2013, 45(7): 146-154.

[4] HIMADRI C, FRANZ D, SUBHASHIS R. Analysis of heat transfer in simultaneously developing pulsating laminar flow in a pipe with constant wall temperature [J]. International Communications in Heat and Mass Transfer, 2006, 33(4): 475-481.

[5] JIN D X, LEE Y P, LEE D Y. Effects of the pulsating flow agitation on the heat transfer in a triangular grooved channel [J]. Heat Mass and Transfer, 2007, 50: 3062-3071.

[6] 何雅玲, 楊衛衛, 趙春鳳, 等. 脈動流動強化換熱的數值研究[J]. 工程熱物理學報, 2005, 26(3): 495-497.
HE Y L, YANG W W, ZHAO C F,. Numerical study on enhancing heat transfer by pulsating flow [J]. Journal of Engineering Thermophysics, 2005, 26(3): 495-497.

[7] 路慧霞, 馬曉建, 趙凌. 脈動流動強化傳熱的研究進展[J]. 節能技術, 2008, 26(148): 168-172.
LU H X, MA X J, ZHAO L. Research advancements of heat transfer enhancement by pulsating flow [J]. Energy Conservation Technology, 2008, 26(148): 168-172.

[8] 賈寶菊, 孫發明, 卞永寧, 等. 波壁管內的脈動流動及傳質強化的數值模擬[J]. 化工學報, 2009, 60(1): 6-12.
JIA B J, SUN F M, BIAN Y N,. Numerical simulation of pulsatile flow and mass transfer enhancement in a wavy-walled tube [J]. CIESC Journal, 2009, 60(1): 6-12.

[9] 陰繼翔, 楊剛, 郭瑞, 等. 非對稱波紋通道內流動與傳熱的數值計算[J]. 太原理工大學學報, 2011, 42(1): 20-24.
YIN J X, YANG G, GUO R,. Numerical calculation of fluid flow and heat transfer in asymmetrical wavy channels [J]. Journal of Taiyuan University of Technology, 2011, 42(1): 20-24.

[10] LI Y, JIN D X, JING Y Q,. An experiment investigation of heat transfer enhancement by pulsating laminar flow in rectangular grooved channels [J]. Advanced Materials Research, 2013, 732(733): 74-77.

[11] YANG B C, JIN D X. An experimental investigation of heat transfer enhancement by pulsating laminar flow in a triangular grooved channel [J]. Advanced Materials Research, 2012, 516(517): 249-252.

[12] 謝公南, 王秋旺, 曾敏, 等. 漸擴漸縮波紋通道內脈動流的傳熱強化[J]. 高?;瘜W工程學報, 2006, 20(1): 31-35.
XIE G N, WANG Q W, ZENG M,. Heat transfer enhancement inside converging-diverging wavy channel by pulsating flow [J]. Journal of Chemical Engineering of Chinese Universities, 2006, 20(1): 31-35.

[13] ESAM M A, RAED I B. The role of periodic vortex shedding in heat transfer enhancement for transient pulsatile flow inside wavy channels [J]. Proceedings of World Academy of Science EngineeringTechnology, 2008, (22): 446-452.

[14] 過增元. 對流換熱的物理機制及其控制[J]. 科學通報, 2001, 45(19): 2118-2122.
GUO Z Y. The convective heat transfer of physical mechanism [J]. Science in China, 2001, 45(19): 2118-2122.

[15] 過增元, 魏澍, 程新廣, 等. 換熱器強化的場協同原則[J]. 科學通報, 2003, 48(22): 2324-2327.
GUO Z Y, WEI S, CHENG X G,. The field synergy principle in heat exchanger [J]. Science in China, 2003, 48(22): 2324-2327.

[16] 張健華, 李思文, 沈忠良, 等. PIV技術在脈動流傳熱實驗中的應用[J]. 浙江工業大學學報, 2014, 42(1): 109-113.
ZHANG J H, LI S W, SHENG Z L,. The application of PIV technology on the pulsating heat transfer experiment [J]. Journal of Zhejiang University of Technology, 2014, 42(1): 109-113.

[17] 陶文銓. 數值傳熱學[M]. 西安: 西安交通大學出版社, 2001.
TAO W Q. Numerical Heat Transfer [J]. Xi’an: Xi’an Jiaotong University Press, 2001.

[18] 王永慶, 董其伍, 劉敏珊. 混沌對流強化傳熱的場協同分析[J]. 鄭州大學學報, 2011, 32(3): 6-9.
WANG Y Q, DONG Q W, LIU M S. Analysis of field synergy on heat transfer enhancement in chaotic advection [J]. Journal of Zhengzhou University, 2011, 32(3): 6-9.

[19] 李志信, 過增元. 對流傳熱優化的場協同理論[M]. 北京: 科學出版社, 2010.
LI Z X, GUO Z Y. The Theory of Field Synergy on Convection Heat Transfer Optimization [M]. Beijing: Science Press, 2010.

[20] 何雅玲, 雷勇剛, 田麗亭, 等. 高效低阻強化換熱技術的三場協同性探討[J]. 工程熱物理學報, 2009, 30(11): 1904-1906.
HE Y L, LEI Y G, TIAN L T,. An analysis of three-field synergy on heat transfer augmentation with low penalty of pressure drop [J]. Journal of Engineering Thermophysics, 2009, 30(11): 1904-1906.

Influence of vortex on heat transfer enhancement in triangular grooved channel by pulsating flow

HUANG Qi, WANG Xunting, YANG Zhichao, ZHONG Yingjie

(Engineering Research Center of Pulse Technology, Institute of Energy and Power Engineering, Zhejiang University of Technology, Hangzhou 310014, Zhejiang, China)

In this paper, heat transfer enhancement in the triangular grooved channel by a laminar pulsating flow is studied. The influence of several main parameters on heat transfer enhancement is analyzed. The parameters are Reynolds number, Strouhal number and pulsation amplitude. The experimental results show that the enhancement of heat transfer rate increases with the Reynolds number and pulsation amplitude, and there exists an optimal Strouhal number for the greatest enhancement of heat transfer in the triangular grooved channel. To analyze the correlation between the pulsating flow behaviors and the heat transfer enhancement characteristics, the PIV investigation is performed. The PIV results show that the heat transfer enhancement results from the strong mixing caused by the repeating sequence of vortex generation, growth, expansion and ejection from the groove to the main stream by the pulsating flow. The repeating sequence of vortex variation changes the flow pattern, which leads to destroy the boundary layer and make the mixing speed faster in the different zones. What’s more, the numerical research has been conducted to investigate the synergy of the temperature, velocity and pressure fields on the laminar pulsating flow in a triangular grooved channel. The numerical results indicate that an increase of the intersection angle between velocity and pressure gradient improves the synergy between the velocity and pressure fields with an equal heat transfer enhancement, resulting in a reduction of penalty of pressure drop. Therefore, the improvement of three-field synergy is the basic mechanism for the heat transfer enhancement in the triangular grooved channel by a laminar pulsating flow.

heat transfer; pulsating flow; laminar flow; vortex; synergy; experimental measurement; numerical analysis

supported by the Funded Project of Science Technology Department in Zhejiang Province (2014C31034).

date: 2016-03-02.

Prof. ZHONG Yingjie, zhong_yingjie@zjut. edu.cn

TK 11

A

0438—1157(2016)09—3616—09

10.11949/j.issn.0438-1157.20160231

浙江省科技廳資助項目(2014C31034)。

2016-03-02收到初稿,2016-05-11收到修改稿。

聯系人:鐘英杰。第一作者:黃其(1988—),男,博士研究生。

猜你喜歡
溫度梯度壓力梯度脈動
壓力梯度對湍流邊界層壁面脈動壓力影響的數值模擬分析
RBI在超期服役脈動真空滅菌器定檢中的應用
致密-低滲透油藏兩相啟動壓力梯度變化規律
嚴寒地區混凝土箱梁實測溫度梯度分析
溫度梯度場對聲表面波器件影響研究
基于概率需求的高速鐵路無砟軌道板溫度荷載取值研究Ⅱ:溫度梯度作用
高速鐵路CRTSⅢ型板式無砟軌道溫度梯度試驗研究
疊加原理不能求解含啟動壓力梯度滲流方程
有限水域水中爆炸氣泡脈動的數值模擬
海上稠油砂巖油藏啟動壓力梯度測定方法及應用——以秦皇島32-6油田為例
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合