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

        ?

        風(fēng)扇激波相關(guān)噪聲預(yù)測(cè)方法及試驗(yàn)驗(yàn)證

        2021-08-27 06:46:32蔣永松董洪瑞鄭文濤王詠梅
        航空發(fā)動(dòng)機(jī) 2021年4期
        關(guān)鍵詞:模態(tài)

        蔣永松,董洪瑞,鄭文濤,王詠梅

        (中國(guó)航發(fā)沈陽(yáng)發(fā)動(dòng)機(jī)研究所,沈陽(yáng)110015)

        0 引言

        隨著民用航空發(fā)動(dòng)機(jī)涵道比的不斷增大,噴流噪聲大幅降低,而風(fēng)扇噪聲所占比重卻隨之增加,因此,風(fēng)扇噪聲的產(chǎn)生機(jī)理和傳播特性及其控制仍然是民用航空重要的研究?jī)?nèi)容[1-3]。針對(duì)風(fēng)扇噪聲源的控制,最常用的方法是根據(jù)“Cut-off”準(zhǔn)則[4]選取轉(zhuǎn)子和靜子葉片數(shù),減少可傳播的聲模態(tài)。另外降低風(fēng)扇速度也可有效降低噪聲,然而為了減輕風(fēng)扇設(shè)計(jì)的難度,風(fēng)扇速度不可能取得很低,設(shè)計(jì)時(shí)不可避免地會(huì)使風(fēng)扇轉(zhuǎn)子部分葉高處于超聲速工況,此時(shí)會(huì)在葉片前緣或通道內(nèi)形成激波,并隨之產(chǎn)生與激波相關(guān)的噪聲。

        為了研究激波噪聲的傳播規(guī)律,Morfey和Fish?er[5]根據(jù)弱激波理論,提出了可以預(yù)測(cè)理想葉柵的1維傳播模型,該模型存在大量的簡(jiǎn)化和假設(shè),并且需要給定激波初始強(qiáng)度;Fisher[6]從弱激波的黎曼方程出發(fā),提出了激波在硬壁管道中傳播的時(shí)域計(jì)算方法(Time Domain Numerical Solution,TDNS);隨后McAl?pine和Fisher[7-8]又提出了可以考慮聲襯等軟壁面影響的頻域計(jì)算方法(Frequency Domain Numerical So?lution,F(xiàn)DNS)。但無(wú)論是TDNS還是FDNS,都需要給定初始激波強(qiáng)度,而這又很難精確得到。實(shí)際上,計(jì)算流體力學(xué)(Computational Fluid Dynamics,CFD)已在風(fēng)扇/壓氣機(jī)領(lǐng)域廣泛運(yùn)用,通過(guò)對(duì)風(fēng)扇/壓氣機(jī)進(jìn)行CFD仿真,可以對(duì)其內(nèi)部流場(chǎng)進(jìn)行精細(xì)解讀,詳細(xì)分析激波產(chǎn)生過(guò)程和向上游傳播的特征,進(jìn)而研究激波相關(guān)噪聲。Prasad[9]就曾對(duì)比分析了CFD仿真中使用Navier-Stokes方程和Euler方程在2維和3維算例中捕捉激波的能力。對(duì)于風(fēng)扇轉(zhuǎn)子而言,相對(duì)坐標(biāo)系下靜壓的周向分布與絕對(duì)坐標(biāo)系下固定觀測(cè)點(diǎn)的靜壓脈動(dòng)存在等價(jià)關(guān)系,而前傳噪聲的大小在一定程度上直接取決于風(fēng)扇前靜壓的周向分布,因此,可以用定常CFD仿真來(lái)分析激波相關(guān)噪聲[10]。唐慧敏等[11]根據(jù)PW公司Prasad[12]的思路,發(fā)展了1套風(fēng)扇激波噪聲程序,分析了不同差分格式和網(wǎng)格密度對(duì)噪聲預(yù)測(cè)結(jié)果的影響。RR公司的Wilson則將自己發(fā)展的波分解方法[13]應(yīng)用到歐盟SILENCER項(xiàng)目的跨聲速轉(zhuǎn)子多學(xué)科優(yōu)化設(shè)計(jì)中。

        本文以大涵道比渦扇發(fā)動(dòng)機(jī)風(fēng)扇/增壓級(jí)部件試驗(yàn)件作為研究對(duì)象,采用CFD定常仿真和波分解方法分析激波相關(guān)噪聲的產(chǎn)生和傳播特征,并與試驗(yàn)結(jié)果進(jìn)行對(duì)比,驗(yàn)證方法的可靠性和工程適用性。

        1 聲學(xué)評(píng)估方法

        Morfey[14]提出了聲強(qiáng)的計(jì)算公式

        穿過(guò)某個(gè)截面S的聲功率為

        針對(duì)本文研究的問(wèn)題,在靜止坐標(biāo)系下的時(shí)間脈動(dòng)完全由轉(zhuǎn)子的轉(zhuǎn)動(dòng)引起,任意變量的時(shí)間平均值等價(jià)于變量的周向平均值。因此,風(fēng)扇進(jìn)口某軸向截面S的聲功率可寫(xiě)為

        式中:B為葉片數(shù)為x方向平均流速度;rh和rt分別為流道內(nèi)徑和外徑。

        如果直接利用式(3)對(duì)風(fēng)扇噪聲進(jìn)行評(píng)估,會(huì)存在短波長(zhǎng)的震蕩,并且震蕩幅度隨激波強(qiáng)度增大而增大。Wilson[13]指出這是由進(jìn)口處的數(shù)值反射產(chǎn)生了短波長(zhǎng)的數(shù)值震蕩導(dǎo)致的,為了消除該非物理震蕩,本文采用Wilson提出的波分解法提取由激波引起的相關(guān)單音噪聲,詳細(xì)分析其產(chǎn)生和傳播特征,其基本原理如下。

        在沿x方向均勻流條件下,小擾動(dòng)聲波傳播控制方程為

        式中:c0為當(dāng)?shù)仄骄髀曀俸蛈均由c0無(wú)量綱化;ρ為脈動(dòng)密度,由平均流密度無(wú)量綱化;p由γp0無(wú)量綱化,其中p0為平均大氣壓。

        由方程(4)求時(shí)間的導(dǎo)數(shù)減去對(duì)方程(5)求散度,并利用線(xiàn)化的物態(tài)方程得到對(duì)流波動(dòng)方程

        對(duì)于在內(nèi)、外徑分別為rh和rt的平行環(huán)形管道內(nèi)以頻率ω傳播的聲擾動(dòng),在圓柱坐標(biāo)系下方程(6)的通解為

        式中:m和n分別為周向和徑向模態(tài)數(shù)分別為向下游和向上游傳播的軸向波數(shù);特征函數(shù)分 別 為 待 定 系數(shù)分別為第1、2類(lèi)貝塞爾函數(shù),且特征函數(shù)ψmn和特征值μmn、軸向波數(shù)λmn滿(mǎn)足方程

        以及頻散關(guān)系

        式中:k=ω c0;λmn為軸向波數(shù)

        對(duì)于硬壁圓環(huán)管道,可由對(duì)應(yīng)的Bessel方程和邊界條件求出對(duì)應(yīng)特征值μmn。為方便起見(jiàn),令特征函數(shù)ψmn中的系數(shù)A=1,則系數(shù)因此有

        求得徑向特征值μmn后,可以求出其軸向波數(shù)

        式中:

        實(shí)際上,對(duì)于無(wú)旋流,無(wú)論是Cut-on模態(tài)還是Cut-off模態(tài),式(13)所代表的2個(gè)根都分別與前傳波和后傳波對(duì)應(yīng),熵波和渦波對(duì)于壓力都沒(méi)有貢獻(xiàn)。在平行圓環(huán)管道內(nèi),對(duì)于給定的頻率ω、周向模態(tài)和徑向模態(tài)(m,n),可將壓力和軸向速度脈動(dòng)分解為向上游和下游傳播2部分[13],即

        其中對(duì)于給定的軸向截面,如果已知p( )x,r,θ的分布,則可以利用傅里葉級(jí)數(shù)和貝塞爾函數(shù)的正交性,通過(guò)逆變換得到Pmn(x),即

        式中:Γ為特征函數(shù)ψmn的正交值,即

        相應(yīng)地,速度ux的傅里葉和貝塞爾變換系數(shù)

        2 數(shù)值分析

        本文以大涵道比渦扇發(fā)動(dòng)機(jī)風(fēng)扇/增壓級(jí)部件試驗(yàn)件作為研究對(duì)象,該試驗(yàn)件由風(fēng)扇轉(zhuǎn)子、外涵出口導(dǎo)向葉片(Outlet Guide Vane,OGV)和帶有進(jìn)口導(dǎo)向葉片的4級(jí)增壓級(jí)組成,其中風(fēng)扇轉(zhuǎn)子葉片數(shù)為22,風(fēng)扇轉(zhuǎn)子葉尖進(jìn)口相對(duì)馬赫數(shù)為1.48。

        本文使用多塊結(jié)構(gòu)化氣動(dòng)數(shù)值模擬程序Map[16]進(jìn)行3維定常流場(chǎng)的計(jì)算,控制方程組的空間離散采用以單元為中心的有限體積方法,對(duì)流通量采用低耗散的通量分裂格式求解,并利用MUSCL提高插值精度,采用中心差分格式求解黏性通量,使用隱式方法進(jìn)行時(shí)間離散,并運(yùn)用Gauss-Seidel算法求解,在對(duì)湍流的計(jì)算中,采用Spalart-Allmaras湍流模型對(duì)控制方程進(jìn)行封閉。計(jì)算網(wǎng)格如圖1所示。網(wǎng)格點(diǎn)總數(shù)為787萬(wàn),其中風(fēng)扇轉(zhuǎn)子域網(wǎng)格點(diǎn)數(shù)為597萬(wàn),軸向每個(gè)波長(zhǎng)的網(wǎng)格點(diǎn)約為20個(gè),以保證分辨率[12],為減少計(jì)算域邊界反射效應(yīng)對(duì)計(jì)算結(jié)果的影響,計(jì)算域進(jìn)口選在風(fēng)扇轉(zhuǎn)子上游約8倍軸向弦長(zhǎng)處,并且網(wǎng)格的密度隨著遠(yuǎn)離風(fēng)扇前緣逐漸稀疏。

        圖1 計(jì)算網(wǎng)格

        風(fēng)扇/增壓級(jí)的外涵特性如圖2所示。圖中除給出了Map的仿真結(jié)果外,還給出了對(duì)應(yīng)風(fēng)扇/增壓級(jí)部件試驗(yàn)件的試驗(yàn)結(jié)果。在風(fēng)扇噪聲分析過(guò)程中,主要針對(duì)相對(duì)轉(zhuǎn)速n=0.700、0.800、0.935、0.957、1.000進(jìn)行,首先選取靠近共同工作線(xiàn)的狀態(tài)點(diǎn)進(jìn)行噪聲分析。共同工作線(xiàn)上風(fēng)扇噪聲向上游傳播時(shí)總聲功率的變化如圖3所示。圖中總聲功率只考慮了第1、2階葉片通過(guò)頻率(Blade Passing Frequency,BPF)上噪聲的貢獻(xiàn),聲功率用最大值進(jìn)行了無(wú)量綱化,長(zhǎng)度的單位用葉尖弦長(zhǎng)Btip進(jìn)行了無(wú)量綱化。眾所周知,對(duì)于非線(xiàn)性聲波,特別是由激波引起的噪聲,在傳播過(guò)程中聲波的幅值和聲功率都會(huì)隨著傳播距離的增加而逐漸衰減。從圖中可見(jiàn),除n=0.7外,所有狀態(tài)的聲功率近似與風(fēng)扇前緣的軸向距離呈線(xiàn)性變化;當(dāng)n=0.8時(shí),噪聲衰減最慢;且從轉(zhuǎn)速n=0.8升到n=1.0,衰減率逐漸增大。在n=0.7時(shí),聲功率快速衰減并發(fā)生大幅振蕩,主要原因在于此時(shí)葉尖處于亞聲速狀態(tài),所有聲模態(tài)均被截止,壓力波以指數(shù)規(guī)律衰減。

        圖2 風(fēng)扇/增壓級(jí)外涵特性

        圖3 共同工作線(xiàn)上風(fēng)扇噪聲總聲功率

        風(fēng)扇前3個(gè)不同軸向位置(與風(fēng)扇前緣的距離分別為1.0、1.8、3.0倍葉尖軸向弦長(zhǎng))不同狀態(tài)下,風(fēng)扇聲功率隨風(fēng)扇進(jìn)口葉尖相對(duì)馬赫數(shù)的變化曲線(xiàn)如圖4所示。圖中進(jìn)口葉尖相對(duì)馬赫數(shù)Max為軸向平均馬赫數(shù),Matip為轉(zhuǎn)子葉尖馬赫數(shù)。從圖中可見(jiàn),在距離風(fēng)扇前緣1.8和3.0倍弦長(zhǎng)處,在n=0.957和n=1.000狀態(tài)下,風(fēng)扇噪聲分別比在距離風(fēng)扇前緣1.0倍弦長(zhǎng)處的噪聲降低約5 dB和10 dB。從這些曲線(xiàn)可以得出與激波相關(guān)的傳播特性,當(dāng)MaFrel≈1.12時(shí),噪聲達(dá)到峰值;當(dāng)MaFrel<1.12時(shí),隨著進(jìn)口葉尖相對(duì)馬赫數(shù)的增大,噪聲值急劇增大;當(dāng)MaFrel>1.12時(shí),隨著MaFrel的增大噪聲反而降低,只是降低速率較為平緩??傊?,風(fēng)扇的噪聲與來(lái)流相對(duì)馬赫數(shù)并不成正比關(guān)系,在葉尖馬赫數(shù)略大于1.1時(shí)達(dá)到最大值,之后隨著葉尖馬赫數(shù)進(jìn)一步增大反而減小。

        圖4 不同軸向位置聲功率隨風(fēng)扇進(jìn)口葉尖相對(duì)馬赫數(shù)的變化

        為了對(duì)各狀態(tài)的風(fēng)扇噪聲特性進(jìn)行精細(xì)分析,給出了各轉(zhuǎn)速工作點(diǎn)上風(fēng)扇典型流場(chǎng)的計(jì)算結(jié)果,如圖5~9所示。從轉(zhuǎn)子前無(wú)量綱脈動(dòng)壓力圖可見(jiàn),隨著轉(zhuǎn)速的降低,相應(yīng)BPF頻率降低,對(duì)應(yīng)脈動(dòng)壓力波的波長(zhǎng)逐漸變大,而且當(dāng)n=0.8時(shí),壓力波以第1階徑向模態(tài)形式在管道內(nèi)傳播,此時(shí)壓力波能向上游傳播更遠(yuǎn)距離,這也可以從圖3中看出,此時(shí)風(fēng)扇噪聲最強(qiáng)。

        圖5 在n=1.0工作點(diǎn)上的風(fēng)扇噪聲分析

        圖6 在n=0.957工作點(diǎn)上的風(fēng)扇噪聲分析

        圖7 在n=0.935工作點(diǎn)上的風(fēng)扇噪聲分析

        圖8 在n=0.8工作點(diǎn)上的風(fēng)扇噪聲分析

        圖9 在n=0.7工作點(diǎn)上風(fēng)扇噪聲分析

        當(dāng)風(fēng)扇工作在設(shè)計(jì)轉(zhuǎn)速工作點(diǎn)時(shí),轉(zhuǎn)子葉尖形成了壓氣機(jī)典型的雙波系結(jié)構(gòu),進(jìn)口一道斜激波,緊隨其后一道通道激波。在起飛轉(zhuǎn)速下,葉尖的激波恰好被推出槽道,隨著轉(zhuǎn)速的降低,激波進(jìn)一步靠近葉片上游位置,并且激波與軸向的夾角變小。從圖4中可見(jiàn),盡管隨著轉(zhuǎn)速的降低,風(fēng)扇的增壓能力和負(fù)荷水平減小,但由于激波被推出槽道,噪聲水平反而更高。實(shí)際上,正如Prasad等[15]指出的,由于風(fēng)扇葉尖區(qū)域的激波系結(jié)構(gòu)與來(lái)流的相對(duì)馬赫數(shù)相關(guān),當(dāng)來(lái)流馬赫數(shù)達(dá)到約1.2時(shí),風(fēng)扇前緣的激波越斜,與軸向夾角增大,激波強(qiáng)度降低,其引起的風(fēng)扇前傳噪聲也相應(yīng)降低,這也可從各轉(zhuǎn)速下葉尖相對(duì)馬赫數(shù)分布中看出。

        3 試驗(yàn)分析

        為了進(jìn)行激波相關(guān)噪聲數(shù)值仿真結(jié)果與試驗(yàn)的對(duì)比分析,開(kāi)展了風(fēng)扇/增壓級(jí)在全消聲環(huán)境條件下的聲學(xué)試驗(yàn)。試驗(yàn)時(shí)在風(fēng)扇前1倍風(fēng)扇直徑距離內(nèi)布置了2列聲學(xué)傳感器,如圖10所示。結(jié)合試驗(yàn)時(shí)試驗(yàn)件的運(yùn)行狀況,并參考典型聲學(xué)分析狀態(tài)點(diǎn)(如圖11所示),分別錄取了當(dāng)內(nèi)涵處于對(duì)應(yīng)轉(zhuǎn)速工作點(diǎn)時(shí),n=0.800、0.935外涵堵點(diǎn)和工作點(diǎn)的聲學(xué)性能,n=0.957、1.000外涵堵點(diǎn)、工作點(diǎn)和近喘點(diǎn)的聲學(xué)性能。

        圖10 風(fēng)扇前管道內(nèi)聲學(xué)傳感器測(cè)點(diǎn)分布

        圖11 聲學(xué)分析狀態(tài)點(diǎn)

        不同轉(zhuǎn)速下不同狀態(tài)點(diǎn)計(jì)算得到的聲壓級(jí)沿軸向的分布與試驗(yàn)結(jié)果的對(duì)比如圖12所示。圖中聲壓級(jí)用最大總聲壓級(jí)進(jìn)行了無(wú)量綱化。總體而言,除風(fēng)扇轉(zhuǎn)子前緣附近外,計(jì)算的聲壓級(jí)在低轉(zhuǎn)速與試驗(yàn)值吻合得非常好,當(dāng)n=0.800時(shí)誤差約為2 dB,但隨著轉(zhuǎn)速的升高誤差增大,這主要是由于隨著轉(zhuǎn)速的升高,更高階的徑向模態(tài)被激發(fā)。

        圖12 總聲壓級(jí)沿軸向?qū)Ρ?/p>

        n=0.8和n=1.0各階模態(tài)的聲壓級(jí)沿軸向?qū)Ρ热鐖D13所示。從圖中可見(jiàn),當(dāng)n=0.8時(shí),只有第1階模態(tài)處于傳播狀態(tài);而當(dāng)n=1.0時(shí),前3階模態(tài)均被激發(fā)。正如McAlpine等[7-8]所指出的,高階徑向模態(tài)被激發(fā)時(shí)將影響壁面聲壓的測(cè)量,這在一定程度上也增大了計(jì)算結(jié)果與試驗(yàn)結(jié)果之間的誤差。

        圖13 各階模態(tài)聲壓級(jí)沿軸向?qū)Ρ?/p>

        為了進(jìn)一步考察在同一轉(zhuǎn)速線(xiàn)上不同工作狀態(tài)下的聲學(xué)特性,給出了轉(zhuǎn)速n=0.957時(shí),不同狀態(tài)下風(fēng)扇轉(zhuǎn)子前無(wú)量綱脈動(dòng)壓力,如圖14所示;各階模態(tài)聲壓級(jí)沿軸向的對(duì)比如圖15所示。結(jié)合圖14、15可見(jiàn),在同一轉(zhuǎn)速線(xiàn)上,隨著風(fēng)扇工作狀態(tài)的提高,脈動(dòng)壓力的分布逐漸由第1、2階模態(tài)占主導(dǎo)變?yōu)榈?階徑向模態(tài)占主導(dǎo),此時(shí)計(jì)算的聲壓級(jí)與試驗(yàn)值誤差變小。

        圖14 在n=0.957不同狀態(tài)下風(fēng)扇轉(zhuǎn)子前無(wú)量綱脈動(dòng)壓力

        圖15 在n=0.957各階模態(tài)聲壓級(jí)沿軸向?qū)Ρ?/p>

        需指出的是,跨聲速風(fēng)扇的噪聲源主要包含激波噪聲、葉片定常載荷噪聲、轉(zhuǎn)/靜干涉噪聲和寬頻噪聲等,當(dāng)轉(zhuǎn)子葉尖馬赫數(shù)超過(guò)1.2時(shí),與激波相關(guān)的噪聲將占主導(dǎo)地位,因此對(duì)激波相關(guān)噪聲的預(yù)測(cè)和控制將是風(fēng)扇設(shè)計(jì)過(guò)程中的主要關(guān)注點(diǎn),盡管激波的產(chǎn)生是強(qiáng)烈的非線(xiàn)性過(guò)程,而波分解方法的前提條件是基于線(xiàn)性的小擾動(dòng)假設(shè),因此,不可避免地存在誤差,但從本文的結(jié)果來(lái)看,使用波分解方法分析風(fēng)扇前的單音噪聲水平,在一定程度上可以滿(mǎn)足工程使用需求,特別是在風(fēng)扇氣動(dòng)設(shè)計(jì)階段,可以采用該方法在對(duì)設(shè)計(jì)方案進(jìn)行性能評(píng)估的同時(shí)開(kāi)展相關(guān)噪聲的評(píng)估,指導(dǎo)風(fēng)扇的低噪聲設(shè)計(jì)。實(shí)際上,本文在前期提出了針對(duì)轉(zhuǎn)/靜干涉噪聲的氣動(dòng)/聲學(xué)一體化設(shè)計(jì)方法[17-18],該方法以聲學(xué)快速評(píng)估方法和詳細(xì)評(píng)估方法為基礎(chǔ),分層次保證了風(fēng)扇OGV的低噪聲設(shè)計(jì),遵循同樣的技術(shù)路線(xiàn),本文作為詳細(xì)評(píng)估方法,結(jié)合正在發(fā)展的激波相關(guān)噪聲快速評(píng)估方法,未來(lái)將以此為基礎(chǔ)開(kāi)展針對(duì)風(fēng)扇轉(zhuǎn)子的氣動(dòng)/聲學(xué)一體化設(shè)計(jì)。

        4 結(jié)論

        (1)針對(duì)跨聲速風(fēng)扇,前傳噪聲與來(lái)流相對(duì)馬赫數(shù)并不成正比關(guān)系,其在葉尖馬赫數(shù)處于1.1~1.2時(shí)達(dá)到最大值,之后隨著葉尖馬赫數(shù)進(jìn)一步增大反而減小;

        (2)隨著風(fēng)扇工作狀態(tài)的提升,前緣的激波越斜,與軸向夾角增大,激波強(qiáng)度降低,其引起的風(fēng)扇前傳噪聲也相應(yīng)降低;

        (3)計(jì)算結(jié)果與試驗(yàn)結(jié)果的對(duì)比分析表明,采用CFD仿真結(jié)合波分解的方法能夠?qū)げㄏ嚓P(guān)噪聲進(jìn)行有效的量化評(píng)估,當(dāng)只有低階徑向模態(tài)被激發(fā)時(shí),最小誤差能達(dá)到2 dB。

        猜你喜歡
        模態(tài)
        基于BERT-VGG16的多模態(tài)情感分析模型
        跨模態(tài)通信理論及關(guān)鍵技術(shù)初探
        一種新的基于模態(tài)信息的梁結(jié)構(gòu)損傷識(shí)別方法
        多跨彈性支撐Timoshenko梁的模態(tài)分析
        車(chē)輛CAE分析中自由模態(tài)和約束模態(tài)的應(yīng)用與對(duì)比
        國(guó)內(nèi)多模態(tài)教學(xué)研究回顧與展望
        基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識(shí)別
        由單個(gè)模態(tài)構(gòu)造對(duì)稱(chēng)簡(jiǎn)支梁的抗彎剛度
        利用源強(qiáng)聲輻射模態(tài)識(shí)別噪聲源
        日版《午夜兇鈴》多模態(tài)隱喻的認(rèn)知研究
        電影新作(2014年1期)2014-02-27 09:07:36
        亚洲色国产欧美日韩| 国产伦精品一区二区三区| 日韩精品视频在线观看无| 中文字幕乱码高清完整版| 醉酒后少妇被疯狂内射视频| 日韩欧美在线观看成人| 日本高清一区在线你懂得| 四虎永久在线精品免费一区二区| 欧美另类高清zo欧美| 曰韩精品无码一区二区三区| 国产极品嫩模大尺度在线播放| 中出人妻希奇杰卡西av| 性色av闺蜜一区二区三区| 久久这里只有精品9| 久久综合久中文字幕青草| 午夜天堂一区人妻| 国产精品igao视频网| 全部免费国产潢色一级| 少妇人妻一区二区三飞| 亚洲一区二区女搞男| 国产av电影区二区三区曰曰骚网| 国产成人av在线影院无毒| 国产精品女同一区二区免| 亚洲va韩国va欧美va| 亚洲精品自产拍在线观看| 国产日韩亚洲中文字幕| 国内自拍速发福利免费在线观看| 国产午夜鲁丝片av无码| 亚洲人成精品久久久久| 亚洲av免费看一区二区三区| 天天躁夜夜躁av天天爽| 国产精品人妻一区夜夜爱| 午夜无码亚| 午夜精品免费视频一区二区三区| 白丝兔女郎m开腿sm调教室| 国产成人一区二区三中文| 日韩亚洲国产中文字幕| 乱老年女人伦免费视频| 国产av无码专区亚洲av手机麻豆| 国产精品三级av一区二区| 蜜桃视频网站在线观看一区|