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

    基于OpenFOAM的液態(tài)金屬鉛鉍三維流動換熱特性數(shù)值模擬研究

    2021-06-30 13:56:42何少鵬王明軍田文喜蘇光輝秋穗正
    原子能科學(xué)技術(shù) 2021年6期
    關(guān)鍵詞:邊界層液態(tài)湍流

    何少鵬,王明軍,*,章 靜,田文喜,蘇光輝,秋穗正

    (1.西安交通大學(xué) 動力工程多相流國家重點實驗室,陜西 西安 710049;2.西安交通大學(xué) 核科學(xué)與技術(shù)學(xué)院,陜西 西安 710049)

    液態(tài)金屬鉛鉍快堆作為第4代核能系統(tǒng)的主要堆型之一,有自然循環(huán)能力強、中子經(jīng)濟(jì)性好等諸多優(yōu)點。但液態(tài)金屬鉛鉍由于特殊的熱物性質(zhì)(如低普朗特數(shù)、高熱導(dǎo)率等),其流動邊界層與熱邊界層已不在一個量級,流動換熱特性已與水、空氣等常見工質(zhì)有很大不同[1],這給液態(tài)金屬鉛鉍反應(yīng)堆的研究設(shè)計帶來了許多困難。故亟需針對液態(tài)金屬鉛鉍開發(fā)更為精確的湍流換熱模型以進(jìn)行高精度數(shù)值模擬。

    在主流商用CFD軟件中,通常采取定值湍流普朗特數(shù)的SED模型對液態(tài)金屬流動換熱問題進(jìn)行模擬求解。該模型基于雷諾比擬原理,即假設(shè)熱邊界層與流動邊界層具有相似性,其成立的條件為Pr接近于1。對于空氣、水等常見工質(zhì),使用湍流普朗特數(shù)為0.85~1的SED模型往往能得到工程上可接受的結(jié)果。但針對液態(tài)金屬等低普朗特數(shù)流體開展的一些CFD模擬研究表明,定值湍流普朗特數(shù)的SED模型已不能達(dá)到很好的預(yù)測精度。Kawamura等[2]利用直接數(shù)值模擬研究了雷諾數(shù)與普朗特數(shù)對液態(tài)金屬流動換熱的影響,結(jié)果表明采用定值湍流普朗特數(shù)的SED模型已不能很好地與實驗關(guān)聯(lián)式吻合。Cheng等[3]通過CFD模擬研究了單通道內(nèi)液態(tài)金屬的熱工水力現(xiàn)象,發(fā)現(xiàn)當(dāng)SED模型的湍流普朗特數(shù)由0.9增加到1.5時,預(yù)測努塞爾數(shù)將會下降約30%。Chen等[4]使用一種新網(wǎng)格標(biāo)記分離方法研究了繞絲棒束組件三維流動換熱特性,結(jié)果表明繞絲會引起橫流以及軸向速度周期振蕩,影響熱工水力特性。Manservisi等[5]針對液態(tài)金屬開發(fā)了一種更為精確的k-ε-kθ-εθ四因子模型。該模型是一種顯式的代數(shù)熱通量模型,充分考慮了熱邊界層與流動邊界層的差異性,引入脈動溫度方均值及其耗散率兩個新參數(shù),建立了湍流熱擴散率與脈動溫度方均值間的數(shù)值關(guān)系,能實現(xiàn)對液態(tài)金屬湍流換熱的精確模擬。Manservisi等[6-7]使用四因子模型在三角形柵元通道和四邊形柵元通道內(nèi)進(jìn)行了數(shù)值模擬與驗證,均取得了較好的預(yù)測結(jié)果。

    考慮到主流商業(yè)軟件代碼閉源,封裝嚴(yán)格,不利于用戶對模型進(jìn)行修改和開發(fā),而開源計算流體力學(xué)軟件OpenFOAM采用C++語言編程,具有編程環(huán)境自由、模型開發(fā)簡便等優(yōu)點,廣泛應(yīng)用于能源化工、海洋船舶等領(lǐng)域?;陂_源平臺進(jìn)行模型添加和求解器開發(fā)對實現(xiàn)更高精度的液態(tài)金屬流動換熱CFD模擬有重要意義。故本文基于k-ε-kθ-εθ四因子模型在開源程序OpenFOAM中開發(fā)針對液態(tài)金屬的自定義求解器,闡述數(shù)值模擬中所使用的數(shù)學(xué)物理模型,包括守恒方程與四因子模型方程,并對帶繞絲棒束通道中液態(tài)金屬鉛鉍的流動換熱問題進(jìn)行模擬研究。

    1 數(shù)學(xué)物理模型

    1.1 守恒方程

    液態(tài)金屬反應(yīng)堆中液態(tài)金屬冷卻劑在正常工況下的流動可視為單相、不可壓縮的湍流流動。因此,經(jīng)過不可壓縮流簡化及雷諾時均處理后,動量、質(zhì)量、能量守恒方程如式(1)所示。

    (1)

    式中:U為速度平均值;T為溫度平均值;U′為湍流脈動速度;T′為湍流脈動溫度;cp為比定壓熱容;λ為熱導(dǎo)率。

    (2)

    Prt=νt/αt

    (3)

    湍流普朗特數(shù)給定后(本文分別取針對液態(tài)金屬推薦的湍流普朗特數(shù)1.5和2.0進(jìn)行模擬計算),借由湍流模型求得的局部湍流黏度νt,便可進(jìn)一步求解出湍流熱擴散系數(shù)αt。至此,方程組(1)得以封閉,方程組理論可解。

    1.2 k-ε-kθ-εθ四因子模型

    1) 低雷諾數(shù)修正的k-ε模型

    標(biāo)準(zhǔn)k-ε模型是一種較為常見的二方程湍流模型,適合模擬充分發(fā)展的湍流流動,對低雷諾數(shù)過渡流和近壁區(qū)域計算結(jié)果不夠理想。在四因子模型中,應(yīng)用了針對低雷諾數(shù)修正的k-ε方程。與標(biāo)準(zhǔn)k-ε模型相似,引入了Boussinesq假設(shè),認(rèn)為湍流雷諾應(yīng)力與牛頓剪切引力相似,均與應(yīng)變呈正比,雷諾應(yīng)力張量可寫為湍流黏度與速度梯度的乘積,即:

    (4)

    式中,τt為湍流雷諾應(yīng)力。

    與標(biāo)準(zhǔn)k-ε方程相似,低雷諾數(shù)修正的k-ε方程也定義了2個新的變量,即湍動能k及其耗散率ε,如式(5)所示。

    (5)

    與標(biāo)準(zhǔn)k-ε模型不同的是,本模型定義的湍流黏度為:

    νt=Cμkτlu

    (6)

    式中,Cμ為模型常數(shù)。為更加準(zhǔn)確地模擬流體的近壁面行為,該模型定義了局部動力學(xué)特征時間τlu[9-10],由4個參數(shù)計算得到,即:

    τlu=f1μA1μ+f2μA2μ

    (7)

    式(7)中4個參數(shù)均模化為歸一化壁面距離Rδ與湍流雷諾數(shù)Ret的函數(shù),如式(8)所示。

    (8)

    歸一化壁面距離Rδ和湍流雷諾數(shù)Ret定義如下:

    (9)

    式中:δ為各網(wǎng)格單元到壁面的距離;ν為層流運動黏性系數(shù)。

    為計算局部湍流黏度νt,需得到式(7)中定義的局部動力學(xué)特征時間τlu,即需要解出歸一化壁面距離和湍流雷諾數(shù),進(jìn)而需解出湍動能k及其耗散率ε。修正后的k-ε方程如式(10)所示。

    (10)

    式中:

    (11)

    (12)

    關(guān)于低雷諾數(shù)修正的k-ε模型的更多闡釋可參考文獻(xiàn)[11]。

    2)kθ-εθ模型

    與k-ε模型類似,kθ-εθ模型同樣定義了2個新變量:脈動溫度方均值kθ及其耗散率εθ。

    (13)

    與SED模型不同的是,在四因子模型中將湍流熱擴散率αt定義為式(14)所示的3個參數(shù)的乘積。

    αt=Cθkτlθ

    (14)

    式中:Cθ為模型常數(shù);τlθ為局部熱力學(xué)特征時間,用來更加精確地描述湍流流動中的換熱特性[5]。

    τlθ=f1θB1θ+f2θB2θ

    (15)

    式中,f1θ、f2θ、B1θ與B2θ4個參數(shù)被定義為歸一化壁面距離Rδ、湍流雷諾數(shù)Ret等變量的函數(shù)。

    (16)

    (17)

    式中,R為熱力學(xué)特征時間與動力學(xué)特征時間的比值。

    (18)

    綜上所述,為計算湍流熱擴散系數(shù)αt,需得到局部熱力學(xué)特征時間τlθ,進(jìn)一步需解出脈動溫度方均值kθ及其耗散率εθ,總結(jié)kθ-εθ方程組如下:

    (19)

    式中:Cd1=1、CP1=0.925、CP2=0.9均為模型常數(shù);Pk定義見式(11);Pθ與Cd2定義見式(20)。

    (20)

    2 帶繞絲棒束通道CFD模擬

    2.1 幾何模型

    本文選取德國卡爾斯魯厄理工學(xué)院(KIT)盒間流實驗[13]中的帶繞絲棒束組件作為CFD模擬和模型驗證對象。該棒束組件為六邊形,包含有7根帶繞絲燃料棒,每根棒表面施以均勻熱流密度以模擬燃料棒對鉛鉍合金流體加熱。從組件進(jìn)口向出口看,7根繞絲順時針旋轉(zhuǎn)并纏繞在7根棒上。該帶繞絲棒束組件具體幾何模型與截面視圖如圖1所示,主要幾何尺寸如下:P=20.5 mm、D=16.0 mm、d=4.4 mm、W=20.8 mm、FF=61.05 mm、L=262 mm、Dh=10.34 mm。

    圖1 帶繞絲棒束幾何模型及其截面視圖Fig.1 Geometric model and sectional view of wire-wrapped rod bundle

    2.2 計算方法

    本文選取繞絲螺距長度為262 mm的1個通道作為物理模型。進(jìn)口溫度均為573 K,施加速度入口邊界條件。出口為固定壓力條件,組件壁和繞絲表面設(shè)為絕熱壁面,棒束表面施加均勻的熱流密度。模擬中考慮了液態(tài)鉛鉍物性隨溫度的變化,梯度項離散采用Gauss linear格式,擴散項采用Bounded Gauss linear格式,其余使用一階迎風(fēng)差分格式,分別使用SED模型與四因子模型進(jìn)行計算,具體模擬工況列于表1。

    表1 模擬工況Table 1 Simulated condition

    2.3 網(wǎng)格無關(guān)性分析

    如圖1所示,本文選取組件出口截面上1條垂直線段A,以沿線段A的冷卻劑溫度變化特性及速度分布為對象獲得網(wǎng)格獨立性解,如圖2所示。從計算結(jié)果可看出,當(dāng)網(wǎng)格數(shù)為166萬、196萬和234萬時,沿線段A的冷卻劑速度、溫度變化相差不大,故最終選取166萬網(wǎng)格數(shù)的網(wǎng)格為網(wǎng)格獨立解開展后續(xù)數(shù)值模擬研究。

    圖2 網(wǎng)格無關(guān)性分析Fig.2 Grid independence analysis

    2.4 計算結(jié)果分析

    帶繞絲棒束通道內(nèi)各截面處的速度、溫度云圖如圖3所示。結(jié)果表明,由于繞絲對流體的攪混作用,各截面處速度分布呈現(xiàn)非對稱、不均勻特征,且沿軸向發(fā)生劇烈變化。由溫度云圖還可看出,通道內(nèi)不對稱速度分布引起了不對稱溫度分布。通道中心截面橫流速度場矢量圖如圖4所示,由于繞絲的攪混作用,橫流速度場呈現(xiàn)出旋轉(zhuǎn)特點,且橫流旋轉(zhuǎn)方向與繞絲旋轉(zhuǎn)方向一致。

    圖3 各截面速度與溫度云圖Fig.3 Velocity and temperature contour of several sections

    圖4 中間截面(z=131 mm)橫流速度矢量圖Fig.4 Cross flow velocity vector graph of middle section (z=131 mm)

    本文共選取3個典型局部通道(圖5),即中心通道、邊通道、角通道進(jìn)行模型驗證。選取截面位置為z=262 mm,即出口截面。針對液態(tài)金屬建立的Mikityuk換熱實驗關(guān)聯(lián)式[14]和Ushakov換熱實驗關(guān)聯(lián)式[15]如式(21)、(22)所示,式中X為柵距棒徑比,X=P/D。

    圖5 局部通道示意圖Fig.5 Local channel schematic diagram

    1) Mikityuk實驗關(guān)聯(lián)式:

    Nu=0.047(1-e-3.8(X-1))(Pe0.77+250)

    30≤Pe≤5 000;1.1≤X≤1.95

    (21)

    2) Ushakov實驗關(guān)聯(lián)式:

    1≤Pe≤4 000;1.3≤X≤2.0

    (22)

    關(guān)聯(lián)式(21)、(22)的實驗條件與本文模擬條件之間的對比列于表2。由表2可知,本文模擬條件與關(guān)聯(lián)式在棒束柵格幾何形狀、柵距棒徑比X、Pe等方面具有較好的一致性。盡管本文模擬幾何的X未在Ushakov關(guān)聯(lián)式推薦范圍內(nèi),但文獻(xiàn)[14]將X為1.15~1.95的實驗數(shù)據(jù)與Ushakov關(guān)聯(lián)式進(jìn)行了對比,結(jié)果表明,對于X<1.3的實驗數(shù)據(jù),Ushakov關(guān)系式也能符合很好。因此在本文模擬條件下,Mikityuk關(guān)聯(lián)式與Ushakov關(guān)聯(lián)式具有一定的可靠度。

    表2 關(guān)聯(lián)式實驗條件與模擬條件的對比Table 2 Comparison of correlations’ experimental condition and simulation condition

    使用湍流普朗特數(shù)取1.5和2.0的SED模型以及k-ε-kθ-εθ四因子模型的計算結(jié)果示于圖6,圖中M關(guān)聯(lián)式為Mikityuk關(guān)聯(lián)式預(yù)測值,U關(guān)聯(lián)式為Ushakov關(guān)聯(lián)式預(yù)測值。對比結(jié)果表明,四因子模型與實驗關(guān)聯(lián)式符合良好,在本文所模擬條件下能更為精確地預(yù)測繞絲棒束通道中的Nu。

    圖6 中心通道、邊通道、角通道內(nèi)努塞爾數(shù)Fig.6 Nusselt numbers of center, edge and corner channels

    圖7為使用SED模型(Prt=1.5)和四因子模型得到的沿圖1中直線A的湍流黏度νt與湍流熱擴散率αt的分布曲線。

    a——SED模型(Prt=1.5);b——四因子模型圖7 沿線段A的湍流黏度和湍流熱擴散率Fig.7 Turbulent viscosity and turbulent thermal diffusivity along line A

    從圖7a可看出,使用SED模型時,νt與αt的變化曲線具有高度一致性,兩者呈現(xiàn)一定的等比例關(guān)系。這是由于SED模型基于雷諾比擬原理,直接使用式(3)進(jìn)行湍流熱擴散率αt的計算,且由式(3)可知,νt與αt間的比例系數(shù)即為定值的湍流普朗特數(shù)1.5。然而對于低普朗特數(shù)的液態(tài)金屬鉛鉍,其動量擴散與能量擴散具有不同的行為特征[16],此時使用假設(shè)熱邊界層與流動邊界層相似的雷諾比擬原理往往會造成較大的誤差。

    與SED模型不同,四因子模型考慮了熱邊界層與流動邊界層的差異,引入了脈動溫度方均值及其耗散率2個參數(shù)來描述熱邊界層的能量耗散行為,構(gòu)造控制方程(式(19))求解湍流熱擴散率αt,從而避免了使用雷諾比擬原理。從圖7b可看出,νt和αt間并不呈完全的等比例關(guān)系,尤其在棒束壁面附近,表征動量擴散的νt與表征能量擴散的αt呈不同的變化趨勢,說明四因子模型能較好地模擬動量擴散與能量擴散的不一致性,更加準(zhǔn)確地描述液態(tài)金屬鉛鉍流場中的熱擴散行為。

    上述結(jié)果表明:四因子模型在預(yù)測繞絲棒束通道內(nèi)的Nu、描述液態(tài)鉛鉍流場動量和能量擴散行為方面相比傳統(tǒng)SED模型表現(xiàn)出了一定優(yōu)勢,本文基于四因子模型開發(fā)的適用于低普朗特數(shù)液態(tài)金屬流動傳熱模擬的OpenFOAM求解器可較好地預(yù)測該模擬工況下繞絲棒束通道中液態(tài)金屬鉛鉍流動換熱現(xiàn)象。

    3 結(jié)論

    本文基于開源CFD平臺OpenFOAM開發(fā)了適用于液態(tài)金屬模擬的k-ε-kθ-εθ四因子模型自定義求解器,對帶繞絲棒束通道中液態(tài)鉛鉍流動換熱現(xiàn)象進(jìn)行了模擬研究。計算中充分考慮了熱邊界層與流動邊界層的差異性,采用針對低雷諾數(shù)修正的k-ε模型求解湍流黏度,使用1組kθ-εθ偏微分方程求解湍流熱擴散率,實現(xiàn)了對液態(tài)金屬流動換熱問題較為精確的數(shù)值模擬。計算得到了帶繞絲棒束通道中液態(tài)金屬鉛鉍速度、溫度等重要熱工水力參數(shù)三維分布,深入分析了繞絲對棒束通道內(nèi)液態(tài)金屬流動換熱過程的影響規(guī)律,對比分析了四因子模型和傳統(tǒng)SED模型的模擬結(jié)果,證明了四因子模型在預(yù)測液態(tài)金屬鉛鉍流動傳熱過程中具有一定優(yōu)勢,上述模型和求解器能較好地模擬液態(tài)金屬鉛鉍在繞絲棒束中的流動換熱問題。

    猜你喜歡
    邊界層液態(tài)湍流
    液態(tài)金屬
    基于HIFiRE-2超燃發(fā)動機內(nèi)流道的激波邊界層干擾分析
    重氣瞬時泄漏擴散的湍流模型驗證
    2017年中外液態(tài)食品機械行業(yè)大事記
    淺談液態(tài)渣的顯熱利用和工藝技術(shù)
    資源再生(2017年3期)2017-06-01 12:20:59
    內(nèi)陸核電廠放射性液態(tài)流出物“近零排放”探討
    一類具有邊界層性質(zhì)的二次奇攝動邊值問題
    非特征邊界的MHD方程的邊界層
    “青春期”湍流中的智慧引渡(三)
    “青春期”湍流中的智慧引渡(二)
    色综合站精品国产| 亚洲四区av| 在线播放无遮挡| 日韩欧美三级三区| 97超级碰碰碰精品色视频在线观看| 精品久久久久久久末码| 亚洲一级一片aⅴ在线观看| 国产精品1区2区在线观看.| 国产v大片淫在线免费观看| 极品教师在线免费播放| 亚洲av成人av| 直男gayav资源| 久久精品夜夜夜夜夜久久蜜豆| 免费看日本二区| ponron亚洲| 亚洲最大成人手机在线| 一级黄色大片毛片| 久久6这里有精品| 国产精品亚洲美女久久久| 琪琪午夜伦伦电影理论片6080| 精品乱码久久久久久99久播| 老熟妇仑乱视频hdxx| 色综合婷婷激情| 麻豆精品久久久久久蜜桃| 久久久精品欧美日韩精品| 久久人人精品亚洲av| 啦啦啦啦在线视频资源| 久久久午夜欧美精品| 亚洲av第一区精品v没综合| 亚洲国产日韩欧美精品在线观看| 日韩强制内射视频| 成人精品一区二区免费| 国产一区二区激情短视频| 无人区码免费观看不卡| 日韩中文字幕欧美一区二区| 嫩草影院入口| 波野结衣二区三区在线| 欧美成人a在线观看| 国产午夜精品论理片| 免费大片18禁| 不卡视频在线观看欧美| 一夜夜www| 又黄又爽又免费观看的视频| 成人午夜高清在线视频| 国产精品国产高清国产av| 日韩av在线大香蕉| 精品福利观看| videossex国产| 亚洲国产高清在线一区二区三| 禁无遮挡网站| 国内精品美女久久久久久| 真人一进一出gif抽搐免费| 不卡一级毛片| 日韩一本色道免费dvd| 变态另类成人亚洲欧美熟女| 两性午夜刺激爽爽歪歪视频在线观看| 九九热线精品视视频播放| 国产亚洲精品av在线| 观看美女的网站| 日韩精品中文字幕看吧| 毛片一级片免费看久久久久 | 黄色一级大片看看| 岛国在线免费视频观看| 老女人水多毛片| 久久婷婷人人爽人人干人人爱| www.www免费av| 成人综合一区亚洲| 亚洲国产欧美人成| 免费av观看视频| 欧美日韩乱码在线| 看免费成人av毛片| 少妇猛男粗大的猛烈进出视频 | videossex国产| 在线观看免费视频日本深夜| 高清毛片免费观看视频网站| av在线蜜桃| 午夜激情福利司机影院| 国产精品福利在线免费观看| 成年女人毛片免费观看观看9| eeuss影院久久| av黄色大香蕉| 亚洲18禁久久av| 国产激情偷乱视频一区二区| 欧美一区二区亚洲| 色综合色国产| 亚洲最大成人手机在线| 又粗又爽又猛毛片免费看| 最新在线观看一区二区三区| 国产不卡一卡二| av在线老鸭窝| 九九热线精品视视频播放| 少妇熟女aⅴ在线视频| 97热精品久久久久久| 亚洲成人久久性| 97超级碰碰碰精品色视频在线观看| 国产一区二区在线观看日韩| 中文在线观看免费www的网站| 哪里可以看免费的av片| 亚洲国产日韩欧美精品在线观看| 欧美日韩国产亚洲二区| 91久久精品国产一区二区成人| 中文亚洲av片在线观看爽| 99热这里只有是精品在线观看| 久久久久久久精品吃奶| 91麻豆精品激情在线观看国产| 久久久午夜欧美精品| 国产亚洲精品久久久久久毛片| 老熟妇乱子伦视频在线观看| 男插女下体视频免费在线播放| 日日啪夜夜撸| 亚洲精品乱码久久久v下载方式| 琪琪午夜伦伦电影理论片6080| 色综合色国产| 国产成人aa在线观看| 在线观看一区二区三区| 欧美潮喷喷水| 69人妻影院| 一区二区三区高清视频在线| 22中文网久久字幕| 亚洲人成网站在线播放欧美日韩| 色综合婷婷激情| 老师上课跳d突然被开到最大视频| 51国产日韩欧美| 性色avwww在线观看| 麻豆国产97在线/欧美| 一个人看的www免费观看视频| 三级毛片av免费| 高清日韩中文字幕在线| 搡老熟女国产l中国老女人| 一本久久中文字幕| 免费看a级黄色片| 又黄又爽又刺激的免费视频.| 亚洲av第一区精品v没综合| 身体一侧抽搐| 国产精品久久电影中文字幕| 丰满乱子伦码专区| 99视频精品全部免费 在线| 欧美中文日本在线观看视频| 免费高清视频大片| 成年版毛片免费区| 在线看三级毛片| 欧美黑人欧美精品刺激| 亚洲国产精品合色在线| 乱系列少妇在线播放| 一进一出好大好爽视频| 日韩国内少妇激情av| 18禁黄网站禁片午夜丰满| 精品一区二区免费观看| 内射极品少妇av片p| 国产淫片久久久久久久久| 久久99热这里只有精品18| 久久久久免费精品人妻一区二区| 给我免费播放毛片高清在线观看| 欧美bdsm另类| 国产午夜精品论理片| 乱系列少妇在线播放| 97超视频在线观看视频| 看十八女毛片水多多多| 久久久国产成人精品二区| 亚洲国产欧美人成| 精品一区二区三区视频在线观看免费| 亚洲avbb在线观看| 国产精品自产拍在线观看55亚洲| 夜夜看夜夜爽夜夜摸| 在线播放国产精品三级| 日日干狠狠操夜夜爽| 别揉我奶头~嗯~啊~动态视频| 精品无人区乱码1区二区| 一本精品99久久精品77| 国产精品女同一区二区软件 | 国产一区二区三区在线臀色熟女| 在线观看舔阴道视频| 午夜福利在线观看吧| 日韩中字成人| 亚洲最大成人中文| 国产v大片淫在线免费观看| 在线观看美女被高潮喷水网站| 久久久精品欧美日韩精品| 91久久精品电影网| 99热精品在线国产| 久久国内精品自在自线图片| 日日撸夜夜添| 欧美日韩综合久久久久久 | 最近最新中文字幕大全电影3| 国产精品久久视频播放| 国产精品无大码| 尾随美女入室| 男人狂女人下面高潮的视频| 超碰av人人做人人爽久久| 午夜福利18| 久久精品国产亚洲av涩爱 | 小蜜桃在线观看免费完整版高清| 男人和女人高潮做爰伦理| 国产欧美日韩精品亚洲av| 国产免费一级a男人的天堂| 男女那种视频在线观看| 国内精品宾馆在线| 国产亚洲91精品色在线| 免费在线观看影片大全网站| 老女人水多毛片| 久久午夜亚洲精品久久| 国产成年人精品一区二区| 一区二区三区激情视频| av在线老鸭窝| 亚洲国产色片| 成年女人永久免费观看视频| 男人舔奶头视频| 特级一级黄色大片| 女的被弄到高潮叫床怎么办 | 又爽又黄a免费视频| 国产精品久久久久久精品电影| 亚洲精品在线观看二区| 成人国产麻豆网| 午夜福利在线观看免费完整高清在 | 亚洲欧美日韩卡通动漫| 99热网站在线观看| 麻豆久久精品国产亚洲av| 成人欧美大片| 在线看三级毛片| 久久久久久久精品吃奶| 丰满的人妻完整版| 亚洲中文字幕一区二区三区有码在线看| 成人综合一区亚洲| 欧美xxxx性猛交bbbb| 18禁在线播放成人免费| 午夜福利欧美成人| 一级黄色大片毛片| 午夜福利欧美成人| 22中文网久久字幕| 亚洲精品亚洲一区二区| 国产爱豆传媒在线观看| 久久久成人免费电影| 久久人人爽人人爽人人片va| 国产精品日韩av在线免费观看| 亚洲18禁久久av| 十八禁网站免费在线| 久久久久久国产a免费观看| 91久久精品国产一区二区成人| 亚洲熟妇中文字幕五十中出| 久久欧美精品欧美久久欧美| 搞女人的毛片| 国产高清三级在线| 麻豆精品久久久久久蜜桃| 丝袜美腿在线中文| 欧美精品国产亚洲| 精品久久久久久,| 网址你懂的国产日韩在线| 久久久精品大字幕| 成人鲁丝片一二三区免费| 嫩草影院入口| 波多野结衣高清作品| 日日啪夜夜撸| 97人妻精品一区二区三区麻豆| 成年女人永久免费观看视频| 色精品久久人妻99蜜桃| 狠狠狠狠99中文字幕| 国产私拍福利视频在线观看| 精品人妻偷拍中文字幕| 日本熟妇午夜| 欧美日韩精品成人综合77777| 国产午夜福利久久久久久| 欧美人与善性xxx| 国产精品一区二区免费欧美| 亚洲成a人片在线一区二区| 亚洲va在线va天堂va国产| 久久久久久大精品| 日韩欧美一区二区三区在线观看| 一区二区三区高清视频在线| 国产毛片a区久久久久| 国产av麻豆久久久久久久| 成熟少妇高潮喷水视频| 简卡轻食公司| av视频在线观看入口| 十八禁网站免费在线| 97人妻精品一区二区三区麻豆| 午夜福利在线观看免费完整高清在 | 老司机午夜福利在线观看视频| 日本一二三区视频观看| 18禁在线播放成人免费| 亚洲av美国av| 美女免费视频网站| 1000部很黄的大片| 少妇熟女aⅴ在线视频| 欧美最黄视频在线播放免费| 校园春色视频在线观看| 午夜激情欧美在线| 久久天躁狠狠躁夜夜2o2o| 有码 亚洲区| 在现免费观看毛片| 国产aⅴ精品一区二区三区波| 麻豆成人av在线观看| 亚洲最大成人中文| a级一级毛片免费在线观看| 欧洲精品卡2卡3卡4卡5卡区| 日韩 亚洲 欧美在线| 久久久久国内视频| 嫩草影院新地址| 国产一区二区亚洲精品在线观看| 99热网站在线观看| www日本黄色视频网| 欧美+亚洲+日韩+国产| 女同久久另类99精品国产91| 两性午夜刺激爽爽歪歪视频在线观看| a级一级毛片免费在线观看| 真实男女啪啪啪动态图| 欧美另类亚洲清纯唯美| x7x7x7水蜜桃| 麻豆成人午夜福利视频| 国产成人影院久久av| 亚洲欧美日韩卡通动漫| 久久精品国产鲁丝片午夜精品 | 国产亚洲精品久久久久久毛片| 少妇熟女aⅴ在线视频| 尾随美女入室| 成人精品一区二区免费| 亚洲在线观看片| 又爽又黄a免费视频| 中文字幕免费在线视频6| 国产单亲对白刺激| 亚洲人成伊人成综合网2020| 欧美色视频一区免费| 最后的刺客免费高清国语| 成人高潮视频无遮挡免费网站| 国产探花在线观看一区二区| 国产大屁股一区二区在线视频| 日本三级黄在线观看| 欧美日韩亚洲国产一区二区在线观看| 免费看av在线观看网站| 在线看三级毛片| 非洲黑人性xxxx精品又粗又长| 午夜爱爱视频在线播放| 亚洲欧美日韩高清专用| 在现免费观看毛片| 男人和女人高潮做爰伦理| 国产亚洲精品久久久com| 最后的刺客免费高清国语| 国产精品久久视频播放| 亚洲精品影视一区二区三区av| 夜夜爽天天搞| а√天堂www在线а√下载| 精品久久久久久久久av| 亚洲乱码一区二区免费版| 色在线成人网| 热99在线观看视频| 最后的刺客免费高清国语| 在线播放国产精品三级| 少妇高潮的动态图| 欧美日本亚洲视频在线播放| 亚洲aⅴ乱码一区二区在线播放| 国产美女午夜福利| 精品一区二区三区av网在线观看| 国产精品久久久久久亚洲av鲁大| av国产免费在线观看| 国产蜜桃级精品一区二区三区| 97超视频在线观看视频| 国产男人的电影天堂91| 精品久久久久久久久av| 国产精品精品国产色婷婷| 少妇裸体淫交视频免费看高清| 日本免费a在线| 少妇猛男粗大的猛烈进出视频 | 网址你懂的国产日韩在线| 波野结衣二区三区在线| 99热网站在线观看| 麻豆成人午夜福利视频| 国产免费又黄又爽又色| 国产精品久久久久久久电影| 久久久久精品久久久久真实原创| 国产综合精华液| 免费在线观看成人毛片| 亚洲中文av在线| 国产成人免费观看mmmm| 青春草国产在线视频| 久久久欧美国产精品| 纵有疾风起免费观看全集完整版| 99久久中文字幕三级久久日本| 亚洲内射少妇av| 综合色丁香网| 丰满少妇做爰视频| 国产淫片久久久久久久久| 韩国av在线不卡| 亚洲国产成人一精品久久久| 97在线人人人人妻| 联通29元200g的流量卡| 国产毛片在线视频| 三级经典国产精品| 国产成人精品久久久久久| 亚洲一区二区三区欧美精品| av国产免费在线观看| 99久久中文字幕三级久久日本| 多毛熟女@视频| 最近中文字幕高清免费大全6| av国产久精品久网站免费入址| 欧美成人一区二区免费高清观看| 噜噜噜噜噜久久久久久91| 国产精品精品国产色婷婷| 日韩成人av中文字幕在线观看| 国产精品一及| 狂野欧美激情性xxxx在线观看| 亚洲国产色片| 日韩,欧美,国产一区二区三区| 日本午夜av视频| 国产极品天堂在线| 2021少妇久久久久久久久久久| 91狼人影院| 国产男人的电影天堂91| 老熟女久久久| av免费观看日本| 一级黄片播放器| 天堂中文最新版在线下载| 精品少妇久久久久久888优播| 亚洲丝袜综合中文字幕| 久久精品国产自在天天线| 狂野欧美激情性bbbbbb| 欧美三级亚洲精品| 久久午夜福利片| 91精品一卡2卡3卡4卡| 日韩中字成人| 久久人人爽人人片av| 夜夜看夜夜爽夜夜摸| av卡一久久| 99国产精品免费福利视频| 亚洲国产欧美人成| 麻豆精品久久久久久蜜桃| 五月开心婷婷网| 国产精品人妻久久久影院| 亚洲精品第二区| 成年美女黄网站色视频大全免费 | 亚洲内射少妇av| 亚洲国产欧美在线一区| 乱系列少妇在线播放| 精品国产一区二区三区久久久樱花 | 人妻 亚洲 视频| 大陆偷拍与自拍| av播播在线观看一区| 最近最新中文字幕免费大全7| 视频中文字幕在线观看| 99久国产av精品国产电影| 婷婷色av中文字幕| 免费人成在线观看视频色| 性高湖久久久久久久久免费观看| 国内少妇人妻偷人精品xxx网站| 又大又黄又爽视频免费| 日韩强制内射视频| 亚洲欧美清纯卡通| 晚上一个人看的免费电影| 亚洲丝袜综合中文字幕| 黑丝袜美女国产一区| 久久99精品国语久久久| 中文精品一卡2卡3卡4更新| 99久久综合免费| 91精品国产九色| 久久亚洲国产成人精品v| 两个人的视频大全免费| 九九久久精品国产亚洲av麻豆| 一级毛片久久久久久久久女| 国产片特级美女逼逼视频| 国产在线一区二区三区精| 亚洲欧美日韩东京热| 丝袜喷水一区| 国产亚洲av片在线观看秒播厂| 国产欧美亚洲国产| 午夜福利在线在线| av卡一久久| 十八禁网站网址无遮挡 | 国产综合精华液| 99久久人妻综合| 国产精品一区二区三区四区免费观看| 蜜臀久久99精品久久宅男| 2018国产大陆天天弄谢| 亚洲内射少妇av| 亚洲精品国产成人久久av| 在线观看国产h片| 久久久久久久精品精品| 国精品久久久久久国模美| 日韩 亚洲 欧美在线| 久久久欧美国产精品| 又大又黄又爽视频免费| 最近最新中文字幕免费大全7| 亚州av有码| 亚洲欧洲日产国产| 国产成人91sexporn| 亚洲国产av新网站| 中国三级夫妇交换| 欧美精品国产亚洲| 狂野欧美激情性xxxx在线观看| 久久久久人妻精品一区果冻| 综合色丁香网| 91在线精品国自产拍蜜月| 观看美女的网站| 久久 成人 亚洲| 黄色视频在线播放观看不卡| 黄色欧美视频在线观看| 大话2 男鬼变身卡| 久久久成人免费电影| 美女内射精品一级片tv| 妹子高潮喷水视频| 国产又色又爽无遮挡免| 久久热精品热| 精品午夜福利在线看| 精品久久久久久久末码| 99久久人妻综合| 日本黄色日本黄色录像| 亚洲精品,欧美精品| 中国国产av一级| 最近2019中文字幕mv第一页| 涩涩av久久男人的天堂| 国产精品久久久久久久电影| 精华霜和精华液先用哪个| 18禁在线播放成人免费| 久热久热在线精品观看| 少妇的逼水好多| 永久免费av网站大全| 久久久精品94久久精品| 老司机影院成人| 免费黄频网站在线观看国产| 久久久久久久国产电影| 亚洲在久久综合| 国产黄频视频在线观看| av在线播放精品| 亚洲av不卡在线观看| 久久国产精品大桥未久av | 日韩欧美一区视频在线观看 | 人妻制服诱惑在线中文字幕| 国产成人精品婷婷| 亚洲国产精品国产精品| 久久久久久伊人网av| 在线亚洲精品国产二区图片欧美 | 老司机影院成人| 欧美最新免费一区二区三区| av免费观看日本| 搡老乐熟女国产| 国产一级毛片在线| 97超碰精品成人国产| 我的女老师完整版在线观看| 丝瓜视频免费看黄片| av在线观看视频网站免费| a 毛片基地| 三级经典国产精品| 国产高清国产精品国产三级 | 免费大片黄手机在线观看| 亚洲美女黄色视频免费看| 久久久久久久精品精品| 尤物成人国产欧美一区二区三区| 性高湖久久久久久久久免费观看| 日韩不卡一区二区三区视频在线| 一本久久精品| 国产日韩欧美在线精品| 两个人的视频大全免费| 下体分泌物呈黄色| 国产男女超爽视频在线观看| 亚洲成人手机| 国产毛片在线视频| 天堂8中文在线网| 蜜桃在线观看..| 成人18禁高潮啪啪吃奶动态图 | 国产黄色免费在线视频| 好男人视频免费观看在线| 成人亚洲欧美一区二区av| 国产av国产精品国产| 国产成人freesex在线| 中文精品一卡2卡3卡4更新| 日日啪夜夜爽| 国产午夜精品久久久久久一区二区三区| 毛片女人毛片| 国产一区二区三区综合在线观看 | 久久精品国产自在天天线| 亚洲欧美日韩另类电影网站 | 亚洲成人av在线免费| 亚洲av中文av极速乱| 2018国产大陆天天弄谢| 中文字幕制服av| 国产av精品麻豆| 一级二级三级毛片免费看| 极品教师在线视频| 国产在线免费精品| 国产精品成人在线| 久久久久久久精品精品| 成年av动漫网址| 亚洲av日韩在线播放| 欧美zozozo另类| 国产亚洲av片在线观看秒播厂| 国产亚洲5aaaaa淫片| 亚洲精品久久午夜乱码| 国产在线男女| 亚洲精品久久午夜乱码| 老女人水多毛片| 国产精品人妻久久久影院| 国产大屁股一区二区在线视频| 黑人猛操日本美女一级片| 久久精品国产亚洲av涩爱| 91在线精品国自产拍蜜月| 国产白丝娇喘喷水9色精品| 国产成人a∨麻豆精品| 久久精品国产亚洲av涩爱| 18禁在线播放成人免费| 在线精品无人区一区二区三 | 国产精品蜜桃在线观看| 国产黄色免费在线视频| 久久久久性生活片| 男女下面进入的视频免费午夜| 九九爱精品视频在线观看| 亚洲成人中文字幕在线播放| 中文在线观看免费www的网站| 美女国产视频在线观看| 免费人成在线观看视频色| 亚洲国产精品一区三区| 日韩电影二区| 色视频在线一区二区三区| 人妻一区二区av| 欧美成人一区二区免费高清观看| 国产高潮美女av| 亚洲精品日本国产第一区| av在线观看视频网站免费| 女性被躁到高潮视频|