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

    基于改進LRNk-ε模型的汽車氣動特性研究

    2015-12-29 02:59:56谷正氣,陳陣,黃泰明
    中國機械工程 2015年18期
    關鍵詞:風洞試驗

    基于改進LRNk-ε模型的汽車氣動特性研究

    谷正氣1,2陳陣1黃泰明1豐成杰1宗軼琦1

    1.湖南大學汽車車身先進設計制造國家重點實驗室,長沙,4100822.湖南工業(yè)大學,株洲,412007

    摘要:針對由Jones和Launder提出的LRN(低雷諾數(shù))k-ε湍流模型對轉捩預測不足等問題,引入全應力限制(TSL)方法及流線曲率因子以改進湍動能耗散率及湍流黏性系數(shù),最終獲得了一種改進的低雷諾數(shù)湍流模型,并將其應用于汽車外部流場計算仿真中,利用ISIS-CFD求解器實現(xiàn)計算。將計算結果與其他常用于汽車外部流場仿真的湍流模型以及風洞試驗得到的結果進行了對比分析,結果表明,在同等計算條件下改進LRN k-ε湍流模型能更加準確且高效地模擬車身表面的氣流分離以及尾部流場。

    關鍵詞:湍流模型;轉捩;汽車外流場;風洞試驗;Ahmed

    中圖分類號:U461.1

    收稿日期:2014-11-26

    基金項目:國家自然科學基金資助項目(50975083);交通運輸部新世紀十百千人才培養(yǎng)資助項目(20120222);“中國高水平汽車自主創(chuàng)新能力建設”資助項目;湖南大學汽車車身先進設計與制造國家重點實驗室自主課題資助項目(734215002);財政部創(chuàng)新團隊資助項目(0420036017)

    作者簡介:谷正氣,男,1963年生。湖南大學機械與運載工程學院教授、博士研究生導師,湖南工業(yè)大學教授。主要研究方向為汽車空氣動力學。獲部級科技一等獎2項,發(fā)表專著7部,論文150余篇。陳陣,男,1988年生。湖南大學機械與運載工程學院博士研究生。黃泰明,男,1982年生。湖南大學機械與運載工程學院博士研究生。豐成杰,男,1991年生。湖南大學機械與運載工程學院碩士研究生。宗軼琦,男,1987年生。湖南大學機械與運載工程學院博士研究生。

    Research on Aerodynamic Characteristics of Vehicle with an

    Improved LRNk-εTurbulence Model

    Gu Zhengqi1,2Chen Zhen1Huang Taiming1Feng Chengjie1Zong Yiqi1

    1.State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body,

    Hunan University,Changsha,410082

    2.Hunan University of Technology,Zhuzhou,Hunan,412007

    Abstract:Since the LRN(low Reynolds number) k-ε turbulence model presented by Jones and Launder failed to predict the laminar-turbulent transition,the streamline curvature factor and total stress limitation(TSL) were incorporated to improve the turbulent dissipation rate and turbulence viscosity coefficient.Thus, a new modified LRN k-ε turbulence model was obtained for the prediction of the aerodynamic characteristics of vehicles.Computations were carried out using the ISIS-CFD flow solver. The comparisons among the results obtained by the modified LRN k-ε model with ones by other commonly used turbulence models and wind tunnel experiments show that the new model is more accurate and efficient to simulate the flow separation around the body surface and the wake vortexes.

    Key words:turbulence model;transition;vehicle outflow field;windtunnel test;Ahmed

    0引言

    汽車的氣動特性直接影響其操縱穩(wěn)定性、動力性和燃油經(jīng)濟性[1-2]。作為研究汽車氣動特性的主要方法之一,數(shù)值模擬憑借其便捷性在汽車設計中得到了廣泛的應用。目前,對汽車外部流場的數(shù)值模擬大多采用RANS(時均)模式,其中k-ε湍流模型(k為湍動能,ε為湍動能耗散率)由于占用計算資源少、數(shù)值收斂快,并且能有效捕捉汽車復雜表面流場的物理特性,被廣泛應用于汽車外部流場的數(shù)值模擬[3-6]。k-ε湍流模型最顯著的兩個缺陷是:①即使引入壁面函數(shù),也無法彌補模型對邊界層中未充分發(fā)展湍流模擬的不足[7-8];②無法對轉捩現(xiàn)象進行準確模擬。針對無法準確模擬邊界層流動的問題,Jones等[9]將衰減函數(shù)和經(jīng)驗系數(shù)分別引入湍流黏度系數(shù)和生成項,提出了一種低雷諾數(shù)湍流模型;Abe等[10]通過引入Kolmogorov速度尺度消除了近壁面處低雷諾數(shù)影響;Chang等[11]提出了一種針對突擴流中奇點預測的低雷諾數(shù)湍流模型。低雷諾數(shù)湍流模型可避免使用經(jīng)驗式預測邊界層流動,且模型有一定的預測轉捩的能力。考慮到低雷諾數(shù)湍流模型并未深入剖析轉捩過程的物理機理,一般認為這種預測能力僅僅是一個巧合[12]。

    為了提高模型對轉捩的預測能力,Menter等[13]通過分析平板繞流等典型流動的轉捩機理,利用經(jīng)驗因子的輸運模型,提出了γ-Reθ轉捩模型。1997年,通過觀察到轉捩區(qū)域的漣漪,Mayle等[14]第一次建立了基于層流湍動能的轉捩模型。Cutrone等[15]對比了6種轉捩模型,結果表明層流湍動能轉捩模型最有效。Walters等[16-18]成功建立了“splat mechanism”模型,認為轉捩過程是壓力應變由壓制到激發(fā)的過程。文獻[19]將壓力應變效應引入Walters等的模型,得到了較好的轉捩模擬效果。文獻[20]通過對全應力作出限制,提高了低雷諾數(shù)k-ω模型對轉捩的模擬能力。

    對于汽車這種周邊充斥著極其復雜的高應變率及大曲率氣流的模型,要準確預測轉捩現(xiàn)象,需要進行更為行之有效的改進。鑒于文獻[20]提出的全應力限制(TSL)方法能使低雷諾數(shù)k-ω模型獲得良好的模擬轉捩的能力,本文將此方法應用于Jones等[9]的低雷諾數(shù)k-ε模型,同時,針對此方法在近壁面分離區(qū)對轉捩預測不足的缺陷,引入流線曲率因子加以改進,消除對湍流黏度的過高估計,使模型能在汽車流場計算中更好地捕捉到轉捩現(xiàn)象,進而更好地模擬汽車氣動特性。

    1改進低雷諾數(shù)湍流模型

    由Jones等[9]提出的基于時均納維-斯托克斯方程的低雷諾數(shù)湍流模型為

    ?(ρuik)/?xi=?[(μ+μt/σk)?k/?xi]/?xi+

    ?(ρuiε)/?xi=?[(μ+μt/σε)?ε/?xi]/?xi+cε1f1Pkε/k-

    μt=ρCμfμk2/ε

    其中

    fμ=exp[-2.5/(1+ReT/50)]

    ReT=ρk2/(με)

    式中,ρ為流體密度;xi、xj為坐標軸分量;ui、uj為速度分量;μ為分子黏性系數(shù);μt為湍流黏性系數(shù);Pk為湍動能產(chǎn)生項;Sij為時均應變率;常數(shù)η0=4.377,β=0.012,Cμ=0.09,f1=1,cε2=1.92。

    1.1基于TSL方法的改進

    由實驗結果及直接數(shù)值計算(DNS)結果可知[20],兩方程湍流模型的一個共同缺陷是過早地預測出再附著和過晚地預測出分離。究其機理,即是湍流模型對湍動能耗散率預測不足,下面說明其具體改進過程。

    根據(jù)湍動能的定義

    可得

    (2k)2=

    引入剪切應力的影響,由雷諾應力的定義可得

    可得不等式

    根據(jù)Boussinesq假設,可得

    2ρk≥{2μt[Sij-δij?uk/(3?xk)]-2δijρk/3}2

    進一步變換,可得

    φ=2.5

    綜上所述,對原始低雷諾數(shù)湍流模型作出以下限制:

    ε0是由原始湍流模型得出的湍動能耗散率。這項改進消除了湍流模型對湍動能耗散率的過強依賴性,避免了由原始湍流模型對湍動能耗散率的過低估計而帶來的負面影響。

    1.2近壁面分離區(qū)改進

    以上改進可以提高模型對轉捩的模擬能力,然而,和普通兩方程模型一樣,近壁面分離區(qū)的影響被忽略。根據(jù)兩方程模型的模擬,分離區(qū)與非分離區(qū)速度應變率的大小在同一個量級,然而直接數(shù)值計算結果表明:分離區(qū)應變率應該比非分離區(qū)應變率大很多[20]。在分離區(qū),即使很小的分離泡,擾動也會被快速放大,且湍流生成急劇增加。為了提高模型在近壁面區(qū)模擬分離流以及分離流轉捩的能力,應使湍流黏度根據(jù)彎曲曲率的變化快速作出響應,即曲率一旦增大,湍流黏度迅速響應并減小,使分離提前,再附著遲滯,使轉捩特性更加接近實驗值。

    邊界層分離以大的彎曲曲率為特征,因此流線曲率因子CS被引入作為分離的感應量。CS作為流線S的曲率,可表示為

    其中,CS是流線曲率的大小,也可以被當作是渦的長度尺度,Cx、Cy、Cz是流線S在三個方向偏差的應變率:

    Cx=q·grad(u)-(u/q)(q·grad(q))

    Cy=q·grad(v)-(v/q)(q·grad(q))

    Cz=q·grad(w)-(w/q)(q·grad(q))

    這里我們定義:

    nq=(u/q,v/q,w/q)

    q=(u,v,w)

    CS通過一個量綱一量被引入湍流黏度,此量綱一參數(shù)為

    從物理角度來看,Sl可被當作是湍動能生成項和耗散項的比率。通過參數(shù)α,Sl被引入湍流黏度:

    fμ=exp[-2.5/(1+α+ReT/50)]

    α=C0Rl/(Rl+0.1)

    上式中的0.1是選取的與Rl同數(shù)量級的數(shù)值,可推得

    這樣的改進方法可以減小分離區(qū)內(nèi)湍流黏性,由于湍流黏性作用于湍動能生成項,故可減少分離區(qū)內(nèi)湍流生成,且對非分離區(qū)的湍流黏性幾乎沒有影響。通過這樣的改進方法,在邊界層中,ReT較小,模型具備模擬轉捩的能力。一旦湍流得到充分發(fā)展,ReT劇增,Sl的影響就微乎其微了。

    這種方法結合上述TSL方法即得到改進的對轉捩敏感的低雷諾數(shù)k-ε模型,下文記作new LRNk-ε模型。

    1.3new LRNk-ε模型對轉捩的預測能力

    為證實以上改進的合理性,利用ERCOFTAC(歐洲流動、湍流與燃燒研究協(xié)會)實驗數(shù)據(jù)庫T3A零壓力梯度平板實驗結果[21]驗證 new LRNk-ε對轉捩的預測能力。T3A實驗來流速度及湍流度分別為5.2m/s和3.3%?;趏riginal(原始) LRNk-ε模型和new LRNk-ε模型的計算結果如圖1所示,通過與文獻[22]的實驗數(shù)據(jù)對比,可以看出改進后的模型對轉捩位置及轉捩區(qū)長度的預測均與實驗吻合較好。

    圖1 T3A平板摩擦因數(shù)流向變化曲線

    2模型驗證

    由于Ahmed模型具備接近真實汽車的特征,故其引入該模型以驗證new LRNk-ε模型的準確性,Ahmed模型長1044mm,寬389mm,高288mm,背部斜邊長222mm,斜角為25°,底面距離地面50mm,如圖2所示。圖3為Lienhart等[22]在LSTM(德國流體力學研究所)低速風洞對Ahmed 模型的氣動阻力、表面壓力系數(shù)以及尾部流場進行的測量試驗。

    圖2 Ahmed模型

    圖3 Ahmed風洞試驗(LSTM低速風洞)

    對車輛進行數(shù)值模擬時的外部計算域為長方體,綜合考慮計算效率及阻塞效應的影響,選取計算域?qū)挒?倍模型寬(左右各3個車寬),高為5倍模型高,長為11倍模型長(前面3個車長,后面7個車長)[23]。

    選取new LRNk-ε、original LRNk-ε、可實現(xiàn)k-ε模型(realizablek-ε)、基于剪切應力運輸?shù)姆蛛x渦模型(DES-SST)以及可實現(xiàn)k-ε及大渦模擬混合模型(RANS/LES)對Ahmed模型的外流場進行計算分析。

    選取來流速度為40m/s,計算域上壁面和側壁面邊界設定為自由滑移壁面。計算域下壁面(即地面)設定為移動壁面邊界,移動速度與來流速度相同。汽車模型壁面邊界滿足固壁無滑移條件,出口邊界給定壓力邊界條件,相對壓力為零。

    控制方程的離散采用二階迎風格式,壓力和速度耦合策略為SIMPLE,采用隱式求解器。

    2.1網(wǎng)格無關性

    基于本文提出的低雷諾數(shù)湍流模型,采用三種非結構性網(wǎng)格以驗證網(wǎng)格無關性,網(wǎng)格數(shù)量分別為8.5×106,1.34×107和2.05×107,其第一層網(wǎng)格質(zhì)心到壁面的量綱一距離y+均小于1,滿足低雷諾數(shù)湍流模型對網(wǎng)格的要求。利用DELL Power 32 處理器計算所得的阻力系數(shù)Cd、仿真與實驗的誤差e及計算時間結果如表1所示,由表1可知,1.34×107個網(wǎng)格能在保證計算精度的條件下消耗較少的計算資源,因此,本文的仿真中采用1.34×107個網(wǎng)格。

    表1 各種網(wǎng)格的阻力系數(shù)及響應的計算時間

    2.2結果與分析

    圖4所示為Ahmed模型表面縱對稱面壓力系數(shù)分布,通過與文獻[22]實驗結果及Realizablek-ε等湍流模型對比,可以發(fā)現(xiàn)在模型斜背處,各個湍流模型顯示出較大差異。new LRNk-ε即使是在最復雜的邊緣處依然給出與實驗較為接近的結果,Realizablek-ε模型及RANS/LES模型由

    于引入了壁面函數(shù),而這對復雜的分離流和再附著流來說,必然引入較大的誤差;DES-SST需要計算較多的標量方程,消耗大量計算資源,同時對分離流的預測出現(xiàn)較大波動。由以上結果可知,合理地捕捉轉捩是準確模擬氣流分離及再附著的基礎。

    圖4 表面縱對稱面壓力系數(shù)分布

    如圖5a實驗結果所示,模型尾部出現(xiàn)兩個典型的渦,這里標記為vortex1與vortex2。vortex1是氣流流經(jīng)尾部斜面邊緣時失去附著,形成一股強大的剪切流,并受車尾負壓卷吸而成;同理,vortex2是流經(jīng)車底的剪切流與車尾負壓聯(lián)合作用形成的。圖5將基于不同湍流模型得出的尾部流場與實驗結果進行了比較,結果表明new LRNk-ε模型得到的尾部流場與粒子成像測速測量結果最接近。Realizablek-ε與original LRNk-ε模型計算得到的兩個渦在大小及位置上均與實驗結果有所出入。

    (a)實驗 (b)new LRN k-ε(c)Realizable k-ε

    (d)DES-SST (e)original LRN k-ε(f)RANS/LES

    圖5 模型尾部流場流線圖及壓力系數(shù)分布

    不同湍流模型計算得出的Ahmed模型尾部速度剖面圖見圖6,由new LRNk-ε模型模擬所得來流在頂部邊緣分離,再附著于斜面,與實驗結果相符,original LRNk-ε模型也得到了較為吻合的結果,然而不能準確計算速度大小,DES-SST與RANS/LES模型所得結果與實驗相差較大。

    (a)水平分量

    (b)垂直分量

    表6Ahmed模型尾部速度剖面圖

    表2所示為計算與實驗得到的Ahmed模型阻力系數(shù),從表中可以看出new LRNk-ε模型得到的結果與實驗值之間的誤差最小。

    表2 Ahmed模型的阻力系數(shù)

    圖7所示為基于new LRNk-ε模型且由ISO壓力曲面描繪出的渦結構,車尾左右各有一個旋轉方向相反的拖曳縱向旋轉渦,將其記為counter-rotating vortex pair。該旋轉渦是氣流流經(jīng)后側圍邊緣時失去附著,形成一股強大的剪切流,并受車尾負壓卷吸而成的,該漩渦一方面繞自身渦核旋轉,另一方面受車頂高速剪切流的影響向著車尾方向拖曳而去,渦核位置隨著拖曳距離的增加向著地面靠近。counter-rotating vortex pair下方為類似于符號“Ω”的渦,與圖5尾部流線圖相符,記為Ω shape vortex。模型尾流中靠近地面為一馬蹄形渦,記為horseshoe vortex,其形成過程為:counter-rotating vortex pair的兩個縱向漩渦相互作用并靠近,渦的不穩(wěn)定性使其中間上部流體形成低壓區(qū),下部流體因積聚使流速減慢形成高壓區(qū)。上下部分不同流速的流體相互剪切造成一個與counter-rotating vortex pair相連接的橫向渦,從而形成馬蹄形渦,模型底部支柱后的渦街也被較好地捕捉。綜上所述, new LRNk-ε模型能較為準確地捕捉Ahmed模型的流場信息。

    圖7 Ahmed模型ISO壓力曲面 (基于new LRN k-ε模型)

    3結束語

    針對原始低雷諾數(shù)k-ε模型無法準確預測轉捩的問題,引入了TSL方法及流線曲率因子以提高模型對轉捩的敏感性,進而獲得了能準確預測轉捩的new LRNk-ε模型,并應用于汽車外部流場計算。

    在相同的網(wǎng)格條件下,相比其他湍流模型,本文提出的new LRNk-ε模型能更加準確地模擬汽車外部氣流的分離以及再附著、縱對稱面壓力系數(shù)、汽車氣動阻力系數(shù)及尾部流場分布,提高了計算精度,為類似于汽車外部流場的復雜流動的仿真分析打下了堅實的基礎。

    參考文獻:

    [1]谷正氣. 汽車空氣動力學[M].北京:人民交通出版社,2005.

    [2]傅立敏, 吳允柱, 賀寶琴. 超車過程的車輛氣動特性仿真研究[J]. 中國機械工程, 2007,18(5):621-625.

    Fu Limin, Wu Yunzhu, He Baoqin. Simulation Research on Automotive Aerodynamic Characteristics during the Overtaking Process[J].China Mechanical Engineering,2007,18(5):621-625.

    [3]容江磊, 谷正氣, 楊易, 等. 基于kriging模型的跑車尾翼斷面形狀的氣動優(yōu)化[J]. 中國機械工程, 2011,22(2):243-247.

    Rong Jianglei,Gu Zhengqi,Yang Yi,et al.Aerodynamic Optimization of Cross Sectional Shape for a Sports Car’s Rear Wing Based on Kriging Surrogate Model[J].China Mechanical Engineering,2011,22(2):243-247.

    [4]楊易,徐永康,聶云,等.非定常來流對汽車氣動升力瞬態(tài)特性的影響[J].中國機械工程,2014,25(12):1681-1686.

    Yang Yi,Xu Yongkang,Nie Yun,et al.Effects of Unsteady Stream on Transient Characteristics of Automotive Lift[J].China Mechanical Engineering, 2014, 25(12):1681-1686.

    [5]何憶斌, 李偉平, 劉孟祥,等. F1方程式賽車后升力翼的氣動特性[J]. 航空動力學報, 2013, 28(10):2343-2347.

    He Yibin,Li Weiping,Liu Mengxiang,et al.Rear Lifting Airfoil Aerodynamic Characteristic of F1 Racing Car[J]. Journal of Aerospace Power, 2013, 28(10):2343-2347.

    [6]田紅旗.中國列車空氣動力學研究進展[J].交通運輸工程學報,2006,6(1):1-9.

    Tian Hongqi.Study Evolvement of Train Aerodynamics in China[J].Journal of Traffic and Transportation Engineering,2006,6(1):1-9.

    [7]Versteeg H K,Malalasekera W.An Introduction to Computational Fluid Dynamics:The Finite Volume Method(2nd Edition)[M].New York: Prentice Hall,2007.

    [8]汪怡平,谷正氣,鄧亞東.基于準k-ε-v2/LES模型的汽車外流場數(shù)值模擬[J].機械工程學報,2012,48(14):97-103.

    Wang Yiping,Gu Zhengqi,Deng Yadong.Aerodynamic Simulations of Vehicle by Using the Hybrid Semik-ε-v2/LES Model[J].Journal of Mechanical Engineering,2012,48(14):97-103.

    [9]Jones W P,Launder B E.The Calculation of Low-Reynolds-number Phenomena with a Two-equation Model of Turbulence[J].International Journal of Heat and Mass Transfer,1973,16:1119-1130.

    [10]Abe K,Kondoh T,Nagano Y.A New Turbulence Model for Predicting Fluid Flow and Heat Transfer in Separating and Reattaching Fflows I:Low Fleld Calculations[J].International Journal of Heat and Mass Transfer,1994,37(1):139-151.

    [11]Chang K C,Hsieh W D,Chen C S.A Modified Low-Reynolds-number Turbulence Model Applicable to Recirculating Flow in Pipe Expansion[J].Transactions of the ASME,Journal of Fluids Engineering,1995,117: 417-423.

    [12]符松,王亮.湍流轉捩模式研究進展[J].力學進展,2007,37(3):409-416.

    Fu Song,Wang Liang.Progress in Turbulence/Transition Modeling[J].Advances in Mechanics,2007,37(3):409-416.

    [13]Menter F,Langtry R,Likki S,et al.Correlatio-based Transition Model Using Local Variables-Part 1:Model Formulation[J].Journal of Turbomachinery,2006,128(3):413-422.

    [14]Mayle R E, Schulz A.The Path to Predicting Bypass Transtion[J].Journal of Turbomachinery, 1997,119(3):405-411.

    [15]Cutrone L,Palma P D,Pascazio G,et al.An Evaluation of Bypass Transiton Models for Turbomachinery Flows[J].Heat and Fluid Flow,2006,28(1):161-177.

    [16]Walters D K,James H L.A New Model for Boundary Layer Transition Using a Single-point RANS Approach[J].Journal of Turbomachinery,2004,126(1):193-202.

    [17]Walters D K,Cokljat D.A Three Equation Eddy Viscosity Model for Reynolds Averaged Navie-stokes Simulations of Transitional Flow[J].Journal of Fluids Engineering,2008,130(12):121401.1-14.

    [18]Walters D K.Physical Interpretation of Transition-sensitive RANS Models Employing the Laminar Kinetic Energy Concept[J].ERCOFTAC Bulletin,2009,80:61-67.

    [19]Taghavi-Zenouz R,Salari M,Etemadi M.Prediction of Laminar,Transitional and Turbulent Flow Regimes,Based on Three Equationk-ωTurbulence Model[J].The Aeronautical Journal,2008,112(1134):469-476.

    [20]Zheng X,Liu C,Liu F,et al.Turbulent Transition Simulation Using thek-ωModel[J]. International Journal for Numerical Methods in Engineering, 1996,42(10):907-926.

    [21]Coupland J.Special Interest Group on Laminar to Turbulence Transition and Re Transition,T3A and T3B Test Cases[R].London:ERCOFTAC,1990.

    [22]Lienhart H,Becker S.Flow and Turbulence in the Wake of A Simplified Car Model[R].SAE,2003-01-0656.

    [23]李學武.某微車的氣動特性分析及優(yōu)化[D].長沙:湖南大學,2008.

    (編輯袁興玲)

    猜你喜歡
    風洞試驗
    羽毛球風洞試驗
    直升機前飛狀態(tài)旋翼結冰風洞試驗研究
    飛翼布局飛機陣風減緩主動控制風洞試驗
    低風壓架空導線的風洞試驗
    電線電纜(2017年5期)2017-10-18 00:52:03
    滾轉機動載荷減緩風洞試驗
    特種風洞試驗中氣動伺服彈性失穩(wěn)故障分析
    弧形內(nèi)凹大跨屋蓋結構風荷載特性的風洞試驗與數(shù)值模擬
    民機高速風洞試驗的阻力雷諾數(shù)效應修正
    飛機全動平尾顫振特性風洞試驗
    航空學報(2015年4期)2015-05-07 06:43:33
    港口地貌大氣邊界層風洞試驗模擬
    水道港口(2015年1期)2015-02-06 01:25:44
    国产欧美日韩一区二区精品| 亚洲人与动物交配视频| 国产高清三级在线| 一个人看的www免费观看视频| 国产精品国产高清国产av| 亚洲一区二区三区色噜噜| 少妇被粗大猛烈的视频| 欧美一区二区精品小视频在线| 欧美一区二区亚洲| 亚洲专区中文字幕在线| 白带黄色成豆腐渣| 欧美日本视频| 色综合婷婷激情| av专区在线播放| 又爽又黄a免费视频| aaaaa片日本免费| 精品熟女少妇八av免费久了| 国模一区二区三区四区视频| 国产精品自产拍在线观看55亚洲| 亚洲激情在线av| 成人三级黄色视频| 精品人妻1区二区| 1000部很黄的大片| 日本黄色视频三级网站网址| 在线天堂最新版资源| 国产真实伦视频高清在线观看 | 少妇被粗大猛烈的视频| 亚洲黑人精品在线| www.熟女人妻精品国产| 亚洲精品一卡2卡三卡4卡5卡| 一级a爱片免费观看的视频| 淫秽高清视频在线观看| 国产精品亚洲av一区麻豆| 国产伦精品一区二区三区四那| 国产精品1区2区在线观看.| 婷婷色综合大香蕉| 夜夜爽天天搞| 国产在视频线在精品| 国产乱人伦免费视频| 中文字幕人成人乱码亚洲影| 免费大片18禁| 久久精品国产99精品国产亚洲性色| 波多野结衣巨乳人妻| 一进一出抽搐gif免费好疼| 色吧在线观看| 欧美xxxx黑人xx丫x性爽| 亚洲精品456在线播放app | 成人国产一区最新在线观看| 久久久精品欧美日韩精品| 禁无遮挡网站| 高清毛片免费观看视频网站| 老司机深夜福利视频在线观看| 国产v大片淫在线免费观看| 久久久国产成人精品二区| 日本黄色片子视频| 美女高潮的动态| 极品教师在线视频| 桃色一区二区三区在线观看| 欧美黄色淫秽网站| 男人狂女人下面高潮的视频| 直男gayav资源| 啦啦啦观看免费观看视频高清| 亚洲三级黄色毛片| 国产精品乱码一区二三区的特点| 国产高潮美女av| 长腿黑丝高跟| 成人午夜高清在线视频| 精品熟女少妇八av免费久了| 久久伊人香网站| 成人一区二区视频在线观看| 欧美+日韩+精品| 每晚都被弄得嗷嗷叫到高潮| 亚洲国产精品sss在线观看| 日韩欧美三级三区| av视频在线观看入口| 三级毛片av免费| 如何舔出高潮| 天堂网av新在线| 国产色爽女视频免费观看| 日日夜夜操网爽| 国产亚洲精品综合一区在线观看| 色5月婷婷丁香| 变态另类成人亚洲欧美熟女| 精品国产亚洲在线| 久久精品人妻少妇| 丁香欧美五月| 精品午夜福利在线看| 一二三四社区在线视频社区8| 国产精品精品国产色婷婷| 亚洲欧美激情综合另类| 精品福利观看| 久久久久久久久大av| 国产成人影院久久av| 久久国产精品人妻蜜桃| 日本免费一区二区三区高清不卡| 午夜亚洲福利在线播放| 高清在线国产一区| 精品一区二区三区视频在线观看免费| 黄片小视频在线播放| АⅤ资源中文在线天堂| 日日摸夜夜添夜夜添小说| 亚洲av电影不卡..在线观看| 国产 一区 欧美 日韩| 精品国内亚洲2022精品成人| 精品久久久久久久人妻蜜臀av| 午夜激情福利司机影院| 特大巨黑吊av在线直播| 国产欧美日韩精品一区二区| 麻豆成人av在线观看| 国产午夜精品久久久久久一区二区三区 | 搡女人真爽免费视频火全软件 | 国产精品久久视频播放| 一本一本综合久久| aaaaa片日本免费| 一本精品99久久精品77| 亚洲av成人精品一区久久| 成人精品一区二区免费| 欧美激情久久久久久爽电影| 动漫黄色视频在线观看| 国产精品久久久久久亚洲av鲁大| 婷婷色综合大香蕉| 精品乱码久久久久久99久播| 亚洲狠狠婷婷综合久久图片| 欧美成狂野欧美在线观看| 久久草成人影院| 国产白丝娇喘喷水9色精品| 色综合亚洲欧美另类图片| 一本一本综合久久| 久久精品国产亚洲av涩爱 | 啦啦啦韩国在线观看视频| 亚洲 国产 在线| 18禁在线播放成人免费| 国产麻豆成人av免费视频| 免费黄网站久久成人精品 | 亚洲自拍偷在线| 91麻豆精品激情在线观看国产| 琪琪午夜伦伦电影理论片6080| 亚洲av电影在线进入| 亚洲五月天丁香| 能在线免费观看的黄片| 久久久久国产精品人妻aⅴ院| 久久这里只有精品中国| 国产人妻一区二区三区在| 久久久久免费精品人妻一区二区| 久久精品国产自在天天线| 中国美女看黄片| 国产色婷婷99| www日本黄色视频网| 成人无遮挡网站| 精品国内亚洲2022精品成人| 搡老岳熟女国产| www.www免费av| 日韩大尺度精品在线看网址| 色哟哟哟哟哟哟| 亚洲片人在线观看| 国产国拍精品亚洲av在线观看| 久久精品综合一区二区三区| 亚洲人与动物交配视频| 最新在线观看一区二区三区| 亚洲在线观看片| 欧美乱色亚洲激情| 日韩免费av在线播放| www.色视频.com| 我要看日韩黄色一级片| 国产亚洲精品久久久com| 别揉我奶头 嗯啊视频| 内射极品少妇av片p| 日韩精品青青久久久久久| 亚洲专区中文字幕在线| 国产精品久久久久久亚洲av鲁大| 久久香蕉精品热| 成人av一区二区三区在线看| 黄色日韩在线| 久久久久亚洲av毛片大全| 亚洲真实伦在线观看| 国产精品免费一区二区三区在线| 18禁在线播放成人免费| 亚洲人成电影免费在线| 亚洲美女搞黄在线观看 | 婷婷六月久久综合丁香| 成人午夜高清在线视频| 日韩亚洲欧美综合| 99热这里只有精品一区| 啦啦啦观看免费观看视频高清| 亚洲无线观看免费| 在线a可以看的网站| 国产精品久久久久久亚洲av鲁大| 日韩欧美三级三区| 亚洲精品影视一区二区三区av| 在线免费观看不下载黄p国产 | 成人永久免费在线观看视频| 少妇人妻一区二区三区视频| 禁无遮挡网站| 白带黄色成豆腐渣| 听说在线观看完整版免费高清| www.熟女人妻精品国产| 神马国产精品三级电影在线观看| 国内毛片毛片毛片毛片毛片| xxxwww97欧美| 免费观看精品视频网站| 特级一级黄色大片| 99久久九九国产精品国产免费| 日韩欧美国产一区二区入口| 少妇丰满av| 国产精品一区二区三区四区久久| 国模一区二区三区四区视频| 99久久精品国产亚洲精品| 日本在线视频免费播放| 9191精品国产免费久久| 我的老师免费观看完整版| 国产91精品成人一区二区三区| 成人无遮挡网站| a级一级毛片免费在线观看| 日韩中字成人| 美女 人体艺术 gogo| 国产精品,欧美在线| 精品一区二区三区视频在线| 女同久久另类99精品国产91| 综合色av麻豆| 中文字幕免费在线视频6| 免费看光身美女| 精品人妻1区二区| 丰满乱子伦码专区| 日本一二三区视频观看| 性欧美人与动物交配| 好看av亚洲va欧美ⅴa在| 欧美区成人在线视频| 精品人妻视频免费看| 亚洲狠狠婷婷综合久久图片| 无人区码免费观看不卡| 丁香六月欧美| 久久精品夜夜夜夜夜久久蜜豆| 能在线免费观看的黄片| 一区二区三区激情视频| 久久6这里有精品| 99热精品在线国产| 99热只有精品国产| 一级av片app| 一进一出抽搐gif免费好疼| 久久6这里有精品| 波多野结衣巨乳人妻| 久久久成人免费电影| 非洲黑人性xxxx精品又粗又长| 波野结衣二区三区在线| 岛国在线免费视频观看| 国产伦在线观看视频一区| 在线免费观看的www视频| 一级av片app| 老熟妇仑乱视频hdxx| 一个人观看的视频www高清免费观看| 少妇的逼水好多| 亚洲激情在线av| 两人在一起打扑克的视频| 午夜两性在线视频| 久久久久久久精品吃奶| 两性午夜刺激爽爽歪歪视频在线观看| 成人特级av手机在线观看| 亚洲第一电影网av| 亚洲乱码一区二区免费版| 国产亚洲精品久久久com| 欧美日韩瑟瑟在线播放| 制服丝袜大香蕉在线| 99热精品在线国产| 精品久久久久久久久亚洲 | 久久精品国产亚洲av香蕉五月| 欧美性感艳星| 欧美色欧美亚洲另类二区| a级一级毛片免费在线观看| 欧美日韩黄片免| 久久久久久久久大av| 狠狠狠狠99中文字幕| 日韩 亚洲 欧美在线| 99热这里只有是精品50| 全区人妻精品视频| 婷婷六月久久综合丁香| 一边摸一边抽搐一进一小说| 一本精品99久久精品77| 成人av一区二区三区在线看| 欧美性猛交╳xxx乱大交人| 欧美乱妇无乱码| 美女黄网站色视频| 非洲黑人性xxxx精品又粗又长| 一个人看视频在线观看www免费| 久久精品夜夜夜夜夜久久蜜豆| 国产成人影院久久av| 在线国产一区二区在线| 在线观看一区二区三区| 久久午夜亚洲精品久久| 亚洲中文日韩欧美视频| 日韩中字成人| 亚州av有码| 欧美乱妇无乱码| 日韩av在线大香蕉| 国产精品野战在线观看| 欧美日本亚洲视频在线播放| 色播亚洲综合网| 两人在一起打扑克的视频| 一个人看视频在线观看www免费| 亚洲片人在线观看| 色精品久久人妻99蜜桃| 中文在线观看免费www的网站| 亚洲中文字幕一区二区三区有码在线看| 久久久久久久久久黄片| 禁无遮挡网站| 99热精品在线国产| 神马国产精品三级电影在线观看| 日日干狠狠操夜夜爽| 别揉我奶头 嗯啊视频| 精品人妻视频免费看| x7x7x7水蜜桃| 欧美+亚洲+日韩+国产| 好男人在线观看高清免费视频| 国产老妇女一区| 搡老妇女老女人老熟妇| 久久伊人香网站| 熟女电影av网| 久久精品综合一区二区三区| 国产精品亚洲av一区麻豆| 少妇高潮的动态图| 欧美黄色片欧美黄色片| 欧美乱妇无乱码| a级一级毛片免费在线观看| 色在线成人网| 丝袜美腿在线中文| 亚洲第一电影网av| 欧美一区二区精品小视频在线| 成人特级av手机在线观看| 中出人妻视频一区二区| 黄色日韩在线| 亚洲熟妇熟女久久| 亚洲av成人精品一区久久| 国产欧美日韩一区二区三| 欧美成人免费av一区二区三区| 精品人妻1区二区| 欧美xxxx性猛交bbbb| 亚洲国产高清在线一区二区三| 中亚洲国语对白在线视频| 少妇高潮的动态图| 天天躁日日操中文字幕| 波野结衣二区三区在线| 久久久久国内视频| 夜夜看夜夜爽夜夜摸| 欧美黑人巨大hd| 老熟妇乱子伦视频在线观看| 欧美激情在线99| 国产精品乱码一区二三区的特点| 国产精品女同一区二区软件 | 韩国av一区二区三区四区| 亚洲无线在线观看| 人人妻人人看人人澡| 国产黄a三级三级三级人| 啦啦啦观看免费观看视频高清| 88av欧美| 国产精品电影一区二区三区| 少妇的逼好多水| 人妻丰满熟妇av一区二区三区| 久久中文看片网| 直男gayav资源| 少妇丰满av| 欧美成人免费av一区二区三区| 亚洲国产精品999在线| 夜夜夜夜夜久久久久| 免费看a级黄色片| xxxwww97欧美| 丁香欧美五月| 欧美高清性xxxxhd video| 精华霜和精华液先用哪个| 中出人妻视频一区二区| 国产主播在线观看一区二区| 免费在线观看影片大全网站| 国产亚洲精品久久久com| 亚洲第一区二区三区不卡| 国产精华一区二区三区| 欧美日韩黄片免| 国产精品嫩草影院av在线观看 | 亚洲成人久久性| 日韩亚洲欧美综合| 好男人电影高清在线观看| 简卡轻食公司| 精品人妻1区二区| 日本黄色片子视频| 久久精品国产自在天天线| 亚洲人成网站在线播| 国产午夜精品久久久久久一区二区三区 | 波多野结衣高清无吗| 99热6这里只有精品| av天堂在线播放| 久久久久性生活片| 亚洲精品成人久久久久久| 中出人妻视频一区二区| 99久久精品热视频| 淫妇啪啪啪对白视频| 欧美黑人欧美精品刺激| 99精品在免费线老司机午夜| 久久精品综合一区二区三区| 一本综合久久免费| 国产精品自产拍在线观看55亚洲| 国产aⅴ精品一区二区三区波| 好男人电影高清在线观看| 国产精品久久视频播放| 国产av麻豆久久久久久久| 亚洲av.av天堂| 又爽又黄无遮挡网站| 亚洲中文字幕日韩| or卡值多少钱| 色尼玛亚洲综合影院| 国产av在哪里看| a在线观看视频网站| netflix在线观看网站| 我的女老师完整版在线观看| 日韩 亚洲 欧美在线| h日本视频在线播放| 精品久久久久久久久av| 欧美3d第一页| 在线观看av片永久免费下载| 国产亚洲欧美98| 五月伊人婷婷丁香| 在线播放国产精品三级| 黄片小视频在线播放| 男女那种视频在线观看| 免费高清视频大片| 国产精品野战在线观看| 少妇裸体淫交视频免费看高清| 国产伦一二天堂av在线观看| 久久6这里有精品| 变态另类成人亚洲欧美熟女| 97热精品久久久久久| 午夜两性在线视频| 一进一出抽搐动态| 久久人人精品亚洲av| 免费看日本二区| 国产一区二区在线观看日韩| 国产大屁股一区二区在线视频| 亚洲av免费在线观看| 国产高清三级在线| 国产三级黄色录像| netflix在线观看网站| 久久久精品欧美日韩精品| 午夜视频国产福利| 日日干狠狠操夜夜爽| 久久国产精品影院| 亚洲欧美日韩高清专用| 最好的美女福利视频网| 欧美另类亚洲清纯唯美| 五月玫瑰六月丁香| 亚洲av中文字字幕乱码综合| 男女之事视频高清在线观看| 免费av不卡在线播放| 男人和女人高潮做爰伦理| 好男人电影高清在线观看| 又爽又黄无遮挡网站| 亚洲美女视频黄频| 精品久久久久久久久久免费视频| 亚洲av日韩精品久久久久久密| 特级一级黄色大片| 日韩高清综合在线| 成人国产一区最新在线观看| 最新在线观看一区二区三区| 精品一区二区免费观看| 看免费av毛片| 亚洲av一区综合| .国产精品久久| 最近在线观看免费完整版| 又粗又爽又猛毛片免费看| 伦理电影大哥的女人| 狠狠狠狠99中文字幕| 国产视频一区二区在线看| 欧美日韩福利视频一区二区| 美女免费视频网站| av黄色大香蕉| 亚洲av成人av| 欧美绝顶高潮抽搐喷水| 18美女黄网站色大片免费观看| 悠悠久久av| xxxwww97欧美| 亚洲美女视频黄频| 美女被艹到高潮喷水动态| 亚洲av一区综合| 一级黄色大片毛片| 高清日韩中文字幕在线| 亚洲人成电影免费在线| 在线观看免费视频日本深夜| 我要看日韩黄色一级片| 日韩大尺度精品在线看网址| 国产精品,欧美在线| 精品午夜福利视频在线观看一区| 国内精品美女久久久久久| 久久亚洲真实| 日韩国内少妇激情av| 超碰av人人做人人爽久久| 在线播放无遮挡| 波多野结衣高清作品| 欧美一级a爱片免费观看看| 五月玫瑰六月丁香| 怎么达到女性高潮| 少妇被粗大猛烈的视频| 精品久久久久久久久亚洲 | 欧美午夜高清在线| 又爽又黄无遮挡网站| 嫩草影视91久久| 熟女人妻精品中文字幕| 欧美又色又爽又黄视频| 非洲黑人性xxxx精品又粗又长| 2021天堂中文幕一二区在线观| 99久久精品国产亚洲精品| 亚洲成人中文字幕在线播放| 首页视频小说图片口味搜索| 色在线成人网| 国产成人啪精品午夜网站| 久久这里只有精品中国| 国产久久久一区二区三区| 日韩欧美 国产精品| 国产精华一区二区三区| 99视频精品全部免费 在线| 最近最新免费中文字幕在线| 欧美绝顶高潮抽搐喷水| 午夜亚洲福利在线播放| 久久久久久久久久成人| 国产一区二区在线av高清观看| 国产精品永久免费网站| av国产免费在线观看| 欧美黄色淫秽网站| 亚洲自拍偷在线| 亚洲第一欧美日韩一区二区三区| 亚洲美女黄片视频| 桃色一区二区三区在线观看| 中文字幕av成人在线电影| 1024手机看黄色片| 久久国产精品影院| 嫩草影院入口| 欧美不卡视频在线免费观看| 欧美zozozo另类| 在线十欧美十亚洲十日本专区| 51午夜福利影视在线观看| 亚洲国产日韩欧美精品在线观看| 国产免费男女视频| 在线观看美女被高潮喷水网站 | 午夜影院日韩av| 99久久精品国产亚洲精品| 午夜福利18| 国产在视频线在精品| 国产精品精品国产色婷婷| 三级男女做爰猛烈吃奶摸视频| 国产美女午夜福利| 久久精品国产自在天天线| 欧美潮喷喷水| 日本 欧美在线| 国产精品影院久久| 尤物成人国产欧美一区二区三区| 黄色视频,在线免费观看| 国产精品久久久久久久电影| 国产成人欧美在线观看| 哪里可以看免费的av片| 亚洲av五月六月丁香网| 国产精品美女特级片免费视频播放器| 色噜噜av男人的天堂激情| 脱女人内裤的视频| 欧美zozozo另类| 国产乱人视频| 国产在视频线在精品| 免费大片18禁| 久久香蕉精品热| 又爽又黄无遮挡网站| 久久人人爽人人爽人人片va | 精品久久久久久久久av| 亚洲,欧美,日韩| 十八禁人妻一区二区| av中文乱码字幕在线| 午夜久久久久精精品| 制服丝袜大香蕉在线| 1000部很黄的大片| 欧美日韩亚洲国产一区二区在线观看| 色综合站精品国产| 国产白丝娇喘喷水9色精品| 国内精品久久久久久久电影| 国产91精品成人一区二区三区| 特级一级黄色大片| 国产人妻一区二区三区在| 久久婷婷人人爽人人干人人爱| 少妇的逼水好多| 桃红色精品国产亚洲av| 深夜a级毛片| av专区在线播放| 欧美成人一区二区免费高清观看| 最好的美女福利视频网| 欧美一级a爱片免费观看看| 日韩av在线大香蕉| 婷婷精品国产亚洲av在线| 亚洲在线自拍视频| 一个人观看的视频www高清免费观看| 永久网站在线| 久久久久性生活片| 一区福利在线观看| 又爽又黄a免费视频| 亚洲熟妇中文字幕五十中出| 色综合婷婷激情| а√天堂www在线а√下载| 国产精品98久久久久久宅男小说| 国产亚洲精品久久久com| 神马国产精品三级电影在线观看| 夜夜夜夜夜久久久久| 午夜免费成人在线视频| 99久国产av精品| 国产91精品成人一区二区三区| 亚洲欧美日韩高清在线视频| 国产一区二区在线av高清观看| 日日摸夜夜添夜夜添av毛片 | 中亚洲国语对白在线视频| av在线观看视频网站免费| 国产探花在线观看一区二区| 久久午夜福利片| 亚洲乱码一区二区免费版| 成人欧美大片| 午夜免费激情av| 五月伊人婷婷丁香|