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

地鐵列車雙錐內嵌隔板矩形管的耐撞性優化

2021-03-17 01:24:14尚昱煌姜士鴻姚曙光胡正維
振動與沖擊 2021年5期
關鍵詞:模型設計

尚昱煌,許 平,姜士鴻,姚曙光,胡正維

(1. 中南大學 交通運輸工程學院,長沙 410075; 2. 中南大學 軌道交通安全教育部重點實驗室,長沙 410075;3. 中南大學 軌道交通安全國際合作聯合實驗室,長沙 410075; 4. 中國中車長春軌道客車股份有限公司,長春 130062)

隨著軌道交通行業的快速發展,乘員安全問題所受關注度日益增加。吸能裝置能在列車發生碰撞事故時,通過自身的塑性變形來耗散車輛的動能,從而有效減輕乘員所受傷害,因此國內外學者對其進行了大量研究[1-4]。其中金屬薄壁方管具有易于制造、吸能效率高等優點,因此被廣泛應用于列車的吸能裝置中。

Wierzbicki等[5]對薄壁矩形管在軸向壓縮下的變形模式進行了分析,并建立了折疊單元壓縮模型,通過超折疊單元理論對其平臺力MCF(mean crushing force)進行了預測。Zhou等[6-8]將折紋引入薄壁方管,并對折紋方管的耐撞性進行了深入研究,通過多種工況的試驗分析了軸向長度與厚度比值、缺陷形狀等因素對折紋方管的變形模式、MCF以及比吸能的影響。Qi等[9]將錐度引入薄壁方管,對錐形方管和直方管進行多個角度的斜向加載,發現錐形方管的比吸能SEA(specific energy absorption)均高于直方管,且峰值力PCF(peak crushing force)均低于直方管。Reid等[10]對錐度矩形管的MCF進行了理論推導,并通過準靜態壓縮試驗驗證了MCF公式預測的準確性。Song等[11]在薄壁方管上開設窗孔誘導結構,并將其與傳統薄壁方管分別在正向和斜向加載下進行了吸能性能的對比研究,發現在正向加載時窗孔誘導結構能有效降低薄壁方管的PCF且增加MCF,在斜向加載下窗孔誘導結構能防止薄壁方管發生彎曲壓潰。Zhang等[12]將縱向凹槽引入薄壁方管,發現適當尺寸的縱向凹槽能夠有效提高薄壁方管的SEA,并在一定程度上減小PCF。Asanjarani等[13]將橫向凹槽引入錐形方管,并以錐度、截面邊長、壁厚、凹槽尺寸和數量為變量進行了多目標耐撞性優化。李健等[14]將多個隔板嵌于直方管內,發現隔板的引入能改變直方管的變形模式,通過形成更多的褶皺而增加吸能量。Xu等[15-16]通過進行有限元分析和準靜態軸向壓縮試驗,對內嵌隔板方管進行了研究,并以吸能結構不同部分的厚度為設計變量,進行了多目標優化。Xing等[17]將隔板引入嵌套方管,并從幾何尺寸和厚度等多種設計角度對吸能結構進行了耐撞性優化。

以往研究通過引入錐度、凹槽、隔板、窗孔等途徑來提高金屬薄壁方管的耐撞性,尚未將錐度與隔板同時引入矩形截面的方管。本文以實際的工程項目為背景,對某地鐵車輛所使用的雙錐內嵌隔板矩形管進行研究,該吸能結構中矩形截面管的一組對稱面具有錐度,內部嵌入了多個隔板。通過建立有限元模型,對雙錐內嵌隔板矩形管和傳統矩形管的壓縮過程進行數值仿真,從而對比兩者耐撞性,并進行靜態壓縮試驗以驗證有限元分析的準確性。為進一步提高雙錐內嵌隔板矩形管的耐撞性,實現SEA最大化和PCF最小化,分別將有錐度和無錐度對稱面的厚度以及隔板厚度作為三個設計變量,通過拉丁超立方采樣進行試驗設計,利用徑向基函數方法構建出代理模型,采用多目標遺傳算法進行優化而得到更佳的設計方案。

1 模型建立

1.1 幾何模型

如圖1所示,一個包含了四個吸能單元的吸能組件安裝于某地鐵列車底架前端,其中一個吸能單元為雙錐內嵌隔板矩形管,該吸能單元僅有一組對稱面具有錐度。圖2展示了雙錐內嵌隔板矩形管的具體結構及尺寸。圖2(a)為雙錐內嵌隔板矩形管的內部結構,其由一個外管、兩個端板和五個隔板組成,外管的高度為146 mm,長度為639 mm,被五個隔板不均等劃分。

圖2(b)展示了各個隔板所在外管截面的寬度,且外管的一組對稱面具有錐度,錐度為2°。由于外管不同截面的寬度各不相同,圖2(c)中各個隔板W和W0的數值也各不相同。各個隔板W的數值等于其所在矩形截面的寬度數值,各個隔板W0的數值比其W的數值小130。此外,外管和隔板厚度均為3 mm,兩個端板的厚度為5 mm。

圖1 雙錐內嵌隔板矩形管在某地鐵列車的應用Fig.1 The application of DTRTD in a certain subway train

(a) 內部結構

(b) 俯視圖(c) 隔板圖2 雙錐內嵌隔板矩形管的結構及尺寸Fig.2 Detailed geometry of DTRTD

1.2 有限元模型

1.2.1 有限元模型的定義

為了研究雙錐內嵌隔板矩形管的耐撞性,本文使用非線性有限元軟件 LS-DYNA進行準靜態壓縮數值仿真分析,并將雙錐內嵌隔板矩形管與傳統矩形管進行對比。其中,傳統矩形管為各面無錐度、內部未嵌入隔板的矩形截面管,其總長度、高度、外管厚度以及端板厚度均與雙錐內嵌隔板矩形管相同,寬度為194 mm即雙錐內嵌隔板矩形管前后兩端寬度的平均值。

如圖3所示,有限元模型包括剛性墻、連接板、吸能結構、壓頭四個部分。兩塊連接板焊于吸能結構兩端,左側連接板與剛性墻接觸,右側連接板與壓頭接觸,設置剛性墻固定不動,壓頭沿軸向對右側連接板加載使吸能結構壓縮。為了準確模擬吸能結構的變形,吸能結構和連接板采用四邊形殼單元,單元厚度方向采用五點積分,面內采用單點積分,剛性墻和壓頭采用六面體實體單元。帶連接板的吸能結構的自接觸采用“Automatic_Single_Surface”算法,與剛性墻和壓頭的接觸均采用“Automatic_Surface_To_Surface”算法。接觸的靜摩擦因數設為0.3,動摩擦因數設為0.2[18],使用第四種沙漏黏度類型來控制沙漏能量。通過采用不同邊長的四邊形殼單元觀察建模時吸能結構EA的變化,發現當邊長小于或等于7 mm時EA趨向于平穩。因此,為減少計算時間且保證計算精度,將四邊形殼單元的邊長設置為7 mm 來對吸能結構和連接板進行建模。

圖3 有限元模型Fig.3 Finite element model

1.2.2 材料模型

為獲得吸能結構材料的力學性能,如圖4所示,使用MTS 647液壓萬能材料試驗機對加工好的標準試驗件進行準靜態拉伸試驗,表1和圖5(a)分別為測得的材料參數和工程應力-應變曲線。采用分段線性塑性材料模型Mat.024定義吸能結構的材料屬性,Mat.024可錄入由材料真實應力應變曲線所得到的等效塑性應力及應變,使材料屬性定義更加精確;連接板與吸能結構的材料相同,故同樣采用Mat.024材料模型;剛性墻和壓頭使用Mat.020剛性材料模型定義其屬性。

圖4 材料拉伸試驗Fig.4 Material tensile test

表1 吸能結構材料參數

真實應力應變曲線可分為兩段,可通過式(1)~(5)進行求解[19]。標準試驗件頸縮之前的部分為第一段,真實應變εt和頸縮之前的真實應力σt可如下求得

εt=ln(εe+1)

(1)

σt=σe(εe+1)

(2)

式中:σe為工程應力;εe為工程應變。

(a) 工程應力-應變曲線

(b) 真實應力-應變曲線圖5 材料應力-應變曲線Fig.5 Stress-strain curves of material

標準試驗件頸縮之后的部分為第二段,頸縮之后的真實應力σt求解如下

(3)

式中:n為應變硬化指數。C和n可如下求得

C=σu(e/n)n

(4)

n=ln[1+1/(0.24+0.013 95σu]

(5)

式中:σu為極限強度;e為自然常數。圖5(b)為最終求得的真實應力應變曲線。

2 試驗驗證

2.1 耐撞性的評價指標

吸能量EA(energy absorption)是衡量吸能結構耐撞性的關鍵指標,可以由吸能結構壓縮過程的力-位移曲線積分而得[20]

(6)

式中:F(s) 為瞬態撞擊力,是瞬態位移s的函數。

平臺力MCF為吸能結構變形過程中的平均壓縮力,可以通過吸能量EA 除以壓縮位移S得到

MCF=EA/S

(7)

比吸能SEA用來表示吸能結構單位質量的吸能量,可以通過吸能量EA除以吸能結構的總質量m得到[21]

SEA=EA/m

(8)

峰值力PCF對碰撞事故發生時乘員所受傷害有著重要影響,在數學上定義如下

PCF=max[F(s)]

(9)

2.2 準靜態軸向壓縮試驗

雙錐內嵌隔板矩形管的準靜態軸向壓縮試驗配置如圖6所示,由剛性墻、試件、液壓加載裝置以及多種測試儀器構成。試件加工過程中,其外管由垂向軸對稱的兩部分焊接而成,隔板和連接板均與外管焊接。液壓加載裝置受電磁系統控制,可提供穩定而均勻的位移或壓力。此試驗使用位移控制模式,以大約20 mm/min的速度對試件進行準靜態軸向加載,從而確保不會有動態影響產生[22],壓縮總行程設定為400 mm。測力傳感器安裝在剛性墻上,實時采集試件所受壓力;位移計測量剛性墻與壓頭的間距變化,進而得到試件的實時壓縮量。多臺數碼攝像機從不同角度記錄試件的準靜態軸向壓縮過程。

圖6 準靜態軸向壓縮的試驗配置Fig.6 Test configuration of quasi-static axial compression

2.3 結果分析及驗證

將雙錐內嵌隔板矩形管的數值仿真結果與試驗數據進行對比,從而分析其吸能特性并驗證有限元模型的準確性。圖7(a)將數值仿真所得雙錐內嵌隔板矩形管的力-位移曲線和試驗進行了比較,可明顯看出兩者具有較高的吻合度。由于現場試驗情況的復雜性,數值仿真的力-位移曲線無法與試驗完全一致。雙錐內嵌隔板矩形管的外管在整個壓縮過程中吸收了80.12 kJ的能量,而其余部分僅吸收了3.04 kJ的能量(隔板和端板分別吸收了2.99 kJ和0.05 kJ),如圖7(b)所示。

(a) 仿真和試驗的力-位移曲線

(b) 外管和其余部分吸能量圖7 數值仿真與試驗結果曲線Fig.7 Result curves of numerical simulation and test

表2將雙錐內嵌隔板矩形管數值仿真和試驗的EA和PCF進行了對比,誤差均不超過2%。其中,雙錐內嵌隔板矩形管數值仿真的吸能量為83.16 kJ,而傳統矩形管數值仿真的吸能量為62.31 kJ。因此,當兩者壓縮至行程終點時,雙錐內嵌隔板矩形管較傳統矩形管而言,能夠吸收更多的能量。雙錐內嵌隔板矩形管和傳統矩形管的質量分別為16.1 kg和12.5 kg,可根據式(8)求得其比吸能分別為5.17 kJ/kg和4.98 kJ/kg,故雙錐內嵌隔板矩形管的比吸能高于傳統矩形管。

表2 仿真與試驗的耐撞性指標對比

圖8中前兩行對比了雙錐內嵌隔板矩形管在試驗和數值仿真過程中三個時間點的壓縮變形,兩者呈現出很好的一致性。由圖中可看出雙錐內嵌隔板矩形管在壓頭的軸向加載下,每相鄰的隔板與端板或相鄰隔板之間形成一個褶皺,褶皺自壓頭一側向左依次形成,呈現穩定而規律的變形模式。當壓縮至行程終點即400 mm時,雙錐內嵌隔板矩形管共形成了五個褶皺。圖8中第三行為傳統矩形管在數值仿真過程中的壓縮變形,可發現傳統矩形管的變形不如雙錐內嵌隔板矩形管緊湊,在壓縮至行程終點時傳統矩形管僅形成了四個褶皺。

從以上的對比分析可得出,無論是力-位移曲線和耐撞性指標,還是變形模式,雙錐內嵌隔板矩形管準靜態壓縮的試驗結果均驗證了有限元模型的準確性。相較于傳統矩形管,雙錐內嵌隔板矩形管變形更為緊湊,能夠在相同壓縮行程下形成更多的褶皺,具有更高的吸能量和比吸能,因此雙錐內嵌隔板矩形管多方面的耐撞性能均優于傳統矩形管。

圖8 試驗和數值仿真的變形過程對比Fig.8 Comparison of deformation processes betweenexperiment and numerical simulation

3 多目標優化

3.1 優化問題的提出

在耐撞性設計中,要求吸能結構具有盡可能高的SEA,因此優化的第一個目標為最大化SEA[23]。同時,PCF需要盡可能減小,以避免碰撞發生時減速度過大而導致乘員傷亡慘重,故優化的第二個目標為最小化PCF。此外,根據該地鐵車輛的能量吸收分配方案,在優化中約束EA須不低于75 kJ,以確保雙錐內嵌隔板矩形管在規定壓縮行程內吸收足夠的能量。以上多目標優化問題在數學上定義如下:

(10)

3.2 試驗設計

外管由有錐度的一組對稱面和無錐度的一組對稱面組成,有錐度的對稱面的厚度定義為設計變量A,無錐度的對稱面的厚度定義為設計變量B,變化范圍均為1~5 mm。隔板的厚度定義為設計變量C,變化范圍為2~6 mm。三個設計變量的變化范圍覆蓋了此類吸能結構在實際工程中的常用厚度區間[24]。EA、SEA和PCF被設為三個輸出響應,壓縮行程與試驗和仿真相同,設定為400 mm。

試驗設計基于拉丁超立方采樣而建立。拉丁超立方采樣是一種從多元參數分布中近似隨機抽樣的方法,屬于分層抽樣技術,常用于計算機實驗或蒙特卡洛積分等。在統計抽樣中,拉丁方陣是指每行、每列僅包含一個樣本的方陣。拉丁超立方則是拉丁方陣在多維中的推廣,每個與軸垂直的超平面最多含有一個樣本。表3列出了應用拉丁超立方采樣而獲得的三組設計變量的不同取值組合,并基于已驗證的有限元模型求出對應的輸出響應。

表3 試驗設計矩陣

3.3 代理模型與誤差分析

由于實際工程中復雜產品的物理模型往往非常復雜,設計變量與目標性能之間通常不具有顯式的函數關系式,且表現為多參數、高維數、強非線性問題。為了能高效、準確地獲得優化結果,利用代理方法對離散的數據進行擬合,在不降低模型精度情況下建立高效的模型來替代實際模型,這類模型也稱作是代理模型。徑向基函數是一種高效的代理方法,尤其是在對SEA的預測上,擁有良好的精確性,故被選用于本文以建立代理模型。

平均相對誤差ARE(average relative error),最大相對誤差MRE(maximum relative error)和確定系數R2用于衡量代理模型的精確度,在數學上定義如下

(11)

(12)

(13)

表4 代理模型的誤差分析結果

3.4 參數研究

分別以三個設計變量中的一個作為自變量,將其余兩個設計變量設定為各自變化范圍的中間值,根據代理模型得到各輸出響應的變化情況;再將各個自變量按比例換算為厚度級別,其中自變量的中間值的厚度級別設為0,最小值和最大值分別設為-1和+1,從而得到圖9所示的輸出響應與設計變量的關系曲線圖。為得到各輸出響應的三維云圖,將一個設計變量設定為其變化范圍的中間值,另外兩個設計變量均作為自變量,根據代理模型構建圖10中的各響應面。

(a) 設計變量與吸能量(b) 設計變量與比吸能(c) 設計變量與峰值力

由圖9(a)可見:EA隨設計變量A和B的增加而單調遞增;EA隨著設計變量C的增加,先增長而后略微有下降的趨勢。由圖10(a)可看出,在B設定為3 mm時,EA隨A的增加而增加,隨C的增加而先增后減;在A為5 mm且C大約為4 mm時,EA達到最大。

與EA相似,圖9(b)中SEA隨設計變量A和B的增加而單調遞增;SEA隨著C的增加而先增后減。由圖10(b)可看出,在A設定為3 mm時,SEA隨B的增加而增加,隨C的增加而先增后減;在B為5 mm且C大約為3.5 mm時,SEA達到最大。

PCF隨三個設計變量的增加均保持單調遞增,如圖9(c)所示,但C增長至5.5 mm附近時,PCF的增長趨勢明顯放緩而趨于平穩。由圖10(c)可看出,在C設定為4 mm時,PCF隨A和B的增加而增加,在A和B均為5 mm時,PCF達到最大。

從圖9可看出,根據曲線斜率的絕對值,設計變量對EA、SEA以及PCF的影響程度從大到小依次排序均為B、A、C。

3.5 優化算法

遺傳算法是當前較為流行的一種優化方法,派生于自然選擇法則和生物進化的遺傳機制。遺傳算法可以確定空間里的所有解,并且不會在局部優化時陷入快速衰退的陷阱,具有很強的全局搜索能力。本文所采用的多目標遺傳算法是為解決多目標優化問題而對遺傳算法進行的擴展。在有多個目標時,由于存在目標之間的沖突和無法比較的現象,一個解在某個目標上是最好的,在其他的目標上可能是最差的。這些在改進任何目標函數的同時,必然會削弱至少一個其他目標函數的解稱為非支配解或帕累托解,一組目標函數最優解的集合稱為帕累托前沿。在多目標優化問題中,有不止一個目標函數需要最大化或者最小化,因此最終目的不是尋求單個最優解而是尋求帕累托前沿。表5羅列了本文所使用多目標遺傳算法的參數定義,圖11為多目標優化算法的流程圖。

(a) 吸能量三維云圖(b) 比吸能三維云圖(c) 峰值力三維云圖

表5 多目標遺傳算法的參數定義

圖11 多目標遺傳算法流程圖Fig.11 Flowchart of multi-objective genetic algorithm process

3.6 結果驗證及討論

多目標優化所得的帕累托前沿如圖12所示,雙錐內嵌隔板矩形管的比吸能與峰值力呈正相關,故最大化比吸能和最小化峰值力這兩個目標具有互斥性[15]。盡管如此,圖中圈內的點相較于優化前結果而言,均實現了更高的比吸能以及更低的峰值力,而圈外的其余點所對應的設計方案則是對兩個目標之一起到了明顯的優化效果。

圖12 多目標優化帕累托前沿Fig.12 Pareto front of MOO

對于帕累托前沿中具體設計方案的選擇,則可根據實際工程中的需求對各目標進行權重分配,式(14)中的權函數取最大值時所對應的設計方案即為該權重分配準則下的最優解。

f(x)=

(14)

式中:SEA(x)和PCF(x)為帕累托前沿中各設計方案對應的SEA和PCF的值;SEAU和SEAL分別為試驗設計中SEA的上限和下限;PCFU和PCFL分別為試驗設計中PCF的上限和下限;W1和W2分別為分配給SEA和PCF的權重因子。

本文對三種權重分配準則的最優設計方案進行求取,第一種權重分配準則更偏向于提升SEA,故設定W1=0.7且W2=0.3;第二種權重分配準則無偏向性,故設定W1=0.5且W2=0.5;第三種權重分配準則更偏向于降低PCF,故設定W1=0.3且W2=0.7。求取的各權重分配準則對應的最優設計方案如表6所示,并根據各方案對應的設計變量取值進行數值仿真,結果表明三種權重分配準則所對應最優設計方案的比吸能和峰值力的預測誤差均小于百分之四,從而驗證了代理模型的準確性。

表6 各權重分配準則的最優設計方案及驗證

4 結 論

本文提出一種用于某地鐵列車碰撞吸能的雙錐內嵌隔板矩形管。通過準靜態軸向壓縮試驗驗證有限元分析的準確性,采用拉丁超立方采樣進行試驗設計,根據徑向基函數獲得代理模型,從而分析各變量對響應的影響,并以最大化比吸能和最小化峰值力為目標,對雙錐內嵌隔板矩形管進行多目標優化設計。研究結果表明:

(1) 通過數值仿真對雙錐內嵌隔板矩形管進行耐撞性分析,所得到的力-位移曲線、耐撞性響應以及變形模式等結果均與準靜態軸向壓縮試驗有著良好的一致性,驗證了所建立的有限元模型的準確性。雙錐內嵌隔板矩形管在準靜態軸向壓縮下呈現出規律而穩定的變形,具有較為理想的耐撞性響應。

(2) 對比雙錐內嵌隔板矩形管與傳統矩形管準靜態壓縮的數值仿真結果,發現雙錐內嵌隔板矩形管較傳統矩形管而言,變形更為緊湊,在相同壓縮行程下能形成更多褶皺,吸能量和比吸能均更高,故雙錐內嵌隔板矩形管具有更優的耐撞性能。

(3) 各設計變量對吸能量、比吸能以及峰值力的影響程度從大到小依次排序均為無錐度外管對稱面厚度、有錐度外管對稱面厚度、隔板厚度。其中無錐度外管對稱面厚度和有錐度外管對稱面厚度對三個響應均產生正相關的影響。隨著隔板厚度的增加,吸能量和比吸能先增加而后呈現減小的趨勢,峰值力先增加而后趨于穩定。

(4) 多目標優化后得到的帕累托前沿表明,最大化比吸能和最小化峰值力這兩個優化目標具有互斥性,因而無法同時達到最優解。盡管如此,帕累托前沿中的部分設計方案相較于優化前而言,能夠同時實現更高的比吸能以及更低的峰值力。根據實際工程中的不同需求,可為各目標進行相應的權重分配,并通過求取權函數最大值,能夠實現不同權重分配準則的最優設計方案的選擇。因此,本文多目標優化有效提升了雙錐內嵌隔板矩形管的耐撞性,具有較大的工程價值。

猜你喜歡
模型設計
一半模型
重要模型『一線三等角』
何為設計的守護之道?
現代裝飾(2020年7期)2020-07-27 01:27:42
重尾非線性自回歸模型自加權M-估計的漸近分布
《豐收的喜悅展示設計》
流行色(2020年1期)2020-04-28 11:16:38
瞞天過海——仿生設計萌到家
藝術啟蒙(2018年7期)2018-08-23 09:14:18
設計秀
海峽姐妹(2017年7期)2017-07-31 19:08:17
有種設計叫而專
Coco薇(2017年5期)2017-06-05 08:53:16
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 国产大片喷水在线在线视频| 激情无码字幕综合| 日韩欧美中文在线| 制服丝袜一区二区三区在线| 日韩欧美国产另类| aa级毛片毛片免费观看久| 中文字幕免费视频| 99热最新网址| 40岁成熟女人牲交片免费| 91偷拍一区| 国产欧美日韩一区二区视频在线| 欧美日韩在线国产| 亚洲欧美另类中文字幕| 国产精品久久自在自线观看| 欧美五月婷婷| 欧美成人午夜影院| 丰满人妻一区二区三区视频| 欧美第一页在线| 视频二区中文无码| 国产人成在线视频| 老司国产精品视频91| 好紧好深好大乳无码中文字幕| 不卡无码网| 国产精品短篇二区| 黄色网在线| 国产精品999在线| 美女啪啪无遮挡| 国产簧片免费在线播放| 国产精品太粉嫩高中在线观看| 在线日韩日本国产亚洲| 亚洲第一中文字幕| 玖玖精品在线| 无码视频国产精品一区二区 | 免费啪啪网址| 黄色国产在线| 亚洲第一区在线| 国产成人乱无码视频| 国内精品免费| 日韩经典精品无码一区二区| 婷婷激情亚洲| 一级黄色片网| 成年人国产网站| a色毛片免费视频| 日本亚洲成高清一区二区三区| 国产资源免费观看| 亚洲精品无码久久久久苍井空| 99精品欧美一区| 欧美另类视频一区二区三区| 在线观看免费AV网| 国产精品亚欧美一区二区| 成人一级黄色毛片| 粗大猛烈进出高潮视频无码| 国产亚洲欧美在线视频| 色综合中文字幕| 亚洲侵犯无码网址在线观看| 激情在线网| 久久精品这里只有精99品| 一区二区三区四区精品视频 | 国产精品流白浆在线观看| 日韩精品无码免费一区二区三区 | 无码aⅴ精品一区二区三区| 日韩在线2020专区| 国产高清色视频免费看的网址| 亚洲精品国产综合99| 女人天堂av免费| 亚洲综合一区国产精品| 99ri国产在线| 亚洲国产第一区二区香蕉| 中文纯内无码H| 亚洲中文字幕无码爆乳| 不卡无码h在线观看| 天天综合网在线| 男人天堂伊人网| 亚洲不卡影院| 亚洲欧美不卡| hezyo加勒比一区二区三区| 亚洲丝袜中文字幕| 婷婷丁香色| 亚洲va视频| 精品国产三级在线观看| 伊人久久福利中文字幕| 亚洲AⅤ永久无码精品毛片|