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

    彎曲網(wǎng)格上的間斷有限元湍流數(shù)值解法研究

    2014-04-30 07:24:44秦望龍呂宏強伍貽兆
    空氣動力學學報 2014年5期
    關(guān)鍵詞:物面粘性湍流

    秦望龍,呂宏強,伍貽兆

    (南京航空航天大學 航空宇航學院,江蘇 南京 210016)

    彎曲網(wǎng)格上的間斷有限元湍流數(shù)值解法研究

    秦望龍,呂宏強,伍貽兆

    (南京航空航天大學 航空宇航學院,江蘇 南京 210016)

    采用間斷有限元方法對雷諾平均Navier-Stokes(RANS)方程進行了數(shù)值求解,對Spalart-Allmaras單方程湍流模型進行了部分修正,使得求解器更加魯棒。構(gòu)造了分段高次多項式來逼近真實物面,同時物面附近采用多層彎曲網(wǎng)格來避免網(wǎng)格交叉,此外提出了一種快速計算積分點的曲面物面距的方法。采用混合網(wǎng)格對NACA0012翼型以及RAE翼型進行了數(shù)值計算,并與實驗數(shù)據(jù)以及前人數(shù)據(jù)進行了對比。計算結(jié)果表明,采用物面彎曲網(wǎng)格結(jié)合修正的湍流模型方法在相對稀疏的網(wǎng)格上就能得到比較好的數(shù)值解。

    間斷有限元;RANS;Spalart-Allmaras湍流模型;物面距;混合網(wǎng)格

    0 引 言

    近年來,隨著眾多學者的深入研究,間斷Galerkin有限元(Discontinous Galerkin,DG)方法求解雷諾平均Navier-Stokes(Reynolds-averaged Navier-Stokes,RANS)方程取得了迅速發(fā)展[1-11]。目前采用DG方法求解RANS方程的湍流模型主要有一方程Spalart-Allmaras模型[4-10]、兩方程Wilcoxk-ω模型[1-3]以及兩方程剪切力運輸(SST)模型[11]。

    由于RANS方程本身的非光順性質(zhì),間斷有限元方法對其求解比較困難。尤其當采用稀疏網(wǎng)格,單元內(nèi)會出現(xiàn)吉布斯(Gibbs)現(xiàn)象導致求解失敗。針對SA一方程模型的穩(wěn)定求解問題,一些學者提出了相應(yīng)的方法[4-7]。文中采用一種修正的S-A一方程模型較好的抑制了Gibbs現(xiàn)象,順利求解了RANS方程。

    間斷有限元方法對物面形狀的表達精度非常敏感,如果物面采用分段線性網(wǎng)格進行計算會導致求解不收斂或數(shù)值結(jié)果不光順[12-13]。為了解決以上問題,一般采用二次及以上多項式進行高階近似以逼近真實物面,文獻[14-18]給出了不同的高階近似方法。然而在求解粘性流動,尤其是雷諾數(shù)較大時,由于第一層網(wǎng)格厚度很小,物面高階近似會帶來網(wǎng)格的交叉纏繞,從而出現(xiàn)負體積。為了使求解順利進行,文中采用多層彎曲網(wǎng)格遞推的方法,以很小的計算代價解決了以上問題。

    Spalart-Allmaras單方程模型及SST兩方程模型的求解需要計算積分點的物面距。對于彎曲物面,積分點的物面距的求解比較困難,Liu等人采用有限元方法對Eikonal方程進行數(shù)值離散求解了內(nèi)流問題中的物面距[19]。然而,該方法在求解外流問題時會不穩(wěn)定,Schoenawa采用流線-擴散方法對Eikonal方程進行離散求解,并引入人工粘性方法使其穩(wěn)定魯棒[11]。為了不額外引入計算方程,本文采用一種簡單的數(shù)值方法對物面距進行計算。

    1 控制方程

    可壓縮的雷諾平均Navier-Stokes方程耦合修正的一方程Spalart-Allmaras湍流模型可以寫成如下的守恒形式:

    其中,Ω是有界域。守恒變量u、無粘數(shù)值通量Fc、粘性數(shù)值通量Fv可以寫成如下矢量形式:

    其中,ρ、p、e0、h分別是密度、壓強、單位總能和單位總焓,ui=(vx,vy)是笛卡爾坐標系下的速度分量。壓強p滿足理想氣體方程:

    γ是定壓比熱容Cp與定容比熱容Cv的比值,文中γ取1.4。σij是粘性應(yīng)力張量,對于牛頓流體,

    其中,μ是動力粘性系數(shù),可以通過Sutherland公式得出。μT是湍流粘性系數(shù):

    湍流生成項和消散項中的參數(shù)分別如下:

    2 數(shù)值方法

    2.1 數(shù)值離散

    將計算域表示成單元集合Ω=∪kΩk,文中Ωk是四邊形單元或三角形單元。定義Vh,p是單元h上采用多項式基函數(shù)表示的p階有限元空間,p≥0。單元變量uh可以在Vh,p下表示,即uh∈Vh,p。對方程(1)離散,有

    式(10)中,F(xiàn)(uh,?uh)是對流通量和粘性通量的和。邊界通量H=F(uh,?uh)·n的求解需要引入合適的數(shù)值通量函數(shù)。文中采用LLF數(shù)值通量求解邊界無粘數(shù)值通量[21]。對于粘性項,采用BR2格式(the second scheme of Bassi and Rebay)[1]進行數(shù)值離散。

    2.2 時間積分

    數(shù)值離散之后,方程最終可以寫成一個常微分方程系統(tǒng):

    式(11)中,M是質(zhì)量矩陣,R(Uh)是殘值,Uh是未知變量。文中采用隱式向后歐拉方法對式(11)進行數(shù)值求解。每個迭代步產(chǎn)生一個線性方程組:

    對于式(12)的求解文中采用ILU預(yù)處理的廣義最小余量方法[22]。

    3 網(wǎng)格彎曲處理

    高階間斷有限元方法對物面形狀的表達精度非常敏感,文中采用三次多項式對物面邊進行分段高階近似,從而逼近真實物面。求解高雷諾數(shù)粘性流動時,在物面進行高階近似的同時需要對相鄰外層網(wǎng)格進行高階近似才能避免網(wǎng)格交叉纏繞。

    區(qū)別于文獻[18]中的映射方法,文中采用參數(shù)方程實現(xiàn)物面高階近似和曲面網(wǎng)格陣推。如圖1,將曲線AB表示成參數(shù)方程形式:

    根據(jù)邊界點A、B的坐標及加權(quán)方法算出的法矢nA、nB,可以求出上述方程的八個參數(shù),從而實現(xiàn)物面彎曲。在此基礎(chǔ)求出AB線段1/3、2/3點a、b處的彎曲距離α、β。根據(jù)每一層的陣推距離K1α、K2β以及邊界點A1、B1,我們可以求出所有彎曲邊的高階參數(shù)方程形式。這里,K1、K2是調(diào)節(jié)因子,用于調(diào)節(jié)該面彎曲的程度,文中均取1.0。采用以上方法可以求出任意彎曲層數(shù)下彎曲邊的表達形式,進一步提高了網(wǎng)格彎曲的自動化,并且為計算曲面物面距中搜索點的法矢提供了保證。需要指出的是,在實現(xiàn)高階近似過程中,對于物面一些特殊點,如翼型后緣尖點,由于過該點的相鄰邊不是圓滑過渡,所以該點的法矢不能簡單采用加權(quán)平均的方法進行計算,否則逼近的外形會偏離真實外形。

    圖1 網(wǎng)格單元高階近似示意圖Fig.1 Illustration of high order approximation

    4 曲面物面距

    為了不額外引入方程而準確計算積分點的曲面物面距,我們對物面邊進行遍歷,求出積分點到物面的最小距離d。

    對于固定曲面物面,根據(jù)積分點的位置可以將其分為三類(圖2),然后根據(jù)式(14)分別求解。

    當積分點處于b情況時,對該曲面物面進行二分法遍歷,直到在曲面上找到點g滿足n(g)·l(bg)=cosα=0。n(g)表示曲面上g點的單位法矢,l(bg)表示b點到g點的單位向量,α是兩向量的夾角。為了減少遍歷次數(shù),文中α取不小于80°。采用以上方法對RAE翼型的物面距的計算結(jié)果如圖3,等勢線顯示范圍為[0.001,0.5]。該算例計算了2695個網(wǎng)格單元中心點到物面的距離,其中1320個點需要計算case b情形。由于對α的設(shè)置,該算例case b情形的遍歷次數(shù)基本在2次及以下。

    圖2 高斯積分點的分類Fig.2 Possible cases of Gaussian quadrature points

    圖3 RAE翼型的物面距分布Fig.3 Wall distance solution of RAE airfoil

    5 數(shù)值結(jié)果與分析

    采用上述高階間斷有限元方法結(jié)合物面彎曲混合網(wǎng)格單元對RANS方程進行數(shù)值求解。

    5.1 NACA0012翼型算例

    計算來流條件為Ma=0.3,α=0°,Re=1.85×106。高雷諾數(shù)條件決定了第一層網(wǎng)格的厚度很薄,文中第一層網(wǎng)格取弦長的2×10-5,對應(yīng)Y+=4。計算網(wǎng)格單元如圖4,其中非結(jié)構(gòu)網(wǎng)格單元由結(jié)構(gòu)網(wǎng)格單元剖分而來。由于物面采用大長寬比結(jié)構(gòu)網(wǎng)格單元,該算例計算網(wǎng)格總數(shù)控制在2592。

    圖4 NACA0012翼型計算網(wǎng)格Fig.4 Mesh around NACA0012 airfoil

    圖5給出了計算空間點的物面距的等勢線,顯示范圍為[0.001,0.5]。

    圖5 NACA0012翼型的物面距分布Fig.5 Wall distance solution of NACA0012 airfoil

    圖6給出了本算例的馬赫云圖及湍流粘性系數(shù)云圖,圖7給出了本算例四階精度的表面壓力系數(shù)Cp分布及摩擦阻力系數(shù)Cf分布曲線。可以看出,Cp計算結(jié)果與實驗結(jié)果吻合較好。由于Cf沒有實驗數(shù)據(jù),這里我們采用Fluent軟件在14000個單元上計算得到的Cf數(shù)據(jù)作為比較,可以看出曲線走勢大致吻合。

    圖6 NACA0012翼型馬赫云圖及湍流粘性系數(shù)云圖Fig.6 Mach contours and eddy viscosity contours of NACA0012 airfoil

    圖7 Cp和Cf分布曲線Fig.7 Cpand Cfdistribution of NACA0012 airfoil

    5.2 RAE翼型算例

    計算來流條件為Ma=0.4,α=2.79°,Re=6.5×106。采用混合網(wǎng)格單元,總數(shù)為2695(圖8)。計算精度從1階(p=0)到4階(p=3),為使p階初值盡可能的靠近最終收斂解,采用p-1階的收斂解作為其初值。圖9給出了本算例p=3時的馬赫云圖及湍流粘性系數(shù)云圖。圖10給出了本算例的表面壓力系數(shù)Cp分布及摩擦阻力系數(shù)Cf分布曲線,計算結(jié)果與參考文獻[8]的計算結(jié)果吻合較好。由于文中曲線數(shù)據(jù)均未采用任何光順方法,所以在單元邊界處會有局部不光滑現(xiàn)象。

    圖8 RAE翼型計算網(wǎng)格Fig.8 Mesh around RAE airfoil

    圖9 RAE翼型馬赫云圖及湍流粘性系數(shù)云圖Fig.9 Mach contours and eddy viscosity contours of RAE airfoil

    圖10 Cp和Cf分布曲線Fig.10 Cpand Cfdistribution of RAE airfoil

    6 結(jié) 論

    (1)采用高階間斷有限元方法對雷諾平均Navier-Stokes方程進行了數(shù)值離散和求解,湍流模型采用Spalart-Allmaras單方程湍流模型,并對其進行了部分修正,保證了數(shù)值求解的穩(wěn)定性和魯棒性。

    (2)采用混合網(wǎng)格單元,物面采用參數(shù)方程實現(xiàn)高階物面近似及曲面網(wǎng)格陣推,由于物面采用多層大長寬比彎曲結(jié)構(gòu)網(wǎng)格單元,外部流場計算使用非結(jié)構(gòu)網(wǎng)格單元,在較少的網(wǎng)格上就可以得到相當精度的數(shù)值解。

    (3)給出了一種計算曲面物面距的方法,以很少的計算時間得到了滿足計算需要的曲面物面距。

    (4)采用以上方法對NACA0012翼型和RAE翼型高雷諾數(shù)算例進行了數(shù)值計算,驗證了文中的物面距求解方法和湍流模型方法。

    [1]BASSI F,CRIVELLINI A,REBAY S,et al.Discontinuous Galerkin solution of the Reynolds-averaged Navier-Stokes andkωturbulence model equations[J].Computers&Fluids,2005,34(2):507-540.

    [2]BASSI F,CRIVELLINI A,GHIDONI A,et al.High-order discontinuous Galerkin discretization of transonic turbulent flows[C].In Proceedings of the 47th Aerospace Sciences Meeting,Oriando;AIAA-2009-180.

    [3]HARTMANN R,HELD J,LEICHT T.Adjoint-based error estimation and adaptive mesh refinement for the RANS andk-ω turbulence model equations[J].Journal of Computational Physics,2011,230(11):4268-4284.

    [4]NGUYEN N C,PERSSON P O,PERAIRE J.RANS solutions using high order discontinuous Galerkin methods[C].The 45th Aerospace Sciences Meeting,Reno,2007;AIAA-2007-914.

    [5]OLIVER T,DARMOFAL D.An unsteady adaptation algorithm for discontinuous Galerkin discretizations of the RANS equations[C].The 18th AIAA Computational Fluid Dynamics Conference,Miami,2007;AIAA-2007-3940.

    [6]LANDMANN B,KESSLER M,WAGNER S,et al.A parallel,high-order discontinuous Galerkin codes for laminar and turbulent flows[J].Computers&Fluids,2008,37(2):427-438.

    [7]MORO D,NGUYEN N C,PERAIRE J.Navier-Stokes solution using hybridizable discontinuous Galerkin methods[C].The 20th Computational Fluid Dynamics Conference,Honolulu,2011;AIAA-2011-3407.

    [8]BURGESS N K,MAVRIPLIS D J.Robust computation of turbulent flows using a discontinuous Galerkin method[C].The 50th Aerospace Sciences Meeting,Nashville,2012;AIAA-2012-0457.

    [9]BURGESS N K,MAVRIPLIS D J.High-Order Discontinuous Galerkin Methods for Turbulent High-lift Flows[C].The 7th International Conference on Computational Fluid Dynamics(ICCFD7),Big Island,2012,ICCFD7-4202.

    [10]WANG L,ANDERSON W K,ERWIN J T,et al.Solutions of high-order methods for three-dimensional compressible viscous flows[C].The 42nd AIAA Fluid Dynamics Conference and Exhibit,New Orleans,2012.

    [11]SCHOENAWA S,HARTMANN R.Discontinuous Galerkin discretization of the Reynolds-averaged Navier-Stokes equations with the shear-stress transport model[J].Journal of Computational Physics,2014.

    [12]LüH Q,WU Y Z,ZHOU C H,et al.High resolution of subsonic flows on coarse grids[J].Acta Aeronautica et Astronautica Sinca,2009,30(2):200-204(in Chinese).

    呂宏強,伍貽兆,周春華等.稀疏非結(jié)構(gòu)網(wǎng)格上的亞音速流高精度數(shù)值模擬[J].航空學報,2009,30(2):200-204.

    [13]LU H.High-order discontinuous Galerkin solution of low-Re viscous flows[J].Modern Physics Letters B,2009,23(03):309-312.

    [14]LüBON C,et al.High-order boundary discretization for discontinuous Galerkin codes[C].In Proceedings of the 24th Applied Aerodynamics Conference,F(xiàn)rancisco,2006;AIAA-2006-2822.

    [15]HARTMANN R,HELD J,LEICHT T,et al.Discontinuous Galerkin methods for computational aerodynamics-3D adaptive flow simulation with the DLR PADGE code[J].Aerospace Science and Technology,2010,14(7):512-519.

    [16]PERSSON P,PERAIRE J.Curved mesh generation and mesh refinement using Lagrangian solid mechanics[C].In Proceedings of the 47th Aerospace Sciences Meeting,Oriando,2009;AIAA-2009-949.

    [17]YU J,YAN C.Study on discontinuous Galerkin method for Navier-Stokes equations[J].Chinese Journal of Theoretical and Applied Mechanics,2010,42(5):962-969(in Chinese)

    于劍,閻超.Navier-Stokes方程間斷Galerkin有限元方法研究[J].力學學報,2010,42(5):962-969.

    [18]QIN W L,LüH Q,WU Y Z.High-order discontinuous Galerkin solution of N-S equations on hybrid mesh[J].Chinese Journal of Theoretical and Applied Mechani.,2013,45(6):987-991.(in Chinese)

    秦望龍,呂宏強,伍貽兆.基于混合網(wǎng)格的高階間斷有限元黏流數(shù)值解法[J].力學學報,2013,45(6):987-991.

    [19]LIU C B,NITHIARASU P,TUCKER P G.Wall distance calculation using the Eikonal/Hamilton-Jacobi equations on unstructured meshes:a finite element approach[J].Engineering Computations,2010,27(5):645-657.

    [20]SPALART P,ALLMARASS.A one-equation turbulence model for aerodynamic flows[C].In Proceedings of the 30th Aerospace Sciences Meeting and Exhibit,Reno,1992;AIAA-92-0439.

    [21]XIA Y D,WU Y Z,LüH Q,et al.Parallel computation of a high-order discontinuous galerkin method on unstructured grids[J].ACTA Aerodynamica Sinica,2011,29(5):537-541.(in Chinese)

    夏軼棟,伍貽兆,呂宏強,等.高階間斷有限元法的并行計算研究[J].空氣動力學學報,2011,29(5):537-541.

    [22]SAAD Y,SCHULTZ M H.GMRES:A generalized minimal residual algorithm for solving nonsymmetric linear systems[J].SIAM Journal on Scientific and Statistical Computing,1986,7(3):856-869.

    Discontinuous Galerkin solution of RAN Sequations on curved mesh

    QIN Wanglong,Lü Hongqiang,WU Yizhao
    (College of Aerospace Engineering,Nanjing University of Aeronautics and Astronautics,Nanjing 210016,China)

    Discontinuous Galerkin(DG)finite element method was adopted for the numerical approximation of the Reynolds-averaged Navier-Stokes(RANS)equations with the Spalart-Allmaras turbulence model.In order to make the solver robust,the original turbulence model equation were modified accordingly.Furthermore,high order approximation of the real solid boundary was used and several layers of curved meshes were constructed to avoid inconsistent mesh cross-overs.For the computation of the distance of each quadrature point to the nearest curved wall boundaries,a fast straightforward numerical method was proposed.The DG discretization of the RANS equations were demonstrated for turbulent flows past a NACA0012 airfoil and RAE airfoil based on hybrid mesh.Numerical results indicate that highly accurate solutions can be obtained with the modified turbulent equation on coarse curved hybrid mesh.

    Discontinuous Galerkin(DG);RANS;Spalart-Allmaras turbulent model;high order approximation;hybrid mesh

    V211.3

    Adoi:10.7638/kqdlxxb-2014.0058

    0258-1825(2014)05-0581-06

    2014-06-24;

    2014-08-10

    國家自然科學基金(11272152);航空科學基金(20101552018);江蘇高校優(yōu)勢學科建設(shè)工程資助項目

    秦望龍(1988-),男,江蘇南通人,博士研究生,主要研究方向:計算流體力學,高階間斷有限元.E-mail:qinwanglong@126.com

    呂宏強,博士,副教授.E-mail:hongqiang.lu@126.com

    秦望龍,呂宏強,伍貽兆.彎曲網(wǎng)格上的間斷有限元湍流數(shù)值解法研究[J].空氣動力學學報,2014,32(5):581-586.

    10.7638/kqdlxxb-2014.0058. QIN W L,LU H Q,WU Y Z.Discontinuous Galerkin solution of RAN Sequations on curved mesh[J].ACTA Aerodynamica Sinica,2014,32(5):581-586.

    猜你喜歡
    物面粘性湍流
    一類具有粘性項的擬線性拋物型方程組
    激波/湍流邊界層干擾壓力脈動特性數(shù)值研究1)
    力學學報(2021年7期)2021-11-09 06:26:08
    帶粘性的波動方程組解的逐點估計
    重氣瞬時泄漏擴散的湍流模型驗證
    粘性非等熵流體方程平衡解的穩(wěn)定性
    讓吸盤掛鉤更牢固
    新型單面陣自由曲面光學測量方法成像特性仿真
    家庭醫(yī)生增強基層首診粘性
    賽車空動優(yōu)化的秘密
    汽車之友(2014年21期)2014-11-03 17:33:45
    “青春期”湍流中的智慧引渡(三)
    九九久久精品国产亚洲av麻豆| 全区人妻精品视频| 亚洲久久久久久中文字幕| 国产黄色小视频在线观看| av天堂中文字幕网| 两个人看的免费小视频| 非洲黑人性xxxx精品又粗又长| 久久久久久大精品| 中文字幕人成人乱码亚洲影| 99国产精品一区二区三区| 免费在线观看日本一区| 狂野欧美白嫩少妇大欣赏| 成年女人永久免费观看视频| 一级a爱片免费观看的视频| 老熟妇乱子伦视频在线观看| 国产成人a区在线观看| 精品午夜福利视频在线观看一区| 别揉我奶头~嗯~啊~动态视频| а√天堂www在线а√下载| 美女高潮喷水抽搐中文字幕| 国产av在哪里看| 欧美最新免费一区二区三区 | 欧美黑人欧美精品刺激| 91久久精品电影网| 国产精品亚洲一级av第二区| 国产亚洲精品一区二区www| 亚洲第一电影网av| tocl精华| 精品久久久久久久末码| 亚洲性夜色夜夜综合| 在线免费观看的www视频| 99久久久亚洲精品蜜臀av| 国产成+人综合+亚洲专区| 亚洲不卡免费看| 国产精品久久久久久久久免 | 热99在线观看视频| 村上凉子中文字幕在线| 此物有八面人人有两片| 久久精品亚洲精品国产色婷小说| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 搡老岳熟女国产| 午夜免费成人在线视频| 国产高潮美女av| 大型黄色视频在线免费观看| 国产99白浆流出| 1000部很黄的大片| 国产高潮美女av| 一本综合久久免费| 看免费av毛片| 99热这里只有是精品50| 亚洲avbb在线观看| 亚洲国产中文字幕在线视频| 中出人妻视频一区二区| 桃色一区二区三区在线观看| 精品无人区乱码1区二区| 精品一区二区三区视频在线 | 两性午夜刺激爽爽歪歪视频在线观看| 亚洲国产日韩欧美精品在线观看 | 欧美午夜高清在线| 黄色女人牲交| 久久婷婷人人爽人人干人人爱| e午夜精品久久久久久久| 午夜精品在线福利| 成人无遮挡网站| 久久亚洲真实| 国产伦人伦偷精品视频| 久久久国产精品麻豆| 久久中文看片网| 久久99热这里只有精品18| 一区福利在线观看| 少妇丰满av| 欧美bdsm另类| 久久精品国产自在天天线| 欧美在线黄色| 亚洲片人在线观看| 亚洲片人在线观看| 国产精品久久久久久精品电影| 一级a爱片免费观看的视频| 国产真实伦视频高清在线观看 | 午夜激情欧美在线| 女警被强在线播放| 午夜激情福利司机影院| 久99久视频精品免费| 国产淫片久久久久久久久 | 国产伦精品一区二区三区视频9 | 亚洲七黄色美女视频| 午夜老司机福利剧场| 国产亚洲av嫩草精品影院| 高清毛片免费观看视频网站| 欧美性猛交╳xxx乱大交人| 久久久国产精品麻豆| 一个人免费在线观看的高清视频| 国产亚洲欧美在线一区二区| 国产一区二区激情短视频| 淫秽高清视频在线观看| 99久久精品热视频| 亚洲电影在线观看av| 亚洲av成人精品一区久久| 亚洲国产色片| 亚洲黑人精品在线| 美女高潮喷水抽搐中文字幕| 亚洲精品一区av在线观看| 老熟妇乱子伦视频在线观看| 亚洲第一电影网av| 国产精品嫩草影院av在线观看 | 每晚都被弄得嗷嗷叫到高潮| 亚洲天堂国产精品一区在线| 午夜老司机福利剧场| 欧美午夜高清在线| 亚洲中文字幕一区二区三区有码在线看| 精品一区二区三区av网在线观看| 禁无遮挡网站| 变态另类丝袜制服| 亚洲精品在线观看二区| 少妇高潮的动态图| 激情在线观看视频在线高清| 亚洲色图av天堂| 免费av不卡在线播放| 国产99白浆流出| 精品人妻偷拍中文字幕| 久久精品人妻少妇| 国产真实伦视频高清在线观看 | 国产精品美女特级片免费视频播放器| 99在线视频只有这里精品首页| 亚洲成人久久爱视频| 亚洲成人久久爱视频| 欧美3d第一页| 高清在线国产一区| 亚洲在线观看片| 久久精品国产综合久久久| 色综合婷婷激情| 欧美乱码精品一区二区三区| 免费看光身美女| 色老头精品视频在线观看| 欧美乱色亚洲激情| 在线十欧美十亚洲十日本专区| 亚洲欧美日韩东京热| 国产免费男女视频| 免费看光身美女| 99视频精品全部免费 在线| 一进一出抽搐gif免费好疼| 国产av在哪里看| 啦啦啦免费观看视频1| 国内揄拍国产精品人妻在线| 黄色女人牲交| 99精品欧美一区二区三区四区| 哪里可以看免费的av片| 亚洲国产欧美人成| 国产黄色小视频在线观看| 免费人成视频x8x8入口观看| 天天一区二区日本电影三级| 日韩精品中文字幕看吧| av福利片在线观看| 欧美三级亚洲精品| 悠悠久久av| 午夜免费成人在线视频| 亚洲精品日韩av片在线观看 | 90打野战视频偷拍视频| 日韩欧美免费精品| 欧美激情在线99| 十八禁人妻一区二区| 久99久视频精品免费| www国产在线视频色| 99久久成人亚洲精品观看| 啦啦啦观看免费观看视频高清| 午夜免费男女啪啪视频观看 | 久久精品国产综合久久久| 免费无遮挡裸体视频| 国产精品亚洲一级av第二区| 1000部很黄的大片| 久久人人精品亚洲av| 久久久久久久亚洲中文字幕 | 国产精品亚洲一级av第二区| 免费在线观看日本一区| 欧美中文日本在线观看视频| 欧美成人a在线观看| 欧美性猛交黑人性爽| 他把我摸到了高潮在线观看| 亚洲美女黄片视频| 欧美区成人在线视频| www日本在线高清视频| 亚洲天堂国产精品一区在线| 免费在线观看亚洲国产| 亚洲精品久久国产高清桃花| 偷拍熟女少妇极品色| 99久久精品热视频| 日本与韩国留学比较| 欧美日本亚洲视频在线播放| 国产伦精品一区二区三区视频9 | 母亲3免费完整高清在线观看| 亚洲欧美日韩高清专用| av片东京热男人的天堂| 高清在线国产一区| 好看av亚洲va欧美ⅴa在| h日本视频在线播放| 成人国产一区最新在线观看| 一a级毛片在线观看| 免费高清视频大片| 男女之事视频高清在线观看| 在线观看午夜福利视频| 亚洲欧美激情综合另类| 丰满人妻熟妇乱又伦精品不卡| 少妇的逼水好多| 精品福利观看| 久久精品人妻少妇| 九九久久精品国产亚洲av麻豆| 乱人视频在线观看| 日本 欧美在线| 国内精品久久久久久久电影| 国产aⅴ精品一区二区三区波| 18+在线观看网站| 国产亚洲精品久久久久久毛片| av中文乱码字幕在线| 岛国视频午夜一区免费看| 少妇裸体淫交视频免费看高清| 俺也久久电影网| 国产精品亚洲av一区麻豆| 真实男女啪啪啪动态图| 99国产精品一区二区三区| 一区二区三区高清视频在线| 最新中文字幕久久久久| 久久久久久国产a免费观看| 嫩草影院入口| 桃色一区二区三区在线观看| 一级毛片女人18水好多| 欧美日韩瑟瑟在线播放| 精品电影一区二区在线| 精品国产美女av久久久久小说| 欧美性感艳星| 午夜老司机福利剧场| 特大巨黑吊av在线直播| 欧美最黄视频在线播放免费| 亚洲精品一卡2卡三卡4卡5卡| 神马国产精品三级电影在线观看| 国产精品,欧美在线| 一夜夜www| 国产真实乱freesex| 级片在线观看| 欧美精品啪啪一区二区三区| 成人国产一区最新在线观看| 美女黄网站色视频| 51午夜福利影视在线观看| 桃红色精品国产亚洲av| 国产成人啪精品午夜网站| 国产精品亚洲美女久久久| 亚洲人成网站在线播放欧美日韩| 男人舔奶头视频| 五月玫瑰六月丁香| 午夜福利在线观看吧| 99久国产av精品| 免费观看精品视频网站| 老司机在亚洲福利影院| 99热这里只有精品一区| 国产高清三级在线| 99热这里只有精品一区| 成人精品一区二区免费| 99视频精品全部免费 在线| 亚洲片人在线观看| 黄色视频,在线免费观看| 最近在线观看免费完整版| svipshipincom国产片| 老司机深夜福利视频在线观看| 精品电影一区二区在线| 中文字幕人妻丝袜一区二区| 啪啪无遮挡十八禁网站| 性欧美人与动物交配| 国产伦精品一区二区三区四那| 午夜两性在线视频| 国产一区在线观看成人免费| 噜噜噜噜噜久久久久久91| 午夜福利在线观看吧| 首页视频小说图片口味搜索| 午夜日韩欧美国产| 色哟哟哟哟哟哟| 国产成年人精品一区二区| 69人妻影院| 又爽又黄无遮挡网站| 国产av麻豆久久久久久久| 亚洲精品在线观看二区| 国产一区二区在线av高清观看| 热99re8久久精品国产| 天堂网av新在线| 午夜两性在线视频| 国产高清激情床上av| 久久精品综合一区二区三区| 午夜免费激情av| 99久久成人亚洲精品观看| 99久久九九国产精品国产免费| 99国产极品粉嫩在线观看| 国产三级黄色录像| 热99re8久久精品国产| 啦啦啦观看免费观看视频高清| 身体一侧抽搐| 国产精品一及| 午夜免费观看网址| 少妇裸体淫交视频免费看高清| 无人区码免费观看不卡| 给我免费播放毛片高清在线观看| 国产成人aa在线观看| 国产 一区 欧美 日韩| 国产色爽女视频免费观看| 国产精品98久久久久久宅男小说| 国产一区在线观看成人免费| 一夜夜www| 欧美中文日本在线观看视频| 9191精品国产免费久久| av福利片在线观看| 午夜a级毛片| 两个人视频免费观看高清| 欧美一级毛片孕妇| 韩国av一区二区三区四区| 国产爱豆传媒在线观看| 亚洲内射少妇av| 在线观看日韩欧美| 一级毛片女人18水好多| 三级国产精品欧美在线观看| 日韩欧美三级三区| 夜夜爽天天搞| 国内久久婷婷六月综合欲色啪| 国产乱人伦免费视频| 国产高清激情床上av| 夜夜爽天天搞| 校园春色视频在线观看| 亚洲国产精品合色在线| 久久精品国产自在天天线| 女人十人毛片免费观看3o分钟| 午夜福利成人在线免费观看| 在线十欧美十亚洲十日本专区| 欧美在线黄色| 亚洲黑人精品在线| 法律面前人人平等表现在哪些方面| 老司机深夜福利视频在线观看| 日本与韩国留学比较| 国产精品女同一区二区软件 | 免费高清视频大片| 成人无遮挡网站| 久久伊人香网站| 在线观看免费视频日本深夜| 久久精品国产清高在天天线| 99riav亚洲国产免费| 中文在线观看免费www的网站| 国产老妇女一区| 久久中文看片网| 国产亚洲精品久久久com| 国产精品野战在线观看| 国产成人系列免费观看| 天天躁日日操中文字幕| 精品久久久久久久久久久久久| 在线天堂最新版资源| 波多野结衣高清作品| 久久午夜亚洲精品久久| 脱女人内裤的视频| 欧美+亚洲+日韩+国产| 中文字幕熟女人妻在线| 成年免费大片在线观看| 天天一区二区日本电影三级| 久久精品国产自在天天线| 国模一区二区三区四区视频| 中文在线观看免费www的网站| 国产精品亚洲一级av第二区| 免费av毛片视频| 日本一二三区视频观看| xxx96com| 一区二区三区高清视频在线| 久久精品综合一区二区三区| aaaaa片日本免费| 757午夜福利合集在线观看| 欧美bdsm另类| 久久香蕉国产精品| 亚洲人成电影免费在线| 日日摸夜夜添夜夜添小说| 午夜福利在线在线| 一个人免费在线观看电影| 女人被狂操c到高潮| 老熟妇乱子伦视频在线观看| www.色视频.com| 国产淫片久久久久久久久 | 亚洲精品影视一区二区三区av| svipshipincom国产片| 亚洲电影在线观看av| АⅤ资源中文在线天堂| 成人av在线播放网站| 欧美日本视频| 精华霜和精华液先用哪个| 一个人观看的视频www高清免费观看| 亚洲男人的天堂狠狠| 国产精品久久久久久久电影 | 国内久久婷婷六月综合欲色啪| 国产精品久久久久久久电影 | 国产色爽女视频免费观看| 一区二区三区高清视频在线| 中文字幕熟女人妻在线| 中出人妻视频一区二区| 国产精品99久久99久久久不卡| 欧美黑人欧美精品刺激| 99热6这里只有精品| 国产精品国产高清国产av| 色综合亚洲欧美另类图片| 一本一本综合久久| 又黄又粗又硬又大视频| 午夜a级毛片| 欧美成狂野欧美在线观看| 在线国产一区二区在线| 男女午夜视频在线观看| 国产精品久久久久久人妻精品电影| 国产精品 欧美亚洲| 1000部很黄的大片| 99国产综合亚洲精品| 色尼玛亚洲综合影院| 在线观看av片永久免费下载| 俄罗斯特黄特色一大片| av黄色大香蕉| ponron亚洲| 天天躁日日操中文字幕| 精品久久久久久久久久免费视频| 国产免费一级a男人的天堂| 国产一区二区三区在线臀色熟女| 亚洲av电影不卡..在线观看| 欧美三级亚洲精品| 俄罗斯特黄特色一大片| 国产在视频线在精品| 日本免费一区二区三区高清不卡| 性色av乱码一区二区三区2| 中文字幕久久专区| 草草在线视频免费看| 91麻豆av在线| 村上凉子中文字幕在线| 黄色片一级片一级黄色片| 欧美日韩瑟瑟在线播放| 黄片小视频在线播放| 99热6这里只有精品| 色老头精品视频在线观看| 波多野结衣高清作品| 成年免费大片在线观看| 一级黄色大片毛片| 久久久久久国产a免费观看| 少妇的丰满在线观看| 在线天堂最新版资源| 国产精品精品国产色婷婷| 精品99又大又爽又粗少妇毛片 | 在线观看午夜福利视频| 日日摸夜夜添夜夜添小说| 最近最新免费中文字幕在线| 久久精品国产99精品国产亚洲性色| 黄色成人免费大全| 淫秽高清视频在线观看| 亚洲avbb在线观看| 日韩欧美 国产精品| 丝袜美腿在线中文| 脱女人内裤的视频| 国产主播在线观看一区二区| 男人和女人高潮做爰伦理| 欧美+日韩+精品| 久久久久久久久久黄片| 亚洲av电影不卡..在线观看| 欧美+亚洲+日韩+国产| 好男人在线观看高清免费视频| 又紧又爽又黄一区二区| 久久久久久久精品吃奶| 看黄色毛片网站| 一区福利在线观看| 18禁黄网站禁片午夜丰满| 日韩国内少妇激情av| 久久精品影院6| 欧美另类亚洲清纯唯美| 一区二区三区高清视频在线| 欧美在线一区亚洲| 成人一区二区视频在线观看| 亚洲片人在线观看| 在线免费观看不下载黄p国产 | 午夜亚洲福利在线播放| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 蜜桃亚洲精品一区二区三区| 老汉色∧v一级毛片| 噜噜噜噜噜久久久久久91| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 桃色一区二区三区在线观看| 色哟哟哟哟哟哟| 亚洲五月天丁香| 9191精品国产免费久久| 久久久久久久久久黄片| 国产精品亚洲av一区麻豆| 1000部很黄的大片| 午夜福利在线观看吧| 成人永久免费在线观看视频| 久久中文看片网| 欧美日韩亚洲国产一区二区在线观看| 国产爱豆传媒在线观看| 欧美xxxx黑人xx丫x性爽| 国产午夜精品论理片| 天堂动漫精品| 国产精品久久电影中文字幕| 神马国产精品三级电影在线观看| 日本一二三区视频观看| 一个人看视频在线观看www免费 | 欧美丝袜亚洲另类 | 中国美女看黄片| 91九色精品人成在线观看| 波多野结衣高清作品| 精品久久久久久成人av| av天堂在线播放| 美女cb高潮喷水在线观看| 日韩精品中文字幕看吧| 露出奶头的视频| 欧美大码av| 色噜噜av男人的天堂激情| 麻豆国产av国片精品| 欧美日韩综合久久久久久 | 三级毛片av免费| av天堂中文字幕网| 亚洲欧美一区二区三区黑人| 免费一级毛片在线播放高清视频| 久久久久久久久大av| 国产精品久久久久久久久免 | 国产精品一区二区免费欧美| 亚洲精品成人久久久久久| 成人性生交大片免费视频hd| 久久国产精品影院| 最近最新中文字幕大全免费视频| 国产精品野战在线观看| 中文字幕人妻丝袜一区二区| 色尼玛亚洲综合影院| 3wmmmm亚洲av在线观看| 非洲黑人性xxxx精品又粗又长| 色哟哟哟哟哟哟| 久久久久九九精品影院| 免费看光身美女| 黄色视频,在线免费观看| 免费在线观看日本一区| 久久中文看片网| 国内久久婷婷六月综合欲色啪| 亚洲成人精品中文字幕电影| 国产成人av教育| 人妻夜夜爽99麻豆av| 黄片大片在线免费观看| 网址你懂的国产日韩在线| 欧美日韩乱码在线| 韩国av一区二区三区四区| 床上黄色一级片| 国产在线精品亚洲第一网站| 国产高清激情床上av| 亚洲男人的天堂狠狠| 亚洲成人精品中文字幕电影| 亚洲欧美激情综合另类| 午夜福利高清视频| 给我免费播放毛片高清在线观看| 变态另类成人亚洲欧美熟女| 日韩av在线大香蕉| 激情在线观看视频在线高清| 欧美激情在线99| 又紧又爽又黄一区二区| 午夜福利视频1000在线观看| 动漫黄色视频在线观看| 欧美xxxx黑人xx丫x性爽| 国产精品99久久99久久久不卡| 国产伦人伦偷精品视频| 亚洲中文日韩欧美视频| 久久久久久国产a免费观看| 脱女人内裤的视频| 亚洲精品色激情综合| 亚洲天堂国产精品一区在线| 男人舔女人下体高潮全视频| 亚洲第一欧美日韩一区二区三区| 最好的美女福利视频网| aaaaa片日本免费| 国产精品影院久久| bbb黄色大片| 久久久久久人人人人人| 99在线视频只有这里精品首页| 午夜两性在线视频| 99riav亚洲国产免费| 国产老妇女一区| 最后的刺客免费高清国语| av福利片在线观看| 天天躁日日操中文字幕| 亚洲精品一卡2卡三卡4卡5卡| 中亚洲国语对白在线视频| 国产乱人伦免费视频| 久久6这里有精品| 又黄又爽又免费观看的视频| 欧美激情在线99| 激情在线观看视频在线高清| 欧美在线一区亚洲| 成年版毛片免费区| 亚洲精品影视一区二区三区av| 国产精品99久久99久久久不卡| 亚洲欧美日韩卡通动漫| 亚洲国产精品合色在线| 亚洲国产精品成人综合色| 亚洲最大成人中文| 狠狠狠狠99中文字幕| 特级一级黄色大片| 国产精品女同一区二区软件 | 少妇的逼好多水| 国产成人av激情在线播放| 亚洲熟妇中文字幕五十中出| 精品国内亚洲2022精品成人| 国产一区二区在线av高清观看| 日日摸夜夜添夜夜添小说| 丰满乱子伦码专区| 9191精品国产免费久久| 99国产综合亚洲精品| 国产免费男女视频| 日日干狠狠操夜夜爽| 色老头精品视频在线观看| 夜夜爽天天搞| 欧美日韩综合久久久久久 | 成人亚洲精品av一区二区| 五月伊人婷婷丁香| 老司机福利观看| 国产主播在线观看一区二区| 国产av一区在线观看免费| 男女下面进入的视频免费午夜| 欧美国产日韩亚洲一区| 久久久久久久久久黄片|