• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    交界面剪切力對(duì)傾斜扁平管內(nèi)蒸汽凝結(jié)影響的CFD模擬

    2016-07-04 03:43:20鄧慧白焰李欣欣張東明華北電力大學(xué)控制與計(jì)算機(jī)工程學(xué)院北京006中國(guó)電能成套設(shè)備有限公司北京00080北京國(guó)電智深控制技術(shù)公司北京000
    化工學(xué)報(bào) 2016年4期
    關(guān)鍵詞:數(shù)值模擬湍流

    鄧慧,白焰,李欣欣,張東明(華北電力大學(xué)控制與計(jì)算機(jī)工程學(xué)院,北京 006;中國(guó)電能成套設(shè)備有限公司,北京 00080;北京國(guó)電智深控制技術(shù)公司,北京 000)

    ?

    交界面剪切力對(duì)傾斜扁平管內(nèi)蒸汽凝結(jié)影響的CFD模擬

    鄧慧1,白焰1,李欣欣2,張東明3
    (1華北電力大學(xué)控制與計(jì)算機(jī)工程學(xué)院,北京 102206;2中國(guó)電能成套設(shè)備有限公司,北京 100080;3北京國(guó)電智深控制技術(shù)公司,北京 102200)

    摘要:針對(duì)傾斜扁平管內(nèi)的湍流蒸汽,考慮交界面剪切力,建立蒸汽凝結(jié)的CFD模型。在一根實(shí)驗(yàn)樣管上(長(zhǎng)寬高尺寸:2600 mm×3 mm×50 mm,傾斜角:60°),對(duì)比管內(nèi)平均凝結(jié)傳熱系數(shù)和蒸汽凝結(jié)率的CFD解與實(shí)驗(yàn)值,驗(yàn)證模型的有效性。計(jì)算不同蒸汽入口流速下的交界面剪切力值,發(fā)現(xiàn)剪切力的大小由蒸汽入口流速?zèng)Q定,其數(shù)值隨蒸汽流動(dòng)持續(xù)遞減。對(duì)剪切力對(duì)蒸汽凝結(jié)影響進(jìn)行CFD模擬,發(fā)現(xiàn),剪切力增大管內(nèi)局部凝結(jié)傳熱系數(shù),減小蒸汽凝結(jié)質(zhì)量;在0~0.8 m管段,剪切力明顯破壞了重力對(duì)液膜的累積,削薄了液膜厚度,但從1.0 m到蒸汽出口,液膜厚度由重力控制,剪切力的影響可忽略不計(jì);剪切力迫使液膜加速流動(dòng),其加速作用0~0.2 m管段表現(xiàn)突出。

    關(guān)鍵詞:湍流;凝結(jié);界面剪切力;傾斜扁平管;數(shù)值模擬

    2015-06-24收到初稿,2015-10-14收到修改稿。

    聯(lián)系人:白焰。第一作者:鄧慧(1979—),男,博士研究生。

    Received date: 2015-06-24.

    引 言

    空冷凝汽器由于節(jié)水效果好,在我國(guó)“富煤缺水”地區(qū)的火電機(jī)組上得到了廣泛應(yīng)用[1-2]??绽淠髌毡槭褂帽馄缴咝纬峁茏鰮Q熱器的基管,由于基管蒸汽入口Reynolds數(shù)較高,蒸汽在管內(nèi)處于湍流狀態(tài),而汽液交界面剪切力對(duì)湍流狀態(tài)下的蒸汽凝結(jié)有重要影響。

    Nusselt[3]研究純飽和蒸汽在垂直平板上的膜狀凝結(jié),忽略交界面剪切力,提出邊界層模型。Chen 等[4]研究?jī)上嚅]式虹吸熱管內(nèi)的蒸汽凝結(jié),發(fā)現(xiàn)交界面剪切力增加了液膜厚度,減小了管內(nèi)蒸汽凝結(jié)傳熱系數(shù)(HTC)。Fiedler等[5]實(shí)驗(yàn)研究R134a制冷劑在傾斜圓管內(nèi)的逆流凝結(jié),結(jié)果表明,剪切力對(duì)R134a凝結(jié)的影響可忽略不計(jì)。Rohsenow等[6]研究垂直平板上的蒸汽凝結(jié),發(fā)現(xiàn)剪切力加劇了蒸汽凝結(jié),由此提出:用量綱1剪切力系數(shù)修正蒸汽凝結(jié)HTC計(jì)算。Lee等[7]實(shí)驗(yàn)研究垂直圓管內(nèi)的混合氣體凝結(jié),結(jié)果表明,剪切力會(huì)削弱不可凝氣體對(duì)蒸汽凝結(jié)的不利影響。張俊霞等[8]用VOF (volume of fluid)模型[9]研究垂直圓管內(nèi)的蒸汽凝結(jié),發(fā)現(xiàn)在蒸汽入口區(qū),剪切力對(duì)凝結(jié)有明顯的強(qiáng)化作用。

    在交界面剪切力的計(jì)算上,基于Couette流理論[10]的近似經(jīng)驗(yàn)公式[11-13]得到了廣泛應(yīng)用,Oh等[14]和Zhou等[15]還分別設(shè)計(jì)迭代算法計(jì)算剪切力。

    在以上文獻(xiàn)中,流體的凝結(jié)發(fā)生在垂直圓管內(nèi)或二維平壁上,圓管和平壁都具有標(biāo)準(zhǔn)的幾何特征。針對(duì)非標(biāo)準(zhǔn)幾何特征的傾斜扁平管,缺少利用CFD模擬剪切力對(duì)蒸汽凝結(jié)和管內(nèi)流場(chǎng)影響的相關(guān)報(bào)道。

    本文考慮蒸汽湍流狀態(tài)下的氣液交界面剪切力,提出模擬傾斜扁平管內(nèi)蒸汽凝結(jié)的數(shù)學(xué)模型;分別用Ansys Fluent?和Wolfram Mathematica?求解數(shù)學(xué)模型中的蒸汽控制方程和液膜控制方程;設(shè)計(jì)迭代算法,計(jì)算兩相控制方程間的耦合變量即交界面剪切力;用Ansys Fluent?數(shù)值模擬剪切力對(duì)蒸汽凝結(jié)和管內(nèi)流場(chǎng)的影響。

    1 研究對(duì)象

    杜小澤等[16]用單排大扁管釬焊蛇形翅片做管束基管,設(shè)計(jì)2排管束(各155根翅管) 以60°角對(duì)立傾斜,組成Λ形換熱器,換熱器下布置軸流風(fēng)機(jī),構(gòu)成空冷凝汽器樣機(jī)。樣機(jī)的最大設(shè)計(jì)蒸汽流量為0.28 kg·s-1,正常工作狀態(tài)下蒸汽流量0.17 kg·s-1。

    本文以文獻(xiàn)[16]空冷凝汽器樣機(jī)的翅管為數(shù)值模擬對(duì)象,翅管及其管束的形狀如圖1所示。

    圖1 翅管和管束的形狀Fig.1 Pictorial view of plate wave-finned tube and bundle

    2 數(shù)值計(jì)算方法

    2.1凝結(jié)假設(shè)條件

    在Nusselt理論[3]基礎(chǔ)上,通過以下假設(shè)條件,建立液膜和蒸汽的流體控制方程,得到蒸汽凝結(jié)的數(shù)學(xué)模型。

    (1)液相假設(shè)。凝結(jié)液膜為不可壓縮黏性流體,定常層流表面有波動(dòng);忽略液膜運(yùn)動(dòng)方程中的慣性項(xiàng)及能量方程中的對(duì)流項(xiàng)。

    (2)氣相假設(shè)。蒸汽為穩(wěn)態(tài)可壓縮流體,在管內(nèi)定常湍流,湍流特性滿足標(biāo)準(zhǔn)k-ε方程;管內(nèi)蒸汽凝結(jié)為純蒸汽的膜狀凝結(jié),無(wú)不可凝氣體影響;在流動(dòng)凝結(jié)過程中,忽略蒸汽管內(nèi)流動(dòng)的壓損,認(rèn)為蒸汽狀態(tài)處于均勻的飽和溫度Ts,蒸汽中無(wú)溫度梯度,不計(jì)蒸汽凝結(jié)過程中釋放的顯熱。

    (3)液膜熱邊界條件假設(shè)。液膜表面溫度等于蒸汽飽和溫度;蒸汽的膜狀凝結(jié)只在大扁管的平壁上發(fā)生,平壁為溫度待定的等溫壁,大扁管上下小半橢圓內(nèi)壁上無(wú)蒸汽凝結(jié),假設(shè)為絕熱壁。

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

    基于圖1中管束的傾斜角和大扁管的幾何形狀,將大扁管的管高、管寬和軸向管深方向,定義為局部坐標(biāo)系的x、y和z方向,得到大扁管平壁上的凝結(jié)液膜微分體如圖2所示。

    圖2 局部坐標(biāo)系下的液膜微分控制體Fig.2 Differential control volume of condensate film under local Cartesian coordinates

    由2.1節(jié)假設(shè)條件,得液膜的Navier-Stokes方程和邊界條件

    蒸汽的連續(xù)和Navier-Stokes方程[17]為

    以標(biāo)準(zhǔn)k-ε方程為蒸汽的湍流方程(8)

    蒸汽在液膜表面上釋放的潛熱,由液膜的導(dǎo)熱傳至管外,蒸汽凝結(jié)失去的質(zhì)量由液膜獲得,蒸汽和液膜在交界面上的剪切力相等。管內(nèi)氣液兩相流體在交界面(即液膜表面)上傳熱傳質(zhì)的連續(xù)和剪切力耦合,用以下液膜表面微分體的質(zhì)量、剪切力和熱量連續(xù)條件描述。

    質(zhì)量

    剪切力

    由于管內(nèi)蒸汽在管內(nèi)的軸向流速遠(yuǎn)大于橫向流速,忽略剪切力管高方向上的分量,式(9)簡(jiǎn)化為

    熱量

    將連續(xù)條件式(10)代入液膜動(dòng)量方程的邊界條件,得液膜速度分布函數(shù)

    式(12)反映了蒸汽壓力、重力和交界面剪切力對(duì)液膜流速的影響。由式(12)可計(jì)算液膜沿管高和管深方向的質(zhì)量流量

    將連續(xù)條件式(11)中的交界面溫度代入液膜熱量方程的邊界條件,可得液膜溫度梯度,由Fourier導(dǎo)熱定律,考慮液膜表面波動(dòng)效應(yīng)[18],對(duì)熱導(dǎo)率修正,得液膜微分體的熱通量

    將式(15)代入連續(xù)條件式(11)中的等熱流方程,得氣液交界面換熱平衡方程

    聯(lián)立連續(xù)條件式(9)和式(13)~式(16),得傾斜大扁管內(nèi)壁上的液膜厚度偏微分方程

    由連續(xù)條件式(9),聯(lián)立式(13)和式(14),可計(jì)算蒸汽在微分體表面上的冷凝量,即為蒸汽控制方程的質(zhì)量源項(xiàng)

    在CFD計(jì)算中,式(18)中A取大扁管內(nèi)壁網(wǎng)格面積,V取與內(nèi)壁貼面的蒸汽網(wǎng)格體積。由質(zhì)量源項(xiàng),可計(jì)算蒸汽方程的動(dòng)量和熱量源項(xiàng)

    2.3交界面剪切力的計(jì)算

    蒸汽和液膜的控制方程之間,存在耦合變量即交界面剪切力,兩套方程須聯(lián)立求解,設(shè)計(jì)迭代算法求解交界面剪切力。

    (1)輸入大扁管內(nèi)壁溫度。

    (2)假設(shè)交界面剪切力分布初值,構(gòu)建剪切力的插值函數(shù),用Wolfram Mathematica?求解液膜厚度偏微分方程式(17)。

    (3)將大扁管內(nèi)壁設(shè)置成絕熱壁,由式(18)~式(20)計(jì)算蒸汽方程的源項(xiàng),用Ansys Fluent?求解蒸汽的控制方程和湍流方程式(4)~式(8)。

    (4)用Ansys Fluent?的UDF(User Defined File)提取蒸汽在幾何邊界上的剪切力,與步驟(2)初值比較,循環(huán)執(zhí)行步驟(2)~步驟(4)直至剪切力差值收斂。

    (5)由液膜厚度偏微分方程的數(shù)值解,計(jì)算蒸汽凝結(jié)放熱量和液膜流量。

    為簡(jiǎn)化計(jì)算,計(jì)算蒸汽方程時(shí),忽略液膜厚度,將大扁管內(nèi)壁設(shè)置為蒸汽流場(chǎng)的幾何邊界,且蒸汽只在內(nèi)壁的貼面網(wǎng)格內(nèi)冷凝。

    Ansys Fluent?使用標(biāo)準(zhǔn)笛卡爾全局坐標(biāo)系,在Wolfram Mathematica?中,使用圖2所示的局部坐標(biāo)系,兩個(gè)坐標(biāo)系下的向量場(chǎng)通過歐拉旋轉(zhuǎn)矩陣轉(zhuǎn)換[19]。迭代計(jì)算中,用Mathlink?在計(jì)算軟件之間交換界面剪切力。

    2.4翅管的翅側(cè)計(jì)算

    由于管壁很薄(厚度為1 mm的不銹鋼),可認(rèn)為管內(nèi)外壁溫相等。本文假設(shè)大扁管為溫度未知的等溫壁,異步計(jì)算翅管兩側(cè)換熱平衡,以估算管壁溫度。具體方法為:設(shè)置壁溫條件從Ts-0.25 K到Ts-0.05 K,以0.001 K為步長(zhǎng),用Ansys Fluent?計(jì)算翅側(cè)換熱,用2.2節(jié)中的迭代算法計(jì)算管內(nèi)換熱,比較翅側(cè)和管內(nèi)的計(jì)算結(jié)果,求得兩側(cè)換熱差值的最小值,當(dāng)最小值滿足誤差要求時(shí),管內(nèi)和管外達(dá)到熱平衡狀態(tài)。

    在用Ansys Fluent?計(jì)算翅側(cè)換熱時(shí),假設(shè)所有翅片空氣入口條件一致,根據(jù)實(shí)驗(yàn)樣機(jī)入口空氣流量計(jì)算并設(shè)置空氣入口為速度入口。取翅管單側(cè)的2個(gè)翅片做數(shù)值模擬,得到空氣換熱量后在乘以564作為單根翅管的翅側(cè)換熱值(單根翅管單側(cè)共564個(gè)翅片)。

    驗(yàn)證翅側(cè)計(jì)算網(wǎng)格的獨(dú)立性時(shí)發(fā)現(xiàn),隨著網(wǎng)格數(shù)從260642到781926變動(dòng),空氣在翅片出口處溫度和速度的偏差在2%內(nèi),CFD最后采用260642個(gè)網(wǎng)格計(jì)算。翅側(cè)數(shù)值計(jì)算的控制方程、邊界條件和方法可參考文獻(xiàn)[20],數(shù)值仿真得到翅側(cè)冷卻空氣的溫度場(chǎng)如圖3所示。

    圖3 翅側(cè)換熱空氣的溫度瀑布圖Fig.3 Numerical flood illustration for air temperature of fin-side on tube

    2.5蒸汽側(cè)網(wǎng)格劃分

    用Ansys ICEM?建立大扁管的幾何模型并劃分計(jì)算網(wǎng)格。管內(nèi)蒸汽區(qū)采用O網(wǎng)格如圖4所示,為精確模擬蒸汽在管壁上的凝結(jié),網(wǎng)格在管壁貼面區(qū)域加密,形成邊界層。采用Bigeometric劃分大扁管平壁網(wǎng)格如圖5所示,網(wǎng)格單元的尺寸沿管深和管寬方向單向增大。

    圖4 蒸汽區(qū)域O網(wǎng)格沿管深剖面Fig.4 Plate view in tube axial direction for O-grid of vapour zone

    圖5 大扁管平壁的Bigeometric網(wǎng)格Fig.5 Bigeometric grid of condensation wall on flat tube

    隨著網(wǎng)格數(shù)從292200、686928到925600變化,蒸汽在管壁上剪切力和速度的偏差在2%內(nèi)變化,可認(rèn)為CFD的數(shù)值解獨(dú)立于網(wǎng)格數(shù),CFD最后采用292200個(gè)網(wǎng)格計(jì)算。

    2.6蒸汽側(cè)計(jì)算的條件和方法

    2.6.1邊界條件設(shè)置蒸汽進(jìn)口velocity-inlet,蒸汽出口outflow,其他邊界均設(shè)置為絕熱壁。

    2.6.2湍流設(shè)置蒸汽湍流方程設(shè)置如下

    2.6.3計(jì)算方法壓力與速度耦合采用SIMPLE算法,壓力項(xiàng)離散采用Standard格式,其他項(xiàng)離散為二階Upwind格式。數(shù)值結(jié)果判斂標(biāo)準(zhǔn)為:交界面剪切力殘差小于10-3,方程殘差小于10-4。

    3 結(jié)果與討論

    3.1模型驗(yàn)證

    針對(duì)空冷樣機(jī)蒸汽凝結(jié)率和管內(nèi)凝結(jié)平均傳熱系數(shù)(HTC),對(duì)比凝結(jié)數(shù)學(xué)模型的CFD解與文獻(xiàn)[16]實(shí)驗(yàn)數(shù)據(jù),驗(yàn)證模型的有效性。

    空冷凝汽器樣機(jī)的蒸汽凝結(jié)量用單根基管凝結(jié)量乘以基管總數(shù)計(jì)算。單根基管的蒸汽凝結(jié)量,為凝結(jié)液膜沿管高和管深方向質(zhì)量流量之和,取Wolfram Mathematica?對(duì)以下兩積分式的數(shù)值解

    管內(nèi)平均凝結(jié)HTC,取Wolfram Mathematica?對(duì)以下積分式的數(shù)值解

    表1是凝汽器樣機(jī)蒸汽凝結(jié)率的CFD解與文獻(xiàn)[16]實(shí)驗(yàn)表1結(jié)果對(duì)比,由表可見,數(shù)值解與實(shí)驗(yàn)值的誤差不超過4%,數(shù)值解與實(shí)驗(yàn)結(jié)果吻合良好,精確度較高。

    圖6是管內(nèi)平均凝結(jié)HTC的CFD數(shù)值解與文獻(xiàn)[16]實(shí)驗(yàn)表3結(jié)果對(duì)比,數(shù)值解與實(shí)驗(yàn)值的最大誤差不超過12%,數(shù)值解與實(shí)驗(yàn)結(jié)果匹配良好。

    圖6 蒸汽平均凝結(jié)HTC的CFD解與文獻(xiàn)[16]結(jié)果對(duì)比Fig.6 Comparison of mean condensation HTC between the CFD results and experimental results of Ref.[16]

    表1 蒸汽凝結(jié)率CFD解與文獻(xiàn)[16]實(shí)驗(yàn)結(jié)果對(duì)比Table 1 Comparison of condensate rate between CFD results and experimental results of Ref.[16]

    3.2交界面剪切力的數(shù)值模擬

    取文獻(xiàn)[16]表3實(shí)驗(yàn)5的工況條件,由空冷凝汽器正常工作蒸汽流量0.17 kg·s-1,計(jì)算對(duì)應(yīng)的基管蒸汽入口流速為4.86 m·s-1,用CFD數(shù)值模擬蒸汽流速4.86 m·s-1下的Z方向剪切力(剪切力在軸向管深方向上的分量),其瀑布線如圖7所示。由圖7可見,Z向剪切力分布基本不受湍流影響,數(shù)值沿蒸汽流動(dòng)持續(xù)遞減。在大扁管尾部,剪切力有很小的負(fù)值出現(xiàn),表示在蒸汽出口區(qū),由于大部分蒸汽凝結(jié)成水,蒸汽流速已低于液膜流速,同時(shí)受出口未凝結(jié)完蒸汽的回流影響,交界面剪切力對(duì)液膜有微弱的逆向剪切作用。

    在3組蒸汽流速下(4.86、7.29和9.72 m·s-1),分別計(jì)算剪切力Z分量沿管高的平均值,對(duì)比3組CFD解如圖8所示。由圖可見,蒸汽入口流速?zèng)Q定管內(nèi)交界面剪切力的大小,流速越高,剪切力越大。在蒸汽入口區(qū)非常短的距離內(nèi),3組剪切力都急劇下降一半,然后隨著蒸汽的流動(dòng),剪切力基本以線性速度遞減。

    圖7 蒸汽流速4.86 m·s-1下Z方向剪切力的瀑布線Fig.7 Flood and line of Z shear stress predicted by proposed model in vapour velocity 4.86 m·s-1

    圖8 不同蒸汽流速的Z方向剪切力平均值Fig.8 Average Z-shear stress vs axial distance from inlet by variation of vapour velocity

    3.3剪切力對(duì)蒸汽凝結(jié)的影響

    取文獻(xiàn)[16]表3實(shí)驗(yàn)5的工況條件,用4組CFD計(jì)算蒸汽的局部凝結(jié)HTC和局部凝結(jié)量,第1組計(jì)算忽略交界面剪切力,第2~4組計(jì)算剪切力,蒸汽速度為4.86、7.29和9.72 m·s-1。4組數(shù)值計(jì)算結(jié)果的對(duì)比,分別見圖9和圖10。

    圖9 對(duì)數(shù)-線性標(biāo)尺下不同蒸汽流速的局部凝結(jié)HTCFig.9 Local condensate HTC vs axial distance from inlet in log-line scale by variation of vapour velocity

    圖10 不同蒸汽流速的局部凝結(jié)量Fig.10 Local condensate rate vs axial distance from inlet by variation of vapour velocity

    由圖9可見,剪切力對(duì)凝結(jié)的強(qiáng)化作用在0~1 m管深區(qū)表現(xiàn)明顯,隨著蒸汽凝結(jié),剪切力對(duì)局部凝結(jié)HTC的放大作用逐漸減弱。由圖9還可見,蒸汽流速越大,局部凝結(jié)HTC越高,蒸汽流速7.29 m·s-1(凝汽器最大設(shè)計(jì)蒸汽流量對(duì)應(yīng)流速)和9.72 m·s-1下的局部凝結(jié)HTC曲線基本吻合,這表示:當(dāng)蒸汽流量超過設(shè)計(jì)流量時(shí),蒸汽流速(剪切力)對(duì)凝結(jié)的強(qiáng)化作用達(dá)到極限。

    由圖10可見,剪切力作用減小了蒸汽的局部凝結(jié)量。蒸汽流速越大,蒸汽在大扁管內(nèi)的凝結(jié)量越小(原因分析見3.4節(jié)),4組局部凝結(jié)量在管深0.5 m后,均呈明顯的線性增長(zhǎng)。

    3.4剪切力對(duì)管內(nèi)流場(chǎng)的影響

    3.4.1剪切力對(duì)液膜厚度的影響取文獻(xiàn)[16]表3實(shí)驗(yàn)5的工況條件,用兩組CFD數(shù)值模擬0~0.5 m管深液膜厚度,第1組忽略剪切力;第2組計(jì)算剪切力,蒸汽入口流速4.86 m·s-1。兩組數(shù)值模擬結(jié)果的對(duì)比見圖11。

    圖11 線性-對(duì)數(shù)標(biāo)尺下不同蒸汽流速的液膜平均厚度Fig.11 Average film thickness vs axial distance from inlet in line-log scale by variation of vapour velocity

    圖12 兩組CFD對(duì)管深0~1 m區(qū)域液膜厚度的數(shù)值模擬Fig.12 Numerical simulation of film thickness on plate from 0 to 0.5 m in tube axial length under global Cartesian coordinates

    由圖11(a)可見,由于翅管傾斜60°布置,重力對(duì)管壁上凝結(jié)液膜的體積力,被分解成沿管高和管深方向上的兩個(gè)分量,分別驅(qū)動(dòng)液膜沿管高和管深單向累積,使得液膜厚度在管高和管深方向上,呈明顯的梯田狀分布。

    由圖11(b)可見,0~0.2 m管段,剪切力是液膜的決定性驅(qū)動(dòng)力,剪切力迫使液膜厚度沿管高和管深趨于平均,與水泥刮刀抹平水泥的效果非常類似,在此區(qū)域,液膜厚度呈明顯的平原狀分布。

    取文獻(xiàn)[16]表3實(shí)驗(yàn)5的工況條件,用4組CFD計(jì)算全管段的液膜厚度,第1組忽略剪切力,第2~4組計(jì)算剪切力,蒸汽入口流速4.86、7.29和9.72 m·s-1,在線性-對(duì)數(shù)標(biāo)尺下,對(duì)比4組數(shù)值計(jì)算結(jié)果如圖12所示。

    由圖12可見,從蒸汽入口到管深0.2 m,剪切力作用壓倒重力作用,明顯削薄了液膜厚度。從管深0.5 m開始到0.8 m,剪切力和重力共同作用于液膜,重力加速度累積液膜,剪切力迫使液膜變薄。隨著蒸汽流動(dòng)凝結(jié),剪切力減小,剪切力與重力對(duì)液膜的作用此消彼長(zhǎng)。從管深1 m到蒸汽出口,液膜完全受重力作用控制,4組CFD計(jì)算的液膜厚度穩(wěn)定不變,表示液膜厚度不再受蒸汽入口流速(交界面剪切力)影響。

    由圖12還可見,從蒸汽入口到管深1.0 m,蒸汽流速越大,液膜越薄。由式(22)、式(23)的計(jì)算結(jié)果可知,由于大扁管軸向深度(2600 mm)遠(yuǎn)遠(yuǎn)大于管高尺(50 mm),蒸汽凝結(jié)量主要來(lái)自對(duì)液膜沿管高方向匯流量的積分,而由式(13)、式(14)可知,匯流量取決于液膜管高方向平均流速和液膜厚度,因?yàn)榧羟辛?duì)液膜管高方向流速的影響可忽略不計(jì),蒸汽凝結(jié)量變成液膜厚度的單元積分函數(shù)。這就是圖10中蒸汽流速越高,蒸汽凝結(jié)量越小的原因。

    3.4.2剪切力對(duì)交界面速度的影響取文獻(xiàn)[16]表3實(shí)驗(yàn)5的工況條件,用4組CFD數(shù)值模擬全管段交界面Z速度(沿軸向管深方向的分量),第1組忽略剪切力;第2~4組計(jì)算剪切力,蒸汽入口流速4.86、7.29和9.72 m·s-1。4組計(jì)算結(jié)果在對(duì)數(shù)-線性標(biāo)尺下的對(duì)比如圖13所示,前兩組數(shù)值模擬結(jié)果的對(duì)比見圖14。

    由圖13和圖14(a)可見,當(dāng)忽略交界面剪切力時(shí),液膜由重力驅(qū)動(dòng),重力加速度迫使液膜沿管深單向加速,同時(shí)隨著蒸汽不斷凝結(jié),液膜厚度增加導(dǎo)致液膜移動(dòng)減緩,重力加速度和蒸汽凝結(jié)量對(duì)交界面速度的影響至管深0.1 m處達(dá)到平衡,從此處一直到蒸汽出口,交界面速度保持穩(wěn)定不變。

    由圖13和圖14(b)可見,剪切力對(duì)液膜加速作用明顯,忽略剪切力計(jì)算交界面速度會(huì)導(dǎo)致其數(shù)值減小1個(gè)數(shù)量級(jí);蒸汽流速越高,交界面速度越大。剪切力對(duì)液膜的加速作用在管深0~0.2 m區(qū)域達(dá)到頂峰,然后由于蒸汽持續(xù)相變凝結(jié),蒸汽質(zhì)量喪失,能量和動(dòng)量也隨之凈減,剪切力對(duì)液膜的加速作用減弱,液膜流動(dòng)速度不斷減慢。

    圖13 對(duì)數(shù)-線性標(biāo)尺下不同蒸汽流速的交界面平均速度Fig.13 Average W interfacial velocity vs axial distance in log-line scale from inlet by variation of vapour velocity

    圖14 兩組CFD對(duì)全管深區(qū)域交界面W速度的數(shù)值模擬Fig.14 Numerical simulation of W interfacial velocity under global Cartesian coordinates

    由圖13和圖14(b)還發(fā)現(xiàn),蒸汽入口流速4.86 m·s-1時(shí),液膜流速?gòu)墓苌?.0 m處急劇下降,表示在蒸汽出口區(qū),交界面剪切力減小至0.15 N·m-2以下時(shí)(圖8),液膜流動(dòng)受為凝結(jié)完蒸汽回流的影響較大,導(dǎo)致液膜速度急劇下降。而蒸汽流速7.29和9.72 m·s-1時(shí),交界面剪切力在蒸汽出口區(qū)大于0.15 N·m-2,液膜流動(dòng)不受蒸汽回流影響,不會(huì)發(fā)生流速急劇降低的現(xiàn)象。

    4 結(jié) 論

    (1)管內(nèi)交界面剪切力的大小由蒸汽入口流速?zèng)Q定,蒸汽流速越高,剪切力越大,剪切力隨著蒸汽流動(dòng)凝結(jié)持續(xù)遞減。

    (2)剪切力提高了管內(nèi)局部凝結(jié)HTC,減小了局部凝結(jié)量,當(dāng)蒸汽流速超過凝汽器最大蒸汽流量對(duì)應(yīng)的流速時(shí),剪切力對(duì)凝結(jié)的影響達(dá)到極限。

    (3)剪切力在0~0.8 m管段內(nèi)削薄液膜厚度;從1.0 m開始,重力成為控制液膜厚度的主要驅(qū)動(dòng)力,剪切力對(duì)液膜厚度的影響可忽略不計(jì);液膜厚度變薄,是導(dǎo)致蒸汽流速越大蒸汽凝結(jié)量越小的原因。

    (4)在0~0.2 m管段,剪切力對(duì)液膜的加速作用明顯,隨著蒸汽流動(dòng),剪切力的加速作用逐漸減弱;蒸汽流速4.86 m·s-1時(shí),在蒸汽出口區(qū),液膜流速受蒸汽回流影響而急劇下降。

    符號(hào)說明

    A ——面積,m2

    cp——比定壓熱容,J·kg-1·K-1

    g ——重力加速度,m·s-2

    H ——比焓,J·kg-1

    h ——傳熱系數(shù),W·m-2·K-1

    k ——湍動(dòng)能,J·kg-1

    m ——質(zhì)量流量,kg·s-1

    n ——翅管數(shù)

    p ——壓力,Pa

    q″ ——熱通量,kW·m-2

    S,s ——控制方程源項(xiàng)

    T ——溫度,K

    ΔT ——對(duì)數(shù)平均溫差,K

    u,v,w ——速度分量,m·s-1

    V ——體積,m3

    v ——速度,m·s-1

    x,y,z ——坐標(biāo),m

    β ——翅管傾斜角,rad

    δ ——液膜厚度,m

    ε ——能量耗散率,m

    η ——蛇形翅片肋效率

    κ ——熱導(dǎo)率,W·m-1·K-1

    λ——汽化潛熱,kJ·kg-1

    μ——?jiǎng)恿︷ざ龋琸g·m-1·s-1

    ρ——密度,kg·m-3

    τ,-τ——分別為剪切力、剪切力張量,N·m-2

    Φ——熱量,J

    下角標(biāo)

    a——空氣

    con——凝結(jié)

    e——能量方程

    eff——有效值

    f——液膜

    fin——翅片

    I——交界面

    i——管內(nèi)

    m——質(zhì)量方程

    o——翅側(cè)

    s——蒸汽飽和狀態(tài)

    t——大扁管與翅片

    tube——大扁管

    v——能量方程

    w——管壁

    x,y,z——沿x,y,z坐標(biāo)方向

    1,2——分別為翅側(cè)空氣進(jìn)口、翅側(cè)空氣出口

    References

    [1] 郭民臣,紀(jì)執(zhí)琴,安廣然,等. 干-濕冷卻系統(tǒng)對(duì)空冷機(jī)組熱經(jīng)濟(jì)性影響的分析[J]. 化工學(xué)報(bào),2015,66(1): 433-440. DOI: 10.11949/j.issn.0438-1157.20140990. GUO M C,JI Z Q,AN G R,et al. Thermal economy analysis of power unit with wet-dry hybrid cooling system [J]. CIESC Journal,2015,66(1): 433-440. DOI: 10.11949/j.issn.0438-1157.20140990.

    [2] 明廷臻,黨艷輝,劉偉,等. 橢圓管矩形翅片空冷器流體流動(dòng)與傳熱特性數(shù)值分析 [J]. 化工學(xué)報(bào),2009,60(6): 1380-1384. DOI: 10.3321/j.issn:0438-1157.2009.06.007. MING T Z,DANG Y H,LIU W,et al. Numerical analysis of fluid flow and heat transfer characteristics on elliptical tube with rectangular fins of air cooler [J]. CIESC Journal,2009,60(6): 1380-1384. DOI: 10.3321/j.issn:0438-1157.2009.06.007.

    [3] NUSSELT W. Die Oberflachenkondesation des Wasserdamffes the surface condensation of water [J]. Zetrschr. Ver. Deutch. Ing.,1916,60: 541-546.

    [4] CHEN S J,REED J G,TIEN C L. Reflux condensation in a two-phase closed thermosyphon [J]. Int. J. Heat Mass Tran.,1984,27(9): 1587-1594. DOI: 10.1016/0017-9310(84)90271-0.

    [5] FIEDLER S,AURACHER H. Experimental and theoretical investigation of reflux condensation in an inclined small diameter tube [J]. Int. J. Heat Mass Tran.,2004,47(19): 4031-4043. DOI: 10.1016/j.ijheatmasstransfer.2004.06.005.

    [6] ROHSENOW W M,HARTNETT J P,CHO Y I. Handbook of Heat Transfer [M]. 3rd ed. New York,USA: McGraw-Hill,1998: 14.7-14.8.

    [7] LEE K Y,KIM M H. Effect of an interfacial shear stress on steamcondensation in the presence of a noncondensable gas in a vertical tube [J]. Int. J. Heat Mass Tran.,2008,51(21): 5333-5343. DOI: 10.1016/j.ijheatmasstransfer.2008.03.017.

    [8] 張俊霞,王立,李運(yùn)剛,等. 界面剪切力對(duì)蒸汽垂直下流膜狀凝結(jié)傳熱的影響分析 [J]. 化工學(xué)報(bào),2011,62(10): 2733-2739. ZHANG J X,WANG L,LI Y G,et al. Analysis of interfacial shear stress effects on vapor downward condensation heat transfer [J]. CIESC Journal,2011,62(10): 2733-2739.

    [9] HIRT C W,NICHOLS B D. Volume of fluid (VOF) method for the dynamics of free boundaries [J]. J. Comput. Phys.,1981,39(81): 201-225. DOI: 10.1016/0021-9991(81)90145-5.

    [10] KAYS W M,CRAWFORD M E. Convective Heat and Mass Transfer,[M]. 3rd ed. New York,USA: McGraw-Hill,1993.

    [11] MUNOZ-COBO J L,HERRANZ L,SANCHO J,et al. Turbulent vapor condensation with noncondensable gases in vertical tubes [J]. Int. J. Heat Mass Tran.,1996,39(96): 3249-3260. DOI: 10.1016/0017-9310(96)00012-9.

    [12] KIM D E,YANG K H,HWANG K W,et al. Pure steam condensation model with laminar film in a vertical tube [J]. Int. J. Multiphase Flow,2011,37(8): 941-946. DOI: 10.1016/j.ijmultiphaseflow.2011.04.006.

    [13] MUNSON B R,YOUNG D F,OKIISHI T H. Fundamentals of Fluid Mechanics [M]. 4th ed. New York,USA: John Wiley & Sons,2002.

    [14] OH S,Revankar S T. Analysis of the complete condensation in a vertical tube passive condenser [J]. Int. Commun. Heat Mass,2005,32(2): 716-727. DOI: 10.1016/j.icheatmasstransfer.2004.10.013.

    [15] ZHOU W,HENDERSON G,REVANKAR S T. Condensation in a vertical tube bundle passive condenser (Ⅰ): Through flow condensation [J]. Int. J. Heat Mass Tran.,2010,53(5): 1146-1155. DOI: 10.1016/j.ijheatmasstransfer.2009.10.039.

    [16] 杜小澤,楊立軍,金衍勝,等.火電站直接空冷凝汽器傳熱系數(shù)實(shí)驗(yàn)關(guān)聯(lián)式 [J].中國(guó)電機(jī)工程學(xué)報(bào),2008,28(14): 32-37. DU X Z,YANG L J,JIN Y S,et al.Development of experimental correlation for heat transfer coefficient of direct air-cooled condenser in power plant [J]. Proceedings of the CSEE,2008,28(14): 32-37.

    [17] BATCHELOR G K.An Introduction to Fluid Dynamics[M]. Cambridge,UK: Cambridge University Press,1967.

    [18] MCADAMS W H. Heat Transmission[M]. New York,USA: McGraw-Hill,1954.

    [19] WEISSTEIN. Euler angles,Wolfram MathWorld[EB/OL]. Champaign,USA: Wolfram Research,Inc.,[2015-05-29]. http://mathworld. wolfram.com/EulerAngle.html

    [20] HARUN B,STEPHEN D,Edward K L. Numerical modeling of finned heat exchangers [J]. Appl. Therm. Eng.,2013,61(2): 278-288. DOI: 10.1016/j.applthermaleng.2013.08.002.

    CFD simulation on effect of interfacial shear force on water vapor condensation in inclined flat tube

    DENG Hui1,BAI Yan1,LI Xinxin2,ZHANG Dongming3
    (1School of Control and Computer Engineering,North China Electric Power University,Beijing 102206,China;2China Power Complete Equipment Co.,Ltd.,Beijing 100080,China;3Beijing GDZS Control Technology Co.,Beijing 102200,China)

    Abstract:Considering the interfacial shear force,a mathematical model to the condensation of turbulent vapor flowing downward in an inclined flat tube is proposed and implemented in computational fluid dynamics (CFD). The predicted results from the CFD model are compared with the experimental results from the literature for the vapor condensation in a prototype tube (2600 mm length,3 mm width and 50 mm altitude in 60° inclination to vertical). It is found that the condensate rate and mean condensation heat transfer coefficient (HTC) from CFD simulation agree very well with the experimental quantities. Using CFD model to calculate the interfacial shear stress by varying vapour velocity,the results demonstrate that the value of shear force depends on the vapor velocity at the tube inlet,and shear force decreases continuously with the vapor flow and condensation. Simulating the interfacial shear effect on the condensation,it shows that the interfacial shear increases the local condensation HTC,and meanwhile,reduces the local condensate rate. The simulation results also shows that the interfacial shear weakens the gravitational effect on the film accumulation and obviously decreases film thickness from 0 to 0.8 m in the tube axial length. However,from 1.0 m to the tube outlet,the gravitational force dominates over the shear force,and thus the shear effects can be completely neglected. It is also found that the condensate film isspeeded up particularly from 0 to 0.2 m in axial length thanking for the interfacial shear.

    Key words:turbulent flow;condensation;interfacial shear stress;inclined flat tube;numerical simulation

    DOI:10.11949/j.issn.0438-1157. 20150976

    中圖分類號(hào):TK 124

    文獻(xiàn)標(biāo)志碼:A

    文章編號(hào):0438—1157(2016)04—1215—10

    Corresponding author:Prof. BAI Yan,by@ncupu.edu.cn

    猜你喜歡
    數(shù)值模擬湍流
    “湍流結(jié)構(gòu)研究”專欄簡(jiǎn)介
    重氣瞬時(shí)泄漏擴(kuò)散的湍流模型驗(yàn)證
    張家灣煤礦巷道無(wú)支護(hù)條件下位移的數(shù)值模擬
    科技視界(2016年18期)2016-11-03 23:14:27
    張家灣煤礦開切眼錨桿支護(hù)參數(shù)確定的數(shù)值模擬
    科技視界(2016年18期)2016-11-03 22:57:21
    跨音速飛行中機(jī)翼水汽凝結(jié)的數(shù)值模擬研究
    科技視界(2016年18期)2016-11-03 20:38:17
    姚橋煤礦采空區(qū)CO2防滅火的數(shù)值模擬分析
    雙螺桿膨脹機(jī)的流場(chǎng)數(shù)值模擬研究
    科技視界(2016年22期)2016-10-18 14:53:19
    一種基于液壓緩沖的減震管卡設(shè)計(jì)與性能分析
    科技視界(2016年20期)2016-09-29 11:08:27
    湍流十章
    “青春期”湍流中的智慧引渡(三)
    欧美激情久久久久久爽电影 | 国精品久久久久久国模美| 亚洲精品乱久久久久久| 操出白浆在线播放| 久久性视频一级片| 欧美激情 高清一区二区三区| 天天躁日日躁夜夜躁夜夜| 久久精品国产亚洲av高清一级| 在线观看免费视频网站a站| 精品一区二区三区视频在线观看免费 | 一级a爱视频在线免费观看| 久久影院123| 国产欧美日韩一区二区精品| 国产精品一区二区在线观看99| 天天影视国产精品| 成人特级黄色片久久久久久久| 亚洲一区二区三区不卡视频| 欧美丝袜亚洲另类 | 亚洲自偷自拍图片 自拍| 亚洲一区高清亚洲精品| 女警被强在线播放| 欧美日韩福利视频一区二区| 亚洲一卡2卡3卡4卡5卡精品中文| 国产成人影院久久av| 亚洲七黄色美女视频| avwww免费| 久热爱精品视频在线9| 精品第一国产精品| 国产精品欧美亚洲77777| 国产成人欧美在线观看 | 欧美最黄视频在线播放免费 | 日本a在线网址| 精品免费久久久久久久清纯 | 人妻丰满熟妇av一区二区三区 | 不卡av一区二区三区| 韩国av一区二区三区四区| 亚洲av片天天在线观看| 狠狠婷婷综合久久久久久88av| 日韩免费av在线播放| 亚洲情色 制服丝袜| 两性夫妻黄色片| 91精品三级在线观看| 黄网站色视频无遮挡免费观看| 国产高清videossex| 很黄的视频免费| 麻豆av在线久日| 18在线观看网站| 一级,二级,三级黄色视频| 国产精品久久久久久人妻精品电影| 国产av精品麻豆| 叶爱在线成人免费视频播放| 欧美日韩中文字幕国产精品一区二区三区 | 大片电影免费在线观看免费| 男女之事视频高清在线观看| a级片在线免费高清观看视频| 丝瓜视频免费看黄片| 国产精品一区二区精品视频观看| 久久 成人 亚洲| www.精华液| 日本vs欧美在线观看视频| 18禁裸乳无遮挡免费网站照片 | 成年人午夜在线观看视频| 欧美大码av| 一进一出抽搐动态| 人人妻人人澡人人爽人人夜夜| 久久狼人影院| 日韩欧美一区二区三区在线观看 | 99在线人妻在线中文字幕 | 在线观看日韩欧美| 波多野结衣一区麻豆| 男女免费视频国产| 日本精品一区二区三区蜜桃| 91字幕亚洲| 少妇粗大呻吟视频| 国产在视频线精品| 搡老乐熟女国产| 丁香六月欧美| 丁香六月欧美| 午夜免费鲁丝| 久久国产亚洲av麻豆专区| 亚洲午夜理论影院| 国产精品乱码一区二三区的特点 | 亚洲中文av在线| 日韩熟女老妇一区二区性免费视频| 国产精品九九99| 老司机影院毛片| 精品国产一区二区三区久久久樱花| 嫩草影视91久久| 国产精品98久久久久久宅男小说| 最近最新中文字幕大全电影3 | 国产高清激情床上av| 黄网站色视频无遮挡免费观看| 757午夜福利合集在线观看| 国产人伦9x9x在线观看| 少妇裸体淫交视频免费看高清 | 亚洲av日韩在线播放| 精品国产乱码久久久久久男人| 精品免费久久久久久久清纯 | 国产在线精品亚洲第一网站| 久久午夜亚洲精品久久| 大香蕉久久网| 巨乳人妻的诱惑在线观看| 国产亚洲一区二区精品| 欧美精品啪啪一区二区三区| 一级毛片高清免费大全| 最新美女视频免费是黄的| 久久中文字幕人妻熟女| 精品国产一区二区三区四区第35| 三上悠亚av全集在线观看| 国产精品国产高清国产av | 久久狼人影院| 亚洲国产毛片av蜜桃av| 午夜免费观看网址| 国产欧美亚洲国产| 亚洲人成77777在线视频| 日本a在线网址| 成人手机av| 国产成人av教育| 久久精品成人免费网站| 国产免费现黄频在线看| 亚洲精品一二三| 搡老岳熟女国产| 18禁国产床啪视频网站| www.999成人在线观看| 国产高清激情床上av| 一a级毛片在线观看| 日韩欧美一区二区三区在线观看 | 丝袜美腿诱惑在线| av网站免费在线观看视频| 成人手机av| bbb黄色大片| 久久久久国产一级毛片高清牌| 校园春色视频在线观看| 国产精品永久免费网站| 亚洲aⅴ乱码一区二区在线播放 | 国产亚洲欧美98| 制服诱惑二区| 欧美激情高清一区二区三区| 国产精品久久久人人做人人爽| 一进一出好大好爽视频| 久久国产精品男人的天堂亚洲| 国产精品免费一区二区三区在线 | 中亚洲国语对白在线视频| 日本一区二区免费在线视频| 精品无人区乱码1区二区| 美女扒开内裤让男人捅视频| 久久久久久久久免费视频了| 亚洲欧美一区二区三区久久| 99精品久久久久人妻精品| 另类亚洲欧美激情| 久久精品国产a三级三级三级| 国产高清视频在线播放一区| 亚洲欧美激情在线| 黄网站色视频无遮挡免费观看| 日韩欧美国产一区二区入口| 黄色成人免费大全| 欧美日韩瑟瑟在线播放| 欧美国产精品va在线观看不卡| 91九色精品人成在线观看| 男男h啪啪无遮挡| 亚洲av片天天在线观看| 91在线观看av| 亚洲欧美一区二区三区久久| 国产亚洲精品一区二区www | 1024香蕉在线观看| 亚洲成av片中文字幕在线观看| 岛国毛片在线播放| 国产三级黄色录像| 捣出白浆h1v1| 国产精品av久久久久免费| 久久久国产精品麻豆| 亚洲av美国av| 黄色视频,在线免费观看| 精品国产一区二区三区久久久樱花| 一区二区三区精品91| 国产精品亚洲av一区麻豆| 制服人妻中文乱码| 国产日韩一区二区三区精品不卡| 男女高潮啪啪啪动态图| av国产精品久久久久影院| 精品国产国语对白av| 女同久久另类99精品国产91| 大型黄色视频在线免费观看| 99久久精品国产亚洲精品| 午夜精品在线福利| 欧美丝袜亚洲另类 | 国产伦人伦偷精品视频| 午夜福利乱码中文字幕| 在线视频色国产色| 黑人巨大精品欧美一区二区蜜桃| 少妇 在线观看| 欧美人与性动交α欧美软件| 精品卡一卡二卡四卡免费| 精品第一国产精品| 高清av免费在线| 国产人伦9x9x在线观看| 亚洲情色 制服丝袜| 黄色毛片三级朝国网站| 激情在线观看视频在线高清 | 久久亚洲真实| 国产精品亚洲av一区麻豆| 免费观看精品视频网站| 97人妻天天添夜夜摸| 亚洲精品在线观看二区| 十八禁网站免费在线| 99国产综合亚洲精品| 国产成人精品久久二区二区91| 午夜福利在线免费观看网站| 亚洲性夜色夜夜综合| 99精国产麻豆久久婷婷| 香蕉丝袜av| 欧美在线一区亚洲| 国产黄色免费在线视频| 一级作爱视频免费观看| 美女 人体艺术 gogo| 女人被狂操c到高潮| 丰满人妻熟妇乱又伦精品不卡| 精品一区二区三区av网在线观看| 色综合婷婷激情| 黑人巨大精品欧美一区二区mp4| 国产高清视频在线播放一区| 女人高潮潮喷娇喘18禁视频| av视频免费观看在线观看| 欧美黑人欧美精品刺激| 国产一区二区三区综合在线观看| 啦啦啦 在线观看视频| av有码第一页| 啪啪无遮挡十八禁网站| 十分钟在线观看高清视频www| 一级毛片精品| 身体一侧抽搐| 热99久久久久精品小说推荐| 精品免费久久久久久久清纯 | 久久国产精品影院| 日韩成人在线观看一区二区三区| 国产主播在线观看一区二区| 国产亚洲欧美在线一区二区| 国产成人精品在线电影| 老鸭窝网址在线观看| 精品熟女少妇八av免费久了| 高清在线国产一区| 正在播放国产对白刺激| 成人影院久久| 91九色精品人成在线观看| 999精品在线视频| 精品国产美女av久久久久小说| 满18在线观看网站| 欧美色视频一区免费| 91大片在线观看| 777久久人妻少妇嫩草av网站| 亚洲伊人色综图| 男女床上黄色一级片免费看| 在线观看午夜福利视频| 老司机福利观看| 国产在视频线精品| av天堂在线播放| 十八禁网站免费在线| 欧美黑人欧美精品刺激| 国产免费现黄频在线看| 成人亚洲精品一区在线观看| 欧美日韩中文字幕国产精品一区二区三区 | 侵犯人妻中文字幕一二三四区| 妹子高潮喷水视频| 国产亚洲精品第一综合不卡| 国产一卡二卡三卡精品| 叶爱在线成人免费视频播放| 欧美日韩福利视频一区二区| 精品久久久久久久久久免费视频 | 夫妻午夜视频| 欧美激情久久久久久爽电影 | 欧美丝袜亚洲另类 | 国产蜜桃级精品一区二区三区 | 免费不卡黄色视频| 国产精品 欧美亚洲| 亚洲国产欧美网| 麻豆乱淫一区二区| 91成人精品电影| 亚洲精品美女久久久久99蜜臀| 欧美日韩乱码在线| 别揉我奶头~嗯~啊~动态视频| 老司机福利观看| 在线观看免费高清a一片| 女同久久另类99精品国产91| 国产亚洲精品第一综合不卡| www.999成人在线观看| 亚洲精华国产精华精| 免费不卡黄色视频| 久久久久久久久免费视频了| 在线十欧美十亚洲十日本专区| 亚洲伊人色综图| 亚洲熟妇熟女久久| 天天操日日干夜夜撸| 很黄的视频免费| 久久国产精品人妻蜜桃| av免费在线观看网站| 亚洲国产精品sss在线观看 | 亚洲国产欧美日韩在线播放| 亚洲熟妇中文字幕五十中出 | 青草久久国产| 亚洲第一欧美日韩一区二区三区| 欧美一级毛片孕妇| 黄色丝袜av网址大全| 日日夜夜操网爽| 黄色成人免费大全| 1024香蕉在线观看| 久久精品91无色码中文字幕| 久久久久视频综合| 国产成人av激情在线播放| 在线观看www视频免费| 亚洲成人国产一区在线观看| 亚洲aⅴ乱码一区二区在线播放 | 国产淫语在线视频| 免费人成视频x8x8入口观看| 老司机福利观看| 亚洲一区中文字幕在线| 亚洲片人在线观看| av超薄肉色丝袜交足视频| 成人18禁在线播放| 久久热在线av| 国产精品99久久99久久久不卡| 999久久久精品免费观看国产| 嫩草影视91久久| 日韩精品免费视频一区二区三区| 人妻一区二区av| 欧美日韩乱码在线| 欧美日韩成人在线一区二区| 日本一区二区免费在线视频| 91大片在线观看| 不卡av一区二区三区| av天堂在线播放| 狠狠狠狠99中文字幕| 91麻豆av在线| 亚洲精品在线美女| 国产欧美亚洲国产| av片东京热男人的天堂| 热99国产精品久久久久久7| 国产又色又爽无遮挡免费看| 夜夜爽天天搞| 国产人伦9x9x在线观看| 精品国产超薄肉色丝袜足j| 久久久久精品人妻al黑| 男人的好看免费观看在线视频 | 在线观看免费日韩欧美大片| 操出白浆在线播放| 免费在线观看日本一区| 亚洲成人免费电影在线观看| 80岁老熟妇乱子伦牲交| 精品福利永久在线观看| 男女高潮啪啪啪动态图| 成人特级黄色片久久久久久久| 宅男免费午夜| 男女免费视频国产| svipshipincom国产片| 欧美日韩精品网址| 久久亚洲精品不卡| 免费看十八禁软件| 多毛熟女@视频| 黄色 视频免费看| 国产男女超爽视频在线观看| 可以免费在线观看a视频的电影网站| 色婷婷久久久亚洲欧美| 人妻 亚洲 视频| 如日韩欧美国产精品一区二区三区| 亚洲精品久久午夜乱码| 日本五十路高清| 久久人妻福利社区极品人妻图片| 无人区码免费观看不卡| 少妇猛男粗大的猛烈进出视频| 久久久国产欧美日韩av| 夜夜爽天天搞| 99精品在免费线老司机午夜| 亚洲美女黄片视频| 精品福利永久在线观看| 色综合欧美亚洲国产小说| 国产精品98久久久久久宅男小说| 亚洲精品美女久久久久99蜜臀| 欧美日韩福利视频一区二区| 韩国av一区二区三区四区| 国产在线观看jvid| 亚洲精品中文字幕一二三四区| 久久久国产成人精品二区 | 9色porny在线观看| 黑人巨大精品欧美一区二区蜜桃| 18禁美女被吸乳视频| 正在播放国产对白刺激| 国产亚洲av高清不卡| 亚洲精品一二三| a在线观看视频网站| 久久久国产一区二区| www日本在线高清视频| 美女午夜性视频免费| 18禁美女被吸乳视频| 成熟少妇高潮喷水视频| 国产精品免费视频内射| 99在线人妻在线中文字幕 | 国产一区二区三区视频了| 国产精品自产拍在线观看55亚洲 | 亚洲成人免费av在线播放| 中亚洲国语对白在线视频| 成人影院久久| 12—13女人毛片做爰片一| 国产欧美日韩精品亚洲av| 亚洲人成电影免费在线| 看片在线看免费视频| 在线观看一区二区三区激情| 国产精华一区二区三区| 91在线观看av| 国产高清激情床上av| 99riav亚洲国产免费| 国产激情欧美一区二区| 精品久久久久久,| 99riav亚洲国产免费| 中文字幕高清在线视频| 80岁老熟妇乱子伦牲交| 日本一区二区免费在线视频| 午夜激情av网站| 99在线人妻在线中文字幕 | 飞空精品影院首页| 久久人妻熟女aⅴ| 欧美大码av| 亚洲少妇的诱惑av| 18在线观看网站| 久久精品国产综合久久久| 亚洲视频免费观看视频| 午夜福利在线观看吧| 法律面前人人平等表现在哪些方面| 又黄又爽又免费观看的视频| 日本欧美视频一区| 男女下面插进去视频免费观看| 欧美一级毛片孕妇| 欧美不卡视频在线免费观看 | 别揉我奶头~嗯~啊~动态视频| 又黄又粗又硬又大视频| 亚洲片人在线观看| 嫁个100分男人电影在线观看| 丁香六月欧美| av欧美777| 男人操女人黄网站| 久久精品亚洲熟妇少妇任你| 亚洲精品中文字幕一二三四区| 国产av又大| 一区二区三区国产精品乱码| 亚洲第一av免费看| 国产成人影院久久av| 日日爽夜夜爽网站| 亚洲美女黄片视频| 一边摸一边抽搐一进一小说 | 日韩 欧美 亚洲 中文字幕| av超薄肉色丝袜交足视频| 精品国内亚洲2022精品成人 | 国产精品一区二区在线观看99| 国产成人啪精品午夜网站| 午夜视频精品福利| 亚洲一码二码三码区别大吗| 成人影院久久| 无遮挡黄片免费观看| 日韩人妻精品一区2区三区| 高潮久久久久久久久久久不卡| 女人高潮潮喷娇喘18禁视频| 久久久水蜜桃国产精品网| 亚洲人成电影免费在线| 免费一级毛片在线播放高清视频 | 少妇裸体淫交视频免费看高清 | 男女午夜视频在线观看| 久久九九热精品免费| 一二三四在线观看免费中文在| 中文字幕人妻丝袜制服| 91精品国产国语对白视频| 精品一品国产午夜福利视频| 精品国产一区二区久久| 老司机在亚洲福利影院| 日韩熟女老妇一区二区性免费视频| 精品午夜福利视频在线观看一区| 建设人人有责人人尽责人人享有的| 亚洲中文av在线| 黄色丝袜av网址大全| 国产一区有黄有色的免费视频| 精品熟女少妇八av免费久了| 在线观看66精品国产| 精品国产乱码久久久久久男人| 久久国产精品人妻蜜桃| 波多野结衣av一区二区av| 国产精品偷伦视频观看了| 欧美国产精品va在线观看不卡| 亚洲中文av在线| 欧洲精品卡2卡3卡4卡5卡区| 嫩草影视91久久| 在线观看免费日韩欧美大片| 岛国在线观看网站| 国产精品.久久久| 伦理电影免费视频| 日韩精品免费视频一区二区三区| 热99re8久久精品国产| 一区二区三区激情视频| 亚洲av欧美aⅴ国产| 一本大道久久a久久精品| 亚洲精品粉嫩美女一区| 欧美大码av| 国产蜜桃级精品一区二区三区 | 热re99久久国产66热| 精品视频人人做人人爽| 如日韩欧美国产精品一区二区三区| 一本综合久久免费| 国产欧美日韩一区二区三| 视频区图区小说| 又黄又爽又免费观看的视频| 欧美黄色片欧美黄色片| 免费不卡黄色视频| 久久国产精品影院| 免费女性裸体啪啪无遮挡网站| 成人影院久久| 亚洲片人在线观看| 亚洲成a人片在线一区二区| 婷婷成人精品国产| 老鸭窝网址在线观看| 免费在线观看视频国产中文字幕亚洲| 视频区欧美日本亚洲| 亚洲欧美日韩高清在线视频| 国产av一区二区精品久久| 婷婷丁香在线五月| 精品人妻1区二区| 日本一区二区免费在线视频| 天堂俺去俺来也www色官网| 夜夜爽天天搞| 欧美日韩瑟瑟在线播放| 俄罗斯特黄特色一大片| 一区二区三区国产精品乱码| 日本黄色日本黄色录像| 正在播放国产对白刺激| 在线视频色国产色| 国产熟女午夜一区二区三区| 亚洲三区欧美一区| 天天操日日干夜夜撸| 一本大道久久a久久精品| 亚洲av日韩精品久久久久久密| 欧美一级毛片孕妇| 亚洲精品在线观看二区| 久久精品91无色码中文字幕| 高清欧美精品videossex| 成人免费观看视频高清| 精品久久久久久久久久免费视频 | 中文字幕高清在线视频| 亚洲成av片中文字幕在线观看| 多毛熟女@视频| 国产一区在线观看成人免费| 性少妇av在线| 大型黄色视频在线免费观看| 欧美成狂野欧美在线观看| 黑人巨大精品欧美一区二区mp4| 亚洲一码二码三码区别大吗| 91九色精品人成在线观看| 岛国毛片在线播放| 丰满迷人的少妇在线观看| 亚洲片人在线观看| av网站在线播放免费| 在线视频色国产色| 99热网站在线观看| av中文乱码字幕在线| 欧美日韩一级在线毛片| 日本黄色视频三级网站网址 | 色播在线永久视频| 亚洲中文av在线| 久久 成人 亚洲| 黄色视频不卡| 91精品国产国语对白视频| 岛国毛片在线播放| 啦啦啦在线免费观看视频4| 亚洲精品国产精品久久久不卡| 久久 成人 亚洲| 女人久久www免费人成看片| 在线国产一区二区在线| 免费高清在线观看日韩| 精品国产国语对白av| 黑丝袜美女国产一区| 亚洲午夜理论影院| 俄罗斯特黄特色一大片| 99久久99久久久精品蜜桃| 19禁男女啪啪无遮挡网站| 淫妇啪啪啪对白视频| 18禁国产床啪视频网站| 国产高清videossex| 久久精品aⅴ一区二区三区四区| 人妻久久中文字幕网| 久久国产精品男人的天堂亚洲| 自拍欧美九色日韩亚洲蝌蚪91| 巨乳人妻的诱惑在线观看| 欧美一级毛片孕妇| 欧美av亚洲av综合av国产av| 女人精品久久久久毛片| 久久久国产欧美日韩av| 久久草成人影院| 国产淫语在线视频| www.精华液| 久久青草综合色| 下体分泌物呈黄色| 欧美在线一区亚洲| 成年人午夜在线观看视频| 美女高潮到喷水免费观看| 成人精品一区二区免费| 99re6热这里在线精品视频| 欧美大码av| 亚洲精品粉嫩美女一区| netflix在线观看网站| 国产一区二区三区在线臀色熟女 | 成熟少妇高潮喷水视频| 国产精品香港三级国产av潘金莲| 91麻豆av在线| 日韩 欧美 亚洲 中文字幕| 国产激情久久老熟女| 久久天堂一区二区三区四区| 淫妇啪啪啪对白视频| 亚洲精品粉嫩美女一区| 丝袜美足系列| 性色av乱码一区二区三区2| 欧美 亚洲 国产 日韩一|