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

        ?

        二維增升裝置前緣縫翼的遠(yuǎn)場噪聲分析

        2012-11-09 00:49:52劉志仁王福新宋文濱李亞林
        空氣動力學(xué)學(xué)報 2012年3期
        關(guān)鍵詞:聲壓級圓柱氣動

        劉志仁,王福新,宋文濱,李亞林

        (1.上海交通大學(xué)航空航天學(xué)院,上海200240;2.中國商用飛機(jī)有限責(zé)任公司上海飛機(jī)設(shè)計研究院,上海200235)

        0 引 言

        隨著社會的發(fā)展和工業(yè)技術(shù)的進(jìn)步,國際社會對民用航空業(yè)的環(huán)保要求越來越苛刻,如何進(jìn)一步降低飛機(jī)的噪聲是民用航空業(yè)目前面臨的一個重要問題。飛機(jī)的噪聲主要包括兩大類:發(fā)動機(jī)噪聲和機(jī)體噪聲。在航空業(yè)發(fā)展的初期發(fā)動機(jī)噪聲強(qiáng)度遠(yuǎn)高于機(jī)體噪聲,但經(jīng)過幾十年的發(fā)展,發(fā)動機(jī)的噴流噪聲已經(jīng)得到了顯著的降低。這使得飛機(jī)的機(jī)體噪聲達(dá)到和發(fā)動機(jī)噪聲同等的水平,航空界也越來越重視機(jī)體噪聲的研究[1]。

        已有的研究表明,在飛機(jī)的進(jìn)場和著陸過程中增升裝置和起落架是兩個主要的噪聲源。其中,前緣縫翼是增升裝置在飛機(jī)起飛和降落過程中最主要的低中頻率噪聲源。對于一些較新的客機(jī),比如波音777等,前緣縫翼的噪聲甚至超過了機(jī)翼其它部分產(chǎn)生的噪聲[2-4]。因此,研究前緣縫翼的噪聲機(jī)理以及減噪措施對降低機(jī)體噪聲具有重要的意義。目前,國外一些學(xué)者已經(jīng)通過風(fēng)洞實驗和數(shù)值計算等手段對縫翼噪聲的機(jī)理、特性進(jìn)行了研究,研究結(jié)果表明前緣縫翼噪聲主要來源于縫翼空腔處的紊流及尾緣脫落渦[5-7]。Ma Z和 Smith M 等[8]通過數(shù)值計算模擬了安裝在縫翼縫道上的聲襯的減噪效果,并對聲襯的位置和尺寸進(jìn)行了優(yōu)化。在國內(nèi),對飛機(jī)氣動噪聲的研究還處于起步階段,研究手段也比較簡單,但發(fā)展態(tài)勢良好。游亞飛等[9]采用基于 CFD(Computational Fluid Dynamics,計算流體力學(xué))和FW-H(Ffcows Williams-Hall)理論的方法研究了干凈機(jī)翼后緣的氣動噪聲特性,他們認(rèn)為翼型形狀、后緣厚度等對噪聲強(qiáng)度存在影響。

        另外,增升裝置在飛機(jī)起飛著陸階段發(fā)揮著重要作用,前緣縫翼的縫道參數(shù)是增升裝置的主要設(shè)計參數(shù),對增升裝置的最大升力系數(shù)有著很大的影響。因此,本文從工業(yè)實際應(yīng)用出發(fā),通過研究前緣縫翼縫道參數(shù)對噪聲的影響來分析增升裝置的氣動噪聲和氣動性能之間的關(guān)系,為工業(yè)設(shè)計提供參考。Emunds和Fischer[10]采用在平均流場中注入渦的方法分析了帶前緣縫翼的二維二段型縫道參數(shù)變化對縫翼噪聲的輻射方向、頻譜的影響,他們認(rèn)為盡早地在設(shè)計過程中將噪聲納入設(shè)計參數(shù),能夠更好地平衡降噪和氣動性能的需求。本文以通用計算流體力學(xué)軟件FLUENT(v12.1)為平臺,先由LES模型計算得到噪聲源,然后采用FW-H積分方法獲取遠(yuǎn)場噪聲總聲壓級,分析了縫道參數(shù)對遠(yuǎn)場噪聲的影響,并使用響應(yīng)面對原翼型的縫道參數(shù)進(jìn)行了優(yōu)化。

        1 縫翼噪聲計算原理簡介

        1.1 前緣縫翼噪聲源

        目前,國外一些研究機(jī)構(gòu)通過風(fēng)洞試驗,測量縫翼產(chǎn)生的流場和噪聲頻譜等信息,提出了縫翼噪聲產(chǎn)生的機(jī)理。Takeda和Zhang等[5]通過風(fēng)洞實驗觀測了縫翼區(qū)域的非定常流場和該區(qū)域產(chǎn)生的噪聲特性,并成功地捕捉到了縫翼尾緣渦脫落的現(xiàn)象,他們的實驗表明,尾緣渦脫落是縫翼尾部噪聲的主要來源。這一結(jié)論與Khorrami等[7]進(jìn)行的數(shù)值計算結(jié)果是一致的。Takeda等[6]人針對縫翼空腔的實驗表明空腔表面與剪切層之間的相互作用是低頻噪聲的來源。一般認(rèn)為,縫翼噪聲源分為低中頻和高頻兩個部分:高頻的窄頻噪聲是由氣流經(jīng)過縫翼與主翼之間的縫道造成的嘯叫引起的,而縫翼尾緣渦的脫落是嘯叫最主要的原因;低中頻率的寬頻噪聲是由縫翼空腔內(nèi)各類湍流流動與壁面之間的相互作用造成的。

        1.2 氣動聲學(xué)數(shù)值計算方法

        經(jīng)過近二十年的發(fā)展,計算氣動聲學(xué)已經(jīng)取得了很大的進(jìn)步,對航空氣動噪聲的預(yù)測也有了多種方法。計算氣動聲學(xué)方法大致可以分為兩類。一類是直接求解法,通過解NS方程直接求取計算域內(nèi)的聲場,這類方法需要的計算成本相當(dāng)大,以目前的計算機(jī)能力來說,還無法進(jìn)行工程應(yīng)用。另一類是混合求解法,先用非定常RANS、LES/DES等方法計算聲源或者采用RANS加入擾動因子模擬聲源,然后再用LEE方法或者FW-H等積分方法計算聲波的傳播或者輻射[11]。本文采用混合方法,先用RANS求得流場的穩(wěn)態(tài)解,然后以此為初始值用LES方法求解流場的瞬態(tài)解,得到噪聲聲源,最后使用FW-H方程進(jìn)行積分得到遠(yuǎn)場觀測點(diǎn)的聲壓隨時間的波動。

        2 二維增升裝置縫翼遠(yuǎn)場噪聲計算

        2.1 流場求解器及FW-H聲學(xué)模型

        本文采用的FLUENT商業(yè)軟件以雷諾平均N-S(Reynolds Averaged Navier Stokes,RANS)方程作為控制方程,使用有限體積法(Finite Volume Method,F(xiàn)VM)進(jìn)行離散,采用SIMPLEC方法進(jìn)行速度與壓力耦合求解,時間和空間離散精度均為二階。(1)采用SA(Spalart-Allmaras)湍流模型求解流場穩(wěn)態(tài)解,SA一方程模型其容錯功能好,處理復(fù)雜流動的能力強(qiáng),方程中包含對流項,擴(kuò)散項和源項,以非守恒形式建立。(2)采用LES方法來求解流場中的非定常壓力脈動,大渦模擬(LES)采用非穩(wěn)態(tài)的N-S方程直接模擬大尺度渦,但不計算小尺度渦,小渦對大渦的影響通過亞網(wǎng)格雷諾壓力模型來近似地模擬。

        在FLUENT噪聲模塊中,用于預(yù)測遠(yuǎn)場觀測點(diǎn)噪聲預(yù)報的方法是基于FW-H方程的聲類比積分方法。其基本原理是根據(jù)瞬態(tài)流場模擬計算結(jié)果,獲得積分表面脈動壓力隨時間的變化,然后根據(jù)FW-H方程計算得到觀測點(diǎn)隨時間變化的聲壓,再通過快速傅利葉變換(FFT)就可以得到觀測點(diǎn)的噪聲頻譜[12]。

        2.2 方法驗證

        由縫翼噪聲機(jī)理可知,縫翼噪聲源主要是偶極子聲源,其產(chǎn)生的遠(yuǎn)場噪聲指向性圖也明顯呈偶極子態(tài)。與縫翼噪聲類似的是,圓柱繞流的噪聲是由渦周期性脫落造成的,是一個典型的偶極子聲源[13]。因此,首先計算圓柱繞流的噪聲,以驗證基于FLUENT平臺的遠(yuǎn)場噪聲計算方法是否能有效預(yù)測偶極子噪聲。

        在圓柱繞流噪聲計算中,自由來流速度為69.2 m/s,雷諾數(shù)為89000,圓柱直徑為19.05mm,計算流域為二維矩形區(qū)域,坐標(biāo)原點(diǎn)位于圓柱圓心,計算域的上下和前邊界距原點(diǎn)距離為5倍圓柱直徑長度,后邊界為20倍圓柱直徑長度。采用全結(jié)構(gòu)網(wǎng)格,網(wǎng)格結(jié)點(diǎn)數(shù)為80492,第一層近壁網(wǎng)格高度為0.0005mm,計算的時間步長為5×10-6s。圖1顯示的是圓柱上方局部的瞬態(tài)渦量圖,從圖可以清晰地觀察到渦的脫落過程。

        圖1 圓柱繞流中渦的脫落Fig.1 Vortex shedding in flow around a circular cylinder

        卡門渦周期性的脫落,會引起圓柱表面的壓力按渦的脫落頻率進(jìn)行波動。為了驗證噪聲計算結(jié)果是否正確反應(yīng)了噪聲產(chǎn)生的機(jī)理,首先對圓柱隨時間變化的升力系數(shù)進(jìn)行快速傅利葉變換,得到其頻譜如圖2所示,圖中的峰值頻率即為渦的脫落頻率。

        圖2 圓柱繞流升力系數(shù)的功率譜密度Fig.2 Power spectral density of lift coefficient in flow around a circular cylinder

        卡門渦街脫落的斯德魯哈爾數(shù)一般認(rèn)為是0.2左右,根據(jù)實驗的雷諾數(shù)不同而有小幅度的不同。表1中列出了計算得到渦脫落頻率的斯德魯哈爾數(shù)與物理實驗值、經(jīng)驗值的對比??梢姳敬斡嬎爿^為準(zhǔn)確地預(yù)測了亞臨界雷諾數(shù)情況下的圓柱繞流渦放頻率。根據(jù)圓柱的噪聲機(jī)理可知,其噪聲頻譜峰值頻率應(yīng)該與渦脫落頻率一致。在圓柱正上方0.66675m(35倍圓柱直徑)的位置設(shè)置一個觀測點(diǎn),采用圓柱表面為積分面,使用FW-H噪聲模型積分得到該點(diǎn)的聲壓,通過快速傅利葉轉(zhuǎn)換可得到其頻率如圖3所示,其峰值頻率的斯德魯哈爾數(shù)同樣為0.195,與渦脫落頻率一致。這一結(jié)果說明了,基于FLUENT平臺的噪聲計算方法是能夠正確模擬偶極子聲源的噪聲產(chǎn)生機(jī)理。

        表1 圓柱繞流渦脫落頻率的斯特魯哈爾數(shù)比較Table 1 Comparison of strouhal number for vortex shedding in flow around a circular cylinder

        圖3 觀測點(diǎn)的噪聲功率譜密度Fig.3 Power spectral density of noise at observation point

        2.3 幾何模型及網(wǎng)格

        本文計算中采用UK National High Lift Programme中的L1T2模型[15]。幾何外形如圖4所示,翼型初始弦長為0.7635m,前緣縫翼的偏角為25°,后緣襟翼的偏角20°,前緣縫道寬度(Gap)為0.0212,縫道重疊量(Overlap)為-0.01,無量綱縫道參數(shù)采用的參考長度為翼型弦長0.7635m。該模型的縫翼尾緣尖端與真實縫翼一樣,是有一定厚度的。

        圖4 L1T2增升裝置模型的幾何形狀Fig.4 Geometry shape of L1T2high lift device model

        本文中計算采用的網(wǎng)格為混合網(wǎng)格,整個計算域中絕大部分網(wǎng)格為結(jié)構(gòu)網(wǎng)格,縫翼附近的網(wǎng)格用四邊形的非結(jié)構(gòu)網(wǎng)格代替結(jié)構(gòu)網(wǎng)格,以提高網(wǎng)格對不同前緣縫道參數(shù)的適應(yīng)性,具體網(wǎng)格如圖5和圖6所示,網(wǎng)格結(jié)點(diǎn)數(shù)為429,438,計算域為50倍翼型弦長的圓(圖中顯示的網(wǎng)格進(jìn)行了稀疏,密度為真實計算網(wǎng)格的1/4)。

        圖5 增升裝置附近的網(wǎng)格劃分Fig.5 Grid distribution near high lift device

        圖6 縫翼附近的網(wǎng)格劃分Fig.6 Grid distribution near slat

        2.4 遠(yuǎn)場噪聲計算及其結(jié)果

        對L1T2模型的噪聲計算分為以下幾步:第一步,采用SA湍流模型求解穩(wěn)態(tài)流場;第二步以流場穩(wěn)態(tài)解為初始值采用LES方法求解瞬態(tài)流場;第三步,瞬態(tài)流場計算穩(wěn)定后,加入FW-H聲學(xué)模型,聲源修正長度(Source Correlation Length)設(shè)為4m,以縫翼固壁為積分面求解遠(yuǎn)場噪聲。計算的條件為:攻角20.18°,來流速度0.197馬赫,計算域半徑50倍弦長的圓。

        1)圖7所示的是SA湍流模型計算得到的L1T2模型表面壓力系數(shù)與實驗值[15]的比較,計算得到的升力系數(shù)為3.997(實驗值為4.071),計算結(jié)果基本與風(fēng)洞實驗值相吻合。

        圖7 增升裝置的表面壓力系數(shù)分布Fig.7 Distribution of surface pressure coefficient of high lift device

        2)非定常計算的時間步長為5×10-6s,迭代步數(shù)先為12000步,此時流場結(jié)果已經(jīng)收斂,然后加入FW-H聲學(xué)模型再迭代9000步,得到觀測點(diǎn)的噪聲聲壓。圖8和圖9是計算得到的瞬態(tài)流場中渦量的分布圖,從圖中可以看出,計算成功捕捉到了縫翼尾緣和空腔內(nèi)的渦的運(yùn)動,這些渦的運(yùn)動是造成縫翼噪聲的主要聲源。圖10顯示的是距模型正下方50m遠(yuǎn)處觀測點(diǎn)的噪聲頻譜。

        圖8 縫翼區(qū)域的渦量云圖Fig.8 Vorticity contour near slat

        圖9 縫翼末端的脫落渦Fig.9 Vortex shedding at slat trailing edge

        3)在距模型50m遠(yuǎn)處,每隔5°取一個觀測點(diǎn),計算觀測點(diǎn)的噪聲總聲壓級,得到縫翼噪聲輻射的指向性圖如圖11所示。從圖中可以看出縫翼噪聲的輻射形態(tài)是一個典型的偶極子聲源的樣式,這個結(jié)果與縫翼噪聲原理是相吻合的。

        圖10 1/3倍頻程的聲壓級頻譜Fig.10 Acoustic spectrum in 1/3-octave band

        圖11 縫翼遠(yuǎn)場噪聲的指向圖Fig.11 Far-field directivity pattern of slat

        3 二維增升裝置的縫翼噪聲特性研究

        為了研究前緣縫道參數(shù)對前緣縫翼遠(yuǎn)場噪聲的影響,首先以L1T2模型為基準(zhǔn)模型使用不同的縫道參數(shù)構(gòu)造了24個不同的二維增升裝置模型。這些模型的縫道參數(shù)分別由5個無量綱縫道寬度值(0.01,0.02,0.03,0.04,0.05)和6個無量綱重疊量(-0.04,-0.03,-0.02,-0.01,0,0.01)組合而成(去除6組不合理的縫道參數(shù))。按照前面計算L1T2模型縫翼噪聲的方法,分別計算這24個幾何模型的縫翼噪聲,噪聲的觀測點(diǎn)統(tǒng)一設(shè)在模型垂直下方的50m處。計算的條件為來流攻角12°,速度0.2馬赫。以計算得到的24個總聲壓級數(shù)據(jù)為基礎(chǔ),使用優(yōu)化設(shè)計工具OPSYS建立縫翼噪聲總聲壓級關(guān)于縫道參數(shù)的響應(yīng)面如圖12所示。從圖中可以看出,當(dāng)縫道寬度和縫道重疊量同時增大時,觀測點(diǎn)的總聲壓級先減小再增大,縫道寬度和重疊量都很小或者很大時噪聲總聲壓級顯著增大;當(dāng)縫道寬度減小重疊量增加時,噪聲總聲壓級持續(xù)減小。

        當(dāng)然,只分析縫道參數(shù)對噪聲的影響是不夠的,在設(shè)計過程中增升裝置首先要滿足氣動性能的要求,因此還應(yīng)該分析噪聲與氣動性能之間的關(guān)系。為了分析多段翼型氣動性能與噪聲之間的相關(guān)性,同時計算了不同縫道參數(shù)對應(yīng)的最大升力系數(shù),圖13中是對應(yīng)縫道參數(shù)組合的最大升力系數(shù)云圖。將圖12和圖13比較可知,最大升力系數(shù)與總聲壓級的最優(yōu)區(qū)域并不完全重疊,但是在一定的設(shè)計范圍內(nèi)兩者還是能夠得到較好的平衡。將24個不同增升裝置的最大升力系數(shù)與噪聲總聲壓級通過縫道參數(shù)對應(yīng)起來,得到圖14,從圖中也可以看出好的設(shè)計是能夠很好地滿足氣動性能和噪聲等級兩方面要求的。

        圖12 縫翼噪聲的OASPL云圖Fig.12 OASPL contour of slat noise

        圖13 最大升力系數(shù)云圖Fig.13 Maximum lift coefficient contour

        圖14 總聲壓級與最大升力系數(shù)的關(guān)系Fig.14 Relationships between OASPL and CLmax

        4 二維增升裝置噪聲性能優(yōu)化

        采用兩種優(yōu)化方案,使用OPSYS工具[16]以降低縫翼噪聲為目標(biāo)對L1T2模型的縫道參數(shù)進(jìn)行優(yōu)化。優(yōu)化方案一:設(shè)定增升裝置的最大升力系數(shù)必須不小于4.0,才能滿足氣動性能要求,得到的有效設(shè)計參數(shù)范圍如圖15所示。優(yōu)化方案二:設(shè)定增升裝置的最大升力系數(shù)必須不小于原模型的最大升力系數(shù)。優(yōu)化的結(jié)果如表2所示,方案一和方案二得到的噪聲總聲壓級分別較原翼型下降1.07%和0.09%,可以看出原翼型的噪聲輻射強(qiáng)度本來就比較低。

        圖15 縫翼噪聲的優(yōu)化Fig.15 Optimization of slat noise

        表2 優(yōu)化結(jié)果比較Table 2 Comparison of optimized results

        5 結(jié) 論

        本文以計算流體力學(xué)軟件FLUENT為平臺,以L1T2為基準(zhǔn)模型,采用LES模型計算了增升裝置瞬態(tài)流場,并利用FW-H方法計算了遠(yuǎn)場觀測點(diǎn)的噪聲輻射。分析了縫翼縫道參數(shù)對噪聲的影響以及縫翼噪聲與氣動性能之間的關(guān)系,并對L1T2模型的縫道參數(shù)進(jìn)行了優(yōu)化。得到如下結(jié)論:

        (1)縫道寬度和重疊量大小適中時,噪聲總聲壓級較??;縫道寬度和重疊量都很小或者很大時噪聲總聲壓級顯著增大;當(dāng)縫道寬度減小重疊量增加時,噪聲總聲壓級持續(xù)減小。

        (2)最大升力系數(shù)和噪聲總聲壓級的最優(yōu)參數(shù)值范圍并不完全重合,但好的設(shè)計參數(shù)是能夠使多段翼型同時在氣動性能和噪聲強(qiáng)度兩方面得到優(yōu)化。

        (3)L1T2模型是一個具有較好噪聲特性的增升裝置模型,但通過優(yōu)化仍然能夠得到少量的改進(jìn)。

        致謝:本研究工作得到中國商用飛機(jī)有限責(zé)任公司上海飛機(jī)設(shè)計研究院提供專項資助。

        [1]GEOFFREY THOMAS,GUY NORRIS,SMITH C F,et al.Plane simple truth[M].U.S.:Aerospace Techinical Publications International Oty Ltd,2008:1-208.

        [2]DOBRZYNSKI W.Almost 40years of airframe noise research-what did we achieve[J].Journalofaircraft,2010,47(2):353-367.

        [3]SMITH M,CHOW L.Aerodynamic noise sources on high lift slats and flaps[R].AIAA Paper,2003-3226.

        [4]GUO Y,YAMAMOTO K,STOKER R.Componentbased empirical model for high-lift system noise prediction[J].JournalofAircraft,2003,40(5):914-922.

        [5]TAKEDA K,ZHANG X,NELSON P.Unsteady aerodynamics and aeroacoustics of a high-lift device configuration[A].AIAA.40th Aerospace Sciences Meeting &Exhibit[C].Reno,NV,2002:14-17.

        [6]TAKEDA K,ASHCROFT G,ZHANG X,et al.Unsteady aerodynamics of flap cove flow in a high-lift device configuration[A].AIAA 39th Aerospace Sciences Meeting and Exhibit[C].Reno,NV,2001.

        [7]KHORRAMI M R,BERKMAN M E,CHOUDHARI M.Unsteady flow computations of a slat with a blunt trailing edge[J].AIAAJournal,2000;38(11):2050-2058.

        [8]MA Z,SMITH M,RICHARDS S K,et al.Slat noise attenuation using acoustic liner[A].AIAA/CEAS 11th Aeroacoustics Conference(26th Aeroacoustics Conference)[C].2005:1-14.

        [9]游亞飛,宋文萍,郝海兵.機(jī)翼后緣噪聲預(yù)測研究[J].應(yīng)用聲學(xué),2008,27(2):140-147.(YOU Y F,SONG W P,HAO H Z B.Prediction of trailing edge noise of a clean wing[J].AppliedAcoustics,2008,27(2):140-147.)

        [10]EMUNDS R,F(xiàn)ISCHER M.Effect of slat settings(gap and overlap)on slat noise based on a test vortex injected upstream of the slat hook[A].AIAA/CEAS.12th Aeroacoustics Conference(27th AIAA Aeroacoustics Conference)[C].Cambridge,MA:American Institute of Aeronautics and Astronautics,2006.

        [11]MA Z.Numerical investigation of slat noise attenuation using acoustic liners[D].Southampton,UK:University of Southampton;2008.

        [12]ANSYS Inc.Ansys fluent 12.1user's guide[EB/OL].http://www.ansys.com/.2009.

        [13]ORSELLI R M,MENEGHINI J R,SALTARA F.Two and three-dimensional simulation of sound generated by flow around a circular cylinder[A].AIAA 30th AIAA Aeroacoustics Conference[C].2009.

        [14]REVELL J D,HAYS A P,PRYDZ A,et al.Experimental study of aerodynamic noise vs drag relationships for circular cylinders[J].AIAAJournal,1978,16:889-897.

        [15]MOIR I.Measurements on a two-dimensional aerofoil with high-lift devices[R].1994.

        [16]SONG W,KEANE A.Surrogate-based aerodynamic shape optimization of a civil aircraft engine nacelle[J].AIAA Journal,2007,45:2565-2574.

        猜你喜歡
        聲壓級圓柱氣動
        工程學(xué)和圓柱
        機(jī)器噪聲平均聲壓級計算方法差異性實證研究
        電動工具(2024年1期)2024-02-29 01:40:24
        中寰氣動執(zhí)行機(jī)構(gòu)
        圓柱的體積計算
        基于NACA0030的波紋狀翼型氣動特性探索
        一種計算消聲室聲壓級的新方法
        全新DXR mkll有源揚(yáng)聲器
        演藝科技(2019年4期)2019-03-30 03:21:46
        基于反饋線性化的RLV氣動控制一體化設(shè)計
        削法不同 體積有異
        Diodes1.9W D類音頻放大器提供高聲壓級水平并延長電池壽命
        丰满人妻一区二区三区免费视频| 日本精品国产1区2区3区| 自拍情爱视频在线观看| 国产乱码人妻一区二区三区| 国产成人精品一区二区不卡| 久久精品国产亚洲AV成人公司| 亚洲一区二区三区在线更新| 国产情侣自拍在线视频| 风韵多水的老熟妇| 国产日产精品久久久久久| 一区二区三区精彩视频在线观看 | 乱人伦人妻中文字幕不卡| 日本一区二区三区爱爱视频| 人妻少妇久久中文字幕| 曰本极品少妇videossexhd| 亚洲中文无码精品久久不卡| 精品久久一区二区三区av制服 | 在线a亚洲视频播放在线观看| 麻豆av在线免费观看精品| 欧美又大粗又爽又黄大片视频| 无码人妻精品丰满熟妇区| 91精品国产91久久久久久青草| 国产精品一品二区三区| 精品综合久久久久久888蜜芽| 亚洲 高清 成人 动漫| 国产亚洲无码1024| 亚洲精品视频中文字幕| 成人免费无码大片a毛片软件| 亚洲国产精品久久久久久网站| 亚洲综合在线一区二区三区| 亚洲成av人在线观看网址| 少妇的肉体k8经典| 国产不卡一区二区三区视频| 后入丝袜美腿在线观看| 饥渴的熟妇张开腿呻吟视频| 一区二区三区婷婷在线| 亚洲综合一区二区三区在线观看 | 亚洲国产一区二区三区精品| 特级无码毛片免费视频尤物| 无码不卡一区二区三区在线观看| 亚洲美女一区二区三区三州|