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

高速列車頭型長細(xì)比對(duì)氣動(dòng)噪聲的影響1)

2017-11-11 01:54:16莫晃銳劉青泉
力學(xué)學(xué)報(bào) 2017年5期

安 翼 莫晃銳 劉青泉

?(中國科學(xué)院力學(xué)研究所流固耦合系統(tǒng)力學(xué)重點(diǎn)試驗(yàn)室,北京100190)

?(北京理工大學(xué)宇航學(xué)院力學(xué)系,北京100081)

高速列車頭型長細(xì)比對(duì)氣動(dòng)噪聲的影響1)

安 翼?莫晃銳?劉青泉?,2)

?(中國科學(xué)院力學(xué)研究所流固耦合系統(tǒng)力學(xué)重點(diǎn)試驗(yàn)室,北京100190)

?(北京理工大學(xué)宇航學(xué)院力學(xué)系,北京100081)

高速列車的頭尾車外形對(duì)氣動(dòng)噪聲具有重要的影響.工程實(shí)踐中隨著車速的增加,車輛頭部越來越細(xì)長,日本高速磁懸浮列車實(shí)踐中甚至出現(xiàn)了具有極端長細(xì)比的頭部形狀.本文以討論頭型長細(xì)比對(duì)列車氣動(dòng)噪聲的影響規(guī)律為出發(fā)點(diǎn),應(yīng)用非線性聲學(xué)求解器(NLAS)和FW–H聲學(xué)比擬法的混合算法,在3種運(yùn)行速度下對(duì)基于CRH380A高速列車頭型概化的4種不同頭型長細(xì)比的模型車的氣動(dòng)噪聲進(jìn)行了數(shù)值模擬.給出了不同頭型長細(xì)比列車的流場特征、氣動(dòng)阻力和氣動(dòng)噪聲.結(jié)果表明,列車的氣動(dòng)總阻力隨頭型長細(xì)比的增大而減小,且頭型長細(xì)比對(duì)列車總氣動(dòng)阻力的影響隨運(yùn)行速度的增加而增強(qiáng).而頭型長細(xì)比對(duì)氣動(dòng)噪聲的影響呈現(xiàn)出較為復(fù)雜的影響,并不存在單調(diào)的影響關(guān)系;綜合考慮氣動(dòng)阻力和氣動(dòng)噪聲,長細(xì)比最大的頭型綜合性能較優(yōu),但差異并不顯著,因此在不考慮微氣壓波等因素的條件下,簡單增加車頭長細(xì)比并不一定能帶來明顯的氣動(dòng)噪聲性能提升.

高速列車,氣動(dòng)噪聲,氣動(dòng)阻力,頭型長細(xì)比

引言

近年來,我國的高速鐵路迅速發(fā)展,已成為我國最主要的城際客運(yùn)系統(tǒng)之一,更高速的磁懸浮列車也正在研發(fā)中.隨著列車運(yùn)行速度的不斷提升,噪聲問題日顯突出,成為影響高速列車可持續(xù)發(fā)展的關(guān)鍵問題之一.高速列車的噪聲主要由機(jī)械噪聲和氣動(dòng)噪聲組成[1],氣動(dòng)聲學(xué)理論指出,氣動(dòng)噪聲的聲功率與速度的6~8次方成正比[2],而機(jī)械噪聲則與速度的低次冪相關(guān).研究表明[34],當(dāng)列車運(yùn)行速度超過300km/h時(shí),氣動(dòng)噪聲將顯著增強(qiáng),并主導(dǎo)列車的總體噪聲.高速列車的氣動(dòng)噪聲主要來自于頭尾車、轉(zhuǎn)向架、受電弓和車體[5],其中頭尾車產(chǎn)生的氣動(dòng)噪聲是其主要來源之一[6].Mellet等[7]分析了不同時(shí)速下的TGV-Duplex和ICE3高速列車的大量噪聲實(shí)測數(shù)據(jù),發(fā)現(xiàn)頭尾車噪聲占全車噪聲的比重隨著列車速度的提高而快速增長,當(dāng)運(yùn)行速度達(dá)到300km/h以上時(shí),頭尾車輻射的噪聲超過其余八節(jié)車廂輻射的噪聲,且頭車的噪聲比尾車噪聲還要顯著.由于高速列車頭尾車的幾何外形決定著周圍流動(dòng)的附著、邊界層的發(fā)展和分離,以及列車尾部的流動(dòng)分離和所產(chǎn)生的非定常尾流[8],頭尾車氣動(dòng)噪聲的產(chǎn)生與其幾何形狀密切相關(guān)[9].

Kitagawa和 Nagakura[10]分析了日本新干線高速列車的氣動(dòng)噪聲組成以及聲源位置,發(fā)現(xiàn)光滑的車體表面可以有效地減少車體上部產(chǎn)生的氣動(dòng)噪聲.Torii和Ito[11]對(duì)新干線列車噪聲源的研究發(fā)現(xiàn),對(duì)列車鼻形的改進(jìn)可以降低標(biāo)準(zhǔn)測點(diǎn)處(距離軌道中心線25m,距地面高3.5m)約2dB(A)的噪聲級(jí),同時(shí)可有效減少列車在隧道中的壓力波.Maeda等[12]和Ido等[13]通過風(fēng)洞實(shí)驗(yàn)進(jìn)行了一系列長細(xì)比下的高速列車頭型的氣動(dòng)阻力測試,發(fā)現(xiàn)列車的氣動(dòng)阻力隨著頭型長細(xì)比的增大而有效降低.喻華華[14]曾在不同來流速度條件下,對(duì)CRH380高速列車的5種備選頭型的氣動(dòng)噪聲進(jìn)行了風(fēng)洞測試,結(jié)果表明,在相同長細(xì)比條件下,當(dāng)頭型滿足流線型設(shè)計(jì)要求時(shí),其不同橫截面形狀的車鼻對(duì)列車總體氣動(dòng)噪聲的影響較為有限.王成強(qiáng)等[15]應(yīng)用基于NLAS(nonlinear acoustics solver)的CAA模擬方法對(duì)高速列車的氣動(dòng)噪聲進(jìn)行了數(shù)值模擬研究.潘忠和陸森林[16]發(fā)現(xiàn)表面聲功率級(jí)和脈動(dòng)壓力級(jí)最大值都出現(xiàn)在鼻錐、雨刷器等表面曲率變化較大的部位.

高速列車的頭車和尾車具有一致的外形設(shè)計(jì),一般為復(fù)雜的三維曲面[17],其橫截面存在明顯形狀或面積變化的區(qū)段稱為車鼻,通常由此確定了頭型最主要的幾何特性[18].車鼻外形由眾多參數(shù)決定,為頭型特征結(jié)構(gòu)的氣動(dòng)噪聲特性研究帶來了困難,實(shí)際研究中,常定義車頭鼻形部位長度與后部車身斷面等效半徑之比為頭型的長細(xì)比[19],其與車鼻橫截面形狀分布一起,成為頭型設(shè)計(jì)的重要參數(shù).實(shí)踐中,日本在 2015年試驗(yàn)的下一代磁懸浮列車L0系采用了長達(dá)15m的車鼻,而其車廂斷面僅為3.1m×2.9m,長細(xì)比高達(dá)8.8,對(duì)列車功能和使用模式的設(shè)計(jì)都產(chǎn)生了影響.而在我國高速磁懸浮列車發(fā)展中,長細(xì)比在氣動(dòng)外形設(shè)計(jì)中的位置也是一個(gè)值得思考的問題.

高速列車頭型的長細(xì)比對(duì)列車氣動(dòng)性能有著顯著的作用和影響,然而,至今關(guān)于頭型長細(xì)比對(duì)列車氣動(dòng)噪聲影響的研究還較少,仍缺乏深刻的規(guī)律性認(rèn)識(shí).為此,本文將通過數(shù)值模擬方法,以我國自行設(shè)計(jì)的CRH380A高速列車為對(duì)象,參考飛行器設(shè)計(jì)中的優(yōu)化技術(shù)[2022],針對(duì)特定的車鼻橫截面形狀函數(shù),探討其分布區(qū)間的改變,即頭型長細(xì)比的變化,對(duì)不同運(yùn)行速度下高速列車氣動(dòng)噪聲的影響,以期為合理進(jìn)行頭型降噪設(shè)計(jì)提供科學(xué)依據(jù).

1 數(shù)值模擬方法

1.1 總體方法

本文采用計(jì)算流體力學(xué)/計(jì)算氣動(dòng)聲學(xué) (hybrid CFD/CAA)的混合方法對(duì)高速列車的氣動(dòng)噪聲進(jìn)行數(shù)值模擬.將計(jì)算區(qū)域分為非線性聲源區(qū)(近場流動(dòng)區(qū))和線性聲傳播區(qū),分別采用相應(yīng)的計(jì)算方法求解.雖然大渦模擬(LES)和脫體渦模擬(DES)類方法以及FW-H方法已經(jīng)逐漸用于噪聲計(jì)算[23],但其計(jì)算量很大.考慮到本文研究對(duì)象為細(xì)長流線體,其流動(dòng)的不穩(wěn)定性較弱,本文選用由Batten等[24]提出的計(jì)算量較低的非線性聲學(xué)求解器(non-linear acoustics solver,NLAS)作為近場流動(dòng)的求解算法.即首先使用cubic k-ε RANS模型求解Navier-Stokes方程,獲得流動(dòng)的統(tǒng)計(jì)定常解,然后運(yùn)用非線性聲學(xué)求解器(NLAS)求解流動(dòng)的非定常時(shí)空演化和壓力脈動(dòng),獲得近場流場和噪聲聲源信息.而對(duì)于遠(yuǎn)場噪聲的預(yù)測則采用聲比擬法實(shí)現(xiàn),即采用FW-H方程[25],通過在控制面的積分得到遠(yuǎn)場噪聲.由于近場部分求解的是脈動(dòng)量,避免了LES等方法直接使用流動(dòng)量造成的數(shù)值誤差,同時(shí)對(duì)流項(xiàng)影響較小,NLAS對(duì)網(wǎng)格的需求也遠(yuǎn)低于LES和DES等方法.

1.2 非線性聲學(xué)求解器(NLAS)

非線性聲學(xué)求解器是由Batten等[2627]提出的一種求解統(tǒng)計(jì)定常狀態(tài)流動(dòng)中聲的產(chǎn)生與傳播的數(shù)值算法,其控制方程是從Navier-Stokes方程的擾動(dòng)推導(dǎo)而來,稱之為非線性擾動(dòng)方程(NLDE)[24],其形式為

式中,ρ為密度,u為速度,e為能量,τ為切應(yīng)力,p為壓強(qiáng),δ為delta函數(shù).

忽略密度擾動(dòng),對(duì)上述方程組取時(shí)間平均,可以消去時(shí)間演變項(xiàng)和所有線性通量項(xiàng),得到

其中,Ri中的物理量對(duì)應(yīng)于標(biāo)準(zhǔn)雷諾應(yīng)力張量和湍流的熱通量項(xiàng),通過RANS方法可以求得這些未知項(xiàng).

1.3 FW-H聲學(xué)比擬法

遠(yuǎn)場聲壓的預(yù)測基于Farassat[28]給出的FW-H方程的時(shí)域積分解

這里,Qi=(ρ∞?ρ)vi+ρui,Li=p?ni+ρui(uj?vj)nj;ρ∞和c∞分別為遠(yuǎn)場未受擾動(dòng)流體介質(zhì)的密度和聲速,ui和vi分別表示當(dāng)?shù)亓黧w速度和物體表面速度;?ni和?ri分別為物面單位法向矢量n和單位發(fā)射矢量(x?y)/r在3個(gè)方向的分量,r=|x?y|為觀測點(diǎn)與聲源之間的距離,其中x和y分別表示觀測點(diǎn)和聲源的位置矢量.符號(hào)[]ret代表在延遲時(shí)間τ=t?r/c∞下取值,其中t和τ分別為聲源發(fā)出聲波的時(shí)間和聲波到達(dá)觀測點(diǎn)的時(shí)間;Mr=vi?ri/c∞為聲源與觀測點(diǎn)方向上的馬赫數(shù).

1.4 數(shù)值方法驗(yàn)證

首先用經(jīng)典的汽車后視鏡噪聲算例做驗(yàn)證.Hold等[29]和Siegert等[30]對(duì)放置于平板上的汽車后視鏡簡化模型的氣動(dòng)噪聲進(jìn)行了風(fēng)洞實(shí)驗(yàn).如圖1所示,后視鏡簡化模型由1/2圓柱與1/4球體拼接組成,豎直放置于平板上,圓柱直徑和高度以及1/4球體的直徑均為0.2m.將圓柱的下底面圓心設(shè)為坐標(biāo)原點(diǎn),流場中設(shè)置了兩個(gè)壓力探點(diǎn),探點(diǎn)a位于半圓柱下游表面的邊緣,坐標(biāo)為(0.0,0.117,0.085),探點(diǎn)b位于下游尾流的平板表面處,坐標(biāo)為(0.3978,0.0,0.14181).本文采用與Siegert[30]實(shí)驗(yàn)相同的后視鏡模型及幾何配置,運(yùn)用非線性聲學(xué)求解器(NLAS)數(shù)值計(jì)算后視鏡簡化模型的近場流動(dòng)與壓力脈動(dòng).

圖1 放置在平板上的汽車后視鏡簡化模型及壓力探點(diǎn)分布Fig.1 The geometrical model and probe points location of the simpli fi ed wind mirror test case

采用六面體結(jié)構(gòu)網(wǎng)格離散求解空間,流動(dòng)求解域范圍為 x ∈ [?5D,15D],y ∈ [0,10D],z ∈[?5D,5D],網(wǎng)格量約 4.9×105. 來流速度 U∞=200km/h,瞬態(tài)計(jì)算的時(shí)間步長△t=2×10?5s,約為5.6×10?3T0,其中,T0=D/U0≈ 3.6× 10?3s為流場以平均流速傳播一個(gè)特征長度D所需的時(shí)間.數(shù)值求解時(shí),后視鏡模型和底部平板采用絕熱的無滑移固壁條件,其它邊界處為來流速度200km/h、參考溫度298.5K和參考?jí)毫?9530Pa的遠(yuǎn)場邊界條件.計(jì)算可解析的信號(hào)最高頻率 fmax=1/(2?t)=25000Hz.

圖2為探點(diǎn)a與探點(diǎn)b處聲壓級(jí)的數(shù)值計(jì)算結(jié)果與Siegert實(shí)驗(yàn)結(jié)果的對(duì)比.結(jié)果顯示,在探點(diǎn)a處,數(shù)值結(jié)果低估了約40~400Hz頻段內(nèi)的聲壓級(jí),但其余頻段的聲壓級(jí)與實(shí)驗(yàn)結(jié)果吻合很好;而在探點(diǎn)b處,數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果在40~2000Hz的頻段內(nèi)整體上有著較好的吻合,在60~100Hz頻段處,NLAS較為理想地預(yù)測了探點(diǎn)b處聲壓譜的峰值特征,但峰值頻率的預(yù)測結(jié)果較實(shí)驗(yàn)結(jié)果偏低,整體上數(shù)值結(jié)果和實(shí)驗(yàn)測量結(jié)果符合良好.

圖2 探點(diǎn)a與b處聲壓級(jí)的數(shù)值結(jié)果與實(shí)驗(yàn)結(jié)果的對(duì)比Fig.2 Comparison of the calculated and experimental SPL values at the probe points

其次應(yīng)用在同濟(jì)大學(xué)的氣動(dòng)--聲學(xué)風(fēng)洞中開展的1:8模型試驗(yàn)結(jié)果做驗(yàn)證[14].該風(fēng)洞在噴口速度160km/h試驗(yàn)段背景噪聲SPL(A)為61dB,截止頻率為50Hz.模型試驗(yàn)使用CRH 380A三編組模型,為突出頭型影響,去除轉(zhuǎn)向架并將其和車廂連接處填充光滑.試驗(yàn)中測點(diǎn)在距模型7.5m處平行于車長方向布置,有4組測點(diǎn)采集了噪聲數(shù)據(jù),試驗(yàn)共研究了200,230和250km/h三種風(fēng)速.本文對(duì)三種風(fēng)速條件下的噪聲頻譜和分布進(jìn)行了研究,典型工況(250km/h,中間測點(diǎn))的計(jì)算和試驗(yàn)測量結(jié)果的比較如圖3所示.圖中可以看出,模擬結(jié)果和試驗(yàn)測量結(jié)果在200Hz以上區(qū)域吻合良好,低頻部分有差異.差異可能主要由風(fēng)洞低頻背景噪聲所引起,但總體上在列車噪聲所關(guān)注的頻率范圍,本文的數(shù)值方法給出了較好的結(jié)果.總體上驗(yàn)證了所采用數(shù)值方法的有效性和準(zhǔn)確性.

圖3 250km/h條件下CRH380A縮尺模型典型測點(diǎn)噪聲頻譜數(shù)值模擬和測量結(jié)果的比較Fig.3 Comparison of the calculated and experimental SPL values for the 1:8 scaled CRH380A model in the wind of 250km/h

2 物理模型及計(jì)算條件

2.1 物理模型

基于CRH380A高速列車的基本頭型,抽象簡化出代表列車主體結(jié)構(gòu)特征的細(xì)長結(jié)構(gòu)體,作為本文研究的列車幾何模型,模型列車為頭車、單節(jié)中間車廂和尾車組成的三節(jié)車體編組.由于車體復(fù)雜部件對(duì)流場有一定影響[3132],本文去掉車廂間隙、受電弓和轉(zhuǎn)向架等非光順曲面結(jié)構(gòu)以突出頭部形狀的影響,并將CRH380A高速列車的車鼻橫截面形狀作為初始的分布函數(shù).

將列車車鼻看作是其二維橫截面沿長度方向的分布,則可使用函數(shù)S(l)表征其形狀,S(l)為距離車頭鼻尖l位置處的橫截面形狀函數(shù),設(shè)其相對(duì)應(yīng)的面積為A(l),其中,l∈[0,ln],ln為車鼻長度,其一般應(yīng)小于單節(jié)車體的總長度.當(dāng)確定車鼻橫截面形狀函數(shù)S(l)和其分布區(qū)間[0,ln]后,頭型的最主要幾何外形即可確定,其長細(xì)比λ=ln/r,其中為頭型后部斷面的等效半徑值.

為研究頭型長細(xì)比對(duì)氣動(dòng)噪聲的影響,采用了4種具有不同長細(xì)比的頭型,它們對(duì)應(yīng)的列車幾何模型如圖4.CRH380A的頭型ln為12m,長細(xì)比λ0約為6.36,A,B,C,D四種列車模型對(duì)應(yīng)的頭型長細(xì)比分別為:λA=0.75λ0,λB=1.0λ0,λC=1.25λ0,λD=1.5λ0.

圖4 4種頭型細(xì)長比對(duì)應(yīng)的列車幾何模型Fig.4 Models of four trains with di ff erent slenderness ratio

2.2 計(jì)算區(qū)域與網(wǎng)格

簡化列車模型的單節(jié)車體長Ls=26m,寬W≈3m,高 H≈3.5m,取列車下底面面心為計(jì)算域原點(diǎn),從而車體長度在x軸上的范圍為[?39m,39m].數(shù)值計(jì)算區(qū)域長度范圍為 x∈[?4.5Ls,10.5Ls],寬y∈ [?4Ls,4Ls],高z∈ [?h,4Ls],其中h=0.371m為車體底面距離地面的高度.

計(jì)算網(wǎng)格為四面體非結(jié)構(gòu)網(wǎng)格,壁面附近區(qū)域的網(wǎng)格分辨率約為,在車體的尾流區(qū)域?qū)W(wǎng)格進(jìn)行了加密,當(dāng)?shù)鼐W(wǎng)格分辨率約為.列車車體壁面邊界層采用了棱柱網(wǎng)格,按照NLAS計(jì)算原理,第一層網(wǎng)格取在對(duì)數(shù)區(qū),其高度對(duì)應(yīng)的y+≈150,模型的計(jì)算網(wǎng)格總量約為5.8×106(見圖5).

圖5 求解空間的四面體非結(jié)構(gòu)網(wǎng)格Fig.5 Unstructured mesh of the solving region

2.3 計(jì)算條件

數(shù)值求解時(shí),列車車體為絕熱的無滑移固壁條件,地面采用不可滑移的運(yùn)動(dòng)固壁條件,其運(yùn)動(dòng)速度與來流速度一致,其他邊界處為參考溫度298.5K和參考?jí)毫?101325Pa的均勻來流的遠(yuǎn)場邊界條件.近場流動(dòng)的非定常計(jì)算使用隱式的雙重時(shí)間步(dualtime-stepping)方法,瞬態(tài)計(jì)算的時(shí)間步長?t=5×10?5s. 計(jì)算可解析的信號(hào)最高頻率 fmax=1/(2?t)=10000Hz.

計(jì)算采用了250km/h,350km/h,500km/h三種運(yùn)行速度,對(duì)應(yīng)的馬赫數(shù) Ma分別為0.204,0.286和0.408.由于上述馬赫數(shù)之間的差異較大,并處于弱可壓縮區(qū)間,近場流動(dòng)的數(shù)值求解統(tǒng)一采用了可壓縮形式的控制方程,同時(shí),通過預(yù)處理方法調(diào)節(jié)控制方程的Jacobi系數(shù)矩陣特征值u?c,u+c和u的大小,以減小聲速u±c和流體質(zhì)點(diǎn)速度u之間的差異,減少控制方程系數(shù)矩陣的特征值分散,使問題的剛性降低以提高收斂速度.

2.4 計(jì)算流程

求解過程中,首先使用cubic k-ε RANS模型求得流動(dòng)的定常解,進(jìn)行脈動(dòng)重建后運(yùn)用非線性聲學(xué)求解器(NLAS)進(jìn)行非定常流動(dòng)演進(jìn),并在預(yù)設(shè)噪聲面上采集壓力脈動(dòng).在這一求解步,開始收集聲源信息之前流動(dòng)經(jīng)過了額外的0.5s的非定常演變,以消除脈動(dòng)重建所造成的人為影響,獲得統(tǒng)計(jì)穩(wěn)定的非定常流動(dòng).隨后,在時(shí)間上繼續(xù)推進(jìn)0.5s以便在噪聲面上記錄近場流動(dòng)的壓力脈動(dòng),脈動(dòng)信號(hào)的采集時(shí)長足夠FW-H遠(yuǎn)場積分的需求.最后進(jìn)行壓力脈動(dòng)的FW-H積分以獲得遠(yuǎn)場興趣點(diǎn)處的噪聲信息.

3 計(jì)算結(jié)果及討論

3.1 流速場與流動(dòng)形態(tài)特征

流場特征決定了列車的氣動(dòng)噪聲,為此我們首先分析了列車周圍的流場特征.結(jié)果顯示不同模型在不同速度下的流場表現(xiàn)出基本相似的特征.圖6所示為B頭型在350km/h運(yùn)行速度下,列車周圍流場不同剖面的流線.流線使用當(dāng)?shù)亓飨蛩俣萓與自由來流速度U∞的差值?U進(jìn)行渲染,為清晰對(duì)比當(dāng)?shù)亓飨蛩俣扰c自由來流速度的相對(duì)大小,將速度差?U的取值范圍限定在[?10m/s,10m/s]的區(qū)間中.

圖6 列車車體周圍流場的流線形態(tài),由?U=U?U∞著色渲染Fig.6 Flow fi eld around the train,rendered with?U=U ?U∞

可見,在列車頭車位置,流體從鼻尖駐點(diǎn)沿車鼻表面開始加速,且車鼻近地表面處的流體加速較快,在較短的距離內(nèi),當(dāng)?shù)亓黧w的流向速度便超過了來流速度.而在車鼻表面的其他位置,流體隨著曲面截面的擴(kuò)張不斷加速,并在車鼻截面達(dá)到最大時(shí)達(dá)到了局部的最高速度,隨后,流體沿著具有固定截面形狀的車體發(fā)展,當(dāng)?shù)亓飨蛩俣融吔趤砹魉俣?當(dāng)流動(dòng)發(fā)展到尾車附近時(shí),車體上表面流體在尾車車鼻頂部存在加速段,隨后速度快速下降,進(jìn)入尾流區(qū).在車體后部的兩側(cè)位置,靠近車體表面的流體流向速度在較長距離內(nèi)都低于來流速度,它們與車體上表面的流體一起,匯聚進(jìn)入車體尾部的尾流.

圖 7為運(yùn)行速度 350km/h時(shí),不同頭型長細(xì)比的列車周圍流場的 Q判據(jù)等值面,Q的取值為(U∞/H)2,其中U∞=350km/h為來流速度,H≈3.5m為車體高度.Q判據(jù)等值面由當(dāng)?shù)貕毫與來流壓力p∞=101325Pa的差值?p渲染和著色,為更清楚對(duì)比,?p的取值限制在區(qū)間[?100Pa,100Pa]內(nèi).

根據(jù)圖7所示的計(jì)算結(jié)果,4種車型的Q判據(jù)等值面具有相似的空間外形,主要存在于車體表面附近和尾流區(qū)域,由此可知,細(xì)長列車體流動(dòng)中的主要聲音產(chǎn)生區(qū)域?yàn)檐圀w表面和尾流區(qū).不同車型的Q判據(jù)等值面上的壓力分布規(guī)律也具有一致性:頭車車鼻位置具有高于來流壓力的壓力狀態(tài),隨著流體沿著車鼻表面的加速流動(dòng),車鼻后部的流體速度高于自由來流速度,導(dǎo)致其當(dāng)?shù)貕毫Φ陀谧杂蓙砹鞯膲毫?隨后,中間車體的壓力逐漸趨于自由來流的壓力值,當(dāng)流體運(yùn)動(dòng)到列車尾部時(shí),局部的加速導(dǎo)致當(dāng)?shù)貕毫焖傧陆担?dāng)加速的流體進(jìn)入尾流區(qū)域,其壓力再次迅速上升.不同長細(xì)比的頭型具有相異的車鼻長度和不同的截面變化率,導(dǎo)致當(dāng)?shù)亓黧w的加速和減速狀態(tài)存在差異,在壓力分布上,表現(xiàn)為車頭和車尾位置壓力狀態(tài)的空間差異.

圖7 350km/h時(shí)的Q判據(jù)等值面,Q=(U∞/H)2,并由?p=p?p∞著色渲染.自上而下分別為車型A,車型B,車型C,車型DFig.7 The iso-surface of the Q criterion(Q=(U∞/H)2)at 350km/h,rendered with?p=p? p∞.From top to bottom:shape A,B,C and D

為了更清楚地了解列車周圍的渦結(jié)構(gòu),圖8給出了模型B尾流區(qū)域的Q判據(jù)等值面的局部視圖.從中可以看出,列車尾流中存在明顯的渦列對(duì),并伴隨著一系列的環(huán)狀渦,環(huán)狀渦從車體尾部開始發(fā)展,沿著尾流逐步擴(kuò)張并隨之破碎,尾渦結(jié)構(gòu)的發(fā)展及破碎行為將產(chǎn)生可觀的氣動(dòng)噪聲.而尾渦中的壓力分布整體上是從車體尾部向下游方向降低,逐步恢復(fù)到自由來流的壓力狀態(tài),但其局部的壓力狀態(tài)是十分復(fù)雜的,比如,同一環(huán)狀渦的不同位置,其壓力值存在較大的差異.

圖8 來流速度為350km/h時(shí),Q判據(jù)等值面描述的模型B的尾渦渦列,Q=(U∞/H)2,由?p=p?p∞著色Fig.8 The tail vortex of the shape B at 350km/h,the iso-surface is determined by Q=(U∞/H)2,rendered with?p=p? p∞

近場流動(dòng)中壓力脈動(dòng) p′的空間分布特征如圖9所示,p′為各位置的瞬時(shí)壓力p與時(shí)均壓力ˉp的差,即p′=p?ˉp.從車體中縱剖面處的壓力脈動(dòng)灰度描述(圖9(a))可知,列車表面的壓力脈動(dòng)呈現(xiàn)出沿車體表面正負(fù)壓力交替分布的偶極子聲源特性,而尾流區(qū)的壓力脈動(dòng)是在體空間內(nèi)正負(fù)交替分布的,表現(xiàn)出四極子聲源特性.圖9(b)給出了近場壓力脈動(dòng)的三維空間描述,從中可以清晰地看出,列車表面和尾流中近場壓力脈動(dòng)在空間分布上的正負(fù)交替特征.

圖9 列車(a)中縱剖面壓力脈動(dòng)p′=p?ˉp的灰度描述,(b)Q=(U∞/H)2判據(jù)等值面處的壓力脈動(dòng)Fig.9(a)Pressure perturbation p′=p?ˉp at the middle plane,(b)pressure perturbation at the iso-surface of Q=(U∞/H)2

3.2 氣動(dòng)阻力

圖10 4種長細(xì)比頭型的列車模型的阻力系數(shù)Fig.10 Drag coefficient of four trains with di ff erent slenderness ratio

高速列車頭尾車外形的變化將影響周圍流動(dòng)的發(fā)展和分離,不僅影響列車的氣動(dòng)噪聲,同樣影響其氣動(dòng)性能,特別是阻力性能,列車氣動(dòng)噪聲的優(yōu)化不應(yīng)以氣動(dòng)阻力的過度增加為代價(jià).為此,分析了A,B,C,D四種列車模型在不同運(yùn)行速度下的阻力系數(shù),如圖10所示.在三種運(yùn)行速度下,頭車和尾車阻力系數(shù)都隨著長細(xì)比λ的增大而減小,這可能是由于長細(xì)比的增大,導(dǎo)致車鼻處的橫截面變化率降低,從而使周圍流動(dòng)的穩(wěn)定性增強(qiáng),分離減弱.但頭、尾車氣動(dòng)阻力的明顯下降段發(fā)生在[λA,λC]區(qū)間,當(dāng)λ>λC時(shí),頭尾車阻力系數(shù)的下降趨勢將顯著減弱.中間車廂的阻力系數(shù)與頭型長細(xì)比之間沒有呈現(xiàn)出明確的規(guī)律,由于C型車與D型車在頭尾車阻力上差異很小,而前者的中間車廂阻力系數(shù)較小于后者,從而C型車有著最小的總阻力.

圖10 4種長細(xì)比頭型的列車模型的阻力系數(shù)(續(xù))Fig.10 Drag coefficient of four trains with di ff erent slenderness ratio(continued)

與模型A相比,在250km/h的運(yùn)行速度下,模型B,C,D的總阻力分別降低2.31%,8.56%,7.70%,在350km/h的運(yùn)行速度下,總阻力分別下降2.50%,8.95%,8.01%,而在500km/h時(shí),其總阻力分別下降2.81%,9.36%,8.45%.可見,在一定程度上,頭型長細(xì)比增大有利于氣動(dòng)阻力的減小,但當(dāng)長細(xì)比增大到一定程度時(shí),氣動(dòng)阻力反而有所增大;而且頭型長細(xì)比對(duì)列車氣動(dòng)阻力的影響隨著運(yùn)行速度的增加而增強(qiáng).

3.3 遠(yuǎn)場噪聲

遠(yuǎn)場噪聲的求解基于近場聲源面的 FW-H積分,由于湍動(dòng)能的主要集中區(qū)域是氣動(dòng)噪聲的主要聲源區(qū),本文選取最大湍動(dòng)能的1%等值面作為近場聲源面的參照,最終所選取的聲源面位置如圖11所示,聲源面長度范圍為[?45m,60m],寬度為[?3m,3m],高度區(qū)間為[0,5.5m].

高速列車鐵路沿線噪聲的評(píng)估通常基于距離軌道中心線25m、高度3.5m處的聲壓級(jí),為探討列車沿線的噪聲特性,沿列車長度方向設(shè)置了如圖12所示的聲壓探點(diǎn),探點(diǎn)總數(shù)為13個(gè),各探點(diǎn)在x方向上的間距為10m,從車頭向車尾方向依次編號(hào),其中P7探點(diǎn)與列車車體中點(diǎn)對(duì)齊.

圖11 FW-H聲源積分面空間位置的選取.上方彩圖:1%最大湍動(dòng)能的等值面(由流向速度U渲染著色).下方虛線線框:作為FW-H積分邊界的聲源面的空間位置Fig.11 The integral surface for FW-H integral,top:the iso-surface of 1%turbulence kinetic energy rendered with stream velocity U,bottom:the wire frame shows the location of the actual used integral surface in the simulation

圖12 列車長度方向的聲壓探點(diǎn)分布示意圖Fig.12 The location of the sound pressure probes in the stream direction

圖13為模型B的3個(gè)沿線探點(diǎn)P1,P7和P13在3種運(yùn)行速度下的聲壓級(jí)頻譜.可以看出,其氣動(dòng)噪聲分布在很寬的頻率范圍內(nèi),并不存在明顯的主峰,隨著運(yùn)行速度的提高,各探點(diǎn)位置的聲壓級(jí)在整個(gè)頻段內(nèi)都相應(yīng)地增加.

圖13 模型B的3個(gè)沿線探點(diǎn)的聲壓級(jí)Fig.13 SPL of three probes for the case of shape B

4種車型的沿線噪聲總聲壓級(jí) (overall sound pressure level,OASPL)隨運(yùn)行速度的變化如圖14,各型列車沿線噪聲在車長方向上的分布具有相似的特征,列車的沿線噪聲從車體前部向尾部方向上升,并且前部探點(diǎn)的總聲壓級(jí)趨于線性增長.隨著運(yùn)行速度的增大,各探點(diǎn)的總聲壓級(jí)相應(yīng)增大,同時(shí),越靠近車體尾部,探點(diǎn)的總聲壓級(jí)的增加幅度越大,這可能是由于速度的增大,使得列車后部產(chǎn)生了更高強(qiáng)度的尾流流動(dòng).

圖14 4種車型沿線噪聲總聲壓級(jí)隨運(yùn)行速度的變化Fig.14 The OASPL at di ff erent locations

圖14 4種車型沿線噪聲總聲壓級(jí)隨運(yùn)行速度的變化(續(xù))Fig.14 The OASPL at di ff erent locations(continued)

圖15給出了列車沿線噪聲總聲壓級(jí)與頭型長細(xì)比之間的關(guān)系,總體上講,4種頭型中列車沿線噪聲總聲壓級(jí)從高到低依次為A,C,B,D,且A與C具有相近的沿線噪聲聲壓級(jí),而B和D具有相近的聲壓級(jí).同時(shí),4種頭型沿線噪聲的差異隨著運(yùn)行速度的提高而減少,也就是說,頭型對(duì)列車氣動(dòng)噪聲的影響隨著列車速度的增大而減弱.

圖15 4種長細(xì)比頭型的列車沿線噪聲總聲壓級(jí)對(duì)比Fig.15 The comparison of the OASPL for di ff erent shapes

圖16 沿列車周向分布的遠(yuǎn)場壓力探點(diǎn)及相應(yīng)的總聲壓級(jí)分布Fig.16 The OASPL distribution around the train

為更全面地比較頭型對(duì)列車遠(yuǎn)場噪聲聲壓級(jí)和方向性的影響,圍繞列車車體設(shè)置了如圖16所示的遠(yuǎn)場壓力探點(diǎn),探點(diǎn)位于距離地面3.5m的水平面內(nèi),沿半徑為150m的圓周等間隔分布,間隔角度?θ=15°,圓周的圓心與列車的中心具有相同的x和y坐標(biāo)值.

沿列車周向分布遠(yuǎn)場壓力探點(diǎn)的總聲壓級(jí)分布如圖16所示,總聲壓級(jí)的圓周分布在列車兩側(cè)是對(duì)稱的,聲壓級(jí)從高到低依次為A,C,B,D,其中,B和D的聲壓級(jí)在車頭位置(θ=0°)存在小的差異,但在整個(gè)圓周上都是十分接近的,這與前面的沿線噪聲結(jié)果一致.

4 結(jié)論

本文基于簡化的列車細(xì)長體模型,在250km/h,350km/h,500km/h三種運(yùn)行速度下,探討了頭型長細(xì)比分別為 0.75λ0,1.0λ0,1.25λ0和 1.5λ0四種頭型的列車氣動(dòng)性能和氣動(dòng)噪聲特性,得到如下主要結(jié)論:

(1)列車頭尾車的氣動(dòng)阻力隨著頭型長細(xì)比的增大而減少,但當(dāng)頭型長細(xì)比超過其一定值后,其對(duì)氣動(dòng)阻力的影響將減弱;且頭型長細(xì)比對(duì)列車總氣動(dòng)阻力的影響隨運(yùn)行速度的增加而增強(qiáng).在本文研究中,C、D兩種頭型列車的頭尾車氣動(dòng)阻力差異較小,整體上長細(xì)比為1.25λ0的C型車在四種車型中具有最低的氣動(dòng)總阻力.

(2)不同長細(xì)比的各型列車,其沿線噪聲在車長方向上的分布具有相似的特征,列車的沿線噪聲從車體前部到尾部逐漸上升,并且前部探點(diǎn)的總聲壓級(jí)趨于線性增長.隨著運(yùn)行速度的增加,各探點(diǎn)的總聲壓級(jí)相應(yīng)增大,同時(shí),越靠近車體尾部,探點(diǎn)的總聲壓級(jí)的增加幅度越大.

(3)四種頭型沿線噪聲之間的差異隨著運(yùn)行速度的提高而減少,也就是說,頭型對(duì)列車氣動(dòng)噪聲的影響作用隨著列車速度的增大而減弱.

(4)對(duì)列車沿線噪聲和遠(yuǎn)場圓周噪聲的分析表明,長細(xì)比最小的A型車的遠(yuǎn)場聲壓級(jí)最高,其次是車型C,也就是說,具有最優(yōu)氣動(dòng)阻力的C并不具有最低的遠(yuǎn)場噪聲.長細(xì)比為1.0λ0的B和1.5λ0的D在空間上具有十分接近的聲壓級(jí)分布,在車頭等局部位置,D的總聲壓級(jí)略微低于B.

(5)綜合考慮氣動(dòng)阻力和氣動(dòng)噪聲,長細(xì)比最大的D型車綜合性能較優(yōu),但差異并不顯著,因此在不考慮微氣壓波等因素的條件下,簡單增加車頭長細(xì)比并不一定能帶來明顯的氣動(dòng)阻力和噪聲性能提升.

1 David Thompson,et al.Railway Noise and Vibration:Mechanisms,Modelling and Means of Control.Amsterdam:Elsevier,2008

2 King WF.III.A precis of developments in the aeroacoustics of fast trains.Journal of Sound and Vibration,1996,193(1):349-358

3 楊國偉,魏宇杰,趙桂林等.高速列車的關(guān)鍵力學(xué)問題.力學(xué)進(jìn)展,2015,45:201507(Yang Guowei,Wei Yujie,Zhao Guilin et al.,Research progress on the mechanics of high speed rails.Advances in Mechanics,2015,45:201507(in Chinese))

4 Talotte C.Aerodynamic noise:a critical survey.Journal of Sound and Vibration,2000,231(3):549-562

5 張軍,孫幫成,郭濤等.高速列車整車氣動(dòng)噪聲及分布規(guī)律研究.鐵道學(xué)報(bào),2015,37(2):10-17(Zhang Jun,Sun Bangcheng,Guo Tao,et al.Research on aerodynamic noise radiated from whole body surface of high-speed train and its distribution.Journal of the China Railway Society,2015,37(2):10-17(in Chinese))

6 張亞東,張繼業(yè),李田.高速列車整車氣動(dòng)噪聲聲源特性分析及降噪研究.鐵道學(xué)報(bào),2016,38(7):40-49(Zhang Yadong,Zhang Jiye,Li Tian,Research on aerodynamic noise source characterization and noise reduction of high-speed trains vehicle.Journal of the China Railway Society,2016,38(7):40-49(in Chinese))

7 Mellet C,L′etourneaux F,Poisson F,et al.High speed train noise emission:Latest investigation of the aerodynamic/rolling noise contribution.Journal of Sound and Vibration,2006,293(3):535-546

8 SunZX,ZhangY,YangGW.SurrogateBasedOptimizationofAerodynamic Noise for Streamlined Shape of High Speed Trains.Applied Sciences-Basel,2017,7(2):192

9 劉加利,田愛琴,杜健等.頭部控制線形狀對(duì)高速列車氣動(dòng)噪聲的影響.中國鐵路,2014,11:58-62(Liu Jiali,Tian Aiqin,Du Jian,et al.In fl uence of the head car control line shape on the aerodynamic noise of high speed train.Chinese Railways,2014,11:58-62(in Chinese))

10 Kitagawa T,Nagakura K.Aerodynamic noise generated by Shinkansen cars.Journal of Sound and Vibration,2000,231(3):913-924

11 Torii A,Ito J.Development of the series 700 Shinkansen train-set(improvement of noise level).Japanese Railway Engineering,2000,40(1)

12 Maeda T,Kinoshita M,Kajiyama H,et al.Aerodynamic drag of shinkansen electric cars(series 0,series 200,series 100).Railway Technical Research Institute,Quarterly Reports,1989,30(1)

13 Ido A,Iida M,Maeda T.Wind tunnel tests for nose and tail of train.RTRI Report JNR,1993,7(7)(in Japanese)

14 喻華華.高速列車氣動(dòng)噪聲產(chǎn)生與控制的機(jī)理研究.[博士論文].北京:中國科學(xué)院力學(xué)研究所,2013(Yu Huahua.On the mechanism of generation and control of aerodynamic noise for highspeed trains.[PhD thesis].Beijing:Institute of Mechanics,Chinese Academy of Sciences,2013(in Chinese))

15 王成強(qiáng),邢海英,鄭繼峰.基于CAA的高速動(dòng)車組氣動(dòng)噪聲仿真研究.華東交通大學(xué)學(xué)報(bào),2015,32(1):9-15(Wang Chengqiang,XingHaiying,ZhengJifeng,Simulationstudyonaerodynamicnoise of the high speed trains based on CAA.Journal of East China Jiaotong University,2015,32(1):9-15(in Chinese))

16 潘忠,陸森林.高速列車外流場氣動(dòng)噪聲數(shù)值模擬研究.重慶理工大學(xué)學(xué)報(bào)(自然科學(xué)),2016,30(5):8-14(Pan Zhong,Lu Senlin.Numerical simulation of aerodynamic noise for high-speed train in exterior fl ow fi eld.Journal of Chongqing University of Technology(Natural Science),2016,30(5):8-14(in Chinese))

17 Ding SS,Li Q,Tian AQ,et al.Aerodynamic design on high-speed trains.Acta Mechanica Sinica,2016,32(2):215-232

18 張亮,張繼業(yè),李田等.超高速列車流線型頭型多目標(biāo)優(yōu)化設(shè)計(jì).機(jī)械工程學(xué)報(bào),2017,53(2):106-114(Zhang Liang,Zhang Jiye,Li Tian,et al.Multi-objective optimization design of the streamlined head shape of super high-speed trains.Journal of Mechanical Engineering,2017,53(2):106-114(in Chinese))

19 Lee JS,Kim JH.Approximate optimization of high-speed train nose shape for reducing micropressure wave.Structural and Multidisciplinary Optimization,2008,35(1):79-87

20 常興華,馬戎,張來平等.基于計(jì)算流體力學(xué)的“虛擬飛行”技術(shù)及初步應(yīng)用.力學(xué)學(xué)報(bào),2015,47(4):596-604(Chang Xinghua,Ma Rong,Zhang Laiping,et al.Study on CFD-based numerical virtual fl ight technology and preliminary application.Chinese Journal of Theoretical and Applied Mechanics,2015,47(4):596-604(in Chinese))

21 胡守超,崔凱,李廣利等.基于實(shí)驗(yàn)設(shè)計(jì)方法的高超聲速飛機(jī)前緣型線優(yōu)化分析.力學(xué)學(xué)報(bào),2016,48(2):290-299(Hu Shouchao,Cui Kai,Li Guangli,et al.Optimization and analysis of the leading edge shape for hypersonic airplanes based on doe methods.Chinese Journal of Theoretical and Applied Mechanics,2016,48(2):290-299(in Chinese))

22 李廣利,崔凱,肖堯等.高壓捕獲翼前緣型線優(yōu)化和分析.力學(xué)學(xué)報(bào),2016,48(4):877-885(Li Guangli,Cui Kai,Xiao Yao,et al.Leading edge optimization and parameter analysis of high pressure capturing wings.Chinese Journal of Theoretical and Applied Mechanics,2016,48(4):877-885(in Chinese))

23 馮峰,郭力,王強(qiáng).高亞聲速湍流噴流氣動(dòng)噪聲數(shù)值分析.力學(xué)學(xué)報(bào),2016,48(5):1049-1060(Feng Feng,Guo Li,Wang Qiang.Numerical investigation of noise of a high subsonic turbulent jet.Chinese Journal of Theoretical and Applied Mechanics,2016,48(5):1049-1060(in Chinese))

24 Batten P,Ribaldone E,Casella M,et al.Towards a generalized non-linear acoustics solver.AIAA-2004-3001//10th AIAA/CEAS Aeroacoustics Conference.Manchester,United Kingdom,May 10-12(2004)

25 Ffowcs Williams JE,Hawkings DL.Sound generation by turbulence and surfaces in arbitrary motion.Philosophical Transactions of the Royal Society of London,Series A,1969,264:321-342

26 Batten P,Goldberg U,Chakravarthy S.Reconstructed sub-grid methods for acoustics predictions at all Reynolds numbers//AIAA Paper,2511,2002

27 Batten P,Enrico R,Chakravarthy S.Towards a generalized nonlinear acoustics solver//AIAA Paper,2004-3001,2004

28 Farassat F.Discontinuities in aerodynamics and aeroacoustics:The concept and applications of generalized derivatives.Journal of Sound and Vibration,1977,55(2):165-193

29 Hold R,Brenneis A,Eberle A,et al.Numerical simulation of aeroacoustic sound generated by generic bodies placed on a plate:Part I-prediction of aeroacoustic sources//AIAA paper,1999,20:40-60

30 Siegert R,Schwarz V,Reichenberger J.Numerical simulation of aeroacoustic sound generated by generic bodies placed on a plate:Part II-prediction of radiated sound pressure//AIAA paper,1999,20:62-72

31 Yu HH,Li JC,Zhang HQ.On aerodynamic noises radiated by the pantograph system of high-speed trains.Acta Mechanica Sinica,2013,29(3):399-410

32 Guo DL,Shang KM,Zhang Y,et al.In fl uences of affiliated components and train length on the train wind.Acta Mechanica Sinica,2016,32(2):191-205

STUDY ON THE INFLUENCE OF THE NOSE SLENDERNESS RATIO OF HIGH-SPEED TRAIN ON THE AERODYNAMIC NOISE1)

An Yi?Mo Huangrui?Liu Qingquan?,2)?(LMFS,Institute of Mechanics,Chinese Academy of Sciences,Beijing 100190,China)
?(School of Aerospace Engineering,Beijing Institute of Technology,Beijing 100081,China)

In the high-speed train design,the nose shape is a crucial control factor in fl uencing not only aerodynamic performance but also the aerodynamic noise.In the engineering practice,the nose shape becomes more and more slender along with the increasing of the design speed,e.g.the Japanese high-speed maglev train L0 series even has a 15m long slender nose(the slenderness ratio reach to 8.8).This study aims to discuss the in fl uence of the slenderness ratio of the nose shape on the aerodynamic noise.The hybrid numerical method of nonlinear acoustics solver(NLAS)and Ffowcs Williams-Hawkings(FW-H)acoustic analogy method is employed to study the aerodynamics noise characteristics.The numerical method is validated with a standard wind mirror test case and a set of acoustics wind tunnel experiments of the CRH380A train.The shape of the CRH380A train is chosen as a bench mark,and four di ff erent nose shapes of di ff erent slenderness ratio under di ff erent running speed situation are studied with numerical simulation.The fl ow fi eld,aerodynamic drag,and the aerodynamic noise are obtained and discussed.The result shows that the total drag decrease with the increase of the slenderness ratio,and this e ff ect enhances when the train speed increases.However,the in fl uence of the slenderness ratio on the aerodynamic noise is much complex as no simple trend is observed.Considering both the aerodynamic and aeroacoustics characteristics,the train with the most slender nose shape is the best while this advantage is not notable compared with the second-best.Thus,simply increase the slenderness does not necessarily result in better aerodynamic noise performance if the e ff ect of tunnel boom is not considered.

high-speed rail,aerodynamic noise,aerodynamic drag,slenderness ratio

O354.1,TB533+.2,U260.16,U266

A

10.6052/0459-1879-17-126

2017–04–17收稿,2017–08–11 錄用,2017–08–20 網(wǎng)絡(luò)版發(fā)表.

1)中國科學(xué)院知識(shí)創(chuàng)新工程方向性資助項(xiàng)目(KJCX2-EW-L02-1).

2)劉青泉,教授,主要研究方向:環(huán)境流體力學(xué).E-mail:liuqq@bit.edu.cn

安翼,莫晃銳,劉青泉.高速列車頭型長細(xì)比對(duì)氣動(dòng)噪聲的影響.力學(xué)學(xué)報(bào),2017,49(5):985-996

An Yi,Mo Huangrui,Liu Qingquan.Study on the in fl uence of the nose slenderness ratio of high-speed train on the aerodynamic noise.Chinese Journal of Theoretical and Applied Mechanics,2017,49(5):985-996

主站蜘蛛池模板: 少妇精品网站| 亚洲精品无码久久毛片波多野吉| 亚洲国产一区在线观看| 狠狠亚洲五月天| 欧美伊人色综合久久天天| 国产色图在线观看| 精品国产三级在线观看| 91精品啪在线观看国产60岁 | 欧美.成人.综合在线| 在线观看热码亚洲av每日更新| 国产亚洲精品97AA片在线播放| 国产网站免费| 久久性妇女精品免费| 91视频首页| 国产成人高清精品免费软件| 成人一区专区在线观看| 亚洲成人播放| 国产第一色| 久久婷婷色综合老司机| 一本大道香蕉久中文在线播放| 亚洲爱婷婷色69堂| 免费无码一区二区| 免费观看亚洲人成网站| 国产尤物jk自慰制服喷水| 激情综合网址| 国产亚洲第一页| 国产一区亚洲一区| 国产精品微拍| 亚洲最大福利网站| 国产精品人人做人人爽人人添| 在线va视频| 亚洲视频a| Aⅴ无码专区在线观看| 国产激爽爽爽大片在线观看| 97视频免费在线观看| 亚洲色中色| 在线色综合| 久久99久久无码毛片一区二区| 美女无遮挡拍拍拍免费视频| 成人在线不卡视频| 97国产成人无码精品久久久| 国产性生大片免费观看性欧美| 精品欧美一区二区三区久久久| 亚洲欧美综合另类图片小说区| 亚洲天天更新| 嫩草国产在线| 亚洲A∨无码精品午夜在线观看| 国内精品自在欧美一区| 亚洲AV一二三区无码AV蜜桃| 国产永久免费视频m3u8| 黄色网在线| 国产人成乱码视频免费观看| 国产免费高清无需播放器| 成人中文在线| 日本免费高清一区| 亚洲福利网址| YW尤物AV无码国产在线观看| 国产v精品成人免费视频71pao| 国产精品99r8在线观看| 精品国产电影久久九九| 亚洲日韩每日更新| 国产农村妇女精品一二区| 国产亚洲精品91| 中文国产成人精品久久| 国产99精品久久| 国产精品一区在线观看你懂的| 久久精品电影| 国产91精品调教在线播放| 国产永久无码观看在线| 欧美在线导航| 亚洲天堂777| 国产在线第二页| 人妻无码中文字幕第一区| 久久毛片免费基地| av在线无码浏览| 狠狠做深爱婷婷综合一区| 国产成人永久免费视频| 色播五月婷婷| 国产高潮视频在线观看| 色网站在线免费观看| 伊人久久大线影院首页| 欧美高清三区|