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

彈性離散支撐鋼軌的振動和聲輻射特性研究

2022-07-12 04:26:54劉清源何遠鵬陳宗平圣小珍
鐵道學報 2022年6期
關鍵詞:振動模型

劉清源,成 功,2,何遠鵬,陳宗平,韓 健,圣小珍

(1.西南交通大學 牽引動力國家重點實驗室,四川 成都 610031;2.華東交通大學 軌道交通基礎設施性能監測與保障國家重點實驗室,江西 南昌 330013;3.西南交通大學 機械工程學院,四川 成都 610031;4.上海工程技術大學 城市軌道交通學院, 上海 201620)

隨著軌道交通的快速發展,鐵路噪聲與振動問題受到了社會的廣泛關注,鋼軌輻射噪聲作為鐵路噪聲的主要來源之一,其振動和聲輻射特性的建模和預測也受到了學界的廣泛關注。

Remington[1]在建立的輪軌噪聲預測模型中使用連續彈性支撐的歐拉梁模型模擬鋼軌,Thompson[2]在Remington的基礎上,使用離散彈性支撐的鐵木辛柯梁模型模擬鋼軌,并將其納入TWINS軟件。由于簡單的鐵木辛柯梁模型在2 000 Hz以上的高頻已不再適用,Wu等[3]提出一種采用雙層鐵木辛柯梁模擬鋼軌的方法,與單梁模型相比,雙層梁模型能反映高頻激勵下的截面內變形。Knothe等[4]采用梁模型模擬軌頭和軌底、板模型模擬軌腰,使用梁和板的組合模型得到鋼軌高頻響應。萬淑敏等[5]建立離散彈性支撐的三維實體有限元鋼軌模型,并得到鋼軌的聲輻射效率。

Remington[1]的噪聲預測模型把鋼軌作為線聲源處理,根據鋼軌的聲輻射效率和振動速度計算鋼軌的輻射聲功率。為了獲得高頻聲輻射,羅斌等[6]使用統計能量分析得到鋼軌的1/3倍頻程輻射噪聲。Sheng等[7-8]基于傅里葉變換,提出一種基于波數域的方法,求解周期性結構在移動和靜止諧載荷作用下的響應,并使用2.5維邊界元法得到鋼軌在移動和靜止諧載荷下的聲輻射特性。Thompson[9]總結了鋼軌噪聲和軌道衰減率的關系式,根據鋼軌截面的輻射效率和振動衰減率,計算無限長鋼軌的輻射聲功率。

研究鋼軌輻射分為兩步,計算鋼軌振動和計算鋼軌輻射。計算鋼軌振動一般采用解析法或者有限元法,可根據需要視鋼軌為梁單元模型、多層梁單元模型或者實體單元模型;計算鋼軌聲輻射主要采用聲學邊界元法。二維邊界元法忽略了振動沿軌道方向的變化,與實際情況相去甚遠,但其算出的鋼軌輻射效率可結合軌道衰減率計算鋼軌噪聲。三維邊界元受模型網格的影響很大,網格過密會造成計算效率過低,網格過疏會導致模型在高頻段的計算結果可信度降低,對研究無限長鋼軌的振動特性產生限制。鋼軌可以視作在一個方向上無限長的同質、均勻且橫截面形狀不變的波導體,利用這一特點可以使用于計算鋼軌聲輻射的三維邊界積分方程得到簡化。通過將無限長方向的空間變量從空間域變換到波數域,就可以實現在波數域中解耦求解,這就是2.5維邊界元法。應用這種方法計算鋼軌的聲輻射可以克服以上兩種模型的缺點,在準確反映振動沿鋼軌傳播的同時提高計算效率[7-8]。

為了探究基于不同模型和方法所計算出的振動和聲輻射特性的差異及其原因,確定它們的適用范圍,為研究鋼軌振動和聲輻射提供參考,并分析剛性軌道板對鋼軌輻射聲功率的影響,本文采用三種鋼軌模型計算無限長鋼軌在靜止和移動諧載荷下的振動特性和軌道衰減率,并使用2.5維聲學邊界元方法和衰減率與噪聲之間的關系公式計算無限長鋼軌在靜止諧載荷下的輻射聲功率級,并考慮剛性軌道板的影響。

1 鋼軌聲振模型建立

1.1 鋼軌振動模型的建立

本文建立的軌道模型包括UIC60鋼軌及軌下扣件系統,其中扣件系統模擬為結構阻尼彈簧,鋼軌分別模擬為有限長鐵木辛柯梁單元、有限長三維實體單元-鐵木辛柯梁單元組合結構及無限長鐵木辛柯梁,軌道結構如圖1所示。

圖1 軌道結構

前兩種結構都是有限長結構,長度為210 m。值得注意的是,第二種模擬方法是在21 m長的三維實體有限元鋼軌兩端分別連接一段梁單元鋼軌,一方面可以更準確地得到鋼軌的振動特性,另一方面在兼顧計算效率的同時盡量減小截斷邊界的反射效應。該方法在Ansys中進行模擬,在梁單元和實體單元之間使用多點約束(MPC)連接,如圖2所示。其原理是將連接截面處的梁單元節點和實體單元節點的自由度進行剛性連接,用于耦合單元之間的載荷傳遞。

圖2 MPC連接實體單元和梁單元

值得注意的是,要研究無限長鋼軌的振動和聲輻射特性,需要探究多長的有限元模型能夠滿足要求。本文使用不同長度鐵木辛柯梁模型,將端部導納與原點導納的比值作為衡量標準,如圖3所示。

圖3 不同長度鋼軌端部與原點導納比值

從圖3可以看出,較大的比值出現在400~3 000 Hz。隨著鋼軌長度的增加,鋼軌端部導納減小,當鋼軌長度為200 m時,導納比值最大不超過20%,此時得到的原點導納與后文采用的無限長鐵木辛柯梁的計算結果基本一致,且不會出現因截斷邊界反射產生的波動,所以導納比值不超過20%時(即相差14 dB)可以將有限長的鋼軌近似為無限長鋼軌。

第三種結構是將鋼軌模擬為一根無限長的鐵木辛柯梁,該梁在周期性的扣件位置被離散支撐,其鋼軌在單位高速移動諧載荷作用下的速度振動響應(t=0)可以寫為[8]

(1)

1.2 鋼軌聲輻射計算方法

Thompson[9]對無限長鋼軌的聲功率計算公式進行了簡化,得到聲功率級與軌道衰減率之間的關系式為

(2)

式中:Wref為參考聲功率,Wref=10-12W;vref為參考速度,vref=10-9m/s;σ為輻射效率;ρ0c0為空氣特性阻抗;P為鋼軌橫截面周長;v(0)為激勵點處的鋼軌速度;δ為軌道衰減率,dB/m。

該聲功率計算方法認為鋼軌下部的支撐是連續支撐,因此認為其振動衰減率與波數β的虛部成正比,但由式(1)易知離散支撐下的鋼軌在每跨不同位置的振動幅值都不相同。因此,本文在文獻[7]和文獻[8]的基礎上將無限長鐵木辛柯梁鋼軌和2.5維邊界元相結合,得到更符合真實情況的鋼軌輻射聲功率,其計算公式為

(3)

(4)

(5)

(6)

式中:K0為修正的0階第二類貝塞爾函數;k0=2πf/c0,c0為聲速,f為頻率。

1.3 模型參數的選取

本文只考慮鋼軌垂向振動特性和垂向振動的輻射聲功率,因此僅設置了軌下墊板的垂向剛度。鋼軌和軌下墊板的參數見表1。

表1 軌道參數

采用梁單元模擬鋼軌受力情況與實際情況有差異,梁模型中集中載荷施加在鋼軌的中心,彈簧單元實際支撐的位置也是鋼軌中心。實體模型受力需要與梁模型受力一致,所以對鋼軌截面劃分網格時,保證截面中心位置有一個節點,實現對鋼軌的中心處施加集中載荷和提供彈簧支撐,如圖4所示。

圖4 有限元實體模型受力示意

使用有限元法計算鋼軌振動和輻射聲功率對單元長度有要求,需要單元的最大邊長不超過鋼軌傳遞最短波長的1/6。以自由鐵木辛柯梁鋼軌為例,無限長鋼軌彎曲波的頻散曲線如圖5所示。由于梁模型只能傳遞彎曲波,因此圖5只有一條曲線,但是考慮鋼軌的其他傳播波(在這些振動波中,鋼軌橫截面發生變形)在固定頻率時的波數要小于彎曲波的波數,當鋼軌網格滿足彎曲波要求時,這些傳播波也天然滿足。

圖5 自由鋼軌彎曲波頻散曲線

2 鋼軌垂向原點位移導納特性分析

原點導納是研究振動特性的一個常用指標,是結構在單位諧載荷下的響應(位移、速度、加速度)。分別對有限長鐵木辛柯梁、有限長三維實體-鐵木辛柯梁組合結構及無限長鐵木辛柯梁施加單位垂向諧載荷,載荷作用于相鄰扣件的跨中位置,載荷頻率范圍為10~5 000 Hz,步長為10 Hz。對實體-梁組合模型,荷載施加在跨中截面的中心軸處。根據文獻[10],使用有限元法計算實體單元的原點導納會產生奇異性,本文為避免這種情況,將與激勵點相隔一個單元長度(50 mm)的點作為響應點,將該響應點的響應作為原點導納,得到原點位移導納如圖6所示。

圖6 原點位移導納

圖6中,導納第一個峰值出現在170 Hz,是鋼軌整體在扣件上的垂向振動特征頻率,是軌道系統的一階共振頻率,其可以計算為

(7)

式中:k為單位長度鋼軌的剛度,N/m2;m為單位長度鋼軌的質量,kg/m。

更高頻率的兩個峰值(1 080 Hz和2 940 Hz)對應鋼軌的兩階pinned-pinned頻率。pinned-pinned頻率主要受鋼軌彎曲剛度和扣件間距的影響,與扣件支撐鋼軌的寬度也有關系但影響不大。在pinned-pinned頻率處,當載荷作用在跨中時,原點導納會出現峰值,作用在扣件上方會出現谷值。

pinned-pinned模態產生的主要原因是駐波效應,即該頻率下的半波長等于支撐間距,駐波節點正處于扣件位置。鐵木辛柯梁的垂向振動pinned-pinned頻率為

(8)

式中:n0為鋼軌pinned-pinned模態階數;L為軌枕間距,m。

通過對比三種模型可以發現:

(1)有限長鐵木辛柯梁模型和無限長鐵木辛柯梁模型的原點導納基本保持一致,說明在計算原點導納時,200 m的有限元梁模型已經可以很好地還原鋼軌的無限長特性。

(2)三種模型在頻率小于3 000 Hz時的垂向位移導納基本一致,當頻率大于3 000 Hz以后,實體-梁組合模型的垂向位移導納開始顯著大于另外兩種梁模型的位移導納。因為鐵木辛柯梁模型的使用需要滿足對截面變形的一些假設,受到高頻激勵后,這些假設不再有效,因此梁的模型已經不再適合。實體-梁模型受到高頻載荷作用,鋼軌根部發生翹曲,橫截面發生變形,不再維持在一個平面,與兩個梁模型相比變形也更大。

3 軌道衰減率分析

軌道衰減率是評價鋼軌振動縱向傳播特性的重要指標,其對鋼軌的聲輻射至關重要[11]。軌道衰減率代表了振動沿鋼軌方向減弱的能力,軌道衰減率大說明該頻率下的振動沿鋼軌方向傳播能力差,軌道系統的動態阻尼大。鋼軌的支撐方式和扣件阻尼都對軌道衰減率有很大影響。

工程上通常根據歐洲標準EN 15641:2008+A1:2010[12]計算軌道衰減率,該方法采用移動激勵的方式,在不同的位置用力錘進行敲擊,測出不同距離的導納,得到軌道衰減率為

(9)

式中:A(xn)為鋼軌原點(0點)施加單位載荷時第n點的導納函數;A(x0)為在0點施加單位載荷時原點的導納函數;Δxn為第n個測點到第(n-1)個測點的距離。值得注意的是,工程測試時錘擊點和測試點都位于鋼軌軌頂。

本文采用該方法對軌道衰減率進行評價,并且將與激勵點相隔一個單元長度(50 mm)的點作為0點,這也更符合工程測試實際。三種模型的軌道衰減率如圖7所示。

圖7 三種模型的軌道衰減率

三種模型計算的軌道衰減率在低于3 000 Hz的范圍內基本一致,低頻段的軌道衰減率接近9 dB/m,該范圍內的振動會沿軌道迅速耗散,在pinned-pinned頻率附近出現峰值。

在高于3 000 Hz的頻域范圍,三種模型的軌道衰減率都有升高的趨勢,但實體-鐵木辛柯梁模型的上升幅度遠大于另外兩種梁模型。原因是在高頻振動下,實體-鐵木辛柯梁模型的激勵原點截面會發生變形,且某些頻率下變形后截面不再符合平截面假定,而梁模型的截面為剛性,所以實體-鐵木辛柯梁模型鋼軌在高頻范圍的動態阻尼較鐵木辛柯梁模型偏大。

通常情況下,軌道衰減率往往通過靜止諧載荷激勵得到,而在實際中,鋼軌受到的是移動輪軌力激勵,移動諧載荷激勵下的鋼軌振動,可以通過文獻[7]中式(42)計算移動諧載荷下的位移導納得到移動諧載荷下的軌道衰減率,如圖8所示。從圖8可以看出,載荷前和載荷后的軌道衰減率有較大區別,前者整體大于后者,同時前者大于靜止諧載荷軌道衰減率,后者小于靜止諧載荷軌道衰減率。靜止諧載荷會在pinned-pinned頻率處(1 080 Hz)出現峰值,而在移動諧載荷作用時pinned-pinned頻率會出現分叉(載荷前變為1 020、1 160 Hz,載荷后變為1 000、1 180 Hz)。

圖8 移動諧載荷與靜止諧載荷的軌道衰減率

4 鋼軌輻射聲功率分析

4.1 實體-鐵木辛柯梁鋼軌和無限長鐵木辛柯梁聲輻射對比

圖7表明實體-鐵木辛柯梁鋼軌和鐵木辛柯梁鋼軌在高頻段的軌道衰減率有較大差異,根據式(2)易知軌道衰減率對鋼軌的輻射聲功率具有較大影響。基于式(2),可以得到單位垂向諧載荷作用下的實體-鐵木辛柯梁鋼軌和無限長鐵木辛柯梁鋼軌輻射聲功率,如圖9所示。載荷作用于相鄰扣件之間的跨中位置,載荷頻率范圍為10~5 000 Hz,步長為10 Hz。對于實體-梁組合模型,載荷施加在跨中截面的中心軸處。

圖9 有限長實體-鐵木辛柯梁組合模型和無限長鐵木辛柯梁模型的輻射聲功率對比

從圖9可以看出,有限長實體-鐵木辛柯梁組合模型的輻射聲功率在高于3 000 Hz時明顯小于無限長鐵木辛柯梁模型所計算的輻射,隨著頻率的增加,這種差異也隨之增大,在5 000 Hz時差異可以達到6 dB。由式(2)可知,當軌道衰減率增加時,輻射聲功率會減小,原點速度導納增加時,輻射聲功率會增大。由圖6和圖7可知,頻率高于3 000 Hz時,實體-鐵木辛柯梁模型的原點速度導納大同時軌道衰減率也大,但軌道衰減率的差異更大,所以實體-鐵木辛柯梁模型的聲功率級更小。

4.2 式(2)與式(3)的聲輻射對比

根據式(2)和式(3)得到的自由場鋼軌輻射聲功率級如圖10所示。

圖10 鋼軌輻射聲功率級對比

由圖10可以看出:

(1)頻率低于150 Hz時,式(3)得到的鋼軌聲功率小于式(2)的聲功率,這是由于在該頻段,鋼軌的衰減率較大,鋼軌趨近于偶極子聲源,與式(2)線聲源的假設不符。

(2)高于150 Hz的大部分頻率,式(3)得到的鋼軌聲功率大于式(2)的聲功率,最大高約2 dB,這是由于式(2)的推導是在連續彈性支撐梁的基礎上進行的,而在扣件等效剛度相等時,連續彈性支撐梁的聲輻射會低于離散支撐梁。差異的峰值出現在第一、第二階pinned-pinned頻率及2 000 Hz附近。將該區域放大可以發現,兩種方法得到的輻射聲功率頻率峰值存在差異,式(2)的峰值在1 840 Hz,式(3)的峰值出現在1 930 Hz,差異出現的原因在于式(2)的方法忽略了聲在軌道方向傳播的波動性,該波動性對聲輻射的影響將在4.3節中進行詳細闡述。

(3)兩個公式均預測了pinned-pinned頻率附近的峰值。

4.3 剛性軌道板對鋼軌輻射聲功率的影響

在實際情況中,鋼軌作為聲源并不存在于自由聲場中。對于高速鐵路無砟軌道而言,鋼軌位于軌道板上方,因此有必要考慮鋼軌位于軌道板上方時的輻射聲功率,由于軌道板表面阻抗很大,可以認為聲波在軌道板表面發生了全反射,通常可以構造鏡像聲源來模擬剛性反射面,如圖11所示。用式(3)可以得到考慮軌道板反射時的鋼軌輻射聲功率,其中的聲壓是考慮了源鋼軌與鏡像鋼軌的相互影響后得到的聲壓。

圖11 鋼軌及鋼軌鏡像

目前的板式無砟軌道中,鋼軌和軌道的間距常常處于40~100 mm之間,所以本文計算了間距分別為40、60、80 mm三種工況的鋼軌輻射聲功率,如圖12所示。

圖12 剛性軌道板上鋼軌的輻射聲功率

從圖12可以看出:相比于自由場下的鋼軌輻射聲功率,考慮軌道板的反射后最顯著的特征是會出現一些峰值,對應于40、60、80 mm三種工況,峰值分別出現在830、4 470 Hz,790、3 060 Hz,730、2 320 Hz。隨著距離的增加,峰值頻率會減小,峰值的幅值也會降低。

為了探明這些峰值產生的原因,分別畫出不同激振頻率和不同軌道方向波數時單位力激勵的鋼軌振動速度譜,如圖13所示,圖13中的色彩代表力作用點處鋼軌的振動速度水平,粗黑線是自由鋼軌的頻散曲線,以及以單位速度振動時鋼軌表面平均聲壓譜(即鋼軌同一截面上,表面各點聲壓的平均值,鋼軌不同截面振動速度相同因此聲壓分布也相同),自由空間和剛性軌道板上方鋼軌的表面平均聲壓如圖14和圖15所示。圖13~圖15中顏色越深紅表示振動速度/聲壓的幅值越大。

圖13 無限長周期離散支撐鋼軌的速度譜

圖14 自由空間鋼軌的表面平均聲壓

圖15 剛性軌道板上方鋼軌(間距為60 mm)的表面平均聲壓

由圖13可以看出,鋼軌的振動能量主要由彎曲波提供,但由于離散支撐的存在,頻率在500 Hz以下的彎曲波與自由鋼軌存在一定差異。易知隨著扣件剛度的增加,該頻率也會隨之增加,170 Hz是一個cut-on頻率。

對于聲壓譜,由圖14可以看到,自由場中的鋼軌聲壓在2 000 Hz附近存在一條亮帶,這條亮帶是由鋼軌本身的聲學邊界形狀所決定的,波數為0時,亮帶的中間頻率為1 840 Hz,這一頻率正是式(2)計算得到的一個鋼軌聲功率峰值,但式(3)得到的聲功率峰值出現在這條亮帶與鋼軌彎曲波頻率曲線的交點上,這一交點的對應頻率為1 930 Hz;由圖15可以看到,具有反射面聲場中的鋼軌聲壓相比于自由場鋼軌聲壓新增了兩條亮帶,這兩條亮帶由鋼軌和軌道板的間距所決定,可以看到兩條亮帶與鋼軌頻散曲線的交點(790、3 060 Hz)正是圖12中的鋼軌輻射噪聲峰值頻率。值得注意的是,鋼軌的聲壓由式(5)決定。鋼軌的聲壓譜為何會產生峰值,筆者認為可能由特征波導致,但鋼軌振動速度為0時,其聲壓譜的峰值在圖14和圖15中并不相同。對此筆者還在進一步研究。

5 結論

(1) 使用有限長鋼軌模型計算振動特性,端部導納和原點導納比值最大不超過20%時,其原點導納與無限長鋼軌模型基本一致,不會出現因截斷邊界產生的波動,此時可以將有限長鋼軌近似為無限長。

(2) 使用MPC連接實體模型和鐵木辛柯梁鋼軌模擬鋼軌,在保留實體鋼軌模型的優點同時,可以大量減少計算消耗。當網格劃分合適,實體-梁鋼軌模型的動力響應結果在高于3 000 Hz的頻率比梁單元鋼軌的可信度更高。無限長鐵木辛柯梁模型相比另外兩種模型計算更為方便,可以集成軌道模型和聲學邊界元模型建立軌道系統噪聲預測軟件。

(3) 由于實體-梁鋼軌模型在高頻會產生較多的變形,所以導納和軌道衰減率會比梁鋼軌模型更大,這也會導致實體-鐵木辛柯梁組合模型在高頻段的輻射聲功率比鐵木辛柯梁模型小了很多,而且這個差異隨著頻率的增加而增大。

(4) 由于多普勒效應,移動諧載荷下的軌道衰減率在pinned-pinned頻率處會出現兩個峰值,同時荷載前的鋼軌衰減率要大于荷載后的鋼軌衰減率。

(5) 相比于基于周期離散支撐鋼軌振動和2.5維邊界元計算鋼軌聲功率的方法,基于軌道衰減率和2維邊界元法的鋼軌聲功率傳統方法的結果在150 Hz以下高估了鋼軌聲功率,在150 Hz以上低估了鋼軌聲功率。

(6) 剛性軌道板會使鋼軌輻射聲功率出現多個峰值,峰值頻率隨鋼軌與軌道距離的增加而減小,其原因是峰值頻率位于鋼軌的速度譜與聲壓譜亮帶的交點上,聲壓譜亮帶頻率隨距離的增加而降低。

鳴謝

在本文的研究過程中,多次與英國南安普敦大學聲與振動研究所(ISVR)的D.J.Thompson教授進行了討論,得到了他的大力指導。

猜你喜歡
振動模型
一半模型
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
噴水推進高速艇尾部振動響應分析
重要模型『一線三等角』
This “Singing Highway”plays music
重尾非線性自回歸模型自加權M-估計的漸近分布
振動攪拌 震動創新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動性
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 一级黄色网站在线免费看| 亚洲第一成年网| 久久人妻xunleige无码| 国产91蝌蚪窝| a毛片免费在线观看| 激情视频综合网| 午夜老司机永久免费看片| 在线看片中文字幕| 国产色伊人| 91精品免费高清在线| 色偷偷av男人的天堂不卡| 91亚洲影院| 成人午夜网址| 国产黄在线观看| 无码精品国产dvd在线观看9久| 欧美成a人片在线观看| 国产在线精彩视频二区| 午夜国产精品视频黄| 日本国产精品一区久久久| 久久综合九九亚洲一区| 国产二级毛片| 99在线国产| 四虎国产在线观看| 国产成人综合久久精品尤物| 国产爽妇精品| 一级全黄毛片| 日韩色图区| 毛片在线播放网址| 久久精品只有这里有| 日本高清在线看免费观看| 欧美日韩精品在线播放| 欧美激情第一欧美在线| 免费观看国产小粉嫩喷水| 国产SUV精品一区二区6| 99re经典视频在线| 99ri精品视频在线观看播放| 亚洲欧美日韩中文字幕一区二区三区 | 亚洲毛片在线看| 亚洲中文字幕在线观看| 国产一区二区影院| 天天做天天爱夜夜爽毛片毛片| 日本AⅤ精品一区二区三区日| 在线无码九区| 干中文字幕| 一级片免费网站| 成人夜夜嗨| 色偷偷一区二区三区| 97se亚洲综合不卡| 国内精品视频| 亚洲精品亚洲人成在线| 最新国产高清在线| 免费aa毛片| 97se亚洲| 91视频免费观看网站| 亚洲色图欧美激情| 国产欧美精品一区aⅴ影院| 无码视频国产精品一区二区| 日韩av无码DVD| 天堂av综合网| 无码一区二区波多野结衣播放搜索| 青青草原国产精品啪啪视频| 亚洲欧美日韩中文字幕一区二区三区 | 亚洲av无码牛牛影视在线二区| 亚洲精品成人7777在线观看| 亚洲AV无码不卡无码| 国产不卡在线看| 亚洲欧美h| 欧洲av毛片| 色播五月婷婷| 亚洲AV无码久久精品色欲| 国产欧美日韩精品综合在线| 小蝌蚪亚洲精品国产| 亚洲h视频在线| 中文字幕久久波多野结衣 | 国产真实乱了在线播放| 99这里只有精品免费视频| 91精品国产无线乱码在线| 美女啪啪无遮挡| 日韩黄色大片免费看| 色婷婷国产精品视频| 性欧美精品xxxx| 国产91蝌蚪窝|