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

地鐵小半徑曲線鋼彈簧浮置板軌道鋼軌波磨成因分析

2024-01-11 15:41:00王志強雷震宇
中國機械工程 2023年24期
關鍵詞:振動

王志強 雷震宇

同濟大學鐵道與城市軌道交通研究院,上海,201804

0 引言

城市軌道交通的不斷發展使得人們對地鐵線路振動和噪聲水平的上限要求越來越嚴格。目前,國內許多地鐵線路采用了各種各樣的新型軌道來實現減振降噪的目的,諸如套靴軌枕軌道、梯形軌枕軌道和浮置板軌道等。雖然改變軌道結構形式在一定程度上緩解了振動噪聲問題,但現場調研卻發現減振軌道上出現了廣泛的鋼軌波磨現象[1-3]。鋼軌波磨發生在鋼軌表面,呈現為周期性的波狀磨耗。鋼軌波磨的存在容易誘發輪軌系統產生高頻振動和惱人噪聲,這有悖于減振軌道的最初目的,并且還會造成車輛軌道部件的疲勞失效[4]。如何有效預防鋼軌波磨一直都是鐵路行業的重要研究課題,而理解鋼軌波磨的成因機制則是實現控制鋼軌波磨的首要前提。

針對鋼軌波磨形成原因,國內外研究人員從輪軌系統的不同角度展開了探索。陳光雄等[5]提出了輪軌系統摩擦自激振動誘發鋼軌波磨的觀點,并結合線路測試驗證了上述觀點的有效性[6-7]。之后,CUI等[8]進一步拓展了摩擦自激振動理論,并研究了摩擦自激振動和軌面不平順反饋振動的聯合作用對鋼軌波磨的影響。結合現場測量和數值仿真,LIU等[9-10]發現曲線軌道上產生的波磨可歸因于摩擦誘發的扭轉振動引起的摩擦功率的周期性波動,并且受扣件影響的垂向動力學決定了扭轉振動的頻率,進而造成了實測線路上波磨波長隨扣件類型變化的現象。摩擦自激觀點能夠解釋現有地鐵系統中的大部分鋼軌波磨現象,并且能夠精確預測鋼軌波磨的波長和位置,但在表征波磨演化的輪軌微觀接觸方面還需進一步補充和完善。鋼軌波磨種類繁多,影響因素復雜,因而使用一種理論并不能完全描述其形成機制。早在1993年,根據波長固定機理和損傷機理,GRASSIE等[11]將當時鐵路系統中的鋼軌波磨現象劃分為六類,分別為pinned-pinned共振、車轍、其他P2共振、重載、輕軌和特定軌道形式波磨,并分別給出了相應的治理對策。基于此,大量針對性的案例研究在波磨領域涌現。通過建立輪軌非穩態滾動接觸有限元模型,YU等[12]研究了輪軌接觸對短波波磨的影響以及列車通過頻率對系統動力學的影響,結果發現當列車通過頻率接近鋼軌pinned-pinned頻率時,pinned-pinned共振會發生,并將促進扣件附近波谷處的滑動磨耗。關慶華等[13]詳細研究了車輛軌道系統P2共振的影響因素,并結合試驗提出和驗證了通過軌道振動屬性反推P2共振頻率的方法。金鋒等[14]對神朔重載鐵路波磨地段進行了長期跟蹤測試,發現其成因與有砟軌道P2共振頻率有關,并將該線路上的波磨發展總結為初期萌生、中期快速發展和末期穩定三個階段。此外,SUN等[15-16]、WANG等[17]、堯輝明等[18]從輪軌黏滑角度探究了波磨的形成機制,認為輪軌微觀接觸界面的不穩定黏滑振動促使了波磨的生成和發展。黏滑理論結合了微觀接觸和宏觀振動兩個尺度,能夠更為全面且直觀地表征鋼軌材料磨耗和系統振動特性,因而更有利于描述鋼軌波磨的發生和發展過程。現有文獻中關于鋼軌波磨黏滑振動的研究主要基于半空間理論和穩態假設,且未考慮輪軌瞬態滾動接觸特性[19],鑒于此,本文采用有限元方法研究鋼軌波磨的形成原因,該方法能夠詳盡地考慮輪軌系統動力學、接觸力學與表面損傷之間的相互作用,這對解釋鋼軌波磨現象至關重要。

目前地鐵系統中鋼彈簧浮置板軌道使用較為廣泛,且一般用于減振要求較高的特殊地鐵區段。鋼彈簧浮置板軌道主要通過浮置板受力后的慣性作用實現振動傳遞的衰減,從而減少浮置板傳遞到道床的振動能量。雖然鋼彈簧浮置板軌道的采用有效降低了輪軌系統周圍環境的振動水平,但它似乎對輪軌界面并不友好,大量的現場調研發現鋼彈簧浮置板軌道上出現了嚴重的鋼軌波磨現象,尤其是在小半徑曲線區間。為了闡明小半徑曲線鋼彈簧浮置板軌道上的波磨成因,本文采用有限元方法,從輪軌黏滑振動角度描述鋼軌波磨的演化過程。首先,根據現場線路情況,建立輪對-浮置板軌道三維有限元模型;然后,利用上述模型分析輪軌接觸黏滑特征和鋼軌縱向磨耗特性,從而量化波磨的發展趨勢;最后,借助復模態理論,研究輪軌系統的固有振動屬性與波磨發展的關系,以期理解浮置板軌道上的波磨機制。

1 模型的建立

1.1 波磨調研

實測波磨區間位于天津某地鐵小半徑曲線線路,曲線半徑為350 m,外軌超高為95 mm,軌距為1435 mm,鋼軌類型為CHN60,軌下結構為鋼彈簧浮置板,其上部通過ZX-2扣件與鋼軌相連,下部則直接與道床連接。線路運營車輛為地鐵B型車,運行速度為60 km/h。現場調研發現,小半徑曲線區間內鋼軌波磨較為嚴重,且主要表現為短波波磨,波長約為25 mm。鋼軌波磨現場照片見圖1。

(a)整體圖

1.2 有限元模型

參照波磨區間線路情況,利用有限元軟件ABAQUS建立了輪對-浮置板小半徑曲線軌道三維實體模型,如圖2所示,主要包括輪對、鋼軌、ZX-2扣件、浮置板、鋼彈簧隔振器和道床,其中扣件連接鋼軌和浮置板,鋼彈簧隔振器連接浮置板與道床。有限元模型中軌道線型為圓曲線,曲線半徑為350 m,外軌超高參考實際線路設為95 mm,軌距為1435 mm,軌道長度為25 m,其中鋼軌和道床沿模型全長布置,而浮置板共計4塊,每塊長度為6 m并沿軌道中心線連續布置。車輪型面為LM磨耗型,鋼軌型面為CHN60,外輪外軌接觸角設為18.6°,內輪內軌接觸角設為1.50°[20]。扣件間距為0.6 m,鋼彈簧隔振器間距為1.2 m,道床底部均勻布置接地彈簧與地基相連,以上連接部件均通過彈簧阻尼單元進行模擬,能夠考慮三個方向上的剛度和阻尼特性。模擬各部件的實體結構均采用C3D8R單元進行離散,輪軌接觸區域網格尺寸設置為1 mm,其余輪軌區域網格尺寸逐漸過渡至30 mm,浮置板和道床結構網格尺寸均勻設置為100 mm,整個模型共計279.4516萬個單元,325.3375萬個節點。有限元模型材料及結構參數取值見表1[21-22],其中材料屬性均考慮為線彈性。

表1 模型結構參數取值Tab.1 Values of model structural parameters

(a)等軸視圖

有限元模型中輪軌接觸部分采用Coulomb摩擦模型進行表征,法向接觸為硬接觸,切向接觸采用罰函數法描述,輪軌摩擦因數設為0.35[23]。模型中車輪平移速度為60 km/h,轉動速度為39.68 rad/s(車輪名義滾動圓半徑為420 mm),且輪對兩端軸箱位置處承受幅值為69 kN的垂向力[24],以模擬一系懸掛力。模型邊界條件設置如下:軌道結構縱向兩端面固定約束,浮置板和道床板橫向兩端面約束橫向位移,其他表面及結構無任何約束。

本文主要采用隱式積分方法分析輪軌瞬態滾動接觸行為,該方法在求解t+Δt時刻的控制方程時需要基于t和t+Δt兩個時刻的動態參量。因為上述時刻的動態參量隱含于控制方程,所以需要進行非線性計算。對于一般性的結構問題,隱式積分方法在能夠給出可接受解時,通常需要更多的時間步,即需要更長的計算時間。此外,結構響應預測精度將隨著時間步長Δt相對典型響應振型周期T的增大而減小,因此,在選擇最大允許時間步長時,應著重考慮如下三個方面:①施加荷載的變化速率;②非線性阻尼和剛度屬性的復雜性;③結構的典型振動周期。一般情況下,當Δt/T<1/10時,隱式積分方法能夠提供可靠的計算結果[25]。在ABAQUS軟件中,隱式積分方法的時間步能夠在殘差半增量[25]的基礎上進行自動選擇,當獲得t+Δt時刻的解后,通過對比平衡殘差在t+Δt/2時刻的值,便可以評估解的精確性并適當地調整時間步。

隱式積分方法的動態求解過程主要由位移和速度積分的Newmark公式定義:

u|t+Δt=

(1)

(2)

基于以上計算方法,下面將著重探究輪軌接觸黏滑與鋼軌縱向磨耗特征及輪軌系統固有振動屬性,以期解釋地鐵小半徑曲線鋼彈簧浮置板軌道上的鋼軌波磨形成原因。

2 接觸黏滑與磨耗特征分析

2.1 輪軌接觸黏滑特性

通過運用輪對-浮置板小曲線軌道有限元模型并計算,可以獲得輪對運行過程中各個參量的動態響應。模型中輪對運行距離為6.25 m,即運行時間為0.375(6.25/60×3.6)s,其中包括1 m(5.25~6.25 m)的輪軌接觸求解區域;計算分析步長設為0.02 ms。上述輪對運行距離、輪軌接觸求解區域長度和計算分析步長參數的選擇主要依據文獻[26-27]。由于輪軌接觸求解區域的網格尺寸對模型數值解的精度具有重要影響,因此,在分析輪軌接觸黏滑特性之前,本節評估了網格尺寸大小的合理性。對于輪軌接觸求解區域,分別選擇大小為0.25 mm、0.50 mm和1.00 mm三種網格進行相同工況計算,相應的結果(對應接觸求解區域鋼軌中間位置)匯總于表2。由表2可得,三種網格尺寸的計算結果誤差較小,均在10%以內,但計算時間卻隨著網格尺寸減小而劇增,因此,考慮到計算效率,選擇1.00 mm的網格尺寸是合理的。以上計算所使用的處理器為Inter(R) Xeon(R) W-2145 CPU @ 3.70 GHz,機帶RAM為256.0 GB,處理器個數為8。

表2 不同網格尺寸下的數值結果匯總Tab.2 Summary of numerical results for different mesh sizes

輪軌接觸狀態是表征黏滑特性的重要參數,因此,本節主要從輪軌接觸狀態變化出發分析內外側輪軌界面的黏滑特性。輪對在接觸求解區域穩態運行過程中5個代表性時刻(代表一個完整的接觸狀態變化過程),即0.350 28 s(時刻1)、0.350 70 s(時刻2)、0.351 12 s(時刻3)、0.351 54 s(時刻4)和0.351 96 s(時刻5)對應的內外側輪軌接觸狀態如圖3、圖4所示。需要補充的是,選擇上述時刻跨度的依據是其表征了內外軌接觸狀態的最小變化周期。

(a)時刻1 (b)時刻2 (c)時刻3

由圖3、圖4可得,在輪對運行過程中,內側輪軌接觸狀態呈現周期性變化,即輪軌黏滑分布呈現周期性變化(部分黏著-部分黏著-滑移-部分黏著-部分黏著),說明內側輪軌系統發生了周期性的黏滑振動,進而可能導致內軌表面出現周期性的波狀磨耗,最終形成鋼軌波磨;外側輪軌接觸狀態總體而言未有明顯變化,表明外側輪軌系統不會出現黏滑振動,從而不易形成鋼軌波磨。由于鋼軌波磨與微觀材料磨耗密切相關,因此,為進一步量化鋼軌波磨的形成過程,以下將對鋼軌縱向磨耗特征進行分析。

2.2 鋼軌縱向磨耗特征

根據Archard磨耗模型[28]

(3)

可以得到磨耗體積Vw與法向載荷FN和滑移距離S成正比。其中,k為磨耗比例系數,其大小與輪軌法向接觸應力及相對滑移速度有關;H為兩接觸物體中較軟材料的硬度。法向載荷決定了材料的磨耗程度,而滑移則決定了磨耗是否發生,即在黏著區,輪軌間沒有相對滑移,故不會產生磨耗,而在滑移區,由于節點之間存在相對滑移,故會造成材料的磨耗。然而實際上,黏著區會存在黏著磨損,且黏著磨損的存在也會對鋼軌波磨的形成產生影響,因此,為使計算結果更為精確,本文在建模過程中指定切向接觸罰函數法中的最大彈性滑移為表面特征尺寸的0.005%,以考慮黏著磨損的作用。鑒于此,本節主要根據鋼軌滑移云圖分析鋼軌縱向磨耗特征,內外側鋼軌的滑移云圖見圖5、圖6。

(a)鋼軌縱向

對比圖5a和圖6a可發現,在鋼軌縱向上,內側鋼軌滑移量明顯大于外側鋼軌滑移量,表明內側鋼軌縱向磨耗程度更大。同時,通過對比圖5b和圖6b可以發現,在鋼軌橫向上,內外側鋼軌滑移量相差不大,表明內外側鋼軌橫向磨耗程度相近。由于鋼軌磨耗的發生是縱向和橫向滑移共同作用的結果,因此,通過提取軌面接觸區域節點縱橫向滑移量,借助數學軟件MATLAB,可繪制縱向和橫向滑移綜合作用下的滑移云圖,如圖7a所示。根據式(3)計算得到的鋼軌磨耗曲線如圖7b所示。其中滑移距離S包括黏著區中的彈性滑移(對應黏著磨損),k和H的取值參考文獻[29]。由圖7可得,內側鋼軌總滑移量稍大于外側鋼軌總滑移量,這是內外側鋼軌較大且相近的橫向滑移量的補償作用所致。此外,由圖5~圖7容易看出,內外側鋼軌表面縱橫向滑移云圖分布及磨耗曲線均表現出了一定程度的周期性,尤其是內軌縱向滑移云圖,其波長28 mm(圖8)與實測波磨波長25 mm相近,這與內側輪軌界面接觸黏滑狀態的周期特性相呼應,從而進一步說明了內軌更容易發生嚴重的周期性磨耗,即鋼軌波磨。上述分析表明輪軌系統動態響應始終存在周期特性,而周期特性如何形成還需進一步研究。

(a)鋼軌總滑移云圖

圖8 內軌縱向滑移局部云圖Fig.8 Local nephogram of longitudinal slip of inner rail

3 輪軌系統振動特性及參數敏感性分析

3.1 振動特性分析

通過運用復模態理論[5-6]并進行復特征值分析,提取了輪軌系統在0~1200 Hz范圍內的不穩定振型,如圖9、圖10所示。其中穩定/不穩定振型主要以復特征值實部的正負進行反映,如果復特征值實部大于0,則表明系統振型是不穩定的,且數值越大越不穩定,反之則相反。需要說明的是,輪軌接觸斑沿鋼軌縱向長度一般約為20 mm,由于接觸斑對波磨波長具有接觸過濾作用[30],所以在實測區間車輛運行速度為60 km/h的情況下,根據頻率計算公式

圖9 664.1 Hz對應的系統不穩定振型Fig.9 Unstable mode of the system corresponding to 664.1 Hz

圖10 665.8 Hz對應的系統不穩定振型Fig.10 Unstable mode of the system corresponding to 665.8 Hz

(4)

可得波磨波長對應的最高頻率為833.3 Hz,由此可知復特征值分析選取的頻率范圍0~1200 Hz是足夠的。式(4)中,ω為波磨通過頻率,Hz;v為車輛運行速度,km/h;λ為波磨波長,mm。

由圖9、圖10可知,在0~1200 Hz范圍內,輪軌系統共出現了兩個不穩定振型,分別對應頻率為664.1 Hz和665.8 Hz,振型均表現為內輪相對于軌道的彎曲振動,且兩個不穩定振型對應的復特征值實部分別為15.9和16.2,這說明內側輪軌系統在此頻率下發生了失穩。同時,由1.1節可得實測波磨區間車輛運行速度為60 km/h,內軌波磨波長約為25 mm,因此,由式(4)可得內軌波磨通過頻率約為666.7 Hz。因為復特征值分析所得不穩定振型對應頻率與內軌波磨通過頻率均相近,所以可以得出系統響應周期特性的形成機制為輪軌系統664.1 Hz和665.8 Hz對應的不穩定振型所致,且該振型均表現為內輪相對于軌道的彎曲振動。

綜合輪軌接觸黏滑與磨耗特征及輪軌系統振動特性分析結果,可以得出實測地鐵小半徑曲線浮置板軌道上的內軌波磨是由664.1 Hz和665.8 Hz對應的輪軌系統不穩定振型誘發的黏滑振動造成的,且該不穩定振型均表現為內輪相對于軌道的彎曲振動。

3.2 參數敏感性分析

由于鋼軌波磨是輪軌系統參數的一種集中體現,因此為確定不穩定黏滑振動的產生條件以實現鋼軌波磨的有效控制,本節分析了系統參數對不穩定振動的影響。選取扣件和鋼彈簧三個方向上的剛度和阻尼特性作為代表性輸入參數,各參數分布類型均為均勻分布,均值與表1中數值一致,最小/最大極限值如表3所示[31-33],其中Kfv、Kft、Kfl分別表示扣件垂向、橫向、縱向剛度;Dfv、Dft、Dfl分別表示扣件垂向、橫向、縱向阻尼;Ksv、Kst、Ksl分別表示鋼彈簧垂向、橫向、縱向剛度;Dsv、Dst、Dsl分別表示鋼彈簧垂向、橫向、縱向阻尼。基于3.1節中的數值模型并以復特征值實部作為分析目標,利用有限元隨機分析模塊并執行概率設計分析,可以獲得反映鋼軌波磨發生趨勢的不穩定振動的評價指標-復特征值實部的統計特性,以作為參數敏感性分析的基礎。

表3 輸入參數最小/最大極限值Tab.3 Minimum/maximum limit values of input parameters

采用Spearman秩相關系數rsj[34]間接表示參數的敏感性,其計算公式如下:

(5)

j=1,2,…,12 -1≤rsj≤1

式中,Pkj為第j個隨機參數的秩;Qkj為與Pkj對應的響應量的秩;n為總樣本數,本節取為800。

Spearman秩相關系數可用于度量變量之間聯系的強弱,如果計算所得系數接近1或-1,則認為輸入變量對輸出變量影響顯著;如果計算所得系數趨近0,則認為影響微弱。如果系數為正值,表明增大輸入變量的值,輸出變量的值也增大;如果系數為負值,表明增大輸入變量的值,輸出變量的值反而減小。

將復特征值實部對各系統參數的敏感度匯總于表4,其中ξi表示復特征值實部。由表4可得,對復特征值實部影響較大的前兩個參數分別為扣件垂向剛度和鋼彈簧垂向剛度,且兩者對應的敏感度為負,表明它與復特征值實部成負相關,這也可從圖11、圖12所示的復特征值實部隨扣件垂向剛度和鋼彈簧垂向剛度變化的散點圖中看出。由于減小復特征值實部可以有效抑制系統的不穩定振動,從而控制鋼軌波磨的發生,因此,適當地增大扣件和鋼彈簧的垂向剛度能夠對鋼軌波磨的控制起到積極作用,這也從系統參數角度說明鋼彈簧較低的垂向剛度特性誘發了實測線路上鋼軌波磨的形成(實測線路上ZX-2扣件為普通減振扣件,即非低剛度扣件)。其他參數對復特征值實部影響較小,說明這些參數對實測線路上的波磨形成影響不大。

表4 復特征值實部ξi的敏感度Tab.4 Sensitivities of real part ξi of complex eigenvalues

圖11 ξi隨Kfv變化散點圖Fig.11 Scatter diagrams of ξi varying with Kfv

圖12 ξi隨Ksv變化散點圖Fig.12 Scatter diagrams of ξi varying with Ksv

4 結論

為解釋地鐵小半徑曲線鋼彈簧浮置板軌道上的波磨成因,本文從輪軌黏滑振動角度詳細地闡述了鋼軌波磨的形成過程。首先,依據現場線路情況完成了輪對-浮置板軌道三維有限元模型的構建;然后利用上述模型分析了輪軌接觸黏滑特征和鋼軌縱向磨耗特性,從而量化了波磨的發生趨勢;最后,基于復模態理論研究了輪軌系統的固有振動屬性與波磨形成的關系。主要得到以下結論:

(1)在輪對運行階段,內側輪軌黏滑分布呈現周期性變化,說明內側輪軌系統發生了周期性的黏滑振動,進而可能導致內軌表面出現周期性的波狀磨耗,最終形成鋼軌波磨;外側輪軌黏滑分布總體上未有明顯變化,說明外側輪軌系統未出現黏滑振動,從而不易形成鋼軌波磨。

(2)在鋼軌縱向上,內軌表面滑移量明顯大于外軌表面滑移量,表明內軌縱向磨耗程度更大;在鋼軌橫向上,內外側鋼軌表面滑移量相差不大,表明內外側鋼軌橫向磨耗程度相近。

(3)內外側鋼軌表面縱橫向滑移云圖分布均表現出了一定程度的周期性,尤其是內軌縱向滑移云圖,其波磨波長為28 mm,與實測波磨波長25 mm相近,這與內側輪軌界面接觸黏滑狀態的周期特性相呼應,從而進一步說明了內軌更容易發生嚴重的周期性磨耗,即鋼軌波磨。

(4)復模態分析顯示,內側輪軌系統在與實測波磨通過頻率666.7 Hz相近頻率處發生了失穩,結合輪軌接觸黏滑與磨耗特征分析結果,可以得出實測地鐵小半徑曲線浮置板軌道上的內軌波磨是由664.1 Hz和665.8 Hz對應的輪軌系統不穩定振型誘發的黏滑振動造成的,且該不穩定振型均表現為內輪相對于軌道的彎曲振動。

(5)參數敏感性分析顯示扣件和鋼彈簧垂向剛度對復特征值實部影響較大,且與復特征值實部負相關。由于復特征值實部的減小能夠有效地抑制系統的不穩定振動,從而控制波磨的發生,因此,適當地增大扣件和鋼彈簧的垂向剛度能夠對波磨控制起到積極作用。

猜你喜歡
振動
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
某調相機振動異常診斷分析與處理
大電機技術(2022年5期)2022-11-17 08:12:48
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
This “Singing Highway”plays music
具非線性中立項的廣義Emden-Fowler微分方程的振動性
中立型Emden-Fowler微分方程的振動性
基于ANSYS的高速艇艉軸架軸系振動響應分析
船海工程(2015年4期)2016-01-05 15:53:26
主回路泵致聲振動分析
UF6振動激發態分子的振動-振動馳豫
計算物理(2014年2期)2014-03-11 17:01:44
帶有強迫項的高階差分方程解的振動性
主站蜘蛛池模板: 538国产在线| 日本精品中文字幕在线不卡 | 欧美综合中文字幕久久| 亚洲国产欧美国产综合久久 | 日本伊人色综合网| 无码免费视频| 亚洲av无码人妻| 免费看美女自慰的网站| 亚洲制服丝袜第一页| 国产精品免费电影| 熟妇人妻无乱码中文字幕真矢织江| 青草午夜精品视频在线观看| 国产性爱网站| 国产亚洲精| 成人亚洲国产| 97视频免费看| 国产网友愉拍精品视频| 激情乱人伦| 亚洲无线观看| 欧美午夜精品| 日韩国产精品无码一区二区三区| 国产黄在线观看| 成人福利在线免费观看| 在线国产91| 亚洲国产欧美自拍| 91午夜福利在线观看| 亚洲六月丁香六月婷婷蜜芽| 狠狠色香婷婷久久亚洲精品| 色亚洲成人| 久久狠狠色噜噜狠狠狠狠97视色| 国产精女同一区二区三区久| 丰满人妻一区二区三区视频| 在线毛片网站| 中文字幕亚洲精品2页| 欧美国产精品不卡在线观看| 亚洲色图在线观看| 免费无码AV片在线观看中文| 亚洲一区网站| 2021国产乱人伦在线播放| 91精品国产一区| 国产一区三区二区中文在线| 国产污视频在线观看| 亚洲乱码视频| 91精品国产综合久久香蕉922 | 日韩不卡高清视频| 成人va亚洲va欧美天堂| 国产一区自拍视频| 999国产精品永久免费视频精品久久| 人妻无码一区二区视频| 四虎影视库国产精品一区| 国产午夜看片| 欧美精品啪啪| 亚洲国产系列| 亚洲一级毛片在线观播放| 欧美成一级| WWW丫丫国产成人精品| 免费看美女毛片| 国产成人AV综合久久| 亚洲福利一区二区三区| 久久综合色天堂av| www欧美在线观看| 国产在线自揄拍揄视频网站| 日韩精品亚洲一区中文字幕| 嫩草在线视频| 成人一级黄色毛片| 蝴蝶伊人久久中文娱乐网| 欧美精品啪啪一区二区三区| 91福利免费| 一区二区三区四区精品视频| 丁香婷婷综合激情| 又爽又大又黄a级毛片在线视频 | 中国国产A一级毛片| 亚洲午夜国产片在线观看| 久久国产精品波多野结衣| 精品福利视频网| 国产在线一二三区| 国产乱人激情H在线观看| 超清无码一区二区三区| 无码福利日韩神码福利片| 国产特级毛片aaaaaaa高清| 欧美日韩资源| 无码福利日韩神码福利片|