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

    DSMC計算中碰撞對取樣和時間推進環(huán)節(jié)的高效處理方法

    2018-03-09 07:54:54王江峰中國飛行試驗研究院陜西西安710089南京航空航天大學(xué)航空宇航學(xué)院江蘇南京10016
    空氣動力學(xué)學(xué)報 2018年1期
    關(guān)鍵詞:乘數(shù)步長流場

    屈 程, 王江峰(1. 中國飛行試驗研究院, 陜西 西安 710089; . 南京航空航天大學(xué) 航空宇航學(xué)院, 江蘇 南京 10016)

    0 引 言

    DSMC方法是Bird從流動的物理模擬出發(fā)發(fā)展的一種直接模擬蒙特卡洛方法[1],是模擬過渡領(lǐng)域真實氣體流動問題最成功的方法之一[2-5]。然而,DSMC方法的較大缺點是其對較高密度稀薄流問題的模擬能力不強,并且計算耗時太長,這大大限制了該方法的大規(guī)模應(yīng)用。為克服上述缺點,研究者們發(fā)展了多種方法,如采用多線程并行計算[6-7]、設(shè)計高效率分子搜索算法[7]、更改DSMC方法中的數(shù)據(jù)結(jié)構(gòu)[8]、引入碰撞距離思想[9-10]、發(fā)展自適應(yīng)當?shù)貢r間步長方法[11-12]等。

    DSMC方法中,計算網(wǎng)格起兩種作用[13]。第一種是對流場宏觀流動參數(shù)進行空間離散,第二種是促進碰撞分子碰撞對的選擇,使其滿足基本的幾何接近。為了保證模擬結(jié)果的正確性,DSMC數(shù)值模擬中需要保證網(wǎng)格尺寸不超過當?shù)胤肿悠骄杂沙痰娜种弧8叱曀傧”×鲌鲋屑げ▍^(qū)域密度相對較大,分子平均自由程相對較小,而其它區(qū)域密度相對較低,分子平均自由程相對較大,在數(shù)值模擬過程中很難保證全流場滿足上述網(wǎng)格尺寸要求[10];同時,為了對密度相對較大的區(qū)域有較高的平均自由程和碰撞頻率的分辨率,流場時間步長必須取得非常小,這會降低流場整體計算效率[11]。

    為了放寬DSMC方法對網(wǎng)格尺寸的限制,文獻[7]在選擇分子碰撞對時引入了碰撞距離的思想,預(yù)先給定分子碰撞距離,以第一個分子的位置為圓心(或球心),以碰撞距離為半徑做圓(或球),在分子所在網(wǎng)格與圓(或球)的公共區(qū)域選擇第二個分子,該方法在每次分子碰撞對選擇時都需要進行公共區(qū)域判定及區(qū)域內(nèi)模擬分子標記,這必然會帶來較大的計算消耗。文獻[10]將自適應(yīng)碰撞距離的思想引入虛擬子網(wǎng)格方法,能夠避免上述處理帶來的計算消耗,但是該文獻采用常數(shù)作為分子碰撞對選取的最大循環(huán)次數(shù),常數(shù)值過小會對計算精度造成影響,而過大又會降低計算效率;并且,當循環(huán)結(jié)束仍未找到滿足要求分子碰撞對時,該文獻直接采用隨機選擇的分子對進行碰撞處理,這會降低較疏網(wǎng)格情況下的計算精度。本文發(fā)展了一種基于自適應(yīng)碰撞距離的分子碰撞對選取方法,能夠在程序中動態(tài)設(shè)定適合于每個網(wǎng)格單元的分子碰撞對最大循環(huán)次數(shù),同時對于循環(huán)結(jié)束仍不能找到滿足碰撞距離要求分子碰撞對的情況也做了針對性處理。

    為了加快流場的時間發(fā)展歷程,研究者們發(fā)展了自適應(yīng)當?shù)貢r間步長方法,通過構(gòu)造流場探測器,在流場計算中動態(tài)更新每個網(wǎng)格單元內(nèi)的時間乘數(shù)因子Ti,網(wǎng)格單元i內(nèi)的模擬分子在Ti倍基時間步長間隔內(nèi)只處理一次。本文結(jié)合自適應(yīng)當?shù)貢r間步長方法和DSMC方法數(shù)據(jù)結(jié)構(gòu)特點,發(fā)展了一種自適應(yīng)分子時間步長方法,能夠以模擬分子為最小時間步長調(diào)整單位,為其分配適當?shù)臅r間步長。

    1 計算方法

    1.1 DSMC方法介紹

    DSMC方法使用大量的模擬分子模擬真實的氣體,用一個模擬分子代表一定數(shù)量的真實氣體分子,整個模擬中分子之間以及分子和物面之間不斷進行能量交換,在充足的模擬時間后,通過采樣統(tǒng)計得到每個網(wǎng)格的宏觀流場結(jié)果[14-15]。本文采用非結(jié)構(gòu)貼體網(wǎng)格,計算中采用可變硬球(Variable Hard Sphere,VHS)分子模型,采用恒溫邊界條件,物面采用完全漫反射條件,使用Larsen- Borgnakke唯象論模型處理模擬分子間的能量交換。

    1.2 自適應(yīng)碰撞距離碰撞對選取方法

    為了在保證計算精度的同時降低DSMC方法對網(wǎng)格尺寸的限制,本文引入碰撞距離的思想,發(fā)展了一種自適應(yīng)碰撞距離的分子碰撞對選取方法,網(wǎng)格單元i中的自適應(yīng)碰撞距離di由下式確定:

    di=min(Hi,ref,λi/3)

    (1)

    其中,Hi,ref為網(wǎng)格單元i的特征尺寸,λi為單元i中的局部分子平均自由程,數(shù)值模擬過程中每隔400步對網(wǎng)格單元i中的λi和di進行一次更新,式(1)使碰撞距離限制在三分之一局部自由程之內(nèi)。

    分子碰撞對選取的主要處理步驟如下:

    1) 給定碰撞對中第二個分子選取的最大循環(huán)次數(shù)pi=max[1,int(NUMi/2)],其中NUMi為網(wǎng)格單元i中的模擬分子數(shù),集合{1,2,3...,NUMi}中的每個正整數(shù)對應(yīng)網(wǎng)格單元i中的一個模擬分子;

    2) 首先在集合{1,2,3...,NUMi}中隨機選取一個正整數(shù),并找到與之對應(yīng)的模擬分子M;

    3) 采用步驟2中的方法隨機選取第二個模擬分子L,計算并存儲M和L之間的距離dML;

    4) 判斷0

    上述處理步驟中采用與網(wǎng)格單元i中分子數(shù)相關(guān)的max[1,int(NUMi/2)]作為碰撞對中第二個分子選取的最大循環(huán)次數(shù),采用這種方式能夠獲得較高的碰撞對選取效率,引入max()函數(shù)是為了保證循環(huán)順利進行。步驟4中,達到pi次循環(huán)后仍未找到滿足需求的分子時,直接選取循環(huán)過程中距離M最近的分子與之構(gòu)成分子碰撞對,這樣處理可以在盡量滿足計算精度的同時保證計算效率。

    1.3 自適應(yīng)分子時間步長方法

    采用DSMC方法進行流場模擬時,需要記錄每個分子的性質(zhì)信息,包括速度、內(nèi)能、位置、所在網(wǎng)格單元以及分子組分等。經(jīng)過分析可以發(fā)現(xiàn),如果在DSMC模擬中額外存儲每一個分子的時間乘數(shù)因子信息,并參照自適應(yīng)當?shù)貢r間步長的思想,利用分子速度大小和與分子所在網(wǎng)格單元相關(guān)的特征長度構(gòu)造時間乘數(shù)因子,則可以發(fā)展出一種以模擬分子為最小時間步長調(diào)整單位的自適應(yīng)分子時間步長方法,主要思想如下:

    首先定義模擬分子所在網(wǎng)格單元的特征長度:

    Hk,locate,ref=min(Hi,ref,λi/3)

    (2)

    特征長度與1.2節(jié)中的自適應(yīng)碰撞距離大小一致,同樣隨流場宏觀參數(shù)變化而改變。這樣可以保證自適應(yīng)分子時間步長方法和自適應(yīng)碰撞距離碰撞對選取方法的兼容,且適用于較疏網(wǎng)格。

    調(diào)整編號為k分子的時間乘數(shù)因子Tk需要確定3個參數(shù):最大允許時間乘數(shù)因子N,各分子初始時間乘數(shù)因子Tk,0,以及基時間步長Δt。最大允許時間乘數(shù)因子由人為設(shè)定,初始時間乘數(shù)因子設(shè)定為1,本文基時間步長選取由下式?jīng)Q定:

    (3)

    其中,Href為網(wǎng)格單元的等效尺寸,V∞為來流速度,V∞為來流分子碰撞頻率。

    k編號分子的時間乘數(shù)因子Tk的調(diào)整方法為:

    1) 保存當前的時間乘數(shù)因子,Tk,old=Tk;

    2) 由分子速度大小Vk及分子所在單元特征長度Hk,locate,ref,確定其適合的時間乘數(shù)因子Tk,id= int(Hk,locate,ref/(VkΔt)),若Tk,id

    如果當前時間步Nk,t滿足mod(Nk,t,Tk)=0,則開始處理編號為k的分子,該分子下一次運動和碰撞的時間步長Δtp通??梢院唵蔚乇硎緸棣p=TkΔt,但由于分子可能在下一次被處理前已經(jīng)產(chǎn)生時間乘數(shù)因子的更新,此時Δtp≠TkΔt,因此,分子k下一次運動的時間步長可以表示為Δtp=[mod(Nk,t,Tk,old)+Tk,new-mod(Nk,t,Tk,new)]Δt,即為時間乘數(shù)因子調(diào)整以前剩余的運動時間步長mod(Nk,t,Tk,old)Δt和時間乘數(shù)因子調(diào)整以后的剩余時間步長[Tk,new-mod(Nk,t,Tk,new)]Δt之和。

    DSMC方法的本質(zhì)是在近似小的時間步長內(nèi),將分子的運動和分子間的碰撞進行解耦,采用自適應(yīng)分子時間步長方法時,同一個單元內(nèi)不同編號分子的時間乘數(shù)因子不盡相同,在Nk,t時間步、i單元中只有滿足mod(Nk,t,Tk)=0的分子才進行了運動模擬,因此,不能直接采用NTC法對i單元中的所有分子進行碰撞對取樣,需要作出適當調(diào)整,本文的做法是在每個基時間步都對i單元中滿足mod(Nk,t,Tk)=0的分子進行標記和數(shù)目統(tǒng)計,根據(jù)滿足上述條件的分子計算碰撞數(shù),并在其中進行分子碰撞對選取,通過上述方式保證單元中碰撞頻率與實際一致。

    2 算例驗證1:自適應(yīng)碰撞距離碰撞對選取

    采用不同網(wǎng)格尺寸下的圓柱外形對發(fā)展的自適應(yīng)碰撞距離分子碰撞對選取方法進行了計算分析。圓柱半徑R=0.08 m,來流氣體摩爾百分比為N2∶O2=0.763∶0.237,分子數(shù)密度為n=3.84×1020/m3,來流速度7500 m/s,來流溫度為Tinf=186 K,物面溫度為Twall=1000 K,采用統(tǒng)一時間步長。

    為了對自適應(yīng)碰撞距離碰撞對選取方法進行驗證,本文選取了兩種不同尺寸的非結(jié)構(gòu)網(wǎng)格進行了計算分析:第一種為較密的網(wǎng)格(物面附近網(wǎng)格尺度為來流分子平均自由程的六分之一,遠場網(wǎng)格尺度為來流分子平均自由程的三分之一),網(wǎng)格單元數(shù)為70890;第二種為較疏的網(wǎng)格,網(wǎng)格單元數(shù)為17710。對第二種網(wǎng)格分別采用自適應(yīng)碰撞距離方法與不采用該方法兩種工況進行了計算分析。

    圖1、圖2和圖3分別給出了兩種網(wǎng)格尺寸下物面熱流、壓力以及摩阻分布,從中可以看出,17710個單元的較疏網(wǎng)格在采用自適應(yīng)碰撞距離方法后能夠得到與70890個單元的較密網(wǎng)格符合較好的物面熱流、壓力以及摩阻分布結(jié)果,不采用自適應(yīng)碰撞距離方法的較疏網(wǎng)格與較密網(wǎng)格結(jié)果吻合得相對較差,圓柱外形算例表明本文發(fā)展的自適應(yīng)碰撞距離碰撞對選取方法能夠一定程度放寬DSMC方法對網(wǎng)格尺寸的限制,在相對較疏的網(wǎng)格下取得相對令人滿意的結(jié)果。

    圖1 不同網(wǎng)格下物面熱流密度分布Fig.1 Heat flux on the surface at different grid sizes

    圖2 不同網(wǎng)格下物面壓力分布Fig.2 Pressure on the surface at different grid sizes

    圖3 不同網(wǎng)格下物面摩阻分布Fig.3 Friction on the surface at different grid sizes

    3 算例驗證2:自適應(yīng)分子時間步長

    以算例2中的圓柱外形為例,選用網(wǎng)格單元數(shù)為17710的較疏網(wǎng)格,在采用自適應(yīng)碰撞距離碰撞對選取方法的基礎(chǔ)上,進行采用自適應(yīng)分子時間步長方法和采用統(tǒng)一時間步長的對比計算,基時間步長取為2×10-7s,計算條件與算例2保持一致,流場模擬分子總數(shù)趨于穩(wěn)定后約為28萬。對比計算是在一臺配置為Intel(R) Core(TM)2 2.80GHz,2.00GBRAM的微機上進行,運動和碰撞模擬共進行100000步,流場參數(shù)統(tǒng)計進行1000步。

    圖4給出了采用不同最大時間乘數(shù)因子每執(zhí)行100步運動和碰撞模擬所用時間的比較,N=1表示采用統(tǒng)一時間步長。從圖4可以看出,隨著N的增大,計算時間下降較快,當N>8時,計算時間趨于平穩(wěn),這主要是由于局部分子碰撞頻率控制了局部分子時間步長的自適應(yīng)過程,致使計算時間不能繼續(xù)下降,該趨勢與文獻[7]一致。本算例中最大時間乘數(shù)因子取N=10的計算時間為3.12 h,采用統(tǒng)一時間步長的計算時間為8.92 h,加速比為2.86。因此,采用本文發(fā)展的自適應(yīng)分子時間步長方法可以顯著縮短模擬所需的計算時間。

    圖4 每100步計算時間比較Fig.4 Comparison of computing time every 100 steps

    圖5為采用統(tǒng)一時間步長(N=1)的結(jié)果和N=10的流場平動溫度和壓力等值線對比圖,等值線圖按完全相同的參數(shù)繪制,從圖中可以看出,兩種方法下流場等值線幾乎完全一致。圖6為N=1和N=10的物面壓力和熱流密度的對比圖,同樣符合很好。圖5、圖6表明采用本文發(fā)展的自適應(yīng)分子時間步長方法對計算結(jié)果幾乎不產(chǎn)生影響。

    圖5 不同最大時間乘數(shù)因子對流場結(jié)構(gòu)影響的比較Fig.5 Comparison of field structure with different N

    (a) 壓力

    (b) 熱流密度

    4 結(jié) 論

    本文研究了DSMC計算中的高效處理方法,在碰撞對取樣環(huán)節(jié)引入碰撞距離的思想,發(fā)展了一種自適應(yīng)碰撞距離的分子碰撞對選取方法;在時間推進環(huán)節(jié)以模擬分子為最小時間步長調(diào)整單位,發(fā)展了一種自適應(yīng)分子時間步長方法。對圓柱外形驗證算例的計算分析表明在DSMC計算中采用本文發(fā)展的高效處理方法能夠在相對較疏的網(wǎng)格下得到與傳統(tǒng)DSMC計算相一致的結(jié)果。本文的方法能夠一定程度放寬DSMC方法對網(wǎng)格尺寸的限制,并且能夠大幅度縮短流場達到穩(wěn)定需要的計算時間。

    [1]沈青. 稀薄氣體動力學(xué)[M]. 北京: 國防工業(yè)出版社, 2003

    [2]樊菁. 稀薄氣體動力學(xué): 進展與應(yīng)用[J]. 力學(xué)進展, 2013, 43(2): 185- 201

    [3]Boyd I D. Computation of hypersonic flows using the direct simulation Monte Carlo method[C]//21st AIAA Computational Fluid Dynamics Conference, San Diego, CA, June 24-27, 2013

    [4]Qu C, Wang J F. Hybrid grid DSMC method for chemical nonequilibrium with rarefied flow heating[J]. Transactions of Nanjing

    University of Aeronautics and Astronautics, 2015, 32(4): 409-414

    [5]Moss J N, Glass C E, Greene F A. DSMC simulation of Apollo capsule aerodynamics for hypersonic rarefied conditions[R]. AIAA 2006-3577

    [6]Kim M G, Kim H S. A parallel cell based DSMC method with dynamic load balancing using unstructured adaptive meshes[R]. AIAA 2003-1033

    [7]王學(xué)德. 高超聲速稀薄氣流非結(jié)構(gòu)網(wǎng)格DSMC及并行算法研究[D]. 南京: 南京航空航天大學(xué), 2006

    [8]Matthias L. Optimization and parallelization of the DSMC method on unstructured grids[R]. AIAA 97-2512, 1997

    [9]Lebeau G J, Boyles K A, Lumpkin F E. Virtual sub-cells for the direct simulation Monte Carlo method[R]. AIAA 2003-1031

    [10]黃飛, 陳智, 程曉麗, 等. 一種基于自適應(yīng)碰撞距離的DSMC虛擬子網(wǎng)格方法[J]. 空氣動力學(xué)學(xué)報, 2014, 32(4): 506-510

    [11]Matthias L. Local time stepping with automatic adaptation for the DSMC method[R]. AIAA-98-2670, 1998

    [12]朱榮麗. 過渡區(qū)DSMC數(shù)值模擬的一種自適應(yīng)當?shù)貢r間步長新方法[J]. 宇航學(xué)報, 2008, 29(4): 1142-1146

    [13]Gallis M A, Torczynski J R, Rader D J, et al. Convergence behavior of a new DSMC algorithm[J]. Journal of Computational Physics, 2009, 228: 4532- 4548

    [14]Bird G A. Molecular gas dynamics and the direct simulation of gas flow[M]. Oxford: Clarendon Press, 1994

    [15]陳偉芳, 吳明巧, 趙玉新, 等. 超聲速平板繞流DSMC/EPSM混合算法研究[J]. 空氣動力學(xué)學(xué)報, 2002, 20(4): 441-445.

    猜你喜歡
    乘數(shù)步長流場
    基于Armijo搜索步長的BFGS與DFP擬牛頓法的比較研究
    大型空冷汽輪發(fā)電機轉(zhuǎn)子三維流場計算
    看錯了數(shù)字
    轉(zhuǎn)杯紡排雜區(qū)流場與排雜性能
    基于HYCOM的斯里蘭卡南部海域溫、鹽、流場統(tǒng)計分析
    理性認知西藏投資乘數(shù)小于1問題:以1996—2014年為例
    西藏研究(2016年4期)2016-06-05 11:31:15
    尋找突破角巧解算式謎
    基于瞬態(tài)流場計算的滑動軸承靜平衡位置求解
    基于逐維改進的自適應(yīng)步長布谷鳥搜索算法
    一種新型光伏系統(tǒng)MPPT變步長滯環(huán)比較P&O法
    電測與儀表(2014年2期)2014-04-04 09:04:00
    国产精品av久久久久免费| 久久人妻福利社区极品人妻图片 | 好男人电影高清在线观看| 欧美久久黑人一区二区| 国产黄频视频在线观看| 97在线人人人人妻| 午夜91福利影院| 国产精品 欧美亚洲| 男女之事视频高清在线观看 | 国产主播在线观看一区二区 | 性色av乱码一区二区三区2| 欧美激情 高清一区二区三区| 久久久久久久国产电影| 一级毛片女人18水好多 | 我要看黄色一级片免费的| 久久这里只有精品19| 久久中文字幕一级| 国产又爽黄色视频| 一二三四在线观看免费中文在| 国产熟女午夜一区二区三区| 人体艺术视频欧美日本| 97人妻天天添夜夜摸| 欧美中文综合在线视频| 国产精品 国内视频| 五月开心婷婷网| 美女午夜性视频免费| av不卡在线播放| 脱女人内裤的视频| 午夜激情久久久久久久| 欧美激情极品国产一区二区三区| 建设人人有责人人尽责人人享有的| 人人澡人人妻人| 啦啦啦在线观看免费高清www| 亚洲熟女精品中文字幕| 亚洲av在线观看美女高潮| 老司机深夜福利视频在线观看 | 国产欧美亚洲国产| 欧美黑人精品巨大| 亚洲图色成人| 久久精品亚洲熟妇少妇任你| 黄色a级毛片大全视频| 婷婷丁香在线五月| 久久狼人影院| 成年av动漫网址| 97人妻天天添夜夜摸| 国产人伦9x9x在线观看| 亚洲av成人精品一二三区| 免费高清在线观看视频在线观看| 欧美+亚洲+日韩+国产| 男男h啪啪无遮挡| 亚洲国产欧美在线一区| 99精品久久久久人妻精品| 丝袜美腿诱惑在线| 国产在线视频一区二区| 人人妻,人人澡人人爽秒播 | 97人妻天天添夜夜摸| 亚洲 欧美一区二区三区| xxx大片免费视频| 后天国语完整版免费观看| 成年美女黄网站色视频大全免费| 国产一卡二卡三卡精品| 97在线人人人人妻| 在线观看www视频免费| 又黄又粗又硬又大视频| 精品少妇久久久久久888优播| 国产高清videossex| 午夜老司机福利片| 大话2 男鬼变身卡| 亚洲精品美女久久av网站| 国产精品香港三级国产av潘金莲 | 亚洲精品自拍成人| 99精国产麻豆久久婷婷| 欧美人与善性xxx| 少妇裸体淫交视频免费看高清 | 伊人亚洲综合成人网| 天天添夜夜摸| 亚洲人成电影观看| 国产精品 国内视频| 97在线人人人人妻| 男男h啪啪无遮挡| 国产一区亚洲一区在线观看| 亚洲七黄色美女视频| 国产91精品成人一区二区三区 | av国产精品久久久久影院| 天天添夜夜摸| 国产真人三级小视频在线观看| 久久人人爽人人片av| 久久精品国产综合久久久| 看免费成人av毛片| 老司机靠b影院| 免费在线观看视频国产中文字幕亚洲 | 亚洲国产精品一区二区三区在线| 欧美 亚洲 国产 日韩一| av网站免费在线观看视频| www.999成人在线观看| 久久国产精品影院| 中文字幕高清在线视频| av网站在线播放免费| 成年美女黄网站色视频大全免费| 宅男免费午夜| 久久人人97超碰香蕉20202| 日日摸夜夜添夜夜爱| 久久久久国产一级毛片高清牌| 成人手机av| 免费在线观看完整版高清| 精品人妻1区二区| 亚洲av综合色区一区| 欧美日韩福利视频一区二区| 亚洲男人天堂网一区| 在线观看人妻少妇| 国产精品99久久99久久久不卡| 91麻豆av在线| 国产成人影院久久av| 美女高潮到喷水免费观看| 丝袜喷水一区| 精品国产超薄肉色丝袜足j| 丰满少妇做爰视频| 婷婷丁香在线五月| 亚洲专区中文字幕在线| 亚洲色图综合在线观看| 亚洲av成人精品一二三区| 高清黄色对白视频在线免费看| 国产欧美日韩一区二区三 | 一级毛片我不卡| 免费在线观看日本一区| 婷婷色麻豆天堂久久| 免费在线观看完整版高清| 日韩 欧美 亚洲 中文字幕| 国语对白做爰xxxⅹ性视频网站| √禁漫天堂资源中文www| 久久九九热精品免费| 国产成人精品久久二区二区91| 成年人免费黄色播放视频| 老司机靠b影院| 亚洲熟女毛片儿| 亚洲精品日本国产第一区| 中文精品一卡2卡3卡4更新| 亚洲av日韩精品久久久久久密 | 91精品伊人久久大香线蕉| 十分钟在线观看高清视频www| 日韩一本色道免费dvd| 两个人免费观看高清视频| 亚洲av综合色区一区| 侵犯人妻中文字幕一二三四区| 最近手机中文字幕大全| 久久久久精品国产欧美久久久 | 2018国产大陆天天弄谢| av网站免费在线观看视频| 国产一区二区激情短视频 | 国产精品99久久99久久久不卡| av电影中文网址| 又粗又硬又长又爽又黄的视频| 亚洲欧美一区二区三区黑人| 91麻豆精品激情在线观看国产 | 激情五月婷婷亚洲| 亚洲av成人精品一二三区| 久久天躁狠狠躁夜夜2o2o | 人体艺术视频欧美日本| 婷婷成人精品国产| 五月天丁香电影| 国产欧美日韩综合在线一区二区| 久久国产精品大桥未久av| 男人添女人高潮全过程视频| 国产片特级美女逼逼视频| 真人做人爱边吃奶动态| 多毛熟女@视频| 日韩制服骚丝袜av| 男的添女的下面高潮视频| 亚洲成色77777| 国产又爽黄色视频| av天堂在线播放| 水蜜桃什么品种好| 王馨瑶露胸无遮挡在线观看| 国产亚洲午夜精品一区二区久久| 国产极品粉嫩免费观看在线| 国产亚洲欧美在线一区二区| 亚洲精品美女久久久久99蜜臀 | 久久久久久久久免费视频了| 爱豆传媒免费全集在线观看| 十分钟在线观看高清视频www| 七月丁香在线播放| 大香蕉久久成人网| 中国美女看黄片| 少妇裸体淫交视频免费看高清 | 一边摸一边做爽爽视频免费| 国产成人啪精品午夜网站| 日韩制服骚丝袜av| 国产精品国产三级专区第一集| 中文字幕另类日韩欧美亚洲嫩草| 在线精品无人区一区二区三| www.自偷自拍.com| 狂野欧美激情性bbbbbb| 午夜福利影视在线免费观看| 国产免费一区二区三区四区乱码| 欧美国产精品一级二级三级| 精品久久蜜臀av无| 在线观看国产h片| 亚洲国产中文字幕在线视频| 国产麻豆69| 成年女人毛片免费观看观看9 | 国产熟女欧美一区二区| av欧美777| 亚洲国产成人一精品久久久| 啦啦啦啦在线视频资源| 老司机靠b影院| 天堂俺去俺来也www色官网| 日韩大片免费观看网站| 亚洲伊人久久精品综合| 精品国产乱码久久久久久小说| 人成视频在线观看免费观看| 国产有黄有色有爽视频| 黄色视频不卡| kizo精华| 欧美日本中文国产一区发布| 黄色视频不卡| 亚洲欧美清纯卡通| 免费一级毛片在线播放高清视频 | 国产高清videossex| 一区二区三区四区激情视频| 色播在线永久视频| 免费日韩欧美在线观看| 侵犯人妻中文字幕一二三四区| 精品国产超薄肉色丝袜足j| 国产成人a∨麻豆精品| 黑人欧美特级aaaaaa片| 一边亲一边摸免费视频| 亚洲成国产人片在线观看| 免费在线观看黄色视频的| 国产麻豆69| 国产熟女午夜一区二区三区| 国产一区二区在线观看av| 人人妻人人澡人人看| 最新在线观看一区二区三区 | 国产xxxxx性猛交| 亚洲欧美中文字幕日韩二区| 久久 成人 亚洲| 日日夜夜操网爽| 美女福利国产在线| 中文字幕av电影在线播放| 精品少妇久久久久久888优播| 一级毛片女人18水好多 | 看免费av毛片| 十分钟在线观看高清视频www| 国产成人av教育| 中文字幕色久视频| 欧美另类一区| 精品久久蜜臀av无| 中文字幕人妻熟女乱码| 欧美日韩成人在线一区二区| 免费高清在线观看日韩| 国产又色又爽无遮挡免| 深夜精品福利| 久久久国产精品麻豆| 亚洲人成77777在线视频| 欧美激情极品国产一区二区三区| 欧美激情极品国产一区二区三区| 尾随美女入室| 十八禁网站网址无遮挡| 精品国产乱码久久久久久男人| 欧美日韩亚洲综合一区二区三区_| 在线精品无人区一区二区三| 91成人精品电影| 久久精品久久久久久久性| 丝袜在线中文字幕| 中文字幕最新亚洲高清| 成人国语在线视频| 久久久久国产一级毛片高清牌| 首页视频小说图片口味搜索 | 精品福利永久在线观看| 久久精品人人爽人人爽视色| 国产成人精品久久久久久| 国产在线视频一区二区| 在线观看一区二区三区激情| 在线观看一区二区三区激情| 亚洲国产毛片av蜜桃av| 久久热在线av| av国产精品久久久久影院| 青青草视频在线视频观看| 久久久久久久久免费视频了| 久久青草综合色| 亚洲色图综合在线观看| 国产精品三级大全| 欧美日韩亚洲国产一区二区在线观看 | 国产男女超爽视频在线观看| 免费观看人在逋| 免费观看av网站的网址| 欧美日韩亚洲综合一区二区三区_| 精品国产一区二区三区久久久樱花| 大陆偷拍与自拍| 亚洲精品美女久久av网站| 十八禁网站网址无遮挡| 亚洲av电影在线进入| 汤姆久久久久久久影院中文字幕| 国产免费现黄频在线看| 在线观看人妻少妇| 男女无遮挡免费网站观看| kizo精华| 久久亚洲国产成人精品v| 久久亚洲精品不卡| 精品视频人人做人人爽| 悠悠久久av| 另类精品久久| 老熟女久久久| 精品亚洲成a人片在线观看| 少妇精品久久久久久久| 国产成人系列免费观看| 国产成人系列免费观看| 午夜福利一区二区在线看| 欧美黑人欧美精品刺激| 国产男女超爽视频在线观看| 色精品久久人妻99蜜桃| 亚洲欧美一区二区三区黑人| 中文字幕人妻丝袜制服| 亚洲九九香蕉| 国产成人精品无人区| 大话2 男鬼变身卡| 免费观看a级毛片全部| 久久久久久久久久久久大奶| 国产精品秋霞免费鲁丝片| 大码成人一级视频| 人人妻人人添人人爽欧美一区卜| 亚洲欧洲国产日韩| 国产精品欧美亚洲77777| 国产成人a∨麻豆精品| 国产高清视频在线播放一区 | 99热国产这里只有精品6| 亚洲av成人精品一二三区| 久久午夜综合久久蜜桃| 欧美乱码精品一区二区三区| www日本在线高清视频| 建设人人有责人人尽责人人享有的| 亚洲精品第二区| 亚洲精品国产av成人精品| 亚洲精品第二区| 晚上一个人看的免费电影| 啦啦啦啦在线视频资源| 韩国精品一区二区三区| 精品一品国产午夜福利视频| 国产精品九九99| 欧美人与性动交α欧美精品济南到| 欧美日韩亚洲国产一区二区在线观看 | 亚洲欧美日韩另类电影网站| 久久久久国产精品人妻一区二区| 国产无遮挡羞羞视频在线观看| 丁香六月天网| 国产午夜精品一二区理论片| 老汉色av国产亚洲站长工具| 高清视频免费观看一区二区| 99精国产麻豆久久婷婷| 国产成人精品久久久久久| 久久精品国产亚洲av高清一级| 免费不卡黄色视频| 久久国产精品人妻蜜桃| 亚洲国产最新在线播放| 侵犯人妻中文字幕一二三四区| 国产欧美日韩一区二区三区在线| 男女午夜视频在线观看| 国产在线一区二区三区精| 80岁老熟妇乱子伦牲交| 日韩,欧美,国产一区二区三区| 国产av精品麻豆| 视频在线观看一区二区三区| bbb黄色大片| 狂野欧美激情性xxxx| 成人18禁高潮啪啪吃奶动态图| 国语对白做爰xxxⅹ性视频网站| 麻豆乱淫一区二区| 91成人精品电影| 亚洲,欧美精品.| 女性生殖器流出的白浆| 免费看不卡的av| 日本黄色日本黄色录像| 国产成人精品久久二区二区91| 人体艺术视频欧美日本| 国产成人影院久久av| 亚洲五月婷婷丁香| 日本wwww免费看| 日本wwww免费看| 久久综合国产亚洲精品| 久久精品亚洲熟妇少妇任你| 亚洲专区中文字幕在线| www.av在线官网国产| 国产熟女午夜一区二区三区| 欧美日韩亚洲国产一区二区在线观看 | 精品一区在线观看国产| 热re99久久国产66热| 91字幕亚洲| 亚洲第一青青草原| 高清黄色对白视频在线免费看| 日韩伦理黄色片| 亚洲国产中文字幕在线视频| 婷婷色综合www| 每晚都被弄得嗷嗷叫到高潮| 久久国产精品男人的天堂亚洲| 免费看不卡的av| 在线看a的网站| 国产成人一区二区三区免费视频网站 | 精品国产国语对白av| 欧美性长视频在线观看| 午夜视频精品福利| 九草在线视频观看| 亚洲精品自拍成人| 黄片播放在线免费| 国产精品一区二区精品视频观看| 免费黄频网站在线观看国产| 啦啦啦 在线观看视频| 亚洲激情五月婷婷啪啪| 久久这里只有精品19| 国产成人a∨麻豆精品| 婷婷成人精品国产| 免费在线观看完整版高清| 日韩一本色道免费dvd| 久久久久视频综合| 国产精品免费视频内射| 91字幕亚洲| 国产精品二区激情视频| 中文欧美无线码| 午夜久久久在线观看| 亚洲欧美日韩高清在线视频 | 黑人欧美特级aaaaaa片| 老司机在亚洲福利影院| 女人被躁到高潮嗷嗷叫费观| avwww免费| 大话2 男鬼变身卡| 亚洲伊人久久精品综合| 中文字幕av电影在线播放| 一二三四在线观看免费中文在| 欧美黑人精品巨大| 国产成人精品在线电影| 国产精品成人在线| xxx大片免费视频| 精品国产国语对白av| 欧美亚洲日本最大视频资源| 极品少妇高潮喷水抽搐| 国产精品九九99| 久久99精品国语久久久| 午夜影院在线不卡| 精品人妻熟女毛片av久久网站| 99精品久久久久人妻精品| 中文字幕高清在线视频| 精品国产超薄肉色丝袜足j| 午夜福利乱码中文字幕| 久久99一区二区三区| 一区二区日韩欧美中文字幕| 满18在线观看网站| 国产成人免费观看mmmm| 超碰成人久久| 国产视频一区二区在线看| 国产午夜精品一二区理论片| 亚洲一区二区三区欧美精品| 菩萨蛮人人尽说江南好唐韦庄| 黄色怎么调成土黄色| 最近手机中文字幕大全| 国产成人欧美| 国产精品国产av在线观看| 99精国产麻豆久久婷婷| 极品人妻少妇av视频| 91精品国产国语对白视频| 亚洲专区国产一区二区| 国产真人三级小视频在线观看| 国产成人精品久久二区二区91| 别揉我奶头~嗯~啊~动态视频 | 亚洲国产精品999| 国产一卡二卡三卡精品| 欧美xxⅹ黑人| 一区二区av电影网| 天天躁日日躁夜夜躁夜夜| 91麻豆精品激情在线观看国产 | 无遮挡黄片免费观看| 老汉色av国产亚洲站长工具| 亚洲av日韩精品久久久久久密 | 欧美日韩亚洲综合一区二区三区_| 99热网站在线观看| 五月天丁香电影| 国产精品久久久av美女十八| 日韩,欧美,国产一区二区三区| 免费高清在线观看视频在线观看| 国产精品亚洲av一区麻豆| 精品卡一卡二卡四卡免费| 国产亚洲精品第一综合不卡| 天堂中文最新版在线下载| 久久久国产一区二区| 亚洲自偷自拍图片 自拍| 嫩草影视91久久| 91麻豆精品激情在线观看国产 | 国产精品偷伦视频观看了| 搡老乐熟女国产| 一级黄片播放器| 免费少妇av软件| 国产熟女欧美一区二区| 可以免费在线观看a视频的电影网站| 热99国产精品久久久久久7| 国产又色又爽无遮挡免| 色精品久久人妻99蜜桃| 黄色毛片三级朝国网站| 国产高清视频在线播放一区 | av天堂久久9| 妹子高潮喷水视频| 久久久精品免费免费高清| 亚洲一码二码三码区别大吗| 国产一区亚洲一区在线观看| 99国产精品一区二区三区| 久久人人97超碰香蕉20202| 色网站视频免费| 国产精品成人在线| 色视频在线一区二区三区| 亚洲精品第二区| 亚洲精品美女久久av网站| 免费在线观看黄色视频的| 啦啦啦中文免费视频观看日本| 两个人看的免费小视频| 黄色 视频免费看| 国产精品二区激情视频| 美女扒开内裤让男人捅视频| 超碰成人久久| 80岁老熟妇乱子伦牲交| 最黄视频免费看| 免费久久久久久久精品成人欧美视频| 午夜激情av网站| 男女无遮挡免费网站观看| 久久人妻熟女aⅴ| 婷婷色综合大香蕉| 无限看片的www在线观看| 精品亚洲乱码少妇综合久久| 夫妻性生交免费视频一级片| 国产精品亚洲av一区麻豆| 国产视频一区二区在线看| 99国产精品免费福利视频| 99热国产这里只有精品6| 日日爽夜夜爽网站| 成年av动漫网址| 亚洲欧美精品自产自拍| 男女国产视频网站| 男人爽女人下面视频在线观看| 日本91视频免费播放| 国产成人精品久久久久久| 亚洲欧洲国产日韩| 无限看片的www在线观看| 老司机影院成人| 亚洲,欧美精品.| 久久99精品国语久久久| 丝袜美足系列| 韩国精品一区二区三区| 免费高清在线观看日韩| 大码成人一级视频| 丰满少妇做爰视频| 搡老岳熟女国产| 久久天躁狠狠躁夜夜2o2o | 免费久久久久久久精品成人欧美视频| 老司机深夜福利视频在线观看 | 国产一区二区在线观看av| 手机成人av网站| 黄片小视频在线播放| 免费不卡黄色视频| 女性被躁到高潮视频| 美女午夜性视频免费| 欧美激情 高清一区二区三区| 日韩大码丰满熟妇| 91精品国产国语对白视频| 久久亚洲国产成人精品v| 国产97色在线日韩免费| 夜夜骑夜夜射夜夜干| 日韩电影二区| 亚洲人成77777在线视频| 国产国语露脸激情在线看| 91精品国产国语对白视频| 亚洲专区中文字幕在线| 免费看十八禁软件| 久久精品国产亚洲av涩爱| 国产欧美日韩综合在线一区二区| 制服人妻中文乱码| 在线观看www视频免费| 老司机影院毛片| 欧美国产精品一级二级三级| 韩国高清视频一区二区三区| 精品少妇黑人巨大在线播放| 日本vs欧美在线观看视频| 亚洲精品国产av蜜桃| 99热国产这里只有精品6| 午夜日韩欧美国产| 亚洲专区中文字幕在线| 国产熟女午夜一区二区三区| 性色av一级| 亚洲五月婷婷丁香| 国产日韩欧美视频二区| 国产成人精品久久久久久| 精品国产国语对白av| 国产成人精品在线电影| 国产成人a∨麻豆精品| 女性生殖器流出的白浆| 丁香六月欧美| 国产精品.久久久| 精品第一国产精品| 欧美黑人精品巨大| 国产亚洲欧美精品永久| 一个人免费看片子| 精品一区二区三区四区五区乱码 | 好男人电影高清在线观看| 久9热在线精品视频| 午夜两性在线视频| 免费高清在线观看视频在线观看| 又紧又爽又黄一区二区| 亚洲人成电影观看| 午夜两性在线视频| 国产黄频视频在线观看| 婷婷成人精品国产| 69精品国产乱码久久久| 久久久久久久国产电影| 99精品久久久久人妻精品| 日韩制服丝袜自拍偷拍| 少妇人妻 视频| 18禁观看日本| 久久久精品国产亚洲av高清涩受| 伊人久久大香线蕉亚洲五| 久久久精品国产亚洲av高清涩受|