亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        基于DDES算法的有擾流片腔體氣動噪聲分析

        2015-11-09 00:51:36童明波ZhiweiHu
        空氣動力學(xué)學(xué)報 2015年5期
        關(guān)鍵詞:流片空腔腔體

        劉 瑜,童明波,*,Zhiwei Hu

        (1.南京航空航天大學(xué) 航空宇航學(xué)院,江蘇 南京 210016;2.University of Southampton,Southampton UK SO17 1BJ)

        0 引 言

        隨著戰(zhàn)斗機和無人機的內(nèi)埋彈艙的應(yīng)用,腔體結(jié)構(gòu)的氣動噪聲問題成為研究熱點。腔體的聲學(xué)共振現(xiàn)象[1-2]產(chǎn)生于腔體前緣剪切層的不穩(wěn)定增長,振蕩的高速氣流撞擊到腔體后壁產(chǎn)生壓強波形成聲源,壓強波向前傳播并與前緣剪切層進行能量交換增強了自由剪切層的不穩(wěn)定性,將在腔體后緣再次引發(fā)強烈的氣流撞擊,從而形成周期性的反饋回路[3]。當(dāng)反饋壓強波的頻率與相位同自由剪切層一致時會導(dǎo)致共振[4]。腔體的這種自持振蕩引發(fā)的氣動噪聲量級高達170dB,會產(chǎn)生聲疲勞和聲輻射,造成結(jié)構(gòu)破壞和艙內(nèi)電子儀器失靈[5]。

        國外關(guān)于開式腔體的氣動噪聲分析從1966年Rossiter[6]進行方形腔體噪聲試驗和半經(jīng)驗公式開始到現(xiàn)在已經(jīng)發(fā)展十分成熟,現(xiàn)階段主要研究重點在腔體結(jié)構(gòu)的噪聲抑制方法上,主要包括擾流棒[7-8]、前緣擾流片[9]、質(zhì)量噴流[10]等方法。國內(nèi)從1996年羅柏華等[11]開始關(guān)注腔體噪聲問題,現(xiàn)階段研究主要關(guān)于噪聲機理和計算方法。司海清、王同光[12]等對空腔振蕩頻率估算方程進行了改進,郝宗瑞[13]、馬明生[14]和宋文萍等對空腔噪聲的計算做了研究。噪聲控制方法研究尚屬起步階段,張林[15]、羅新福[16]、楊黨國[17]、范召林等采用試驗的方法研究了腔體振蕩影響因素以及質(zhì)量噴流和后壁傾角等控制手段對腔體振蕩的抑制效果。國外關(guān)于前緣擾流片的噪聲控制方法研究已經(jīng)成熟并被廣泛應(yīng)用在F-22和F-111[18]等飛機上,但國內(nèi)關(guān)于前緣擾流片等被動控制方法的研究成果較少。

        為了使簡單有效的擾流片被動控制方法應(yīng)用于國內(nèi)內(nèi)埋彈艙設(shè)計中,本文基于DDES的CFD方法來對安裝有前緣鋸齒形擾流片的腔體進行流場分析和聲學(xué)計算,以驗證這種被動控制方法對腔體氣動噪聲的抑制效果,并與試驗數(shù)據(jù)進行了對比。

        1 DDES方法

        單方程Spalart-Allmaras(S-A)[19]湍流模型可表示為:

        耗散項表示為:

        源項定義為:

        基于S-A單方程湍流模型的混合RANS/LES方法本質(zhì)上是將方程(2)耗散項中的距離用下式代替:

        其中,CDES是常數(shù),常取為0.65。Lg是網(wǎng)格尺度,定義為:

        Δx、Δy、Δz為某網(wǎng)格點三個方向上的網(wǎng)格尺寸。因此,當(dāng)網(wǎng)格在物面附近時=d,這時流場求解使用基于S-A單方程湍流模型的RANS方法。當(dāng)網(wǎng)格遠離物面時=CDESLg,這時湍流渦粘性系數(shù)的衰減便由當(dāng)?shù)氐木W(wǎng)格尺寸確定。通過式(1)可知,當(dāng)源項的湍流渦流性系數(shù)與耗散項的衰減達到平衡時,這時的渦粘性系數(shù)與2成正比關(guān)系,此時=CDESLg,則有:

        變成了Smagorinsky[20]亞格子模型。即在遠離物面的地方,流場求解表現(xiàn)為LES所需要的亞格子模型。這種混合RANS/LES的方法稱為分離渦模擬(Detached Eddy Simulation,DES)[21]方法。

        由于在物面附近平行于物面方向的網(wǎng)格尺寸通常大于邊界層的厚度,依據(jù)式(4)的判定條件可以保證在邊界層內(nèi)使用RANS算法。但是,若物面網(wǎng)格各個方向上都較為細密則可能使LES算法在邊界層內(nèi)被提前啟動。這種情況下,邊界層內(nèi)的網(wǎng)格是無法滿足LES計算要求的,并且計算效率將會下降。因此,一種稱為延遲分離渦模擬(Delayed Detached Eddy Simulation,DDES)的改進DES算法被提出來以保證邊界層內(nèi)完全使用RANS算法[22]。

        DDES算法是將式(4)的判據(jù)做如下改動:

        其中:

        rd被視為是當(dāng)?shù)赝牧鞒叨扰c距離壁面距離的比值。在遠離壁面區(qū)域rd?1,使得fd=1,該區(qū)域使用LES算法。其它區(qū)域中fd=0,使用RANS算法。不難想象,當(dāng)渦粘度系數(shù)從數(shù)值較大的區(qū)域轉(zhuǎn)向相對較小的區(qū)域時,將使得物面邊界層附近的DDES模型從LES模式向RANS轉(zhuǎn)變。

        2 試驗與計算

        用來與計算結(jié)果做比較的試驗由Ross、Foster等[23]于1991年11月在Bedford英國皇家飛機研究院的ARA風(fēng)洞中進行[24]。試驗腔體長深比L/D=5,寬深比W/D=1。腔體模型的底部均布了10個Kulite壓強傳感器來測量腔體內(nèi)部的非定常壓強變化。試驗的采樣頻率為6kHz,時長3.2s。試驗臺、腔體與擾流片的尺寸與結(jié)構(gòu)如圖1所示。

        計算模型是基于實驗臺基本尺寸向上延伸10倍腔體深度。上邊界和流動方向的前后邊界采用壓強遠場邊界條件,兩側(cè)邊界采用對稱條件,下邊界為無滑移的壁面。來流馬赫數(shù)Ma=0.85,壓強P=62 940Pa,溫度T=270.25K,渦粘比μt/μ0=10。為了準確捕捉到產(chǎn)生聲源的渦結(jié)構(gòu),腔體內(nèi)部以及開口剪切層位置處網(wǎng)格布置十分細密,第一層網(wǎng)格y+值控制在2左右,計算網(wǎng)格總數(shù)為在450萬左右。腔體與擾流片計算網(wǎng)格如圖2所示。

        圖1 試驗腔體與計算擾流片尺寸(mm)Fig.1 Cavity and spoiler model(unit:mm)

        圖2 計算腔體與擾流片網(wǎng)格Fig.2 Mesh for cavity and spoiler

        計算基于Fluent計算平臺,采用二階隱式時間積分的格心有限體積法,無粘通量采用Roe-FDS的通量差分分裂格式,粘性通量采用三階MUSCL空間離散格式。計算時間步長為10-5s,每步進行30次迭代計算,總計算時間持續(xù)了1s。在進行瞬態(tài)計算之前首先進行基于S-A單方程模型的穩(wěn)態(tài)計算直到流場收斂,前0.5s時間的瞬態(tài)計算舍去,最后0.5s的流場信息以100kHz的采樣頻率做聲學(xué)分析。

        計算在英國University of Southampton的Iridis3服務(wù)器上進行,采用36個4核2.27GHz的處理器來進行腔體流場瞬態(tài)計算??涨缓蛶_流片腔體均在30 000步以上的穩(wěn)態(tài)計算后開始瞬態(tài)計算,空腔工況用時55天,帶擾流片腔體用時57天完成。

        3 結(jié)果與分析

        圖3 腔底壓強均方根結(jié)果比較Fig.3 Comparison of experimentand DDES on prms

        圖3給出了DDES計算空腔結(jié)果與試驗結(jié)果的對比,以及計算得到的安裝有擾流片的腔體底面10個壓力測量點的壓強均方根值。

        腔底壓強的波動是腔體自激振蕩的現(xiàn)象,如前文所介紹,剪切層撞擊腔體后壁面使得腔體后部的渦嚴重振蕩導(dǎo)致壓強變化劇烈,如圖3所示。通過與試驗進行對比,本文所用的DDES方法很好的計算出了空腔底部壓強的振蕩變化規(guī)律。安裝前緣鋸齒形擾流片的DDES結(jié)果表明腔體內(nèi)部的壓強波動得到了較大的抑制。其抑制機理是前緣擾流片可以將剪切層向上抬離腔體上方,如圖4(b)所示,使得剪切層在后緣撞擊腔體后壁的強度和進入腔體的流量大為降低。高速氣流流過鋸齒形擾流片產(chǎn)生許多小尺度的渦結(jié)構(gòu),這些渦結(jié)構(gòu)向后傳播起到了穩(wěn)定剪切層緩解其振蕩幅度的作用,這使得擾流片能有效地降低腔體內(nèi)部的壓強脈動水平。

        腔體流動的頻域分析方法可以給出壓強脈動的模態(tài)特征。使用聲壓級的指標對腔底不同位置進行分析。SPL=20lg,其中pref=2×1Pa是國際公認的人耳可聽聲壓閾值。圖5顯示了腔體底面從前到后部四個典型測壓點位置聲壓級變化的空腔試驗和計算結(jié)果以及帶有擾流片的DDES結(jié)果。

        圖4 空腔(上)和裝有擾流片腔體(下)渦量圖Fig.4 Vorticity magnitudes of clean cavity(up)and cavity with spoiler(down)

        圖5 x/L=0.05、0.35、0.65、0.95位置的聲壓級Fig.5 SPLat x/L=0.05,0.35,0.65,0.95

        通過聲譜圖可以看出腔體后部的噪聲等級顯著高于腔體前部,聲壓級譜中的第二階與第三階模態(tài)較為顯著的成為主導(dǎo)模態(tài)。計算結(jié)果對于第二階和第三階主導(dǎo)模態(tài)的預(yù)測與試驗結(jié)果吻合較好,高頻段的預(yù)測值出現(xiàn)多點峰值,這與Xiaoxian Chen等[25]的結(jié)果是一致的。第四階模態(tài)的具體位置并不是很清晰,但是幅值與試驗值差距較小。第一階模態(tài)在腔體前部位置的預(yù)測同樣也不是很清晰,而在腔體后部噪聲等級較大時被清晰的預(yù)測出來。通過在腔體前緣添加鋸齒形擾流片的控制手段后,計算結(jié)果表明腔體內(nèi)部的全局噪聲等級總體上有了明顯的下降,尤其是對于前兩階聲壓級模態(tài)的抑制效果顯著,對于第二階主導(dǎo)模態(tài)幅值的降低達到10dB以上。本文分析的鋸齒形擾流片裝置對聲壓級的第三階主導(dǎo)模態(tài)的抑制效果微弱,但是改變了該模態(tài)的發(fā)生頻率。

        通過Rossiter[15]關(guān)于腔體振蕩模態(tài)的半經(jīng)驗公式的預(yù)測f=,其中m和v分別表示模態(tài)數(shù)和自由來流速度,常數(shù)γ=0.25,κ=0.57分別表示相位延遲、平均擾動對流速度與自由來流的比值。分析腔體最后部的測壓點x/L=0.95位置處的各階模態(tài)頻率與幅值,如表1所示。

        由表1可知Rossiter公式可以大體預(yù)測出各階模態(tài)的頻率位置,與試驗值相比前兩階模態(tài)發(fā)生頻率高17Hz左右,后兩階模態(tài)頻率低于試驗值10~20 Hz,總體上認為誤差較小,說明Rossiter半經(jīng)驗公式對于空腔結(jié)構(gòu)聲調(diào)頻率的預(yù)測精度較高。DDES預(yù)測的各階頻率除了第一階模態(tài)預(yù)測提前外,其余誤差不超過4%。關(guān)于噪聲幅值的預(yù)測DDES在所有四階位置均十分準確,誤差均控制在2dB以內(nèi)。由于腔體流動呈現(xiàn)較強的非定常特性,為了捕捉剪切層的作用和各種尺度的渦的瞬時特性,需要數(shù)量巨大的網(wǎng)格,而且第一階與第四階模態(tài)的幅值相比于主導(dǎo)模態(tài)較小,使得非主導(dǎo)模態(tài)的計算準確度不及主導(dǎo)模態(tài)。計算安裝有擾流片以后通過DDES計算結(jié)果可以看出各階模態(tài)的頻率位置發(fā)生變化。與空腔的試驗值和DDES仿真值相比,所有四階模態(tài)頻率均向高值移動,前兩階模態(tài)增幅為20Hz左右,后兩階模態(tài)增幅為高達50Hz。原因是由于擾流片的作用,除了將腔口處的剪切層抬升遠離腔體外,它可以將大渦打碎成小渦,改變了腔體內(nèi)部的渦結(jié)構(gòu),從而改變了由渦碰撞引發(fā)的聲學(xué)共振現(xiàn)象的各階模態(tài),使得各階模態(tài)發(fā)生頻率均被增大。而與Rossiter公式預(yù)測的頻率相比,更大的差異性說明Rossiter公式在僅考慮腔體長度和來流速度的前提下,只適合預(yù)測空腔結(jié)構(gòu)而對于帶有控制措施的腔體預(yù)測結(jié)果略差。DDES結(jié)果表明擾流片結(jié)構(gòu)在對腔體內(nèi)部總體聲壓級水平進行抑制的同時,對前兩階模態(tài)的幅值也有較大幅度降低。處于主導(dǎo)地位的第二階模態(tài)在添加擾流片以后被降低了13dB,第一階模態(tài)被降低了8dB,降噪效果較為顯著。

        表1 x/L=0.95處各階模態(tài)結(jié)果對比Table 1 Mode comparisons at x/L=0.95

        4 結(jié) 論

        通過對來流Ma數(shù)為0.85下的長深寬比為5∶1∶1的方形空腔和帶有前緣擾流片的腔體使用DDES方法進行氣動噪聲分析后得到如下結(jié)論:

        (1)腔體內(nèi)部聲調(diào)噪聲的產(chǎn)生是由剪切層的不穩(wěn)定振蕩引發(fā),通過剪切層撞擊腔體后壁產(chǎn)生向上游傳播的壓強聲波進一步對剪切層干擾從而形成反饋回路。前緣鋸齒形擾流片可以有效將高速剪切層抬離腔體上方,并且其產(chǎn)生的小渦具有穩(wěn)定剪切層的作用,使得腔內(nèi)后壁的撞擊減弱達到抑制噪聲幅值的效果。

        (2)延遲分離渦模擬(DDES)方法對于腔體氣動噪聲分析的仿真結(jié)果較為準確。腔體所有模態(tài)的聲調(diào)噪聲等級預(yù)測結(jié)果誤差控制在2dB以內(nèi),對于第二三階主導(dǎo)模態(tài)的發(fā)生頻率誤差也控制在4%以內(nèi)。仿真認為前緣鋸齒形擾流片對第二階主導(dǎo)模態(tài)的降噪幅度在10dB以上,對腔體內(nèi)部整體聲壓級降幅在5dB左右。

        本文首次將DDES方法用于腔體氣動噪聲分析中,準確預(yù)測出腔體聲調(diào)噪聲的主導(dǎo)模態(tài)。但是對于占非主導(dǎo)地位的第一階模態(tài)發(fā)生頻率預(yù)測誤差較大,同時對于第四階模態(tài)的仿真結(jié)果出現(xiàn)多個峰值,難以準確確定發(fā)生位置。使用前緣鋸齒形擾流片的被動控制方法來抑制腔體氣動噪聲,可以引發(fā)對不同形式的擾流片(如平板或圓柱形)的降噪效果進行分析和探討。

        [1]Yang Dangguo,Li Jianqiang,Liang Jinmin.Sound generation induced by self-sustained oscillations inside cavities based on CFD and aeroacoustic theory[J].Acta Aerodynamica Sinica,2010,28(6):724-730.(in Chinese)楊黨國,李建強,梁錦敏.基于CFD和氣動聲學(xué)理論的空腔自激振蕩發(fā)聲機理[J].空氣動力學(xué)學(xué)報,2010,28(6):724-730.

        [2]Vakili A D,Wolfe R,Nagle P A.An experimental investigation of cavity aeroacoustics in high speed flows[R].Air Force Office of Scientific Research,Aerospace & Materials Sciences Directorate,1995.

        [3]Yang Dangguo,Luo Xinfu,Li Jianqiang,et al.Analysis of aeroacoustic characteristics in open cavities influenced by boundary-layer thickness[J].Acta Aerodynamica Sinica,2011,29(4):486-490.(in Chinese)楊黨國,羅新福,李建強,等.來流邊界層厚度對開式空腔氣動聲學(xué)特性的影響分析[J].空氣動力學(xué)學(xué)報,2011,29(4):486-490.

        [4]Lawson S J,Barakos G N.Assessment of passive flow control for transonic cavity flow using detached-eddy simulation[J].Journal of Aircraft,2009,46(3):1009-1029.doi:10.2514/1.39894

        [5]Morton M H,Hampson C D,Alexander R A.Final vibration and acoustic loads development for certification of the F-22advanced tactical fighter[C].Schaumburg,IL:49th AIAA/ASME/ASCE/ AHS/ASC Structures,Structural Dynamics,and Materials Conference,2008.

        [6]Rossiter J E.Wind-tunnel experiments on the flow over rectangular cavity at subsonic and transonic speeds[R].London:Ministry of Aviation Aeronautical Research Council,1966.

        [7]Lackey S,Tramel R W,Landrum D B.Weapons bay acoustic suppression using a novel rod in crossflow configuration[C].49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition,2011.

        [8]Thiemann C L,Milne G J,Vakili A D.An experimental investigation of supersonic cavity flow control with vertical cylinders[C].43rd AIAA Fluid Dynamics Conference,2013.

        [9]Nichols R H.A comparison of hybrid RANS LES turbulence models for a generic weapons bay with and without a spoiler[C].26th AIAA Applied Aerodynamics Conference,2008.

        [10]Zhuang N,Alvi F S,Alkislar M B,et al.Supersonic cavity flows and their control[J].AIAA Journal,2006,44(9):2118-2128.doi:10.2514/1.14879

        [11]Luo Baihua,Hu Zhangwei.Experimental study of flow induced cavity oscillationand its suppression by sound excitation[J].Journal of Nanjing University of Aeronautics and Astronautics,1996,28(3):331-336.(in Chinese)羅柏華,胡章偉.流動誘導(dǎo)空腔振蕩及其聲激勵抑制的實驗研究[J].南京航空航天大學(xué)學(xué)報,1996,28(3):331-336.

        [12]Si Haiqing,Wang Tongguang,Zong Huiying.Influence of the plate on the cavity flow-induced oscillationsand the modification of the oscillation frequency equation[J].Journal of Aerospace Power,2006,21(6):6-11.(in Chinese)司海青,王同光,宗慧英.腔內(nèi)平板對空腔自激勵振蕩的影響及預(yù)估振蕩頻率方程的改進[J].航空動力學(xué)報,2006,21(6):6-11.

        [13]Hao Zongrui,Wang Leqin,Zhou Zhonghai,et al.Numerical simulation of cavity flow field and aeroacoustic[J].Journal of Zhejiang University(Engineering Science),2013,47(1):131-138.(in Chinese)郝宗睿,王樂勤,周忠海,等.空腔流場及氣動噪聲數(shù)值模擬[J].浙江大學(xué)學(xué)報(工學(xué)版),2013,47(1):131-138.doi:10.3785/j.issn.1008-973X.2013.01.019

        [14]Ma Mingsheng,Zhang Peihong,Deng Youqi,et al.Numerical simulation investigation of supersonic cavity flow[J].Acta Aerodynamica Sinica,2008,26(3):388-393.(in Chinese)馬明生,張培紅,鄧有奇,等.超聲速空腔流動數(shù)值模擬研[J].空氣動力學(xué)學(xué)報,2008,26(3):388-393.

        [15]Zhang Lin.Aeroacoustic experimental investigation of the cavity flow fields in high speed wind tunnel[D].Changsha:National University of Defense Technology,2006.(in Chinese)張林.高速風(fēng)洞彈艙流場氣動聲學(xué)特性試驗研究[D].長沙:國防科技大學(xué),2006.

        [16]Yang Dangguo,Wu Jifei,Luo Xinfu.Investigation on suppression effect of zero-net-mass-flux jet on aerodynamic noise inside open cavities[J].Acta Aeronauticaet Astronautica Sinica,2011,32(6):1007-1014.(in Chinese)楊黨國,吳繼飛,羅新福.零質(zhì)量射流對開式空腔氣動噪聲抑制效果分析[J].航空學(xué)報,2011,32(6):1007-1014.

        [17]Yang Dangguo.Studies on aeroacoustic characteristics and noise suppressions for internal weapon bays[D].Mianyang:China Aerodynamics Research and Development Center,2010.(in Chinese)楊黨國.內(nèi)埋武器艙氣動聲學(xué)特性與噪聲抑制研究[D].綿陽:中國空氣動力研究與發(fā)展中心,2010.

        [18]Leonard L Shaw,Rodney Clark,Disk Talmadge.F-111generic weapons bay acoustic environment[J].Journal of Aircraft,1988,25(2):147-153.

        [19]Travin A,Shur M,Strelets M,Spalart P.Physical and numerical upgrades in the detached-eddysimulation of complex turbulent flows[M].Advances in LES of Complex Flows,Kluwer Academic Publishers:P Friedrich,W.Rodi,2002:239-254.

        [20]Smagorinsky J.General circulation experiments with the primitive equations[J].Monthly Weather Review,1963,91:99-165.

        [21]Spalart P R.Detached eddy simulation[J].Annual Review of Fluid Mechanics,2009,41:181-202.doi:10.1146/annurev.fluid.010908.165130

        [22]Rodriquez G,Velez C,Ilie M.Numerical studies of high-speed cavity flows using LES,DDES and IDDES[C].Grapevine,Texas:51th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition,2013.

        [23]Foster G W,Ross J A,Ashworth R M.Weapon bay aerodynamics wind tunnel trials and CFD modeling by QinetiQ UK[C]//Proceeding of the RTO/AVT Symposium on Flow-Induced Unsteady Loads and the Impact on Military Applications,Budapest,Hungary:2005.

        [24]Henshaw M J de C.M219cavity case-verification and validation data for computational unsteady aerodynamics[R].Rep.RTOTR-26,AC/323(AVT)TP/19,QinetiQ,UK:2000.

        [25]Chen Xiaoxian,Sandham N D,Zhang Xin.Cavity flow noise predictions[R].MSTARR DARP,Southampton:2007.

        猜你喜歡
        流片空腔腔體
        三擺臂式擾流片矢量氣動力數(shù)值仿真研究
        基于邊光滑有限元法的二維復(fù)合彈性空腔聲振特性分析
        一種多擾流片裝置的推力矢量特性數(shù)值研究 ①
        高鐵復(fù)雜腔體鑄造數(shù)值仿真及控制技術(shù)研究
        高鐵制動系統(tǒng)復(fù)雜腔體鑄造成形數(shù)值模擬
        擾流片式推力矢量控制的氣動力學(xué)研究
        擾流片式推力矢量噴管氣動特性數(shù)值模擬研究
        橡膠擠出裝置
        空腔參數(shù)對重力壩穩(wěn)定的影響分析
        前置污水去油池
        丰满爆乳在线播放| 日韩黄色大片免费网站| 日本本土精品午夜视频| 国内精品视频一区二区三区八戒| a级毛片100部免费看| 98在线视频噜噜噜国产| 日韩精品夜色二区91久久久| 精品国产亚洲第一区二区三区| 国产综合久久久久久鬼色| 中文字幕精品无码一区二区 | 国产真实强被迫伦姧女在线观看| 亚洲欧美日韩国产色另类| 日本五十路熟女在线视频| 国产免费二区三区视频| 国产网红主播无码精品| 亚洲自偷自拍另类图片小说| 中文乱码字幕在线中文乱码| 日韩精品极品免费视频观看| 国产福利精品一区二区| 亚洲男人精品| 国产一区二区三区经典| 蜜桃av在线免费网站| 欧美性猛交xxxx乱大交3| 91精品国产乱码久久久| 亚洲男人天堂一区二区| 毛多水多www偷窥小便| 综合激情网站| 久亚洲一线产区二线产区三线麻豆| 五月综合激情婷婷六月| 成人做爰69片免费看网站| 国产精品无码久久AⅤ人妖| 国产精品日韩av一区二区三区| 人妻无码一区二区不卡无码av| 欧美日韩在线观看免费| 国产一区二区av男人| 亚洲国产精品久久久av| 亚洲欧洲日本综合aⅴ在线| 成人永久福利在线观看不卡| 青青草手机视频免费在线播放| 精品无码国产自产拍在线观看| 久久久国产精品ⅤA麻豆|