• <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第二区| 又黄又爽又免费观看的视频| 国产激情欧美一区二区| 精品日产1卡2卡| 国产精品秋霞免费鲁丝片| 一区在线观看完整版| 一边摸一边抽搐一进一小说| 日韩视频一区二区在线观看| 国产精品一区二区三区四区久久 | 国内精品久久久久精免费| 精品欧美一区二区三区在线| 日韩一卡2卡3卡4卡2021年| 免费看十八禁软件| 日本a在线网址| 久久香蕉国产精品| 99国产精品免费福利视频| АⅤ资源中文在线天堂| 久久中文看片网| 1024香蕉在线观看| 国产av一区二区精品久久| av在线播放免费不卡| 亚洲欧美日韩无卡精品| 91字幕亚洲| 12—13女人毛片做爰片一| 丝袜在线中文字幕| 女同久久另类99精品国产91| 久热爱精品视频在线9| 欧美日韩瑟瑟在线播放| 欧美乱码精品一区二区三区| 一级毛片高清免费大全| 午夜免费鲁丝| 两个人看的免费小视频| 久久影院123| 精品国产乱子伦一区二区三区| 国产人伦9x9x在线观看| 久久久久久久久免费视频了| 久久国产精品人妻蜜桃| 中文字幕av电影在线播放| 亚洲在线自拍视频| 国产精品一区二区在线不卡| www.熟女人妻精品国产| 国产亚洲欧美98| 无人区码免费观看不卡| 国产成人av教育| 国产极品粉嫩免费观看在线| 午夜免费激情av| 久久精品成人免费网站| 国产精品 国内视频| 亚洲国产欧美一区二区综合| 99国产精品99久久久久| 91大片在线观看| 午夜久久久久精精品| 国产成人啪精品午夜网站| 精品久久久久久成人av| 熟女少妇亚洲综合色aaa.| 国产精品久久电影中文字幕| 婷婷丁香在线五月| 中文字幕色久视频| 精品欧美一区二区三区在线| 国产成人欧美在线观看| 国产精品亚洲av一区麻豆| 国产亚洲欧美在线一区二区| 天堂动漫精品| 美女高潮喷水抽搐中文字幕| 亚洲男人天堂网一区| 亚洲av成人av| 亚洲精品粉嫩美女一区| 亚洲精品美女久久av网站| 中文字幕人妻丝袜一区二区| 国产成年人精品一区二区| 午夜福利,免费看| 91麻豆精品激情在线观看国产| 午夜福利一区二区在线看| 亚洲激情在线av| 亚洲av美国av| 国产精品一区二区三区四区久久 | 一级黄色大片毛片| 老熟妇仑乱视频hdxx| 97碰自拍视频| 亚洲国产精品sss在线观看| 日韩精品免费视频一区二区三区| 久久狼人影院| 精品人妻在线不人妻| 少妇被粗大的猛进出69影院| 亚洲成人国产一区在线观看| 国产成年人精品一区二区| 一二三四社区在线视频社区8| 亚洲精品美女久久久久99蜜臀| 黄色视频不卡| 亚洲无线在线观看| 制服丝袜大香蕉在线| 久久久久久久久久久久大奶| 日韩高清综合在线| 在线观看免费日韩欧美大片| 欧美日韩福利视频一区二区| 亚洲人成电影免费在线| 亚洲成a人片在线一区二区| 精品日产1卡2卡| 天堂影院成人在线观看| 国产在线精品亚洲第一网站| 日韩大尺度精品在线看网址 | 国产精品久久久av美女十八| 婷婷六月久久综合丁香| 免费在线观看完整版高清| 国产欧美日韩精品亚洲av| 日韩av在线大香蕉| 国产成人精品在线电影| 日本免费a在线| 色播亚洲综合网| 十八禁人妻一区二区| 亚洲色图 男人天堂 中文字幕| 三级毛片av免费| 最新在线观看一区二区三区| 精品欧美一区二区三区在线| 日本在线视频免费播放| 啦啦啦观看免费观看视频高清 | 久久精品国产亚洲av香蕉五月| 国产私拍福利视频在线观看| 天堂动漫精品| 嫩草影院精品99| 亚洲av熟女| 国产精华一区二区三区| av福利片在线| 精品国产一区二区久久| 久久精品国产99精品国产亚洲性色 | 国产成人精品在线电影| 性色av乱码一区二区三区2| 夜夜爽天天搞| 香蕉国产在线看| 国产成人一区二区三区免费视频网站| 女人精品久久久久毛片| 日韩欧美国产在线观看| 国产精品野战在线观看| 欧美乱色亚洲激情| 黄色成人免费大全| 久久婷婷成人综合色麻豆| 99国产极品粉嫩在线观看| 亚洲第一av免费看| 中文字幕人妻丝袜一区二区| av在线天堂中文字幕| 久久精品91无色码中文字幕| 99香蕉大伊视频| 国产成人啪精品午夜网站| 国产精品野战在线观看| 757午夜福利合集在线观看| 精品少妇一区二区三区视频日本电影| 精品电影一区二区在线| 国产精品久久久av美女十八| 亚洲激情在线av| 淫秽高清视频在线观看| 日日干狠狠操夜夜爽| 久久性视频一级片| 三级毛片av免费| 久久久久精品国产欧美久久久| 久久 成人 亚洲| 久久久久国产精品人妻aⅴ院| 日本黄色视频三级网站网址| 亚洲 国产 在线| 啦啦啦免费观看视频1| 宅男免费午夜| 国产精品av久久久久免费| 满18在线观看网站| av超薄肉色丝袜交足视频| 十八禁人妻一区二区| 国产区一区二久久| 伦理电影免费视频| 男男h啪啪无遮挡| 九色亚洲精品在线播放| 亚洲成人久久性| 日韩成人在线观看一区二区三区| av视频免费观看在线观看| 99在线人妻在线中文字幕| 久久中文字幕一级| 又大又爽又粗| 国产成人一区二区三区免费视频网站| 国产精品久久电影中文字幕| 99国产精品99久久久久| 曰老女人黄片| 日本免费一区二区三区高清不卡 | 桃色一区二区三区在线观看| 天堂影院成人在线观看| 韩国av一区二区三区四区| 老熟妇乱子伦视频在线观看| 亚洲专区中文字幕在线| 日韩欧美一区二区三区在线观看| 91成年电影在线观看| 一个人观看的视频www高清免费观看 | 少妇 在线观看| 男人舔女人的私密视频| 在线观看免费视频网站a站| АⅤ资源中文在线天堂| 大型av网站在线播放| 免费在线观看黄色视频的| www.精华液| 成人永久免费在线观看视频| 一区二区三区激情视频| 久久久国产欧美日韩av| 一二三四社区在线视频社区8| 色精品久久人妻99蜜桃| 精品一区二区三区四区五区乱码| 日本在线视频免费播放| 精品国产国语对白av| 一边摸一边做爽爽视频免费| 亚洲欧美日韩无卡精品| 91成年电影在线观看| 大香蕉久久成人网| 精品久久蜜臀av无| 亚洲精品av麻豆狂野| 中文亚洲av片在线观看爽| 亚洲第一av免费看| 黄频高清免费视频| 亚洲人成电影免费在线| 级片在线观看| 女性被躁到高潮视频| 亚洲少妇的诱惑av| 国产精品日韩av在线免费观看 | 欧美乱色亚洲激情| 亚洲国产欧美日韩在线播放| 在线观看一区二区三区| 19禁男女啪啪无遮挡网站| 午夜福利在线观看吧| 自拍欧美九色日韩亚洲蝌蚪91| 色综合婷婷激情| 动漫黄色视频在线观看| 精品国产美女av久久久久小说| 久久久久久大精品| 一区二区三区精品91| 多毛熟女@视频| 一区二区三区国产精品乱码| 精品国产超薄肉色丝袜足j| 啦啦啦免费观看视频1| 俄罗斯特黄特色一大片| 欧美色视频一区免费| 欧洲精品卡2卡3卡4卡5卡区| 91九色精品人成在线观看| 动漫黄色视频在线观看| 久久精品国产亚洲av香蕉五月| 亚洲性夜色夜夜综合| 在线观看免费午夜福利视频| 日本在线视频免费播放| 成人永久免费在线观看视频| 999久久久国产精品视频| 美女高潮喷水抽搐中文字幕| 亚洲国产欧美一区二区综合| 久久草成人影院| 国产1区2区3区精品| 在线观看www视频免费| 中文字幕人成人乱码亚洲影| 精品国产乱码久久久久久男人| 69精品国产乱码久久久| 国产aⅴ精品一区二区三区波| 中文字幕精品免费在线观看视频| 国产欧美日韩一区二区三区在线| 色综合亚洲欧美另类图片| 搡老妇女老女人老熟妇| 在线观看日韩欧美| 国产高清videossex| 欧美日韩精品网址| 国产野战对白在线观看| 亚洲性夜色夜夜综合| 给我免费播放毛片高清在线观看| 久久久久精品国产欧美久久久| 女性被躁到高潮视频| 搡老岳熟女国产| videosex国产| 国产欧美日韩一区二区精品| 91大片在线观看| 亚洲精品一卡2卡三卡4卡5卡| 村上凉子中文字幕在线| 久久久水蜜桃国产精品网| 嫁个100分男人电影在线观看| 成人18禁在线播放| 国产成人免费无遮挡视频| 麻豆久久精品国产亚洲av| 国产不卡一卡二| 黄色视频,在线免费观看| 日韩有码中文字幕| 高清在线国产一区| 亚洲午夜精品一区,二区,三区| 男女午夜视频在线观看| 女性生殖器流出的白浆| 1024视频免费在线观看| 成人亚洲精品av一区二区| 在线十欧美十亚洲十日本专区| 欧美日本亚洲视频在线播放| 琪琪午夜伦伦电影理论片6080| 亚洲精品久久国产高清桃花| 国产单亲对白刺激| 国产免费男女视频| 最好的美女福利视频网| 亚洲av日韩精品久久久久久密| 欧美日韩福利视频一区二区| 黄片播放在线免费| 999精品在线视频| 久久中文看片网| 一边摸一边抽搐一进一小说| 国内久久婷婷六月综合欲色啪| 欧美激情久久久久久爽电影 | 真人一进一出gif抽搐免费| 九色国产91popny在线| 亚洲欧美日韩另类电影网站| 一进一出好大好爽视频| 亚洲精品久久国产高清桃花| 久久久久精品国产欧美久久久| 一边摸一边抽搐一进一出视频| 天天躁狠狠躁夜夜躁狠狠躁| 亚洲一码二码三码区别大吗| 久久精品国产99精品国产亚洲性色 | 一级毛片精品| 男女午夜视频在线观看| 满18在线观看网站| 欧美另类亚洲清纯唯美| 国产亚洲精品久久久久久毛片| 黄色女人牲交| 亚洲精品久久国产高清桃花| 午夜久久久久精精品| 亚洲 欧美 日韩 在线 免费| 久久亚洲精品不卡| 给我免费播放毛片高清在线观看| 50天的宝宝边吃奶边哭怎么回事| 亚洲精品在线美女| 亚洲成av片中文字幕在线观看| 香蕉国产在线看| 国内精品久久久久精免费| 亚洲中文字幕日韩| 99在线人妻在线中文字幕| x7x7x7水蜜桃| 午夜激情av网站| 老汉色av国产亚洲站长工具| 国内精品久久久久久久电影| 国产精品亚洲美女久久久| 美女午夜性视频免费| 免费搜索国产男女视频| 国产日韩一区二区三区精品不卡| 妹子高潮喷水视频| 亚洲男人天堂网一区| 精品欧美国产一区二区三| 真人做人爱边吃奶动态| 午夜影院日韩av| 男女下面插进去视频免费观看| 国产成人啪精品午夜网站| 久久国产精品影院| 在线免费观看的www视频| tocl精华| 美国免费a级毛片| 黄色女人牲交| 在线观看日韩欧美| 咕卡用的链子| 国产又爽黄色视频| 午夜福利在线观看吧| 两个人看的免费小视频| 悠悠久久av| 91成年电影在线观看| 久久久国产成人精品二区| 老司机深夜福利视频在线观看| 国产主播在线观看一区二区| 18禁裸乳无遮挡免费网站照片 | 免费人成视频x8x8入口观看| 日本三级黄在线观看| 国产精品1区2区在线观看.| 久久中文看片网| 18禁黄网站禁片午夜丰满| 亚洲专区字幕在线| 在线观看免费视频网站a站| 日韩大码丰满熟妇| 亚洲色图综合在线观看| 亚洲av成人不卡在线观看播放网| а√天堂www在线а√下载| 法律面前人人平等表现在哪些方面| 亚洲欧美激情在线| 深夜精品福利| 18禁黄网站禁片午夜丰满| 亚洲国产欧美网| 亚洲自拍偷在线| 久久久久久久精品吃奶| 国产人伦9x9x在线观看| 亚洲一区二区三区色噜噜| 长腿黑丝高跟| 69精品国产乱码久久久| 亚洲中文日韩欧美视频| 色在线成人网| 国产成人精品无人区| 成人18禁在线播放| 国产亚洲精品久久久久久毛片| 成人免费观看视频高清| 欧美黄色淫秽网站| 一边摸一边抽搐一进一出视频| 美女高潮喷水抽搐中文字幕| 免费在线观看黄色视频的| 国产伦一二天堂av在线观看| 亚洲片人在线观看| 成人手机av| 一级毛片高清免费大全| 搡老岳熟女国产| 久久久久九九精品影院| 日日摸夜夜添夜夜添小说| 村上凉子中文字幕在线| 日本 欧美在线| 国产精品一区二区免费欧美| 91字幕亚洲| 性色av乱码一区二区三区2| 午夜久久久在线观看| 99久久久亚洲精品蜜臀av| 黄频高清免费视频| 国产人伦9x9x在线观看| 一进一出抽搐gif免费好疼| 俄罗斯特黄特色一大片| 久久伊人香网站| 一区二区三区高清视频在线| 51午夜福利影视在线观看| 高清在线国产一区| 伊人久久大香线蕉亚洲五| 男女下面插进去视频免费观看| 黄色a级毛片大全视频| 亚洲精品国产一区二区精华液| 国产一区二区三区视频了| 亚洲狠狠婷婷综合久久图片| 午夜福利成人在线免费观看| 国产亚洲精品久久久久久毛片| 91成年电影在线观看| 亚洲第一电影网av| 日本免费一区二区三区高清不卡 | 成人精品一区二区免费| 国产又色又爽无遮挡免费看| 嫁个100分男人电影在线观看| 91精品国产国语对白视频| 久久国产精品男人的天堂亚洲| 亚洲国产欧美一区二区综合| 99国产精品99久久久久| 久久亚洲精品不卡| 禁无遮挡网站| 欧美绝顶高潮抽搐喷水| 99热只有精品国产| av片东京热男人的天堂| 亚洲色图av天堂| 国产欧美日韩精品亚洲av| 精品欧美一区二区三区在线| 亚洲电影在线观看av| 日本免费a在线| 国内精品久久久久精免费| 国产一区在线观看成人免费| 国产一区在线观看成人免费| 免费看十八禁软件| 怎么达到女性高潮| 美女大奶头视频| 99riav亚洲国产免费| 看黄色毛片网站| 日本 av在线| 免费在线观看亚洲国产| 女同久久另类99精品国产91| 黄色视频,在线免费观看| 中文字幕最新亚洲高清| 99国产极品粉嫩在线观看| 91国产中文字幕| 搡老熟女国产l中国老女人| 精品乱码久久久久久99久播| 极品人妻少妇av视频| 亚洲人成网站在线播放欧美日韩| 国产伦人伦偷精品视频| 一区二区三区国产精品乱码| 黄色女人牲交| 成人国语在线视频| 国产高清视频在线播放一区| 电影成人av| 91在线观看av| 天堂√8在线中文| 免费久久久久久久精品成人欧美视频| 国产欧美日韩一区二区三区在线| 久久亚洲真实| 亚洲av片天天在线观看| 91麻豆精品激情在线观看国产| 午夜精品久久久久久毛片777| 久久国产精品影院| 免费高清在线观看日韩| 高潮久久久久久久久久久不卡| 男人舔女人下体高潮全视频| 国产极品粉嫩免费观看在线| 亚洲人成网站在线播放欧美日韩| 欧美日韩乱码在线| 成人欧美大片| 欧美激情 高清一区二区三区| 久久久国产成人精品二区| 成人特级黄色片久久久久久久| 久久人妻福利社区极品人妻图片| 50天的宝宝边吃奶边哭怎么回事| 一区二区日韩欧美中文字幕| 操美女的视频在线观看| 亚洲精品在线美女| 国产精品 国内视频| 老司机深夜福利视频在线观看| 美女高潮到喷水免费观看| 亚洲天堂国产精品一区在线| 国产精品久久久av美女十八| x7x7x7水蜜桃| 色哟哟哟哟哟哟| 老鸭窝网址在线观看| 亚洲最大成人中文| 91成年电影在线观看| 国产精品一区二区免费欧美| 精品一品国产午夜福利视频| 国产麻豆69| ponron亚洲| 大香蕉久久成人网| 成人国产综合亚洲| 久久人妻熟女aⅴ| 真人做人爱边吃奶动态| 人人妻人人爽人人添夜夜欢视频| 国产亚洲欧美精品永久| 国产亚洲精品第一综合不卡| 在线观看免费午夜福利视频| 麻豆成人av在线观看| 欧美在线一区亚洲| 亚洲一区二区三区不卡视频| 久久香蕉国产精品| 国产在线精品亚洲第一网站| 色综合亚洲欧美另类图片| www国产在线视频色| 国产精品久久久久久亚洲av鲁大| 国产av又大| 午夜两性在线视频| 日韩av在线大香蕉| 两性午夜刺激爽爽歪歪视频在线观看 | 最近最新免费中文字幕在线| 老司机深夜福利视频在线观看| 变态另类丝袜制服| 9色porny在线观看| 国产精品久久久人人做人人爽| 国产成人欧美在线观看| 乱人伦中国视频| 色精品久久人妻99蜜桃| 岛国在线观看网站| 亚洲欧美一区二区三区黑人| 国产亚洲精品一区二区www| 亚洲精品美女久久久久99蜜臀| av天堂久久9| 男女下面进入的视频免费午夜 | 女性生殖器流出的白浆| 久久人妻av系列| 两个人免费观看高清视频| 99久久综合精品五月天人人| 神马国产精品三级电影在线观看 | 女同久久另类99精品国产91| 最近最新中文字幕大全免费视频| 黑人操中国人逼视频| 国产成人精品在线电影| 亚洲中文日韩欧美视频| 久久香蕉国产精品| 日日爽夜夜爽网站| 天堂动漫精品| 久久人妻av系列| 亚洲第一青青草原| 在线播放国产精品三级| 巨乳人妻的诱惑在线观看| 欧美色欧美亚洲另类二区 | 亚洲黑人精品在线| 一边摸一边抽搐一进一出视频| 性色av乱码一区二区三区2| 91精品国产国语对白视频| 欧美激情极品国产一区二区三区| 亚洲精品国产一区二区精华液| 国产又爽黄色视频| 久久精品国产99精品国产亚洲性色 | 亚洲中文字幕一区二区三区有码在线看 | av网站免费在线观看视频| 欧美国产精品va在线观看不卡| 嫩草影院精品99| 99国产精品99久久久久| 国产亚洲欧美在线一区二区| 可以在线观看毛片的网站| 久久中文看片网| 伦理电影免费视频| 国产91精品成人一区二区三区| 日韩欧美国产在线观看| 91麻豆精品激情在线观看国产| 中文字幕最新亚洲高清| 在线永久观看黄色视频| 久久伊人香网站| 丰满的人妻完整版| 国产精品一区二区在线不卡| 国产一卡二卡三卡精品| 成人18禁高潮啪啪吃奶动态图| 1024香蕉在线观看| √禁漫天堂资源中文www| 日韩精品免费视频一区二区三区| 国产亚洲欧美精品永久| 男人舔女人下体高潮全视频| 亚洲国产中文字幕在线视频| 久久婷婷人人爽人人干人人爱 | 99久久久亚洲精品蜜臀av| 99香蕉大伊视频| 好男人电影高清在线观看| x7x7x7水蜜桃| 色在线成人网| 亚洲欧美精品综合久久99| 一区二区三区高清视频在线| 99精品久久久久人妻精品| 91字幕亚洲| 99久久精品国产亚洲精品| 91成人精品电影| 精品久久蜜臀av无| 两个人免费观看高清视频| 久久午夜综合久久蜜桃| АⅤ资源中文在线天堂| 免费一级毛片在线播放高清视频 | 亚洲 欧美一区二区三区| 99在线人妻在线中文字幕| 中文字幕久久专区| АⅤ资源中文在线天堂| 黄色毛片三级朝国网站| 嫩草影院精品99| 亚洲最大成人中文| 欧美成人一区二区免费高清观看 | 精品国产超薄肉色丝袜足j| 国产亚洲av嫩草精品影院| 午夜视频精品福利|