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

有限長圓柱繞流氣動噪聲源特性分析

2019-04-11 08:32:28楊志剛劉洋王毅剛
聲學技術 2019年1期
關鍵詞:區域

楊志剛,劉洋,王毅剛

?

有限長圓柱繞流氣動噪聲源特性分析

楊志剛1,2,劉洋1,王毅剛1

(1. 同濟大學上海地面交通工具中心,上海 201804;2. 北京民用飛機技術研究中心,北京 102211)

氣動噪聲源的能級、分布特性及其產生根源還不夠清晰。以有限長三維圓柱繞流為研究對象,基于聲源方程分析氣動噪聲源的種類構成及其與氣動參數的關系,通過數值計算得到可壓非定常流場,利用氣動參數定量計算圓柱頂部、中部和底部的聲源大小分布,研究聲源的分布特性和產生根源。結果表明,在有限長圓柱繞流場中,以偶極子聲源為主,單極子聲源可以忽略不計,四極子源項的值比偶極子小1~2個數量級。偶極子主要分布在來流分離點及圓柱后壁面湍流渦二次碰撞區域,四極子主要分布在來流分離點及其向后拖曳區域。偶極子聲源主要由于圓柱兩側渦脫落處的脈動壓力在橫向(y方向)上的二階梯度引起。以上結果為氣動噪聲控制的進一步研究提供了借鑒和參考。

氣動噪聲;有限長圓柱;聲源方程;偶極子源

0 引言

氣動聲學是指研究流體自身及流體與固體邊界相互作用發聲機理的一門學科[1],氣動噪聲廣泛存在于航空、高速列車及汽車等領域,研究氣動噪聲具有重要的實際意義,探討、分析氣動噪聲源特性,是控制氣動噪聲至關重要的關鍵之一。Lighthill[2-3]提出了聲類比理論,并沒有考慮流動中可能存在的固壁影響,只適用于自由空間,所以Lighthill方程右側項只有四極子。Curle[4]運用基爾霍夫積分公式改進了Lighthill聲類比理論,考慮了靜止固體壁面對流動發聲的影響,固壁的影響相當于在固壁表面分布有一層偶極子源,所以方程右側的源項為偶極子和四極子的疊加。Ffowcs Williams等[5]應用廣義函數法將Curle的結果擴展到考慮運動固體邊界對聲音的影響,即物體在流體中運動的發聲問題,提出Ffowcs Williams-Hawkings方程(FW-H方程)。從FW-H方程可以看出,運動物體與流體相互作用產生的聲場是由四極子源、偶極子源以及由于質量移動效應所產生的單極子源的疊加組成的。Howe[6-7]簡化了Lighthill方程,得到了適用于低馬赫數、高雷諾數的流場的Howe方程,方程的右端項表明流場中的漩渦可以作為聲源。

關于圓柱繞流氣動噪聲的數值計算方法和噪聲特性分析的相關研究已經比較廣泛,而對于圓柱繞流氣動噪聲源的研究比較缺乏,尤其是三維氣動噪聲源的研究。Cox等[8-9]通過雷諾時均法(Reynolds Average Navier-Stokes, RANS)計算二維非定常圓柱繞流,并基于FW-H方程計算了偶極子源和四極子源,發現四極子源對總噪聲的貢獻相對于偶極子源可以忽略。Takaishi等[10]通過大渦模擬得到瞬態流場特征,運用廣義格林公式準確預測了偶極子聲源在高速列車轉向架區域的分布。A. Iida等[11]通過廣義格林公式和相干輸出功率(Coherent Output Power, COP)兩種方法預測了偶極子聲源的分布,兩種方法的結果吻合得很好,結果顯示偶極子聲源是在卡門渦街形成區域附近產生的,分離的剪切層對偶極子聲源有著重要的作用。鄭朝榮等[12]采用可實現-湍流模型和寬頻帶噪聲源模型相結合的方法研究鈍體繞流的氣動噪聲源,結果表明鈍體繞流氣動噪聲主要分布在來流分離和湍流運動比較劇烈的位置,四極子源項的值相對于偶極子源項小很多。

從上述的研究文獻中發現,對于三種聲源的具體數量級大小、分布還需要深入研究,聲源的產生根源也不是很清楚。本文以有限長圓柱繞流為例,通過數值計算定量研究單極子、偶極子和四極子聲源的具體大小和分布特性,并分析了聲源的產生根源。

1 氣動噪聲理論

1.1 基爾霍夫定理

當反射、吸收甚至發射聲波的固體邊界存在時,我們假設存在一個包圍源和固體邊界的假想面,觀測點在表面的外部,如圖1所示,推導得

該結果稱為基爾霍夫定理[13]。Curle運用這一定理將Lighthill定理推廣到靜止固壁存在情況下聲場的求解,并應用到后續的FW-H方程中。此后在流體噪聲和氣動噪聲的各種推廣應用中,有關聲場的解均基于基爾霍夫積分公式。

1.2 聲源方程及聲源項的確定

由靜止介質中的連續方程和動量方程,可得聲波方程:

聲波方程描述了聲波在靜止介質中的傳播,但是不包含任何有關聲源的信息?,F在假定在域中有質量源,單位為kg·m-3,則連續性方程為

將式(4)對微分,將式(5)對微分,然后兩式相減得

該式為聲源的通用方程[14]。方程的左邊為傳播項,方程的右邊三項依次為單極子聲源、偶極子聲源和四極子聲源。

式(7)的物理意義為力點源可能由于壓力梯度產生,也可能由于速度的旋度產生,或者是兩者的組合。故偶極子聲源為

四極子聲源是由于物體周圍紊流流體內的應力變化產生的,本質上是一種流體湍流應力源。-方程(動量方程)為

將式(10)代入式(9),得:

式(11)與連續性方程聯立,解得

故四極子源項為

2 圓柱繞流數值計算

2.1 數值計算方法

以氣動-聲學風洞試驗的圓柱繞流模型為研究對象,如圖2所示,圓柱直徑=0.1 m,高度和直徑之比為18.5:1,高=1.85 m,在圓柱頂部有和試驗模型一致的半球,使圓柱頂部流動分離不過于劇烈,圓柱底部有一個錐度很小的底座,使圓柱固定更加穩定。流場計算區域示意圖如圖1所示,以圓柱底部中心點為坐標原點,流向長度取為-10~25,橫向為-8~8,展向為0~28。聲源區為如圖所示的矩形區域,尺寸為8.2×4×20,聲源區前端距圓柱軸線1.2,后端距圓柱軸線7,基本包含了圓柱繞流渦脫落影響的核心區域。

計算域分塊劃分六面體網格,用O-Grid型網格劃分圓柱的邊界層和加密區,網格在圓柱近壁面沿周向等分,徑向在近壁面密布邊界層,大渦模擬通常要求+(第一層網格質心到壁面的無量綱距離)值小于1,+也可以增加到5[15],由于計算資源的限制,本次計算+取2,第一層邊界層厚度為0.018 mm,增長率為1.05,邊界層一共39層。聲源區最大網格尺寸為10 mm,最終網格總數約為4 086萬。

圖2 圓柱繞流模型圖

表1 定常流動計算的邊界條件設置

2.2 氣動噪聲數值計算及風洞試驗

氣動噪聲數值計算也是在Fluent的非定常流場中進行,采用Ffowcs Williams-Hawkings(FW-H)聲學模型,并設置選取聲源區,將大渦模擬非定常流動計算中獲得的流動參數,如壓力和速度等,代入FW-H方程的聲源項中,通過對遠場的積分,得到遠場接收點的聲壓數據。

為了驗證上述圓柱繞流仿真計算的正確性及氣動噪聲控制措施,本文開展了風洞實驗研究。試驗在同濟大學上海地面交通風洞中心的整車氣動-聲學風洞中進行。試驗用圓柱為實心鋼圓柱,頂部設計成一個半球結構以避免強的氣動噪聲產生,其根部采用圓臺結構(增大其與地面的接觸面積)固定于試驗段的地面上,以保證圓柱在較高風速下不晃動。圓柱直徑為0.1 m,高為1.85 m。圖3標明了試驗圓柱在風洞中的具體位置和遠場傳聲器測點位置。3個傳聲器軸線距圓柱5 m,傳聲器間距1.3 m,高為1.2 m。傳聲器為丹麥GRAS公司的產品,數據采集及分析系統為HEAD公司的產品。試驗內容包括光滑圓柱氣動噪聲試驗和覆蓋多孔介質圓柱氣動噪聲試驗。試驗風速為120 km·h-1。

圖3 圓柱和麥克風位置俯視圖

圖4為風速為120 km·h-1下,距圓柱5 m遠處傳聲器測量和仿真計算的A計權頻譜。從圖4中的測量結果可以看出,在64 Hz附近出現了明顯峰值,聲壓級為78.18 dB,是圓柱典型的渦脫落引起的。由于圓柱頂部和地面影響,在峰值附近出現較低的另外峰值并有一定的帶寬,但并沒有掩蓋繞流圓柱渦脫落的峰值特征。從仿真結果可以看出,在67 Hz附近出現明顯峰值,聲壓級為72.74 dB。本圓柱渦脫落頻率為66.67 Hz,測量結果和仿真結果都接近該值,峰值聲壓級相差5.44 dB,表明仿真結果和試驗結果兩者峰值頻率和大小基本一致,說明了仿真計算的正確性。

圖4 仿真和實測的遠場噪聲頻譜圖

3 圓柱繞流數值計算結果分析

3.1 氣動噪聲源分布特性

氣動噪聲源的分布特性與流場各項氣動參數有關,對于有限長度的圓柱,圓柱頂部、中部和根部的流態不盡相同,需要分別加以分析。

對于圓柱頂部,選取=0截面進行分析。圖5為=0截面圓柱頂部的單極子分布圖,單極子源項的數量級為0~10,數量級很小,可以忽略不計。

圖5 圓柱頂部單極子分布

圖6為=0截面圓柱頂部的偶極子分布圖。在圓柱頂部-90°~27°,偶極子緊貼圓柱頂部分布,且分布區域逐漸變厚;在27°位置處,偶極子的分布區域開始從圓柱表面分離,形成一個拖曳區域;在90°位置處也存在一個小的分離,同時形成一個小的拖曳區域;在圓柱截面來流后壁面1.63~1.76 m處,存在一個貼著壁面分布的偶極子區域,厚度約為10 mm。圓柱頂部偶極子的最大值數量級為109~1010。

圖7為=0截面圓柱頂部的四極子分布圖,在圓柱頂部-90°~0°,緊貼圓柱表面分布著一層很薄的四極子;從0°~43°四極子區域逐漸變厚,并在43°處開始從圓柱表面分離,產生一個較大的拖曳區域;在圓柱截面來流后壁面1.62~1.80 m處,存在一個分布較為廣泛的四極子區域,在1.68~1.77 m處,四極子的分布范圍更加向后拉長,長度約為85 mm。圓柱頂部四極子的最大值數量級為108~109。

圖7 圓柱頂部四極子分布

對于圓柱中部,分別選取=0.8、1.0、1.2 m處的橫截面加以分析。圖8顯示了圓柱在=0.8、1.0、1.2 m處橫截面上單極子聲源在圓柱繞流場中的大小分布情況,單極子主要分布在圓柱分離點附近及尾渦區域,并且圓柱來流側的單極子較大。單極子源項的值分布在0~10左右,數量級很小,可以忽略不計。

圖9依次顯示的是圓柱=0.8、1.0、1.2 m處橫截面上偶極子聲源在圓柱繞流場中的大小分布情況。在圓柱中部三個不同高度截面的偶極子分布情況基本相同,且左右基本對稱。從0°~100°附近,緊貼圓柱表面存在一層較薄的偶極子區域,隨著角度的增大,偶極子區域逐漸變厚;在100°附近偶極子分布區域開始從圓柱表面脫離,并向后形成一個短小的拖曳區域;在100°~180°區域,緊貼圓柱壁面形成一層較薄且不連續的偶極子分布區域。圓柱中部偶極子的最大值數量級為109~1010。

圖8 圓柱中部單極子分布

圖9 圓柱中部偶極子分布

圖10分別顯示的是圓柱=0.8、1.0、1.2 m處橫截面上四極子聲源在圓柱繞流場中的大小分布情況。圓柱中部三個不同高度截面的四極子分布基本相同,且左右基本對稱。在0°~98°緊貼圓柱壁面分布有一層較薄的四極子,在0°和90°附近的四極子相對較薄、較少;在98°~142°區域形成一個比偶極子分布區域更寬更長的四極子分布區域;在142°~180°圓柱表面和后面分布有一些值較小且雜亂的四極子區域。圓柱中部四極子的最大值數量級為108~109。

經過上述分析,在圓柱中部偶極子的強度最大,單極子比偶極子小7~8個數量級,四極子比偶極子小1個數量級。為了更具體地對比圓柱中部偶極子和四極子的大小分布情況,在=1.0 m的橫截面上取=0、25、50、75 mm,從-100~300 mm每隔10 mm取一個點,圓柱區域內的點除外,如圖11所示。對比每個點上偶極子和四極子源項的值,如圖12~15所示。

由圖12~15可知,在=1.0 m截面上緊貼圓柱表面0°和180°處,偶極子源項的值比四極子大1~2個數量級,在圓柱兩側的渦脫落區域從-10~20 mm的范圍內,偶極子比四極子大1~2個數量級,在其他區域,偶極子和四極子大小基本相同。

圖10 圓柱中部四極子分布

圖11 圓柱中部z=1 m截面測點分布

對于圓柱底部,分別選取=0截面和=60 mm截面進行觀察分析。圖16為圓柱底部單極子的大小分布圖,單極子量級為0~10,可以忽略不計。圖17、18分別為圓柱底部偶極子和四極子的大小分布圖,在=60 mm截面處的偶極子和四極子分布與圓柱中部截面分布類似,從=0 mm截面觀察,圓柱底部偶極子和四極子分布與圓柱中部也并沒有什么不同。因此,圓柱底部的偶極子、四極子分布特性與圓柱中部的偶極子、四極子分布特性相同。

圖12 y=0 mm測點的偶極子和四極子大小分布

圖13 y=25 mm測點的偶極子和四極子大小分布

圖14 y=50 mm測點的偶極子和四極子大小分布

圖15 y=75 mm測點的偶極子和四極子大小分布

圖16 圓柱底部單極子大小分布

圖17 圓柱底部偶極子大小分布

圖18 圓柱底部四極子大小分布

3.2 偶極子源分析

偶極子是圓柱繞流的主要聲源,需要對偶極子的分布特征和產生根源做進一步探究。圖19為圓柱頂部速度流線圖,偶極子聲源是由于圓柱對其邊界上的流體的作用力產生的,在圓柱頂部-90°~27°范圍內,來流與圓柱表面相互作用,產生一個緊貼圓柱頂部表面的偶極子區域;接著氣流渦向下回旋,與圓柱后壁面1.63~1.76 m范圍發生碰撞作用,故該區域分布有一層不連續且較厚的偶極子區域。

圖20為=1.7 m截面處的流線圖,來流經過圓柱頂部向下形成湍流渦,二次碰撞圓柱后壁,形成偶極子區域。

圖21為=1.0 m處橫截面上的速度流線圖,在速度流線圖中可以明顯地看到,在圓柱表面0°~100°范圍,來流與圓柱表面碰撞產生作用力,故存在一層緊貼圓柱表面的偶極子區域;在圓柱后側100°~180°范圍,氣流渦二次碰撞圓柱后壁,故產生一層緊貼圓柱后壁的不連續的偶極子區域。

在有限長圓柱繞流中,偶極子作為主要噪聲源。在聲源方程中,偶極子表示為

偶極子源由脈動壓力在x、y、z三個方向的二階梯度構成,為了進一步分析產生偶極子的根源,分別計算脈動壓力在x、y、z三個方向的二階梯度,進行比較分析。圖22為z=1.0 m截面上的脈動壓力云圖,脈動壓力在圓柱兩側的來流分離點附近及其向后形成的拖曳區域沿y方向變化比較明顯。圖23依次為z=1.0 m截面上dp/dx、dp/dy和dp/dz的云圖,dp/dy的值大于dp/dx和dp/dz,且主要分布在圓柱兩側的來流分離點附近及其向后形成的拖曳區域。圖24依次為z=1.0 m截面上d2p/dx2、d2p/dy2和d2p/dz2的云圖,d2p/dy2的云圖與圖9中z=1.0 m 截面上的偶極子分布圖基本相同;d2p/dx2主要分布在緊貼圓柱的來流側,緊貼圓柱后側表面也分布有一層不連續的區域,且分布區域很??;d2p/dz2的值相對于d2p/dx2和d2p/dy2較小,對?2p影響不顯著。因此,??f0主要來源于d2p/dy2,偶極子主要由來流分離處的脈動壓力沿y方向的二階梯度引起。

圖20 z=1.7 m截面流線圖

圖21 z=1.0 m截面流線圖

圖22 圓柱體z=1.0 m截面上的脈動壓力云圖

3.3 四極子源分析

圖24 圓柱體z=1.0 m截面上的d2p/dy2、d2p/dy2、d2p/dy2云圖

圖25 圓柱體z=1.0 m截面上的分布情況

圖26 圓柱體z=1.0 m截面上的分布情況

4 結論

本文通過大渦模擬計算得到三維可壓非定常圓柱繞流場,基于聲源方程從新的角度定量分析圓柱繞流場中單極子、偶極子和四極子聲源的分布特性,更加清晰地認識到聲源的產生根源,為氣動噪聲的控制提供了理論依據。得到的主要結論如下:

(1) 通過大渦模擬和FW-H方程得到的遠場噪聲與風洞試驗結果比較吻合。表明大渦模擬得到的圓柱繞流可壓非定常流場可以為氣動噪聲源的計算提供可靠的聲源信息。

(2) 對于有限長圓柱,偶極子氣動噪聲源占主導地位,單極子源項的值比偶極子低8~9個數量級,四極子源項的值比偶極子低1個數量級左右。

(3) 偶極子源主要分布在圓柱頂部27°來流分離點附近和圓柱中部、底部兩側的分離點附近及圓柱來流側表面,偶極子源主要由脈動壓力沿方向的二階導數d2/d2產生。

(5) 由于計算過程中運用的聲源方程右側源項僅表征了源項的幅值大小,并未考慮聲傳播及聲源輻射效率。關于不同聲源對接收點處總聲場的貢獻率有待進一步研究。

[1] 孫曉峰, 周盛. 氣動聲學[M]. 北京: 國防工業出版社, 1994.

[2] LIGHTHILL M J. On sound generated aerodynamically. i. general theory[J]. Proceedings of the Royal Society of London, 1952, 211(1107): 564-587.

[3] 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.

[4] CURLE N.The Influence of solid boundaries upon aerodynamic sound[C]//Proceedings of the Royal Society of London, 1955, Series A(231): 506-514.

[5] WILLIAMS J E F, HAWKINGS D L. Sound generation by turbulence and surfaces in arbitrary motion[J]. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 1969, 264(1151): 321-342.

[6] HOWE M S. Edge, cavity and aperture tones at very low mach numbers[J]. J. Fluid. Mech, 1997, 330(4): 61-84.

[7] HOWE M S. Acoustics of fluid-structure interactions[M]. Cambridge University Press, 1998.

[8] COX J S, BRENTNER K S, RUMSEY C L. Computation of vortex shedding and radiated sound for a circular cylinder: subcritical to transcritical reynolds numbers[J]. Theoretical & Computational Fluid Dynamics, 1998, 12(4): 233-253.

[9] BRENTNER K S, COX J S, RUMSEY C L, et al. Computation of sound generated by flow over a circular cylinder: an acoustic analogy approach[C]//Computational Aeroacoustics Workshop on Benchmark Problems NASA. 1997: 1-9.

[10] TAKAISHI T, SAGAWA A, NAGAKURA K, et al. Numerical analysis of dipole sound source around high speed trains[J]. J. Acoust. Soc. Am., 2002, 111(6): 2601-2608.

[11] IIDA A, MIZUNO A, KATO C. Visualization of aerodynamic sound source with compact green's function[C]//Aiaa/ceas Aeroacoustics Conference & Exhibit. 2002.

[12] 鄭朝榮, 王笑寒, 武岳. 鈍體繞流氣動噪聲源特性數值研究[J]. 哈爾濱工業大學學報, 2017, 49(12): 146-151.

ZHENG Chaorong, WANG Xiaohan, WU Yue. Numerical investigation on the characteristics of aerodynamic noise sources induced by flows around bluff bodies[J]. Journal of Harbin Institute of Technology, 2017, 49(12): 146-151.

[13] [美]戈德斯坦. 氣動聲學[M]. 閆再友譯. 北京: 國防工業出版社, 2014.

[America]Goldstein. Aeroacoustics[M]. RAN Zaiyou translate. Beijing: National Defense Industry Press, 2014.

[14] 唐狄毅, 李文蘭, 喬渭陽. 飛機噪聲基礎[M]. 西安: 西北工業大學出版社, 1995.

TANG Diyi, LI Wenlan, QIAO Weiyang. Aircraft Noise[M]. Xi’an: Northwestern Polytechnical University Press, 1995.

[15] LEE A H, CAMPBELL R L, HAMBRIC S A. Coupled delayed-detached-eddy simulation and structural vibration of a selt-oscillating cylinder due to vortex-shedding[J]. Journal of Fluids and Structural, 2014, 48(7): 216-234.

[16] 蔡建程, 潘杰, 鄂世舉, 等. 圓柱繞流氣動聲的偶極子及四極子源法定量研究[J]. 聲學學報, 2016, 41(3): 420-427.

CAI Jiancheng, PAN Jie, E Shiju, et al.Quantitive study of the aerodynamic sound induced by the flow past a cylinder based on dipole and quadrupole m odels[J]. Acta Acustica, 2016, 41(3): 420-427.

Study of aeroacoustic noise source induced by a cylindrical flow offinite length

YANG Zhi-gang1,2, LIU Yang1, WANG Yi-gang2

(1. Tongji University, Shanghai Automotive Wind Tunnel Center, Shanghai 201804, China; 2. Beijing Aeronautical Science & Technology Research Institute, Beijing 102211, China)

The energy level, distribution characteristics and origin of aerodynamic noise source are not clear enough. The three-dimensional flow around a finite-length cylinder is taken as the research object. Based on acoustic source equations, the types of aerodynamic noise sources and their relationship with aerodynamic parameters are analyzed. Compressible unsteady flow is calculated numerically and then the aerodynamic parameters are used to quantitatively calculate the sound sources distribution at the top, middle and bottom of the cylinder in order to study the distribution characteristics and origin of the sound sources. Results show that the dipole sound source is dominant in the flow field and the monopole sound source is negligible, the quadrupole sound source is 1-2 orders of magnitude smaller than the dipole sound source, the dipole are mainly distributed at the separation points of incoming flow and the secondary collision zone of eddy flow on the back of cylinder, and the quadrupole is mainly distributed at the separation points of incoming flow and its drag zone. The dipole sound source is mainly caused by the two-step gradient of dynamic pressure in the transverse direction (y-direction) at both sides of the cylinder. The above conclusions can provide a theoretical basis for the aerodynamic noise control.

aeroacoustic noise; cylinder with limited length; sound source equation; dipole sound source

O422.8

A

1000-3630(2019)-01-0005-10

10.16300/j.cnki.1000-3630.2019.01.002

2018-07-12;

2018-08-28

上海市地面交通工具空氣動力與熱環境模擬重點實驗室資助項目、國家自然科學基金資助項目(51375342)。

楊志剛(1961-), 男, 上海人, 博士, 研究方向為汽車空氣動力學。

王毅剛,E-mail: yigang.wang@sawtc.com

猜你喜歡
區域
分割區域
探尋區域創新的密碼
科學(2020年5期)2020-11-26 08:19:22
基于BM3D的復雜紋理區域圖像去噪
軟件(2020年3期)2020-04-20 01:45:18
小區域、大發展
商周刊(2018年15期)2018-07-27 01:41:20
論“戎”的活動區域
敦煌學輯刊(2018年1期)2018-07-09 05:46:42
區域發展篇
區域經濟
關于四色猜想
分區域
公司治理與技術創新:分區域比較
主站蜘蛛池模板: 永久在线精品免费视频观看| 亚洲日本精品一区二区| 免费毛片在线| 国产性爱网站| 成年人国产视频| 中日韩一区二区三区中文免费视频| 精品国产福利在线| 看国产毛片| 一级毛片在线播放| 九九热在线视频| 亚洲一区精品视频在线| 影音先锋亚洲无码| 欧美精品在线看| 国产午夜小视频| 91无码网站| 精品丝袜美腿国产一区| 久久综合亚洲色一区二区三区| 亚洲欧洲日韩久久狠狠爱| 永久免费精品视频| 色偷偷一区二区三区| 精品伊人久久久大香线蕉欧美| 国产视频只有无码精品| 91精品国产丝袜| 国产精品99久久久久久董美香| 免费人欧美成又黄又爽的视频| 亚洲愉拍一区二区精品| 亚洲视频欧美不卡| 白浆免费视频国产精品视频| 國產尤物AV尤物在線觀看| 中文字幕在线欧美| 亚洲精品手机在线| 日韩不卡高清视频| аⅴ资源中文在线天堂| 久久人妻系列无码一区| 久久性视频| 天天综合网色| 亚洲免费毛片| 亚洲日本精品一区二区| 福利姬国产精品一区在线| 狠狠综合久久久久综| 国产91高清视频| 日本成人在线不卡视频| 日韩高清无码免费| 亚洲午夜福利精品无码| 热九九精品| 亚洲欧美h| 国产精品女熟高潮视频| 这里只有精品免费视频| 九色最新网址| 国产一级在线观看www色 | 亚洲人成网7777777国产| 一本久道久久综合多人| 欧美一级黄片一区2区| 亚洲日韩高清在线亚洲专区| 国产chinese男男gay视频网| 亚洲欧美一区二区三区蜜芽| 欧美成人在线免费| 国产成人禁片在线观看| 精品无码国产一区二区三区AV| 国产中文一区a级毛片视频 | 18黑白丝水手服自慰喷水网站| 色欲色欲久久综合网| 日本不卡免费高清视频| 色婷婷丁香| 在线一级毛片| a级毛片一区二区免费视频| 亚洲三级成人| 亚洲综合一区国产精品| 72种姿势欧美久久久大黄蕉| 91免费国产高清观看| 国产精品福利在线观看无码卡| 91原创视频在线| 97色婷婷成人综合在线观看| 日韩欧美中文在线| 国产无遮挡裸体免费视频| 国产免费自拍视频| 亚洲热线99精品视频| 亚洲第七页| 国产网站一区二区三区| 亚洲品质国产精品无码| 婷婷亚洲视频| 妇女自拍偷自拍亚洲精品|