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

風雨場條件下汽車乘員艙氣動噪聲數值模擬

2024-01-11 17:08:09宗軼琦陶海沈輝楊易羅澤敏
時代汽車 2023年24期
關鍵詞:模型

宗軼琦 陶海 沈輝 楊易 羅澤敏

摘 要:本研究以某汽車為研究對象,基于數值模擬探討不同降雨量工況下的汽車乘員艙氣動噪聲聲壓級水平。采用Realizable /LES湍流模型來模擬無雨時的單相流流場,再添加離散相模型(DPM)來模擬有雨時的兩相流流場,以聲類比(FW-H)方法獲得了不同降雨量下車身表面各子系統的1/3倍頻程平均輸入激勵,采用混合有限元-統計能量分析(FE-SEA)方法獲得了駕駛員耳旁氣動噪聲聲壓級水平。仿真結果表明:在20-1000Hz頻段內,有雨工況下的駕駛員耳旁的聲壓級在各中心頻率處都低于無雨工況,小雨工況和中雨工況下的駕駛員耳旁總聲壓級較為接近,大雨工況下的駕駛員耳旁總聲壓級最低。

關鍵詞:氣液兩相流 乘員艙氣動噪聲 Realizable /LES 離散相模型 有限元-統計能量分析法 數值模擬

汽車在低速行駛時,車內噪聲主要是發動機噪聲和路面輪胎噪聲,當汽車速度超過80km/h時,風噪占主導地位[1]。風噪是一種空氣動力性噪聲,封閉乘員室內部的氣動噪聲聲源項主要是偶極子聲源,偶極子聲源是是由車身表面湍流邊界層內的擾動、表面脈動壓力共同引起的。

如今,越來越多學者、專家致力于對風噪的研究,他們從實驗、理論分析、數值模擬這三個方面出發,在討論汽車流場、汽車風噪分析技術和降低汽車風噪方面提供了許多新思路和要點。鄒銳[2]運用CFD方法對某車型進行了外流場瞬穩態仿真,穩態上分析了外流場氣流流動狀況及氣流分離情況,機艙蓋尾渦、A柱渦、后視鏡尾渦的形成、發展以及對車內噪聲的影響,瞬態上在A柱、后視鏡和側窗玻璃上選取了若干監測點,從流場與聲場上具體分析了車外湍流對該區域的影響。宗軼琦[3]運用LES與FE-SEA方法對車內噪聲進行了研究,發現了FE-SEA模型在20-100Hz能夠較為準確的捕捉車內噪聲響應峰值,但與實車道路試驗對比,計算精度略遜于FEM模型;在200-500Hz區域,FE-SEA模型相比于FEM模型、SEA模型、BEM模型,計算精度最高;在500Hz以后的高頻區域內,FE-SEA模型也能保證較高的計算精度。然而這些研究都僅限于研究汽車由于氣流分離產生的氣動噪聲,也即只考慮了由單相流工況下的氣動噪聲,沒有考慮到多相流工況下的氣動噪聲,如汽車在雨天行駛時,就屬于氣液兩相流工況,因為此時的環境變量既包括空氣,又包括雨滴。這里例舉一些其他機械在氣液兩相流工況下的響應情況。曾廣志[4]對風雨環境下橋上城際列車的運行安全性做了研究,研究結果表明:列車和橋梁迎風側表面附近的雨滴密度隨著側風風速和風向角的增加而增加,較之于無雨工況下,在有雨條件下列車的表面壓力、側向力和傾覆力矩系數有增大的趨勢。張坻[5]等對輸流管道的兩相流噪聲進行了研究,研究結果表明:由于管道中的氣泡生成與發展和兩相流產生的壓力脈動和速度脈動是兩相流噪聲產生的根本原因,低馬赫數下,偶極子聲源為主要聲源。楊顯鋒[6]使用CFX和Virtual.Lab Acoustic 模擬發動機排氣管內聲場,獲得了管內聲場在低液相體積分數下的分布規律,并搭建了發動機排氣噴淋冷卻模擬實驗臺,驗證了冷卻水的噴入對降低管內排氣噪聲的積極作用。

綜上所述,這些研究只考慮到了汽車在單相流工況下的氣動噪聲特性,而氣液兩相流中的流場與聲場特性相對于單相流是有變化的。本研究探討了汽車在氣液兩相流工況下的流場與聲場特性。本研究以30m/s行駛的某汽車為研究對象,以20-1000Hz范圍內1/3倍頻程為研究范圍,選擇了合適的湍流模型和兩相流模型來分別計算汽車在無雨工況時的單相流流場和有雨工況時的兩相流流場,構建與車型尺寸相適應的計算域,對無雨工況和有雨工況各劃分一種網格,做了網格無關性驗證以保證網格的精度。對比分析了汽車在無雨、小雨、中雨、大雨工況下的流場特性,采用聲類比的方法并在車身表面選取合適的監測點以得到車內噪聲的輸入激勵。構建了整車FE-SEA混合模型,計算出模型需要的關鍵參數,獲得駕駛員耳旁在不同降雨量下的聲壓級水平。

1 汽車流場數值計算方案

1.1 CFD數值模擬模型

1.1.1 湍流模型

如今較為常用的湍流數值模擬方法有三種:直接數值模擬(DNS)、雷諾時均模擬(RANS)、大渦數值模擬(LES)。DNS一般只適用于雷諾數較低的湍流運動,且計算量大,需要消耗較多的CPU時間和內存消耗。RANS是當今較為熟知的湍流模擬方法,其對應的湍流模型有標準模型、RNG 模型、Realizable 模型和其他湍流模型。Realizable 模型相對于前面的兩種模型精度更高,適用于旋轉流動、邊界層流動、流動分離等,即適用于汽車高速行駛時的流場。LES不但能夠精確求解某個尺度以上所有湍流的運動,捕捉RANS方法無法實現的許多非穩態,非平衡過程中出現的大尺度效應和擬序結構,而且克服了直接數值模擬計算量巨大的問題[7]。宗軼琦[3]以國際標模MIRA模型為基礎,以車身縱對稱面和側窗表面監測點靜壓系數為參考對象,評估了各種湍流模型,并與試驗對比,結果表明Realizable /LES湍流模型計算精度最高。

因此,對于無雨工況,本文選用Realizable 模型作穩態計算,以穩態計算的結果作為LES瞬態計算的初始值。

1.1.2 多相流模型

根據參與流動的項的數目,多相流可分為兩相流、三相流、四相流等,其中兩相流最為常見[8],本文研究的多相流流場包括空氣和雨滴,所以屬于氣液兩相流問題。FLUENT中的模擬多相流的模型包括歐拉-歐拉類多相流模型和歐拉-拉格朗日類多相流模型,前者連續相和離散相都采用歐拉法進行求解,后者連續相采用歐拉法,離散相采用拉格朗日法求解。

本研究流場域中的雨滴體積占有率遠小于10%。對于體積分數小于10%的氣泡、液滴、和粒子負載流動,應采用離散相模型。FLUENT中離散相模型采用的就是歐拉-拉格朗日法的計算思路。在離散相模型中,連續相介質的運動仍然由經典的N-S方程控制,離散相介質由獨立的動量方程所控制。

因此,對于有雨工況,選用Realizable 模型作穩態計算,并以穩態計算的結果作為初始值,采用LES模型與離散相模型進行風雨兩相流場的同步迭代計算。

1.2 實車模型

本文所采用的實車模型如圖1所示,該模型長5.016m,寬1.866m,高1.509m。為了提高計算效率,在保證計算精度的同時,簡化車身主體,省去車門把手及雨刷器等附件。

1.3 計算域及網格劃分

計算域設置如圖2所示,域為11倍的車長,7倍車寬,5倍車高。計算域入口距車頭3倍車長,出口距車尾7倍車長。

車身是一個復雜的幾何體,其包含眾多曲面。四面體網格適用于復雜的幾何體,因此選取四面體網格作為體網格,選取三角形網格作為面網格。

汽車高速行駛時,車身周圍的流場常伴有渦的分離與脫落,并在車身表面形成湍流邊界層,流場極其復雜,因此對車身周圍500mm范圍內的網格進行適當的加密,這可以提高湍流的計算精度。為了準確獲取邊界層內部流動的信息,車身近壁面因采用精細的六面體網格,考慮到兩相流流場計算時采用了離散相模型,應該滿足網格尺寸要大于顆粒尺寸,這里的顆粒尺寸指的是雨滴直徑。因此對無雨工況和有雨工況的邊界層設置了不同尺寸的網格,無雨工況邊界層初始高度1mm,層數4層,網格增長比例為1.2;有雨工況邊界層初始高度4mm,層數2層,網格增長比例為1.2。為了準確的獲取車身表面壓力分布狀況,對重點表面進行適當的加密,如后視鏡、A柱、前側窗等。最終無雨工況網格總數為1200萬,有雨工況網格總數為905萬。

1.4 求解器及邊界條件

對于無雨和有雨工況,如圖2所示,計算域入口設置為速度入口,出口設置為壓力出口;來流速度30m/s;在計算域的側面、頂面、底面采用滑移壁面,車身采用無滑移壁面。穩態計算時使用SIMPLE算法對速度場和壓力場進行耦合求解,動量選擇二階迎風離散格式;瞬態計算時使用PISO算法對速度場和壓力場進行耦合求解,瞬態方程選擇有界二階隱式。無雨和有雨工況瞬態設置的采樣時間都為1s,時間步長為0.0005s,最大迭代次數為20次。

對于有雨工況,離散相模型中的顆粒相邊界條件包括顆粒直徑、速度、雨滴釋放方式、質量流率。下面說明這些邊界條件的推導過程。

1.5 雨場參數

1.5.1 降雨強度

表1給出了降雨強度分類情況,其中小時雨強更能直觀的反應一個地區的實時氣候條件[9],因此本文采用小時雨強,分析在小雨、中雨、大雨工況下的汽車流場特性和駕駛員頭部氣動噪聲聲壓級水平。

1.5.2 雨滴譜分布

根據已有的測量結果,發現天然的雨滴直徑一般在0.1-6mm之內,且服從馬歇爾-帕爾默普分布(簡稱M-P譜)[10]:

式中:為直徑為的雨滴數量;,為降雨強度,單位mm/h;為濃度,取常數值8000。

采用0.5-3.5mm范圍內7種直徑的雨滴來模擬連續直徑分布的降雨,見表2。

1.5.3 雨滴釋放速度

雨滴的釋放速度包括水平速度和豎直速度,水平速度等于空氣流速,為30m/s,豎直速度為雨滴降落時的末速度[11],公式如下:

1.5.4 雨滴釋放方式

雨滴以包裹面的方式釋放。根據式(2)、(3)、(4),直徑小于1mm的雨滴的末速度遠小于直徑大于1mm的雨滴的末速度,為了使采樣時間后0.5s內雨滴、流場、汽車三者充分耦合,因此設置了兩種釋放表面位置,對直徑小于等于1mm直徑的雨滴在速度入口處釋放,直徑大于1mm的雨滴在距離底面2m高處釋放。這里例舉直徑為1mm和2mm雨滴的包裹面,如圖3所示。

1.5.5 雨滴質量流率

雨滴質量流率()可按下式計算。

式中:單位kg/s;為雨滴末速度;為釋放表面面積。

1.6 有效性驗證

在進行流場穩態計算時,做了網格無關性驗證,以氣動阻力系數為評價指標,如圖4所示。當網格數達到905萬后,氣動阻力系數為0.283,當網格數達到1200萬后,氣動阻力系數為0.275,網格數達到1360萬后,氣動阻力系數為0.277,有雨工況與無雨工況的氣動阻力系數相對于網格數為1360萬的氣動阻力系數誤差分別為2.17%、0.7%,滿足工程允許誤差值5%[12]。

2 汽車流場分析

2.1 汽車縱對稱面流場分析

為了得到汽車在無雨、小雨、中雨和大雨時外部流場的流動狀態,考慮到汽車的對稱性,選取汽車縱對稱面瞬時等值壓力云圖和速度云圖進行分析,如圖5所示,汽車行駛速度為30m/s,T=1s。

在無雨時的汽車縱對稱面壓力云圖和速度云圖中,車頭前面正壓力很大,壓力梯度變化明顯,在接近車頭的區域,氣流速度接近為零,此時壓力達到最大值。這是因為車頭處于正面迎風區域,氣流不斷沖擊車頭,速度越大,正壓力越大。氣流向車頭部上面流動時,由于車頭上圓角曲率大,產生流動分離,速度提高形成負壓區。之后縱對稱面上的氣流沿著發動機蓋流動,在發動機蓋與前風擋夾角處發生氣流分離,并形成正壓區。氣流到達風擋上邊緣時,由于結構變化,速度增加,并在車頂面形成負壓區。之后氣流沿著后檔玻璃流動,由于后備箱蓋的阻礙,在后檔與后備箱蓋上部之間形成正壓區,由于后備箱蓋結構的變化,氣流流速降低,在車尾附近又形成負壓區。

比較各降雨量下的縱對稱面等值壓力云圖可知,在有雨情況下,汽車頭部、發動機蓋、前檔、頂部區域壓力梯度大小趨勢與無雨情況相似,但在后擋風玻璃與后備箱蓋之間的區域,該區域正壓區的面積隨著降雨量的增加逐漸減少。比較各降雨量下的等值速度云圖可知,汽車縱對稱面速度梯度大小趨勢與無雨情況相似,有區別的地方是,汽車尾部“真空區”(速度接近為零的區域)的面積隨著降雨量的增加逐漸減小,真空區會在車尾端產生吸力作用,增大模型表面脈動壓力,因此會導致汽車表面的脈動壓力隨著降雨量的增加而降低。

2.2 汽車表面流場分析

圖6為T=1s時不同降雨量下的車身表面瞬態靜壓云圖,汽車有雨與無雨工況下表面的壓力梯度大小趨勢相似。在車頭、發動機蓋后部、前風擋前端與后視鏡前端存在著較大的正壓力,這與圖 5中的汽車縱對稱面正壓力區域一致。在不同降雨量下,汽車后視鏡、A柱與前側窗玻璃區域,大都處于負壓狀態,且這些負壓區的范圍基本一致,負壓是產生風噪的重要原因。

2.3 汽車表面子系統平均氣動壓力譜的計算

在進行混合FE-SEA模型計算之前,需要獲取車身表面各子系統在不同降雨量下的平均氣動壓力譜,以此作為混合FE-SEA模型的輸入激勵。這里以圖6(a)做為參考,在各子系統表面選取若干監測點。選取的原則是:在靜壓分布較為密集的地方適當的增加監測點數量,在靜壓分布較為緩和的地方適當的減少監測點數量[13]。這里以左前側窗為例,選取的監測點如圖7所示。

瞬態計算完成后取后0.5秒的數據,得到監測點的脈動壓力譜后,通過傅里葉變換,得到各監測點的三分之一倍頻程聲壓級,并將各個子系統上不同監測點的聲壓級做平均,得到各子系統的平均氣動壓力譜。左前側窗在不同降雨量下的平均氣動壓力譜如圖8所示。用同樣的方法,可以求出其余子系統在不同降雨量下的平均氣動壓力譜。

3 車內氣動噪聲數值計算

3.1 模型建立

為獲取駕駛員耳旁氣動噪聲聲壓級水平,建立整車混合FE-SEA模型,如圖9所示。該模型忽略了不影響數值仿真結果的后視鏡、輪胎、門把手等部位。

依據模態相似原則,將整車混合模型劃分為FE子系統114個,SEA平板子系統206個,其中FE子系統主要由汽車A、B、C柱、H柱、前后門門檻梁、后輪弧、橫梁組成,如圖10所示,SEA子系統主要由車門、側窗、前后擋風玻璃、底板、防火墻、發動機艙蓋、行李箱蓋、防火墻、頂棚、車燈、儀表臺、中控、座椅組成,半車SEA模型如圖11所示。

整車各子系統的物理屬性見表3。

整車聲腔子系統共分為36個,半車則為18個,如圖12。對于駕駛員室,從上至下分別為駕駛員頭部聲腔,腰部聲腔,腿部聲腔,如圖13。建立好FE-SEA模型和聲腔子系統后,創建各子系統的點、線、面連接,檢查線連接是否斷開,面連接是否正常顯示,確保各子系統之間能夠實現正常的能量流動。

3.2 建模中的關鍵參數

混合FE-SEA模型三個基本的參數為模態密度、內損耗因子和耦合損耗因子[14]。

3.2.1 模態密度的計算

模態密度是反映子系統在某一頻段內模態數密集度的一個物理量,它表征子系統從外界接收能量并引發振動的一種能力。模態密度越高,SEA方法就越能發揮其優勢。

可以通過試驗或者理論計法來計算子系統的模態密度,但由于試驗條件的限制,對于汽車模型所有的子系統,將其簡化為幾何形狀規則,厚度均勻的二維平板,其模態密度計算公式如下:

式中:為二維平板的縱向波數;為彈性模量;為材料的密度;為平板的表面積;為平板厚度;為泊松比。

左前側窗、前擋風玻璃和左前門板的模態密度如圖14所示。

聲腔子系統的模態密度可表示為:

式中:為聲腔的體積;為聲腔的表面積為聲腔的棱邊長度;為聲速。由公式可見,聲腔子系統的模態密度是頻率的函數,其主要由聲腔體積、表面積及棱邊邊長決定,其受邊界條件、阻尼和吸聲影響不大。

圖15為駕駛員頭部聲腔模態密度曲線。

3.2.2 內損耗因子

內損耗因子指系統在單位頻率、單位時間損耗的能量與平均存儲的能量之比,其計算公式如下:

式中:為結構損耗因子,其取值如表4;為邊界阻尼損耗因子,可忽略不計;為聲輻射損耗因子,其可通過下式求得:

式中:為結構的輻射比;為聲速;為空氣密度;為中心圓頻率;為臨界頻率對應的臨界波長;為平板的周長;為輻射面積;為臨界頻率;為平板的邊界條件系數,簡支邊取1,固支邊取2,一般邊界條件取。

左前側窗、前擋風玻璃和左前門板的內損耗因子如圖16所示。

聲腔內損耗因子是通過試驗測量聲場混響時間計算出來,公式為:

由于試驗條件限制,引用陳鑫[15]的數據,取平均吸聲系數為0.009,計算得到,從而繪出駕駛員頭部聲腔內損耗因子,如圖17所示。

3.2.3 耦合損耗因子

耦合損耗因子大小反映了子系統之間耦合能力的強弱,可通過試驗或者理論推導的方法獲取,也可借助VAONE軟件,因為該軟件采用了先進的波傳遞理論,將各個子系統自動連接后,便可求出子系統之間的耦合損耗因子,這大大減少了計算量,并且具有很高的精度,本研究采取的就是這種方法,計算出駕駛員頭部聲腔與左前側窗之間的耦合損耗因子,如圖18所示。

3.3 結果分析

將不同降雨量下的車身表面各子系統的平均氣動壓力譜分別施加到車身混合FE-SEA模型上,并輸入各子系統和聲腔的模態密度、內損耗因子和耦合損耗因子,數值計算出不同降雨量下且車速為30m/s時的駕駛員耳旁(即駕駛員頭部聲腔)的聲壓級(以下簡稱聲壓級),如圖19所示。

隨著頻率的增加,聲壓級水平整體都呈逐漸降低的趨勢,整體降低幅度在50%左右,聲波能量集中在20-400Hz范圍內。有雨工況下的聲壓級在各中心頻率處都低于無雨工況;小雨和中雨工況下的聲壓級變化趨勢較為吻合;大雨工況下的聲壓級在絕大多數中心頻率處都低于無雨、小雨、中雨工況。

4 結論

(1)采用Realizable /LES湍流模型來模擬無雨時的單相流流場,添加離散相模型(DPM)來模擬有雨時的兩相流流場,做了網格無關性驗證以保證汽車流場分析精度,從仿真結果中可以清楚看到整車在不同降雨量下的流場結構。

(2)在全頻段內,有雨工況下的駕駛員耳旁聲壓級在各中心頻率處都低于無雨工況,小雨和中雨工況下的駕駛員耳旁總聲壓級較為一致,大雨工況下的駕駛員耳旁總聲壓級最低。

項目基金:國家自然科學基金(51875186,51975197)。

參考文獻:

[1]龐劍,諶剛,何華.汽車噪聲與振動——理論與應用[M].北京:北京理工大學出版社,2006:353.

[2]鄒銳.乘用車高速工況車內風噪聲預測研究[D].河北工業大學,2021.

[3]宗軼琦.基于LES-FE-SEA與ALE方法的汽車氣動噪聲品質研究[D].湖南大學,2018.

[4]曾廣志.風雨環境對橋上城際列車運行安全性影響研究[D].五邑大學,2021.

[5]張坻,李孔清,王嘉,洪娜.氣液兩相流噪聲數值模擬[J].礦業工程研究,2017,32(01):71-78.

[6]楊顯鋒.管內氣液兩相流場和聲場研究[D].哈爾濱工程大學,2013.

[7]楊曉濤.汽車乘員艙氣動噪聲研究與控制[D].湖南大學,2013.

[8]車得福,李會熊.多相流及其應用[M].西安.西安交通大學出版社,2007:1.

[9]余文林,柯世堂.風雨耦合下大型冷卻塔流場特性與表面氣動力[J].南京航空航天大學學報,2020,52(04):666-674.

[10]高乾豐,董輝,鄧宗偉,朱志祥,彭文春.大型風力機風雨結構三場耦合分析[J].中南大學學報(自然科學版),2016,47(03):1011-1016.

[11]付興,林友新,李宏男.風雨共同作用下高壓輸電塔的風洞試驗及反應分析[J].工程力學,2014,31(01):72-78.

[12]張甫仁,夏文艷.邊界層網格參數對汽車外流場模擬結果的影響[J].重慶交通大學學報(自然科學版),2020,39(04):110-115.

[13]宗軼琦,張乾坤,楊易,江財茂,羅澤敏.基于LBM-FE-SEA方法汽車風窗噪聲數值模擬研究[J].噪聲與振動控制,2022,42(04):184-189+207.

[14]陳剛. 基于混合FE-SEA方法的轎車車內中頻噪聲分析與優化[D].吉林大學,2016.

[15]陳鑫. 基于SEA方法的轎車車內噪聲分析與控制研究[D].吉林大學,2008.

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产一区免费在线观看| 日韩精品一区二区三区视频免费看| 国产精品专区第1页| 九九久久99精品| 精品福利国产| 色婷婷在线影院| 国产精品成人观看视频国产| 亚洲国产成人精品青青草原| 国产人成乱码视频免费观看| 人妻中文久热无码丝袜| 999国产精品永久免费视频精品久久| 天堂岛国av无码免费无禁网站 | 亚洲欧美日韩另类在线一| 青草午夜精品视频在线观看| 91福利免费视频| 啦啦啦网站在线观看a毛片| 91九色视频网| 国产精品无码久久久久久| 国产aaaaa一级毛片| 在线播放91| 丁香综合在线| 国产视频自拍一区| 999国内精品久久免费视频| 4虎影视国产在线观看精品| 丰满人妻久久中文字幕| 免费精品一区二区h| 九色在线观看视频| 国产专区综合另类日韩一区| 国产在线拍偷自揄拍精品| 自拍偷拍欧美| 亚洲综合片| 亚洲成年人片| a亚洲视频| 国产亚洲精品97在线观看| 成人福利在线看| 免费看a级毛片| 婷婷激情五月网| 免费av一区二区三区在线| 日本不卡在线视频| 毛片大全免费观看| 久草性视频| 露脸真实国语乱在线观看| 午夜不卡视频| 亚洲国产成人自拍| 91区国产福利在线观看午夜| 美女免费黄网站| 青草精品视频| 日韩无码精品人妻| 精品久久久久无码| 在线观看免费黄色网址| 精品少妇人妻无码久久| 国产一区二区三区精品久久呦| 国产真实乱人视频| 一级看片免费视频| 伊人中文网| 国产成人乱无码视频| 午夜福利视频一区| 亚洲欧美另类久久久精品播放的| 国产日韩欧美黄色片免费观看| 伊人色在线视频| 全部免费毛片免费播放| m男亚洲一区中文字幕| 国产精品无码一二三视频| 日本精品视频一区二区| 伊人色综合久久天天| 国产精品不卡片视频免费观看| 国产在线自揄拍揄视频网站| 国产黄视频网站| 沈阳少妇高潮在线| 91精品啪在线观看国产91九色| 欧美日韩国产高清一区二区三区| 国产呦视频免费视频在线观看| 日韩无码一二三区| 波多野结衣中文字幕一区| 成人免费视频一区| 午夜成人在线视频| 色婷婷视频在线| 成人日韩视频| 波多野结衣爽到高潮漏水大喷| 色网站在线视频| 中文字幕亚洲专区第19页| 国产第二十一页|