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

基于周期結構法的彈性短軌枕軌道聲振特性分析

2020-10-17 03:28:24徐涆文劉曉龍梁樹林陳秉智
中國鐵道科學 2020年5期
關鍵詞:有限元振動結構

徐涆文,韓 健,劉曉龍,田 彩,梁樹林,陳秉智

(1.西南交通大學 牽引動力國家重點實驗室,四川 成都 610031;2.大連交通大學 機車車輛工程學院,遼寧 大連 116028)

隨著國內列車運行速度的不斷提高,以及城市軌道交通的迅猛發展,人們對于軌道的減振降噪性能要求也越來越高。軌道作為主要的噪聲和振動源之一,國內外研究人員對其振動和聲輻射特性也進行了大量研究[1-2]。大量的減振型軌道也隨之產生,并被鋪設于振動噪聲要求較高的路段,減振型無砟軌道由于其良好的減振降噪性能,運用在我國的地鐵、輕軌以及城際鐵路上。彈性短軌枕軌道即為減振型軌道的1 種。在軌枕下加入軌下彈性墊層,提高了彈性短軌枕軌道結構的減振和降噪性能,同時其具有結構簡單,施工方便,軌枕下支承剛度易于調整的優點。但是目前對于彈性短軌枕振動和聲輻射特性的了解較少,對相關參數選取認知還不夠深入。若參數選取不合理則可能無法發揮出其減振作用,甚至會增加軌道振動和聲輻射響應,導致車內噪聲增大。因此,分析彈性短軌枕軌道不同結構參數對軌道的振動聲輻射影響就顯得尤為重要。

楊俊斌[3]對2 種彈性支撐軌道進行了各個結構參數對結構穩定性以及高速鐵路上的適應性影響分析,確定了高速列車荷載作用下彈性支承軌道各結構參數的合理取值。方銳等[4-5]基于有限元—邊界元法建立了有砟軌道垂向振動和聲輻射有限元模型,分析了鋼軌、軌枕和道床等相關結構參數對軌道振動及其聲輻射的影響。王根平[6]基于FEM和BEM 方法對減振型雙塊式無砟軌道進行了垂向振動特性以及軌道各結構層聲輻射效率的分析。

受計算能力限制,在對周期結構進行動態特性分析時,往往需要對周期結構進行有限長度截取,為了消除截取的周期結構邊界波動反射對于仿真結果的影響,需要對截取界面選擇合理的邊界條件。劉晶波,谷音[7-8]基于矩陣等效原理采用有限單元來模擬二維以及三維黏彈性邊界單元,以消除邊界反射波干擾。而這種邊界結構僅適用于材料均勻連續的結構,并不適用于離散支撐的軌道結構。金浩[9]基于譜單元法建立了軌道的元胞模型從而能夠較好地表征軌道的“無限長”結構特性,利用模型分析了鋼軌動力吸振器對于軌道振動特性的影響。馬龍祥[10]利用軌道結構的周期結構特性建立了軌道周期結構模型,分析了軌道結構在移動載荷激勵下的振動特性。但是上述2 種軌道結構的鋼軌均采用的是經典歐拉梁模型,其不能較好地模擬出實際的鋼軌截面形狀。同樣各種商用有限元軟件中也提出了相應的無反射邊界單元,但是對于離散支撐的有限元軌道結構其并不能有效地衰減有限截取長度軌道邊界的反射波干擾。目前國內消除實體有限元軌道邊界反射波的方法多為加長結構計算模型。這達不到降低計算效率的建模初衷。

為了消除有限截取模型的邊界彈性波反射的影響,本文基于周期結構原理建立了無反射邊界的有限元軌道模型,分析彈性短軌枕軌道的隔振性能、軌道各結構的聲輻射特性,并研究不同軌枕下支承剛度對其振動和聲輻射特性的影響。

1 軌道周期結構模型

1.1 模型建立

周期結構是由元胞按一定的規律重復排列而形成的結構形式。當組成周期結構的元胞按規律在空間中是無限重復時即形成理想的周期結構[10]。軌道可看作為1種縱向方向的理想周期結構,其中每1 塊道床板長度的軌道結構即可看作周期結構中的元胞。圖1為運用有限元軟件ANSYS 建立的軌道結構元胞有限元模型。

圖1 軌道結構元胞有限元模型

在軌道結構元胞有限元模型中,鋼軌、軌枕和道床采用實體單元(solid185)進行離散,扣件系統采用線性的彈簧阻尼單元(combine14)進行模擬,其垂向剛度為80 MN·m-1,橫向剛度為50 MN·m-1,縱向剛度為40 MN·m-1。模型中,鋼軌為UIC60 軌,其彈性模量為210 GPa,泊松比為0.3,質量為60 kg·m-1;軌枕選用SK-1 型軌枕,其彈性模量為35 GPa,泊松比為0.2,質量為90 kg;道床板混凝土設計標號為C55,其彈性模量為36 GPa,泊松比0.2,密度為2 500 kg·m-3。

軌道結構元胞強迫振動動力學方程為[11]

其中,

式中:D為軌道結構元胞總動剛度矩陣;S為軌道結構元胞位移矩陣;F為軌道結構元胞受到的作用力向量;C,K和M分別為軌道結構元胞的阻尼矩陣、剛度矩陣和質量矩陣;ω為角頻率;SI,SL和SR分別為軌道結構元胞中鋼軌中部,左邊界,右邊界的節點位移;FI,FL和FR分別為軌道結構元胞中鋼軌中部,左邊界和右邊界節點受到的作用力;Dij為對應于元胞結構中部、左邊界和右邊界節點位移下的局部剛度矩陣,下標i,j=I,L,R。

由于軌道結構元胞的中部不受外激勵力的作用,即FI=0。因此可以通過等式SI=-(DILSL+DIRSR)對上述的動剛度矩陣D進行縮聚處理[11]

其中,

式中:D(1)為經過縮聚后的單個元胞動剛度矩陣。

下面將2 個軌道結構元胞進行組合,有限元模型如圖2所示。

圖2 元胞組合結構圖

其強迫振動方程為

此時2 個軌道結構元胞集合成1 個整體,對新生成的結構剛度矩陣進行縮聚得

其中,

因此對于具有周期性軌道結構,可以通過上述方法進行多個相同元胞疊加,得到理想長度的軌道周期結構總動剛度矩陣與初始激勵輸入和最終響應輸出的關系。最終疊加后的理想周期結構總動剛度矩陣D(n)維度大小同單塊道床板長度的軌道結構總動剛度矩陣相同,其中n=1,2,…,N,N為迭代的次數。因此將對應項系數進行替換即可通過單塊道床板長度的軌道結構模擬無限長軌道結構的振動衰減特性。

1.2 模型驗證

本文研究單位垂向簡諧力作用F下的鋼軌垂向振動固有特性,激勵力作用位置為鋼軌跨中的鋼軌頂面。為提高計算效率,根據軌道結構的橫向對稱特性,截取軌道結構橫向的一半模型作為計算模型。

軌道是1 種沿縱向方向的周期結構,鋼軌頂面在單位簡諧力的作用下產生的振動波沿著鋼軌縱向向兩側傳播,實際中的軌道結構可看作理想的周期結構,彈性波沿著鋼軌縱向傳播不斷衰減最終趨于零[12]。通過周期結構原理對單塊道床板長度的軌道模型的邊界進行處理以模擬左右兩側無限延長的軌道結構,結構如圖3所示。

提取單塊道床板長度的軌道有限元模型的動剛度矩陣,根據上一節的方法疊加出理想周期結構軌道的總動剛度,由于此總動剛度矩陣縮聚后矩陣維度同單塊軌道板的總動剛度矩陣相同,將多次疊加后的總動剛度矩陣內的對應項參數賦值給單塊道床板長度的有限元軌道結構,最后將其加到截取后的軌道有限元結構的邊界,從而模擬無限長軌道振動的縱向衰減特性,消除結構邊界反射波對于原點響應分析的影響。

圖3 軌道周期結構模型

為了驗證計算模型的準確性,對某地鐵非彈性軌枕軌道進行力錘敲擊試驗,獲取其頻率響應函數,與仿真計算結果進行對比,進而修正有限元模型。試驗中,鋼軌類型為UIC60,軌枕以及道床板均為混凝土結構,鋼軌扣件間距為0.625 m。

軌道結構較復雜,并且軌道頻響測試結果會受到較多因素影響,如扣件裝配位置精度,鋼軌,道床板以及扣件彈性墊等材料的實際屬性,測試環境溫度等,而仿真模型則是抓住主要的影響因素,采用簡化的模型模擬出軌道的實際振動特性。

圖4給出了采用3 塊和5 塊道床板長度有限元軌道模型及周期結構軌道模型計算得到的鋼軌垂向加速度導納。

圖4 不同軌道模型下的鋼軌垂向加速度導納

由圖4可知:加長軌道模型能夠較好地緩解邊界反射波對鋼軌頻響函數的干擾,但是過長的軌道模型會給仿真分析帶來極大的計算量,而建立的軌道模型周期結構邊界能夠有效地消除邊界反射波的影響,從而大大提升模型的計算效率。

圖5給出了鋼軌垂向加速度導納曲線。由圖5可知:2 條曲線整體符合程度較高,且曲線隨頻率變化的趨勢也是基本相同;特別是在230和1 020 Hz處仿真和實測結果均出現振動峰,并且振動峰的位置以及大小基本相同。由有限元軌道的模態分析可得這2種峰值對應的模態振型,分別為230 Hz的鋼軌垂向1階振動和1 020 Hz鋼軌的Pinned-Pinned振動,與實際的軌道振動模態也是完全符合。因此建立的模型能夠較為準確地反映軌道的振動特性,可以利用此模型對軌道的振動以及聲輻射特性進行分析。

圖5 周期結構邊界的軌道模型計算和實測的鋼軌垂向加速度導納

2 彈性短軌枕軌道聲振特性

彈性短軌枕軌道是1 種減振型軌道結構,同普通的非彈性軌枕軌道相比,彈性短軌枕軌道為了緩解軌道垂向和橫向的振動傳遞,在軌枕下方和軌枕四周加入了軌下彈性墊層和橡膠套結構[3],如圖6所示。

圖6 彈性短軌枕結構

為了分析在軌枕下加入軌下彈性墊層后對鋼軌振動以及聲輻射的影響,在上節驗證的剛性道床軌道模型的基礎上,在軌枕和道床板之間通過彈簧和阻尼單元進行連接,以模擬實際軌道中的橡膠套和軌下彈性墊層結構,研究軌枕下加入軌下彈性墊層后鋼軌的振動特性。根據文獻[3],取軌枕的相關結構參數為:軌枕彈性模量35 GPa,密度2 500 kg·m-3,單塊軌枕質量140 kg,泊松比0.2。軌枕的長、寬和高分別為700,320 和250 mm。扣件剛度同上述剛性道床軌道。軌下彈性墊層剛度為120 MN·m-1,軌枕橡膠套剛度為5 000 kN·m-2。彈性短軌枕有限元模型如圖7所示。

圖7 彈性短軌枕有限元模型

2.1 彈性短軌枕軌道振動響應特性

軌道結構的振動特性分析是其聲輻射分析的基礎,同理基于上述周期結構原理建立具有無反射邊界的彈性短軌枕軌道有限元模型,分析其在單位激勵力下的振動響應。

圖8為彈性短軌枕和剛性道床軌道上激勵力作用點處鋼軌的加速度導納。

圖8 彈性短軌枕軌道激勵力作用點處鋼軌加速度導納

由圖8可知:2 種軌道上鋼軌的加速度導納曲線除了在0~400 Hz 位置處振動峰有較大差異外,在Pinned-Pinned 振動峰以及其余頻率位置鋼軌加速度導納曲線并未受到軌枕下結構改變的影響;在頻率為0~400 Hz 的范圍內,彈性短軌枕軌道的鋼軌激勵點垂向加速度導納曲線出現2 個峰值,峰值頻率分別為124和270 Hz。

由軌道結構進行模態分析可知:彈性短軌枕軌道的第1 個峰值為鋼軌的1 階垂向彎曲振動模態,如圖9(a)所示,彈性短軌枕軌道中鋼軌1階垂向彎曲的峰值頻率較剛性道床軌道低,這是由于低頻的鋼軌振動主要受軌道垂向剛度影響,在軌枕下方加入軌下彈性墊層將降低軌道垂向的總剛度值,使得鋼軌的1 階彎曲變形頻率向低頻移動,幅值減小;頻率為270 Hz 的峰值則是由軌枕相對于鋼軌的反相振動引起的,如圖9(b)所示,這是因為在軌枕下加入軌下彈性墊層增加了軌道的垂向自由度,在對應頻率的簡諧單位力激勵下被激發出來。

圖9 彈性短軌枕軌道振型模態

彈性短軌枕軌道軌下結構的振動響應如圖10所示。由圖10可知:軌枕和道床板在124,270 和1 460 Hz頻率位置均出現了振動峰,由軌道模態可知,124 和270 Hz 處的振動峰分別由鋼軌的1 階垂向彎曲振動和鋼軌軌枕間的反向振動模態引起,而1 460 Hz處的振動峰則是由在簡諧激勵力作用下軌枕結構的1階垂向彎曲振動模態激發。

圖10 彈性短軌枕軌道各結構的位移導納

圖11給出了2 種軌道垂向振動的傳遞情況。由圖11可知:相對于彈性短軌枕軌道,剛性道床軌道的道床板垂向振動位移高出將近2 個數量級,而2 種軌道的鋼軌垂向振動差異較小,這表明彈性短軌枕軌道良好的垂向振動衰減能力。

圖11 彈性短軌枕軌道和剛性道床軌道的位移導納

為了更直接地對比2種軌道的減振性能,圖12給出了2 種軌道的隔振效率對比圖。由圖12可知:彈性短軌枕軌道的隔振效率在整個研究的頻率范圍內均高于剛性道床軌道,彈性短軌枕軌道的隔振效率均保持在1.000 0,表明彈性短軌枕軌道具有較好地隔振性能,并且隔振能力幾乎不會受到外激勵頻率的影響;相對彈性短軌枕軌道,剛性道床軌道隔振效率則要小,最大的隔振效率出現在0~140 Hz 的頻率范圍,為0.999 6,同時軌道的隔振效率會隨著外激勵頻率有較大變化,當外激勵頻率增加時,剛性道床軌道的隔振效率有降低的趨勢,但是在高頻有較大的激增波動,在頻率為1 020 Hz位置處,剛性道床軌道的隔振效率突然增大為0.998 0,此頻率對應的是鋼軌的垂向Pinned-Pinned 振動,表明鋼軌的Pinned-Pinned振動的垂向傳遞較差。

圖12 彈性短軌枕軌道和剛性道床板軌道隔振效率對比

可見,于剛性道床軌道相比,彈性短軌枕軌道具有較為優良的隔振性能。

2.2 彈性短軌枕軌道聲輻射特性

采用直接邊界元法計算彈性短軌枕軌道的聲輻射[4]。對于三維振動體的外部聲輻射問題,邊界單位法向量指向單元外部,當x點位于區域邊界上時,Helmholz積分方程為

式中:p(x)為邊界點聲壓;ρ為空氣密度;an(Q)為邊界區域任意點Q的法向加速度;G(Q,x)為自由空間的格林函數;A為邊界區域面積。

對邊界積分方程進行離散,可得邊界元求解方程

式中:H和L均為系數矩陣;p為表面點的壓力矩陣;an為模型表面點法向方向的速度矩陣。

本文通過指定模型表面上的節點法向速度(Neumann邊界條件)求解模型表面上的聲壓。

根據所求得的結構表面各點的聲壓,并通過結構的振動響應求得結構表面的法向速度,進而可求得軌道結構的輻射聲功率Wrad為

對于軌道聲輻射研究,進行軌道鋼軌和軌下結構聲輻射貢獻量分析對于噪聲的控制具有較大的價值。對聲功率級進行A 計權,從而能夠更好地模擬人耳對于聲音的感受程度。

基于聲學邊界元理論,在軌道有限元結構的表面進行單元劃分就可以得到所需要的聲學網格,基于上述邊界元理論可求得軌道結構的輻射聲功率,軌道結構的聲學網格模型如圖13所示。

圖13 軌道結構的聲學網格模型

圖14給出了單位垂向簡諧激勵力下的軌道整體和軌道各結構A 計權的聲功率級。計算頻率范圍為0~3 500 Hz,在軌道的垂向振動聲輻射分析中,在頻率低于240 Hz 時道床板和軌枕較鋼軌對噪聲的貢獻量大,隨著外激勵頻率增加,軌下結構總的貢獻量下降幅度較大,鋼軌此時成為主要的噪聲輸出源。從鋼軌的聲功率曲線能夠看出:在1 020 Hz 頻率位置出現了聲功率的極值,這是由鋼軌的Pinned-Pinned 振動引起的;軌下結構在1 460 Hz處也出現了聲功率級的峰值,這是由軌枕在激勵力下產生的1階垂向彎曲振動所致,模態振型如圖15所示。

由軌道各結構的聲功率級曲線可知:對軌道進行噪聲控制時,在低頻時控制鋼軌下結構的振動聲輻射,而在高于240 Hz 的頻率時,主要對鋼軌采取相應的降噪措施以降低軌道整體的振動聲輻射。

圖14 軌道各結構聲功率級(A計權)

圖15 軌枕垂向1階彎曲模態

2.3 軌下支承剛度對其聲振特性的影響

通過上述的分析可得,相對于剛性道床軌道,彈性短軌枕軌道在軌枕下加入了軌下彈性墊層大大提升了軌道的垂向減振性能,那么分析軌下彈性墊層剛度對于軌道聲振特性的影響就顯得尤為重要。因此下文對不同軌枕支承剛度下的軌道振動和聲輻射特性進行分析。

軌枕下支承剛度需要根據實際列車運行情況選取合理的參數。對軌枕下支承剛度分別為40,160和200 MN·m-1時的彈性短軌枕軌道振動和聲輻射特性進行分析。

圖16給出了不同軌枕支承剛度下鋼軌和軌下結構的垂向位移導納。由圖16可知:軌枕支承剛度主要影響400 Hz 以下的鋼軌振動,即剛度控制區;當軌枕下剛度為40 MN·m-1時,軌道的垂向剛度較低,鋼軌在外激勵作用下振動幅度較大,即能量更多地停留在鋼軌和軌枕結構上,道床板的振動較小,此時彈性短軌枕軌道減振性能較好,但是較大的鋼軌振動位移使得軌道的穩定性較差;隨著軌枕下支承剛度增加,軌枕與道床板間的耦合作用增加,鋼軌上的振動能量更多地傳遞到軌枕下結構,振動能量衰減效果較好,因此加大軌枕下剛度能提升鋼軌的垂向穩定性;由于鋼軌的振動能量更多地傳遞到軌下結構,道床板的垂向振動位移也隨之增加,軌道的減振性能降低。因此軌枕下墊板垂向剛度的選取需要在保證鋼軌穩定性的同時兼顧軌道的減振性能。

圖16 軌道垂向振動位移導納

對軌道的各部分在不同軌枕支承剛度下的總聲功率進行分析,結果見表1。

表1 軌道各結構聲功率級

對軌道整體結構進行總聲功率級分析可知:在0~3 500 Hz 頻帶范圍內鋼軌的總聲功率級高于軌下結構的總聲功率級,為軌道整體結構聲輻射的主要貢獻源;隨著軌枕下支承剛度增加,鋼軌在0~3 500 Hz的聲功率級有微幅增加,從62.304 5 dB(A)增加到62.355 3 dB(A),這是由于隨著軌枕下剛度增加鋼軌1階垂向彎曲峰向高頻移動,同時鋼軌的聲輻射效率隨著頻率增加而增加,因此鋼軌在1階垂向彎曲振動模態下會輻射出更多的噪聲;軌枕下支承剛度增加使得振動能量更多地傳到軌下結構上,導致軌下結構總聲功率級增大;當軌枕下支承剛度從40 MN·m-1增加到200 MN·m-1時,軌下結構的總聲功率級增加約3 dB(A);隨著軌枕下支承剛度的增加,軌道整體的總聲功率級有微幅增高,從62.926 dB(A)升高到62.928 dB(A),表明降低軌枕支承剛度對降低軌道整體的噪聲輻射作用小。

圖17 彈性短軌枕軌道輻射聲功率

圖17給出了鋼軌、軌下結構和軌道整體結構的聲功率級。由圖17可見:軌枕下支承剛度的改變也僅對結構低頻位置的聲功率級產生影響,當軌枕下剛度改變時,軌下結構的輻射聲功率受到剛度改變影響較大,對應于軌道1階彎曲振動位置的聲功率級受軌枕下支承剛度影響較大,因此在增加軌枕下剛度時軌下結構聲功率級增加較大;對于鋼軌而言,由于鋼軌顯著的聲功率峰值出現在較高的頻率段,軌枕下支承剛度的改變對整體的聲輻射影響也就相對較小,因此對于鋼軌在0~3 500 Hz 的頻段內的聲功率級影響就較小。

通過振動分析可知:隨著軌枕下支承剛度增加,鋼軌的振動向軌下結構傳遞效率加大,使得鋼軌振動減小,然而軌下結構振動加劇降低了軌道的減振性能。從軌道聲輻射來看,隨著枕下支承剛度增加鋼軌的垂向彎曲振動峰向高頻移動,而隨著頻率增加鋼軌的輻射效率增加,最終也將導致鋼軌的聲功率級增加。從軌道的整體結構聲功率級來看,降低軌枕下支承剛度能夠降低軌道整體結構的聲功率級,因此軌下剛度需要選取適中的參數同時保證良好的減振和降噪性能。

3 結 論

(1)具有周期結構邊界的軌道模型能夠較好地模擬軌道的振動特性,消除邊界反射波對激勵力作用點鋼軌響應分析的干擾。

(2)對于彈性短軌枕軌道的聲輻射而言,鋼軌是主要的噪聲輻射源,尤其在高頻位置。而軌枕和道床板結構在低頻位置的聲貢獻量也是不可忽略的。

(3)軌枕下支承剛度主要對軌道結構0~400 Hz 范圍內的振動聲輻射有影響。減小軌枕下支承剛度能夠提升軌道的減振和降噪性能,但是過小的垂向剛度將導致鋼軌穩定性降低,增加輪軌磨耗。因此綜合考慮軌道各項性能后,軌枕下支承剛度應當在40~160 MN·m-1范圍內進行合理選取。

猜你喜歡
有限元振動結構
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
中立型Emden-Fowler微分方程的振動性
論《日出》的結構
創新治理結構促進中小企業持續成長
現代企業(2015年9期)2015-02-28 18:56:50
磨削淬硬殘余應力的有限元分析
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: 国产一在线观看| 欧美精品一二三区| 婷婷伊人五月| 精品国产成人a在线观看| 国产黄色免费看| 成人亚洲视频| 国产丝袜第一页| 内射人妻无码色AV天堂| 国产三级国产精品国产普男人| 天天色天天综合| 狠狠色噜噜狠狠狠狠奇米777| 五月婷婷激情四射| 欧美a网站| 一级成人欧美一区在线观看| 中文字幕亚洲精品2页| 一级毛片在线播放免费| 日韩a级毛片| 国产精品99在线观看| 不卡国产视频第一页| 亚洲三级成人| 美女视频黄又黄又免费高清| 成人日韩精品| 精品国产自在在线在线观看| 婷婷开心中文字幕| 五月天综合婷婷| 欧洲av毛片| 午夜国产精品视频| 久久国产精品影院| 久久永久免费人妻精品| 2024av在线无码中文最新| 欧美午夜视频在线| 亚洲天堂成人在线观看| 国产三级成人| 国产91无码福利在线| 五月激情综合网| 国产原创自拍不卡第一页| 国产精品性| 国产精品污视频| 韩国自拍偷自拍亚洲精品| 国产又大又粗又猛又爽的视频| 日韩 欧美 小说 综合网 另类| 久久综合九九亚洲一区 | 手机精品视频在线观看免费| 青青操视频在线| 欧美视频在线第一页| 国产成人精品优优av| 国产欧美日韩视频怡春院| 国产精品美乳| 国产玖玖视频| 色综合久久88色综合天天提莫 | 无码又爽又刺激的高潮视频| 国产精品亚洲欧美日韩久久| 久久一本日韩精品中文字幕屁孩| 91麻豆精品国产91久久久久| 毛片免费在线视频| 免费看美女自慰的网站| 丁香婷婷激情网| 青草娱乐极品免费视频| 极品私人尤物在线精品首页| 四虎免费视频网站| 亚洲精品国产综合99| 日韩精品毛片| 毛片免费视频| 美女国内精品自产拍在线播放| 漂亮人妻被中出中文字幕久久| a色毛片免费视频| 亚洲男人的天堂网| 蝌蚪国产精品视频第一页| 亚洲AV无码乱码在线观看裸奔 | 三上悠亚一区二区| 少妇精品网站| AV不卡在线永久免费观看| 欧美精品二区| 国产人成在线视频| 福利在线一区| 久久久久久久久亚洲精品| 亚洲h视频在线| 亚洲色图欧美在线| 久久精品国产精品国产一区| 亚洲视频四区| 欧美成人区| 亚洲系列无码专区偷窥无码|