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

水蒸氣核化及單體生長過程的分子動力學研究

2017-11-28 01:33:51馮靖伊張燕平
動力工程學報 2017年11期
關(guān)鍵詞:模型

馮靖伊, 張燕平, 高 偉

(華中科技大學 能源與動力工程學院, 武漢 430074)

水蒸氣核化及單體生長過程的分子動力學研究

馮靖伊, 張燕平, 高 偉

(華中科技大學 能源與動力工程學院, 武漢 430074)

通過分子技術(shù)模擬水蒸氣凝結(jié)過程,采用非平衡分子動力學(NEMD)方法模擬并分析了凝結(jié)相變中核化及單體生長過程.計算相變過程中的凝結(jié)速率,并研究了核化和單體生長初期液態(tài)水的直徑與凝結(jié)速率之間的關(guān)系,同時分析了模擬壓力對凝結(jié)速率的影響.結(jié)果表明:在相同過冷度下,凝結(jié)速率隨壓力的降低而減慢;在核化及單體生長過程中,相同的壓力和溫度條件下,凝結(jié)相變的傳熱傳質(zhì)效果隨著液體直徑的增長而改善.

水蒸氣; 分子動力學; 凝結(jié)速率; 液體尺寸; 傳熱傳質(zhì)

作為一種高效的傳熱方式,凝結(jié)換熱一直是眾多學者的重要研究方向.其中水蒸氣的凝結(jié)換熱更是被廣泛應(yīng)用于制冷[1]制熱、發(fā)電和航空航天等工程生產(chǎn)實踐中.幾十年來無數(shù)學者對凝結(jié)換熱進行了實驗和理論研究[2-4],但由于凝結(jié)相變本質(zhì)的復(fù)雜性,目前對其機理研究尚不全面.

近年來,隨著計算機性能的不斷進步,分子模擬技術(shù)在傳熱傳質(zhì)方面的研究也迅速發(fā)展[5].Ikeshoji等[6]提出了適用于周期性邊界條件的非平衡分子動力學算法來計算氬的導(dǎo)熱系數(shù),并通過該方法模擬氬蒸氣的凝結(jié)相變過程,發(fā)現(xiàn)分子間能量和動能的傳遞分別在液相和氣相的熱量傳遞中起主要作用.Krasnochtchekov等[7]采用分子動力學方法模擬了氣態(tài)納米金屬粒子的凝結(jié)相變過程,并在此基礎(chǔ)上分析了表面偏析對金屬粒子結(jié)構(gòu)的影響.Nagayama等[8]用分子動力學(MD)方法模擬了直鏈烷烴分子(如丁烷、辛烷等)在氣液界面上的凝結(jié)/蒸發(fā)行為,分析了相變的影響因素,并總結(jié)出凝結(jié)/蒸發(fā)系數(shù)的預(yù)測方法.Wang等[9]提出了修正過渡態(tài)理論,并采用分子動力學方法依據(jù)修正的過渡態(tài)理論模擬不同溫度下氬蒸氣的凝結(jié)系數(shù).Desbiens等[10]研究了水蒸氣在硅沸石材料中凝結(jié)吸附和液態(tài)水在疏水納米孔洞中的不均勻流動過程.傳熱發(fā)生時傳熱物質(zhì)形態(tài)和尺寸[11]的不同會導(dǎo)致?lián)Q熱效果差異較大,Diao等[12]通過構(gòu)造冷源、熱源以形成溫度梯度,模擬分析對碳納米管傳熱系數(shù)產(chǎn)生影響的因素,發(fā)現(xiàn)納米管封口的存在會大大降低傳熱性能.而水蒸氣凝結(jié)相變發(fā)生時傳熱物質(zhì)形態(tài)也會影響相變效果,如膜狀凝結(jié)的換熱效率遠低于珠狀凝結(jié),在珠狀凝結(jié)實驗數(shù)據(jù)[13]中也顯示,水蒸氣凝結(jié)時液體微團的尺寸對凝結(jié)效果也有影響,但是該宏觀實驗結(jié)果在微觀方面并沒有得到充分的模擬證實.

近年來針對凝結(jié)相變初期(核化及單體生長)的模擬也有一些研究,Mokshin等[14]通過使用粗粒度水分子模型模擬穩(wěn)態(tài)均勻蒸氣凝結(jié)的成核和單體生長過程.Diemand等[15]模擬了大規(guī)模的氬蒸氣均勻成核現(xiàn)象,在證明模擬結(jié)果與實驗值一致的基礎(chǔ)上提出了新的基于自由能函數(shù)的經(jīng)驗成核模型.Zipoli等[16]提出了一種用于模擬水蒸氣凝結(jié)核化的粗粒度模型,將粗粒度技術(shù)的適用范圍擴展到成核現(xiàn)象,通過對比發(fā)現(xiàn)該模型在得到準確結(jié)果的前提下能夠減少計算量.McGrath等[17]采用不同分子間相互作用勢來模擬氬的核化過程,發(fā)現(xiàn)作用勢之間臨界核形成的自由能差異不能僅根據(jù)勢的宏觀性質(zhì)差異來解釋,并對經(jīng)典成核理論的不準確進行了解釋.

雖然已有研究者對凝結(jié)核化和單體生長過程進行了研究,但這些研究還不夠充分,且對不同壓力下凝結(jié)相變影響因素的分析也不夠全面.筆者采用非平衡分子動力學(NEMD)方法,研究汽水換熱過程中水蒸氣在凝結(jié)過程中的核化及單體生長,模擬不同壓力下水蒸氣的凝結(jié)相變過程和傳熱傳質(zhì)現(xiàn)象,并研究了核化和單體生長過程中影響凝結(jié)換熱效果的相關(guān)因素(如液滴尺寸),為提高凝結(jié)換熱效果提供參考.

1 模擬體系建立

采用分子模擬技術(shù)模擬汽水混合換熱過程中水蒸氣的凝結(jié)換熱過程.模擬對象水分子選用CHARMM力場下的TIP3P[18]非剛體模型,水分子模型間的相互作用如下:

(1)

Uangle=Kangle(θ-θ0)2

(2)

Ubond=Kbond(r-r0)2

(3)

式(1)描述的是水分子間的作用勢能,其中右側(cè)第一項為長程靜電作用項,第二項為短程L-J靜電作用項.

整個體系初始狀態(tài)為充滿水分子的長方體模型,XYZ3個方向均采用周期性邊界條件.整個模擬流程圖見圖1.

圖1 模擬流程圖Fig.1 Simulation flow chart

模擬過程分為2個階段:第一階段為預(yù)處理階段,第二階段為凝結(jié)相變階段.第一階段模型的初始狀態(tài)為過熱水蒸氣.整個過程步長設(shè)為1 fs.為了節(jié)省計算時間采用勢能截斷方法對L-J勢能進行處理,截斷半徑取1 nm.采用Ewald加和法(Ewald Sums)處理庫侖力.分子初始速度按高斯分布取樣,求解采用Verlet velocity算法.首先對模型進行能量最小化處理,使盒子中所有分子能量最小化;然后模擬退火處理,控制體系溫度從470 K加熱到520 K然后冷卻至初始狀態(tài),處理總時間為200 ps.最后選用NVT系綜模擬200 ps,過程中使用Nose-Hoover控溫法將模擬溫度控制為470 K,使模型體系達到穩(wěn)態(tài).在此基礎(chǔ)上進行第二階段的凝結(jié)相變模擬.

第二階段采用非平衡態(tài)模擬方法.將模型中心置于點(0,0,0)處并按邊長分為100等份,其中邊長為L,每份長度為b,將立方體盒子的中心區(qū)域-6blt;xlt;6b,-6blt;ylt;6b部分設(shè)為冷源區(qū)域(符號C),對該區(qū)域分子降溫使水分子核化凝結(jié).將x-1/2Lxgt;-15b,x+1/2Lxlt;15b,y-1/2Lygt;-15b,y+1/2Lylt;15b區(qū)域設(shè)為熱源區(qū)域(符號H),加熱區(qū)域中的水分子,使該區(qū)域維持在過熱蒸汽狀態(tài).中間區(qū)域(符號M)為傳遞熱流區(qū)域并出現(xiàn)溫度梯度,區(qū)域劃分如圖2所示.在非平衡過程中整個體系選用NVE系綜,選用Langevin恒溫器控制冷源溫度在435 K,熱源溫度在475 K,模擬200萬步共2.0 ns.

圖2 冷源、熱源及中間區(qū)域劃分示意圖Fig.2 Division of cold source, heat source and middle region

2 結(jié)果與討論

2.1壓力對凝結(jié)相變的影響

在0.4 ns的平衡態(tài)模擬后開始模擬凝結(jié)相變,圖3給出了初始壓力為1.7 MPa的模擬體系發(fā)生凝結(jié)相變的過程.非平衡態(tài)模擬進行到1.1 ns時,水分子相互聚集形成密度約為920 kg/m3的液滴,由成核階段進入單體生長階段,1.45 ns后液化速度減慢,分子數(shù)基本達到穩(wěn)定.整個模擬過程中液相部分的形狀有輕微變化,難以維持標準的球形.這是由液態(tài)水體積較小,氣液界面波動的存在對液體形狀影響較為明顯造成的.

(a) 0.40 ns

(b) 0.55 ns

(c) 0.70 ns

(d) 0.85 ns

(f) 1.30 ns

(g) 1.45 ns

(h) 1.90 ns圖3 水蒸氣凝結(jié)相變過程Fig.3 Water vapor condensation process

通過模擬2.2 MPa、1.7 MPa和1.0 MPa 3個不同初始壓力下水蒸氣的凝結(jié)過程,并對發(fā)生凝結(jié)的水分子數(shù)目和模型壓力進行記錄,結(jié)果如圖4和圖5所示.3個不同初始壓力的模型包含相同的原子數(shù)目、冷源過冷度和冷熱源溫差.圖4為3個模型從非平衡態(tài)模擬開始到結(jié)束發(fā)生凝結(jié)相變的液相原子數(shù)隨時間的增長情況.圖5為各模型在非平衡態(tài)的壓力變化示意圖.從圖5可以看出,相同條件下,壓力降低,水蒸氣的凝結(jié)速率也隨之減慢.

圖6給出了初始壓力1.7 MPa,模型從非平衡態(tài)模擬開始,液態(tài)水的直徑及直徑增長速率隨時間的變化情況.從圖6可以看出,在液化過程開始到液滴直徑達到穩(wěn)定的過程中,液態(tài)水的直徑增長速率呈現(xiàn)先增大后減小的趨勢.模擬過程中,隨著水蒸氣凝結(jié)空間內(nèi)的壓力降低,直徑增長速率應(yīng)呈現(xiàn)減小的趨勢,而0~0.6 ns時液態(tài)水的直徑增長趨勢與該規(guī)律并不相符,表明模擬過程中其他因素影響了凝結(jié)相變的速率.根據(jù)模擬結(jié)果,推測該影響因素為液滴尺寸:在核化和單體生長初期,隨著凝結(jié)液相部分直徑的增長速率加快.

圖4 液相原子數(shù)變化Fig.4 Number of atoms in liquid phase

圖5 模型壓力變化Fig.5 Variation of simulation pressure

圖6 核化及單體生長過程中液體直徑的變化Fig.6 Variation of droplet size during nucleation and monomeric growth

2.2凝結(jié)相變的傳熱分析

以長12 nm、寬12 nm、高6 nm,初始壓力1.7 MPa的模型為例,記錄從非平衡態(tài)模擬階段開始到結(jié)束的整個過程中冷源和熱源的吸熱量、放熱量和模型內(nèi)分子的能量狀態(tài).一方面通過冷源、熱源的熱流量累積可得到凝結(jié)相變過程中的換熱量;另一方面通過記錄模型中分子的動能、勢能的變化情況也可以衡量凝結(jié)相變過程中的換熱量.從相變開始到結(jié)束整個過程中累積的換熱量如圖7所示.

由圖7可知,在0.5~0.8 ns內(nèi)平均壓力為0.921 MPa,該段時間內(nèi)液化了的水分子數(shù)為30,總換熱量為3 422.553 kJ/mol.圖8給出了系統(tǒng)中分子總能量的變化.從圖8可以看出,同樣在0.5~0.8 ns氣液界面的總能量的變化量(即總換熱量)為3 018.293 kJ/mol.2種記錄熱流方法的總換熱量值相差不大,證明了熱流記錄結(jié)果的合理性和準確性.

圖7 換熱量變化圖Fig.7 Changes of heat exchange

圖8 系統(tǒng)總能量變化Fig.8 Variation of total molecular energy in the system

為進一步分析模擬過程中的凝結(jié)換熱效果,需計算單位時間內(nèi)單位換熱面積液滴的傳熱傳質(zhì)情況.通過記錄各個時刻水分子的位置就可以得到各時刻水分子的空間分布以及在平面上的密度分布情況,從而得出液態(tài)水的位置及直徑,并求得換熱面積.在0.5~0.8 ns內(nèi),發(fā)生凝結(jié)的水分子數(shù)為30,經(jīng)過計算可知該統(tǒng)計時段氣液界面上平均凝結(jié)14.755個水分子/(nm2·ns),釋放潛熱130.95 kJ/mol.整個模型內(nèi)部各位置的溫度可通過計算原子的動能得到,具體公式如下:

Ek=dim/2NkT

(4)

式中:Ek為原子總動能,J;dim為模擬的維度,取為3;N為原子個數(shù);k為玻爾茲曼常數(shù),取1.38×10-23J/K;T為溫度,K.

根據(jù)記錄和分析所得熱流量總和Q、換熱面積A和換熱溫差(通過統(tǒng)計溫度分布可得)等條件,可以計算模擬過程中的平均凝結(jié)傳熱系數(shù),0.5~0.8 ns內(nèi)傳熱系數(shù)約為107W/(m2·K),可以推測出在核化和單體生長初期階段凝結(jié)相變的傳熱性能很好.

2.3尺寸對凝結(jié)速率的影響

由第2.1節(jié)分析得出,在凝結(jié)相變的核化和單體生長初期,隨著凝結(jié)液態(tài)水直徑的增長凝結(jié)換熱效果越好.針對該推測構(gòu)造邊長分別為12 nm、18 nm和24 nm,高為6 nm的3個長方體模型.這3個模型具有相同的初始壓力和溫度分布;冷源、熱源及中間區(qū)域的分布比例一致,且初始狀態(tài)的分子密度相同,經(jīng)過平衡態(tài)和非平衡態(tài)模擬使模型中心發(fā)生凝結(jié)相變形成液態(tài)水.模擬過程中監(jiān)測體系的壓力發(fā)生變化,在模型液態(tài)分子數(shù)趨于穩(wěn)定前的0.3 ns,3個模型平均壓力約為0.575 MPa,統(tǒng)計這0.3 ns內(nèi)液相部分的平均換熱面積,記錄該時段水分子凝結(jié)數(shù),得出3個模型的不同尺寸液體微團的水分子凝結(jié)速率,結(jié)果見表1.

表1 凝結(jié)速率與液體尺寸的關(guān)系

由表1可知,液體凝結(jié)速率隨液體直徑增長而加快.表明微觀條件下凝結(jié)相變的核化及單體生長階段凝結(jié)換熱效果與液體的尺寸有關(guān),并且凝結(jié)相變傳熱傳質(zhì)效果隨液態(tài)水直徑的增長而增強.Mokshin等[14]通過模擬不同溫度下穩(wěn)態(tài)均勻蒸氣凝結(jié),發(fā)現(xiàn)核化和單體生長過程中的液體微團直徑大小與時間呈指數(shù)關(guān)系,王愛麗[13]關(guān)于滴狀凝結(jié)的實驗也得出在溫度和壓力不變的情況下凝結(jié)相變初期的凝結(jié)速率有一個快速增長時期,與本文非平衡態(tài)模擬的結(jié)果一致.

為進一步驗證結(jié)論的正確性,針對3個模型在不同壓力(0.678 MPa、0.575 MPa和0.472 MPa)下相同時長內(nèi)的水分子凝結(jié)數(shù)進行統(tǒng)計,并換算成凝結(jié)相變釋放的潛熱量(見圖9).由圖9可知,相同壓力下隨著液滴直徑的增長凝結(jié)相變傳熱傳質(zhì)效果更優(yōu).

表1和圖9從不同的角度驗證了“在凝結(jié)相變的核化和單體生長初期階段,隨著凝結(jié)液態(tài)水直徑的增長凝結(jié)換熱效果更好” 的推測(在不同壓力下)的合理性,這是目前宏觀實驗和模擬中還沒得到過的結(jié)果.

圖9 不同壓力下凝結(jié)釋放的潛熱量與液滴直徑的關(guān)系Fig.9 Latent heat of condensation vs. droplet size at different pressures

3 結(jié) 論

(1) 相同過冷度條件下,隨著壓力的升高凝結(jié)速率越來越快,表明提高壓力可以促進凝結(jié)相變的發(fā)生.

(2) 在凝結(jié)相變的核化和單體生長初期階段,凝結(jié)相變的傳熱系數(shù)約為107W/(m2·K) ,表明凝結(jié)的這段過程具有良好的傳熱傳質(zhì)性能.

(3) 在凝結(jié)相變的核化和單體生長初期,隨著液滴直徑的增長凝結(jié)速率和凝結(jié)相變傳熱傳質(zhì)效果更好,這也是目前宏觀實驗和模擬還未觀察到的結(jié)果.

[2] 胡申華, 嚴俊杰, 王進仕. 非等溫表面Marangoni凝結(jié)特性的研究[J].動力工程學報, 2012, 32(1): 27-30.

HU Shenhua, YAN Junjie, WANG Jinshi. Study of Marangoni condensation characteristics on non-isothermal surface[J].JournalofChineseSocietyofPowerEngineering,2012, 32(1): 27-30.

[3] DAVIES J F, MILES R E H, HADDRELL A E, et al. Influence of organic films on the evaporation and condensation of water in aerosol[J].ProceedingsoftheNationalAcademyofSciencesoftheUnitedStatesofAmerica, 2013, 110(22): 8807-8812.

[4] MILJKOVIC N, WANG E N. Condensation heat transfer on superhydrophobic surfaces[J].MRSBulletin, 2013, 38(5): 397-406.

[5] ACEVEDO O, JORGENSEN W L. Quantum and molecular mechanical Monte Carlo techniques for modeling condensed-phase reactions[J].WileyInterdisciplinaryReviews:ComputationalMolecularScience, 2014, 4(5): 422-435.

[6] IKESHOJI T, HAFSKJOLD B. Non-equilibrium molecular dynamics calculation of heat conduction in liquid and through liquid-gas interface[J].MolecularPhysics, 1994, 81(2): 251-261.

[7] KRASNOCHTCHEKOV P, ALBE K, AVERBACK R S. Simulations of the inert gas condensation processes[J].ZeitschriftfürMetallkunde, 2003, 94(10): 1098-1105.

[8] NAGAYAMA G, TAKEMATSU M, MIZUGUCHI H, et al. Molecular dynamics study on condensation/evaporation coefficients of chain molecules at liquid-vapor interface[J].TheJournalofChemicalPhysics, 2015, 143(1): 014706.

[9] WANG Z J, MIN C, GUO Z Y. Modified transition state theory for evaporation and condensation[J].ChinesePhysicsLetters, 2002, 19(4): 537-539.

[10] DESBIENS N, BOUTIN A, DEMACHY I. Water condensation in hydrophobic silicalite-1 zeolite: a molecular simulation study[J].TheJournalofPhysicalChemistryB, 2005, 109(50): 24071-24076.

[11] QI B J, WEI J J, ZHANG L, et al. A fractal dropwise condensation heat transfer model including the effects of contact angle and drop size distribution[J].InternationalJournalofHeatandMassTransfer, 2015, 83: 259-272.

[12] DIAO J K, SRIVASTAVA D, MENON M. Molecular dynamics simulations of carbon nanotube/silicon interfacial thermal conductance[J].TheJournalofChemicalPhysics, 2008, 128(16): 164708.

[13] 王愛麗. 滴狀冷凝液滴微觀特征及傳熱機制[D]. 大連: 大連理工大學, 2010.

[14] MOKSHIN A V, GALIMZYANOV B N. Steady-state homogeneous nucleation and growth of water droplets: extended numerical treatment[J].TheJournalofPhysicalChemistryB, 2012, 116(39): 11959-11967.

[15] DIEMAND J, ANGéLIL R, TANAKA K K, et al. Large scale molecular dynamics simulations of homogeneous nucleation[J].TheJournalofChemicalPhysics, 2013, 139(7): 074309.

[16] ZIPOLI F, LAINO T, STOLZ S, et al. Improved coarse-grained model for molecular-dynamics simulations of water nucleation[J].TheJournalofChemicalPhysics, 2013, 139(9): 094501.

[17] McGRATH M J, GHOGOMU J N, TSONA N T, et al. Vapor-liquid nucleation of argon: exploration of various intermolecular potentials[J].TheJournalofChemicalPhysics, 2010, 133(8): 084106.

[18] MARK P, NILSSON L. Structure and dynamics of the TIP3P, SPC, and SPC/E water models at 298 K[J].TheJournalofPhysicalChemistryA, 2001, 105(43): 9954-9960.

MolecularDynamicsInvestigationonNucleationandMonomericGrowthofWaterVapor

FENGJingyi,ZHANGYanping,GAOWei

(School of Energy and Power Engineering, Huazhong University of Science and Technology, Wuhan 430074, China)

The process of water vapor condensation was simulated by molecular techniques, while the nucleation and monomeric growth of water vapor in phase change period were simulated and analyzed using non-equilibrium molecular dynamics (NEMD) method. A calculation was conducted on the steam condensation rate in the process of phase change, with a study on the effects of droplet size on the condensation rate in the initial period of nucleation and monomeric growth, and with an analysis on the influence of simulation pressure on the condensation performance. Results show that in the case of same degree of supercooling, the condensation rate decreases with the reduction of pressure. During the process of nucleation and monomeric growth, both the heat and mass transfer effectiveness are relatively good at the same pressure and temperature, which increase with the rise of droplet sizes.

water vapor; molecular dynamics; condensation rate; droplet size; heat and mass transfer

2016-11-07

2017-01-05

國家重點基礎(chǔ)研究發(fā)展計劃(973計劃)資助項目(2015CB251504)

馮靖伊(1991-),女,吉林長春人,碩士研究生,研究方向為分子動力學.

張燕平(通信作者),女,副教授,博士,電話(Tel.):15927592300;E-mail:zyp2817@hust.edu.cn.

1674-7607(2017)11-0912-06

TK124

A

470.10

猜你喜歡
模型
一半模型
一種去中心化的域名服務(wù)本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數(shù)模型及應(yīng)用
p150Glued在帕金森病模型中的表達及分布
函數(shù)模型及應(yīng)用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产农村精品一级毛片视频| 国产一区二区三区精品久久呦| 婷婷激情亚洲| 国产成+人+综合+亚洲欧美| 亚洲无码高清免费视频亚洲| 欧美亚洲国产日韩电影在线| 中文字幕永久视频| 激情综合图区| 全裸无码专区| 国产成年女人特黄特色毛片免| 麻豆国产精品视频| 欧美精品啪啪| 九色免费视频| 香蕉久人久人青草青草| 久久伊伊香蕉综合精品| 99性视频| 亚洲欧美极品| 亚洲欧美综合在线观看| 久久精品人妻中文视频| 人人91人人澡人人妻人人爽| 精品欧美一区二区三区在线| 97精品久久久大香线焦| 欧美成人精品高清在线下载| 激情六月丁香婷婷四房播| 亚洲人成网站色7799在线播放| 欧美日一级片| 欧美日在线观看| 黄色网页在线播放| 欧美亚洲国产精品久久蜜芽| 亚洲精品国产成人7777| 欧美亚洲国产精品久久蜜芽| 狠狠躁天天躁夜夜躁婷婷| 91精品专区| 亚洲欧美日韩成人高清在线一区| 亚洲性日韩精品一区二区| 国产麻豆aⅴ精品无码| 亚洲黄色网站视频| 六月婷婷激情综合| 成人蜜桃网| 欧美人在线一区二区三区| 女高中生自慰污污网站| 天天摸天天操免费播放小视频| 婷婷亚洲视频| 麻豆精品在线视频| 国产91色| 久热精品免费| 午夜福利免费视频| 男人的天堂久久精品激情| 久久精品视频一| 国产精品第5页| 色屁屁一区二区三区视频国产| 欧美a网站| 91毛片网| 久久综合国产乱子免费| 国产成人资源| 国产美女在线观看| 亚洲无码一区在线观看| 国产欧美中文字幕| 欧美亚洲日韩中文| 动漫精品啪啪一区二区三区| 欧美精品成人| 亚洲欧美日本国产综合在线| 欧美伦理一区| 91在线播放国产| 一区二区无码在线视频| 欧美a级完整在线观看| 巨熟乳波霸若妻中文观看免费 | 激情综合网址| 久久精品无码中文字幕| 国产三级a| 精品成人一区二区三区电影| 亚洲无码高清视频在线观看| 亚洲美女视频一区| 免费在线看黄网址| 国产精品9| 无码福利视频| 福利片91| 日本亚洲国产一区二区三区| 一区二区理伦视频| 亚洲天堂区| 中文字幕色站| 午夜视频免费一区二区在线看|