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

    傾斜多孔介質方腔內納米流體自然對流的格子Boltzmann 方法模擬*

    2020-08-29 07:32:58張貝豪鄭林
    物理學報 2020年16期
    關鍵詞:熱效率對流壁面

    張貝豪 鄭林

    (南京理工大學能源與動力工程學院, 南京 210094)

    1 引 言

    近幾十年來人們在電力電子設備、太陽能集熱器、核反應堆等領域的強化傳熱方面開展了廣泛的研究, 而自然對流是該研究領域所涉及到的經典流動與傳熱問題之一. 在這些設備當中, 大多采用水、冷凍液和乙二醇等傳統(tǒng)換熱工質來進行強化換熱, 與之相比, 在傳統(tǒng)的純液體中添加納米顆粒,所得到的納米流體具有較高的導熱系數(shù), 并且能滿足一些特殊條件下的傳熱與冷卻要求, 因而其強化傳熱特性引起了研究者的廣泛關注[1].

    國內外許多學者已經對方腔內納米流體自然對流現(xiàn)象進行了研究. 例如: 李新芳和朱冬生[2]采用CFD 方法對影響納米流體強化傳熱的因素進行了分析, 結果表明, 隨著納米顆粒體積分數(shù)的增大,流體的能量傳輸?shù)玫綇娀? 張晶等[3]采用4 種不同的粘度模型計算納米流體粘度, 并對二維U 型腔內Al2O3-H2O 納米流體的自然對流換熱進行了數(shù)值模擬. Khanafer 等[4]數(shù)值模擬研究了二維方形腔內的納米流體自然對流換熱, 結果表明, 隨著納米顆粒體積分數(shù)的增加, 熱效率增強. Hatami[5]采用有限元方法研究了不同納米流體(TiO2和Al2O3)在加熱肋片的矩形腔內的自然對流問題,結果表明, 隨著體積分數(shù)的增大, 雖然Al2O3-H2O 納米流體在熱肋片處Nuave數(shù)緩慢增大, 但是TiO2-H2O 納米流體在熱肋片處Nuave數(shù)先增大后減小. Chen 和Du[6], Selimefendigil 和Oztop[7],Hatami[8]研究了納米流體在方腔內的自然對流換熱, 結果表明, 隨著Ra 數(shù)的增加, 傳熱速率得到提高. Jahanshahi 等[9]采用瞬態(tài)熱線法和哈密頓模型得到了不同體積分數(shù)下的流體導熱系數(shù), 分別研究了SiO2-H2O 納米流體在熱壁面處不同Ra 數(shù)下Nuave數(shù)隨體積分數(shù)增加的變化趨勢, 結果表明采用哈密頓模型計算得到的Nuave數(shù)呈現(xiàn)輕微減小的趨勢.

    多孔介質中的自然對流與換熱問題普遍存在于自然界和工業(yè)應用中, 得到了學者們的廣泛研究. 目前, 對于封閉多孔介質腔體內自然對流的研究大致可以分為兩大類: 一類為通過水平多孔介質層的傳熱問題, 特別是水平多孔層存在高低溫壁面, 即經典的Benard 流動問題[10], 另一類是具有兩個不同溫度的垂直壁面以及兩個絕熱的水平壁面的多孔介質傳熱問題[11]. Lan 和Prakash 等[12]采用Darcy-Brinkman-Forchheimer (DBF) 模型探究了不同達西數(shù)、Ra數(shù)和畢渥數(shù)下達西和非達西流體的流動機制. Kaluri 和Basak[13]研究了二維多孔介質方腔內納米流體的自然對流換熱問題.Alsabery 等[14]研究了傾斜方腔多孔層內納米流體的自然對流問題, 結果表明施加較小的傾斜角可以顯著提高壁面的Nuave數(shù). Toosi 和Siavashi[15]對二維多孔介質腔內納米流體自然對流現(xiàn)象進行了數(shù)值模擬, 結果表明同時采用多孔介質和納米流體可以獲得最佳的傳熱性能. 相比于均勻受熱的邊界條件, 近年來, 學者們對于壁面處于非均勻溫度加熱的研究不斷增多. Wu 等[16]采用非達西模型的多孔介質控制方程, 分析了左右壁面均為正弦溫度分布的自然對流傳熱特性. Basak 等[17]對比了多孔介質方腔內均勻邊界條件和非均勻邊界條件影響下自然對流強弱的變化情況, 研究指出, 底部中心加熱壁面處為非均勻溫度分布時, 換熱效率更高. 前面的研究主要分析純流體的自然對流問題,由于納米流體具有強化傳熱特征, 考慮納米顆粒影響的多孔介質方腔內自然對流的研究同樣受到一些學者的關注. 比如: Ghasemi 和Siavashi[18]采用LBM 方法研究了不同線性溫度分布的邊界條件對多孔介質內納米流體換熱效率的影響. 此外,Sivasankaran[19]利用交替方向隱式的有限差分方法分析了下壁面為非均勻溫度分布的傾斜多孔介質方腔內Cu-H2O 納米流體的自然對流問題, 發(fā)現(xiàn)在Ra= 103時Nuave數(shù)隨納米顆粒體積分數(shù)的增加而減少.

    綜上所述, 以往研究表明納米流體在多孔介質方腔內具有強化作用, 但是對于不同邊界條件和計算模型, 其流動和換熱過程表現(xiàn)出不同結論, 有的結果表現(xiàn)為強化換熱效果[20], 有的則表現(xiàn)為削弱換熱作用[21]. 值得注意的是, 非均勻溫度邊界條件對傾斜多孔介質方腔內的納米流體流動與傳熱過程產生了重要影響. 與此同時, 由于Al2O3納米顆粒在分散系中化學性質較為穩(wěn)定且價格低廉, 故被廣泛應用于工業(yè)強化傳熱領域中, 因此, 本文對非均勻溫度邊界條件下傾斜多孔介質方腔內Al2O3-H2O 納米流體自然對流進行數(shù)值模擬, 分析了?,Ra數(shù),?和g等參數(shù)對納米流體流動傳熱特性的影響.

    2 問題描述

    2.1 物理模型

    本文主要研究多孔介質方腔內納米流體的自然對流問題. 如圖1 所示, 方腔內納米流體受到重力場的作用, 其下壁面與水平方向的夾角為g, 左壁面溫度分布為T=sin(πy/H) , 而右壁面則保持低溫Tc, 上下壁面絕熱且不可滲透, 流體在固體壁面無滑移. 由于多孔介質內復雜的物理結構, 本文將通過以下4 種假設來簡化物理問題: 1) 假設多孔介質為均質且各向同性; 2) 粘性耗散在能量方程中可以忽略不計; 3) 該流體被認為是不可壓縮的牛頓流體; 4) 可以采用標準的 Boussinesq 假設.

    圖1 物理模型示意圖Fig. 1. Schematic diagram of the physical model.

    2.2 無量綱控制方程

    基于上述基本假設, 本文采用DBF 模型, 對多孔介質方腔內納米流體的自然對流問題模型的無因次形式的連續(xù)性方程、動量方程和能量方程進行描述[22]:

    其中,X和Y是無因次橫向縱向坐標,U和V是沿X和Y方向的無因次速度,Ra是瑞利數(shù),Pr是普朗特數(shù),Da是達西數(shù),T是無因次流體溫度,P是無因次壓力,r是流體密度,μ是流體粘度,k是流體導熱系數(shù),f是納米顆粒體積分數(shù),g是方腔旋轉角度,?是多孔介質孔隙率. 下標f代表基液, 下標nf代表納米流體, 下標m代表考慮多孔介質的有效物理量.

    相應地, 無因次邊界條件由以下表達式來確定:

    上壁面:Y= 1

    下壁面:Y= 0

    左壁面:X= 0

    右壁面:X= 1

    2.3 納米流體的物性參數(shù)

    本文所研究的納米流體是在基液水中加入一定體積分數(shù)的Al2O3納米顆粒而形成的懸浮液. 假設納米流體中的水和Al2O3處于熱平衡狀態(tài), 并且納米顆粒為球形結構且均勻分布在基液中, 不存在納米顆粒的團聚和沉淀現(xiàn)象. 表1 所示為基液、Al2O3納米顆粒和多孔介質的物性參數(shù).

    利用這些已知的熱物性參數(shù), 可以得到不同體積分數(shù)以及不同孔隙率下系統(tǒng)內的熱物性參數(shù)[25,26],具體計算公式如表2 所示:

    為描述熱壁面處自然對流強度, 本研究組分析了壁面處局部Nu數(shù)和Nuave數(shù), 其計算公式分別為[26]

    表1 H2O, Al2O3 和玻璃纖維的熱物理性質Table 1. Thermophysical properties of water,Al2O3 and glass fibers.

    表2 納米流體的熱物性參數(shù)計算公式Table 2. Calculation formula for thermodynamic properties of nanofluids.

    3 網格選擇和程序驗證

    與傳統(tǒng)的數(shù)值方法相比, 格子玻爾茲曼方法(lattice Boltzmann method, LBM)的介觀性質和獨特的計算特點使其具有演化過程清晰、計算效率高、易實現(xiàn)等優(yōu)點[27], 已經被成功地應用到傳熱傳質、湍流、多相流、多孔介質等復雜流體系統(tǒng)上[28?32].因此, 本文采用LBM 開展傾斜多孔介質方腔內納米流體自然對流的數(shù)值模擬研究. 為了驗證程序的可靠性, 首先對經典的多孔介質方腔自然對流問題進行網格無關性驗證, 其中多孔介質方腔的上下壁面為絕熱壁面, 左壁面為高溫壁面, 右壁面為低溫壁面, 不考慮傾斜角變化, 流體在固體壁面處無滑移. 在該模擬研究中, 保持Da= 107,?= 0.9999,Ra= 106,Pr= 6.2 不變, 并采用四種均勻網格(80 × 80, 100 × 100, 120 × 120, 140 × 140)進行網格無關性驗證, 得到了不同網格數(shù)下熱壁面處的Nuave數(shù). 與文獻[33]對比的結果如表3 所示,隨著網格數(shù)的增大,Nuave數(shù)的相對誤差逐漸減小,當網格數(shù)為120 × 120 時,Nuave數(shù)的相對誤差為0.83%, 綜合考慮結果準確性以及計算效率, 本文所有算例將選用120 × 120 作為計算網格數(shù).

    為了進一步驗證所采用LBM 的可行性, 對二維方腔內的自然對流問題進行數(shù)值模擬, 其中, 方腔上下壁面為絕熱壁面, 左壁面為高溫壁面, 右壁面為低溫壁面, 流體在固體壁面處無滑移. 首先驗證忽略多孔介質的影響, 即?= 0.9999, 并在模擬時取Pr= 0.72 和Da= 107, 計算熱壁面處的Nuave數(shù), 并與已有文獻[33]的結果進行對比. 從表4 可知, 本文的計算結果與文獻[33]的計算結果吻合得非常好.

    表3 不同網格數(shù)與文獻[33]的Nuave 數(shù)比較Table 3. Comparison of Nuave number with literature[33] in different grids number.

    表4 本文與文獻[33]的Nuave 數(shù)值結果的比較Table 4. Comparison of Nuave number with previous literature[33].

    其次, 考慮多孔介質的影響, 采用DBF 模型進行數(shù)值模擬計算. 在模擬過程中保持?= 0.6,Pr= 1.0 和Da= 10–2不變, 同樣, 在表5 中對比了不同Ra數(shù)情況下熱壁面處的Nuave數(shù). 結果表明, 本文的結果與文獻[28]的吻合得比較好, 誤差大約在3%以內.

    表5 本文與文獻[34]的Nuave 數(shù)值結果的對比Table 5. Comparison of Nuave number with previous literature[34].

    4 結果與討論

    4.1 孔隙率 ? 對方腔內流動換熱的影響

    本文數(shù)值模擬了多孔介質方腔內Al2O3-H2O 納米流體的自然對流問題, 分析了不同物理參數(shù)對納米流體換熱效率的影響. 在模擬過程中,無量綱參數(shù)的選取為Da= 0.01,f= 0.01,Ra=105,Pr= 6.2,?= 0.9 和g= 0°.

    由于孔隙率與多孔介質滲透率有關, 它的大小會影響多孔介質方腔內納米流體的流動狀態(tài)和傳熱效率. 如圖2 所示, 當?= 0.3 時, 流場中心處形成了一個不規(guī)則的橢圓形流胞, 方腔中心區(qū)域內等溫線與水平壁面存在一定傾斜, 并且左右壁面溫度邊界層較厚, 因此其換熱方式以導熱為主; 隨著孔隙率的增大, 當?= 0.5 時, 橢圓形主流胞變得更加扁平, 溫度邊界層逐漸被擠壓到左右壁面, 其對流換熱強度增大; 當?= 0.7 時, 流體受到的約束阻力進一步變小, 使得其對流換熱區(qū)域擴大; 當?= 0.9 時, 即在高孔隙率條件下, 流胞在方腔中心處轉變?yōu)楸馄降臋E圓形流胞, 并且主流胞中心存在分裂成兩個小流胞的趨勢, 方腔中心區(qū)域存在大量與水平壁面平行的等溫線, 因此方腔內對流換熱占據主導作用. 值得注意的是, 由于左側垂直壁面為正弦溫度分布, 在方腔左上角位置存在壁面吸收熱的現(xiàn)象, 導致在左上角處形成了一個相反于主對流胞的逆時針小漩渦.

    為了分析孔隙率對流體流速的影響, 本研究組還比較了不同?下熱壁面處局部V和上壁面處局部U的分布圖象. 如圖3 所示, 隨著孔隙率的增大, 熱壁面附近的豎直流速明顯增大; 同樣地, 上壁面附近的水平速度也存在相同的增長趨勢. 由于為正弦溫度分布, 熱浮升力在熱壁面中點位置最強, 因此豎直速度的峰值在Y = 0.5 附近; 與此同時, 由于左上角小流胞的阻礙作用, 上壁面附近水平速度的峰值出現(xiàn)在X = 0.8 附近. 值得注意的是, 由于方腔的左上角存在一個逆時針旋轉的小流胞, 因此, 在其附近的豎直和水平速度均為負值.綜上所述, 隨著孔隙率的增大, 流體受到多孔介質的阻力減弱導致流速增強, 提升了熱壁面處的換熱效率.

    圖2 不同 ? 下溫度場和流場的分布圖像 (a) ? = 0.3; (b) ? = 0.5; (c) ? = 0.7; (d) ? = 0.9Fig. 2. Streamlines, isotherms contours for different ? : (a) ?= 0.3; (b) ?= 0.5; (c) ?= 0.7; (d) ? = 0.9.

    圖3 (a)不同 ? 下X = 0 處的豎直速度分布; (b) Y = 1 處的水平速度分布圖Fig. 3. (a) Vertical velocity distribution at X = 0; (b) horizontal velocity distribution at Y = 1 for different ? .

    圖4 (a)不同 ? 下熱壁面處Nuave 數(shù)分布曲線; (b)熱壁面處局部Nu 數(shù)分布曲線Fig. 4. (a) At the heated wall Nuave number; (b) local Nu number for different ? .

    如圖4(a)所示, 隨著孔隙率的增大, Nuave數(shù)呈現(xiàn)上升趨勢, 然而這種增長幅度逐漸減弱, 這一變化規(guī)律與以往均勻溫度分布所得的結論相符[35].這是由于在較高孔隙率的情況下, 方腔內部多孔介質的阻礙作用較小, 孔隙率對方腔內傳熱效率的提升不明顯. 為了進一步分析孔隙率對熱壁面附近局部換熱特性的影響, 本研究組還分析了熱壁面處的局部Nu 數(shù)隨孔隙率變化的情況, 如圖4(b)所示,熱壁面0.2 < Y < 0.6 為主換熱區(qū)域, 不同孔隙率之間的變化差異并不明顯. 與此同時, 在熱壁面0.6 < Y < 0.8, 隨著 ? 的增大, 局部Nu 數(shù)逐漸減小, 這是由于位于方腔中心的主流胞對流逐漸增強, 流胞形態(tài)逐步轉變?yōu)槠綑E圓形, 使得溫度邊界層沿著主流胞流動方向擴散, 因此, 在該區(qū)域內溫度邊界層逐漸變厚, 局部換熱效率下降. 但是在0 < Y < 0.2 和0.9 < Y < 1 靠近上下壁面區(qū)域附近, 其變化趨勢恰好相反, 這是由于最大溫度在熱壁面中點位置, 向上或向下壁面延伸時熱壁面溫度逐漸減小, 在上下壁面局部形成了左側低溫右側高溫的現(xiàn)象, 因此在方腔左上角和左下角出現(xiàn)了熱壁面向方腔內部吸熱的現(xiàn)象, 并且形成了逆時針方向的小漩渦. 隨著孔隙率的變大, 主流胞逐漸擴展至幾乎整個方腔區(qū)域, 上下兩個小流胞逐漸被擠壓變小, 向左換熱的區(qū)域逐漸被壓縮并且換熱效率變弱, 局部Nu 數(shù)的絕對值逐漸變小主導了Nuave數(shù)的發(fā)展趨勢, 因此, 隨著孔隙率的增大, 流胞形態(tài)發(fā)生轉變, 減少了熱壁面從方腔內部的吸熱量, 進一步促進了方腔內的對流換熱效率, 使得熱壁面處的Nuave數(shù)不斷增大.

    4.2 Ra 數(shù)對方腔內自然對流的影響

    在自然對流問題中, Ra 數(shù)是表征自然對流綜合強弱的重要無量綱準則數(shù), 為了研究Ra 數(shù)對方腔內自然對流的影響, 保持 ? =0.9, Da = 0.01, f =0.01, Pr = 6.2 不變, Ra 數(shù)在103—106變化. 如圖5所示, 當Ra = 103時, 對流換熱現(xiàn)象幾乎可以忽略不計, 圓滑的主流胞集中在靠近熱壁面區(qū)域內,大量等溫線呈現(xiàn)上下對稱并垂直于水平方向, 溫度梯度大部分集中在方腔左半側區(qū)域內, 因此該狀態(tài)下流場內主要以導熱為主; 當Ra = 104時, 方腔中心處形成了一個橢圓形狀的對流胞, 由于浮升力的增強, 流體在進入熱壁面下側時促進了對流換熱, 而流體在進入熱壁面上側時等溫線將沿著流動方向擴散, 導致局部換熱能力減弱, 此時方腔內換熱方式仍然以導熱為主; 當Ra = 105時, 方腔中心處存在大量水平方向的等溫線, 方腔中心的流胞開始出現(xiàn)向多個流胞分裂的趨勢, 因此對流區(qū)域內主要以對流換熱為主; 當Ra = 106時, 此時方腔中心處的橢圓形流胞分裂成2 個不規(guī)則小漩渦, 主流胞占據了幾乎整個多孔介質方腔, 邊界層被擠壓變薄, 并且邊界層附近存在很多沿流體流動方向的等溫線. 值得注意的是, 隨著Ra 數(shù)的不斷增大, 方腔左上方的對流胞在逐漸增大, 其中心位置不斷向右側移動, 擠占主流胞的空間.

    本研究組還比較了不同Ra 數(shù)下熱壁面處豎直速度和上壁面處水平速度分布圖像. 如圖6(a)和圖6(b)所示, 在低Ra 數(shù)時熱壁面處的流速十分微弱, 在高Ra 數(shù)時熱壁面處流體流速得到了顯著提升, 這意味著Ra 數(shù)的增加將直接增強方腔內自然對流的強度, 腔體內流胞環(huán)流變得越來越劇烈,流動速率加快, 溫度邊界層變薄, 強化了壁面處的傳熱能力, 使得方腔內部溫度變得越來越均勻. 隨著Ra 數(shù)的增大, 在主換熱區(qū)域流速增大的同時,方腔左上角小漩渦的流動速率也得到增大, 并且速度為負的區(qū)域明顯擴大.

    如圖7(a)所示, 均勻溫度分布模型表明熱壁面的Nuave數(shù)隨Ra 數(shù)增大而增大[20], 但是隨著Ra 數(shù)的增大, 當Ra = 104時, 熱壁面Nuave數(shù)出現(xiàn)了微弱減小現(xiàn)象. 為了進一步解釋這一現(xiàn)象, 本研究組分析了局部Nu 數(shù)的分布曲線. 如圖7(b)所示, 隨著瑞利數(shù)的增大, 局部Nu 數(shù)的峰值位置集中在0.2 < Y < 0.5 附近, 這是由于熱浮升力的增強使得流速加快, 下側熱邊界層變薄, 并且在高瑞利數(shù)時該換熱區(qū)域占據主導地位, 因此, 在Ra =105和106時, 局部Nu 數(shù)產生了顯著的增長趨勢.

    圖5 不同Ra 下溫度場和流場的分布圖像 (a) Ra = 103; (b) Ra = 104; (c) Ra = 105; (d) Ra = 106Fig. 5. Streamlines, isotherms contours for different Ra number: (a) Ra = 103; (b) Ra = 104; (c) Ra = 105; (d) Ra = 106.

    圖6 (a) 不同Ra 下X = 0 處的豎直速度分布; (b) Y = 1 處的水平速度分布圖Fig. 6. (a) Vertical velocity distribution at X = 0; (b) horizontal velocity distribution at Y = 1 for different ? .

    但是, 在0.6 < Y < 0.8, 當Ra = 103和104時局部會出現(xiàn)Nu 數(shù)被削弱的現(xiàn)象, 究其原因, 一方面,是由于流胞逐漸轉變?yōu)楸馄降臋E圓形, 對左上側的熱邊界層擠壓作用減弱; 另一方面, 隨著瑞利數(shù)的增大, 中心流胞流速增大并轉變?yōu)闄E圓形狀, 由于流體流速的增強, 該區(qū)間內等溫線沿流胞主流區(qū)域流動方向擴散, 使得該局部溫度邊界層變厚. 此外,由于熱壁面為非均勻邊界條件, 在上下壁面附近處出現(xiàn)了局部Nu 數(shù)小于0 的情況, 隨著瑞利數(shù)的增加, 中心區(qū)域的換熱效率增強, 同時上壁面附近向左的熱通量也得到加強, 而下壁面受到主流胞擠壓作用向左側的換熱效率微弱下降.

    通過上述分析發(fā)現(xiàn), 當Ra = 104時, 熱壁面上側局部邊界層沿流體流動方向擴散, 局部換熱效率下降顯著, 并且低溫區(qū)域吸熱現(xiàn)象隨Ra 數(shù)增大而增大, 導致左壁面的Nuave數(shù)出現(xiàn)了減小現(xiàn)象.與此同時, 在高Ra 數(shù)情況下, 由于方腔內主換熱區(qū)域換熱效率增長更加明顯, 其壁面Nuave數(shù)增長顯著, 受非均勻溫度邊界層影響較小.

    圖7 (a)不同下Ra 熱壁面處Nuave 數(shù); (b)熱壁面處局部Nu 數(shù)分布曲線Fig. 7. (a) At the heated wall Nuave number; (b) local Nu number for different Ra.

    4.3 傾斜角g 對流體流動換熱的影響

    當多孔介質方腔逆時針旋轉一定角度時, 流體的流動形態(tài)和換熱強度將會發(fā)生明顯變化. 為了分析 傾 斜 角g 的 影 響,保 持 ? = 0.9, Da = 0.01,?=0.01 , Ra = 105, Pr = 6.2 不變, 傾斜角g 在0° ≤ g ≤ 120°變化. 從圖8 可知, 當g ≠ 0°時, 方腔中心流胞的流線更加圓滑, 主流胞也更加向中心集中, 并且方腔中心為單一流胞, 未出現(xiàn)分裂成多個流胞的趨勢. 從溫度場可知, 當方腔施加傾角時,在方腔中心形成了大量沿流胞主流速方向的彎曲等溫線, 但是其溫度邊界層的變化在宏觀云圖中變化不明顯.

    因此, 本研究組計算了水平中心線的局部溫度分布以及熱壁面處Vave/Nuave的對比圖, 其中,Vave為沿Y 方向的平均速度. 如圖9(a)所示, 與g = 0°相比, 熱壁面附近溫度下降更為迅速, 這意味著在有傾斜角的情況下熱壁面中點附近的溫度邊界層變薄, 換熱效率得到一定加強. 深入分析熱壁面整體的換熱效率. 如圖9(b)所示, 當g =40°時, 發(fā)現(xiàn)壁面處的換熱效率得到了較為顯著的提升; 當g = 80°時, Nuave數(shù)存在削弱換熱現(xiàn)象;當g = 120°時, 其強化傳熱作用較為微弱, 這一變化規(guī)律與以往文獻[36]所得結論一致. 由于傾斜角的變化使得熱浮升力的方向發(fā)生改變, 流體在X 方向上存在浮升力的驅動, 影響了熱壁面處Vave的大小, 并且壁面處局部V 會直接影響溫度梯度以及壁面處傳熱強度, 因此Vave/Nuave的變化趨勢基本相同, 并呈現(xiàn)上下波動的趨勢.

    為了進一步探究g 對流體傳熱的影響, 本研究組分析了熱壁面處局部V 的分布曲線, 由于g =120°時旋渦流動方向為逆時針, 本文選取它的相反數(shù)作為參考. 從圖10(a)可以看出, 當g = 40°時,其局部V 的峰值存在最大值; 當g = 80°時, 其局部V 的峰值最小, 并且絕大多數(shù)區(qū)域小于g = 0°時的數(shù)值; 當g = 120°時, 速度峰值數(shù)值與g = 0°時幾乎接近, 因此g = 40°時強化效果最為明顯,而g = 80°時則產生削弱現(xiàn)象. 最后, 為了進一步細致分析g = 120°時局部Nu 數(shù)的分布規(guī)律, 本研究組給出了不同傾角下局部Nu 數(shù)的分布曲線. 如圖10(b)所示, 當g = 0°, 40°和80°時, 非均勻溫度邊界條件在中點以上位置出現(xiàn)了與主換熱區(qū)域相反趨勢的分布情況, 而最大局部Nu 數(shù)存在于熱壁面中點以下區(qū)域. 這是由于當流體流動方向為順時針時, 中點以下的流體受到浮升力和主流包擠壓作用的影響, 對流換熱作用最顯著. 由圖10(b)可知,主換熱區(qū)域內的換熱強度由強到弱依次為40°,0°和80°. 當g = 120°時, 流胞旋轉方向的轉變使得熱壁面最大局部Nu 數(shù)出現(xiàn)在壁面中點以上的位置, 而局部Nu 數(shù)的峰值介于0°和40°之間.

    對方腔施加傾角后, 熱浮升力的有效作用發(fā)生改變, 使得壁面的豎直速度發(fā)生改變, 而壁面Vave與Nuave數(shù)關聯(lián)性增強. 當g = 40°時, 通過提升主換熱區(qū)域的豎直速度可以提高該區(qū)域內局部Nu 數(shù), 使得整體壁面換熱效率增強. 類似地, 當g = 80°時, 壁面豎直速度的削弱對流體傳熱效率的影響較為明顯, 導致壁面換熱效率下降顯著. 這一規(guī)律與文獻[36]所得結論相符, 表明主換熱區(qū)域對壁面Nuave數(shù)影響較為明顯, 壁面低溫區(qū)域從方腔內吸熱的強度受傾斜角的影響較小, 這是由于主換熱區(qū)域和主流區(qū)域近似重疊, 因此, 隨傾斜角的變化, 均勻或非均勻邊界條件下主換熱區(qū)域內Nuave數(shù)強化或削弱作用都將起到主導作用.

    圖8 不同g 下溫度場和流場的分布圖像 (a) g = 0°; (b) g = 40°; (c) g = 80°; (d) g = 120°Fig. 8. Streamlines, isotherms contours for different g number: (a) g = 0°; (b) g = 40°; (c) g = 80°; (d) g = 120°.

    圖9 (a)不同g 時Y = 0.5 處局部溫度分布曲線; (b)熱壁面處Vave/Nuave 的分布曲線Fig. 9. (a) Local temperature distribution along the Y = 0.5; (b) average velocity in the y direction & Nuave number at the heated wall in different g.

    圖10 (a)不同g 下熱壁面處局部豎直速度V; (b)熱壁面處局部Nu 數(shù)的分布曲線Fig. 10. (a) Local velocity in the y direction; (b) local Nuave number at the heated wall in different g.

    本研究組還研究了不同傾斜角下孔隙率的增加對方腔內強化傳熱的影響. 如圖11(a)所示,在施加傾斜角的情況下, 隨著孔隙率的增大, 當? = 0.9 時, Nuave數(shù)出現(xiàn)輕微減小; 而沒有施加傾角時, Nuave數(shù)隨 ? 的增大而增大. 當g = 40°時,Nuave數(shù)的增長幅度最大, 當 ? = 0.7 時, Nuave數(shù)達到 峰 值,與g = 0°的 情 況 相 比, Nuave數(shù) 在 ? =0.7 時提升了17.72%.

    當方腔不施加傾角時, 壁面處的Nuave數(shù)隨孔隙率增加而增加, 但是在施加傾斜角條件時, 方腔熱壁面的Nuave數(shù)出現(xiàn)先增加后減小的趨勢, 為了解釋這一現(xiàn)象, 以g = 40°和0°時孔隙率為0.7 和0.9 的情況為例進行分析. 如圖11(b)所示, 當g =0°時, 在主換熱區(qū)域內(0.2 < Y < 0.6), 孔隙率增大, 局部Nu 數(shù)減小的幅度較小, 而位于上壁面附近區(qū)域內, 局部Nu 數(shù)的增長成為主導. 局部Nu 數(shù)是由導熱系數(shù)比值km/kf和熱壁面溫度梯度所組成, 當g = 0°時, 左壁面溫度梯度的增加幅度相對于km/kf減小幅度更加顯著, 因此, 隨著孔隙率的增加, Nuave數(shù)呈現(xiàn)出單調遞增的趨勢. 當g =40°時, 方腔內對換熱能力起主導作用的依然位于主換熱區(qū)域(0.2 < Y < 0.6), 在此區(qū)域內, 隨著孔隙率的增大, 局部Nu 數(shù)減小的幅度更為明顯, 從而導致了孔隙率為0.9 時Nuave數(shù)的輕微減小. 這是由于隨著孔隙率的進一步增大, 多孔介質導熱對方腔傳熱效率的貢獻逐漸變小, 即km隨孔隙率增大而減小, 而施加不同傾斜角后, 左壁面處的溫度梯度在0.7 < ? < 0.9 增長的幅度減緩, 傾斜角的變化使得km/kf的減小幅度相比于左壁面溫度梯度的增加幅度更加顯著, 因此加熱壁面的局部Nu 數(shù)出現(xiàn)了減小的情況. 值得注意的是, 當g = 80°與g = 0°進行對比時可知, 熱壁面處的Nuave數(shù)在? =0.9 時產生了削弱的作用. 這可能是由于在低孔隙率的情況下, 沿X 方向的作用力可以更有效地促進流場的對流換熱, 而隨著孔隙率的增加, 方腔內流速加快, 在熱壁面沿Y 方向上熱浮升力不再起到主導作用, 從而使得熱壁面的換熱效率低于g = 0°時的換熱效率. 通過分析發(fā)現(xiàn), 在高孔隙率時對方腔施加不同的傾角, 隨孔隙率增大, 壁面溫度梯度增長幅度以及多孔介質對流體流速的影響減弱, 使得熱壁面處的傳熱能力得到抑制, 并且當? =0.7 時左壁面處的Nuave數(shù)達到峰值, 這表明在高孔隙率時, 壁面處的Nuave數(shù)受到km/kf減小的影響將更加明顯.

    圖11 (a)隨著 ? 的增加不同g 時熱壁面Nuave 數(shù)分布曲線; (b)當g = 0°, 40°時, 不同 ? 下局部Nu 數(shù)的分布曲線Fig. 11. (a) Variation of Nuave number as a function of ? in different g at the heated wall; (b) when g = 0°, 40°, variation of local Nu number at the heated wall in different ? .

    圖12 (a)隨著f 的增加不同g 下熱壁面Nuave 數(shù)分布曲線; (b)當g = 0°, 40°時, 不同f 下局部Nu 數(shù)的分布曲線Fig. 12. (a) Variation of Nuave number as a function of f in different g at the heated wall; (b) when g = 0°, 40°, variation of local Nu number at the heated wall in different f.

    最后, 本研究組分析了不同g 下隨著納米顆粒體積分數(shù)的增加, 方腔內部的換熱機理. 如圖12(a)所示, 隨著納米顆粒體積分數(shù)的增加, 當g = 40°,80°和120°時, 熱壁面處的Nuave數(shù)存在微弱增大的趨勢, 且當g = 40°時, Nuave數(shù)最大, 強化換熱現(xiàn)象最明顯, 而當g = 0°時, Nuave數(shù)反而出現(xiàn)了微弱削弱的趨勢, 且 ?=0.04 較 ? =0 時的Nuave數(shù)微弱減小了0.32%. Ho 等[20]采用同樣的導熱系數(shù)和粘度公式進行方腔內數(shù)值模擬, 發(fā)現(xiàn)隨著Al2O3-water 納米流體體積分數(shù)的增加, 熱壁面的Nuave數(shù)緩慢增加. 導致本文的結果與文獻[26]的差異主要由以下兩個方面原因: 一方面, 由于本文采用非均勻溫度邊界, 在增強高溫壁面換熱效率的同時, 由于靠近上下兩個壁面處存在從方腔內部吸熱的現(xiàn)象, 這部分區(qū)域也隨之被加強并在g = 0°時占據主導趨勢; 另一方面, 在g = 0°時, 隨體積分數(shù)的增加流體粘度變大, 熱壁面的平均溫度梯度緩慢減小, 這一減小幅度相比km/kf的增長幅度更為顯著, 導致Nuave數(shù)出現(xiàn)緩慢下降的現(xiàn)象.

    但是, 當對方腔施加任意傾角時, 熱壁面Nuave數(shù)緩慢增長, 為了進一步分析這種不同變化趨勢的原因, 在圖12(b)中, 以g = 0°和40°, 納米顆粒體積分數(shù)為0%和4%的情況為例進行分析,一方面, 當g = 40°時, 在主換熱區(qū)(0.2 < Y < 0.6)附近, 隨著孔隙率的增大, 局部Nu 數(shù)隨之增長;另一方面, 雖然隨體積分數(shù)的增大流體, 有效粘度逐漸變大導致流體流速降低, 對流強度緩慢下降,但是施加不同的傾斜角后, 壁面溫度梯度呈現(xiàn)微弱增加趨勢, 而km/kf增長幅度固定不變, 使得在g = 40°時Nuave數(shù)單調遞增趨勢更加明顯. 這說明通過對方腔施加一定傾斜角, 改變了壁面溫度梯度變化趨勢, 進而改變了其Nuave數(shù). 因此, 熱壁面的換熱效率取決于g 和 ? 對壁面溫度梯度的影響.

    5 結 論

    納米流體作為良好的換熱媒介被廣泛地應用到強化傳熱領域, 很多因素會影響其強化傳熱特性, 其中包括多孔介質孔隙率、納米流體體積分數(shù)、方腔傾斜角度以及邊界條件等因素. 本文以非均勻溫度邊界條件為例, 采用LBM 系統(tǒng)研究方腔內Al2O3-H2O 納米流體的對流傳熱特性, 重點討論了不同g 對腔內換熱效率的作用, 以及 ?和 ? 對不同傾斜角下多孔介質方腔內流動換熱的影響. 通過對非均溫度邊界條件的情況研究表明, Ra 數(shù)、孔隙率和納米顆粒體積分數(shù)并不像均勻溫度邊界條件那樣始終起到強化傳熱的效果.

    主要結論如下:

    1)在低Ra 數(shù)時, 流體受非均勻溫度邊界層低溫區(qū)域影響出現(xiàn)吸熱現(xiàn)象, 熱壁面處的Nuave數(shù)出現(xiàn)輕微減小的趨勢, 這與以往采用均勻溫度分布的文獻[20]所得結論存在差異; 在高Ra 數(shù)時, 由于瑞利數(shù)綜合表現(xiàn)為熱浮升力的強度, 方腔內主換熱區(qū)域換熱效率增長更加明顯, 熱壁面處Nuave數(shù)增大最為顯著, 受到非均勻溫度邊界層影響較小.

    2)隨著g 的增加, 方腔內熱浮升力隨之發(fā)生改變, 通過對壁面流速的改變, 對其壁面強化換熱效果產生不同影響, 當g = 40°時, 左壁面的Nuave數(shù)增長幅度較為明顯; 而當g = 80°時, 左壁面的Nuave數(shù)反而產生削弱作用. 這一規(guī)律與文獻[36]所得結論相符, 這是由于主換熱區(qū)域和主流區(qū)域近似重疊, 而壁面Vave與 Nuave數(shù)關聯(lián)性較強, 因此壁面Nuave數(shù)受主換熱區(qū)域影響較為明顯, 而壁面吸熱區(qū)域受到傾斜角的影響較小.

    3)均勻溫度邊界條件下的研究表明, 壁面處的Nuave數(shù)隨孔隙率增加而增加[12]; 而非均勻溫度邊界條件下的結果表明, 在高孔隙率對方腔施加不同傾角時, 增大孔隙率對壁面溫度梯度的促進作用減弱, 因此, 孔隙率的增加對熱壁面處的傳熱能力存在抑制作用.

    4)對于均勻溫度分布邊界條件下的研究表明,熱壁面的Nuave數(shù)隨Al2O3-H2O 納米流體體積分數(shù)增加而緩慢增加[20], 而本文非均勻溫度邊界條件下的研究表明, 當g= 0°時, 隨?的增大流體受粘性增長的影響較為明顯, 熱壁面溫度梯度緩慢下降,Nuave數(shù)存在微弱減小的趨勢. 而施加傾斜角后, 熱壁面溫度梯度受流體粘性增長的影響減弱,導致隨?的增大熱壁面溫度梯度微弱增大, 使Nuave數(shù)呈現(xiàn)微小增長趨勢.

    通過上述總結歸納發(fā)現(xiàn), 在本文所研究的非均勻溫度邊界物理模型下, 納米顆粒體積分數(shù)所產生的強化換熱作用不再明顯. 給定Ra數(shù)時, 要有效改善納米流體流動換熱效率, 就需要利用多孔介質對有效導熱系數(shù)的提升, 以及傾斜角對系統(tǒng)內的擾動作用, 采用合適的孔隙率和方腔傾斜角度對方腔進行干預.

    猜你喜歡
    熱效率對流壁面
    齊口裂腹魚集群行為對流態(tài)的響應
    二維有限長度柔性壁面上T-S波演化的數(shù)值研究
    壁面溫度對微型內燃機燃燒特性的影響
    提高蒸汽系統(tǒng)熱效率
    基于ANSYS的自然對流換熱系數(shù)計算方法研究
    二元驅油水界面Marangoni對流啟動殘余油機理
    豐田汽車公司的新型高熱效率汽油機
    顆?!诿媾鲎步Ec數(shù)據處理
    考慮裂縫壁面?zhèn)Φ膲毫丫a能計算模型
    基于對流項的不同非線性差分格式的穩(wěn)定性
    91久久精品国产一区二区三区| 国产精品一区www在线观看| 亚洲av男天堂| 少妇人妻一区二区三区视频| 亚洲激情五月婷婷啪啪| 亚洲激情五月婷婷啪啪| 久久久午夜欧美精品| 亚洲自拍偷在线| 欧美日韩一区二区视频在线观看视频在线 | 高清午夜精品一区二区三区| 亚洲成色77777| 人体艺术视频欧美日本| 成人毛片60女人毛片免费| 天天躁夜夜躁狠狠久久av| 亚洲自偷自拍三级| 精品欧美国产一区二区三| 日日啪夜夜爽| 女人久久www免费人成看片| 国产乱人视频| 淫秽高清视频在线观看| 亚洲精品第二区| 99热这里只有是精品50| 日本免费在线观看一区| av在线蜜桃| 老司机影院毛片| 人妻一区二区av| 尾随美女入室| 精品久久久久久久人妻蜜臀av| 亚洲丝袜综合中文字幕| 久久久精品免费免费高清| 中文字幕亚洲精品专区| 九九爱精品视频在线观看| 欧美成人精品欧美一级黄| 99久久精品国产国产毛片| 51国产日韩欧美| 亚洲国产精品成人综合色| 国产精品女同一区二区软件| 免费观看无遮挡的男女| 午夜精品在线福利| 国产视频内射| 免费黄网站久久成人精品| 亚洲国产精品国产精品| 亚洲丝袜综合中文字幕| 80岁老熟妇乱子伦牲交| 国产毛片a区久久久久| 永久网站在线| 91aial.com中文字幕在线观看| 亚洲成人精品中文字幕电影| 国产不卡一卡二| 国产 亚洲一区二区三区 | 国产极品天堂在线| 免费电影在线观看免费观看| 嫩草影院精品99| 69av精品久久久久久| 性插视频无遮挡在线免费观看| 简卡轻食公司| 中文字幕亚洲精品专区| 免费看不卡的av| 成人高潮视频无遮挡免费网站| 在线天堂最新版资源| 不卡视频在线观看欧美| 欧美一级a爱片免费观看看| 久久亚洲国产成人精品v| 久久久久国产网址| 91狼人影院| 搞女人的毛片| 人妻夜夜爽99麻豆av| 两个人的视频大全免费| 亚洲av国产av综合av卡| 国产av码专区亚洲av| 色综合色国产| 在线a可以看的网站| 免费观看在线日韩| 久久鲁丝午夜福利片| 免费黄网站久久成人精品| 亚洲精品日本国产第一区| 国产综合精华液| 久久久久久伊人网av| 国产三级在线视频| 国产黄片美女视频| 久久这里有精品视频免费| 亚洲久久久久久中文字幕| 亚洲精品国产av成人精品| 久久精品久久精品一区二区三区| 高清av免费在线| 久久久久精品久久久久真实原创| 亚洲av国产av综合av卡| 久久精品国产亚洲网站| 成人特级av手机在线观看| 七月丁香在线播放| 男女边摸边吃奶| 搞女人的毛片| 黄色欧美视频在线观看| 一级av片app| 一本久久精品| 久久6这里有精品| 丰满少妇做爰视频| 成人鲁丝片一二三区免费| 成人美女网站在线观看视频| 日本免费在线观看一区| 99久久精品热视频| 国产黄色小视频在线观看| 免费黄色在线免费观看| 国产精品一区www在线观看| 久久国产乱子免费精品| 亚洲欧美日韩无卡精品| 黄片无遮挡物在线观看| 日本黄大片高清| 在线a可以看的网站| 国产成人freesex在线| 九草在线视频观看| 免费观看av网站的网址| 1000部很黄的大片| 国产精品一区二区性色av| 免费看不卡的av| 狠狠精品人妻久久久久久综合| 成年人午夜在线观看视频 | 国产精品一区二区三区四区久久| 免费看a级黄色片| 婷婷色麻豆天堂久久| 亚洲自拍偷在线| 美女国产视频在线观看| 免费电影在线观看免费观看| 免费不卡的大黄色大毛片视频在线观看 | 国产老妇女一区| 又黄又爽又刺激的免费视频.| 又爽又黄无遮挡网站| 免费黄网站久久成人精品| 国产欧美日韩精品一区二区| 日韩av在线免费看完整版不卡| 国产av码专区亚洲av| 国产黄色免费在线视频| 亚洲国产精品专区欧美| 最近视频中文字幕2019在线8| 久久6这里有精品| 亚洲aⅴ乱码一区二区在线播放| 久久韩国三级中文字幕| 综合色丁香网| 久久久久久国产a免费观看| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 亚洲国产成人一精品久久久| 国产久久久一区二区三区| 观看免费一级毛片| 亚洲国产日韩欧美精品在线观看| 精品少妇黑人巨大在线播放| 人妻制服诱惑在线中文字幕| 国产又色又爽无遮挡免| 内地一区二区视频在线| 欧美高清性xxxxhd video| 中文字幕av在线有码专区| 日韩欧美一区视频在线观看 | 亚洲精品日本国产第一区| 久99久视频精品免费| 中文字幕制服av| 大又大粗又爽又黄少妇毛片口| 国产高清不卡午夜福利| 白带黄色成豆腐渣| 美女脱内裤让男人舔精品视频| 欧美日韩在线观看h| 视频中文字幕在线观看| 少妇人妻精品综合一区二区| 免费av观看视频| freevideosex欧美| 夫妻午夜视频| 五月天丁香电影| 久久6这里有精品| 日本午夜av视频| 看非洲黑人一级黄片| 久久这里有精品视频免费| 久久久欧美国产精品| 中文字幕av在线有码专区| 人体艺术视频欧美日本| 亚洲熟女精品中文字幕| 中文精品一卡2卡3卡4更新| 日本免费a在线| 最近的中文字幕免费完整| 久久精品久久久久久久性| 黄色欧美视频在线观看| 国产精品一区二区在线观看99 | 亚洲国产av新网站| 欧美一区二区亚洲| 国产精品一二三区在线看| 真实男女啪啪啪动态图| 亚洲av国产av综合av卡| 欧美成人一区二区免费高清观看| 日韩成人伦理影院| 亚洲自偷自拍三级| 国产91av在线免费观看| 亚洲真实伦在线观看| 国产成人a∨麻豆精品| 国产日韩欧美在线精品| 国产午夜精品久久久久久一区二区三区| 日日摸夜夜添夜夜爱| 国产免费又黄又爽又色| 日韩不卡一区二区三区视频在线| 麻豆精品久久久久久蜜桃| 成人亚洲欧美一区二区av| 九色成人免费人妻av| 国产精品蜜桃在线观看| 人妻夜夜爽99麻豆av| 日韩成人av中文字幕在线观看| av在线老鸭窝| 日本午夜av视频| 日韩亚洲欧美综合| 国产伦一二天堂av在线观看| av免费在线看不卡| 天堂俺去俺来也www色官网 | 嫩草影院新地址| 老女人水多毛片| 2022亚洲国产成人精品| 99久久精品国产国产毛片| 色5月婷婷丁香| 国产精品国产三级国产av玫瑰| 水蜜桃什么品种好| 女人久久www免费人成看片| 一边亲一边摸免费视频| 中文字幕人妻熟人妻熟丝袜美| 午夜日本视频在线| 国产在视频线精品| 成年版毛片免费区| 国产成人福利小说| 日韩欧美国产在线观看| 国产不卡一卡二| 久久久久精品性色| 久久精品国产亚洲av涩爱| 亚洲18禁久久av| 国产精品久久久久久久久免| 午夜福利在线在线| 日韩亚洲欧美综合| 精品国产露脸久久av麻豆 | 青春草视频在线免费观看| av在线蜜桃| 一级片'在线观看视频| 18禁裸乳无遮挡免费网站照片| 国产精品国产三级国产专区5o| 99热这里只有是精品50| av.在线天堂| 久久久久精品久久久久真实原创| 视频中文字幕在线观看| 99热全是精品| av在线亚洲专区| 国产亚洲5aaaaa淫片| 一区二区三区乱码不卡18| 91精品一卡2卡3卡4卡| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 内地一区二区视频在线| 3wmmmm亚洲av在线观看| 久久久久性生活片| 永久免费av网站大全| 性插视频无遮挡在线免费观看| 免费看不卡的av| 国产欧美另类精品又又久久亚洲欧美| 久久精品国产亚洲av天美| freevideosex欧美| 亚洲成人久久爱视频| 国产高清不卡午夜福利| 天美传媒精品一区二区| 久久亚洲国产成人精品v| 日韩三级伦理在线观看| 日日啪夜夜爽| 免费播放大片免费观看视频在线观看| 啦啦啦啦在线视频资源| 亚洲激情五月婷婷啪啪| 成人av在线播放网站| 国产白丝娇喘喷水9色精品| 色哟哟·www| 乱码一卡2卡4卡精品| 亚洲精品日韩av片在线观看| 亚洲精品成人av观看孕妇| 熟女电影av网| 亚洲婷婷狠狠爱综合网| 看免费成人av毛片| 免费大片18禁| 秋霞伦理黄片| 欧美精品国产亚洲| 国产大屁股一区二区在线视频| 哪个播放器可以免费观看大片| 久久久久精品性色| 国产精品一区二区性色av| 亚洲婷婷狠狠爱综合网| 免费无遮挡裸体视频| 国产高清国产精品国产三级 | 成人亚洲精品av一区二区| 亚洲av国产av综合av卡| 国产精品嫩草影院av在线观看| 成人亚洲精品一区在线观看 | 深夜a级毛片| 在线免费十八禁| 欧美不卡视频在线免费观看| 亚洲精品久久午夜乱码| 国产午夜精品一二区理论片| 哪个播放器可以免费观看大片| 大香蕉97超碰在线| 欧美激情国产日韩精品一区| 久久久久久久久大av| 午夜激情欧美在线| 日韩av在线免费看完整版不卡| 大香蕉97超碰在线| 在线天堂最新版资源| 少妇猛男粗大的猛烈进出视频 | 亚洲熟女精品中文字幕| 亚洲成人一二三区av| 成人毛片a级毛片在线播放| 欧美成人一区二区免费高清观看| 亚洲在久久综合| 99热网站在线观看| 国内精品宾馆在线| 80岁老熟妇乱子伦牲交| 中文字幕av成人在线电影| 日韩强制内射视频| 日本免费a在线| 中国国产av一级| 国产精品伦人一区二区| 国产精品一区二区三区四区免费观看| 国产成人91sexporn| 乱系列少妇在线播放| 免费电影在线观看免费观看| 日韩成人av中文字幕在线观看| 中文字幕亚洲精品专区| 日本熟妇午夜| 黄片wwwwww| 高清欧美精品videossex| 欧美一级a爱片免费观看看| 久久久午夜欧美精品| 成人午夜高清在线视频| 在线观看人妻少妇| 永久免费av网站大全| 97在线视频观看| 一级av片app| h日本视频在线播放| 少妇的逼好多水| 成人二区视频| 大陆偷拍与自拍| 国产av不卡久久| 午夜激情欧美在线| 亚洲经典国产精华液单| 亚洲国产欧美在线一区| 夜夜爽夜夜爽视频| 精品午夜福利在线看| 亚洲天堂国产精品一区在线| 丰满人妻一区二区三区视频av| 熟妇人妻不卡中文字幕| 日韩精品有码人妻一区| 亚洲美女搞黄在线观看| 婷婷六月久久综合丁香| 中文字幕人妻熟人妻熟丝袜美| 久久这里只有精品中国| 中文字幕免费在线视频6| 日本午夜av视频| 男人狂女人下面高潮的视频| 亚洲自拍偷在线| 天堂俺去俺来也www色官网 | 啦啦啦中文免费视频观看日本| 国产成人一区二区在线| 少妇裸体淫交视频免费看高清| 尾随美女入室| 观看美女的网站| 国产伦理片在线播放av一区| 久久综合国产亚洲精品| 天天一区二区日本电影三级| 一区二区三区四区激情视频| h日本视频在线播放| 91精品伊人久久大香线蕉| 欧美不卡视频在线免费观看| 午夜福利在线观看吧| 国产精品1区2区在线观看.| 亚洲aⅴ乱码一区二区在线播放| 免费大片黄手机在线观看| 最近中文字幕高清免费大全6| 一个人看视频在线观看www免费| 在线观看美女被高潮喷水网站| 日韩成人伦理影院| 一夜夜www| 欧美+日韩+精品| 日本免费a在线| 久久精品国产亚洲av涩爱| 一夜夜www| 久久鲁丝午夜福利片| 两个人视频免费观看高清| 欧美+日韩+精品| 观看美女的网站| 精品少妇黑人巨大在线播放| 91久久精品国产一区二区成人| 久久久久国产网址| 国产精品久久视频播放| 日本免费在线观看一区| 亚洲国产精品成人久久小说| 神马国产精品三级电影在线观看| 三级国产精品片| 街头女战士在线观看网站| 免费看日本二区| 日韩中字成人| 丝瓜视频免费看黄片| 婷婷色综合大香蕉| 国产片特级美女逼逼视频| 中文字幕久久专区| 日韩视频在线欧美| 亚洲精品日韩在线中文字幕| 欧美zozozo另类| 伦精品一区二区三区| 亚洲熟妇中文字幕五十中出| 日韩欧美国产在线观看| 又大又黄又爽视频免费| 色视频www国产| 人妻夜夜爽99麻豆av| 国产精品久久久久久av不卡| 少妇的逼水好多| 国产精品人妻久久久影院| 内地一区二区视频在线| 深夜a级毛片| 午夜激情福利司机影院| 亚洲欧美一区二区三区黑人 | 久久精品人妻少妇| 亚洲精品国产成人久久av| 99久久人妻综合| 国产亚洲5aaaaa淫片| 亚洲精品日本国产第一区| 色综合色国产| 最近视频中文字幕2019在线8| 欧美97在线视频| 99热网站在线观看| 中文字幕亚洲精品专区| 国产黄频视频在线观看| 欧美xxxx性猛交bbbb| av在线蜜桃| 一区二区三区乱码不卡18| 国产成人精品福利久久| a级毛片免费高清观看在线播放| 极品少妇高潮喷水抽搐| av免费观看日本| 午夜激情欧美在线| 国产精品久久久久久久电影| 最近中文字幕高清免费大全6| a级一级毛片免费在线观看| 国产 亚洲一区二区三区 | 啦啦啦啦在线视频资源| 国产永久视频网站| 国产在视频线精品| 成年版毛片免费区| 亚洲av一区综合| 国产在线男女| 国产精品一区二区三区四区久久| 老女人水多毛片| 国产av不卡久久| 亚洲精华国产精华液的使用体验| av在线播放精品| 老司机影院毛片| 五月天丁香电影| freevideosex欧美| 中文字幕免费在线视频6| 美女国产视频在线观看| 久久这里有精品视频免费| 美女xxoo啪啪120秒动态图| 日本熟妇午夜| 春色校园在线视频观看| 男人舔奶头视频| 欧美 日韩 精品 国产| 精品久久久久久电影网| 成人欧美大片| 校园人妻丝袜中文字幕| 欧美日韩综合久久久久久| 午夜精品在线福利| 久久99精品国语久久久| 亚洲不卡免费看| 国产成年人精品一区二区| 日韩欧美一区视频在线观看 | 亚洲色图av天堂| 久久久久久久久久成人| 夫妻午夜视频| 久久久精品免费免费高清| 国产一区二区亚洲精品在线观看| av在线老鸭窝| 久久99热这里只有精品18| 特大巨黑吊av在线直播| 亚洲熟女精品中文字幕| 超碰97精品在线观看| 国产爱豆传媒在线观看| 成人无遮挡网站| 国产黄片美女视频| 久久精品久久精品一区二区三区| 欧美+日韩+精品| 少妇裸体淫交视频免费看高清| 97超碰精品成人国产| 国产视频首页在线观看| 成人午夜精彩视频在线观看| 高清av免费在线| 中文字幕久久专区| 精品午夜福利在线看| 国产一区有黄有色的免费视频 | 欧美一区二区亚洲| 亚洲精品一二三| 精品一区在线观看国产| 国产伦精品一区二区三区四那| 高清av免费在线| 日韩欧美一区视频在线观看 | 久99久视频精品免费| 国产高潮美女av| 国产精品久久久久久精品电影小说 | 啦啦啦韩国在线观看视频| 青青草视频在线视频观看| 国产精品一区二区在线观看99 | 插逼视频在线观看| 亚洲精品成人av观看孕妇| 热99在线观看视频| 国产伦精品一区二区三区四那| 免费观看性生交大片5| 免费观看av网站的网址| 亚洲一区高清亚洲精品| 亚洲怡红院男人天堂| 日韩欧美一区视频在线观看 | 亚洲精品aⅴ在线观看| 又爽又黄无遮挡网站| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 免费看光身美女| 99久久精品国产国产毛片| 国产亚洲5aaaaa淫片| 尤物成人国产欧美一区二区三区| 一级毛片电影观看| 国产精品一区二区在线观看99 | 成年女人看的毛片在线观看| 麻豆av噜噜一区二区三区| 99热网站在线观看| 亚洲欧美成人精品一区二区| 欧美日韩国产mv在线观看视频 | 国产高潮美女av| 欧美 日韩 精品 国产| xxx大片免费视频| 国产探花极品一区二区| 国产精品熟女久久久久浪| 国产精品无大码| 亚洲人成网站高清观看| 校园人妻丝袜中文字幕| 国产免费福利视频在线观看| 极品教师在线视频| 国产国拍精品亚洲av在线观看| 人妻一区二区av| 日本一二三区视频观看| 中文欧美无线码| 国产精品麻豆人妻色哟哟久久 | 你懂的网址亚洲精品在线观看| 青春草国产在线视频| 色吧在线观看| 国产精品99久久久久久久久| 国产精品熟女久久久久浪| 亚洲精品日韩av片在线观看| 亚洲av成人精品一区久久| 天美传媒精品一区二区| 国产免费一级a男人的天堂| 欧美日韩亚洲高清精品| 一级a做视频免费观看| 白带黄色成豆腐渣| 97精品久久久久久久久久精品| 草草在线视频免费看| 久久久欧美国产精品| 亚洲高清免费不卡视频| 永久免费av网站大全| 亚洲精品成人av观看孕妇| 日本色播在线视频| 日本猛色少妇xxxxx猛交久久| 国产精品国产三级国产专区5o| 日韩av不卡免费在线播放| 亚洲精品日本国产第一区| 国产午夜福利久久久久久| 国产片特级美女逼逼视频| 大陆偷拍与自拍| 两个人视频免费观看高清| 亚洲最大成人中文| 欧美zozozo另类| 欧美日韩一区二区视频在线观看视频在线 | av黄色大香蕉| 亚洲18禁久久av| 美女被艹到高潮喷水动态| 亚洲av男天堂| 精品99又大又爽又粗少妇毛片| 成人国产麻豆网| 免费黄色在线免费观看| 午夜福利在线在线| 狂野欧美白嫩少妇大欣赏| 国产精品国产三级专区第一集| 男人狂女人下面高潮的视频| 亚洲国产av新网站| 免费观看精品视频网站| 一级二级三级毛片免费看| 成人亚洲精品一区在线观看 | 永久网站在线| 亚洲欧美一区二区三区国产| 街头女战士在线观看网站| 久久久久久久久久久丰满| 九色成人免费人妻av| 99re6热这里在线精品视频| 欧美成人一区二区免费高清观看| 欧美精品一区二区大全| 色视频www国产| 免费av毛片视频| 成人一区二区视频在线观看| 搡老妇女老女人老熟妇| 高清午夜精品一区二区三区| 久99久视频精品免费| 亚洲精品日韩av片在线观看| 成年av动漫网址| 五月玫瑰六月丁香| 尾随美女入室| 久久午夜福利片| 亚洲欧美一区二区三区黑人 | 午夜亚洲福利在线播放| 白带黄色成豆腐渣| 精品一区二区三区视频在线| 久久99蜜桃精品久久| 18禁动态无遮挡网站| av免费在线看不卡| 成年女人在线观看亚洲视频 | 欧美性感艳星| h日本视频在线播放| 在线 av 中文字幕| 日本与韩国留学比较|