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

    特征參數(shù)對(duì)DSMC方法模擬聲場(chǎng)中顆粒碰撞的影響

    2016-12-05 05:49:38郭孝武凡鳳仙胡曉紅蘇明旭
    關(guān)鍵詞:方法

    郭孝武, 凡鳳仙, 胡曉紅, 蘇明旭

    (1.上海理工大學(xué) 能源與動(dòng)力工程學(xué)院,上海 200093;2.上海理工大學(xué) 上海市動(dòng)力工程多相流動(dòng)與傳熱重點(diǎn)實(shí)驗(yàn)室,上海 200093)

    ?

    特征參數(shù)對(duì)DSMC方法模擬聲場(chǎng)中顆粒碰撞的影響

    郭孝武1,2, 凡鳳仙1,2, 胡曉紅1,2, 蘇明旭1,2

    (1.上海理工大學(xué) 能源與動(dòng)力工程學(xué)院,上海 200093;2.上海理工大學(xué) 上海市動(dòng)力工程多相流動(dòng)與傳熱重點(diǎn)實(shí)驗(yàn)室,上海 200093)

    直接從顆粒受力和運(yùn)動(dòng)出發(fā),基于直接模擬蒙特卡洛(DSMC)方法建立顆粒碰撞模型,模擬聲場(chǎng)中顆粒的碰撞過(guò)程,通過(guò)改變模擬條件,探討DSMC方法中的特征參數(shù)(取樣顆粒的數(shù)目權(quán)重、網(wǎng)格數(shù)目、時(shí)間步長(zhǎng))對(duì)顆粒碰撞率和計(jì)算時(shí)間的影響.結(jié)果表明:頻率越大,顆粒容積份額變化越迅速,顆??臻g分布越不均勻;數(shù)目權(quán)重的增加對(duì)顆粒碰撞率影響較小,而對(duì)計(jì)算時(shí)間的影響顯著,使其迅速減少;網(wǎng)格數(shù)目增加,碰撞率降低,計(jì)算時(shí)間則先迅速降低,隨后在低頻時(shí)基本不變,高頻時(shí)略有上升.研究還發(fā)現(xiàn):隨著碰撞時(shí)間步長(zhǎng)的增加,低頻聲場(chǎng)中碰撞率單調(diào)增加,高頻聲場(chǎng)中碰撞率先增加而后出現(xiàn)波動(dòng);碰撞時(shí)間步長(zhǎng)的增加將引起計(jì)算時(shí)間減少,減少量在碰撞步長(zhǎng)較小時(shí)最為明顯.

    聲場(chǎng); 顆粒; 碰撞; 直接模擬蒙特卡洛方法; 特征參數(shù)

    目前,用于聲凝并的數(shù)值模擬方法可歸納為區(qū)域算法、矩量法和直接模擬蒙特卡洛(direct simulation Monte Carlo,DSMC)方法[8].其中,DSMC方法立足于氣體分子運(yùn)動(dòng)論,將所計(jì)算的實(shí)際顆粒場(chǎng)用取樣顆粒場(chǎng)進(jìn)行置換,對(duì)每一個(gè)取樣顆粒賦予一個(gè)數(shù)目權(quán)重(即取樣顆粒所代表的真實(shí)顆粒數(shù)目),跟蹤取樣顆粒的運(yùn)動(dòng)軌跡,通過(guò)概率的方法判斷碰撞是否發(fā)生,并對(duì)碰撞后的顆粒凝并事件進(jìn)行處理,該方法可以方便地分析多分散顆粒凝并的動(dòng)態(tài)過(guò)程以及粒徑演變規(guī)律[9-12].2006年,Sheng等[9]利用常體積Monte Carlo方法在初始顆粒及團(tuán)聚體均為球形的前提下,研究了行波聲場(chǎng)中液態(tài)PM2.5的聲凝并,發(fā)現(xiàn)數(shù)值模擬結(jié)果在一些情況下與實(shí)驗(yàn)吻合較好,而當(dāng)顆粒初始粒徑分布發(fā)生改變時(shí),則與實(shí)驗(yàn)存在明顯差異;2007年,他們引入分形維數(shù)以描述顆粒團(tuán)聚體的形狀和結(jié)構(gòu),將上述聲凝并模型推廣到固體顆粒的聲凝并中[10].與Sheng等[9-10]在顆粒通用動(dòng)力學(xué)方程和凝并核函數(shù)的基礎(chǔ)上建立聲凝并模型不同,凡鳳仙等[11-13]從顆粒運(yùn)動(dòng)方程出發(fā),利用DSMC方法研究PM2.5在行波聲場(chǎng)中的碰撞特性[11]、在駐波聲場(chǎng)中的碰撞和凝并規(guī)律[12-13].顆粒碰撞是其發(fā)生凝并的基礎(chǔ)和前提,然而,在利用DSMC方法處理顆粒碰撞時(shí),3個(gè)特征參數(shù)(時(shí)間步長(zhǎng)、取樣顆粒的數(shù)目權(quán)重、計(jì)算區(qū)域內(nèi)的網(wǎng)格數(shù)目)需要人為設(shè)定,這些特征參數(shù)的取值將影響到計(jì)算精度[14-16].時(shí)間步長(zhǎng)與取樣顆粒的數(shù)目權(quán)重越小,計(jì)算越精確,然而隨之出現(xiàn)的是計(jì)算代價(jià)的增加;若在一個(gè)網(wǎng)格范圍內(nèi)判斷顆粒碰撞的發(fā)生,則計(jì)算區(qū)域內(nèi)的網(wǎng)格數(shù)目也將對(duì)計(jì)算結(jié)果和計(jì)算代價(jià)帶來(lái)影響.因此,本文基于DSMC方法,在不同特征參數(shù)條件下對(duì)聲場(chǎng)中顆粒的碰撞過(guò)程進(jìn)行數(shù)值模擬,確定特征參數(shù)對(duì)顆粒碰撞率和計(jì)算時(shí)間的影響,為利用DSMC方法準(zhǔn)確、高效地研究PM10,PM2.5的聲凝并提供參考.

    1 聲場(chǎng)中顆粒動(dòng)力學(xué)模型

    為利用DSMC方法研究聲場(chǎng)中顆粒的碰撞過(guò)程,作出如下簡(jiǎn)化假設(shè):

    a. 聲場(chǎng)為一維平面駐波聲場(chǎng),聲波波動(dòng)方向?yàn)樗椒较?

    b. 由于煙氣的物性參數(shù)和空氣類(lèi)似,且其壓力不太高,將氣體介質(zhì)視為理想空氣;

    c. 顆粒為剛性球體,忽略顆粒的轉(zhuǎn)動(dòng);

    d. 由于聲輻射壓力的作用效果極其微弱[17],不考慮聲輻射壓力;

    e. 顆粒碰撞在各網(wǎng)格內(nèi)進(jìn)行,且為二元碰撞,為著重考慮顆粒的碰撞,暫不考慮其凝并,認(rèn)為顆粒碰撞后彈開(kāi);

    f. 在顆粒碰撞過(guò)程中,顆粒發(fā)生滑移時(shí)所受摩擦力遵守庫(kù)侖摩擦定律.

    1.1 聲波波動(dòng)方程

    由Navier-Stokes方程,可推導(dǎo)出無(wú)旋、無(wú)粘流體中x向駐波聲場(chǎng)的波動(dòng)方程為

    (1)

    式中:ufx為聲波引起的流體介質(zhì)振動(dòng)速度;x為位置坐標(biāo);t為時(shí)間;ua為速度振幅;k為波數(shù),k=ω/c,c為聲速;ω=2πf,f為聲場(chǎng)頻率.

    習(xí)慣上常用聲壓級(jí)和頻率描述聲場(chǎng),聲壓級(jí)與速度振幅ua的關(guān)系為

    (2)

    式中:L為聲壓級(jí),dB;Pr為參考聲壓,Pr=210-5Pa.

    1.2 顆粒運(yùn)動(dòng)方程

    設(shè)重力方向?yàn)閦向,與x和z垂直的方向?yàn)閥向,流體攜帶顆粒沿y向以一定速度通過(guò)聲場(chǎng)空間.根據(jù)牛頓第二定律,顆粒運(yùn)動(dòng)方程可寫(xiě)為

    其中:式(3a)等號(hào)右邊4項(xiàng)依次為Stokes力的x向分量、壓力梯度力、虛擬質(zhì)量力、Basset力;式(3b)等號(hào)右邊為Stokes力的y向分量;式(3c)等號(hào)右邊3項(xiàng)分別為Stokes力的z向分量、重力、浮力.式中:mp為顆粒質(zhì)量;mf為與顆粒等體積的流體質(zhì)量;ufy和ufz分別為流體速度的y,z向分量;upy與upz分別為顆粒速度的y,z向分量;ρp為顆粒密度;g為重力加速度;t′為時(shí)間變量;Cc為Cunningham修正系數(shù),其表達(dá)式為[13]

    式中,Kn為Knudsen數(shù),Kn=2lm/dp,lm為氣體分子平均自由程,m.

    1.3 描述顆粒碰撞的DSMC方法

    對(duì)含有大量PM10,PM2.5的氣固懸浮體系施加駐波聲場(chǎng),顆粒間的碰撞不可避免.為了準(zhǔn)確地描述聲場(chǎng)中顆粒的行為規(guī)律,理想的方法是跟蹤到每一個(gè)顆粒,通過(guò)顆粒運(yùn)動(dòng)軌道來(lái)判斷顆粒碰撞,但當(dāng)顆粒量很大時(shí),這種方法的計(jì)算量將是非常驚人的.由于顆粒數(shù)目龐大,利用顆粒軌道判斷碰撞的方法難以勝任.解決這一問(wèn)題的一個(gè)行之有效的途徑是采用DSMC方法.由于DSMC方法處理的是遠(yuǎn)低于真實(shí)顆粒數(shù)目的取樣顆粒,從而大大降低了計(jì)算量,特別適于在較大規(guī)模顆粒量和較大區(qū)域范圍內(nèi)對(duì)小粒徑顆粒、較稀的氣固兩相流動(dòng)進(jìn)行模擬研究.在DSMC方法的前提下,不同研究者的處理和計(jì)算方法不盡相同,較常用的是將計(jì)算區(qū)域劃分成若干個(gè)網(wǎng)格,在同一網(wǎng)格內(nèi)判斷顆粒之間是否發(fā)生碰撞,本文也采用該方法.

    (5)

    取樣顆粒i與同一網(wǎng)格中其他顆粒的總碰撞概率Pi可表示為

    式中,N為取樣顆粒i所在網(wǎng)格的取樣顆??倲?shù).

    采用修正的Nanbu方法[19]判斷顆粒碰撞的發(fā)生,即對(duì)于任意取樣顆粒i,在滿(mǎn)足Pi<1的前提下,利用[0,1)區(qū)間上均勻分布的隨機(jī)數(shù)R,選擇同一網(wǎng)格內(nèi)的候選被碰取樣顆粒j為

    j=int[RN]+1

    (7)

    式中,int[RN]表示RN的整數(shù)部分.如果R滿(mǎn)足

    (8)

    則認(rèn)為取樣顆粒i與j發(fā)生碰撞.此時(shí),保持顆粒位置不變,顆粒碰撞后速度根據(jù)動(dòng)量定理計(jì)算[19],即

    (9)

    (10)

    (11)

    (12)

    (13)

    (14)

    (15)

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

    2.1 邊界條件與計(jì)算參數(shù)

    取λ2k×1 mm×1 mm的三維空間為計(jì)算區(qū)域,其中λ2k表示頻率f=2 kHz時(shí)對(duì)應(yīng)的聲波波長(zhǎng),沿x向?qū)⒂?jì)算區(qū)域劃分為若干個(gè)大小相等的網(wǎng)格.由聲波的周期性和本文采用的聲波頻率(f=2,4,6 kHz)可知,對(duì)x,y,z這3個(gè)方向的邊界均可采用周期性邊界條件[13]進(jìn)行處理,對(duì)各個(gè)取樣顆粒采用相同的數(shù)目權(quán)重,主要計(jì)算參數(shù)見(jiàn)表1.其中:T為氣體溫度;p為氣體靜態(tài)壓力;np為顆粒數(shù)目濃度;Δt為求解顆粒運(yùn)動(dòng)的時(shí)間步長(zhǎng).本文通過(guò)試算確定Δt,表1中給出的Δt滿(mǎn)足3種計(jì)算頻率下顆粒運(yùn)動(dòng)的計(jì)算精度要求.

    表1 數(shù)值計(jì)算參數(shù)

    2.2 顆粒運(yùn)動(dòng)方程的求解

    采用定步長(zhǎng)四階Runge-Kutta算法對(duì)式(3)進(jìn)行求解.其中,Basset力與顆粒經(jīng)歷的運(yùn)動(dòng)過(guò)程有關(guān).將流體速度和顆粒速度基于經(jīng)歷的時(shí)間步離散,可實(shí)現(xiàn)Basset力的數(shù)值求解.例如,在時(shí)間區(qū)間0~t3離散為0~t1,t1~t2,t2~t3的情況下,如果這3個(gè)時(shí)間區(qū)間對(duì)應(yīng)的流體與顆粒速度增量分別為Δufx0,Δufx1,Δufx2與Δupx0,Δupx1,Δupx2,則t3時(shí)刻顆粒受到的Basset力FB可按下式求解[20]:

    (16)

    顆粒運(yùn)動(dòng)軌跡可由各時(shí)間步長(zhǎng)內(nèi)顆粒的位移疊加得到,根據(jù)二階隱式Adams插值算法,顆粒的位移Sp可表示為

    式中,up為顆粒速度.

    2.3 數(shù)值計(jì)算流程

    圖1 計(jì)算流程

    3 結(jié)果與討論

    3.1 聲場(chǎng)中顆粒容積份額的演變

    圖2 不同頻率的聲場(chǎng)中顆粒容積份額的演變

    3.2 特征參數(shù)對(duì)碰撞率和計(jì)算時(shí)間的影響

    3.2.1 數(shù)目權(quán)重對(duì)碰撞率和計(jì)算時(shí)間的影響

    3.2.2 網(wǎng)格數(shù)目對(duì)碰撞率和計(jì)算時(shí)間的影響

    圖3 數(shù)目權(quán)重對(duì)碰撞率和計(jì)算時(shí)間的影響

    圖4 網(wǎng)格數(shù)目對(duì)碰撞率和計(jì)算時(shí)間的影響

    3.2.3 碰撞時(shí)間步長(zhǎng)對(duì)碰撞率和計(jì)算時(shí)間的影響

    圖5 碰撞時(shí)間步長(zhǎng)對(duì)碰撞率和計(jì)算時(shí)間的影響

    4 結(jié) 論

    綜合考慮顆粒受到的Stokes力、壓力梯度力、虛擬質(zhì)量力、Basset力、重力和浮力,建立顆粒運(yùn)動(dòng)模型,基于DSMC方法模擬聲場(chǎng)中顆粒的碰撞過(guò)程,獲得顆粒碰撞率和計(jì)算時(shí)間受DSMC方法中特征參數(shù)(時(shí)間步長(zhǎng)、數(shù)目權(quán)重、網(wǎng)格數(shù)目)的影響規(guī)律,得到以下結(jié)論:

    a. 相同聲場(chǎng)作用時(shí)間下,頻率越大,顆粒容積份額變化越迅速,顆粒空間分布也越不均勻;聲場(chǎng)頻率對(duì)碰撞率有顯著影響,而對(duì)計(jì)算時(shí)間的影響很小.

    b. 取樣顆粒數(shù)目權(quán)重增加,顆粒碰撞率的變化很小,計(jì)算時(shí)間則迅速降低;網(wǎng)格數(shù)目增加,碰撞率降低,計(jì)算時(shí)間先迅速降低,而后低頻時(shí)基本不變,高頻時(shí)略有上升.

    c. 隨著碰撞時(shí)間步長(zhǎng)的增加,低頻聲場(chǎng)中顆粒碰撞率單調(diào)增加,高頻聲場(chǎng)中顆粒碰撞率先增加而后出現(xiàn)波動(dòng);計(jì)算時(shí)間減少,減少量在碰撞步長(zhǎng)較小時(shí)最為明顯.

    [1] ZHANG L W,CHEN X,XUE X D,et al.Long-term exposure to high particulate matter pollution and cardiovascular mortality:a 12-year cohort study in four cities in northern China[J].Environment International,2014,62:41-47.

    [2] WANG T J,JIANG F,DENG J J,et al.Urban air quality and regional haze weather forecast for Yangtze River Delta region[J].Atmospheric Environment,2012,58:70-83.

    [3] WEN C,XU M H,YU D X,et al.PM10formation during the combustion of N2-char and CO2-char of Chinese coals[J].Proceedings of the Combustion Institute,2013,34(2):2383-2392.

    [4] EHRLICH C,NOLL G,KALKOFF W D,et al.PM10,PM2.5and PM1.0—Emissions from industrial plants—results from measurement programmes in Germany[J].Atmospheric Environment,2007,41(29):6236-6254.

    [5] WANG J,HU Z M,CHEN Y Y,et al.Contamination characteristics and possible sources of PM10and PM2.5in different functional areas of Shanghai,China[J].Atmospheric Environment,2013,68:221-229.

    [6] 楊林軍.燃燒源細(xì)顆粒物污染控制技術(shù)[M].北京:化學(xué)工業(yè)出版社,2011.

    [7] 胡惠敏,李瑞陽(yáng),蔡萌,等.聲波與其他方法聯(lián)合作用脫除細(xì)顆粒物的研究進(jìn)展[J].上海理工大學(xué)學(xué)報(bào),2016,38(1):13-18.

    [8] 張明俊,凡鳳仙.細(xì)顆粒物的聲凝并數(shù)值模擬研究進(jìn)展[J].化工進(jìn)展,2012,31(8):1671-1676.

    [9] SHENG C D,SHEN X L.Modelling of acoustic agglomeration processes using the direct simulation Monte Carlo method[J].Journal of Aerosol Science,2006,37(1):16-36.

    [10] SHENG C D,SHEN X L.Simulation of acoustic agglomeration processes of poly-disperse solid particles[J].Aerosol Science and Technology,2007,41(1):1-13.

    [11] 凡鳳仙,袁竹林.外加聲場(chǎng)對(duì)增加PM2.5碰撞幾率的數(shù)值模擬研究[J].中國(guó)電機(jī)工程學(xué)報(bào),2006,26(11):12-16.

    [12] 凡鳳仙,袁竹林,趙兵,等.駐波聲場(chǎng)中細(xì)微顆粒凝并的數(shù)值模擬[J].燃燒科學(xué)與技術(shù),2008,14(3):253-258.

    [13] FAN F X,YANG X F,KIM C N.Direct simulation of inhalable particle motion and collision in a standing wave field[J].Journal of Mechanical Science and Technology,2013,27(6):1707-1712.

    [14] WEI J M.A fast Monte Carlo method based on an acceptance-rejection scheme for particle coagulation[J].Aerosol and Air Quality Research,2013,13:1273-1281.

    [15] HE Y X,ZHAO H B,WANG H M,et al.Differentially weighted direct simulation Monte Carlo method for particle collision in gas-solid flows[J].Particuology,2015,21:135-145.

    [16] KARCHANI A,MYONG R S.Convergence analysis of the direct simulation Monte Carlo based on the physical laws of conservation[J].Computers & Fluids,2015,115:98-114.

    [17] 宋曉通,凡鳳仙.駐波聲場(chǎng)中可吸入顆粒物漂移的影響因素分析[J].熱能動(dòng)力工程,2016,31(1):81-86.

    [18] CLECKLER J,ELGHOBASHI S,LIU F.On the motion of inertial particles by sound waves[J].Physics of Fluids,2012,24(3):033301.

    [19] TSUJI Y,TANAKA T,YONEMURA S.Cluster patterns in circulating fluidized beds predicted by numerical simulation(discrete particle model versus two-fluid model)[J].Powder Technology,1998,95(3):254-264.

    [20] ROSTAMI M,ARDESHIR A,AHMADI G,et al.On the effect of gravitational and hydrodynamic forces on particle motion in a quiescent fluid at high particle Reynolds numbers[J].Canadian Journal of Physics,2008,86(6):791-799.

    (編輯:丁紅藝)

    Influence of Characteristic Parameters in the Simulation of Acoustic Particle Collision Using DSMC Method

    GUO Xiaowu1,2, FAN Fengxian1,2, HU Xiaohong1,2, SU Mingxu1,2

    (1.School of Energy and Power Engineering,University of Shanghai for Science and Technology,Shanghai 200093,China; 2.Shanghai Key Laboratory of Multiphase Flow and Heat Transfer in Power Engineering,University of Shanghai for Science and Technology,Shanghai 200093,China)

    An inter-particle collision model was established based on the direct simulation Monte Carlo (DSMC) method as well as the force analysis and the motion equation.The collision process of the particles in acoustic field was examined.Through changing the numerical simulation conditions,the influences of the characteristic parameters,such as the weight value of the number of sampling particles,the cell number and the time-step size used in DSMC method,on the particle collision rate and computational time were discussed.The results show that the higher the acoustic frequency is,the more rapidly the particle volume fraction changes and the more unevenly the particles are distributed in space.The increase of the weight value of the number has a small effect on the collision rate,but it significantly reduces the computational time.As the cell number increases,the collision rate decreases,whereas the computational time decreases rapidly at first and then almost keeps constant at lower acoustic frequency and increases slightly at higher acoustic frequency.It is also found that as the time-step size for determining the inter-particle collision increases,the collision rate increases monotonously at lower frequency,while it increases at first and then fluctuates at higher frequency.The increase of the time-step size results in the decrease of the computational time and the decrease is more obvious with a smaller time-step size.

    acoustic field; particle; collision; direct simulation Monte Carlo method; characteristic parameter

    1007-6735(2016)05-0419-08

    10.13255/j.cnki.jusst.2016.05.003

    2016-01-25

    國(guó)家自然科學(xué)基金資助項(xiàng)目(51206113,51176128,51576130);上海市科委科研計(jì)劃資助項(xiàng)目(13DZ2260900)

    郭孝武(1988-),男,碩士研究生.研究方向:氣固兩相流數(shù)值模擬.E-mail:ross_2015@163.com

    凡鳳仙(1982-),女,副教授.研究方向:燃燒源污染物排放控制、氣固流動(dòng)與傳遞.

    E-mail:fanfengxian@hotmail.com

    X 513

    A

    猜你喜歡
    方法
    中醫(yī)特有的急救方法
    中老年保健(2021年9期)2021-08-24 03:52:04
    高中數(shù)學(xué)教學(xué)改革的方法
    化學(xué)反應(yīng)多變幻 “虛擬”方法幫大忙
    變快的方法
    兒童繪本(2020年5期)2020-04-07 17:46:30
    學(xué)習(xí)方法
    可能是方法不對(duì)
    用對(duì)方法才能瘦
    Coco薇(2016年2期)2016-03-22 02:42:52
    最有效的簡(jiǎn)單方法
    山東青年(2016年1期)2016-02-28 14:25:23
    四大方法 教你不再“坐以待病”!
    Coco薇(2015年1期)2015-08-13 02:47:34
    賺錢(qián)方法
    听说在线观看完整版免费高清| 一区二区三区高清视频在线| 亚洲成av人片在线播放无| 国产激情久久老熟女| 夜夜爽天天搞| 19禁男女啪啪无遮挡网站| 2021天堂中文幕一二区在线观| 亚洲av片天天在线观看| 亚洲精品456在线播放app | 国产综合懂色| 可以在线观看的亚洲视频| 俺也久久电影网| 亚洲五月天丁香| 精品久久久久久久毛片微露脸| 大型黄色视频在线免费观看| 小说图片视频综合网站| 国产欧美日韩精品一区二区| 亚洲人成伊人成综合网2020| 亚洲国产色片| 成年人黄色毛片网站| 99国产精品一区二区蜜桃av| 久久久久久久久久黄片| 亚洲国产精品sss在线观看| 欧美精品啪啪一区二区三区| 久久国产精品人妻蜜桃| 少妇的丰满在线观看| 亚洲欧美日韩卡通动漫| 淫秽高清视频在线观看| 美女高潮喷水抽搐中文字幕| 久久久国产欧美日韩av| 欧美+亚洲+日韩+国产| 夜夜看夜夜爽夜夜摸| 在线观看66精品国产| 亚洲午夜理论影院| 一级毛片高清免费大全| 国产精品野战在线观看| 欧美最黄视频在线播放免费| 国产成人av激情在线播放| 级片在线观看| 99久久成人亚洲精品观看| 精品久久久久久成人av| 成人鲁丝片一二三区免费| 免费搜索国产男女视频| а√天堂www在线а√下载| 精品久久久久久久毛片微露脸| 99在线人妻在线中文字幕| 日韩欧美三级三区| 免费av毛片视频| 一区福利在线观看| 美女免费视频网站| 97人妻精品一区二区三区麻豆| 每晚都被弄得嗷嗷叫到高潮| 亚洲中文字幕日韩| 日本成人三级电影网站| 日韩精品中文字幕看吧| 亚洲国产精品合色在线| 成人国产综合亚洲| av天堂中文字幕网| 在线观看免费视频日本深夜| 国产精品亚洲一级av第二区| 97人妻精品一区二区三区麻豆| 在线视频色国产色| 国产成人系列免费观看| 成人三级黄色视频| 一a级毛片在线观看| 最好的美女福利视频网| 精品人妻1区二区| 99视频精品全部免费 在线 | 国产av在哪里看| 日本黄大片高清| 国产乱人伦免费视频| 天堂网av新在线| 在线看三级毛片| 欧美日韩福利视频一区二区| 午夜精品在线福利| 午夜福利视频1000在线观看| 久久九九热精品免费| 狠狠狠狠99中文字幕| 午夜亚洲福利在线播放| 国产一区二区激情短视频| 成年女人永久免费观看视频| 一二三四社区在线视频社区8| 五月玫瑰六月丁香| 久久精品国产亚洲av香蕉五月| 岛国视频午夜一区免费看| 久久久久久国产a免费观看| 五月玫瑰六月丁香| 国产精品久久久人人做人人爽| 神马国产精品三级电影在线观看| 在线观看免费视频日本深夜| 黄色成人免费大全| 老熟妇乱子伦视频在线观看| 少妇熟女aⅴ在线视频| 亚洲无线观看免费| 亚洲无线观看免费| 熟妇人妻久久中文字幕3abv| 天堂网av新在线| 制服人妻中文乱码| 国产精品永久免费网站| 欧美国产日韩亚洲一区| 狂野欧美激情性xxxx| 一进一出抽搐动态| 欧美日韩一级在线毛片| bbb黄色大片| 日韩有码中文字幕| 国产精品美女特级片免费视频播放器 | 狠狠狠狠99中文字幕| 日韩欧美免费精品| cao死你这个sao货| 免费人成视频x8x8入口观看| 久久亚洲真实| 亚洲av熟女| av视频在线观看入口| а√天堂www在线а√下载| 99国产精品一区二区蜜桃av| 99久久无色码亚洲精品果冻| 日韩高清综合在线| 日本一二三区视频观看| 亚洲色图 男人天堂 中文字幕| 99热6这里只有精品| 两人在一起打扑克的视频| 老司机午夜十八禁免费视频| aaaaa片日本免费| 99久国产av精品| 丰满人妻熟妇乱又伦精品不卡| 亚洲人成网站高清观看| 成人欧美大片| 国产午夜精品论理片| 曰老女人黄片| 九九在线视频观看精品| 亚洲中文字幕日韩| 又粗又爽又猛毛片免费看| 天天躁狠狠躁夜夜躁狠狠躁| 婷婷六月久久综合丁香| 在线观看一区二区三区| 在线永久观看黄色视频| 国产黄a三级三级三级人| 亚洲av电影在线进入| a级毛片a级免费在线| 搡老妇女老女人老熟妇| 国产黄色小视频在线观看| 国产黄色小视频在线观看| 亚洲精华国产精华精| 午夜免费激情av| 18禁黄网站禁片午夜丰满| 亚洲人成伊人成综合网2020| 久久久国产成人精品二区| 熟女少妇亚洲综合色aaa.| 国产久久久一区二区三区| 久久久久精品国产欧美久久久| 亚洲av成人精品一区久久| 高潮久久久久久久久久久不卡| 啦啦啦韩国在线观看视频| 高清在线国产一区| 99久久无色码亚洲精品果冻| 熟女电影av网| 亚洲国产欧洲综合997久久,| 国产精品亚洲av一区麻豆| 国产不卡一卡二| av在线蜜桃| 亚洲成人久久性| 精品久久久久久,| 久久久久久久久中文| 可以在线观看的亚洲视频| 在线免费观看不下载黄p国产 | 亚洲国产欧美一区二区综合| 嫩草影院入口| 日本 欧美在线| 午夜免费激情av| 国产成人av教育| 老司机福利观看| 啦啦啦观看免费观看视频高清| 国产欧美日韩精品一区二区| 国产精品亚洲美女久久久| 白带黄色成豆腐渣| 日韩精品中文字幕看吧| 成年女人看的毛片在线观看| 久久精品人妻少妇| 亚洲国产精品sss在线观看| 国产私拍福利视频在线观看| 综合色av麻豆| 欧美成人性av电影在线观看| 好男人在线观看高清免费视频| 亚洲人与动物交配视频| 免费看日本二区| 亚洲成a人片在线一区二区| 天堂影院成人在线观看| 国产1区2区3区精品| 99国产精品一区二区三区| 中文亚洲av片在线观看爽| 欧美+亚洲+日韩+国产| 国产伦精品一区二区三区视频9 | 特级一级黄色大片| 久久国产精品影院| 人人妻,人人澡人人爽秒播| 国产免费av片在线观看野外av| 久久久水蜜桃国产精品网| 日本a在线网址| 免费观看的影片在线观看| 日韩有码中文字幕| 亚洲天堂国产精品一区在线| 黑人欧美特级aaaaaa片| 黑人操中国人逼视频| 岛国在线观看网站| 国产精品综合久久久久久久免费| 亚洲片人在线观看| 男人舔女人下体高潮全视频| 黄片小视频在线播放| 十八禁网站免费在线| 黄频高清免费视频| 毛片女人毛片| 亚洲欧美激情综合另类| 日本黄色片子视频| 日本在线视频免费播放| 国产成人精品久久二区二区免费| 黄色成人免费大全| 国产av在哪里看| 一二三四社区在线视频社区8| 我要搜黄色片| 男人和女人高潮做爰伦理| 桃色一区二区三区在线观看| 日本一本二区三区精品| 亚洲国产精品合色在线| 91字幕亚洲| 88av欧美| 国产成人精品无人区| 国产成人一区二区三区免费视频网站| 真实男女啪啪啪动态图| 少妇熟女aⅴ在线视频| 精品电影一区二区在线| 久久精品国产亚洲av香蕉五月| 舔av片在线| 给我免费播放毛片高清在线观看| 精品99又大又爽又粗少妇毛片 | 欧美国产日韩亚洲一区| 国产av麻豆久久久久久久| 亚洲午夜精品一区,二区,三区| 神马国产精品三级电影在线观看| 一个人免费在线观看的高清视频| 国产日本99.免费观看| 一级毛片精品| 亚洲国产精品合色在线| 国产精品影院久久| 亚洲,欧美精品.| 狠狠狠狠99中文字幕| 亚洲成av人片免费观看| 国产激情久久老熟女| 亚洲av中文字字幕乱码综合| 全区人妻精品视频| 最新中文字幕久久久久 | 男女视频在线观看网站免费| 久久久久性生活片| 天堂网av新在线| 老司机在亚洲福利影院| 五月伊人婷婷丁香| 90打野战视频偷拍视频| 亚洲色图 男人天堂 中文字幕| 法律面前人人平等表现在哪些方面| 欧美绝顶高潮抽搐喷水| 美女免费视频网站| 成在线人永久免费视频| 麻豆国产av国片精品| 可以在线观看毛片的网站| 精品一区二区三区视频在线观看免费| 午夜免费成人在线视频| 久久久久久久午夜电影| 亚洲成av人片免费观看| www.熟女人妻精品国产| 久久久久国产精品人妻aⅴ院| 亚洲精品在线观看二区| 国产在线精品亚洲第一网站| 老司机深夜福利视频在线观看| 两个人的视频大全免费| 国产伦精品一区二区三区视频9 | 成人三级黄色视频| 听说在线观看完整版免费高清| 亚洲欧美日韩高清专用| svipshipincom国产片| 国内精品一区二区在线观看| 男人舔女人的私密视频| 亚洲国产欧洲综合997久久,| 这个男人来自地球电影免费观看| 精品久久蜜臀av无| 看黄色毛片网站| 久久国产精品影院| 男女下面进入的视频免费午夜| 黄频高清免费视频| 亚洲一区高清亚洲精品| 99久久精品一区二区三区| 亚洲中文字幕一区二区三区有码在线看 | 亚洲激情在线av| 国产亚洲精品久久久com| 亚洲欧美激情综合另类| 色尼玛亚洲综合影院| 久久人妻av系列| 淫妇啪啪啪对白视频| 国产熟女xx| 亚洲专区国产一区二区| 三级男女做爰猛烈吃奶摸视频| 全区人妻精品视频| av福利片在线观看| 一a级毛片在线观看| 午夜亚洲福利在线播放| 久久久国产成人免费| 91av网站免费观看| 99久久国产精品久久久| 午夜精品一区二区三区免费看| 网址你懂的国产日韩在线| 中文字幕人成人乱码亚洲影| 午夜福利免费观看在线| 日韩av在线大香蕉| 两人在一起打扑克的视频| 国产探花在线观看一区二区| 亚洲欧美日韩东京热| 色在线成人网| 舔av片在线| 18禁黄网站禁片午夜丰满| 小蜜桃在线观看免费完整版高清| 国产午夜福利久久久久久| 精品国产美女av久久久久小说| 免费大片18禁| 久久这里只有精品19| 精品国产乱子伦一区二区三区| 国产精品98久久久久久宅男小说| 久久久国产欧美日韩av| 色综合站精品国产| 在线观看免费午夜福利视频| 亚洲欧美日韩高清在线视频| 欧美三级亚洲精品| 他把我摸到了高潮在线观看| 亚洲aⅴ乱码一区二区在线播放| 国产亚洲欧美98| 国产欧美日韩精品一区二区| 国产极品精品免费视频能看的| 久久精品国产亚洲av香蕉五月| 少妇裸体淫交视频免费看高清| 久久午夜亚洲精品久久| av黄色大香蕉| 成人国产综合亚洲| 天天躁狠狠躁夜夜躁狠狠躁| 淫秽高清视频在线观看| 国产美女午夜福利| 看黄色毛片网站| 99国产极品粉嫩在线观看| 亚洲精品色激情综合| 午夜免费激情av| 国产激情欧美一区二区| 免费在线观看影片大全网站| av天堂在线播放| 日韩欧美免费精品| 欧美成人免费av一区二区三区| 国内精品一区二区在线观看| 日本 av在线| 午夜福利在线观看吧| 亚洲欧洲精品一区二区精品久久久| 一级黄色大片毛片| 日本精品一区二区三区蜜桃| 成人av在线播放网站| 久久久久久久久中文| 亚洲成人精品中文字幕电影| 免费看美女性在线毛片视频| 精品日产1卡2卡| h日本视频在线播放| 亚洲av美国av| 女警被强在线播放| 亚洲一区二区三区不卡视频| 午夜a级毛片| 久久性视频一级片| 国内精品久久久久精免费| 亚洲一区二区三区色噜噜| 久久久久久国产a免费观看| 黄色片一级片一级黄色片| 国产亚洲精品av在线| 欧美av亚洲av综合av国产av| 日本a在线网址| 国内久久婷婷六月综合欲色啪| 黄色成人免费大全| 国产精品一及| 久久久色成人| 久久久久久久久中文| 国产野战对白在线观看| 亚洲aⅴ乱码一区二区在线播放| 日韩欧美国产在线观看| 精品电影一区二区在线| 国产伦人伦偷精品视频| 伦理电影免费视频| 成人国产一区最新在线观看| 精品国产乱码久久久久久男人| 色综合欧美亚洲国产小说| 国产黄片美女视频| 亚洲美女视频黄频| 男人的好看免费观看在线视频| 国产亚洲欧美在线一区二区| 免费观看人在逋| 成人亚洲精品av一区二区| 两性午夜刺激爽爽歪歪视频在线观看| 国产不卡一卡二| 久久久久九九精品影院| 久久久精品大字幕| 每晚都被弄得嗷嗷叫到高潮| 久久久久久国产a免费观看| 人妻久久中文字幕网| 一级黄色大片毛片| 亚洲国产欧美人成| 首页视频小说图片口味搜索| 欧美绝顶高潮抽搐喷水| 久久久成人免费电影| 桃红色精品国产亚洲av| 免费一级毛片在线播放高清视频| 午夜福利在线观看免费完整高清在 | 精品99又大又爽又粗少妇毛片 | 成人高潮视频无遮挡免费网站| 亚洲国产看品久久| 亚洲五月婷婷丁香| 9191精品国产免费久久| 在线观看一区二区三区| 两性夫妻黄色片| 成人一区二区视频在线观看| 亚洲美女黄片视频| 99国产综合亚洲精品| 日韩欧美在线二视频| 日韩欧美 国产精品| 叶爱在线成人免费视频播放| 亚洲国产中文字幕在线视频| 午夜福利在线观看免费完整高清在 | 国产视频一区二区在线看| 99在线人妻在线中文字幕| 桃色一区二区三区在线观看| 在线a可以看的网站| 十八禁网站免费在线| 欧美黄色淫秽网站| 亚洲国产精品合色在线| 天堂动漫精品| 狂野欧美激情性xxxx| 国产亚洲av高清不卡| 亚洲熟妇熟女久久| 性欧美人与动物交配| 欧美高清成人免费视频www| 欧美日韩黄片免| 91av网站免费观看| 亚洲av成人精品一区久久| 波多野结衣高清无吗| 一边摸一边抽搐一进一小说| 波多野结衣巨乳人妻| 丰满的人妻完整版| 性欧美人与动物交配| 国产午夜福利久久久久久| 啦啦啦观看免费观看视频高清| 精品一区二区三区av网在线观看| 免费看a级黄色片| 久久午夜亚洲精品久久| 天天一区二区日本电影三级| 波多野结衣高清作品| 精品久久蜜臀av无| 欧美三级亚洲精品| 51午夜福利影视在线观看| 中文字幕最新亚洲高清| 亚洲真实伦在线观看| 亚洲五月婷婷丁香| 99国产精品一区二区三区| 国产成人福利小说| 日本在线视频免费播放| 此物有八面人人有两片| 久久久久国内视频| 好看av亚洲va欧美ⅴa在| 久久久久久久久久黄片| 中出人妻视频一区二区| 免费在线观看影片大全网站| 两个人视频免费观看高清| 国产激情偷乱视频一区二区| 久久天堂一区二区三区四区| 97超级碰碰碰精品色视频在线观看| 免费无遮挡裸体视频| 色综合站精品国产| 高清在线国产一区| 男女午夜视频在线观看| 真人一进一出gif抽搐免费| h日本视频在线播放| 欧美性猛交黑人性爽| 真实男女啪啪啪动态图| 美女大奶头视频| 精品熟女少妇八av免费久了| 久久香蕉精品热| 日日干狠狠操夜夜爽| 窝窝影院91人妻| 国产91精品成人一区二区三区| 中文字幕精品亚洲无线码一区| 夜夜看夜夜爽夜夜摸| 狠狠狠狠99中文字幕| 小蜜桃在线观看免费完整版高清| 欧美黑人欧美精品刺激| 黄色丝袜av网址大全| 欧美在线黄色| bbb黄色大片| 麻豆成人av在线观看| 精品免费久久久久久久清纯| 国产精品野战在线观看| 黄色女人牲交| 成人亚洲精品av一区二区| 日韩欧美国产一区二区入口| 男插女下体视频免费在线播放| 日本熟妇午夜| 久久久国产精品麻豆| 在线播放国产精品三级| 欧美在线黄色| 两人在一起打扑克的视频| 免费在线观看影片大全网站| 欧美黑人巨大hd| 国产精品乱码一区二三区的特点| 露出奶头的视频| 亚洲人成网站高清观看| 日本 av在线| 特大巨黑吊av在线直播| 亚洲中文av在线| 国内久久婷婷六月综合欲色啪| 日韩免费av在线播放| 亚洲五月婷婷丁香| 黄色日韩在线| 两个人视频免费观看高清| 亚洲国产日韩欧美精品在线观看 | 18禁黄网站禁片免费观看直播| 麻豆国产av国片精品| 亚洲av片天天在线观看| 国产高清激情床上av| 欧美一级毛片孕妇| 亚洲国产高清在线一区二区三| 久久精品影院6| 九色国产91popny在线| 久久精品aⅴ一区二区三区四区| 免费看日本二区| 精品人妻1区二区| 别揉我奶头~嗯~啊~动态视频| 日本黄色视频三级网站网址| 99久久无色码亚洲精品果冻| 日韩欧美在线乱码| 丰满的人妻完整版| 一级作爱视频免费观看| 观看免费一级毛片| 日韩中文字幕欧美一区二区| 色综合亚洲欧美另类图片| 成人无遮挡网站| 日韩欧美国产在线观看| 999精品在线视频| 国产午夜精品久久久久久| 亚洲在线观看片| 婷婷六月久久综合丁香| 久久精品国产99精品国产亚洲性色| 日韩精品青青久久久久久| 两性夫妻黄色片| 在线观看免费视频日本深夜| 亚洲黑人精品在线| 免费无遮挡裸体视频| 最近最新中文字幕大全免费视频| 国产亚洲精品一区二区www| 国产精品爽爽va在线观看网站| 中文字幕av在线有码专区| 欧美日韩中文字幕国产精品一区二区三区| 天堂网av新在线| 久久久水蜜桃国产精品网| 国内精品久久久久精免费| 精品国产美女av久久久久小说| 午夜成年电影在线免费观看| 日本 av在线| www.熟女人妻精品国产| 欧美午夜高清在线| av黄色大香蕉| 丰满的人妻完整版| 高清在线国产一区| 久久中文看片网| 日本一本二区三区精品| 国产91精品成人一区二区三区| 最新美女视频免费是黄的| 91麻豆精品激情在线观看国产| 午夜久久久久精精品| 日韩欧美一区二区三区在线观看| 91麻豆av在线| 国产又色又爽无遮挡免费看| 日韩欧美在线二视频| 久久精品国产清高在天天线| 国产高清三级在线| 久久久久久久久久黄片| 精品一区二区三区av网在线观看| АⅤ资源中文在线天堂| 中文字幕久久专区| 国产欧美日韩精品亚洲av| 啦啦啦观看免费观看视频高清| 久久久水蜜桃国产精品网| 麻豆一二三区av精品| 亚洲性夜色夜夜综合| 久久精品亚洲精品国产色婷小说| 男女之事视频高清在线观看| 亚洲欧美日韩东京热| 欧美在线一区亚洲| 黄色日韩在线| 丝袜人妻中文字幕| 国产精品亚洲一级av第二区| 青草久久国产| 免费在线观看影片大全网站| 欧美日韩一级在线毛片| 男人舔女人的私密视频| 久久久久性生活片| 国产精品久久久av美女十八| 国产激情久久老熟女| 在线永久观看黄色视频| 亚洲欧洲精品一区二区精品久久久| 婷婷六月久久综合丁香| 级片在线观看| 床上黄色一级片| 久久久久久久久中文| 成在线人永久免费视频| 99国产精品一区二区三区| 嫩草影院入口| 日本黄大片高清| 操出白浆在线播放| 男女午夜视频在线观看|