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

基于Kriging代理模型的城軌車輪踏面外形優化

2018-07-04 06:27:04施振宇邢宗義
鐵道標準設計 2018年7期
關鍵詞:優化模型

施振宇,戴 碩,邢宗義

(南京理工大學自動化學院,南京 210094)

為了減輕交通壓力以及降低環境污染,城市軌道交通的發展開始加速。城市軌道交通載客量大、安全準點,能夠滿足人民群眾的日常出行,而且具備節約土地資源以及環保的優勢,與人民群眾生產生活密切相關[1]。地鐵的運行線路與普通鐵路不同,彎道較多且曲線軌道半徑小,需要經常性的啟動和制動,造成輪軌接觸關系分析變得困難起來。隨著累計行程的逐漸增加,輪軌之間的損耗加大,造成輪軌系統性能變差,從而危害車輛運行的安全行駛,極大阻礙了城市軌道交通的發展[2]。城軌車輪磨耗除了使得輪緣輪徑的尺寸減小之外,同時踏面外形輪廓也會因此發生較大變化。部分地鐵線路運行一段距離之后輪緣發生磨損,出現踏面剝離、擦傷和平輪等形式的損傷[3]。而車輪踏面磨耗嚴重也會導致輪對壽命遠低于正常的使用周期[4]。踏面剝離是指由于持續的運營導致輪軌接觸疲勞,致使踏面表面出現剝落,產生裂紋,此時須進行鏇修,若出現磨耗超限則要及時更換,這將會產生一定的費用并大大減少輪對使用壽命[5]。輪軌外型不匹配將導致列車平穩性與曲線通過性能差、較大的輪軌接觸應力、踏面嚴重磨損和出現裂紋、運行中噪聲等一系列的問題,嚴重時甚至會引起列車脫軌等重大交通事故[6]。因此采用恰當的車輪踏面外形,可以減少輪軌之間的磨損,節約成本和減少維修次數,提高車輪的穩定性,增加使用壽命,改善車輛的動力學性能。

常崇義等[7]采用ABAQUS建立了三維動態輪軌接觸有限元模型,以現場試驗得到的實測數據對數值分析模型和方法進行了驗證。BRAGHIN等[8]應用基于簡化FASTSIM算法對輪軌關系進行了計算分析,研究出當車輛行駛超過20萬km再進行鏇輪作業,能夠將輪對的使用壽命延長一倍以上。員華等[9]對輪對鏇修工作現行采用的等級鏇修制度進行了分析,選擇不同的輪緣厚恢復等級進行鏇修,盡可能減少輪徑的鏇削量。PASCUAL等[10]分析研究城軌車輪輪對磨耗實測歷史數據,認為當輪緣厚度下限為27.5 mm,鏇修后的輪緣厚度值為30.5 mm時產生的輪對鏇修總費用是最小的。許宏等[11]基于高斯過程建立了輪對磨耗模型,采用蒙特卡羅仿真進行求解,實現鏇修參數的優化。王珍[12]通過對動車組歷史數據的分析,確定了單個車輪的最優鏇修參數選擇。Jendel[13]將車輛軌道動力學仿真計算放在時域進行,并在Gensys軟件中進行模擬。

輪對的鏇修工作主要分為以下兩個步驟:一是需要確定輪緣厚的恢復閾值和輪徑的切削量,即鏇床要按照鏇修后的車輪輪緣厚值進行鏇修加工,相應對輪徑進行切削以降低踏面測量的基準點;第二,在輪對的輪緣輪徑參數確定之后,還需要保證鏇修后的輪對能夠保持一個平滑的合理的外形輪廓,使得列車在正常運行過程中,滿足車輛的一系列動力學性能。本文主要對第二個步驟進行重點研究。

本文圍繞城軌車輪鏇修后的踏面外形輪廓展開研究,提出了基于Kriging[14-15]代理模型的單目標優化方案和多目標優化方案,對改進之后的結果與標準踏面的性能進行對比,研究動力學性能的各項指標,并比較兩種模型的優劣。結果表明多目標優化相比單目標而言可以同時對多個相互影響和制約的目標進行優化,在多個優化目標之間尋找一個折中點,使系統的綜合性能得到提高。

1 城軌車輪踏面外形Kriging代理模型變量

Kriging代理模型是在分析設計變量在空間上的相關性的基礎上,根據已知樣本信息來構建輸入變量與輸出變量之間的映射關系,從而預測某一未知待測點的信息的。

Kriging模型是對線性回歸分析的一種改進,Kriging模型需要利用在某一點為中心的鄰域內的已知設計點的信息對待測樣本點進行預測,即利用線性加權法將該點的鄰域內樣本設計點的信息加權來估計該點未知信息。利用獲得信息的最小誤差方差來確定加權的權系數。因此,Kriging模型本質上是最優的線性無偏估計。

Kriging模型的組成:多項式和隨機分布

y(x)=fT(x)β+Z(x)

(1)

式中,fT(x)=[f1(x1),f2(x2),…,fq(xm)]T為已知的多項式回歸模型;β=[β1,β2,…,βq]T為相對應的待估參數;Z(x)為一隨機過程,具有如下統計特征

(2)

式中,xi和xj是訓練樣本中的兩個點;R(θ,xi,xj)是相關性系數,包含待定參數θ,用于表示兩個樣本點之間的空間相關性。因此,Kriging代理模型將任意的響應值都視為一個服從正態分布隨機變量,如此模型具有更好的靈活性,而不限定于某種特定的形式。

2009年10月至12月,廣州地鐵進行了較為集中的鏇輪作業,先后實測結果如表1所示。可以看出,城軌車輪在經過鏇修之后再次上線后,車輪的磨耗以及鋼軌的磨耗會加劇。輪對鏇修會破壞已經進入穩定磨損狀態的輪軌接觸關系,輪軌將會重新回到運行初期的磨合階段,磨耗量平均值增大。因此鏇修作業后車輪保持一個具備良好的輪軌接觸關系的外形十分重要。現有的車輪型面的優化,主要是對踏面外形的優化,本文在TBT449—2003標準中已經給出的幾種LM磨耗型踏面外形的基礎上基于車輛動力學性能做相應的型面優化。

表1 2009年廣州地鐵5號線磨耗狀況對比

選取我國廣泛應用的LM磨耗型踏面進行分析。其幾何特性如圖1所示。根據TBT449—2003分析,LM踏面外形能夠依據輪緣厚的尺寸不同分為4種,分別為26,28,30,32 mm[16],綜合4種磨耗型踏面的幾何特性,有幾個不變的參數可以進行確定。我們可以分析得到需要采用5個參數能夠描述這段圓弧r1,r2,r3,xR3,xR4,從左向右4段圓弧依次為圓OR1,OR2,OR3,OR4,A點為輪緣厚測量點,對于確定了輪緣厚值的預鏇輪對,A點和E點值為固定值A(-39.8,-12),E(42.3,2.52)。

圖1 LM磨耗型踏面幾何特性(單位:mm)

由已知的A點坐標可得到兩個方程

(3)

可得圓OR1的圓心坐標O1(xR1,yR2)。

由于圓OR3必經過踏面基準點,則現已知圓OR3的圓心橫坐標xR3以及半徑r3,則可得出O3的縱坐標yR3。

已知E點坐標,根據圓OR4的特點及圓OR3和圓OR4外切的特性,可得以下兩個方程

(4)

可以算出yR4和r4。

(5)

可得出圓OR2的圓心坐標O2(xR2,yR2)。

至此,4段圓弧均可唯一確定。表1即為使用的設計變量及其范圍的依據。

表2 標準試驗系統結果數據設計變量

2 單目標車輪踏面Kriging優化方案

2.1 目標函數

磨耗性能是用于表征車輛的曲線通過性能,采用英國首創的磨耗指數來表示

WZ=Fxξx+Fyξy

(6)

式中,WZ為列車輪對的磨耗指數;Fx,Fy分別為軌道對輪對的左右軌的縱向和橫向蠕滑力;ξx,ξy分別為輪對的左右輪與鐵軌的縱向和橫向蠕滑率。

目標函數f(x)與輪軌間損耗有關,為

(7)

WZL(t)=FxL(t)ξxL(t)+FyL(t)ξyL(t)

WZR(t)=FxR(t)ξxR(t)+FyR(t)ξyR(t)

(8)

式中:WZL(t),WZR(t)分別為某時間點列車導向輪對雙輪的磨耗指數;FxL(t),FyL(t),FxR(t),FyR(t)分別為某時刻導向輪對的左右軌的縱向和橫向蠕滑力;ξxL(t),ξyL(t),ξxR(t),ξyR(t)分別為某時刻導向輪對的左右輪與軌道的縱向和橫向蠕滑率;S1,S2,S3分別為[t0,t1],[t1,t2],[t2,t3]時間段內車輛運行的距離。

2.2 約束條件

為避免發生事故,減輕輪軌損耗,約束條件采用脫軌系數,輪重減載率,平穩性指標。

(1)約束函數g1(x)與脫軌系數相關,則

(9)

式中,Lv0為車輪的脫軌系數的最大上限值,根據TB/T2360-93取0.8;PL(t),PR(t)分別為某時刻鐵軌作用于導向輪對的垂向力;fir(·,lfir,ufir)為帶通濾波器;lfir,ufir分別為帶通濾波器的截斷低頻率和截斷高頻率。

(2)約束函數g2(x)和輪重減載率相關,則

(10)

式中,ΔP為輪重減載量;P為平均輪軌垂向力;P0=0.6;P0為輪重減載率安全允許限值。

(3)目標函數g3(x)與平穩性相關,則

g3(x)=SPl-SPl0≤0

(11)

按照《客運專線鐵路工程竣工驗收動態檢測指導意見》規定,SPl0=2.5。

2.3 模型精度檢驗

依次對其他指標進行精度判定,并采用均方根誤差值進行誤差計算,該評價指標的最大值為0.2。由表3可以看出,各指標的誤差值均在合格范圍內,因而此模型可以進行仿真并對其進行計算。

表3 單目標優化均方根誤差

2.4 Kriging優化模型的求解過程

模型求解過程見圖2,具體步驟如下。

(1)給定待優化的初始踏面形狀,車體動力學模型,確定樣本點的初始空間E0,中心點x0,步長δ0。

(2)在初始空間進行最優拉丁超立方試驗設計,進而得到樣本點。

(3)利用Kriging插值法對目標函數及約束條件進行插值,建立車輪踏面優化模型。回歸部分采用一階多項式模型。隨機過程采用GAUSS函數。進行模型精度檢驗。

(4)求解優化問題

(12)

利用懲罰因子將優化問題轉化成只帶邊界約束條件的優化問題,將約束條件通過罰因子表示出來,即

(13)

Ωg={i|gi(x,θ)>0,j=1,2,…,N}

(14)

則可得子優化問題

(15)

再利用遺傳算法即可得到優化問題的最優解x(k+1)*以及所對應的代理模型的目標值。

(5)判斷最優解的信息符合收斂準則與否。當第一次迭代時,直接執行步驟6。當迭代兩次之后時,計算優化過程中第k次與第k-1次所求的可能最優解所對應的模型響應值的相對差值,判斷該值符合收斂標準與否。若收斂,過程結束,則第4個步驟獲得的最優解就是真實模型的最優解;若不收斂,即進入第6個步驟。

(6)采用信賴域更新方法。可以根據下式來確定信賴域的更新半徑δ(k+1)

(16)

式中,m1,m2,S1,S2為系數,其中01,0

第k+1次迭代的子空間Ek+1由式確定

n=1,2,…,N}

圖2 單目標優化Kriging模型求解流程

3 多目標車輪踏面Kriging優化方案

3.1 目標函數

與磨耗指數相關的目標函數f1(x)見2.1節。

與脫軌系數相關的目標函數

(17)

式中,QL(t),QR(t)分別為鋼軌受到車輪作用產生的橫向力;PL(t),PR(t)分別為左右車輪的垂向力;fir(·,lfir,ufir)為帶通濾波器,lfir,ufir分別為帶通濾波器的截斷低頻率和截斷高頻率。

與平穩性指標相關的目標函數

f3(x)=SPl

(18)

3.2 約束條件

與輪軌橫向力有關的約束函數

g1(x)=max{|fir(QL(t),lfir,ufir)|,

|QR(t),lfir,ufir|}-Q0≤0

(19)

與輪軸橫向力相關的約束函數

g2(x)=max{|fir(HL(t),lfir,ufir)|,

|HR(t),lfir,ufir|}-H0≤0

(20)

3.3 模型精度檢驗

采用均方根誤差值進行誤差計算,結果如表4所示。由表4可以看出,各指標的誤差值均在合格范圍內,因而可以對該模型進行仿真。

表4 多目標優化均方根誤差

3.4 優化模型的求解

以磨耗指數、脫軌系數和平穩性為目標函數的多目標Kriging優化問題

(21)

基于多目標的Kriging優化模型的求解,與單目標的區別之處在于目標函數的不同。運用遺傳算法對其進行求解時,對目標函數采用線性加權和法,由于3個目標函數并無明顯的主次關系,因此權重根據函數值本身不同的數量級確定。再引入懲罰因子,則式(18)轉變為

(22)

再利用遺傳算法對式(15)求解即可得到優化問題的最優解x(k+1)*以及所對應的代理模型的目標值,求出的最優解根據信賴域更新方法判斷符合收斂準則與否[17]。

4 仿真結果分析

4.1 單目標優化仿真結果

初始踏面采用TBT449—2003中所規定的輪緣厚為30 mm的磨耗型車輪踏面,初始樣本中心點為x0=[14,100,500,32,15],δ0=[1.5,10,20,4,2]。進行優化設計。

在第5次迭代時,目標函數的變化值小于給定的閥值,優化過程結束。求出的最優踏面與標準踏面的對比如圖3所示。圖4是兩種踏面的等效錐度的對比。

圖3 單目標優化標準踏面和最優踏面的比較

圖4 單目標優化標準踏面和最優踏面等效錐度比較

由圖4可以看出,優化后的等效錐度要顯著大于標準踏面。在偏移量6 mm時,優化的踏面的等效錐度數值為0.208 5,而同樣的LM標準踏面的等效錐度為0.171 9。當橫移量呈現輪緣貼靠的時候,等效錐度迅速變大。

分析曲線通過性能,根據相關規定設置了3種理想曲線通過軌道,半徑分別為400、600、800 m。以條形圖的形式比較不同條件下最優踏面通過不同半徑彎道曲線的動力學各指標的變化情況[18],如圖5所示。

圖5 單目標優化標準踏面和最優踏面曲線通過能力比較

由圖5(a)可以看出,采用最優踏面的列車在通過曲線軌道時,目標函數磨耗指數有較大幅度的降低,經過計算發現列車在通過3種不同曲線半徑軌道時的磨耗指數降幅分別為13.3%、15.2%和16.6%,曲線半徑越小,降幅越大。但對于其他的指標并不一定有著相同的趨勢。

脫軌系數在半徑400 m和600 m的軌道上分別惡化8.9%與6.7%,當曲線半徑增大到800 m時,最優踏面的脫軌指標發生優化,幅度較小,為3.6%。輪軌橫向力指標在半徑400 m和800 m的軌道上,最優踏面進行了較小幅度的優化,優化率僅為2.75%和1.5%。半徑為600 m時,輪軌橫向力甚至惡化1.6%。輪重減載率在3個不同曲線軌道上均發生了一定的惡化,分別為6.2%、5%和0.64%。輪軸橫向力同樣發生微弱的優化與惡化。

由以上分析可以得出,以磨耗指數單個目標函數做優化,對磨耗指數本身有著較好的優化效果,但對于其他指標來說,優化效果不明顯,甚至可能會變差。同時,會造成直線運行能力的惡化。

4.2 多目標優化仿真結果

初始踏面采用TBT449—2003中所規定的輪緣厚為30 mm的磨耗型車輪踏面進行多目標優化設計。Kriging優化模型的初始樣本中心點為x0=[14,100,500,32,15],δ0=[1.5,10,20,4,2]。利用該優化模型求出的最優踏面與初始踏面的對比如圖6所示。其等效錐度的對比如圖7所示。

圖6 多目標優化標準踏面和最優踏面的比較

圖7 多目標優化標準踏面和最優踏面等效錐度比較

由圖7可以看出,橫移量在0~4 mm,優化后的踏面的等效錐度與標準踏面基本持平,橫移量在4~7 mm,等效錐度略小于標準踏面,橫移量大于7 mm后,優化踏面的等效錐度略大于標準踏面。

同樣設置了3種半徑分別為400、600、800 m的理想曲線通過軌道。車輛以適當的速度通過帶有軌道激勵的曲線軌道。表5為最優踏面通過不同半徑的曲線軌道的性能指標值。同時,將優化后的指標值與標準踏面進行對比,計算得到該優化模型的優化率,如表6所示。

表5 裝載最優踏面的列車曲線通過能力指標

表6 裝載最優踏面的列車曲線通過能力的優化率 %

由表6可以看出,在半徑為400 m的曲線軌道上裝載優化后的踏面的列車的磨耗指數優化率為8.24%,較單目標模型在相同條件下的13.3%的優化率要較小一些。另外兩種大半徑的軌道上磨耗指數的優化率同樣比其要小。但是約束條件都有不同程度的優化,并未發生某些性能變差的情況。半徑為400 m的曲線軌道上,輪軌橫向力的優化率為5.52%,輪軸橫向力的優化率為6.07%,脫軌系數的優化率為6.62%,輪重減載率的優化率為5.95%。半徑為600 m的曲線軌道上,磨耗指數的優化率與400 m軌道相比略微增大為8.42%,輪軌橫向力的優化率為9.12%,輪軸橫向力的優化率為5.09%,脫軌系數的優化率為4.74%,輪重減載率的優化率為5.79%。半徑為800 m的較大半徑曲線軌道上,磨耗指數的優化率與400 m軌道相比提高幅度較大,為10.57%,輪軌橫向力的優化率為8.68%,輪軸橫向力的優化率為6.66%,脫軌系數的優化率為5.28%,輪重減載率的優化率為4.2%。總的來說,列車的曲線通過能力在3種曲線軌道上都得到明顯改善。這與等效錐度的走勢大致相同,在橫移量較大時,等效錐度較高有利于曲線通過。

由如上分析可得,采用線性加權和與信賴域相結合的方法對建立的Kriging優化模型進行優化,最優解只需要信賴域變更5次即可獲得,雖然磨耗指數的優化率沒有之前單目標的優化率高,但沒有某個指標出現惡化情況,導致曲線通過能力不平衡。由此說明所采用的多目標優化模型,在平衡不同性能指標方面有著較好的優化效果。相比建立的單目標優化模型,多目標優化具有較好的均衡性。

5 結論

本文圍繞城軌車輪鏇修后的踏面外形輪廓展開研究,提出了基于Kriging代理模型的單目標優化方案和多目標優化方案,以磨耗指數、平穩性和脫軌系數多個指標為目標函數建立Kriging多目標優化模型,采用線性加權和與信賴域更新相結合的優化算法求解,從優化結果可以看出,裝載優化踏面的列車的直線平穩性與曲線通過能力都可以得到一定程度的提高。多目標優化與單目標相比很大的優勢就是均衡性比較好,同時優化多個相互影響的目標,使系統的各方面綜合性能得到提高,從而在多個優化目標之間尋找一個折中點,使系統的綜合性能得到提高。

[1] 王永亮,張星臣,徐彬,等.城市軌道交通網絡化列車開行方案優化方法[J].中國鐵道科學,2012(5):120-126.

[2] 劉燦龍.地鐵輪軌磨耗的初步研究[J].城市軌道交通研究,2008(10):57-59.

[3] 陸縉華.廣州地鐵1號線車輛的磨損情況分析[J].電力機車技術,2001(3):16-17.

[4] 鄧軍.廣州地鐵車輛輪對損傷行駛及原因淺析[J].機車電傳動,2008(1):71-74.

[5] 李海東.基于CRH2的車輪磨耗和經濟型鏇修踏面的研究[D].南昌:華東交通大學,2014.

[6] 李煜.車輪磨耗對車體動力學性能影響研究[D].大連:大連交通大學,2014.

[7] 常崇義,王成國,金鷹.基于三維動態有限元模型的輪軌磨耗數值分析[J].中國鐵道科學,2008(4):89-95.

[8] Braghin F, Lewis R, et al. A mathematical model to predict railway wheel profile due to wear[J]. Wear, 2006,261(2):1253-1264.

[9] 員華,肖勝強,汪洋.基于磨耗量統計的輪對等級鏇修可行性研究[J].城市軌道交通,2006(1):43-45.

[10] Pascual F,Marcos J A.A common playground for design and maintenance engineering in the Taglo engineering cycle[C]. Proceeding of the Joint Rail Conference, 2004:193-199.

[11] 許宏,員華,王凌,等.基于高斯過程的地鐵車輛輪對磨耗建模及其鏇修策略優化[J].機械工程學報,2010(7):88-95.

[12] 王珍.動車組整車輪對鏇修優化決策模型的研究[D].成都:西南交通大學,2013.

[13] Jendel T. Prediction of wheel profile wear-comparisons with field measurements[J]. Wear, 2002,253(4):89-95.

[14] 謝鋒.基于近似模型的車架不確定性多目標優化[D].長沙:湖南大學,2013.

[15] Lophaven S N,Nielsen H B,Sondergaard J.DACE-A Matlab Kriging Toolbox,Version2.0[Z]. Technical Report,2002.

[16] 朱士友,潘麗莎,員華.輪對等級鏇修對車輛平穩性的影響分析[J].城市軌道交通研究,2016(7):46-48.

[17] 陳國棟.基于代理模型的多目標優化方法及其在車身設計中的應用[D].長沙:湖南大學,2012.

[18] 國家標準局.GB5599—85 鐵道車輛動力學性能評定和試驗鑒定規范[S].北京:中國標準出版社,1986.

猜你喜歡
優化模型
一半模型
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 亚洲第一天堂无码专区| 成人在线综合| 国产麻豆永久视频| 国产精品福利尤物youwu | 婷五月综合| 久久这里只有精品国产99| 国产乱人乱偷精品视频a人人澡| 日韩在线视频网站| 伊人色在线视频| 五月天久久综合国产一区二区| 直接黄91麻豆网站| 欧美日韩理论| 欲色天天综合网| 伊人色天堂| 精品国产免费观看一区| 成人亚洲视频| 精品国产免费第一区二区三区日韩| 2020久久国产综合精品swag| 欧美特黄一级大黄录像| 五月天在线网站| 少妇被粗大的猛烈进出免费视频| 欧美视频在线播放观看免费福利资源 | 99re经典视频在线| 美女视频黄又黄又免费高清| 国产91高清视频| 最新国语自产精品视频在| 午夜不卡福利| 亚洲欧美日本国产专区一区| 欧美一区精品| AⅤ色综合久久天堂AV色综合| 国产浮力第一页永久地址| 精品免费在线视频| 高清国产在线| 自慰网址在线观看| 三级欧美在线| 91久久国产综合精品| 国产成人高清精品免费| аv天堂最新中文在线| 欧美在线精品怡红院| 在线中文字幕日韩| 欧美一级在线| 亚洲一级色| 欧美在线中文字幕| 国产精品第页| 欧美午夜在线观看| 日韩精品一区二区三区大桥未久 | 亚洲最大综合网| 亚洲高清免费在线观看| 久久精品最新免费国产成人| 精品无码一区二区三区电影| 乱人伦中文视频在线观看免费| 亚洲无码高清视频在线观看| 成人在线观看不卡| 激情综合网激情综合| 一级香蕉人体视频| 久久亚洲综合伊人| 视频一区视频二区中文精品| 精品久久人人爽人人玩人人妻| 99在线小视频| 91麻豆久久久| 欧美性猛交一区二区三区 | 亚洲人成人无码www| 欧美成人区| 欧美日韩国产高清一区二区三区| 71pao成人国产永久免费视频| 国产亚洲欧美另类一区二区| 欧美成人a∨视频免费观看| 久久午夜夜伦鲁鲁片不卡| 亚洲AV一二三区无码AV蜜桃| 91探花国产综合在线精品| 黄色网站不卡无码| 日本亚洲欧美在线| 人妻精品久久无码区| 日本久久免费| 嫩草影院在线观看精品视频| 国产亚洲第一页| 乱人伦中文视频在线观看免费| 久久成人免费| 国产成人久久综合777777麻豆| 日韩在线播放中文字幕| 国产无码精品在线| 亚洲无码精品在线播放|