999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

方腔流動氣動噪聲邊界散射數值預測方法研究

2015-03-28 11:06:45劉秋洪蔡晉生
空氣動力學學報 2015年5期

王 芳,劉秋洪,蔡晉生

(西北工業大學航空學院,陜西西安 710072)

方腔流動氣動噪聲邊界散射數值預測方法研究

王 芳,劉秋洪,蔡晉生*

(西北工業大學航空學院,陜西西安 710072)

方腔流動噪聲問題因其流場變化復雜且劇烈而倍受關注,本文主要開展高亞聲速方腔流動氣動噪聲數值預測方法研究。基于任意邊界條件的格林函數解和Lighthill聲模擬理論,提出可以考慮空間邊界影響的氣動噪聲積分計算方法。數值模擬包含流動和噪聲計算兩部分,通過二階精度的DDES模型進行流動數值模擬,邊界積分方法計算散射聲場分布。數值結果顯示聲場分布隨時間呈現周期性變化,與流場的脈動及其脈動周期一致。觀察點的聲壓級隨頻率逐漸下降且在諧波頻率突然增大。本文計算結果與高精度計算氣動聲學方法計算結果相符,表明該方法合理、可靠,并且具有較高的計算效率。

氣動噪聲;高亞聲速流動;方腔;格林函數

0 引言

隨著航空運輸量的日益增加,近年來對適航噪聲的標準也逐步提高,方腔繞流廣泛存在于飛機起落架艙、燃燒室、懸掛物儲藏室等構型中,并且方腔流動包含剪切層與方腔內部流動的相互作用、方腔內部流動的自激振蕩等多種復雜形式,也是流體與聲波相互作用的耦合流動,其流體動力學的不穩定對聲場分布產生劇烈影響,方腔噪聲的數值預測研究因而成為空氣動力學和氣動噪聲關注的熱點。

氣動噪聲的數值預測方法包括CAA(Computational Aeroacoustics)和 HCAA(Hybrid Computational Aeroacoustics)兩類,近年來國內外研究學者在方腔噪聲的數值預測方面開展了大量研究工作,Hardin[1]采用DNS(Direct Numerical Simulation)方法研究了二維低速方腔流動噪聲,Golerfelt等[2]結合高精度方法探討了具有層流和湍流邊界層的三維方腔流動噪聲。國內萬振華等[3]也采用DNS方法考察了亞聲速方腔流動流致振蕩產生的氣動噪聲。HCAA[4]克服了CAA方法對于網格和數值格式的高要求及限制,且時間效率較高,被廣泛的用于復雜流動氣動噪聲數值預測,并且逐漸發展為一類適用于實際工程問題的方法,有利推動了航空飛行器氣動噪聲的研究工作。

考慮到HCAA方法在氣動噪聲研究領域的極大優勢,研究學者們提出了諸多性能優良的聲學計算模型,主要包括 Curle理論[5]、FW-H(Ffowcs Williams-Hawkings)方程[6]、渦聲理論(Vortex Sound Theory)[7]和Kirchhoff方程[8]等。渦聲理論闡述了旋渦是聲音的本質,然而不能反映聲場的多極子特征,并且該方法必須在流場渦量被精確計算的前提下才可以實施。Kirchhoff方程僅能在聲波傳播的線性區域內執行,對積分邊界的選取提出了嚴格的要求。Curle理論和FW-H方程是Lighthill聲模擬理論的拓展形式,其顯著優勢在于可以考慮固體邊界對聲場分布的影響,同時從物理角度清晰的解釋聲波傳播的物理本質,因此也成為當前最有效且簡單的氣動噪聲數值計算方法。Golerfelt[9]結合高精度方法和FW-H積分方程研究了二維亞聲速開口方腔流動噪聲,Ask[10]結合Curle理論考察了低馬赫數開口方腔氣動噪聲聲場分布的偶極子特征,Zhang等[11]采用LES模型和FW-H開展了閉合方腔氣動噪聲的數值預測。然而,上述研究工作存在流場高精度計算的限制,不易開展實際工程應用。作者提出了適用于自由空間內非緊致邊界氣動噪聲的數值預測方法[12],該方法克服了高精度流場計算帶來的巨大工作量,不足的是需要計算物理空間的散射聲源。

本文基于任意邊界條件下點聲源的格林函數解和Lighthill聲模擬理論,提出可以考慮物理空間邊界散射影響的氣動噪聲數值預測方法,考察了方腔內點應力源的噪聲并分析方腔邊界的散射影響;同時,深入開展高亞聲速方腔流動氣動噪聲的數值預測,探討流場與噪聲之間的關系。

1 氣動噪聲數值預測方法

任意邊界條件下,單位點聲源的傳播遵循如下波動方程:

此處g(x,t;y,τ)表示格林函數,又稱為聲波基本解。c表示聲波的傳播速度;x、y分別表示觀察點和聲源點位置。方程(1)描述了τ時刻由y發出的聲波經過t-τ時間段到達x的傳播過程。圖1給出了流場脈動和聲波傳播示意圖,其中y1、y2和y3為聲源點,S表示任意結構物體表面,ΓF和ΓA分別表示流場邊界和聲學遠場邊界,ΩF和ΩA分別表示流場計算區域和聲場計算區域。對于包含空間邊界的物理模型,Sb表示其對應的空間邊界。

圖1 氣動噪聲傳播示意圖Fig.1 The diagram of aeroacoustic noise propagation

采用HCAA開展氣動噪聲的數值計算時廣泛使用的Lighthill聲模擬理論[13]為

對于高雷諾數流動,可以忽略流體黏性應力帶來的影響[14],Lighthill應力張量表示為

式中:p'=p-p0,ρ'=ρ-ρ0分別為壓力脈動和密度脈動;p0和ρ0分別表示來流壓力和密度;ui為流體沿i方向的運動速度。結合方程(1)、(2)可得

方程(3)可以描述物理空間邊界對氣動噪聲的影響,詳細推導過程可參閱文獻[12],其中ph和ρh分別表示可壓縮流體脈動引起的壓力和密度脈動,ρ0表示均勻來流密度。假設規定邊界法向方向由物面指向流場區域,vn表示邊界沿法向的運動速度,此外Tij=ρuiuj+(ph-c2ρh)δij。對于包含非緊致邊界的物理模型成立珔S=S∪Sb。采用方程(3)開展氣動噪聲的數值計算,首先求解物體表面上的散射聲壓,然而將觀察點布置在物體表面上會出現聲源點和觀察點距離較近或者相互重合的情形,從而導致數值積分的奇異性。借鑒Khalighi的研究工作,消除積分奇異性可得:

其中z表示物體表面上的觀察點,滿足‖z-y‖2→0,{z}表示邊界珔S去掉z點的剩余部分。對應給出遠場噪聲計算的表達式:

本文重點考察方腔流動氣動噪聲,在邊界靜止的情形下,可以給出方程(4)和(5)的頻域形式:

其中G(x,y,ω)為g(x,t;y,τ)的頻域形式;ω=2πf,表示圓周頻率;f表示頻率,假設聲場計算的時間步長和采樣數目分別為Δt和N,采樣定律[14]指出聲學計算可分辨的最高截止頻率ftop=1/(2NΔt)。以不包含非緊致邊界的半空間問題為例,格林函數滿足如下邊界條件:

此處珔S=Sb,半空間格林函數的表達式可參閱文獻[15]。

2 方腔邊界散射效應

選取如圖2所示的方腔模型,方腔深度2.54mm,方腔長度滿足L/D=3,沿流向設置流場計算域,方腔前端長度和后端長度分別為5D和4.8D,方腔前端頂點為坐標原點。假設x、y、z坐標軸分別沿來流方向,腔體深度方向和腔體展向,流場計算區域為x∈[-5D,7.8D],y∈[-D,6.8D],z∈[0,4/3D]。

圖2 二維方腔流場計算區域Fig.2 The computational field of two-dimensional cavity flow

為了考察方腔邊界引起的散射效應,結合半空間格林函數考察單個應力源的噪聲,為方便起見取T11=1,T22=T12=T21=0。以方腔展向中間截面為參考面,在方腔內放置兩個點源N1和N2,分別位于方腔內部和方腔左端,位置坐標為(2/3L,-1/2D)和(-1/2D,2/3L),給定k0=35,其對應波長λ≈0.18。沿直線l:x+1.4y=0均勻分布100個觀察點,圖3顯示了基于自由空間和半空間格林函數獲得的噪聲分別用G0和Gb表示,考察不同波數下噪聲沿直線的變化規律,其橫縱坐標分別對應于和噪聲幅值,此處縱坐標采用進行量級處理。圖3顯示半空間格林函數和自由空間的噪聲存在顯著差異。與N2相比,由于N1位于方腔內部,邊界的散射影響較為強烈,半空間格林函數獲得的噪聲并不隨著波數的增大呈現相同的變化趨勢,這與文獻[15]獲得的結果一致,邊界的散射效應使得噪聲變化劇烈且復雜。由于N2與方腔邊界距離較遠,邊界的散射影響較小,因而半空間格林函數獲得的噪聲與自由空間幾乎保持相同的變化趨勢。唯一不同的是其聲壓幅值沿直線先增大后減小,這是因為隨著r的增大,聲源N2與觀察點的距離先增大后減小,也符合聲波的傳播規律。半空間格林函數獲得噪聲與自由空間格林函數噪聲的數值差異說明方腔邊界的散射效應不容忽視,并且應力源受到的散射影響與其物理位置密切相關。

圖3 不同波數下自由空間和半空間格林函數獲得的噪聲Fig.3 The noise of free and half space Green's function with different wavenumber

3 方腔流動噪聲數值計算

3.1 方腔流動計算

已知自由來流馬赫數為0.8,來流溫度為T= 298.15 K,聲波傳播速度為c=346.14 m/s,來流速度為276.91 m/s,來流密度為ρ=1.269 6 kg/m3,相對于方腔深度的雷諾數為48 600。非定常流場數值計算的步長為Δt=10-8s,采用基于SA模型的二階精度DDES方法[16]進行計算。采用Block網格劃分技術將計算區域劃分為168塊,方腔內部包含48個網格塊,且沿流向和深度方向分布的網格密度為257× 193,沿深度方向網格距離壁面最小距離為10-5m。展向均勻分布33層網格,腔內網格數目約164萬,腔內與腔外整個流場計算域的三維網格單元總數近280萬。方腔底部沿流向分布選取展向分布的中間截面z=2/3D為參考面,圖4給出流動計算獲得的不同時刻的渦量云圖。起始時刻腔體內部出現流動旋渦,隨著時間的延續旋渦產生的新脫落渦不斷向腔體后端發展,與腔體后端頂點發生碰撞,導致腔體底部及前端產生分離渦,流場結構發生劇烈變化。方腔內部的流線圖分布見圖5,方腔內部流線封閉,表明方腔內部逆向流動和剪切層流動呈現交替變化,在方腔中間x=2L/3出現大渦,方腔底部前端和后端分別出現零碎小渦。

圖4 不同時刻方腔內部渦量云圖Fig.4 The vortex contours inside the cavity at different time

定義速度均方根值:

圖5 方腔內部流線圖Fig.5 The streamlines in the cavity

其中u'=u-U0,v'=v-U0,分別表示x、y方向的速度脈動,U0為來流速度。圖6給出方腔x=2L/ 3的速度型變化,與Gloerfelt等[2]采用層流模型和Smagorinsky模型獲得的結果吻合,其中紅色△和藍色○點分別表示層流和 Smagorinsky模型的結果,藍色曲線為DDES模擬的結果。與來流速度相比,x和 y方向的速度脈動幅值最高達到30%。流體脈動速度沿豎直方向先增大后減小,并且在y =0附近區域最為劇烈。與方腔內部流動相比,方腔上部流體的脈動速度衰減的更快,這是因為方腔上方形成剪切層,而開式方腔內部受上方剪切層擾動較弱,并且始終存在大渦結構、分離渦及零碎小渦。

圖6 方腔內部x=2L/3速度型分布Fig.6 The distribution of velocity profile inside the cavity with x=2L/3

3.2 方腔噪聲計算

考慮到方腔流動的特殊性,假設將方腔底部邊界看作半空間邊界的一部分,可以按照半空間模型研究方腔噪聲。為減少三維模型由于聲源樣本數目過大所需的巨大存儲空間,同時保證聲學計算的分辨率,選取展向中間截面z=2/3D的流場為參考面。聲學計算每隔5個時間步輸出一個樣本,對應的時間步長為Δt'=5Δt,首先在x=(-1.04D,4.64D)設置檢測點,圖7給出了流場壓力脈動隨時間的變化及該點的SPL隨頻率(St=fL/U)的變化,表明方腔內部流動呈現不規則的周期性變化,流場脈動周期為T=70Δt',并且聲壓的能量譜在St數為0.78時達到最大。其次,沿時間方向采集1 050個樣本開展聲場計算,結合時域方法詳細研究方腔流動噪聲的瞬時變化并探討流動與噪聲之間的關系,在計算區域布置4 073個網格點,圖8給出了聲壓云圖的瞬時分布,以半個周期為間隔輸出兩個周期的聲場計算結果。由圖可知,聲場分布隨時間呈現周期性變化并與流場的脈動及脈動周期一致。由于方腔內部的前端和后端區域流體脈動劇烈,這兩個區域的聲場脈動也較為強烈;運動旋渦在腔體后端發生分離導致流場發生劇烈變化,聲場計算結果表明聲波由腔體后端不斷向腔體斜上方傳播,并且聲壓幅值隨著流場脈動發生周期性變化;流動狀態達到穩定之后流場進入下一個脈動周期,聲場也呈現相同的周期性變化。圖7顯示流場的周期大概為 T=70Δt',圖8的聲壓云圖顯示聲場變化周期與流場相同,充分證明流場和聲場的一致性,流場脈動較為劇烈的區域噪聲也相對強烈,這說明聲學計算結果取決于流動計算,同時也驗證了當前計算方法的有效性。為了詳細研究流場內部局部位置聲壓的變化及分布,選取如圖9所示的曲線,考察不同空間位置的聲壓分布。紅色和綠色實線分別表示直線l1和l2,它們的x、y坐標滿足x+1.4y=0和0.85x+y=2.6D,且分別以O和O'為起點,P為終點。在直線l1和l2上均勻布置200個觀察點,圖10給出了聲壓幅值分布曲線,此處c1、c2分別對應于直線線l1和l2,其中橫坐標分別對應于,表示觀察點x與O和O'的距離,例如St=0.78的紅色實線表示頻率為0.78時聲壓沿直線的變化。觀察圖10可知,聲壓幅值隨著頻率的增大顯著降低,然而相同頻率下聲場隨觀察點位置不同而變化,圖10(a)顯示聲壓幅值沿直線l1逐漸降低,這與Gloerfelt采用DNS方法[9]獲得的結論一致,圖10(b)給出的聲壓幅值在r=3.95D迅速下降,此后聲壓幅值沿直線l2表現出與直線l1相同的變化趨勢。對位于方腔正上方的觀察點滿足r<3.95D,其噪聲幅值變化復雜,渦和二次渦的復雜流場特征密切相關。由此可見,這與方腔內部流動聲源脈動劇烈且不斷產生分離小圖8和圖10清晰的揭示了高亞聲速方腔流動氣動噪聲的聲場特征及其原因。

圖7 流體壓力脈動Fig.7 Pressure fluctuation induced by the flow

圖8 方腔聲壓云圖變化Fig.8 The variations of pressure contours in cavity

圖9 聲場觀察點分布Fig.9 The distributions of observers in acoustic field

圖10 聲壓幅值隨觀察點的變化Fig.10 The variations of pressure amplitude with different observers

已知圖7給出觀察點x=(-1.04D,4.64D)的壓力脈動及其PSD隨頻率的變化,聲學計算仍然選取上述檢測點為觀察點,圖11給出了聲壓級隨St數發生的變化,并與Gloerfelt采用高精度方法[2]獲得的結果進行了數值比較,其中藍色點線表示Gloerfelt的結果,粉色實線表示當前計算結果。圖11顯示本文方法獲得的聲壓級在St=0.78達到峰值,與圖6(b)給出的流場壓力脈動PSD曲線的結論一致,已知Gloerfelt采用高精度方法在St=0.75達到相同的峰值,兩者的數值誤差僅為4%,并且均與Rossiter公式[17]的預測結果符合。隨著頻率的增大聲壓級總體上呈下降趨勢,且聲壓級僅在諧波頻率處出現峰值,本文的數值結果與Gloerfelt[2]及Karamcheti[18]獲得的結論一致。在高頻率下本文方法與Gloerfelt方法的結果顯現出數值差異,Lighthill[19]指出低頻噪聲的能量占總聲能的50%以上,可見頻率St≤3的噪聲攜帶了大部分聲能,并且本文結果是基于二維聲場計算修正得到的[20],并不是真實三維流場獲得的遠場噪聲,因而這種數值差異是可以接受的。

圖11 觀察點聲壓級隨頻率的變化Fig.11 The variations of SPL with frequency for observer

4 結論

本文提出了一種可以考慮物理空間邊界影響的氣動噪聲積分計算模型,主要研究高亞聲速方腔流動氣動噪聲并探討流動與噪聲之間的關系,主要結論如下:

(1)采用半空間格林函數可以考慮半空間邊界對聲場的散射影響且無需計算散射聲源,大大節省了計算資源,噪聲計算簡單且高效。

(2)觀察點的聲壓級在諧波頻率處出現明顯的峰值,當前方法獲得的結果與高精度計算方法吻合,說明這種方法可以精確求解方腔流動邊界對氣動噪聲的散射效應,同時可以避免流場高精度計算所需要的巨大工作量。

(3)聲場分布表現出與流場相同的脈動周期和變化特征,流場脈動劇烈的區域其噪聲也相對強烈,觀察點的流場脈動和噪聲脈動表現出相同的變化趨勢,反映了流場決定聲場的物理事實。

[1]Hardin J C.Sound generation by flow over a two-dimensional cavity[J].AIAA Journal,1995,33(3):407-512.

[2]Gloerfelt X,Bogey C,Bailly C,et al.Aerodynamic noise induced by laminar and turbulent boundary layers over rectangular cavities[C]//8th AIAA/CEAS AeroacousticsConference,Colorado: 2002,2002-2476.

[3]Wan Zhenhua,Zhou Lin,Sun Dejun.Numerical investigation of flow-induced oscillations and noise from a rectangular cavity[J].Acta Aerodynamica Sinica,2012,30(3):291-298.(in Chinese)萬振華,周林,孫德軍.方腔流致振蕩及噪聲的數值研究[J].空氣動力學學報,2012,30(3):291-298.

[4]Bailly C,Bogey C,Gloerfelt X.Some useful hybrid approaches for predicting aerodynamic noise[J].Competes Rendus Mechnique,2005,333:666-675.

[5]Curlen N.The influence of solid boundaries on aerodynamic sound[J].Proceedings of the Royal Society of London,Series A,1952,213(1187):505-514.

[6]Ffowcs Williams J E,Hawkings D L.Sound generation by turbulence and surface in arbitrary motion[J].Proceedings of the Royal Society of London,Series A,Mathematics and Physical Sciences,1969,264:321-342.

[7]Howe M S.Theory of vortex sound[M].New York:Cambridge U-niversity Press,2003.

[8]Farassat F,Mayers,M K.Extension of Kirchhoff's formula to radiation from moving surfaces[J].Journal of sound and vibration,1988,123:451-461.

[9]Gloerfelt X,Bailly C,Juve D.Direct computation of the noise radiated by a subsonic cavity and application of integral methods[J].Journal of Sound and Vibration,2003,266:119-146.

[10]Ask J,Davidson L.Sound generation and radiation of an open twodimensional cavity[J].AIAA Journal,47(6):1337-1349.

[11]Zhang N,Shen H C,Yao H Z.Numerical simulation of cavity flow induced noise by LES and FW-H acoustic analogy[J].Shanghai: 9th International Conference on Hydrodynamics,2010,242-247.

[12]Wang Fang,Liu Qiuhong,Cai Jinsheng.An unified aeroacoustic computational integral method of noise radiation and scattering with noncompact bodies[J].Acta Aeronautica et Astronautica Sinica,2013,34(11):2482-2491.(in Chinese)王芳,劉秋洪,蔡晉生.非緊致結構氣動噪聲輻射/散射統一積分計算方法[J].航空學報,2013,34(11):2482-2491.

[13]Lighthill M J.On sound generated aerodynamically.I.General theory[J].Proceedings of the Royal Society of London,Series A: Mathematical and Physical Sciences,1952,211(1107):564-587.

[14]Morfey C L.The role of viscosity in aerodynamic sound generation[J].Acoustics,2003,2(3):225-240.

[15]Wang F,Cai J S,Liu Q H.Aerodynamic noise calculations of ground effect based on tailored Green's function[J].Journal of Aircraft,2015,52(1):21-30.

[16]Cai J S,Pan S C,Li W F,Zhang Z K.Numerical and experimental investigations of a nonslender delta wing with leading-edge vortex flap[J].Computers and Fluids,2014,99:1-17.

[17]Rossiter J E.Wind-tunnel experiments on the flow over rectangular cavities at subsonic and transonic speeds[R].London:Her Majesty's Stationery Office,1964.

[18]Karamcheti M.Acoustic radiation from two-dimensional rectangular cutouts in aerodynamic surfaces[R].Washington:National Aeronautics and Space Administration,1955.

[19]Lighthill M J.On sound generated aerodynamically.II.Turbulence as a source of sound[J].Proceedings of the Royal Society of London,Series A:Mathematical and Physical Sciences,1954,222 (1148):1-32.

[20]Orselli R M,Meneghini J R,Saltara F.Two and three-dimensional simulation of sound generated by flow around a circular cylinder[C]//15th AIAA/CEAS Aeroacoustics Conference,Flordia:2009,2009-3270.

Numerical prediction method of aerodynamic noise including boundary scattering generated by high subsonic cavity flow

Wang Fang,Liu Qiuhong,Cai Jinsheng*
(1.School of Aeronautics,Northwestern Polytechnical University,Xi'an 710072,China)

The aerodynamic noise of cavity flow is a concern problem due to the complex and serious dynamic fluctuation,the goal of this paper is to study the numerical prediction method of aerodynamic noise generated by high subsonic cavity flow.With the help of Green's function satisfied arbitrary boundary conditions and Lighthill's Acoustic Analogy,the aeroacoustic integral computational method taking into account the influence of physical boundary is presented.The numerical calculation includes two parts,flow and noise computations.The two-order DDES model is used to perform the fluid simulation while the present integral computational method calculates the distribution of scattering noise.Numerical results show that the distribution of noise field varies periodically,and it is in agreement with fluid fluctuation.The SPL gradually reduces with frequency and appears higher amplitudes at harmonic frequencies.Numerical results obtained with the present method agree well with that of high-order aeroacoustic computational method,which indicate that this method is valid and efficient to calculate high subsonic cavity noise.

aeroacoustic noise;high subsonic flow;cavity;Green's function

V211.3

:Adoi:10.7638/kqdlxxb-2015.0004

0258-1825(2015)05-0617-07

2015-01-06;

:2015-03-03

國家自然科學基金(11002116);西北工業大學基礎研究基金(GCKY1006)

王芳(1984-),女,寧夏人,博士研究生,研究方向:氣動噪聲算法研究及數值預測.E-mail:fangw1211@163.com

蔡晉生*,西北工業大學航空學院.Tel:13659115379,E-mail:caijsh@nwpu.edu.cn

王芳,劉秋洪,蔡晉生.方腔流動氣動噪聲邊界散射數值預測方法研究[J].空氣動力學學報,2015,33(5):617-623.

10.7638/kqdlxxb-2015.0004 Wang F,Liu Q H,Cai J S.Numerical prediction method of aerodynamic noise including boundary scattering generated by high subsonic cavity flow[J].Acta Aerodynamica Sinica,2015,33(4):617-623.

主站蜘蛛池模板: 国产精品无码翘臀在线看纯欲| 国产一级片网址| 亚洲AV成人一区国产精品| 日韩黄色精品| 午夜啪啪网| 自慰网址在线观看| 成人在线不卡视频| 无码精油按摩潮喷在线播放 | 久久人与动人物A级毛片| 2018日日摸夜夜添狠狠躁| 99精品国产高清一区二区| 国产三级视频网站| 欧美日韩午夜视频在线观看| 久久国产精品电影| 四虎影视8848永久精品| 日本91视频| 在线观看国产网址你懂的| 国产成人AV男人的天堂| 亚洲欧美色中文字幕| 男女性午夜福利网站| 欧美h在线观看| 久久99蜜桃精品久久久久小说| 青青草一区二区免费精品| 国产精品3p视频| 免费在线播放毛片| 精品福利一区二区免费视频| 免费可以看的无遮挡av无码| 午夜国产理论| 亚洲欧美国产五月天综合| 91综合色区亚洲熟妇p| av午夜福利一片免费看| 久久久久亚洲AV成人人电影软件 | a毛片在线免费观看| 国产永久免费视频m3u8| 亚洲国产日韩一区| 国产极品美女在线播放 | 国产新AV天堂| 国产视频大全| 成人小视频网| 精品成人一区二区| 19国产精品麻豆免费观看| 国产精品天干天干在线观看| 黑色丝袜高跟国产在线91| 国产91视频观看| 精品免费在线视频| 国产成熟女人性满足视频| 亚洲无码精品在线播放| 国产精品网拍在线| 午夜视频免费试看| 影音先锋亚洲无码| 午夜国产大片免费观看| 国产手机在线小视频免费观看| 无码在线激情片| 亚洲欧州色色免费AV| 777午夜精品电影免费看| 国产成人高清精品免费5388| 国产精品妖精视频| 国产成人精品视频一区二区电影| 国产精品美女自慰喷水| 亚洲国产中文欧美在线人成大黄瓜| 毛片免费在线视频| 欧美成人精品欧美一级乱黄| 精品偷拍一区二区| 一级做a爰片久久免费| 国产精品福利社| 欧美成人影院亚洲综合图| 青青草91视频| 免费观看亚洲人成网站| 刘亦菲一区二区在线观看| 国产网友愉拍精品视频| 亚洲精品成人福利在线电影| 伊人蕉久影院| 国产一级二级在线观看| 国产打屁股免费区网站| 国产91久久久久久| 国产福利免费观看| 青青青伊人色综合久久| 精品国产成人三级在线观看| 极品av一区二区| 日韩在线欧美在线| 99精品热视频这里只有精品7| 99这里只有精品免费视频|