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

        ?

        汽車后空調(diào)頂棚風(fēng)管氣動噪聲數(shù)值仿真分析及優(yōu)化

        2023-02-01 06:33:16王偉江龍書成付玉樂
        振動與沖擊 2023年1期
        關(guān)鍵詞:模態(tài)

        黃 毅,王偉江,秦 望,龍書成,付玉樂

        (廣州汽車集團(tuán)股份有限公司 汽車工程研究院,廣州 511434)

        大型SUV車輛為保障乘員艙中后排溫度環(huán)境的舒適性,往往布置有單獨(dú)的后空調(diào)通風(fēng)制冷系統(tǒng)[1],受整車外造型和內(nèi)飾造型的空間限制,后空調(diào)系統(tǒng)中的頂棚通風(fēng)管呈現(xiàn)出截面積小,沿程彎曲多,風(fēng)管內(nèi)部及出風(fēng)口結(jié)構(gòu)復(fù)雜的特點(diǎn),風(fēng)管內(nèi)部局部流速較高,風(fēng)管內(nèi)流場表現(xiàn)出湍流特征,進(jìn)而在車內(nèi)感受到因風(fēng)管內(nèi)部紊流壓力脈動產(chǎn)生的氣動噪聲[2],從而影響車內(nèi)乘客和駕駛員的舒適性,甚至是安全性。因此,對后空調(diào)頂棚風(fēng)管引起的車內(nèi)氣動噪聲水平進(jìn)行前期分析評估和控制至關(guān)重要。

        國內(nèi)外學(xué)者對于空調(diào)風(fēng)管引起的氣動噪聲數(shù)值計算問題研究主要圍繞不同數(shù)值方法和湍流模型的計算精度差異比較方面展開。Ayhan等[3]采用CFD和CAA耦合方法對空調(diào)風(fēng)管進(jìn)行氣動噪聲仿真計算,并對不同湍流模型精度,結(jié)果表明采用LES大渦模擬的精度較分離渦DES和RNG k-epsilon模擬具有更小的誤差。Yves等[4]采用混合耦合方法研究了不同壁面邊界條件對帶翻板的簡易空調(diào)風(fēng)管氣動噪聲數(shù)值計算精度差異,研究結(jié)果表明壁面函數(shù)法模型較低雷諾數(shù)模型更具優(yōu)勢。Franck等[5-6]采用基于格子玻爾茲曼(Lattice Boltzmann method,LBM)的數(shù)值模擬方法計算空調(diào)風(fēng)管氣動噪聲,計算精度較高。Perot等[7-8]等利用LBM方法計算空調(diào)氣動噪聲,指出格柵和彎頭管道處的湍流壓力脈動是中高頻氣動噪聲的主要來源。李啟良等[9]通過在前空調(diào)風(fēng)管出風(fēng)口增加可穿透面的方法,同時考慮偶極子和四極子噪聲,提高了風(fēng)管氣動噪聲數(shù)值計算精度。卿宏軍等[10]通過L型簡易風(fēng)管研究比較聲類比、直接模擬和耦合仿真方法對氣動噪聲的計算精度,結(jié)果表明聲類比方法求解精度最優(yōu)。

        上述學(xué)者側(cè)重于不同氣動噪聲數(shù)值計算方法精度和噪聲源位置的研究,本文針對某大型SUV后空調(diào)頂棚風(fēng)管車內(nèi)氣動噪聲問題,先通過客觀測試和主觀評價確認(rèn)后其空調(diào)頂棚風(fēng)管高檔運(yùn)行工作時的主要問題為200~450 Hz低頻氣動噪聲,然后采用LES大渦模擬和FW-H聲類比相結(jié)合的數(shù)值仿真計算整車頂棚風(fēng)管車內(nèi)氣動噪聲,最后通過頂棚風(fēng)管CFD流場分析揭示了該低頻噪聲的產(chǎn)生機(jī)理,并仿真優(yōu)化了該低頻噪聲,為整車空調(diào)噪聲NVH開發(fā)提供了有效途徑。

        1 風(fēng)管噪聲測試及問題評價

        1.1 動力學(xué)模型

        某大型SUV車型后空調(diào)頂棚風(fēng)管和乘員艙簡化處理后的幾何模型如圖1所示,其中風(fēng)管和乘員艙通過四個帶有格柵的出風(fēng)口相通,乘員艙后輪罩布置有泄壓閥出口,后空調(diào)通過單獨(dú)的后鼓風(fēng)機(jī)轉(zhuǎn)速檔位調(diào)節(jié)為頂棚風(fēng)管入口提供不同的風(fēng)量,以保障中后排乘客的溫度和空氣環(huán)境舒適性。

        圖1 后空調(diào)頂棚風(fēng)管及乘員艙幾何模型Fig.1 Rear air conditioning ceiling duct and compartment geometry model

        采用Bruel&Kjaer人工頭和數(shù)采系統(tǒng)在半消聲室進(jìn)行整車后空調(diào)頂棚風(fēng)管噪聲測試,消聲室背景噪聲為18 dB(A),人工頭置于車內(nèi)中排左側(cè)位置如圖2所示,測試時同時采集左耳(MLL)和右耳(MLR)內(nèi)的麥克風(fēng)聲壓數(shù)據(jù)。

        測試工況為整車上電后空調(diào)系統(tǒng)單獨(dú)1到7檔工作,測試時采樣頻率fs設(shè)置為12.8 kHz,采樣時間為10 s,測試過程中車內(nèi)無人且保持車門關(guān)閉。

        1.2 測試結(jié)果分析及問題評價

        對采集的MLL和MLR兩處的時域聲壓數(shù)據(jù)進(jìn)行A計權(quán)FFT頻譜分析,頻率分辨率設(shè)置為10 Hz,由于人感受到的是左右耳的綜合作用,對MLL和MLR聲壓能量平均化處理后的A計權(quán)聲壓級頻譜及主要成分構(gòu)成如圖3所示,由圖可知在各工作檔位下,總聲壓級逐漸增大,且頻譜存在200~450 Hz頻帶峰值,該頻帶峰值能量占比隨著檔位的升高總體呈增大趨勢,其中4檔該頻帶能量占比相對較低,這是由該檔位工作時噪聲頻譜在100~200 Hz和700~1 000 Hz幅值相對更大造成,與鼓風(fēng)機(jī)和電機(jī)本體在該檔位轉(zhuǎn)速下工作噪聲特性相關(guān),最高7檔工作時頻帶聲能量占比到高達(dá)37%,且具有三個尖峰幅值特征,最大幅值達(dá)到51 dB(A),是主要?dú)鈩釉肼曍暙I(xiàn)成分。

        圖2 B&K人工頭聲壓數(shù)據(jù)采集試驗(yàn)布置Fig.2 Sound pressure data acquisition test layout

        7檔工況下車內(nèi)主觀明顯感受到“轟轟”聲,采用B&K Pulse后處理濾掉200~450 Hz噪聲回放,“轟轟”聲明顯減弱,判斷該頻段是后空調(diào)高檔工作時聲學(xué)環(huán)境變差的主要問題頻段,后文通過數(shù)值仿真方法分析該問題的產(chǎn)生機(jī)理和進(jìn)行仿真優(yōu)化。

        2 頂棚風(fēng)管氣動噪聲模型

        2.1 CFD計算域模型

        以封閉幾何模型作為整個CFD計算域在star ccm+進(jìn)行四面體網(wǎng)格劃分,出風(fēng)口Y向截面網(wǎng)格劃分結(jié)果如圖4所示,其中尺寸分布為:核心區(qū)乘員艙3~15 mm,頂棚風(fēng)管1~3 mm,出風(fēng)口和導(dǎo)流片處進(jìn)行局部加密,為0.5~1 mm;壁面邊界采用3層棱柱層網(wǎng)格,各層增長率為1.1,格柵處棱柱層厚度為0.5 mm,其余壁面邊界棱柱層厚度為1 mm,網(wǎng)格總量約為3.5千萬。

        (a) Y向截面

        2.2 LES大渦模擬控制方程

        將Naiver-Stokes方程進(jìn)行物理空間離散,得到不可壓流體的LES大渦模擬控制方程[11-12]如下

        (1)

        (2)

        (3)

        2.3 氣動噪聲數(shù)值計算方法

        CAA計算氣動聲學(xué)是基于 Lighthill的聲類比理論,經(jīng)過Curle、Ffowcs Williams 和 Hawkings擴(kuò)展成為FW-H聲類比方程[13-14]如下

        (4)

        式中:p′為流場聲壓;ni為流體沿物體表面的法向方向;vn為法向速度;c0為聲速;Tij為Lighthill應(yīng)力張量;δ(f)為Dirac delta函數(shù);H(f)為Heaviside函數(shù)。

        FW-H方程中右邊的三項(xiàng)分別對應(yīng)單極子、偶極子和四極子噪聲源,其中風(fēng)管壁面可以看成剛性,單極子噪聲幾乎為零[15];汽車后空調(diào)高檔運(yùn)行風(fēng)管內(nèi)流體仍屬于低馬赫數(shù),而四極子與偶極子噪聲強(qiáng)度之比正比于馬赫數(shù)的平方[16],四極子噪聲強(qiáng)度遠(yuǎn)小于偶極子,偶極子是主要噪聲貢獻(xiàn)。

        3 計算結(jié)果分析及方案優(yōu)化

        3.1 風(fēng)管氣動噪聲數(shù)值計算驗(yàn)證

        穩(wěn)態(tài)計算相關(guān)設(shè)置如表1所示,除進(jìn)出口外都為無滑移壁面邊界,采用Curle偶極子和Proudman四極子穩(wěn)態(tài)噪聲源模型計算噪聲源分布位置和大小,入口流量邊界為圖1所示鼓風(fēng)機(jī)出口對應(yīng)后空調(diào)7檔流量,圖1所示的整車泄壓閥為壓力出口邊界,近壁面采用壁面函數(shù)法結(jié)合低雷諾數(shù)模型的ALLy+Wall Treatment混合方式處理,即根據(jù)y+的大小自適應(yīng)選擇壁面模型[17],經(jīng)計算壁面區(qū)y+值范圍為0.2~19.4,當(dāng)y+值小于5時近壁面采用低雷諾數(shù)模型求解,而大于5時采用壁面函數(shù)法求解。

        穩(wěn)態(tài)寬頻噪聲源分布計算結(jié)果如圖5所示,由圖可知噪聲源主要集中在頂棚風(fēng)管上,其中偶極子最大噪聲源為88.0 dB(A),四極子最大噪聲源為78.9 dB(A),四極子噪聲能量僅為偶極子的2%左右,因此,偶極子噪聲源是風(fēng)管氣動噪聲的主要噪聲貢獻(xiàn)。

        (a) 偶極子噪聲源

        保持穩(wěn)態(tài)計算邊界條件,設(shè)置瞬態(tài)湍流模型和求解器參數(shù)如表2所示,以收斂后的穩(wěn)態(tài)計算結(jié)果為初始條件再進(jìn)行瞬態(tài)計算,以風(fēng)管和格柵為FW-H聲源積分面,以車內(nèi)中排左乘客右耳處為FW-H聲壓監(jiān)測點(diǎn),輸出瞬態(tài)計算穩(wěn)定后0.1~0.3 s監(jiān)測點(diǎn)時域聲壓數(shù)據(jù),并進(jìn)行A計權(quán)平均聲壓級后處理計算,聲壓級數(shù)值計算和試驗(yàn)測試對比結(jié)果如圖6所示。

        表2 瞬態(tài)計算模型及求解參數(shù)設(shè)置Tab.2 Unsteady model and solver settings

        由圖6可知,數(shù)值計算的頻譜在200 Hz以下較試驗(yàn)結(jié)果低,主要是由于數(shù)值計算未考慮鼓風(fēng)機(jī)的機(jī)械等噪聲。20~1 000 Hz和200~450 Hz頻帶計算總聲壓級精度在95%以上,由于仿真計算未考慮實(shí)車測試時存在鼓風(fēng)機(jī)和暖通箱的噪聲,這從側(cè)面說明高檔運(yùn)行時風(fēng)管是主要?dú)鈩釉肼曉床考?。風(fēng)管頻譜能夠有效捕捉到200~450 Hz的突出特征,并且分別有260 Hz、320 Hz和430 Hz三處的峰值頻率成分和試驗(yàn)結(jié)果對應(yīng),三處峰值較試驗(yàn)值誤差在10%以下,計算值偏大的原因主要是數(shù)值計算未考慮風(fēng)管內(nèi)部噪聲的折反射和風(fēng)口吸音棉對噪聲的吸收。

        (a) A計權(quán)噪聲頻譜

        3.2 低頻帶噪聲產(chǎn)生機(jī)理分析

        氣動噪聲主要是由于壓力脈動引起,對流場進(jìn)行壓力脈動FFT分析,200~450 Hz頻帶壓力脈動分布結(jié)果如圖7所示,由圖7可知該頻帶的最大壓力脈動主要集中在左側(cè)P1和右側(cè)P2點(diǎn)附近導(dǎo)流片處,以及左右側(cè)最窄處P3和P4點(diǎn)附近。

        圖7 200~450 Hz壓力脈動分布Fig.7 200-450 Hz pressure fluctuation distribution

        P1~P4處附近在0.2 s時的湍動能如圖8所示,由圖可知該四處明顯存在較大的湍動能集聚區(qū),而湍動能大小與渦的存在和流體流動分離直接相關(guān)。

        圖8 湍動能分布Fig.8 Turbulent kinetic energy distribution

        (a) P1附近Z截面

        進(jìn)一步提取該四處附近Z向截面的速度流線及壓力分布結(jié)果如圖9所示,由圖9(a)可知,左側(cè)的導(dǎo)流片P1處為了保證左側(cè)兩個出風(fēng)口的流量分配要求,在導(dǎo)流片中間增加了一個流向左后側(cè)出風(fēng)口的間隙,導(dǎo)流片被該間隙分成兩段,該間隙沖出較高流速的流體,由于慣性作用大于黏性作用,使得流體在導(dǎo)流片凹側(cè)發(fā)生分離,表現(xiàn)出湍流狀態(tài),形成局部渦,并形成明顯的負(fù)壓區(qū),這必然加劇該處的壓力脈動,進(jìn)而對外表現(xiàn)出強(qiáng)烈的氣動噪聲。同理,由圖9(b)可知,右側(cè)風(fēng)管P2位置處的導(dǎo)流片為保證右側(cè)出風(fēng)口的流量分配在Y向急劇偏轉(zhuǎn),形成倒L形結(jié)構(gòu),使得流體通過后在導(dǎo)流片凹側(cè)形成局部渦,產(chǎn)生大范圍的負(fù)壓區(qū),也必然加劇該處的壓力脈動。與此同時,P1和P2位置處的導(dǎo)流片明顯使得通過后的流體偏向風(fēng)管外側(cè)(如圖9(c),9(d)所示),使得下游P4和P3處的流體明顯向外側(cè),流體出現(xiàn)分離并形成渦,加劇壓力脈動,形成氣動噪聲源。

        圖10 P1~P4點(diǎn)壓力頻域Fig.10 Frequency domain pressure data of P1-P4

        監(jiān)測記錄瞬態(tài)計算過程中P1~P4點(diǎn)的時域壓力數(shù)據(jù)結(jié)果,并進(jìn)行FFT分析結(jié)果如圖10所示,由圖可知在200~450 Hz之間存在相對較大的壓力脈動,其中P2點(diǎn)呈現(xiàn)出260和320 Hz左右的峰值特征,其它三點(diǎn)呈現(xiàn)出寬頻的特征,由上可知由于P1~P4附近渦的存在,使得該四處存在較大的壓力脈動,最大可達(dá)到7 Pa左右,且壓力脈動的頻率在200~450 Hz相對較大,成為低頻氣動噪聲的主要激勵源。另一方面,氣動噪聲的大小也決定于風(fēng)管聲學(xué)模態(tài)的影響,因此通過有限元計算頂棚風(fēng)管單位激勵下的出口聲壓頻響和聲模態(tài),結(jié)果分別如圖11和圖12所示。

        由圖11可知,在1 kHz以內(nèi)存在多個頻響峰值,在峰值處都有對應(yīng)的聲模態(tài),尤其在200~450 Hz之間頻響明顯較大,若壓力脈動激勵頻率及位置和對應(yīng)聲模態(tài)耦合時,必然造成聲模態(tài)頻率處的噪聲峰值。

        由圖12可知,256和274 Hz聲模態(tài)的反節(jié)點(diǎn)主要位于P3和P4處,328和451 Hz聲模態(tài)在兩側(cè)導(dǎo)流片P1及P2處也有反節(jié)點(diǎn)分布,而反節(jié)點(diǎn)位置又存在和聲模態(tài)耦合的壓力脈動激勵頻率,這會在乘員艙輻射產(chǎn)生劇烈的氣動噪聲,并呈現(xiàn)處聲模態(tài)頻率附近的峰值特征,這與風(fēng)管噪聲仿真和測試結(jié)果中低頻帶中三個峰值頻率成分結(jié)果一致。

        圖11 風(fēng)管出口頻響及聲模態(tài)分布Fig.11 Frequency response of duct oulet and acoustic mode distribution

        3.3 方案優(yōu)化

        頂棚風(fēng)管受到整車內(nèi)外造型和周圍線束的布置限制,風(fēng)管外廓造型很難做大的調(diào)整,并且風(fēng)管聲模態(tài)較為密集,通過風(fēng)管造型調(diào)整聲模態(tài)較為困難。

        由于導(dǎo)流片引起的渦流壓力脈動可以通過導(dǎo)流片結(jié)構(gòu)的優(yōu)化來有效降低,實(shí)現(xiàn)流量分配和噪聲的性能平衡,因此對左右側(cè)的導(dǎo)流片優(yōu)化思路為:①左側(cè)導(dǎo)流片去掉中間間隙以避免從間隙沖出的高速流體形成的分離渦,為彌補(bǔ)間隙流向左側(cè)后出風(fēng)口的流量,將導(dǎo)流片的上游縮短以增加左側(cè)后出風(fēng)口的流量,最終成一段式導(dǎo)流片如圖13優(yōu)化方案左側(cè)所示;②右側(cè)L型導(dǎo)流片在X方向的過渡梯度減小,即讓導(dǎo)流片在Y向的彎折變緩,這樣既可以降低導(dǎo)流片在凹側(cè)的流體分離程度,同時盡量不改變原流量分配比例。對優(yōu)化導(dǎo)流片后的風(fēng)管出風(fēng)口流量分配比例和原狀態(tài)對比結(jié)果如圖14所示,流量分配比例均勻程度稍有變差,比例范圍為22.9%~27.1%,但在(25±2.5)%的工程允許容差范圍內(nèi)。

        優(yōu)化前后t=0.2 s時左右側(cè)導(dǎo)流片附近的流線和靜壓分布對比如圖15、16所示,由圖15(a)和15(b)對比可知:左側(cè)由于導(dǎo)流片之間間隙的取消,優(yōu)化方案導(dǎo)流片凹側(cè)的負(fù)壓區(qū)明顯減小,渦的劇烈程度也大大降低,同時流體流向的左下側(cè)45°偏向程度也有一定的削弱,優(yōu)化后左側(cè)P1和P4點(diǎn)0.1~0.3 s內(nèi)壓力脈動時域和200~450 Hz頻域幅值較原狀態(tài)都有明顯降低(如圖17(a)和17(d)所示);由圖16(a)和16(b)對比可知:優(yōu)化方案右側(cè)風(fēng)管導(dǎo)流片凹側(cè)的負(fù)壓區(qū)范圍縮小,較原方案更分散,可以削弱集中渦的作用,有利于壓力脈動的分散化降低,風(fēng)管右側(cè)P2和P3點(diǎn)處0.1~0.3 s的壓力脈動的時域和200~450 Hz頻域幅值也相應(yīng)的隨著降低(如圖17(b)和17(c)所示),降低了壓力脈動激勵源。

        (a) 256 Hz聲模態(tài)振型

        3.4 噪聲優(yōu)化結(jié)果

        優(yōu)化方案的噪聲仿真結(jié)果如圖18所示,由圖可知:優(yōu)化方案的噪聲頻譜在200~450 Hz頻帶明顯降低,20~1 000 Hz和200~450 Hz頻帶的聲壓級分別降低3.2和4.8 dB(A),三處峰值平均降低5.5 dB(A),低頻的轟轟聲實(shí)現(xiàn)了有效的控制。

        圖13 風(fēng)管導(dǎo)流片優(yōu)化方案Fig.13 Optimization scheme of duct deflector

        圖14 風(fēng)管出風(fēng)口流量分配百分比對比Fig.14 Air duct outlet flow distribution percentage

        (a) 原狀態(tài)

        (a) 原狀態(tài)

        (a) P1點(diǎn)壓力時域及頻域

        (b) 頻帶及峰值聲壓級圖18 原狀態(tài)與優(yōu)化方案噪聲對比Fig.18 Noise comparison between original and optimized scheme

        4 結(jié) 論

        針對某大型SUV車型后空調(diào)高檔運(yùn)行下頂棚風(fēng)管的200~450 Hz低頻轟轟聲問題,在噪聲試驗(yàn)的基礎(chǔ)上,通過LES大渦模擬和FW-H聲類比相結(jié)合的數(shù)值仿真計算方法分析了該頻帶噪聲的產(chǎn)生機(jī)理,并在不影響風(fēng)量分配比例的基礎(chǔ)上實(shí)現(xiàn)了優(yōu)化控制,形成如下結(jié)論:

        (1) 后空調(diào)頂棚風(fēng)管是高檔運(yùn)行氣動噪聲的主要貢獻(xiàn)源,且以偶極子噪聲為主,四極子噪聲源強(qiáng)度只占偶極子的2%左右,噪聲仿真時可以不予考慮;

        (2) LES和FW-H相結(jié)合的方法能夠準(zhǔn)確地計算頂棚風(fēng)管車內(nèi)氣動噪聲,其中總聲壓級精度在95%以上,同時能夠準(zhǔn)確捕捉低頻帶“轟轟”聲的噪聲峰值特征,精度也能達(dá)到90%以上,該方法可用于整車空調(diào)風(fēng)管的氣動噪聲預(yù)測和開發(fā)控制;

        (3) 風(fēng)管內(nèi)導(dǎo)流片處產(chǎn)生的渦流壓力脈動是引起200~450 Hz的氣動噪聲的激勵源,當(dāng)激勵頻率與頂棚風(fēng)管的聲模態(tài)耦合,并且渦的位置和聲模態(tài)反節(jié)點(diǎn)重合時,會形成低頻帶內(nèi)對應(yīng)聲模態(tài)頻率的噪聲峰值;

        (4) 頂棚風(fēng)管導(dǎo)流片的合理設(shè)計可以有效降低渦引起的壓力脈動,進(jìn)而降低對應(yīng)激勵頻率氣動噪聲峰值,實(shí)現(xiàn)風(fēng)管出口流量分配和氣動噪聲的性能平衡。

        猜你喜歡
        模態(tài)
        基于BERT-VGG16的多模態(tài)情感分析模型
        跨模態(tài)通信理論及關(guān)鍵技術(shù)初探
        一種新的基于模態(tài)信息的梁結(jié)構(gòu)損傷識別方法
        多跨彈性支撐Timoshenko梁的模態(tài)分析
        車輛CAE分析中自由模態(tài)和約束模態(tài)的應(yīng)用與對比
        國內(nèi)多模態(tài)教學(xué)研究回顧與展望
        基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識別
        由單個模態(tài)構(gòu)造對稱簡支梁的抗彎剛度
        計算物理(2014年2期)2014-03-11 17:01:39
        利用源強(qiáng)聲輻射模態(tài)識別噪聲源
        日版《午夜兇鈴》多模態(tài)隱喻的認(rèn)知研究
        電影新作(2014年1期)2014-02-27 09:07:36
        亚洲国产人成自精在线尤物| 精品人妻伦九区久久aaa片69| 一本色道久久综合狠狠躁 | 亚洲午夜无码久久yy6080 | 韩国三级大全久久网站| 午夜精品男人天堂av| 久久精品中文字幕第一页| 亚洲精品熟女av影院| 成年丰满熟妇午夜免费视频| a级毛片成人网站免费看| 少妇对白露脸打电话系列| 国产精品中文第一字幕| 久久久诱惑一区二区三区| av网页免费在线观看| 国产乱人伦av在线a麻豆| 国产精品区一区第一页| 麻豆91免费视频| 国产精品人人爱一区二区白浆| 一区二区三区四区中文字幕av| 亚洲一区二区三区四区五区六| japanesehd中国产在线看| 国产aⅴ夜夜欢一区二区三区| 一本色道久久综合亚州精品| 国产精品会所一区二区三区| 欧美成人午夜精品久久久| 亚洲AV无码秘 蜜桃1区| 国产精品日本一区二区三区在线| 最美女人体内射精一区二区 | 亚洲成人免费网址| 亚洲高清国产拍精品熟女| 日韩精品久久中文字幕| 亚洲国产精品福利片在线观看| 亚洲AV综合A∨一区二区| 日本在线一区二区三区视频| 日韩av午夜在线观看| 国产亚洲精久久久久久无码77777 丝袜足控一区二区三区 | 亚洲精品成人网站在线播放| 污污污污污污WWW网站免费| 日本中出熟女一区二区| 欧美国产激情18| 亚洲国产成人精品无码区99|