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

    空間自適應(yīng)的快速粒子Level Set方法

    2010-05-10 06:26:26黃筱云李紹武
    關(guān)鍵詞:交界面高分辨率算例

    黃筱云,李紹武

    (天津大學(xué)建筑工程學(xué)院,天津 300072)

    空間自適應(yīng)的快速粒子Level Set方法

    黃筱云,李紹武

    (天津大學(xué)建筑工程學(xué)院,天津 300072)

    Level set作為一種有效的界面追蹤方法,已經(jīng)被廣泛地應(yīng)用到單相流的自由表面追蹤或多相流的交界面追蹤.快速粒子level set方法(FPLS)是一種改進的level set方法,既具有較高的界面追蹤精度,又有較高的計算效率.將其與自適應(yīng)網(wǎng)格結(jié)合起來,提出了一種基于樹型結(jié)構(gòu)空間自適應(yīng)的快速粒子 level set方法(SA-FPLS).通過在交界面附近設(shè)置高分辨率網(wǎng)格帶,且只在高分辨率網(wǎng)格帶內(nèi)求解level set方程,既進一步提高了交界面模擬精度,又可以節(jié)省計算機內(nèi)存,同時保持較高計算效率.

    快速粒子level set;自適應(yīng)網(wǎng)格;空間自適應(yīng)的快速粒子level set

    在用數(shù)值方法求解帶有交界面的流體力學(xué)問題時,必須通過一定方法預(yù)先給出流體交界面的位置.level set方法是近年發(fā)展起來的一種界面追蹤方法[1],其原理是將交界面看作滿足以下控制方程的符號距離函數(shù)φ的零等值線(面).

    式中:φ稱為level set函數(shù),其值在交界面上為零,離開交界面為正或為負;ui為速度張量.

    為使交界面的解在出現(xiàn)界面融合或較大φ值梯度時仍有較高的精度,一般采用重新初始化的方法[2],即在每完成一步level set方程求解過程后,通過迭代求解方程

    的解,來對level set函數(shù)進行初始化,其中0φ為level set方程的計算結(jié)果,xΔ為網(wǎng)格尺寸.這種level set方法在整個計算域上進行計算,稱之為全局 level set方法.

    為了減少計算量,Peng等[3]提出了一種局部的level set方法,它只需在交界面附近區(qū)域求解level set方程和重新初始化.事實上,level set方法中采用迭代法進行重新初始化是影響計算速度的關(guān)鍵因素,Sethian[4]提出了一種fast marching技術(shù)求解level set方程,其特點是無需重新初始化,但計算精度較低.

    與VOF方法相比較,level set方法的缺點是守恒性差.為了彌補這一缺陷,Sussman等[5]將 VOF與level set方法相結(jié)合,提出了混合型的界面追蹤方法——CLSVOF,該方法利用 VOF函數(shù)來校正 level set函數(shù)結(jié)果以達到減少質(zhì)量損失的目的,但該方法并不能完全消除交界面附近的噪聲,追蹤到的交界面與實際值存在相位差.另外一種改進的level set方法是粒子level set(PLS)方法[6],它將MAC法的特性引入到 level set方法中,計算時通過粒子來修正 level set函數(shù),它捕捉到的交界面十分光滑,幾乎無質(zhì)量損失.為了加快計算速度,Enright等[7]提出了一種局部快速粒子level set方法——FPLS,該方法采用一階精度的半Lagrangian格式[8]求解level set方程,用一階精度fast marching技術(shù)[9]進行重新初始化,粒子的運算則采用二階精度Runge-Kutta格式.

    在實際應(yīng)用中,復(fù)雜交界面需要高分辨率網(wǎng)格來模擬,采用均勻網(wǎng)格將消耗大量計算機內(nèi)存,計算速度也將大大降低.Strain[10]在采用半 Lagrangian方法求解level set方程時,首先引入了樹型結(jié)構(gòu)的自適應(yīng)網(wǎng)格,這種網(wǎng)格結(jié)構(gòu)在交界面附近采用高分辨率網(wǎng)格,而在遠離交界面的區(qū)域則保持較粗的分辨率.這樣一方面提高了計算結(jié)果的精度,另一方面減少了無關(guān)區(qū)域所占用的內(nèi)存.但該方法為了保持在 T型節(jié)點插值的連續(xù)性,重新初始化階段在交界面附近區(qū)域采用三角網(wǎng)格[11],計算量仍較大.2002年,Sochnikov等[12]將 fast marching技術(shù)應(yīng)用到四叉樹型結(jié)構(gòu)的自適應(yīng)網(wǎng)格中,使得計算量進一步減少.這種自適應(yīng)level set方法是通過減小交界面附近網(wǎng)格尺寸來保證解的守恒性,可能會導(dǎo)致交界面周圍的網(wǎng)格與最粗網(wǎng)格的尺寸相差很大,形成了過于復(fù)雜的樹型結(jié)構(gòu).

    筆者將樹型結(jié)構(gòu)的自適應(yīng)性與 FPLS方法結(jié)合起來,提出一種改進的自適應(yīng) level set方法——SAFPLS,該方法在交界面附近能自動形成一個具有固定寬度的高分辨率網(wǎng)格帶,并在該網(wǎng)格帶內(nèi)求解level set方程,進行重新初始化運算,且粒子始終包含在最小網(wǎng)格帶內(nèi).該方法保留了自適應(yīng)level set方法的優(yōu)點,同時 FPLS方法的高精度性使得該方法無需在交界面周圍采用過細的網(wǎng)格,這樣在保證計算結(jié)果高精度的同時,計算速度更快,消耗內(nèi)存更少.SAFPLS方法的具體步驟為:

    (1)按照初始化 level set函數(shù),生成在交界面附近加密的樹型結(jié)構(gòu)網(wǎng)格;

    (2) 采用FPLS方法求解level set方程;

    (3) 根據(jù)level set函數(shù)值重新調(diào)整網(wǎng)格;

    (4) 返回第(2)步.

    1 FPLS方法的基本原理

    在 FPLS方法中,正負粒子被隨機布置在交界面周圍窄帶中(初始時刻正值區(qū)的粒子為正粒子;反之為負粒子).每個粒子有各自的半徑 rp,大小被限定在 最 小 值 rmin= 0 .1max ( Δ x, Δy )和 最 大 值 rmax=0.5max(Δ x,Δy)之間,其值可根據(jù)下列公式確定.

    式中sp為粒子的符號(正粒子為+1,負粒子為-1).粒子的φ值將通過雙線性或三線性插值得到.

    這種確定半徑辦法的目的是使粒子盡量與交界面相切.當(dāng)粒子穿過交界面達一個半徑時,即認為是逃逸粒子.在劇烈變化的流場中,通過 level set方程得到的交界面附近的 level set函數(shù)值可能存在較大的誤差.為了消除誤差,可借助逃逸粒子來對 level set函數(shù)值進行修正.為此,定義一個與逃逸粒子有關(guān)的局部level set函數(shù)pφ,即

    如果函數(shù)pφ和φ值不相等,說明level set方程得出的φ函數(shù)可能存在誤差,取兩者中絕對值較小的作為校正后的φ值.這種對φ的校正在重新初始化的過程中同樣也需要.

    當(dāng)界面發(fā)生拉伸和斷裂時,交界面附近網(wǎng)格粒子變得稀疏,需要補充粒子以確保追蹤界面的準確性.而當(dāng)非逃逸粒子遠離交界面時,則需要將這些粒子刪除.

    在FPLS中,采用了一種低階的半Lagrangian算法求解level set方程,該格式無條件穩(wěn)定且收斂于真解.公式為

    在進行重新初始化過程中,式(3)采用以下一階迎風(fēng)差分格式進行離散.

    通過求解式(7),式(3)的解可以快速地構(gòu)造出來,這就是fast marching技術(shù)的主要思想,即從界面附近區(qū)域出發(fā),逐步構(gòu)造出Eikonal方程的解(見圖1).

    圖1 采用快速步進法進行重新初始化Fig.1 Reinitialization with fast marching technique

    2 SA-FPLS的原理

    為簡單起見,這里僅以二維的四叉樹網(wǎng)格說明SA-FPLS的原理.三維的八叉樹網(wǎng)格與之類似,不再贅述.在四叉樹網(wǎng)格結(jié)構(gòu)中,每個四叉樹網(wǎng)格定義4個指向其子網(wǎng)格的指針(如果子網(wǎng)格不存在,指針指向空)和一個指向其父網(wǎng)格(如果父網(wǎng)格不存在,指針指向空)的指針,每個網(wǎng)格另有4個指向角點的指針和4個指向邊的指針.每個節(jié)點有4個指針指向其相鄰葉網(wǎng)格,對于特殊的節(jié)點,如T型節(jié)點,4個指針中有2個指向同一個網(wǎng)格.每個節(jié)點需要保存該節(jié)點的坐標(biāo)及定義在該點的物理量.最粗的網(wǎng)格為 0層網(wǎng)格,略粗的為第1層,以此類推.圖2中網(wǎng)格Ⅰ有4個子網(wǎng)格Ⅱ、Ⅲ、Ⅳ和Ⅴ和 4個節(jié)點 A、B、C 和 D,與子網(wǎng)格Ⅳ有共同節(jié)點C.網(wǎng)格Ⅳ有 4個子網(wǎng)格Ⅵ、Ⅶ、Ⅷ和Ⅸ和1個父網(wǎng)格Ⅰ.網(wǎng)格Ⅱ與網(wǎng)格Ⅲ共享1號邊.節(jié)點E有4個指針分別指向網(wǎng)格Ⅱ、Ⅲ、Ⅳ和Ⅴ.在這種四叉樹網(wǎng)格結(jié)構(gòu)中,網(wǎng)格在每個方位上都有唯一的鄰居,定義為與其相鄰且層數(shù)不大于該網(wǎng)格的網(wǎng)格中層數(shù)最大的一個.如圖3中網(wǎng)格Ⅱ與網(wǎng)格Ⅰ、Ⅳ、Ⅴ和Ⅵ相鄰,但只有網(wǎng)格Ⅳ才是網(wǎng)格Ⅱ的鄰居.

    圖2 樹型網(wǎng)格結(jié)構(gòu)Fig.2 Tree grid structure

    圖3 網(wǎng)格及其鄰居Fig.3 Given grid and its neighboring grid

    在生成樹型結(jié)構(gòu)網(wǎng)格的過程中,首先分裂包含交界面的網(wǎng)格,若一個網(wǎng)格的角點上φ的符號相反,交界面必定穿越該網(wǎng)格,該網(wǎng)格即需要分裂,直至網(wǎng)格的層數(shù)達到規(guī)定的最大層數(shù).因此,交界面附近的網(wǎng)格分辨率是最高的,那些與高分辨率網(wǎng)格相鄰的網(wǎng)格也需要進行分裂,以保證高分辨率網(wǎng)格具有一定帶寬.在交界面周圍生成高分辨率的網(wǎng)格之后,還需要對網(wǎng)格進行平衡處理[13],即保持相鄰網(wǎng)格的大小相差不會超過 1倍.新產(chǎn)生的網(wǎng)格節(jié)點上的物理量通過線性插值得到.在計算過程中,為了提高交界面追蹤精度,交界面附近的網(wǎng)格應(yīng)始終保持高分辨率,且加密網(wǎng)格隨著交界面的變化而變化.這樣可以保證只在網(wǎng)格帶內(nèi)求解level set方程和進行φ值校正,同時防止因在交界面附近分裂網(wǎng)格時,線性插值導(dǎo)致的較大計算誤差.因此,隨著交界面的運動,需要在交界面附近將預(yù)先設(shè)定的高分辨率帶的網(wǎng)格不斷進行分裂,同時也要將脫離高分辨率網(wǎng)格帶的網(wǎng)格進行合并(見圖 4).

    圖4 交界面附近的高分辨率網(wǎng)格帶Fig.4 Grid band of highest resolution around interface

    在四叉樹網(wǎng)格中采用半 Lagrangian方法求解level set方程與在均勻網(wǎng)格中有所不同,網(wǎng)格節(jié)點在上個時刻投影點的位置需要根據(jù)四叉樹網(wǎng)格的結(jié)構(gòu)尋找.首先,尋找投影點所在的 0層網(wǎng)格;然后,從該網(wǎng)格逐步向下找到投影點所在的子網(wǎng)格;最后,根據(jù)所在子網(wǎng)格角點上前一時刻的φ值,線性插值得到該節(jié)點新時刻的φ值.

    在重新初始化的過程中,由于采用 fast marching技術(shù),φ值是從交界面出發(fā),迎風(fēng)地向外構(gòu)造得到,因此不會影響交界面的位置,當(dāng)計算的某個節(jié)點的φ值大于高分辨率帶的寬度時,向外構(gòu)造過程即停止.值得注意的是,在向外構(gòu)造的過程中,若遇到 T型節(jié)點時,則需要分裂該點周圍所有層數(shù)小于規(guī)定最大層數(shù)的網(wǎng)格,以保證高精度網(wǎng)格帶的寬度.

    3 算例及分析

    在二維算例中,采用如下一階近似算法計算交界面的長度L,即

    3.1 算例1

    本算例模擬星型結(jié)構(gòu)的擴張過程.交界面初始形狀如圖 5(a)所示,交界面以單位速度向向外法向擴張.計算區(qū)域大小為[0 ,2]×[0 ,2],網(wǎng)格層數(shù)為4層,0層網(wǎng)格尺寸為0.04×0.04.初始level set函數(shù)表示為

    式中:Ra為平均半徑;d為偏離Ra的距離;n為觸角的個數(shù);θ為半徑與 x正向的夾角.在本算例中,Ra= 0 .3, d = 0 .2, n= 7 .圖 5(b)是交界面分別在 0~0.05 s 6個時刻的界面位置.可以看出,界面在運動中始終保持與初始形狀特征相似的形狀,只是觸角變得越來越厚.

    圖5 七角星結(jié)構(gòu)擴張過程的數(shù)值模擬結(jié)果Fig.5 Numerical results of expansion of seven-pointed star

    3.2 算例2

    本算例目的在于檢驗 SA-FPLS方法在剪切流場中對復(fù)雜交界面追蹤的有效性.計算區(qū)域為[0 ,1]×[0 ,1],網(wǎng)格層數(shù)為4層,0層網(wǎng)格尺寸為0.04×0.04,界面初始形狀為一個半徑為0.15、圓心在[0 .5,0.75]處的圓(見圖6(a)).速度場表示為

    式中,流場周期T為8 s,時間步長為0.005 s.圓的變形過程如圖 6所示.可見,由于在交界面附近采用了高分辨率的網(wǎng)格,尾部絲狀結(jié)構(gòu)保持得相當(dāng)完整.同時,由于采用自適應(yīng)網(wǎng)格,大部分區(qū)域保持較粗的網(wǎng)格,從而節(jié)約了大量內(nèi)存和計算時間.

    3.3 算例3

    本算例模擬圓在1個有16個渦旋的流場中的變形情況.圓的初始狀態(tài)及初始網(wǎng)格結(jié)構(gòu)與算例 2一致,不同的是,該模型采用周期性邊界條件,其周期性速度場[15]表示為

    式中,流場周期T為2 s,計算時間步長為0.005 s.變形過程如圖 7所示,在變形過程中,由于采用了周期性邊界條件,交界面穿過計算區(qū)域上端后,又重新出現(xiàn)在計算區(qū)域下端,且上下邊界的網(wǎng)格與節(jié)點保持一一對應(yīng).可以看出 FPLS與自適應(yīng)網(wǎng)格相結(jié)合,可以適應(yīng)十分復(fù)雜的流場結(jié)構(gòu).

    圖6 圓在簡單旋轉(zhuǎn)流場中變形過程的數(shù)值模擬結(jié)果Fig.6 Numerical results of transformation of circle in single vortex flow

    圖7 圓在多渦流場中變形情況的數(shù)值模擬Fig.7 Numerical results of transformation of circle in multi-vortex flow

    檢驗?zāi)P褪睾阈詢?yōu)劣的方法是對比初始時刻界面所圍圓面積與 1個周期后界面恢復(fù)到原始形狀后面積的差異.界面所圍面積的計算式為

    式中:L為交界面的期望長度;H為Heaviside函數(shù);φe為真解;φc為計算結(jié)果.算例2和算例3計算結(jié)果的誤差分析如表1所示.

    3.4 算例4

    算例4和算例5將討論SA-FPLS方法在三維空間中的應(yīng)用.算例4模擬一個半徑為0.15的Zalesak圓球[16]在流場作用下的變化情況,球的初始位置為[0 .5,0.75,0.5],切口的寬度為 0.05,深度為 0.125.計算區(qū)域尺度為[0 ,1] × [ 0 ,1] × [ 0 ,1],網(wǎng)格層數(shù)為 4,最粗網(wǎng)格大小為0.0625× 0 .0625× 0 .0625,流場表示為

    表1 計算結(jié)果的誤差分析Tab.1 Error analysis of numerical results

    計算時間為 2 s,時間步長為 0.005 s.界面運動過程如圖 8所示,Zalesak圓球經(jīng)過一個周期的運動后又回到初始位置,而缺口仍然保持尖銳.

    圖8 Zalesak圓球的旋轉(zhuǎn)過程數(shù)值模擬結(jié)果Fig.8 Numerical results of rotation of Zalesak sphere

    3.5 算例5

    本算例模擬一個半徑為 0.15、初始位置為[0.35,0.35,0.35]的球體在三維流場作用下的變形過程,如圖 9所示.計算區(qū)域為[0 , 1] × [ 0 , 1] × [ 0 , 1],網(wǎng)格層數(shù)為 3,最粗網(wǎng)格大小為0.04× 0 .04× 0 .04.流場表示為式中,流場周期 T為 2,s,計算時間步長為 0.005,s.結(jié)果表明,圓球在1.0,s時達到最大拉伸后,又恢復(fù)到原樣,而其中各圖為圓球分別在 t為0~2.0,s時的形狀.

    4 結(jié) 語

    將自適應(yīng)網(wǎng)格與 FPLS結(jié)合起來,探討了一種交界面追蹤方法 SA-FPLS.從算例中可以看出,SAFPLS一方面繼承了 FPLS的優(yōu)點,交界面具有良好的光滑性,總體質(zhì)量有良好的守恒性;另一方面由于它只在交界面附近窄帶內(nèi)求解level set函數(shù),包括進行重新初始化,同時基于四叉樹或八叉樹的自適應(yīng)網(wǎng)格只在交界面附近加密網(wǎng)格,而在遠離交界面的區(qū)域保持較粗網(wǎng)格,因而這種方法既可節(jié)約內(nèi)存,又可提高計算效率.

    [1]Osher S,Sethian J. Fronts propagating with curvature dependent speed:Algorithms based on Hamiliton-Jacobi formulations[J].Journal of Computational Physics,1988,79(1):12-49.

    [2]Sussman M,Puckett E G,Osher S. A level set approach for computing solutions to incompressible two-phase flow[J].Journal of Computational Physics,1994,114(2):146-159.

    [3]Peng D,Merriman B,Osher S,et al. A PDE-based fast local level set method[J].Journal of Computational Physics,1999,155(2):410-438.

    [4]Sethian J. A fast marching level set method for monotonically advancing fronts[J].Applied Mathematics,1996,93(4):1591-1595.

    [5]Sussman M,Puckett E G. A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows[J].Journal of Computational Physics,2000,162(2):301-337.

    [6]Enright D,F(xiàn)edkiw R. A hybrid particle level set method for improved interface capturing[J].Journal of Computational Physics,2002,183(1):83-116.

    [7]Enright D,Losassom F. A fast and accurate semi-Lagrangian particle level set method[J].Computers and Structure,2005,83(6):479-490.

    [8]Strain J. Semi-Lagrangian methods for level set equations[J].Journal of Computational Physics,1999,151(2):498-533.

    [9]Sethian J. Fast marching methods[J].SIAM Review,1999,41(2):199-235.

    [10]Strain J. Tree methods for moving interfaces [J].Journal of Computational Physics,1999,151(2):616-648.

    [11]Strain J. Fast tree-based redistancing for level set computations[J].Journal of Computational Physics,1999,152(2):664-686.

    [12]Sochnikov V,Efrim S. Level set calculations of the evolution of boundaries on a dynamically adaptive grid[J].International Journal for Numerical Methods in Engineering,2002,56(13):1913-1929.

    [13]Mark de B,Mark van K,Mark O,et al.Computational Geometry Algorithms and Applications[M]. Berlin:Springer,1999.

    [14]Sussman M,F(xiàn)atemi E. An efficient,interface-preserving level set redistancing algorithm and its application to interfacial incompressible fluid flow[J].SIAM Journal on Scientific Computing,1999,20(6):1165-1191.

    [15]Smolarkiewicz P. The multi-dimensional crowley advection scheme[J].Monthly Weather Review,1982,110(12):1968-1983.

    [16]Zalesak S. Fully multidimensional flux-corrected transport algorithms for fluids[J].Journal of Computational Physics,1979,31(3):335-362.

    Spatially Adaptive Fast Particle Level Set Method

    HUANG Xiao-yun,LI Shao-wu
    (School of Civil Engineering,Tianjin University,Tianjin 30072,China)

    As an effective method of interface tracking,the level set method has been widely applied to surface tracking of single phase flow and interface tracking of multi-phase flow. The fast particle level set(FPLS)method,as an improved version of the level set method,not only is of relatively high accuracy in surface tracking,but also possesses high efficiency of calculation. By combination of the FPLS with tree-based adaptive gridding system,a spatially adaptive fast particle level set(SA-FPLS)approach is proposed,By solving the level-set equation with higher gridding resolution in a narrow mesh band around the interface,the solution accuracy can be further improved with less computer memory consumption and higher computational efficiency.

    fast particle level set;adaptive grid;spatially adaptive fast particle level set

    TV131.2

    A

    0493-2137(2010)11-0981-07

    2009-09-02;

    2010-05-14.

    黃筱云(1980— ),男,博士研究生,td98298@sina.com.

    李紹武,lishaowu@tju.edu.cn.

    猜你喜歡
    交界面高分辨率算例
    鋼-混凝土交界面法向粘結(jié)性能研究
    高速公路機電工程相關(guān)交界面管理組織建設(shè)探討
    高分辨率合成孔徑雷達圖像解譯系統(tǒng)
    雙塊式無砟軌道軌枕與道床交界面損傷特性分析
    中國鐵路(2019年1期)2019-03-23 01:11:58
    基于振蕩能量的低頻振蕩分析與振蕩源定位(二)振蕩源定位方法與算例
    互補問題算例分析
    高分辨率對地觀測系統(tǒng)
    太空探索(2015年8期)2015-07-18 11:04:44
    基于Curvelet-Wavelet變換高分辨率遙感圖像降噪
    基于CYMDIST的配電網(wǎng)運行優(yōu)化技術(shù)及算例分析
    高分辨率遙感相機CCD器件精密熱控制
    亚洲午夜理论影院| а√天堂www在线а√下载| 两性夫妻黄色片| 国产成人aa在线观看| 久久久久精品国产欧美久久久| 91九色精品人成在线观看| 亚洲成人免费电影在线观看| 久久精品国产清高在天天线| 午夜免费激情av| 亚洲avbb在线观看| 婷婷亚洲欧美| 久久精品夜夜夜夜夜久久蜜豆| 变态另类丝袜制服| 午夜激情福利司机影院| 99精品在免费线老司机午夜| 天天一区二区日本电影三级| 亚洲熟女毛片儿| 12—13女人毛片做爰片一| 久久人人精品亚洲av| 国产亚洲欧美在线一区二区| 久久久国产成人免费| 亚洲av成人精品一区久久| 欧美高清成人免费视频www| 一二三四在线观看免费中文在| 国内精品一区二区在线观看| 最新中文字幕久久久久 | 亚洲精品国产精品久久久不卡| 国产精品 国内视频| 亚洲欧美精品综合久久99| 国产熟女xx| 男人舔奶头视频| 国产一区二区激情短视频| 午夜精品在线福利| 国产综合懂色| 午夜视频精品福利| 亚洲精品粉嫩美女一区| 久久这里只有精品19| 久久久精品大字幕| 精品电影一区二区在线| 国产一区二区三区在线臀色熟女| 99在线视频只有这里精品首页| 丰满人妻熟妇乱又伦精品不卡| 两性午夜刺激爽爽歪歪视频在线观看| 午夜视频精品福利| 久久精品国产亚洲av香蕉五月| 极品教师在线免费播放| 人妻久久中文字幕网| 9191精品国产免费久久| 婷婷精品国产亚洲av在线| 久久久久久久午夜电影| 别揉我奶头~嗯~啊~动态视频| 国产av一区在线观看免费| 成人高潮视频无遮挡免费网站| 久久久久亚洲av毛片大全| 90打野战视频偷拍视频| 一级作爱视频免费观看| 久久久久久久久中文| 18禁观看日本| 久久久精品欧美日韩精品| 一进一出抽搐动态| 久久精品亚洲精品国产色婷小说| 日韩精品青青久久久久久| 日韩三级视频一区二区三区| 18禁黄网站禁片午夜丰满| 亚洲欧洲精品一区二区精品久久久| 757午夜福利合集在线观看| 一级a爱片免费观看的视频| av在线天堂中文字幕| 美女 人体艺术 gogo| 一级毛片精品| 波多野结衣高清无吗| 欧美日韩中文字幕国产精品一区二区三区| 美女cb高潮喷水在线观看 | 国产三级中文精品| 久久久久久大精品| 最新在线观看一区二区三区| 岛国在线观看网站| 国产伦精品一区二区三区视频9 | 国产成人精品久久二区二区91| 午夜激情欧美在线| 国产亚洲精品综合一区在线观看| 亚洲精品在线美女| 久久这里只有精品中国| 好男人在线观看高清免费视频| 真实男女啪啪啪动态图| 亚洲成人久久性| 长腿黑丝高跟| 欧美色欧美亚洲另类二区| 国产伦一二天堂av在线观看| 在线观看一区二区三区| 日日干狠狠操夜夜爽| 久久草成人影院| 特级一级黄色大片| 国产精品久久久人人做人人爽| 综合色av麻豆| 亚洲第一欧美日韩一区二区三区| 久久久久久久久免费视频了| 国产熟女xx| 精品国产亚洲在线| 精品一区二区三区av网在线观看| 亚洲七黄色美女视频| 搡老妇女老女人老熟妇| 99久久无色码亚洲精品果冻| 亚洲国产色片| 欧美成人免费av一区二区三区| 真人一进一出gif抽搐免费| 在线看三级毛片| 亚洲国产中文字幕在线视频| 亚洲第一欧美日韩一区二区三区| 亚洲avbb在线观看| 国产精品 欧美亚洲| 首页视频小说图片口味搜索| 久久国产精品人妻蜜桃| 国产亚洲精品一区二区www| 精品一区二区三区视频在线观看免费| 久久婷婷人人爽人人干人人爱| 69av精品久久久久久| 伦理电影免费视频| 麻豆久久精品国产亚洲av| 亚洲中文av在线| 亚洲人成伊人成综合网2020| 久久国产精品影院| 日本精品一区二区三区蜜桃| 最新在线观看一区二区三区| 中文亚洲av片在线观看爽| 亚洲熟妇中文字幕五十中出| 亚洲人成伊人成综合网2020| 免费av不卡在线播放| 在线免费观看的www视频| 精品一区二区三区视频在线 | 韩国av一区二区三区四区| 午夜福利18| 成人午夜高清在线视频| 久久久久精品国产欧美久久久| 九九在线视频观看精品| 久久精品综合一区二区三区| 亚洲男人的天堂狠狠| 亚洲中文字幕日韩| 国产午夜精品久久久久久| 两性夫妻黄色片| 日韩有码中文字幕| 1024香蕉在线观看| 最新中文字幕久久久久 | 国产高清视频在线观看网站| a级毛片a级免费在线| 一个人看视频在线观看www免费 | 搞女人的毛片| 欧美三级亚洲精品| 成人三级黄色视频| www.999成人在线观看| 99国产综合亚洲精品| 国产精品美女特级片免费视频播放器 | 一级毛片高清免费大全| 亚洲国产精品合色在线| 嫩草影视91久久| 精品国产三级普通话版| 制服丝袜大香蕉在线| 免费观看人在逋| avwww免费| 亚洲精品粉嫩美女一区| 91久久精品国产一区二区成人 | 桃红色精品国产亚洲av| 成人鲁丝片一二三区免费| 国产伦在线观看视频一区| 亚洲自拍偷在线| 国产成人精品久久二区二区91| 亚洲第一欧美日韩一区二区三区| 91在线观看av| www.自偷自拍.com| 久久久久久国产a免费观看| 无人区码免费观看不卡| 人人妻人人看人人澡| 国产乱人视频| 亚洲av日韩精品久久久久久密| 国产亚洲av嫩草精品影院| 制服人妻中文乱码| 国产一级毛片七仙女欲春2| 日本精品一区二区三区蜜桃| 最近在线观看免费完整版| 亚洲人成网站高清观看| 亚洲一区二区三区色噜噜| 亚洲片人在线观看| 男女视频在线观看网站免费| 色尼玛亚洲综合影院| 黑人欧美特级aaaaaa片| 午夜两性在线视频| 久久精品夜夜夜夜夜久久蜜豆| 亚洲七黄色美女视频| 国产真实乱freesex| 国产aⅴ精品一区二区三区波| 好看av亚洲va欧美ⅴa在| 亚洲成人精品中文字幕电影| 在线观看美女被高潮喷水网站 | 亚洲国产色片| 日韩大尺度精品在线看网址| 12—13女人毛片做爰片一| 99精品欧美一区二区三区四区| 亚洲av第一区精品v没综合| 狂野欧美白嫩少妇大欣赏| 精品久久久久久久人妻蜜臀av| 全区人妻精品视频| 午夜影院日韩av| 综合色av麻豆| 桃红色精品国产亚洲av| 日韩欧美在线二视频| 好看av亚洲va欧美ⅴa在| 日本五十路高清| www.999成人在线观看| www.www免费av| 久久午夜综合久久蜜桃| 久久国产乱子伦精品免费另类| 欧美在线一区亚洲| 成人欧美大片| 久久久久久久久久黄片| 亚洲国产中文字幕在线视频| 国内少妇人妻偷人精品xxx网站 | 成年女人看的毛片在线观看| 免费人成视频x8x8入口观看| 日韩大尺度精品在线看网址| 成人av一区二区三区在线看| 观看美女的网站| 成年女人毛片免费观看观看9| 日日夜夜操网爽| 色在线成人网| 操出白浆在线播放| 手机成人av网站| 天堂√8在线中文| 成人三级做爰电影| av在线蜜桃| 欧美zozozo另类| 丁香六月欧美| 18禁国产床啪视频网站| 少妇的丰满在线观看| 夜夜看夜夜爽夜夜摸| 俄罗斯特黄特色一大片| 国产av麻豆久久久久久久| 国产精品av视频在线免费观看| 性色avwww在线观看| 亚洲国产欧美网| 国产一区在线观看成人免费| 亚洲成a人片在线一区二区| 亚洲电影在线观看av| 性色av乱码一区二区三区2| 99久久成人亚洲精品观看| 久久久久亚洲av毛片大全| 十八禁人妻一区二区| 色综合婷婷激情| 黄色女人牲交| 日本黄色片子视频| 午夜福利在线在线| 可以在线观看的亚洲视频| 亚洲欧美精品综合一区二区三区| 国产精品av久久久久免费| 久久久久精品国产欧美久久久| 精品国产超薄肉色丝袜足j| 在线观看美女被高潮喷水网站 | 日本一二三区视频观看| 黄色成人免费大全| 亚洲黑人精品在线| 久久精品国产综合久久久| 免费av毛片视频| 又爽又黄无遮挡网站| 12—13女人毛片做爰片一| 精品久久久久久久末码| 久久久久免费精品人妻一区二区| 老汉色av国产亚洲站长工具| av黄色大香蕉| 日本 av在线| aaaaa片日本免费| 国产极品精品免费视频能看的| 制服丝袜大香蕉在线| 一级毛片精品| 亚洲精品色激情综合| 成人三级做爰电影| 日本免费一区二区三区高清不卡| 一个人观看的视频www高清免费观看 | 十八禁网站免费在线| 人妻久久中文字幕网| 午夜福利在线在线| 美女 人体艺术 gogo| 国产一级毛片七仙女欲春2| 亚洲国产欧美一区二区综合| 久久久久免费精品人妻一区二区| 女生性感内裤真人,穿戴方法视频| 熟女少妇亚洲综合色aaa.| 久久久久久九九精品二区国产| 亚洲国产精品成人综合色| 搡老熟女国产l中国老女人| 19禁男女啪啪无遮挡网站| 久久久久性生活片| 欧美+亚洲+日韩+国产| 丝袜人妻中文字幕| 天天一区二区日本电影三级| 国产伦人伦偷精品视频| 精品国内亚洲2022精品成人| 亚洲,欧美精品.| 狂野欧美激情性xxxx| 国内精品美女久久久久久| 麻豆av在线久日| 亚洲自偷自拍图片 自拍| 日韩三级视频一区二区三区| 中文字幕久久专区| 亚洲18禁久久av| 亚洲精品乱码久久久v下载方式 | 国产成人精品久久二区二区91| 天堂av国产一区二区熟女人妻| 久久天堂一区二区三区四区| 精品日产1卡2卡| 午夜免费成人在线视频| 18禁黄网站禁片午夜丰满| 熟女电影av网| 在线观看午夜福利视频| 一级黄色大片毛片| 少妇裸体淫交视频免费看高清| 国产精品 欧美亚洲| 又粗又爽又猛毛片免费看| 久久九九热精品免费| 无限看片的www在线观看| 日韩欧美精品v在线| 午夜免费激情av| 不卡av一区二区三区| 无限看片的www在线观看| 岛国在线观看网站| 99热只有精品国产| 两人在一起打扑克的视频| 两性午夜刺激爽爽歪歪视频在线观看| 91av网一区二区| 国产亚洲精品久久久com| 成人一区二区视频在线观看| 伊人久久大香线蕉亚洲五| 亚洲精品一区av在线观看| 嫩草影院入口| 成人亚洲精品av一区二区| 国产野战对白在线观看| 中国美女看黄片| 麻豆一二三区av精品| 精品国产三级普通话版| 午夜激情福利司机影院| 国产亚洲av高清不卡| 亚洲国产欧美人成| 日韩欧美在线二视频| 嫩草影视91久久| 国产伦人伦偷精品视频| 午夜福利成人在线免费观看| 每晚都被弄得嗷嗷叫到高潮| 久久中文字幕一级| 日本在线视频免费播放| 男人舔奶头视频| 丝袜人妻中文字幕| 亚洲中文日韩欧美视频| 国产欧美日韩一区二区三| xxxwww97欧美| 午夜福利在线观看免费完整高清在 | 搡老岳熟女国产| av视频在线观看入口| 欧美一级a爱片免费观看看| 久久人人精品亚洲av| 日本一二三区视频观看| 亚洲专区字幕在线| 五月玫瑰六月丁香| 欧美极品一区二区三区四区| 不卡一级毛片| а√天堂www在线а√下载| 国产高清三级在线| 国产男靠女视频免费网站| 这个男人来自地球电影免费观看| 两个人看的免费小视频| 18禁观看日本| 天天躁日日操中文字幕| 美女被艹到高潮喷水动态| 男女下面进入的视频免费午夜| 男人和女人高潮做爰伦理| 国产伦一二天堂av在线观看| 日韩欧美三级三区| 亚洲第一欧美日韩一区二区三区| 成人特级黄色片久久久久久久| 90打野战视频偷拍视频| 蜜桃久久精品国产亚洲av| 国产男靠女视频免费网站| 日本撒尿小便嘘嘘汇集6| 神马国产精品三级电影在线观看| 精品无人区乱码1区二区| 1000部很黄的大片| 亚洲人成伊人成综合网2020| 亚洲欧美日韩无卡精品| 色av中文字幕| 91麻豆精品激情在线观看国产| x7x7x7水蜜桃| 亚洲av五月六月丁香网| av天堂中文字幕网| 国产精品日韩av在线免费观看| 欧美色欧美亚洲另类二区| 精品国产超薄肉色丝袜足j| 亚洲九九香蕉| 久久人妻av系列| 美女午夜性视频免费| 精品一区二区三区视频在线 | 91老司机精品| 男人舔奶头视频| 男插女下体视频免费在线播放| 精品国产乱子伦一区二区三区| 男女午夜视频在线观看| 最近最新免费中文字幕在线| 丁香六月欧美| 日韩欧美精品v在线| 别揉我奶头~嗯~啊~动态视频| 色综合婷婷激情| 1024香蕉在线观看| 操出白浆在线播放| 国产精品亚洲美女久久久| 中文字幕人妻丝袜一区二区| 国内精品久久久久久久电影| 亚洲国产高清在线一区二区三| 欧美一区二区国产精品久久精品| 99久国产av精品| 亚洲专区中文字幕在线| 天天躁日日操中文字幕| 欧美日韩国产亚洲二区| 亚洲国产欧洲综合997久久,| 法律面前人人平等表现在哪些方面| 不卡一级毛片| 国产高清视频在线播放一区| 18美女黄网站色大片免费观看| 午夜日韩欧美国产| 国产精品免费一区二区三区在线| 色综合欧美亚洲国产小说| 精品人妻1区二区| 亚洲色图 男人天堂 中文字幕| 在线观看66精品国产| 亚洲av日韩精品久久久久久密| 久99久视频精品免费| 国产午夜福利久久久久久| 亚洲av熟女| 国产成人精品久久二区二区免费| 久久精品夜夜夜夜夜久久蜜豆| 午夜精品久久久久久毛片777| 99国产精品一区二区三区| 久久天堂一区二区三区四区| 成人性生交大片免费视频hd| 中文亚洲av片在线观看爽| 淫妇啪啪啪对白视频| 日韩成人在线观看一区二区三区| 嫁个100分男人电影在线观看| 成人三级做爰电影| 精品久久久久久久久久免费视频| av天堂在线播放| 免费搜索国产男女视频| 成年版毛片免费区| 在线观看美女被高潮喷水网站 | 18禁观看日本| 午夜精品在线福利| 黑人巨大精品欧美一区二区mp4| 国产精品爽爽va在线观看网站| 亚洲国产中文字幕在线视频| 12—13女人毛片做爰片一| 麻豆av在线久日| 国产伦在线观看视频一区| 亚洲中文av在线| 亚洲avbb在线观看| 日韩精品青青久久久久久| 日本一二三区视频观看| 婷婷精品国产亚洲av| 亚洲av熟女| 搞女人的毛片| 最近视频中文字幕2019在线8| 两个人看的免费小视频| 国产一区二区在线观看日韩 | 精品国产乱码久久久久久男人| 怎么达到女性高潮| 欧美黄色片欧美黄色片| 校园春色视频在线观看| 日本与韩国留学比较| 国产单亲对白刺激| 黄片小视频在线播放| 国产 一区 欧美 日韩| 琪琪午夜伦伦电影理论片6080| 日韩欧美在线乱码| 欧美黑人巨大hd| 三级国产精品欧美在线观看 | 久久久久亚洲av毛片大全| 国内精品久久久久久久电影| 久久亚洲精品不卡| 首页视频小说图片口味搜索| av天堂在线播放| АⅤ资源中文在线天堂| 露出奶头的视频| 精品一区二区三区四区五区乱码| 日韩精品中文字幕看吧| 极品教师在线免费播放| 久久久精品大字幕| 亚洲,欧美精品.| 久久伊人香网站| 国产伦精品一区二区三区视频9 | 亚洲成人免费电影在线观看| 亚洲欧美日韩卡通动漫| 成人欧美大片| 久久香蕉国产精品| 国产av麻豆久久久久久久| 国产成人福利小说| 欧美一级a爱片免费观看看| 日本精品一区二区三区蜜桃| 日本一本二区三区精品| 亚洲精品色激情综合| 天堂网av新在线| 一二三四社区在线视频社区8| 成人18禁在线播放| 九九久久精品国产亚洲av麻豆 | 日韩精品中文字幕看吧| 岛国在线免费视频观看| 婷婷亚洲欧美| 亚洲五月天丁香| 国产亚洲av高清不卡| cao死你这个sao货| 国产免费男女视频| 亚洲七黄色美女视频| 免费在线观看视频国产中文字幕亚洲| 亚洲av电影在线进入| 国产成人精品久久二区二区免费| 久久久精品欧美日韩精品| 免费在线观看视频国产中文字幕亚洲| 免费看美女性在线毛片视频| 国产伦一二天堂av在线观看| а√天堂www在线а√下载| 手机成人av网站| 免费搜索国产男女视频| 在线观看66精品国产| 亚洲专区中文字幕在线| 舔av片在线| 在线观看舔阴道视频| 综合色av麻豆| 国产亚洲精品av在线| 精品国产乱码久久久久久男人| 男女那种视频在线观看| 18禁国产床啪视频网站| 国产av一区在线观看免费| 国产探花在线观看一区二区| 国产一区在线观看成人免费| 久久久久国产一级毛片高清牌| 黄色成人免费大全| 丝袜人妻中文字幕| 特级一级黄色大片| 一夜夜www| 99国产精品一区二区三区| 女生性感内裤真人,穿戴方法视频| 欧美中文综合在线视频| 久久精品人妻少妇| 老司机午夜十八禁免费视频| 成人高潮视频无遮挡免费网站| 中文字幕人妻丝袜一区二区| 亚洲精品久久国产高清桃花| 一区二区三区激情视频| 日韩欧美免费精品| 一个人看视频在线观看www免费 | 久久精品影院6| 欧美日韩中文字幕国产精品一区二区三区| 两人在一起打扑克的视频| 黑人巨大精品欧美一区二区mp4| 亚洲精品久久国产高清桃花| 人人妻,人人澡人人爽秒播| 欧美大码av| 久久精品国产亚洲av香蕉五月| 在线观看日韩欧美| 国产乱人伦免费视频| 99国产极品粉嫩在线观看| 黄色片一级片一级黄色片| 白带黄色成豆腐渣| 最近视频中文字幕2019在线8| 欧美3d第一页| 人妻夜夜爽99麻豆av| 国产真实乱freesex| 日韩免费av在线播放| 午夜精品在线福利| 亚洲中文字幕日韩| 成年女人毛片免费观看观看9| 亚洲国产高清在线一区二区三| 1024香蕉在线观看| av在线天堂中文字幕| 在线观看舔阴道视频| 99久久无色码亚洲精品果冻| 一个人看的www免费观看视频| 18禁美女被吸乳视频| 99热这里只有精品一区 | av福利片在线观看| 18禁裸乳无遮挡免费网站照片| 看黄色毛片网站| 久久久久久大精品| 午夜精品久久久久久毛片777| 国产精品野战在线观看| 欧美zozozo另类| 国内精品美女久久久久久| 两人在一起打扑克的视频| h日本视频在线播放| av女优亚洲男人天堂 | 99在线人妻在线中文字幕| 成人精品一区二区免费| or卡值多少钱| 色综合婷婷激情| 日韩av在线大香蕉| 一个人看的www免费观看视频| 久久香蕉国产精品| 一本综合久久免费| 18禁美女被吸乳视频| 老司机深夜福利视频在线观看| 偷拍熟女少妇极品色| 一个人看的www免费观看视频| 欧美zozozo另类| 一个人看视频在线观看www免费 | 波多野结衣高清作品| 国产激情偷乱视频一区二区| 亚洲国产精品成人综合色| 国产精品国产高清国产av| 国产伦人伦偷精品视频| 亚洲无线观看免费| 在线观看午夜福利视频| 婷婷精品国产亚洲av在线|