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

        ?

        非均布軸流風(fēng)扇葉片力特性及其尖峰噪聲特性預(yù)測(cè)

        2019-07-03 12:39:58潘丁浩吳亞?wèn)|彭志剛歐陽(yáng)華杜朝輝
        關(guān)鍵詞:尖峰軸流聲壓級(jí)

        潘丁浩, 吳亞?wèn)|, 彭志剛, 歐陽(yáng)華, 杜朝輝

        (1. 上海交通大學(xué) 航空航天學(xué)院, 上海 200240; 2. 上海交通大學(xué) 機(jī)械與動(dòng)力工程學(xué)院, 上海 200240;3. 先進(jìn)航空發(fā)動(dòng)機(jī)協(xié)同創(chuàng)新中心, 北京 100083)

        在傳統(tǒng)葉片均布的低馬赫數(shù)軸流風(fēng)扇中,尖峰噪聲(離散頻率)是讓人感到不適的最重要因素,因此降低風(fēng)扇尖峰噪聲是軸流風(fēng)扇設(shè)計(jì)中的重要問(wèn)題[1-2].對(duì)低速葉輪機(jī)械采用葉片非均布設(shè)計(jì)可以改變流場(chǎng)的頻域特性,繼而有效降低其尖峰噪聲[3].因此,非均布葉片設(shè)計(jì)方法已被應(yīng)用到低速風(fēng)扇的設(shè)計(jì)中[4],并且基于非均布風(fēng)扇的噪聲特性研究也有很多[5].關(guān)于軸流風(fēng)扇葉片非均布降噪的研究主要可分為3個(gè)階段[6].第1階段,Mellin等[7]最早注意到非均布設(shè)計(jì)的降噪潛力,然而他們的研究因?yàn)閷W⒂趯?duì)稱形狀而缺乏一般性.第2階段,Ewald等[8]提出了被后來(lái)廣泛研究的MBS(Modulated Blade Spacing)設(shè)計(jì)規(guī)律,并基于理想點(diǎn)力級(jí)數(shù)展開提出了尖峰噪聲的預(yù)測(cè)方法.Duncan等[9]基于相位調(diào)制理論推導(dǎo)了MBS分布律下的尖峰噪聲降噪模型,并由實(shí)驗(yàn)發(fā)現(xiàn)應(yīng)用MBS的多級(jí)軸流壓氣機(jī)在葉片通過(guò)頻率(Blade Passing Frequency, BPF)處尖峰噪聲降低約10 dB,然而他沒(méi)有對(duì)流場(chǎng)非定常細(xì)節(jié)進(jìn)行研究,也未將實(shí)驗(yàn)結(jié)果與其降噪模型進(jìn)行對(duì)比分析.孫曉峰[10]基于Wright[11]的BLH理論,將葉片點(diǎn)力使用貝塞爾級(jí)數(shù)展開,并以此推導(dǎo)了MBS設(shè)計(jì)方案的尖峰噪聲預(yù)測(cè)公式.第3階段,由Lewy[12]提出的比MBS更簡(jiǎn)便的非均布設(shè)計(jì)方案,并基于相位調(diào)制理論推導(dǎo)了自由聲場(chǎng)中單級(jí)轉(zhuǎn)子的尖峰噪聲預(yù)測(cè)模型.在他的基礎(chǔ)上,Dobrzynski[13]與Anghinolfi等[14]先后提出基于不同聲學(xué)評(píng)估指標(biāo)的非均布最優(yōu)化方法,他們的研究表明,聲學(xué)最優(yōu)的角度往往是比較一般的非均布策略,而非MBS,這也揭示了Lewy工作的必要性.

        近年來(lái),Lewy模型[12]在工業(yè)設(shè)計(jì)中已被廣泛應(yīng)用.然而作為非均布風(fēng)扇聲學(xué)特性預(yù)測(cè)工具,該模型存在的局限性不僅是應(yīng)用對(duì)象只限于自由聲場(chǎng)中單級(jí)轉(zhuǎn)子,其更大的局限性是為了使預(yù)測(cè)模型具有更好的一般性而假設(shè)了各個(gè)葉片所受脈動(dòng)力僅存在由周向角度差所導(dǎo)致的時(shí)域相位差.這一假設(shè)對(duì)非均布設(shè)計(jì)程度影響并不大,并且對(duì)葉片尺寸較小的軸流風(fēng)扇來(lái)說(shuō)也是成立的.然而,對(duì)一般的軸流風(fēng)扇來(lái)說(shuō),這一假設(shè)成立與否,則明顯會(huì)影響到非均布風(fēng)扇尖峰噪聲預(yù)測(cè)的準(zhǔn)確性.風(fēng)扇在快速旋轉(zhuǎn)過(guò)程中,葉片所受氣流非定常力是風(fēng)扇尖峰噪聲的主要來(lái)源.因此,采用非均布設(shè)計(jì)的風(fēng)扇所受脈動(dòng)力特性進(jìn)行分析,對(duì)于建立更準(zhǔn)確的非均布風(fēng)扇尖峰噪聲預(yù)測(cè)模型十分必要.

        基于此,本文針對(duì)某型帶有后方支撐的低速軸流風(fēng)扇采用均布和非均布葉片分布的模型進(jìn)行了非定常數(shù)值模擬,通過(guò)對(duì)葉片采用片條化處理,分析了與非均布角度相關(guān)的葉片力特性.基于非均布風(fēng)扇葉片力的特性以及Lowson點(diǎn)力發(fā)聲模型[15],提出了非均布風(fēng)扇的噪聲預(yù)測(cè)方法,并建立了基于風(fēng)扇葉片力特性的低速軸流非均布風(fēng)扇尖峰噪聲降噪模型.通過(guò)基于計(jì)算流體力學(xué)(CFD)計(jì)算結(jié)果的Lowson模型直接預(yù)測(cè)的結(jié)果和實(shí)驗(yàn)數(shù)據(jù),驗(yàn)證了該降噪模型的有效性.

        1 數(shù)值模擬

        使用ANSYS CFX平臺(tái)對(duì)一低速軸流風(fēng)扇在設(shè)計(jì)流量下進(jìn)行非定常數(shù)值模擬,計(jì)算域如圖1(a)所示.計(jì)算域沿流向依次為:圓柱形入口段、風(fēng)扇段、支撐段及大尺寸半球形出口段.考察3個(gè)風(fēng)扇模型(1個(gè)均布,2個(gè)非均布),3個(gè)模型的風(fēng)扇葉片形狀及下游支撐結(jié)構(gòu)均一致.表1給出了2種非均布方案中葉片的周向相位角.表中:b1~b7表示葉片編號(hào),并默認(rèn)b1周向相位角為0°.在ICEM中對(duì)各算例進(jìn)行結(jié)構(gòu)化網(wǎng)格劃分,經(jīng)過(guò)網(wǎng)格無(wú)關(guān)性驗(yàn)證,整個(gè)計(jì)算域網(wǎng)格總數(shù)約374萬(wàn),圖1(b)和(c)給出了風(fēng)扇與支撐段在葉根部和葉頂部的網(wǎng)格細(xì)節(jié).

        圖1 計(jì)算域及網(wǎng)格Fig.1 CFD simulation domain and mesh

        方案周向相位角/(°)b1b2b3b4b5b6b71064105.41452202573102056109.6145.6213.7260307

        根據(jù)實(shí)驗(yàn)條件,出口邊界設(shè)置為開放邊界,動(dòng)靜交界面采用瞬態(tài)的轉(zhuǎn)/靜模型,即采用k-ω湍流模型.為了考察湍流模型,圖2給出了方案1非均布風(fēng)扇的性能曲線.圖中:橫坐標(biāo)為無(wú)量綱流量系數(shù)φ;縱坐標(biāo)為無(wú)量綱壓力系數(shù)Ψ.由圖可見(jiàn):k-ε湍流模型計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果誤差較大; RNGk-ε、k-ω以及SSTk-ω湍流模型計(jì)算結(jié)果比較接近;k-ω湍流模型與實(shí)驗(yàn)結(jié)果更接近,所以選擇了k-ω湍流模型進(jìn)行數(shù)值計(jì)算.

        圖2 風(fēng)扇性能曲線Fig.2 Aerodynamic performance of the uneven fan

        圖3 葉片片條Fig.3 Strips of blade

        由于本風(fēng)扇葉片前掠明顯,并且輪轂比較大,為了更準(zhǔn)確地分析葉片受力,將風(fēng)扇葉片沿徑向分為10個(gè)等寬的片條.片條編號(hào)1#~10#依次代表葉根至葉頂10個(gè)片條,如圖3所示.認(rèn)為各片條所受合力在聲源特性上等效為一個(gè)作用在片條形心的運(yùn)動(dòng)點(diǎn)力.將片條在隨體系中受到的等效點(diǎn)力按方向定義為FT、FD、FL,如圖4所示.圖中:黑色圓點(diǎn)代表某繞著x軸旋轉(zhuǎn)的風(fēng)扇片條;ω為旋轉(zhuǎn)角速度;FT為軸向力,指向上游方向?yàn)檎籉L為徑向力,離心方向?yàn)檎?;FD為切向阻力,與旋轉(zhuǎn)方向相反;R為旋轉(zhuǎn)半徑.

        圖4 片條力方向性定義Fig.4 The definition of point force

        2 風(fēng)扇葉片力特性

        2.1 風(fēng)扇葉片力方向特性

        圖5 方案1中的b1葉片各片條脈動(dòng)力時(shí)均值Fig.5 The time averaged of blade force of b1

        圖6 b1葉片脈動(dòng)力脈動(dòng)特性Fig.6 Force fluctuation in blade b1

        2.2 非均布風(fēng)扇葉片力時(shí)均特性

        圖7給出了各片條軸向力FT的時(shí)均值.與均布風(fēng)扇不同,非均布風(fēng)扇各葉片受力時(shí)均值存在差異.將圖7所有片條所受的FT力求和,均布結(jié)果為 13.43 N,非均布結(jié)果為 13.33 N,兩者相差不到1%.由此可見(jiàn),非均布風(fēng)扇所有葉片受力時(shí)均值的總和保持不變.

        圖7 各風(fēng)扇片條軸向力FT時(shí)均值Fig.7 Time-averaged values of FT

        圖8 各風(fēng)扇片條無(wú)量綱時(shí)均值Fig.8 Normalized Time-averaged values of

        (1)

        式中:A=0.239;B=0.007 35;C=0.383.擬合結(jié)果如圖10所示.由圖可見(jiàn),擬合公式對(duì)采用不同角度方案的非均布風(fēng)扇葉片F(xiàn)T的預(yù)測(cè)具備很好的通用性.

        圖9 75%葉高處壓力云圖Fig.9 Pressure contour at 75% span

        圖10 擬合結(jié)果Fig.10 Fitting results

        2.3 非均布風(fēng)扇葉片力脈動(dòng)項(xiàng)特性

        圖11 方案1風(fēng)扇b1葉片7#片條FT脈動(dòng)項(xiàng)Fig.11 FT Fluctuations in the 7# segment

        3 非均布風(fēng)扇尖峰噪聲預(yù)測(cè)

        3.1 Lowson點(diǎn)力發(fā)聲模型

        Lowson模型可有效預(yù)測(cè)低馬赫數(shù)運(yùn)動(dòng)點(diǎn)力在自由聲場(chǎng)中聲場(chǎng)特性[15].本例計(jì)算域中流動(dòng)馬赫數(shù)小于 0.2,同時(shí)支撐管道極短,可看作自由聲場(chǎng)問(wèn)題;其次,本文對(duì)風(fēng)扇采用片條處理,各片條自身受力形式較同一,因此認(rèn)為作用在各片條形心的等效運(yùn)動(dòng)點(diǎn)力的發(fā)聲特性與風(fēng)扇是等效的.另外,本文考察單倍軸頻至2BPF頻域范圍內(nèi)的噪聲(BPF為葉片通過(guò)頻率), 2BPF諧波波長(zhǎng)約為 0.57 m,大于風(fēng)扇的直徑,低頻諧波波長(zhǎng)則更長(zhǎng).考慮緊湊聲源的定義,考察頻域的聲波波長(zhǎng)大于葉輪直徑,認(rèn)為本風(fēng)扇符合緊湊聲源的條件.基于上述討論,本文運(yùn)用Lowson模型進(jìn)行風(fēng)扇的尖峰噪聲預(yù)測(cè).

        Lowson[15]推導(dǎo)了任意低馬赫數(shù)運(yùn)動(dòng)脈動(dòng)點(diǎn)力自由聲場(chǎng)的聲壓

        (2)

        式中:yi、xi、Fi分別為觀察點(diǎn)位置、脈動(dòng)點(diǎn)力位置和脈動(dòng)力;a0為當(dāng)?shù)芈曀?;r=yi-xi;Ma為脈動(dòng)點(diǎn)力在r方向上速度分量的馬赫數(shù).

        由式(2)推得聲諧波的復(fù)數(shù)值

        (3)

        式中:Fr=(yi-xi)Fi/r;Mai為點(diǎn)力運(yùn)動(dòng)各向馬赫數(shù).

        在n階轉(zhuǎn)頻(RF)處的聲諧波聲壓幅值pnRF=|cn|,由聲壓級(jí)公式SPL=20 lg(p/p0)得到諧波聲壓級(jí).

        3.2 基于葉片力特性的非均布風(fēng)扇尖峰噪聲預(yù)測(cè)方法

        傳統(tǒng)非均布風(fēng)扇尖峰噪聲預(yù)測(cè)方法忽略了非均布風(fēng)扇各葉片受力的差異特性.基于上面討論的非均布風(fēng)扇葉片力特性,提出一種利用Lowson點(diǎn)力模型進(jìn)行尖峰噪聲預(yù)測(cè)的方法步驟:

        (3) 對(duì)葉片數(shù)為7的非均布風(fēng)扇,解方程組

        (4) 葉片力時(shí)均值所致聲諧波使用無(wú)量綱力來(lái)修正,而脈動(dòng)項(xiàng)所致聲諧波不變,因此第j個(gè)葉片所受合力輻射的聲諧波復(fù)數(shù)值

        (4)

        可得n倍軸頻處聲壓級(jí)

        (5)

        步驟(1)中進(jìn)行CFD計(jì)算是為了獲得風(fēng)扇葉片力脈動(dòng)項(xiàng)的具體形式,若將風(fēng)扇葉片力在頻域上作傅里葉展開,即F=Fλe-iλωt(F=FT、FD、FL),可推得葉片均布時(shí)單片風(fēng)扇遠(yuǎn)場(chǎng)聲諧波預(yù)測(cè)式

        (6)

        式中:α=ωRy/(a0r1);y1和y2分別為觀測(cè)點(diǎn)的軸向位置及到軸向的距離.需要指出,如果忽略徑向力FL而只考慮FT和FD,則式(6)可以進(jìn)一步推導(dǎo)得與Lowson模型后續(xù)推導(dǎo)一致的包含貝塞爾函數(shù)的形式.

        3.3 基于葉片力特性的非均布風(fēng)扇尖峰噪聲降噪模型

        理論上均布風(fēng)扇的尖峰噪聲僅在葉片通過(guò)頻率倍頻處有分量[15],因此考慮非均布風(fēng)扇在BPF倍頻處相對(duì)均布風(fēng)扇下降的聲壓級(jí),即

        ΔSPL=20 lg(pue/pe)

        其中:pue和pe分別為非均布和均布的聲壓,代入式(5),可得

        ΔSPL=

        (7)

        (8)

        式中:B為葉片數(shù).μk通過(guò)均布風(fēng)扇的CFD計(jì)算獲得,它表示單只葉片分別由時(shí)均載荷及脈動(dòng)載荷輻射聲諧波在某測(cè)點(diǎn)處的幅值強(qiáng)弱及相位關(guān)系.從定義式上來(lái)看,μk受聲諧波測(cè)定值的影響,即受到測(cè)點(diǎn)距離及動(dòng)靜干涉具體形式的影響.可見(jiàn),μk取值的規(guī)律至少受到3個(gè)方面因素的影響:觀測(cè)點(diǎn)位置、風(fēng)扇形狀和支撐特點(diǎn)(支撐形狀及安裝方式).因此,對(duì)特定的風(fēng)扇及支撐組件,有必要通過(guò)CFD確定μk的取值.

        圖12 各倍軸頻μkFig.12 Complex values of μk

        值得一提的是,對(duì)于本例風(fēng)扇以及選定的測(cè)點(diǎn),5倍及以上軸頻μk模值非常小,如圖12所示.其原因是在該測(cè)點(diǎn)定常力輻射的高倍軸頻聲諧波較弱,而支撐間形狀差異導(dǎo)致葉片力脈動(dòng)項(xiàng)在各倍軸頻處均存在明顯成分.這種情況下,式(8)近似等價(jià)于

        (9)

        式(9)與Lewy[12]的分析一致,其優(yōu)勢(shì)在于無(wú)需任何數(shù)值模擬的前提計(jì)算,即可以得到非均布風(fēng)扇在BPF倍頻處的降噪特性.但是,式(9)僅對(duì)高倍軸頻μk是小量時(shí)成立.對(duì)于不同的風(fēng)扇及支撐構(gòu)型,μk取值的規(guī)律并不確定(見(jiàn)圖12).

        基于式(4)推導(dǎo)出的式(8)給出了任意非均布風(fēng)扇相比于均布風(fēng)扇在葉片通過(guò)頻率倍頻的降噪模型.實(shí)際上,對(duì)于軸頻倍頻,也可以通過(guò)式(4)評(píng)估任意兩個(gè)采取特定非均布角度方案的風(fēng)扇在該頻率處的聲壓級(jí)差,即

        ΔSPLn=

        (10)

        由式(10)可知,基于某型風(fēng)扇采用任一非均布角度方案時(shí)的聲學(xué)測(cè)試結(jié)果,即可預(yù)測(cè)該型風(fēng)扇采取其他任意非均布角度方案時(shí)的尖峰頻率聲壓級(jí),即

        SPLue1=SPLue2,exp+ΔSPLn

        其中:SPLue2,exp為由試驗(yàn)測(cè)得的非均布在倍軸頻處的聲壓級(jí),而ΔSPLn由式(10)得到.

        3.4 降噪模型效果評(píng)估

        基于CFD得到的葉片力結(jié)果,使用 3.1 節(jié)中的Lowson模型進(jìn)行直接尖峰噪聲預(yù)測(cè),監(jiān)測(cè)點(diǎn)選取為旋轉(zhuǎn)軸夾角45° 下游方向,距離風(fēng)扇1 m處.圖13所示為該位置預(yù)測(cè)的聲壓級(jí)譜.由圖可見(jiàn),非等距風(fēng)扇相對(duì)等距風(fēng)扇在BPF處的聲壓級(jí)減少了 4.1 dB,而式(8)預(yù)測(cè)的結(jié)果為 3.92 dB,兩者吻合.另外,圖中2BPF處的聲壓級(jí)降低了 15.7 dB,而式(8)預(yù)測(cè)的結(jié)果為 15.45 dB,兩者也吻合.由此可見(jiàn),3.3 節(jié)中考慮葉片力特性的非均布風(fēng)扇降噪預(yù)測(cè)模型與基于CFD結(jié)果的直接聲學(xué)預(yù)測(cè)十分吻合.

        圖13 聲壓級(jí)譜Fig.13 Predicted SPL spectrum

        針對(duì)本文中使用的2套非均布角度方案,在標(biāo)準(zhǔn)聲學(xué)的實(shí)驗(yàn)條件下,得到了BPF及2BPF處尖峰噪聲的聲壓級(jí),并減去均布風(fēng)扇的試驗(yàn)值,得到葉片通過(guò)頻率倍頻降噪量如表2所示.表中,同時(shí)比較了目前應(yīng)用廣泛的離散噪聲預(yù)測(cè)模型得出的結(jié)果[9],通過(guò)將CFD得到的葉片力直接應(yīng)用Lowson模型得到的結(jié)果,以及 3.3 節(jié)中降噪模型所計(jì)算得到的結(jié)果.

        表2 2種非均布方案相對(duì)均布風(fēng)扇的降噪比較Tab.2 Different noise models and experiment results

        相比Duncun[9]模型,基于CFD結(jié)果的混合方法能更好地預(yù)測(cè)非均布風(fēng)扇的尖峰噪聲,而本文提出的基于葉片力特性的尖峰噪聲降噪模型與混合方法吻合很好,因此可以節(jié)省非均布設(shè)計(jì)時(shí)對(duì)任意角度采用混合方法的過(guò)程成本.

        本文模型相比于Lewy模型來(lái)說(shuō),考慮了風(fēng)扇非均布時(shí)葉片力的特性,因此對(duì)于支撐形狀一致性好且動(dòng)靜干涉強(qiáng)烈的風(fēng)扇、或不存在支撐的風(fēng)扇、或測(cè)點(diǎn)位于對(duì)于空間中某些特殊位置時(shí),即μk值較大時(shí),式(8)相對(duì)于Lewy模型更準(zhǔn)確.另外,Lewy的模型是基于自由聲場(chǎng)中單級(jí)轉(zhuǎn)子,可以認(rèn)為此時(shí)各葉片受力為定常力,則式(8)等價(jià)于下式

        (11)

        4 結(jié)論

        本文以低馬赫數(shù)汽車?yán)鋮s軸流風(fēng)扇為例,根據(jù)CFD風(fēng)扇葉片周向非均布時(shí)葉片力的特性,提出了低速軸流非均布風(fēng)扇尖峰噪聲預(yù)測(cè)的方法,并建立了低速軸流非均布風(fēng)扇尖峰頻率降噪模型.通過(guò)分析,得到的主要結(jié)論如下:

        (1) 低速軸流非均布風(fēng)扇葉片力的時(shí)均值與具體非均布角度方案有關(guān);各葉片受力脈動(dòng)值僅存在與風(fēng)扇周向分布角度相關(guān)的時(shí)域上的相位差,不存在脈動(dòng)幅值和形式上的差異,并且各葉片受力時(shí)均值的差異規(guī)律可以通過(guò)CFD結(jié)果擬合得到.

        (2) 提出了非均布葉片噪聲影響的關(guān)鍵參數(shù),即葉片周向分布規(guī)律,以及葉片所受時(shí)均載荷與脈動(dòng)載荷所致聲諧波的比值μk,μk與測(cè)點(diǎn)位置、風(fēng)扇形狀、以及支撐形式有關(guān).μk?1時(shí),基于非均布風(fēng)扇葉片力特性的尖峰預(yù)測(cè)模型與Lewy的非均布風(fēng)扇噪聲預(yù)測(cè)模型一致;此條件不成立時(shí),本文所提的預(yù)測(cè)模型更為準(zhǔn)確.

        (3) 本文對(duì)低速軸流非均布風(fēng)扇提出了考慮葉片力特性的尖峰降噪模型,它可以有效評(píng)估某軸流風(fēng)扇葉片采用任意非均布角度方案時(shí)的尖峰噪聲降噪性能,其結(jié)果與基于CFD的Lowson點(diǎn)力模型得到的預(yù)測(cè)結(jié)果吻合.

        猜你喜歡
        尖峰軸流聲壓級(jí)
        機(jī)器噪聲平均聲壓級(jí)計(jì)算方法差異性實(shí)證研究
        軸流壓氣機(jī)效率評(píng)定方法
        尖峰石陣
        一種計(jì)算消聲室聲壓級(jí)的新方法
        全新DXR mkll有源揚(yáng)聲器
        演藝科技(2019年4期)2019-03-30 03:21:46
        西澳大利亞——尖峰石陣
        超大型軸流風(fēng)機(jī)在十天高速西秦嶺隧道中的應(yīng)用
        微型軸流風(fēng)機(jī)低溫失效分析
        Diodes1.9W D類音頻放大器提供高聲壓級(jí)水平并延長(zhǎng)電池壽命
        尖峰之年:NASA地球科學(xué)探測(cè)進(jìn)入高潮
        太空探索(2014年4期)2014-07-19 10:08:58
        日本av天堂一区二区三区| 久久精品无码专区东京热| 亚洲精品无码成人片久久不卡 | 激情久久无码天堂| av天堂手机一区在线| 蜜桃av在线免费网站| 久精品国产欧美亚洲色aⅴ大片| 国产成人亚洲日韩欧美| 伊人婷婷色香五月综合缴激情| 五月激情在线观看视频| 日本xxxx色视频在线观看免费| 东北妇女肥胖bbwbbwbbw| 免费 无码 国产在线观看不卡| 国产视频在线播放亚洲| 国产女人精品视频国产灰线| 狠狠色噜噜狠狠狠888米奇视频 | 亚洲国产综合精品久久av| 中文字幕日本五十路熟女| 日韩精品熟妇一区二区三区| 中文字幕一区在线观看视频| 国产精品一区二区久久乐下载| 97超碰国产一区二区三区 | 国产白嫩护士被弄高潮| 国产精品一区二区电影| 免费观看在线一区二区| 中文字幕一区二区三区久久网| 久久亚洲精品无码va白人极品| 永久免费观看的毛片手机视频 | 国产一级自拍av播放| 华人在线视频精品在线| 国产乱子伦农村xxxx| 天天做天天躁天天躁| 欧洲AV秘 无码一区二区三| 亚洲中文字幕精品久久吃奶| 久久久精品一区aaa片| 中日韩欧美在线观看| 国产成人久久精品二区三区| 体验区试看120秒啪啪免费| 竹菊影视欧美日韩一区二区三区四区五区| 国产女精品| 日韩精品一区二区三区免费观影 |