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

    基于Additive Runge-Kutta方法的激波聚焦起爆高精度數(shù)值模擬

    2016-11-18 09:21:14王成宋清官
    關(guān)鍵詞:凹面工作頻率橢球

    王成, 宋清官

    (北京理工大學(xué) 爆炸科學(xué)與技術(shù)國家重點(diǎn)實(shí)驗(yàn)室, 北京 100081)

    ?

    基于Additive Runge-Kutta方法的激波聚焦起爆高精度數(shù)值模擬

    王成, 宋清官

    (北京理工大學(xué) 爆炸科學(xué)與技術(shù)國家重點(diǎn)實(shí)驗(yàn)室, 北京 100081)

    基于詳細(xì)氫氧化學(xué)動力學(xué)模型,建立了描述氫氧爆轟的多組分反應(yīng)歐拉方程組. 針對建立的反應(yīng)歐拉方程組,數(shù)值方法上采用3階Additive Runge-Kutta方法對時(shí)間項(xiàng)進(jìn)行積分,采用5階精度的加權(quán)本質(zhì)無振蕩(WENO)格式對空間對流項(xiàng)進(jìn)行離散,自主研發(fā)了大規(guī)模高精度計(jì)算程序. 該程序能夠處理化學(xué)反應(yīng)源項(xiàng)引起的剛性問題,且能節(jié)省計(jì)算時(shí)間和計(jì)算內(nèi)存. 對半球型、半橢球型、圓錐型3種結(jié)構(gòu)形式凹面腔內(nèi)的激波聚焦起爆過程進(jìn)行了數(shù)值模擬,數(shù)值模擬研究得到了不同結(jié)構(gòu)形式凹面腔內(nèi)的激波聚焦起爆過程.

    詳細(xì)化學(xué)動力學(xué)模型; Additive Runge-Kutta方法; WENO格式;激波聚焦

    氣相爆轟是一個超音速燃燒過程,依靠波陣面后反應(yīng)區(qū)能量的快速釋放來自維持,因此氣相爆轟是燃料氧氣混合物釋放燃料化學(xué)能極其充分的方式,因而具有較高的能量利用效率[1]. 脈沖爆轟發(fā)動機(jī)(pulse detonation engine, PDE)利用了氣相爆轟的特點(diǎn),在發(fā)動機(jī)性能提升方面擁有巨大優(yōu)勢,正逐漸成為了國內(nèi)外新型發(fā)動機(jī)研究熱點(diǎn). 爆轟波的點(diǎn)火起爆是脈沖爆轟發(fā)動機(jī)工作的關(guān)鍵條件,激波聚焦能夠在可燃?xì)怏w中產(chǎn)生高溫高壓區(qū),進(jìn)而可以點(diǎn)火起爆形成爆轟波,是脈沖爆轟發(fā)動機(jī)點(diǎn)火起爆研究的一個重要方向.

    國外在20世紀(jì)90年代已經(jīng)開始了激波聚焦方面的研究,且發(fā)展較為成熟. Izumi等[2]對惰性激波在拋物形或楔形凹面腔的反射問題進(jìn)行了實(shí)驗(yàn)和數(shù)值模擬研究;Levin等[3]提出了基于激波聚焦起爆的兩級脈沖爆轟發(fā)動機(jī)概念,解決了脈沖爆轟發(fā)動機(jī)的點(diǎn)火起爆問題,并使得提高脈沖爆轟發(fā)動機(jī)的工作頻率成為可能. 國內(nèi)激波聚焦問題研究起步較晚. 董剛等[4]數(shù)值模擬研究了氣體反應(yīng)活性、入射激波強(qiáng)度及反射壁面形狀對激波聚焦起爆的影響;曾昊等[5]開展了兩級脈沖爆轟發(fā)動機(jī)的相關(guān)研究,數(shù)值模擬了凹面腔內(nèi)的激波聚焦起爆過程,分析了環(huán)形射流噴口位置、環(huán)形射流初始壓力、入射噴口寬度、凹面腔結(jié)構(gòu)對激波聚焦起爆過程的影響,并進(jìn)行了初步的原理性實(shí)驗(yàn)研究,研究了導(dǎo)流角、氣流出口面積、擴(kuò)張角尾噴管、凹面腔與射流入口間距離等對激波聚焦的影響.

    為進(jìn)一步研究凹面腔內(nèi)的激波聚焦起爆過程,本文基于詳細(xì)氫氧化學(xué)動力學(xué)模型[6],數(shù)值方法上采用3階Additive Runge-Kutta方法對時(shí)間項(xiàng)進(jìn)行積分,空間對流項(xiàng)的離散采用5階精度的加權(quán)本質(zhì)無震蕩(WENO)格式[7],分別對半球型、半橢球型和圓錐型3種結(jié)構(gòu)的凹面腔的激波聚焦起爆過程進(jìn)行了數(shù)值模擬,對比研究了不同凹面腔結(jié)構(gòu)的激波聚焦規(guī)律.

    1 控制方程

    本文采用多組分反應(yīng)Euler方程組描述氫氧預(yù)混氣體的爆轟化學(xué)反應(yīng)過程

    (1)

    (2)

    (3)

    (4)

    式中:p,ρ,e,T,u,v,h分別為壓力、密度、單位質(zhì)量的總能、溫度、x和y方向的速度以及單位質(zhì)量的焓;Y1,Y2,…,YN-1分別為基元反應(yīng)各組分的質(zhì)量分?jǐn)?shù);R為普適氣體常數(shù);ω1,ω2,…,ωN-1分別為各組分的質(zhì)量生成率.

    為準(zhǔn)確地描述氫氧預(yù)混氣體的爆轟反應(yīng)過程,化學(xué)反應(yīng)模型采用9組分20步基元反應(yīng)模型[1]. 參與化學(xué)反應(yīng)的各組分分別為 H2,O2,OH,O,H,H2O2,HO2,H2O,此外氫氧預(yù)混氣體用惰性氣體N2稀釋,N2不參與化學(xué)反應(yīng). 化學(xué)反應(yīng)方程可統(tǒng)一寫為

    (5)

    kfm=AmTnmexp-EamR0T,

    kpm=exp∑9i=1υ″mi-υ′mis0iRi-hiRiT,

    kem=kpmexpp0R0T∑9i=1(υ″mi-υ′mi),

    kbm=kfm/(kem).(6)

    (7)

    式中:Wi為第i組分的摩爾質(zhì)量;γTB為第三體效應(yīng)指數(shù).

    2 數(shù)值方法

    在數(shù)值計(jì)算方法方面,基于WENO-LF[8]對物理通量進(jìn)行分裂,方程(1)離散為

    (8)

    由于化學(xué)反應(yīng)模型采用了9組分20步基元反應(yīng)模型,對流項(xiàng)與化學(xué)反應(yīng)源項(xiàng)之間,以及各個化學(xué)反應(yīng)源項(xiàng)之間時(shí)間尺度差別很大,導(dǎo)致歐拉方程組出現(xiàn)嚴(yán)重的剛性問題. 為了解決剛性問題,本文采用了顯隱式Additive Runge-Kutta方法對時(shí)間方向進(jìn)行積分[10]. 該方法能夠解決下述形式的常微分方程組:

    (9)

    式中:F(U)為N項(xiàng)的總和,對于多組分反應(yīng)歐拉方程組,N=2,即方程右端F(U)包含對流項(xiàng)Fns(U)和化學(xué)反應(yīng)源項(xiàng)Fs(U).

    AdditiveRunge-Kutta方法的主要思想是將剛性項(xiàng)Fns(U)和非剛性項(xiàng)Fs(U)分開處理,對非剛性項(xiàng)Fs(U)運(yùn)用顯式Runge-Kutta(ERK)方法進(jìn)行顯式積分,對剛性項(xiàng)Fns(U)運(yùn)用對角隱式Runge-Kutta(ESDIRK)方法進(jìn)行隱式積分. s步AdditiveRunge-Kutta法的具體使用方法如下:

    (10)

    3 計(jì)算模型及結(jié)果分析

    3.1 一維氣相爆轟

    預(yù)混氣體H2/O2/N2的體積比為0.50/0.25/0.25,初始溫度T0=300 K,初始壓力pa=0.101 MPa. 在左側(cè)封閉端設(shè)置一段高溫高壓點(diǎn)火區(qū),點(diǎn)火溫度Tin=3 000 K,點(diǎn)火壓力pa=0.5 MPa. 數(shù)值方法上對于空間對流項(xiàng)的離散采用5階精度的加權(quán)本質(zhì)無震蕩(WENO)格式.

    為研究網(wǎng)格尺寸對計(jì)算結(jié)果的影響,分別取網(wǎng)格尺寸為2.00,1.00,0.50,0.25,0.10,0.05 mm.

    表1 不同網(wǎng)格尺寸下的爆轟參數(shù)

    從表1中可以看出,平均爆速、CJ爆壓和壁面壓力在不同網(wǎng)格尺寸下基本趨于相同值,即網(wǎng)格尺寸對其影響不大,且各物理量值與理論值基本吻合;而Von Neumann峰壓力隨著網(wǎng)格尺寸減小而增大,且隨著網(wǎng)格尺寸的減少而收斂. 綜上所述,網(wǎng)格尺寸取為0.1 mm時(shí)滿足計(jì)算精度,既能準(zhǔn)確的計(jì)算各爆轟參數(shù),又能減少網(wǎng)格數(shù),節(jié)省計(jì)算時(shí)間.

    選取網(wǎng)格尺寸為0.1 mm,通過對比各組分的體積分?jǐn)?shù),容易發(fā)現(xiàn)各組分變化的量級差別很大,H2,O2,H2O變化的量級為10-1,OH,O,H變化的量級為10-2,H2O2和HO2變化的量級分別為10-4和10-6且衰減很快,如圖1(g)~1(h),因此不同組分在氫氧預(yù)混氣體爆轟波的化學(xué)反應(yīng)中起的作用并不相同,甚至差別很大,可以忽略對反應(yīng)影響較小的組分H2O2和HO2,對基元反應(yīng)模型做進(jìn)一步簡化,簡化的化學(xué)反應(yīng)模型應(yīng)用到下面的激波聚焦起爆研究中.

    3.2 激波聚焦起爆

    為了研究不同凹面腔結(jié)構(gòu)內(nèi)的激波聚焦起爆過程,本文選取了半球型、半橢球型和圓錐型3種結(jié)構(gòu)的凹面腔作為計(jì)算模型,如下:

    工況1:直徑為50 mm,開口直徑為25 mm的半球凹面腔;

    工況2:長徑為50 mm,短徑為25 mm的半橢球凹面腔;

    工況3:高為25 mm,底面圓直徑為25 mm的圓錐凹面腔.

    凹面腔和尾噴管為內(nèi)區(qū)域(實(shí)線),噴管外為外區(qū)域(虛線),如圖2所示. 對于內(nèi)區(qū)域,填充體積比0.23/0.16/0.61的H2/O2/N2預(yù)混氣體,等截面尾噴管軸線方向長12 mm,環(huán)形射流入口為壓力入口邊界,環(huán)形射流入射寬度2.8 mm,壓力pin=0.45 MPa,溫度Tin=450 K;對于外區(qū)域,填充空氣,壓力pa=0.101 MPa,溫度Ta=300 K. 凹面腔及尾噴管為剛性、無滑移、絕熱壁面. 數(shù)值方法上采用3階Additive Runge-Kutta方法對時(shí)間項(xiàng)進(jìn)行積分,空間對流項(xiàng)的離散采用5階精度的加權(quán)本質(zhì)無震蕩(WENO)格式. 從初始網(wǎng)格尺寸為0.1 mm. 選取凹面腔中軸線為壓力觀測線,選取凹面腔中軸線上的點(diǎn)M作為壓力觀測點(diǎn),如圖2所示.

    從圖3中可以看出,高溫高壓的環(huán)形射流由入口噴入半球形凹面腔內(nèi),初始的溫度、壓力不足以起爆預(yù)混氣體,在凹面腔內(nèi)形成了兩道激波,如圖3(a)所示. 隨著激波的傳播,激波的強(qiáng)度逐漸減弱,如圖3(b)所示,但是當(dāng)上下激波相遇碰撞時(shí),壓力增加,激波的強(qiáng)度加強(qiáng),此時(shí)仍未能起爆預(yù)混氣體,如圖3(c) 所示.

    激波繼續(xù)向前傳播,在半球型壁面處發(fā)生反射,由于壁面收斂的作用,壓力逐漸增加,形成了更強(qiáng)的激波. 激波繼續(xù)向前傳播,在半球凹面腔中心處發(fā)生匯聚,如圖3(d)所示,當(dāng)激波在球形凹面腔中心處發(fā)生匯聚時(shí),壓力密度等參量急劇增加,達(dá)到了氫氧預(yù)混氣體的臨界起爆條件,發(fā)生起爆,如圖3(e)所示. 起爆之后,爆轟波在球形凹面腔內(nèi)迅速傳播,如圖3(f)~3(h)所示,實(shí)現(xiàn)了脈沖爆轟發(fā)動機(jī)的起爆過程.

    圖3~圖6可以看出,對于半球凹面腔,在t=310.72 μs時(shí),在中軸線上發(fā)生激波聚焦起爆,起爆壓力約為3.0 MPa;對于半橢球凹面腔,在t=330.05 μs時(shí),在中軸線上發(fā)生激波聚焦起爆,起爆壓力約為2.5 MPa;對于圓錐凹面腔,在t=307.77 μs時(shí),在中軸線附近上發(fā)生激波聚焦起爆,起爆壓力約為2.0 MPa,同時(shí)發(fā)現(xiàn)初始起爆點(diǎn)明顯位于中軸線偏上的位置,本文中給出的入流條件和物理模型均是關(guān)于中軸線對稱的,由此推測這種現(xiàn)象是圓錐凹面腔頂點(diǎn)的奇異性導(dǎo)致的數(shù)值誤差造成的.

    通過對比3種結(jié)構(gòu)的凹面腔的激波聚焦過程,發(fā)現(xiàn)在此初始條件下,3種結(jié)構(gòu)的凹面腔均實(shí)現(xiàn)了激波聚焦起爆過程,且聚焦中心基本均位于中軸線上,壁面收斂導(dǎo)致的激波增強(qiáng)效果不足以起爆預(yù)混氣體. 由于半橢球凹面腔的曲率要比半球凹面腔和圓錐凹面腔的曲率要大,因此半橢球凹面腔的激波聚焦效果最為明顯,聚焦中心位置最為靠前,如圖3~7所示. 此外,通過對比發(fā)現(xiàn),3種結(jié)構(gòu)的凹面腔實(shí)現(xiàn)聚焦起爆時(shí)間不同,其中圓錐凹面腔最早實(shí)現(xiàn)聚焦起爆,半球凹面腔次之,半橢球凹面腔最晚實(shí)現(xiàn)聚焦起爆.

    此外本文還研究了不同凹面腔結(jié)構(gòu)對脈沖爆轟發(fā)動機(jī)工作頻率的影響. 由分析可知,凹面腔內(nèi)預(yù)混氣體起爆之后產(chǎn)生大量高溫高壓氣體,氣體由尾噴管噴出形成推力,氣體噴出后在凹面腔內(nèi)形成了低溫低壓區(qū),此時(shí)入口處的氣體總壓要高于凹面腔內(nèi)的氣體總壓,此時(shí)高溫高壓的射流再次由入口進(jìn)入凹面腔內(nèi),再次實(shí)現(xiàn)聚焦爆轟過程,進(jìn)入下一個循環(huán),因此凹面腔內(nèi)形成了周期性的脈沖爆轟波,實(shí)現(xiàn)了循環(huán)起爆. 從圖中容易看出,對于3種結(jié)構(gòu)的凹面腔,觀測點(diǎn)M的壓力均在400 μs后發(fā)生了周期性的變化,因此對于3種結(jié)構(gòu)的凹面腔,起爆之后均實(shí)現(xiàn)了循環(huán)起爆. 其中,半球凹面腔的工作頻率約為16.7 kHz,半橢球凹面腔的工作頻率約為15.0 kHz,圓錐凹面腔的工作頻率約為20.0 kHz,即圓錐凹面腔的工作頻率最高,半球凹面腔的工作頻率次之,半橢球凹面腔的工作頻率最低.

    4 結(jié) 論

    采用3階Additive Runge-Kutta方法對時(shí)間項(xiàng)進(jìn)行積分,空間對流項(xiàng)的離散采用5階精度的加權(quán)本質(zhì)無震蕩(WENO)格式,對脈沖爆轟發(fā)動機(jī)激波聚焦起爆過程進(jìn)行了數(shù)值模擬,重點(diǎn)研究了半球凹面腔、半橢球凹面腔和圓錐凹面腔對激波聚焦起爆過程的影響,并得到了以下結(jié)論:

    ① 3種結(jié)構(gòu)的凹面腔均實(shí)現(xiàn)了激波聚焦起爆過程,且聚焦中心基本均位于中軸線上,其中半橢球凹面腔的激波聚焦效果最為明顯,聚焦中心位置最為靠前,因此為獲得明顯的聚焦效果,應(yīng)使用半橢球凹面腔;

    ② 3種結(jié)構(gòu)的凹面腔實(shí)現(xiàn)聚焦起爆時(shí)間不同,其中圓錐凹面腔最早實(shí)現(xiàn)聚焦起爆,半球凹面腔次之,半橢球凹面腔最晚實(shí)現(xiàn)聚焦起爆,因此為較快實(shí)現(xiàn)聚焦起爆,應(yīng)使用圓錐凹面腔;

    ③ 3種結(jié)構(gòu)的凹面腔均實(shí)現(xiàn)了循環(huán)起爆,其中,圓錐凹面腔的工作頻率最高,半球凹面腔的工作頻率次之,半橢球的工作頻率最低,因此為達(dá)到較高的工作頻率,應(yīng)使用圓錐凹面腔.

    [1] 張博,白春華.氣相爆轟動力學(xué)特征研究進(jìn)展[J].中國科學(xué):物理學(xué) 力學(xué) 天文學(xué),2014,44(7):665-681.

    Zhang Bo, Bai Chunhua. Research progress on the dynamic characteristics of gaseous detonation[J]. Scientia Sinica: Physica, Mechanica & Astronomica, 2014,44(7):665-681. (in Chinese)

    [2] Izumi K, Aso S, Nishida M. Experimental and computational studies focusing processes of shock waves reflected from parabolic reflectors[J]. Shock Waves, 1994,3(3):213-222.

    [3] Levin V A, Nechaev J N, Tarasov A I. A new approach to organizing operation cycles in pulse detonation engines[C]∥Proceedings of High-Speed Deflagration and Detonation: Fundamentals and Control. Moscov: [s.n.], 2001:31-44.

    [4] 董剛,唐敖,葉經(jīng)方,等.激波聚焦誘導(dǎo)點(diǎn)火和爆轟的數(shù)值研究[J].爆炸與沖擊,2006,25(5):437-444.

    Dong Gang, Tang Ao, Ye Jingfang, et al. Numerical studies on initiation and detonation induced by shock wave focusing[J]. Explosion and Shock Waves, 2006,25(5):437-444. (in Chinese)

    [5] 曾昊,何立明,章雄偉,等.環(huán)形射流噴口位置對激波聚焦起爆的影響分析[J].推進(jìn)技術(shù),2011,32(3):437-442.

    Zeng Hao, He Liming, Zhang Xiongwei, et al. Investigation on the influence of jets spout location on detonation initiation via imploding annular shock waves[J]. Journal of Propulsion Technology, 2001,32(3):437-442. (in Chinese)

    [6] 趙慧,李健,郝莉.不同曲率彎管對氣相爆轟波傳播特性的影響[J].北京理工大學(xué)學(xué)報(bào),2014,34(增刊1):176-180.

    Zhao Hui, Li Jian, Hao Li. The influence of different bend curvature on gaseous detonation wave propagation[J]. Transactions of Beijing Institute of Technology, 2014,34(suppl 1):176-180. (in Chinese)

    [7] 藺偉,宋清官,王成,等.濃度梯度對瓦斯爆炸影響的數(shù)值模擬[J].北京理工大學(xué)學(xué)報(bào),2015,35(4):336-340.

    Lin Wei, Song Qingguan, Wang Cheng, et al. Numerical investigation about the effect of concentration gradient on methane explosion[J]. Transactions of Beijing Institute of Technology, 2015,35(4):336-340. (in Chinese)

    [8] Jiang G S, Shu C W. Efficient implementation of weighted ENO schemes[J]. Journal of Computational Physics, 1996,126(1):202-228.

    [9] Chi-Shu W. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[J]. Lecture Notes in Mathematics, 1998,1697(1):325-432.

    [10] Kennedy C A, Carpenter M H. Additive Runge-Kutta schemes for convection diffusion reaction equations[J]. Applied Numerical Mathematics, 2003,44(1-2):139-181.

    (責(zé)任編輯:劉雨)

    Numerical Application of Additive Runge-Kutta Methods on Detonation Initiation with Convergence of Shock Waves

    WANG Cheng, SONG Qing-guan

    (State Key Laboratory of Explosion Science and Technology, Beijing Institute of Technology, Beijing 100081, China)

    Based on detailed hydrogen and oxygen chemical kinetics model, the multi-species reactive Euler equations were established to describe the hydrogen and oxygen detonation. Using third-order Additive Runge-Kutta methods in time discretization, using fifth-order weighted essentially non-oscillatory (WENO) scheme in spatial discretization, a high-solution large-scale program was developed. It was verified that the program could solve stiff source term caused by multispecies chemical reaction, which could also save computation time and computation memory. With this program, the process of detonation initiation with convergence of shock waves was numerically simulated for 3 types of concave cavity, hemispherical concave cavity, semi-ellipsoidal concave cavity and conical concave cavity. The simulation results show the characteristic of detonation initiation in different types of concave cavity.

    detailed chemical kinetics model; Additive Runge-Kutta methods; weighted essentially non-oscillatory scheme; convergence of shock waves

    2015-10-30

    國家自然科學(xué)基金資助項(xiàng)目(11325209,11272056)

    王成(1972—),男,教授,博士生導(dǎo)師,E-mail:wangcheng@bit.edu.cn.

    O 383

    A

    1001-0645(2016)02-0137-07

    10.15918/j.tbit1001-0645.2016.02.006

    猜你喜歡
    凹面工作頻率橢球
    獨(dú)立坐標(biāo)系橢球變換與坐標(biāo)換算
    橢球槽宏程序編制及其Vericut仿真
    智能制造(2021年4期)2021-11-04 08:54:44
    搜集凹面錐體
    橢球精加工軌跡及程序設(shè)計(jì)
    基于外定界橢球集員估計(jì)的純方位目標(biāo)跟蹤
    CPU故障的處理技巧
    中文信息(2016年11期)2017-02-11 16:48:48
    減壓小心機(jī)
    決策探索(2016年21期)2016-11-28 09:30:57
    高頻非對稱EMCCD增益驅(qū)動電路設(shè)計(jì)
    一種工作頻率可變的雙向DC-DC變換器
    新時(shí)期自適應(yīng)通信技術(shù)應(yīng)用發(fā)展現(xiàn)狀
    国产激情久久老熟女| 伦理电影免费视频| 天堂影院成人在线观看| 中文字幕最新亚洲高清| 成人18禁高潮啪啪吃奶动态图| 欧美日本视频| 99久久久亚洲精品蜜臀av| 亚洲成人国产一区在线观看| 亚洲精华国产精华精| 国产精品亚洲av一区麻豆| 脱女人内裤的视频| 午夜福利高清视频| 老汉色av国产亚洲站长工具| 中国美女看黄片| 十八禁人妻一区二区| 久久这里只有精品19| 欧美大码av| 超碰成人久久| 亚洲成人久久性| 亚洲自拍偷在线| 国产一区在线观看成人免费| 啦啦啦免费观看视频1| 久久精品aⅴ一区二区三区四区| 1024香蕉在线观看| 哪里可以看免费的av片| 久久精品国产亚洲av高清一级| 99久久精品国产亚洲精品| 国产亚洲欧美98| 可以在线观看毛片的网站| 国产午夜福利久久久久久| 国产私拍福利视频在线观看| 亚洲精华国产精华精| 国产精品久久电影中文字幕| 在线观看66精品国产| 波多野结衣av一区二区av| 亚洲性夜色夜夜综合| 可以在线观看的亚洲视频| 人妻丰满熟妇av一区二区三区| 在线观看www视频免费| 国产欧美日韩一区二区三| 夜夜爽天天搞| 最新美女视频免费是黄的| 精品不卡国产一区二区三区| 波多野结衣高清作品| 在线观看免费日韩欧美大片| 窝窝影院91人妻| 国产一区二区三区视频了| 精品乱码久久久久久99久播| a在线观看视频网站| 狂野欧美激情性xxxx| 丁香六月欧美| 亚洲熟女毛片儿| 村上凉子中文字幕在线| 久久久久久久久久黄片| 国产精品精品国产色婷婷| 久9热在线精品视频| 免费看日本二区| 黄片播放在线免费| 色婷婷久久久亚洲欧美| 国产99白浆流出| 一a级毛片在线观看| 国产精品二区激情视频| 久久中文字幕一级| 99riav亚洲国产免费| x7x7x7水蜜桃| www.熟女人妻精品国产| 亚洲国产欧美网| 日韩中文字幕欧美一区二区| 男女下面进入的视频免费午夜 | 人妻久久中文字幕网| 欧美精品啪啪一区二区三区| 免费av毛片视频| 1024视频免费在线观看| 国产精品亚洲av一区麻豆| 欧美绝顶高潮抽搐喷水| 日韩欧美一区视频在线观看| 成人手机av| xxx96com| 伊人久久大香线蕉亚洲五| 欧美性猛交黑人性爽| 国产成人av教育| 午夜日韩欧美国产| 国产又色又爽无遮挡免费看| 男人的好看免费观看在线视频 | 妹子高潮喷水视频| 国产一区二区激情短视频| 久99久视频精品免费| 亚洲av中文字字幕乱码综合 | 一卡2卡三卡四卡精品乱码亚洲| 亚洲午夜精品一区,二区,三区| www日本黄色视频网| 最近最新中文字幕大全免费视频| 最新美女视频免费是黄的| 国产一区二区三区视频了| 久久香蕉国产精品| a在线观看视频网站| 亚洲精品一区av在线观看| 久久国产精品男人的天堂亚洲| 亚洲精品美女久久av网站| 男女视频在线观看网站免费 | 久久久久久大精品| a在线观看视频网站| 国产成人av教育| 黄色丝袜av网址大全| 老司机深夜福利视频在线观看| 老汉色∧v一级毛片| avwww免费| 欧美激情高清一区二区三区| 制服诱惑二区| 最新美女视频免费是黄的| 精品国产超薄肉色丝袜足j| 亚洲久久久国产精品| 在线观看日韩欧美| 欧美日韩瑟瑟在线播放| 久久精品aⅴ一区二区三区四区| 男人的好看免费观看在线视频 | 不卡一级毛片| 国产精品av久久久久免费| 首页视频小说图片口味搜索| 亚洲成人免费电影在线观看| av在线天堂中文字幕| 欧美成人免费av一区二区三区| 老熟妇乱子伦视频在线观看| 精品乱码久久久久久99久播| 欧美黄色淫秽网站| 夜夜看夜夜爽夜夜摸| 国产欧美日韩精品亚洲av| 国内精品久久久久久久电影| 久热爱精品视频在线9| 深夜精品福利| 欧美丝袜亚洲另类 | 亚洲一区二区三区不卡视频| 久久99热这里只有精品18| 禁无遮挡网站| 国产免费男女视频| 欧美黑人欧美精品刺激| 成年女人毛片免费观看观看9| 白带黄色成豆腐渣| 女生性感内裤真人,穿戴方法视频| 日本撒尿小便嘘嘘汇集6| 搞女人的毛片| 国产精品一区二区三区四区久久 | 亚洲精品国产区一区二| 91av网站免费观看| 久久热在线av| 搡老岳熟女国产| 欧美成人一区二区免费高清观看 | 在线看三级毛片| 91在线观看av| 久久99热这里只有精品18| 91国产中文字幕| 男男h啪啪无遮挡| 啦啦啦免费观看视频1| a在线观看视频网站| 黄色毛片三级朝国网站| 巨乳人妻的诱惑在线观看| 在线观看午夜福利视频| 中文资源天堂在线| 久久久久久久久久黄片| 精品午夜福利视频在线观看一区| 久久婷婷成人综合色麻豆| 极品教师在线免费播放| 黄片播放在线免费| 午夜福利欧美成人| 操出白浆在线播放| 岛国在线观看网站| 久久久久久大精品| 日韩欧美在线二视频| 欧美亚洲日本最大视频资源| 人成视频在线观看免费观看| 啦啦啦韩国在线观看视频| 无遮挡黄片免费观看| 亚洲精华国产精华精| 国产99久久九九免费精品| 最好的美女福利视频网| 香蕉国产在线看| 国产亚洲精品av在线| 久久狼人影院| 精品一区二区三区四区五区乱码| 亚洲成人免费电影在线观看| 满18在线观看网站| 看片在线看免费视频| 国产在线观看jvid| 人人妻人人澡人人看| 亚洲va日本ⅴa欧美va伊人久久| 黄片小视频在线播放| 免费在线观看亚洲国产| 在线观看www视频免费| 亚洲九九香蕉| 成人18禁高潮啪啪吃奶动态图| 久久欧美精品欧美久久欧美| 一个人观看的视频www高清免费观看 | 久久久久久久精品吃奶| 啦啦啦韩国在线观看视频| 美女大奶头视频| 黄片大片在线免费观看| 色婷婷久久久亚洲欧美| 欧美黄色片欧美黄色片| 青草久久国产| 一级毛片女人18水好多| 欧美丝袜亚洲另类 | 国产av一区二区精品久久| 欧美在线一区亚洲| 91大片在线观看| 一区二区三区激情视频| 日韩欧美一区视频在线观看| а√天堂www在线а√下载| 国产亚洲欧美98| 好看av亚洲va欧美ⅴa在| 国产亚洲精品第一综合不卡| 听说在线观看完整版免费高清| 97超级碰碰碰精品色视频在线观看| 亚洲专区字幕在线| 国产真人三级小视频在线观看| 久久精品91无色码中文字幕| 国产熟女xx| 国产私拍福利视频在线观看| 一卡2卡三卡四卡精品乱码亚洲| ponron亚洲| 母亲3免费完整高清在线观看| 可以在线观看毛片的网站| 午夜a级毛片| 久99久视频精品免费| 午夜激情av网站| 欧美性猛交黑人性爽| 国产片内射在线| 此物有八面人人有两片| 午夜福利在线观看吧| netflix在线观看网站| 国产精品免费一区二区三区在线| 色综合欧美亚洲国产小说| 国产伦在线观看视频一区| 一a级毛片在线观看| 最近在线观看免费完整版| 天天一区二区日本电影三级| 亚洲成av片中文字幕在线观看| 18禁裸乳无遮挡免费网站照片 | 日韩成人在线观看一区二区三区| 国产精品av久久久久免费| 人妻久久中文字幕网| 午夜福利高清视频| 欧美成狂野欧美在线观看| 亚洲va日本ⅴa欧美va伊人久久| 久久精品影院6| 亚洲男人的天堂狠狠| 18禁观看日本| 亚洲精品久久国产高清桃花| 欧美激情久久久久久爽电影| 久久中文字幕人妻熟女| 50天的宝宝边吃奶边哭怎么回事| 国产不卡一卡二| 国产极品粉嫩免费观看在线| 欧美日韩黄片免| 一夜夜www| 成人欧美大片| 亚洲成av片中文字幕在线观看| 人妻久久中文字幕网| 91字幕亚洲| 日韩大码丰满熟妇| 国产精品99久久99久久久不卡| 亚洲av电影在线进入| 国产又黄又爽又无遮挡在线| av欧美777| 国产成人一区二区三区免费视频网站| 12—13女人毛片做爰片一| 男女视频在线观看网站免费 | 无遮挡黄片免费观看| 99在线视频只有这里精品首页| 精品福利观看| 亚洲精品中文字幕一二三四区| 国产激情偷乱视频一区二区| 亚洲国产欧洲综合997久久, | 久久精品影院6| 日韩欧美在线二视频| 国产三级在线视频| 精品欧美国产一区二区三| 手机成人av网站| 色综合站精品国产| 亚洲天堂国产精品一区在线| 久久欧美精品欧美久久欧美| 久久精品国产亚洲av高清一级| xxx96com| 国产精品久久久久久亚洲av鲁大| 制服丝袜大香蕉在线| 国产成人啪精品午夜网站| 精品国产亚洲在线| 中文字幕人妻丝袜一区二区| 99久久国产精品久久久| 中文资源天堂在线| 日韩 欧美 亚洲 中文字幕| 男人的好看免费观看在线视频 | 久久国产精品影院| 可以免费在线观看a视频的电影网站| 精品国产超薄肉色丝袜足j| cao死你这个sao货| 久久精品亚洲精品国产色婷小说| 极品教师在线免费播放| 久久精品国产综合久久久| 国产在线观看jvid| 成年女人毛片免费观看观看9| 12—13女人毛片做爰片一| 久久精品国产亚洲av高清一级| 欧美一级a爱片免费观看看 | 91字幕亚洲| 精品免费久久久久久久清纯| 欧美中文综合在线视频| 12—13女人毛片做爰片一| av在线天堂中文字幕| 美女扒开内裤让男人捅视频| 黄色成人免费大全| 十八禁网站免费在线| 亚洲avbb在线观看| 在线观看一区二区三区| 成人亚洲精品av一区二区| 午夜福利18| 免费在线观看黄色视频的| 国产亚洲精品久久久久久毛片| 国产成人精品无人区| 国产成+人综合+亚洲专区| 在线天堂中文资源库| 日韩欧美国产在线观看| 国产欧美日韩一区二区三| 欧美中文综合在线视频| 日韩欧美一区视频在线观看| 国产精品 欧美亚洲| 最好的美女福利视频网| 国产97色在线日韩免费| 韩国精品一区二区三区| av福利片在线| 欧美一区二区精品小视频在线| 国产精华一区二区三区| 亚洲三区欧美一区| 嫁个100分男人电影在线观看| 国产黄片美女视频| 国产精品久久视频播放| 国产亚洲av高清不卡| 夜夜爽天天搞| 少妇 在线观看| 亚洲aⅴ乱码一区二区在线播放 | 两个人免费观看高清视频| 人妻久久中文字幕网| 日韩成人在线观看一区二区三区| 很黄的视频免费| 少妇粗大呻吟视频| 国产极品粉嫩免费观看在线| 成人国产一区最新在线观看| 国产日本99.免费观看| ponron亚洲| 99在线视频只有这里精品首页| 精品一区二区三区视频在线观看免费| 香蕉国产在线看| 亚洲电影在线观看av| 又大又爽又粗| 亚洲国产高清在线一区二区三 | 国产精品国产高清国产av| 51午夜福利影视在线观看| 日日夜夜操网爽| 久久久水蜜桃国产精品网| 黄片播放在线免费| 久久久水蜜桃国产精品网| 啦啦啦观看免费观看视频高清| 88av欧美| 亚洲精品av麻豆狂野| 国产欧美日韩精品亚洲av| 50天的宝宝边吃奶边哭怎么回事| 哪里可以看免费的av片| 真人做人爱边吃奶动态| 91成年电影在线观看| 国产激情偷乱视频一区二区| 精品久久久久久久毛片微露脸| 热re99久久国产66热| 午夜福利免费观看在线| 一本综合久久免费| 最近最新中文字幕大全免费视频| 欧美精品亚洲一区二区| 久久久久九九精品影院| 麻豆久久精品国产亚洲av| tocl精华| 国内少妇人妻偷人精品xxx网站 | 别揉我奶头~嗯~啊~动态视频| 免费在线观看影片大全网站| 国产极品粉嫩免费观看在线| 日本五十路高清| e午夜精品久久久久久久| 久久人人精品亚洲av| 国产精品亚洲av一区麻豆| 精品久久久久久,| 成年女人毛片免费观看观看9| 97人妻精品一区二区三区麻豆 | 亚洲精品在线观看二区| 国产成人欧美在线观看| avwww免费| 色尼玛亚洲综合影院| 国产蜜桃级精品一区二区三区| 久久久久精品国产欧美久久久| 午夜免费观看网址| 怎么达到女性高潮| 色综合站精品国产| 国产亚洲精品一区二区www| 精品久久久久久久毛片微露脸| 热re99久久国产66热| 国产精品精品国产色婷婷| 中文字幕精品免费在线观看视频| 国内毛片毛片毛片毛片毛片| 亚洲人成网站高清观看| 国产一卡二卡三卡精品| 一级毛片女人18水好多| 嫩草影院精品99| 中文字幕精品免费在线观看视频| 母亲3免费完整高清在线观看| 少妇 在线观看| 青草久久国产| 国产黄片美女视频| 欧美日韩乱码在线| 国产成人精品无人区| 色av中文字幕| 久久久久久九九精品二区国产 | 少妇被粗大的猛进出69影院| 精品熟女少妇八av免费久了| 黄色毛片三级朝国网站| 国产精品免费一区二区三区在线| 日韩欧美在线二视频| 在线十欧美十亚洲十日本专区| 国产精品一区二区三区四区久久 | 精品国产亚洲在线| cao死你这个sao货| 9191精品国产免费久久| 国产一区二区激情短视频| 日韩欧美国产在线观看| 91字幕亚洲| 很黄的视频免费| 神马国产精品三级电影在线观看 | av片东京热男人的天堂| 国内少妇人妻偷人精品xxx网站 | 国产高清激情床上av| 在线av久久热| 丝袜在线中文字幕| 久久午夜亚洲精品久久| 侵犯人妻中文字幕一二三四区| 国产亚洲精品久久久久5区| 50天的宝宝边吃奶边哭怎么回事| 草草在线视频免费看| 99久久精品国产亚洲精品| 国产av在哪里看| 99久久国产精品久久久| 日本免费一区二区三区高清不卡| 这个男人来自地球电影免费观看| 精品久久久久久久人妻蜜臀av| 国产精品日韩av在线免费观看| 看黄色毛片网站| 人妻丰满熟妇av一区二区三区| 久久久久国内视频| 亚洲色图 男人天堂 中文字幕| 精品卡一卡二卡四卡免费| 亚洲无线在线观看| 亚洲av中文字字幕乱码综合 | 黄色a级毛片大全视频| 国产精品免费一区二区三区在线| 中文在线观看免费www的网站 | 日本精品一区二区三区蜜桃| 老司机福利观看| 波多野结衣高清无吗| 美女 人体艺术 gogo| 国内少妇人妻偷人精品xxx网站 | 亚洲激情在线av| 欧美一级毛片孕妇| 亚洲国产毛片av蜜桃av| 国内精品久久久久久久电影| 久久欧美精品欧美久久欧美| 久99久视频精品免费| 成人国产一区最新在线观看| 国产精品美女特级片免费视频播放器 | 波多野结衣高清无吗| 制服诱惑二区| 国产视频一区二区在线看| 18禁裸乳无遮挡免费网站照片 | 免费电影在线观看免费观看| 91麻豆av在线| 欧美av亚洲av综合av国产av| 老司机午夜十八禁免费视频| 久久欧美精品欧美久久欧美| 中文字幕精品亚洲无线码一区 | 最近最新中文字幕大全免费视频| 日日爽夜夜爽网站| 国产男靠女视频免费网站| 成人手机av| 国产精品1区2区在线观看.| 欧美激情高清一区二区三区| 久久国产亚洲av麻豆专区| 成人亚洲精品一区在线观看| 中文字幕av电影在线播放| 啦啦啦 在线观看视频| 窝窝影院91人妻| 视频区欧美日本亚洲| 精品高清国产在线一区| 国产精品av久久久久免费| 欧美性长视频在线观看| 亚洲欧美一区二区三区黑人| 精品高清国产在线一区| 国产精品爽爽va在线观看网站 | 亚洲五月色婷婷综合| 日日干狠狠操夜夜爽| 白带黄色成豆腐渣| 国产精品久久视频播放| 少妇 在线观看| 欧美乱妇无乱码| 国产99白浆流出| 黄色片一级片一级黄色片| 亚洲电影在线观看av| 久9热在线精品视频| 色哟哟哟哟哟哟| tocl精华| 免费在线观看影片大全网站| 欧美人与性动交α欧美精品济南到| 在线观看www视频免费| 桃红色精品国产亚洲av| 999久久久精品免费观看国产| 后天国语完整版免费观看| 18禁观看日本| 亚洲成a人片在线一区二区| 精品久久久久久,| 亚洲精品国产一区二区精华液| 精品国产美女av久久久久小说| 在线av久久热| 国产又爽黄色视频| 女人被狂操c到高潮| 无遮挡黄片免费观看| 国产激情久久老熟女| 久久99热这里只有精品18| 亚洲第一电影网av| 丰满的人妻完整版| 桃色一区二区三区在线观看| 亚洲自拍偷在线| 美女 人体艺术 gogo| 在线观看免费视频日本深夜| 国产黄色小视频在线观看| 成人永久免费在线观看视频| 在线视频色国产色| 国产极品粉嫩免费观看在线| 老汉色∧v一级毛片| 欧美在线一区亚洲| 欧美绝顶高潮抽搐喷水| 免费人成视频x8x8入口观看| 一级黄色大片毛片| 大香蕉久久成人网| 99久久久亚洲精品蜜臀av| 亚洲人成伊人成综合网2020| 村上凉子中文字幕在线| 熟女少妇亚洲综合色aaa.| 男女做爰动态图高潮gif福利片| 99热这里只有精品一区 | 精品免费久久久久久久清纯| 亚洲av片天天在线观看| 人人妻人人澡欧美一区二区| 人妻久久中文字幕网| 亚洲成人国产一区在线观看| 免费在线观看日本一区| 日韩视频一区二区在线观看| 在线十欧美十亚洲十日本专区| 久久精品91无色码中文字幕| 国产一区二区三区视频了| 777久久人妻少妇嫩草av网站| 视频区欧美日本亚洲| 亚洲专区字幕在线| 免费女性裸体啪啪无遮挡网站| 老司机在亚洲福利影院| 国产精品 国内视频| 精品不卡国产一区二区三区| 后天国语完整版免费观看| 好男人电影高清在线观看| 日韩一卡2卡3卡4卡2021年| 午夜福利在线在线| 亚洲成a人片在线一区二区| 国产成人系列免费观看| 国产激情久久老熟女| 亚洲久久久国产精品| 黑人操中国人逼视频| 宅男免费午夜| 日韩欧美免费精品| 国产亚洲精品综合一区在线观看 | 午夜福利18| 日本在线视频免费播放| 十八禁网站免费在线| 国产免费男女视频| 高清毛片免费观看视频网站| 精品国产超薄肉色丝袜足j| 1024手机看黄色片| 男人舔奶头视频| 色av中文字幕| 日韩欧美一区视频在线观看| 级片在线观看| 色综合亚洲欧美另类图片| 欧美久久黑人一区二区| a级毛片a级免费在线| 亚洲最大成人中文| 在线观看午夜福利视频| 黄色 视频免费看| 97超级碰碰碰精品色视频在线观看| av中文乱码字幕在线| 午夜成年电影在线免费观看| 免费av毛片视频| 中文字幕高清在线视频| 免费搜索国产男女视频| 精品一区二区三区视频在线观看免费| 国产麻豆成人av免费视频| 91成人精品电影| 校园春色视频在线观看| 亚洲,欧美精品.| 色播在线永久视频| 黄色 视频免费看| x7x7x7水蜜桃| av在线播放免费不卡| 欧美最黄视频在线播放免费| 免费女性裸体啪啪无遮挡网站| 白带黄色成豆腐渣|