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

天線陣面沸騰換熱的數值分析及優化

2014-08-07 12:17:01殷翔錢吉裕孔祥舉吳華根曹鋒
西安交通大學學報 2014年11期
關鍵詞:研究

殷翔,錢吉裕,孔祥舉,吳華根,曹鋒

(1.西安交通大學能源與動力工程學院, 710049, 西安; 2.中國電子科技集團第14研究所, 210039, 南京)

天線陣面沸騰換熱的數值分析及優化

殷翔1,錢吉裕2,孔祥舉2,吳華根1,曹鋒1

(1.西安交通大學能源與動力工程學院, 710049, 西安; 2.中國電子科技集團第14研究所, 210039, 南京)

為研究雷達天線陣面的冷卻降溫效果,建立了兩相流蒸發冷卻模型,并著重分析了沸騰換熱在局部高熱流密度條件下的換熱能力,局部熱流高達800 kW/m2。采用VOF模型結合用戶自定義控制方程,數值計算三維流道內的沸騰傳熱現象,從氣液流動趨勢上尋找阻礙換熱的因素。結果表明:VOF模型能較好地用于分析氣液兩相流動中的層狀流、泡狀流以及氣液流動趨勢;通過減小彎道處過流面積提高局部流速的方法可以緩解氣相在彎道處的匯聚,消除了0.9 K的局部過熱;對于單側高熱流蒸發的數值分析,蒸發控制方程的調節系數在500左右為宜。溫度計算結果與實驗結果趨勢一致,總體偏差5K左右,數值分析可以用于研究冷板的沸騰換熱。

VOF模型;蒸發冷卻;沸騰換熱;數值計算

在電子設備領域,電流的發熱現象普遍存在且嚴重危害著設備的正常使用,因而對電子設備的冷卻降溫研究顯得非常重要。尤其是在具有較高熱流密度的電子設備上,冷卻不足往往導致設備過熱損壞,而常規熱容冷卻方式受冷卻極限和運行結構的制約,很難在這一領域滿足實際應用[1]。中國科學院電工研究所將蒸發潛熱冷卻方式應用到大型電機的定子冷卻領域,但依舊以實驗為主要研究方式[2]。

在低熱流密度沸騰換熱方面,實驗和數值模擬的研究都有所發展[3-4],但數值計算傳質傳熱的精度總體偏低。姜彭采用CFX相變蒸發模型研究水霧冷卻的換熱性能,數值計算與實驗的溫差在10 K左右[5]。在高熱流領域,數值研究基本限于微細通道內以及微觀結構表面的氣泡動力行為方面。顓瑞等對高熱流微細通道作了數值研究,觀察了氣泡的產生及湮滅過程,但沒有相關實際應用的數值研究[6-7]。國外許多學者采用VOF模型數值研究微觀表面的沸騰現象,以及微細通道內的氣泡生長[8-9]。總體而言,數值研究蒸發冷卻技術在實際三維換熱器中的應用相對較少,而在高熱流領域,傳熱傳質涉及的問題更為復雜,傳熱性能以及流動性能較一般流動差異較大。常規流道較微細通道的數值研究條件更為苛刻,而局部具有高熱流密度的沸騰傳熱特性與一般流動的差異性就更大了。

本文考慮雷達天線的運動特性以及局部具有較高熱流密度的特征,著重采用數值方法研究冷板換熱器的蒸發冷卻技術,主要涉及管道內沸騰過程的流動以及傳熱傳質問題。通過對矩形蛇形管流道內的沸騰特性、傳熱傳質、氣泡特性以及冷卻效果采用FLUENT軟件結合用戶自定義控制方程進行著重研究分析,同時結合實驗驗證,為高熱流電子設備蒸發冷卻換熱器的設計提供了指導。

1 沸騰傳熱控制理論

沸騰傳熱過程因流動以及傳熱傳質的復雜性,使得其發展處于半經驗半理論甚至更多依靠經驗的狀態,而實驗又很難研究流動的內部狀態,數值研究相應地彌補了這方面的不足。其中,最復雜部分是在沸騰過程中氣相與液相的傳熱傳質問題上。沸騰管內流動狀態受管壁面微觀結構、氣泡形態以及氣泡尾流或氣相二次流等因素的影響,本文在研究三維流道的換熱性能時忽略了這些因素,以保證數值研究在三維實體模型中的可實現性,同時假設氣泡為規則球體且氣泡直徑為定值、氣液相界面光滑。本文從微元層面分析某一網格單元體內的傳熱傳質情況,從而模擬計算三維流道內的流動及換熱效果。

1.1 氣泡動力特性

1.1.1 相間傳熱傳質 假設在某一網格區域產生某一球狀氣泡,氣泡的平面截圖如圖1所示,圓弧為相界面。VOF模型采用線性插值取代實際相界面,當網格步長遠小于氣泡直徑時,這種相界面處理方法產生的數值誤差是非常小的。

圖1 VOF模型氣液界面插值

氣相和液相受熱流驅使,質量和熱量的轉換是在球形界面上完成的。初始時刻,流動區域均為溫度不大于飽和溫度的過冷或飽和液體,受熱升溫后,貼近壁面部分流體溫度首先升高,形成熱邊界層。當邊界層內某一微小區域溫度大于對應條件下工質的飽和溫度時,可以認為液體“熱量多余”,這部分“多余熱量”需轉換為氣相的潛熱,從而使得兩相狀態的溫度值在飽和溫度附近,保證數值計算結果與實際情況相符。這部分的轉換方式是通過源項形式實現的,因而問題的復雜性便落到如何控制氣相和液相之間熱量和質量的轉換上。

基于前文假設,某一光滑相界面在特定條件下的傳輸質量由Hertz Knudsen方程給出

(1)

式中:M為流體的相對分子質量;R為氣體常數;P*為相界面上氣相一側局部壓力。

數值計算控制方程一般采用溫度值作為控制變量,式中含有壓力不利于數值研究。這里的傳質問題是發生在飽和狀態點附近,Clausius Claperyron方程描述為

(2)

當P*足夠靠近飽和狀態點時

(3)

此時,三維模型的單位體積質量源項為

(4)

實際流動過程中流道內壓降比較小,為了方便數值研究,將流體飽和溫度設為定值。令

(5)

在同一狀態條件下,γ為定值常數。氣泡直徑往往難以準確獲得,實際操作中發現,針對R134a或R22工質的沸騰傳熱研究,非惡劣工況下γ的取值在100~500之間,數值計算結果更貼近實驗觀察。

此時,微小單元內的潛熱與顯熱的轉換能量SE=ΔhS,SE為能量控制方程內的能量源項。

1.1.2 氣泡動力特性 從傳熱角度看,過冷流體內的氣泡生長過程主要受壁面熱流和Manrangoni熱流影響[10],流體的表面張力隨著溫度的升高而降低,在熱邊界層內會產生表面張力梯度,氣液相界面產生Manrangoni驅動力會抑制氣泡生長。當壁面過熱度進一步提高,壁面熱流遠高于Manrangoni熱流時,氣泡急劇生長,流動演化為核態沸騰狀態,傳熱系數和熱流密度急劇增大。圖2中Fm為Manrangoni驅動力,Fe為表面張力,共同抑制氣泡生長;Fb為氣泡在熱壁面高熱流的驅動下產生的氣泡膨脹力,促使氣泡成長。

圖2 氣泡受力分布圖

此外,氣泡在生長過程受到的其他力包括氣液相界面壓力波動產生的壓力、促使氣泡脫離壁面的力,以及體積力和實際管道流動中產生的一系列暫不可知的復雜力,對于這些力在本文分析氣泡動力特性中均忽略不計。

1.2 數值計算控制方程

為了能更好地數值研究氣泡的流動特性,本文采用追蹤相界面的方法,即利用VOF模型研究沸騰換熱流動。VOF模型處理相界面采用線性插值取代實際界面的方法,這樣的處理方式使得數值計算的相界面與實際界面吻合程度較高,但會有微小程度的數值擴散,尤其當網格單元步長與氣泡尺寸相當或更大時,這種擴散使得計算結果質量不守恒,產生錯誤的結果,因而網格單元步長需嚴格控制。考慮計算機性能限制,一般網格步長取0.5~1 mm為宜,微細通道研究應選得更小。

正是由于這種追蹤相界面的處理方法,使得研究對象質量守恒方程是通過求解每一相的體積分數的連續性方程實現的。本文研究對象為液相和氣相,其中液相和氣相體積分數方程分別為

(6)

(7)

當液相蒸發為氣相時,有一部分流體的顯熱轉化為氣相的潛熱;當氣相冷凝為液相時,勢必有一部分氣相的潛熱轉化為流體的顯熱。計算中流體表征溫度參數體現為顯熱,氣液焓差表征為潛熱。潛熱無法直接體現在能量方程中,即而采用源項形式處理顯熱與潛熱之間的轉化問題。能量控制方程為

(8)

SE的作用體現為潛熱與顯熱的轉化。

VOF模型的動量方程是對混合相求解單一方程實現的,動量的損失體現為壓力的損失,本文研究中因壓降較小,故忽略動量損失。動量方程為

(9)

式中:F為動量源項,體現為表面張力、壁面黏附力對動量的影響。

2 物理模型

實驗研究中,采用發熱電阻貼合在具有蛇形流道的冷板上來模擬雷達天線陣面的發熱及降溫過程。如圖3a所示,冷板與水平面的夾角記為θ,用以模擬雷達天線在轉動過程中重力加速度方向與流動方向的變化關系。本文選取θ=0°,45°,90°,135°,180°共5種工況作針對性研究。圖3b給出了10個加熱單元的位置編號。

冷板流道區域采用1 mm步長的結構體網格,固體冷板區域采用分塊處理的1 mm步長結構體網格,網格總數在250萬左右。工質采用R134a,物性參數以NIST查詢數據為基礎,采用小區間線性插值方式嵌入FLUENT材料庫。工作壓力為1.68 MPa,飽和溫度為333 K,進口過冷度為5K,流量為0.1 m3/h,環境溫度為293 K,每個加熱單元功率為200 W,熱流密度約800 kW/m2。計算模型采用VOF模型,同時考慮Marangoni效應。

(a)冷板結構示意圖

(b)加熱單元編號

3 結果分析

3.1 整體分析

過冷工質從進口流入冷板換熱器,在第一個加熱單元附近溫度升至飽和溫度,溫度梯度較大,繼續吸熱產生相變之后形成兩相沸騰流動。在工況1~3下,當壁面受熱產生氣泡之后,由于相間力不足以克服重力的作用,氣相在重力加速度影響下,在貼近發熱單元一側聚集,形成分層流動,使得高熱流壁面直接與氣相接觸,從一定程度上抑制了蒸發傳熱。圖4a為形成層狀流動的氣相體積分數φ分布圖。在工況4、5下,氣泡在壁面產生之后,隨即上浮至遠離加熱單元一側,使得氣相能與主流充分進行質、熱交換,形成泡狀流,充分保證了高熱流壁面與冷流體(液相)接觸換熱,換熱性能較好。圖4b為形成泡狀流的氣相體積分數分布圖。

圖4b取自蛇形管道某一連續彎道的直管段處。本文研究發現,氣液兩相流體在流經彎道時并不是均勻的層狀流,氣液相微小單元在相間作用力以及表面張力等綜合作用下流經彎道時,氣相會在彎道處匯聚。當通過連續彎道時,氣相主流區域從上一彎道內壁附近流向下一彎道內壁區域。如圖4c所示,陰影部分是氣相主流示意圖,在連續彎道產生的擾流作用下,氣相主流并不是沿某個截面直線流動。圖4b為沿管道流動中心的垂直截面圖,氣相體積分數分布產生了不連續的現象,這種匯聚趨勢也正是下文優化分析的著重點。模擬結果同時顯示,在過冷度同為5K的情況下,工況4、5的初始沸騰時間要比工況1~3長近4 s,工況4、5的熱邊界層厚度要明顯小于工況1~3。

(a)層狀流示意圖

(b)泡狀流示意圖

(c)彎道氣相匯聚示意圖

為分析不同加熱單元附近的換熱效果,取工況1、2為研究對象,分析流經每一個加熱單元的蒸發量情況,即流經加熱單元后氣相體積分數的變化值Δφ。圖5給出了工質流經每個加熱單元時Δφ的變化趨勢,沿流動方向整體呈上升趨勢,在出口附近有微小差異。工質流經每個加熱單元時,Δφ均有所增加,但增加幅度有所不同。

文獻[11]指出,氣液兩相流動的傳熱性能是核態傳熱和對流傳熱的綜合作用。氣相的導熱性能差,阻礙換熱,但其流動速度比液相大,對流換熱性能好,總體換熱性能是這兩種因素的綜合作用。加熱單元10附近對應流道內氣相體積分數較大,氣膜較厚,氣膜直接與高熱流壁面接觸阻礙了液相的蒸發,同時加熱單元10為最后一個彎道,流動下游產生的擾流作用弱于連續彎道流動,故加熱單元10的總體蒸發量小于其他地方。蛇形管流道氣相體積分數分布圖(圖4b)顯示,在彎道附近容易產生氣相的匯聚、漩渦,形成流動死區,從而使得一部分氣相直接與高熱流壁面長時間接觸,產生局部過熱現象,換熱效果變差,過熱程度尤以工況3更為惡劣。

圖5 不同加熱單元的Δφ趨勢圖

3.2 優化分析

本文研究顯示:當熱邊界層內產生的氣相能較充分地與主流進行傳質傳熱時,換熱性能較好;反之,氣相在熱邊界層內的匯聚部分將抑制換熱。為方便研究,截取蛇形管流道某一處彎道作詳細分析,重力加速度方向以工況3為準。

(a)截取彎道氣相體積分數分布圖

(b)截取彎道溫度分布圖

圖6a、6b分別是截取部分的氣相體積分數以及溫度分布圖。通過流動動畫顯示:圖中A區域產生流動漩渦,氣相在漩渦中匯聚,使得這部分高熱流壁面直接與氣相接觸,換熱效果變差;B區域是在重力加速度作用下形成的分層流動的氣相匯聚。A、B區域溫度局部最大值為333.9 K,局部過熱度為0.9 K。氣相微小單元在壁面產生后,在重力、相間力等各種綜合力作用下形成漩渦,產生流動死區,不能很好地匯入主流從而使得換熱效果變差。若能使得主流對氣相的作用力足以克服漩渦產生的綜合力作用,那么此處的匯聚情況便可以得到改善。

為此,對彎道作減小過流面積處理,圖7a、7b分別是改進模型的數值模擬結果,結果顯示圖6a中A區域的流動漩渦消失,且A區域溫度降至飽和溫度,局部過熱消失。A區域緊貼著加熱單元中心區域,熱流密度較高,當這部分區域的壁面與液相充分接觸時,則可改善換熱效果。

(a)改進模型氣相體積分數分布圖

(b)改進模型溫度分布圖

4 實驗對比

本文主要從氣液相流動特性以及流動的宏觀趨勢來研究沸騰換熱效果,由于流態的實驗驗證在目前技術條件下還很難進行,因此從加熱單元的溫度值驗證了整個數值模擬系統的可信性。實驗冷板水平放置于兩相流實驗平臺,工質由工質泵供給,流經干燥過濾器、流量計,進入實驗冷板。冷板進出口均有試液鏡觀察流動情況,實驗觀察到冷板出口確為兩相狀態,且為層狀流,與工況1的計算結果一致,隨后工質流入冷凝器,最終進入儲液罐,完成循環。

實驗驗證工況下采用R22作為工質,流量為0.618 t/h,進口過冷度約為2 K,加熱單元總功率為1 010 W,溫度為常溫,壓力在1 MPa左右。圖8為加熱單元溫度實驗值與計算值對比圖。加熱單元1、2的溫度計算值與實驗值相差甚微,此處近似為單相管內流動;沸騰現象產生之后,計算值與實驗值總體趨勢基本一致,但均低于實驗值5K左右,在兩相蒸發的數值計算中屬于可接受的范圍,計算結果可以用于冷板沸騰換熱研究。數值研究中,加熱單元、擴熱板以及冷板之間的傳熱都是理想化的,而實際實驗過程中傳熱存在接觸熱阻,必然導致加熱單元的總體溫度略高。由于實驗管道中還存在少量不凝結氣體以及雜質等,使得換熱效果變差,對蒸發換熱性能影響比較大,而計算工況都是理想化的純工質,因此計算工況的冷卻效果好于實驗工況,加熱單元溫度值總體偏低。

圖8 加熱單元溫度計算值與實驗值對比圖

5 結 論

(1)沸騰傳熱模型中傳熱傳質控制量歸結為γ的選擇。一般地,對于R134a、R22工質在較高熱流密度下,γ選擇在500左右較為適宜,可依據網格步長以及熱流大小適當地增大或減小γ值。控制方程中,同一狀態條件下,影響傳質的因素有氣泡直徑和調節系數β,這些因素均反映在γ的選取上。

(2)VOF模型能較好地模擬計算三維流道內的分層流動和泡狀流動,從一定程度上反映兩相流動的氣液分、聚情況,為改進流道以提高換熱系數提供指導。網格步長以及時間步長需慎重選擇,否則將會造成質量的數值擴散。步長若與氣泡尺寸相當或者更大,則將產生較大的數值擴散,導致計算結果質量嚴重不守恒,致使計算結果錯誤。

(3)對于彎道處的氣相匯聚問題,通過減小過流面積提高局部流速的方法可以得到緩解,消除漩渦,從而提高換熱性能。加熱單元溫度值的計算與實驗結果差異較小,宏觀趨勢基本一致,數值方法可以用于冷板沸騰換熱研究。

[1] 姚濤, 候哲, 顧國彪.蒸發冷卻技術應用于大型汽輪發電機的技術可行性 [J].電工技術學報, 2008, 23(2): 1-5.

YAO Tao, HOU Zhe, GU Guobiao.Application of evaporative cooling technique to large turbo generator [J].Transactions of China Electrotechnical Society, 2008, 23(2): 1-5.

[2] 熊楠.蒸發冷卻技術面對的幾個重要問題 [J].電工電能新技術, 2001, 3(1): 36-40.

XIONG Nan.Main tasks for research and practice of evaporation-cooling technology [J].Advanced Technology of Electrical Engineering and Energy, 2001, 3(1): 36-40.

[3] 郭雷.微細通道流動沸騰換熱機理及實驗研究 [D].濟南: 山東大學, 2011.

[4] 曾陽.蒸發冷卻換熱器的數值模擬與實驗研究 [D].長沙: 湖南大學, 2011.

[5] 姜澎, 黃洪雁, 馮國泰, 等.水霧/蒸汽相變沖擊冷卻的數值模擬 [J].哈爾濱工程大學學報, 2009, 30(10): 1097-1101.

JIANG Peng, HUANG Hongyan, FENG Guotai, et al.Numerically simulating mist/steam impingement cooling with phase change [J].Journal of Harbin Engineering University, 2009, 30(10): 1097-1101.

[6] ZHUAN Rui, WANG Wen.Simulation of subcooled flow boiling in a micro-channel [J].Refrigeration International Journal of Refrigeration, 2010, 34(3): 781-795.

[7] ZHUAN Rui, WANG Wen.Simulation on nucleate boiling in micro-channel [J].Heat Mass Transfer, 2010, 53(1): 502-512.

[8] THOME J R, DUPONT V, JACOBI A M.Heat transfer model for evaporation in micro channels [J].Heat Mass Transfer, 2004, 47(14): 3375-3385.

[9] LEE W, SON G.Numerical simulation of boiling enhancement on a micro-structured surface [J].Heat Mass Transfer, 2011, 38(2): 168-173.

[10]ZHANG N L, CHAO D F.Models for enhanced boiling heat transfer by unusual Marangoni effects under microgravity conditions [J].Heat Mass Transfer, 1999, 26(8): 1081-1090.

[11]RICARDO J, LIMA D S, THOME J R, et al.Two-phase flow boiling in horizontal smooth tubes: new heat transfer results for R134a at three saturation temperatures [J].Applied Thermal Engineering, 2009, 29(7): 1289-1298.

(編輯 荊樹蓉)

NumericalInvestigationandOptimizationonBoilingHeatExchangeofAntennaSurface

YIN Xiang1,QIAN Jiyu2,KONG Xiangju2,WU Huanggen1,CAO Fen1

(1.School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China;2.The Fourteenth Institution of China Electronics Technology Group Corporation, Nanjing 210039, China)

The evaporative cooling of two-phase model was established to cool the antenna surface with partial high heat flux.The analysis was focused on the ability of cooling under the local high heat flux of 800 kw/m2.Using VOF model, as well as the governing equation defined by the user, the computation of the boiling heat exchange in the three-dimensional pipe was performed to study the reason of hindering the heat exchange from the flow of the vapor and liquid.Results show that, the model could be well used in the analysis of stratified flow, bubble flow and the flow trend of vapor and liquid in the two-phase flow.Decreasing the flow area of the curved tube and increasing the local velocity could reduce the converging of the vapor.And the local superheat of 0.9 K was dismissed.The coefficient in the evaporation control equation should be increased to about 500 to ensure the accuracy of simulation in the partial high flux model.The overall error of the numerical result was about 5K but the tendency of the temperature was similar with the experiment.So, numerical analysis could be used in the study of the boiling heat exchange for the cold plate.

model of VOF; evaporative cooling; boiling heat exchange; numerical investigation

2014-03-23。

殷翔(1991—),男,碩士生;吳華根(通信作者),男,副教授。

時間:2014-09-01

10.7652/xjtuxb201411011

TM216.1

:A

:0253-987X(2014)11-0064-06

網絡出版地址:http:∥www.cnki.net/kcms/detail/61.1069.T.20140901.1009.004.html

猜你喜歡
研究
FMS與YBT相關性的實證研究
2020年國內翻譯研究述評
遼代千人邑研究述論
視錯覺在平面設計中的應用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
關于遼朝“一國兩制”研究的回顧與思考
EMA伺服控制系統研究
基于聲、光、磁、觸摸多功能控制的研究
電子制作(2018年11期)2018-08-04 03:26:04
新版C-NCAP側面碰撞假人損傷研究
關于反傾銷會計研究的思考
焊接膜層脫落的攻關研究
電子制作(2017年23期)2017-02-02 07:17:19
主站蜘蛛池模板: 少妇人妻无码首页| 色老二精品视频在线观看| 国产成人精品亚洲77美色| 国产第一页第二页| 四虎国产精品永久一区| 综合久久五月天| 在线观看av永久| 中国国语毛片免费观看视频| 中文字幕色在线| 国产精品视频猛进猛出| 波多野结衣中文字幕一区| 欧美日韩免费在线视频| 欧美成人综合视频| 亚洲天天更新| 97se亚洲综合| 亚洲中文精品人人永久免费| 亚洲精品日产精品乱码不卡| 无码视频国产精品一区二区| 久久综合丝袜长腿丝袜| 国产乱子伦无码精品小说| 国产又色又刺激高潮免费看 | 幺女国产一级毛片| 午夜限制老子影院888| 欧美成一级| 日韩免费毛片视频| 19国产精品麻豆免费观看| 国产三区二区| 亚洲天堂成人在线观看| 国产大片黄在线观看| 久久99国产精品成人欧美| 四虎永久在线精品国产免费 | 精品综合久久久久久97超人该| 无码在线激情片| 无码内射中文字幕岛国片| 2021精品国产自在现线看| 国产精品久久久久久久久| a级毛片视频免费观看| 国产精品视频公开费视频| 日本人妻一区二区三区不卡影院 | 丁香亚洲综合五月天婷婷| 国产91全国探花系列在线播放| 一本久道热中字伊人| 日本成人精品视频| 国产无码网站在线观看| 91福利在线观看视频| 巨熟乳波霸若妻中文观看免费| 日本免费精品| 欧美日韩午夜| 欧美成一级| 久久香蕉欧美精品| 91国内外精品自在线播放| 人妻21p大胆| 成人一级免费视频| 亚洲一区二区三区麻豆| a亚洲视频| 国产永久在线观看| 亚洲综合色婷婷中文字幕| 日韩A∨精品日韩精品无码| 亚洲日本中文字幕天堂网| 5555国产在线观看| 成人午夜久久| 色国产视频| 国产精彩视频在线观看| 国产va在线观看免费| 亚洲中文字幕久久无码精品A| 亚洲精品无码人妻无码| www中文字幕在线观看| 国产精品网址你懂的| 日韩a在线观看免费观看| 国产成人在线无码免费视频| 久久精品丝袜高跟鞋| 美女一区二区在线观看| 亚洲区第一页| 国产精品久久久精品三级| 在线国产你懂的| 国产精品第一区在线观看| 99国产精品一区二区| 国产精品私拍99pans大尺度| 国产成人AV大片大片在线播放 | 国产91丝袜| 国产99在线观看| 国产伦精品一区二区三区视频优播|