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

考慮打磨量的重載鋼軌打磨廓形優(yōu)化設(shè)計

2022-08-25 03:03:52康彥兵張華鵬
西南交通大學(xué)學(xué)報 2022年4期
關(guān)鍵詞:優(yōu)化

吳 磊 ,康彥兵 ,董 勇 ,張華鵬

(西南交通大學(xué)機械工程學(xué)院,四川 成都 610031)

重載線路因具有貨運頻繁、線路復(fù)雜、牽引重量大等特點,容易發(fā)生嚴重的輪軌磨耗損傷,尤其在重載曲線鋼軌段,其損傷形式主要表現(xiàn)為鋼軌側(cè)磨和鋼軌壓潰[1],造成輪軌表面應(yīng)力分布不良,影響輪軌動力學(xué)性能和列車運行安全.輪軌接觸幾何型面的不平順會造成輪軌表面接觸應(yīng)力的分布差異[2]及車輛結(jié)構(gòu)疲勞等問題[3],進而引起鋼軌損傷.采用鋼軌打磨目標設(shè)計廓形的方式,可有效改善輪軌幾何型面匹配性能,預(yù)防鋼軌接觸疲勞,減輕鋼軌損傷[4].

在鋼軌廓形目標型面優(yōu)化設(shè)計方面:Persson等[5]以列車脫軌系數(shù)、輪軌滾動接觸疲勞、最大輪對橫向位移等動力學(xué)指標作為懲罰因子加權(quán)建立罰函數(shù),采用遺傳算法得到優(yōu)化性能較好的目標鋼軌廓形;Choi等[6]以最小化曲線段鋼軌磨損為優(yōu)化目標,提出一種基于遺傳算法的非對稱式鋼軌廓形設(shè)計方法,采用動力學(xué)理論分析評價輪軌力大小以及脫軌情況,所優(yōu)化后得到的內(nèi)外鋼軌廓形相比初始廓形具有更小的磨耗指數(shù);Wang等[7]引入輪軌接觸幾何關(guān)系建立廓形優(yōu)化函數(shù),運用Sequential Quadratic Programming算法求解得到匹配LM (標準磨耗型)車輪性能更好地打磨鋼軌廓形;毛鑫等[8-9]基于輪徑差函數(shù)曲線,建立直線和曲線鋼軌打磨目標廓形的優(yōu)化模型,得到輪軌接觸分布良好的目標廓形;王亮等[10]利用 NURBS (non-uniform rational b-splines)曲線構(gòu)建以輪軌接觸性能和輪對曲線通過能力Kriging代理模型為目標函數(shù)的非對稱打磨目標型面多目標優(yōu)化模型,采用 NSGA-II (non dominated sorting genetic algorithm-II)算法優(yōu)化使輪對曲線通過能力改善.

目前對于鋼軌打磨廓形優(yōu)化設(shè)計,主要通過改善輪軌接觸分布和減少鋼軌磨損等方式,并未在設(shè)計過程中直接考慮目標廓形打磨量.肖杰靈等[11-12]提出鋼軌打磨廓形設(shè)計應(yīng)考慮打磨量以及打磨量最小原則,但并未涉及打磨量最小化的具體設(shè)計方法.本文建立目標廓形打磨量計算方法,并將打磨量與主導(dǎo)重載線路鋼軌損傷的磨耗損傷、壓潰損傷進行整合,發(fā)展重載鋼軌打磨廓形多目標優(yōu)化設(shè)計方法,設(shè)計鋼軌打磨廓形.

1 鋼軌打磨量計算方法

綜合考慮輪軌磨耗指數(shù)、輪軌接觸應(yīng)力和打磨量因素,進行鋼軌打磨目標廓形優(yōu)化設(shè)計.磨耗指數(shù)和輪軌接觸應(yīng)力可通過Simpack動力學(xué)軟件仿真計算.此外,通過建立任意鋼軌廓形與目標鋼軌廓形的對齊算法計算目標廓形打磨面積,考慮到用鋼軌橫截面形狀即可描述鋼軌廓形幾何特征,因此,用鋼軌橫截面面積變化量表征打磨量,并將打磨量整合到廓形優(yōu)化模型中.其中,目標廓形源于通過遺傳優(yōu)化算法選擇、交叉、變異算子自動創(chuàng)建的鋼軌廓形,優(yōu)化迭代結(jié)束后最后一代廓形中綜合性能最好的個體即為打磨廓形.

1.1 實測廓形分析

測得某線路不同里程段鋼軌廓形,選取具有代表幾何特征的兩種實測鋼軌型面進行分析.如圖1(a)所示,曲線段實測鋼軌廓形1軌距角處的鋼軌磨損相比較實測鋼軌廓形2更為嚴重,A1、A2兩點分別為實測廓形1和2軌頂處的兩點.此外,圖1(b)打磨目標廓形與實測廓形坐標正方向相反,由于打磨目標廓形需要輸入到車輛-軌道動力學(xué)模型中計算,而根據(jù)Simpack建立的車輛-軌道動力學(xué)模型系統(tǒng)以打磨目標廓形坐標正方向為標準.為了計算實測廓形打磨成目標廓形之后的材料去除截面面積,需以目標廓形坐標朝向為基準對實測鋼軌廓形坐標進行旋轉(zhuǎn)變換預(yù)處理,使打磨目標鋼軌廓形與實測廓形進行重疊放置.

圖1(a)中,將A1、A2作為旋轉(zhuǎn)變換中心,通過編寫數(shù)值程序自動找到該旋轉(zhuǎn)中心點,根據(jù)式(1)對兩種實測廓形分別繞A1、A2進行旋轉(zhuǎn)變換預(yù)處理,使其靠近預(yù)期目標廓形,如圖2所示.圖2中:α1、α2、β為通過反三角函數(shù) arctank1、arctank2、arctank0求得的夾角.

圖1 鋼軌廓形幾何型面Fig.1 Geometric profile of rail

式中:Rph為實測廓形h旋轉(zhuǎn)變換處理后的矩陣;Msph為初始實測廓形h的負向平移矩陣,如式(2);Mrot為繞原點旋轉(zhuǎn)的旋轉(zhuǎn)矩陣,如式(3);Msmh為初始實測廓形h的正向平移矩陣,如式(4);Rinih為實測廓形h初始矩陣,如式(5).

式中:(xh,yh)為旋轉(zhuǎn)中心點橫、縱坐標;θ為繞原點的旋轉(zhuǎn)角度,此處值為π;Rxh、Ryh分別為初始實測廓形h數(shù)據(jù)點的橫、縱坐標向量;E為單位列向量.

1.2 鋼軌廓形橫向?qū)R

圖2中兩實測鋼軌廓形通過旋轉(zhuǎn)變換已被置于預(yù)期目標廓形附近,為了使實測廓形與目標廓形軌距角側(cè)斜直線段重合,需要再次旋轉(zhuǎn)實測廓形,使其與目標廓形內(nèi)側(cè)斜直線在誤差允許范圍內(nèi)斜率相等.旋轉(zhuǎn)角度計算方法如下所述.

1) 保持目標廓形不動,繞坐標原點旋轉(zhuǎn)上述預(yù)處理后的實測廓形,使其內(nèi)側(cè)與預(yù)期目標鋼軌廓形平行.根據(jù)圖2可知斜率及旋轉(zhuǎn)角度為

圖2 實測鋼軌廓形旋轉(zhuǎn)預(yù)處理Fig.2 Pre-processing of measured rail profile rotation

式中:k1、k2、k0均為實測廓形1、2以及目標對齊廓形內(nèi)側(cè)斜直線斜率,實際中廓形軌距角側(cè)(即鋼軌內(nèi)側(cè))廓形段經(jīng)常為近似直線的曲線,可通過最小二乘法對鋼軌廓形內(nèi)側(cè)曲線段進行線性擬合得到斜直線及其斜率;δi為旋轉(zhuǎn)角度;jn為選取廓形內(nèi)側(cè)近似斜直線段的離散點編號;n為該斜直線段離散點數(shù)目;分別為該段點的橫、縱坐標;x為該段數(shù)據(jù)橫坐標平均值.

根據(jù)式(1)計算并獲得與目標廓形內(nèi)側(cè)斜直線段平行的實測廓形矩陣,但在該公式中,初始矩陣為Rph,旋轉(zhuǎn)中心為原點(0,0),旋轉(zhuǎn)矩陣Mrot中θ值為δ,其結(jié)果如圖3所示,陰影部分面積Sz為兩廓形內(nèi)側(cè)段某區(qū)間[13]直線圍成的面積,z為大于等于1的正整數(shù).

圖3 實測廓形與目標廓形內(nèi)側(cè)直線段平行示意Fig.3 Parallelling of the straight line inside between the measured profile and the target profile

2) 對上述處理后的實測鋼軌廓形左、右平移,構(gòu)建平移向量,使其與目標廓形內(nèi)側(cè)斜直線段重合,其迭代流程見圖4.采用最小面積法判斷內(nèi)側(cè)斜直線是否重合[14],如式(8).

圖4 鋼軌廓形橫向?qū)R流程Fig.4 Flowchart for horizontal alignment of rail profile

式中:ε為精確度,即面積閾值(rshreld),取為0.001 5 mm.

選取縱坐標區(qū)間長度為5 mm,設(shè)置移動步長為0.001 mm,則可得到實測廓形的單次平移矩陣為

式中:rstep為移動步長,表示沿x軸方向移動rstep距離,0表示沿y軸方向不進行平移.

通過不斷迭代移動步長改變陰影部分面積,使該面積滿足式(8)時,即可判定鋼軌橫向?qū)R完成,對齊效果見圖5.

由圖5可知:在軌距角處,相比較實測廓形2,目標廓形并未被實測廓形1所包裹,實測廓形1并不能打磨成目標廓形的樣子.

圖5 實測廓形與目標廓形橫向?qū)R示意Fig.5 Horizontal alignment between the measured profile and the target profile

1.3 鋼軌廓形垂向?qū)R

在做垂向?qū)R算法設(shè)計時,需根據(jù)不同實測廓形軌距角處磨損差異,基于軌距角和軌頂建立不同的垂向?qū)R方式.基于軌距角處對實測鋼軌廓形1垂向平移,使其在軌距角處對齊,結(jié)果如圖6(a)所示.基于軌頂最高處對實測鋼軌廓形2垂向平移,使其在軌頂處對齊,結(jié)果如圖6(b)所示.

1.4 實測鋼軌廓形打磨量計算方法

圖6(a)和6(b)中的陰影部分均為由鋼軌廓形曲線圍成的不規(guī)則形狀,計算出該面積即可知所要求解的鋼軌總打磨量.由于實際鋼軌廓形數(shù)據(jù)是二維離散坐標點,且相鄰離散點間距離足夠小,可將圖中陰影部分看作封閉的多邊形,即將兩廓形曲線的橫縱坐標分別拼接為組成該多邊形的兩個向量,見式(10),采用求任意多邊形面積的計算方法即可求出該鋼軌廓形打磨面積的近似值,該近似值滿足本文優(yōu)化所需打磨量參數(shù)條件.

式中:X、Y為實測廓形和目標廓形陰影部分橫、縱坐標組成的向量,其中,前c個點為實測鋼軌廓形1或2的坐標點,后nt-1-c個點為目標鋼軌廓形坐標點,最后再添加點(x1,y1)是為了讓陰影面積對應(yīng)兩曲線構(gòu)成封閉多邊形;Sg為鋼軌打磨面積;[xmym]、[xm+1ym+1]為X、Y第m和m+ 1個點組成的向量;nt為圍成打磨面積區(qū)域(即圖6中陰影部分)的廓形點數(shù)目.

圖6 實測廓形與目標廓形垂向?qū)R示意Fig.6 Vertical alignment between the measured profile and the target profile

2 打磨廓形優(yōu)化方法

輪軌磨耗和輪軌接觸應(yīng)力是描述輪軌相互作用的兩個重要指標,將鋼軌打磨為特定的幾何形狀可以顯著改善磨耗和應(yīng)力.此外,鋼軌打磨量過大會縮減鋼軌服役壽命并消耗過多的打磨資源.因此,本文采用輪軌磨耗指數(shù)、輪軌接觸應(yīng)力和鋼軌打磨量為設(shè)計鋼軌打磨廓形的子目標函數(shù),通過加權(quán)法整合,構(gòu)建綜合目標函數(shù).用Simpack建立C80貨車車輛-軌道動力學(xué)模型,聯(lián)合MATLAB數(shù)值計算軟件編寫遺傳數(shù)值優(yōu)化算法求解目標函數(shù)最優(yōu)解,得到鋼軌打磨目標廓形.需要說明的是第1節(jié)中的目標廓形是為了介紹廓形對齊及打磨量計算方法而引出,具有任意代表性,而此處的打磨目標廓形是通過廓形優(yōu)化設(shè)計方法輸入特定的工況參數(shù)得到的目標廓形.

2.1 綜合優(yōu)化目標函數(shù)

根據(jù)文獻[15]設(shè)計鋼軌廓形變量C,分別賦予輪軌磨耗指數(shù)、輪軌接觸應(yīng)力和鋼軌打磨量不同的影響權(quán)重系數(shù)并加權(quán)建立綜合優(yōu)化目標函數(shù),如式(11),其中,輪軌磨耗指數(shù)和接觸應(yīng)力均選擇車輛穩(wěn)定運行部分的平均磨耗指數(shù)和輪軌接觸應(yīng)力作為評判優(yōu)化廓形性能的磨耗和應(yīng)力指標.

式中:F(C)為綜合優(yōu)化目標函數(shù)值;w1id為輪軌磨耗指數(shù)權(quán)重,w2id為輪軌接觸應(yīng)力權(quán)重,w3id為打磨量權(quán)重;f1id(C)為輪軌磨耗指數(shù)歸一化指標,f2id(C)為輪軌接觸應(yīng)力歸一化指標[16],f3id(C)為上述計算的打磨量Sg歸一化指標.

2.2 目標函數(shù)最優(yōu)解

獲得最優(yōu)鋼軌打磨廓形個體就是求取綜合優(yōu)化目標函數(shù)的最優(yōu)解,即目標函數(shù)的最小值.遺傳算法因具有良好的全局適應(yīng)性和全局搜索性經(jīng)常被用于求解優(yōu)化問題.

選取實際重載線路段不同的實測鋼軌廓形作為初始種群,這些廓形分別來自直線位置和曲線位置,它們新鋪設(shè)時均為標準CN60軌廓形、軌底坡1∶40.運用改進的 NURBS (non-uniform rational b-splines)曲線構(gòu)造方法[17]對鋼軌廓形進行參數(shù)化表示,將其導(dǎo)入建立的Simpack車輛動力學(xué)模型中,計算得到輪軌的磨耗指數(shù)、輪軌接觸應(yīng)力.運用第1節(jié)提到的鋼軌打磨量計算方法計算出實測廓形打磨成任意目標種群廓形的打磨量,并進行數(shù)值歸一化計算,得出每代種群不同廓形個體的目標函數(shù)值F(C),挑選每代最優(yōu)個體(目標函數(shù)最小值對應(yīng)個體)進行后續(xù)計算結(jié)果分析,取其倒數(shù)為算法的適應(yīng)度值,根據(jù)此值篩選出符合設(shè)計要求的廓形個體作為父代個體進行選擇、交叉、變異等操作產(chǎn)生下一子代廓形個體,經(jīng)過多次迭代更新,達到算法最大迭代次數(shù)后,目標廓形更新結(jié)束,選取最后的子代目標函數(shù)值最小的個體作為最優(yōu)設(shè)計的鋼軌打磨廓形.

2.3 輪軌磨耗指數(shù)

輪軌磨耗指數(shù)為

式中:Tx、Ty、T?分別為接觸點的縱向蠕滑力、橫向蠕滑力和自旋蠕滑力矩; ξx、 ξy、 ξ?分別為輪軌接觸斑處的縱向蠕滑率、橫向蠕滑率和自旋蠕滑率.

3 優(yōu)化方法的應(yīng)用

選擇第1節(jié)中提到的實測鋼軌廓形1,應(yīng)用第2節(jié)鋼軌打磨廓形優(yōu)化方法,選取不同的權(quán)重參數(shù)建立不同的打磨廓形設(shè)計優(yōu)化策略分析輪軌磨耗指數(shù)、輪軌接觸應(yīng)力、鋼軌打磨量的迭代演化情況.

3.1 優(yōu)化策略

以重載C80車動力學(xué)參數(shù)模型為載荷工況,并保持工況載荷條件不變,對不同優(yōu)化子目標賦予不同的影響權(quán)重系數(shù),分析不同權(quán)重系數(shù)的優(yōu)化策略對優(yōu)化結(jié)果的影響,具體參數(shù)見表1.

表1 優(yōu)化策略Tab.1 Optimization strategies

工況1同時考慮了輪軌磨耗、接觸應(yīng)力以及鋼軌打磨量的子目標,主要對輪軌磨耗指數(shù)、輪軌接觸應(yīng)力和鋼軌打磨量均進行重點優(yōu)化;工況2同時考慮對磨耗指數(shù)和接觸應(yīng)力子目標進行優(yōu)化;工況3對磨耗指數(shù)和接觸應(yīng)力權(quán)重均賦予0,只考慮打磨量因素,主要針對打磨量子目標進行優(yōu)化.

3.2 鋼軌廓形迭代演化分析

動力學(xué)分析結(jié)果表明,鋼軌廓形迭代計算的變化對車輛4個輪對均具有相同的影響趨勢,尤其一位輪對,即車輛運行的導(dǎo)向輪對,變化速率更為明顯.每代遺傳保留下來的鋼軌廓形幾何型面有多種,選取一位導(dǎo)向輪對計算該輪對分別和該多種鋼軌廓形輪軌磨耗指數(shù)及接觸應(yīng)力的平均值,同時計算鋼軌打磨量的最大值,分析每代鋼軌廓形各優(yōu)化目標的迭代演化情況.

1) 輪軌磨耗指數(shù)迭代演化

圖7為一位輪對曲線外軌和內(nèi)軌的輪軌磨耗指數(shù)迭代演化過程,由圖7可知:工況1的輪軌磨耗指數(shù)變化速率最大,收斂時相比初始廓形外軌和內(nèi)軌平均降低了68.9%;3種工況對于外軌的磨耗指數(shù)優(yōu)化效果更明顯,均呈下降趨勢,而工況2、3對于內(nèi)軌的磨耗指數(shù)優(yōu)化前后無明顯變化;同時考慮磨耗、應(yīng)力和打磨量的優(yōu)化策略對輪軌磨耗優(yōu)化效果最好.

圖7 輪對磨耗指數(shù)迭代演化過程Fig.7 Iterative evolution process of wheelset wear index

2) 輪軌接觸應(yīng)力的迭代演化

圖8為一位輪對接觸應(yīng)力指數(shù)迭代演化過程,分析圖8可知:工況1、2對于內(nèi)軌的接觸應(yīng)力優(yōu)化效果更為明顯;工況1下降速率最快,收斂狀態(tài)下比初始廓形內(nèi)軌接觸應(yīng)力下降了39.1%,而工況2下降了37.7%;工況3對于輪軌接觸應(yīng)力指標起到了負優(yōu)化的作用,其外軌和內(nèi)軌接觸應(yīng)力值相比初始廓形增長幅度明顯,表明該工況下得到的優(yōu)化鋼軌廓形與該車輛模型中的車輪匹配,其輪軌接觸應(yīng)力過大,而較大的輪軌接觸應(yīng)力會加速鋼軌接觸疲勞,應(yīng)盡力避免.

此外,圖8工況3接觸應(yīng)力最大的原因是輪對通過曲線時,曲線外軌更容易在軌肩處接觸,軌肩處鋼軌廓形曲率半徑小,使得接觸應(yīng)力大;圖7中工況3磨耗指數(shù)小的原因是曲線外軌軌肩接觸,增加了輪徑差,輪對通過曲線時曲線外軌的輪軌蠕滑率下降,進而導(dǎo)致磨耗指數(shù)下降.綜上所述,同時考慮磨耗和應(yīng)力的工況2對輪軌接觸應(yīng)力優(yōu)化效果最好,而只考慮打磨量的工況3優(yōu)化效果最差.

圖8 一位輪對接觸應(yīng)力指數(shù)迭代演化過程Fig.8 Iterative evolution process of the first wheelset contact stress index

3) 鋼軌打磨量的迭代演化

圖9為鋼軌打磨量指數(shù)迭代演化過程,由圖9可知:工況1、2、3的每代鋼軌廓形最大打磨量下降均很明顯,并都在37代附近收斂;第200代廓形相比初始廓形分別降低了21.8%、11.3%、24.4%,工況3變化速率最大.

圖9 鋼軌打磨量指數(shù)迭代演化過程Fig.9 Iterative evolution process of grinding amount index

綜上可知:工況1由于給三者均分配了權(quán)重系數(shù),因而對于曲線外、內(nèi)軌的輪軌磨耗指數(shù),內(nèi)軌的接觸應(yīng)力以及鋼軌的打磨量均具有較好的優(yōu)化效果;工況2對于外軌磨耗指數(shù),內(nèi)軌接觸應(yīng)力優(yōu)化效果明顯,但由于沒考慮打磨量因素,打磨量下降速率相比較其他兩種工況較慢,優(yōu)化效果不夠明顯;工況3由于打磨量權(quán)值為1,占權(quán)值全部比重,打磨量下降速率最大,根據(jù)打磨量最小化原則,打磨量指標優(yōu)化效果最好,但曲線鋼軌接觸應(yīng)力卻顯著增大,可能加大輪軌損傷,應(yīng)當避免.因此,綜合3種指標迭代演化過程可知,同時考慮磨耗、應(yīng)力和打磨量的工況1對鋼軌廓形優(yōu)化效果最好.

3.3 最優(yōu)廓形性能分析

考慮到重載曲線鋼軌損傷以曲線外軌磨耗、內(nèi)軌壓潰為主的特點,用外軌磨耗指數(shù)和內(nèi)軌接觸應(yīng)力來評價設(shè)計廓形好壞,作為確定鋼軌打磨廓形設(shè)計方案的依據(jù).選取不同工況最優(yōu)廓形,得輪軌磨耗指數(shù)和接觸應(yīng)力如表2所示.

表2 外軌磨耗和內(nèi)軌接觸應(yīng)力Tab.2 High rail wear and low rail contact stress

分析表2可知:與實測廓形1相比,工況1最優(yōu)廓形外軌磨耗及內(nèi)軌接觸應(yīng)力下降最顯著;工況3對于輪軌磨耗優(yōu)化效果明顯,但卻導(dǎo)致接觸應(yīng)力急劇上升,因此對于應(yīng)力呈負優(yōu)化效果.綜上所述,同時考慮磨耗、應(yīng)力和打磨量的工況1效果最好.

4 結(jié) 論

1) 建立了基于目標廓形的鋼軌打磨廓形自動對齊算法和打磨量計算方法,該方法可以計算任意測試廓形與目標廓形的打磨量.將打磨量整合到鋼軌打磨廓形優(yōu)化設(shè)計綜合評價函數(shù)中,發(fā)現(xiàn)優(yōu)化迭代可以顯著降低打磨量.

2) 所建立的鋼軌打磨廓形優(yōu)化設(shè)計模型可以差異化地最小化輪軌磨耗、輪軌接觸應(yīng)力以及鋼軌目標廓形打磨量.在優(yōu)化過程中只考慮打磨量因素時,打磨量和輪軌磨耗會降低,但輪軌接觸應(yīng)力會顯著增加.

3) 同時考慮輪軌磨耗、輪軌接觸應(yīng)力、打磨量的影響時,或只考慮磨耗和應(yīng)力的影響時,磨耗、應(yīng)力顯著下降,但前者所得目標廓形所需的打磨量更少.因此,在鋼軌打磨目標廓形優(yōu)化設(shè)計中,有必要同時考慮磨耗、應(yīng)力和打磨量的影響.

猜你喜歡
優(yōu)化
超限高層建筑結(jié)構(gòu)設(shè)計與優(yōu)化思考
PEMFC流道的多目標優(yōu)化
能源工程(2022年1期)2022-03-29 01:06:28
民用建筑防煙排煙設(shè)計優(yōu)化探討
關(guān)于優(yōu)化消防安全告知承諾的一些思考
一道優(yōu)化題的幾何解法
由“形”啟“數(shù)”優(yōu)化運算——以2021年解析幾何高考題為例
圍繞“地、業(yè)、人”優(yōu)化產(chǎn)業(yè)扶貧
事業(yè)單位中固定資產(chǎn)會計處理的優(yōu)化
4K HDR性能大幅度優(yōu)化 JVC DLA-X8 18 BC
幾種常見的負載均衡算法的優(yōu)化
電子制作(2017年20期)2017-04-26 06:57:45
主站蜘蛛池模板: 国产精品成人一区二区| 狠狠亚洲五月天| 伊人久久大香线蕉成人综合网| 久久鸭综合久久国产| 青青草原偷拍视频| 国产精品护士| 国产在线八区| 久久一本精品久久久ー99| 久久精品无码一区二区日韩免费| 国产美女主播一级成人毛片| 在线精品自拍| 婷婷久久综合九色综合88| 国产成人精品18| 欧美精品黑人粗大| 国产成人精品18| 国产国拍精品视频免费看| 一级一级特黄女人精品毛片| 在线观看亚洲天堂| 免费国产小视频在线观看| 久久精品无码国产一区二区三区 | 欧美日韩综合网| 日韩黄色精品| 国产一在线观看| 日本高清视频在线www色| lhav亚洲精品| 2021国产精品自产拍在线| 中文字幕亚洲无线码一区女同| 青草91视频免费观看| 99青青青精品视频在线| a毛片基地免费大全| 国产真实乱了在线播放| 人妻精品全国免费视频| 成年人福利视频| 欧美另类精品一区二区三区| 视频一本大道香蕉久在线播放 | 国产亚卅精品无码| 亚洲一区毛片| 天天色综网| 免费久久一级欧美特大黄| 试看120秒男女啪啪免费| 久久久成年黄色视频| 波多野结衣在线一区二区| 精品欧美日韩国产日漫一区不卡| 中文字幕永久视频| 国产地址二永久伊甸园| www.精品国产| 欧美成人影院亚洲综合图| 欧美日本在线一区二区三区| 一本大道AV人久久综合| 一级黄色片网| 久久婷婷五月综合色一区二区| 成人91在线| 露脸国产精品自产在线播| 亚洲欧美在线综合一区二区三区| 亚洲激情区| 国产一级毛片yw| 成人福利在线免费观看| 国产在线自揄拍揄视频网站| 91小视频在线| av在线无码浏览| 精品无码专区亚洲| 这里只有精品在线| 久久精品中文无码资源站| 国产精品浪潮Av| 国产麻豆精品手机在线观看| 欧美成人综合视频| 欧美一区二区人人喊爽| 日韩专区第一页| 免费无码又爽又黄又刺激网站 | 成人国产精品网站在线看| 国产网友愉拍精品| 欧美在线中文字幕| 久久久久亚洲精品成人网| 就去色综合| 亚洲va精品中文字幕| 精品亚洲国产成人AV| 欧美不卡视频在线观看| 亚洲欧美成人网| 凹凸精品免费精品视频| 香蕉蕉亚亚洲aav综合| 国产精品视频第一专区| 欧美精品高清|