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

液壓模塊掛車剛柔耦合動態響應特性及擺臂疲勞分析

2014-09-05 10:11:42武俊達董大偉王媛文孫玉華
振動與沖擊 2014年7期
關鍵詞:模態有限元模型

武俊達, 董大偉, 閆 兵, 王媛文, 孫玉華

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

液壓模塊掛車是一種特大、特重貨物陸路運輸工具,該車主要由4部分組成:結構、機構、液壓系統與電氣控制系統。其中車架屬于大型空間結構,由多模塊拼接而成,掛車被牽引行駛過程中受到來自路面不平度的振動激勵時,車架彈性體的模態特征與剛柔耦合現象直接影響整體系統動態響應。液壓系統與懸架結構結合成為液壓懸掛系統,其用于支撐車架與貨物載重、實現同油缸組的車軸載荷相等、調節車架高度以保持車架平面裝載水平等功能。實際運行中,懸架擺臂已多次發生斷裂故障,對掛車運行安全性與運輸效率有很大的影響。

目前,關于懸架結構強度及疲勞方面的研究主要側重于有限元靜態分析方法[1-2],忽略動態交變載荷的影響。為了實現車軸承受載荷相等,可采用平衡桿系模擬液壓懸掛系統[3],而所述的平衡桿系不能全面地反映液壓系統的作用狀態對懸架結構動載荷的影響。由于掛車具有多體系統動力學復雜的特點,已有學者把該車簡化為一縱六軸線式,建立了多剛體系統動力學模型并進行液壓懸掛系統動載荷分析[4]。因此,運用有限元分析與剛柔多體動力學仿真的相結合方法,應用于液壓模塊掛車動態響應及其擺臂疲勞強度研究具有工程意義與學術價值。

基于有限元分析法與多體動力學理論,運用ANSYS與ADAMS之間接口技術,建立15軸線式液壓模塊掛車剛柔耦合多體動力學模型,其中把車架看為柔性體,結合掛車振動試驗驗證模型的準確性。對掛車在不同路面等級與車速,考慮凹凸路面沖擊下的運行工況,進行仿真計算并分析掛車動態響應特性。以多體動力學仿真計算結果為有限元分析的載荷,計算獲得擺臂結構應力時間歷程。選擇合理的疲勞分析方法,進行擺臂疲勞壽命計算,包括:應力譜的循環計數、損傷累積及壽命預測。

1 基礎理論

1.1 有限元模態分析的縮減自由度Guyan法

含n自由度結構有限元模型的動力平衡方程為:

(1)

式中,M、C與K為質量阻尼與剛度矩陣;u是節點位移矢量;f是隨時間變化的載荷矢量。

縮減自由度Guyan法[5]的理論是將u分為主自由度um與從自由度us對應的位移向量,n=m+s。當沒有載荷施加于從自由度上,且阻尼可忽略時,得到:

(2)

考慮靜態情況時,式(2)可寫為:

(3)

由式(3),可得到:

{us}=-[Kss]-1[Ksm]{um}

(4)

[[Kmm]-[Kms][Kss]-1[Ksm]]{um}={fm}

(5)

(6)

式中,T為靜態轉換矩陣;I為m×m階單位矩陣。縮減到主自由度上的剛度矩陣為:

(7)

用相同方法,將轉換矩陣T擴展到式(2)中,得到縮減到主自由度上的質量矩矩陣與載荷矢量:

(8)

方程(2)可變為:

(9)

其自由度數縮減為主自由度數,比式(2)減少s個自由度。

1.2 柔性體動力學的模態綜合Craig-Bampton法

ADAMS中采用模態綜合技術Craig-Bampton法[6]求解柔性體運動方程。該方法將自由度系統分成邊界自由度與內部自由度,對應兩種模態:約束模態和固定邊界正常模態。柔性體變形矢量U可用兩種模態的物理位移及其對應模態坐標來表示,如公式(2)表示:

(10)

式中,UB、UI為邊界與內部自由度的變形矢量;0為零矩陣;ΦIC,ΦIN為約束模態與正常模態中內部自由度的物理位移;qC、qN為約束模態和正常模態的模態坐標矢量;Φ成為模態轉換矩陣。

總體剛度和質量矩陣轉換可通過模態轉換矩陣獲得的:

(11)

(12)

動力平衡方程由模態坐標矢量通過模態轉換矩陣表示如下:

(13)

(14)

在此基礎上,求解柔性體模型在廣義坐標的位移、速度與加速度向量和通過Lagrange方程式建立動力學微分方程[7]。考慮系統沖擊時可采用傳遞矩陣法的剛柔耦合多體系統的沖擊響應分析方法[8],將沖量模型和碰撞接觸模型與多體系統傳遞矩陣法結合,對系統進行動態響應分析。

1.3 液壓系統的油液運動狀態方程描述

目前,工程流體分析最廣泛應用Lagrange方法描述流體的運動狀態,其根據質量守恒定律、能量轉化與守恒定律等及其油液的物理屬性,建立油液力學微分方程組[9]。在此基礎上,ADAMS/Hydraulic模塊中采用Merrit方法的Navier-Stockes油液運動方程、液體連續方程、能量累積方程等來描述油液的運動狀態方程(Equation of state for liquid)[10]

ADAMS中,基于Lagrange方程式,含液壓系統的多體系統廣義坐標由兩個部分組成:第一部分是決定物體位形(位置和方向)和油缸的伸縮量;第二部分是決定油液狀態的油液壓力。系統的廣義力除了外力和慣性力,還存在在油液體積變遷時對外做功。因此含多油缸聯通的多體系統動力學微分方程包括:多體系統動力學方程[7]和油液的動力學方程。

2 液壓模塊掛車剛柔多體動力學模型的建立

2.1 車架有限元模型的建立

運用ANSYS建立車架有限元模型時,考慮鞍座支架對車架動態特性的影響,把鞍座支架看為車架的一部分。采用板殼單元Shell63與實體單元Solid45劃分網格,結合剛性多節點約束方法來處理二者間的連接處。多體動力學模型中,柔性車架與其他剛體連接通過外部節點,該節點采用多點約束與集中質量單元模擬。材料屬性設為普通合金結構鋼16Mn。

由上述的縮減自由度Guyan法,對車架有限元模型進行模態求解并生成ADAMS所需要的模態中性文件。該方法的準確度將取決于車架有限元模型主自由度的合理選取,根據結構特點與外部節點數量,選取288階車架主自由度模態、20階固有頻率模態。

2.2 柔性車架的處理

多體動力學模型的計算精度將取決于選取柔性體模態階數。根據車架振動試驗分析結果,車架模態頻率主要集中在范圍50Hz以內。因此,除了前6階為剛體模態,對柔性車架選取20階模態。與全自由度模態Block-Lanczos法的計算結果對比如表1所示,可看出:選取ANSYS中有限元模型的主自由度與ADAMS中柔性體的模態階數能夠反映車架動態特征,滿足模型的計算精度。

模型采用鞍座支架的裝載方式,為了實現柔性車架的鞍座支架與剛體貨物間的連接,將大剛度的彈簧力元放置在二者間[11]。

2.3 液壓回路系統的建立

液壓懸掛包括懸臂Xi-j、擺臂Si-j、油缸Ci-j等主要部件,如圖1所示,i=1, 2,…, 15為掛車的軸線,j=1, 2, 3表示每軸線橫向有3個液壓懸掛組。

圖1 液壓懸掛組的主要部件

采用ADAMS/Hydraulic模塊建立多油缸聯通的液壓回路系統時,假設;油缸體、活塞桿與油管是剛性體;忽略液壓回路系統中的安全閥、溢流閥、換向閥等液壓元件和熱能效應。油缸間連接通過油管,油缸與油管的工作內涵特性取決于選取油管類型和幾何參數。設計三點支撐方式的液壓回路系統如下:前6個軸線的18個油缸互相連通成為一組,后9個軸線的27個油缸縱橫交錯連通而分成左右兩組。選取油缸類型為Cylinder1f,其內徑D=0.16 m;油管類型為Pipe_level1,其內徑d=0.007 m;油液的體積彈性模量Eh=1 800 MPa、密度ρ=900 kg/m3[10]。

表1 縮減法與完全自由度法的固有頻率對比

2.4 輪胎等效剛度計算

采用線性等效剛度KT的彈簧力元模擬輪胎,由車輪承受載荷Fw與輪胎變形量δ確定。用Komandi[12]提出經驗公式計算輪胎變形量,每組輪軸有四個8.25R15型輪胎,計算得到KT=6.5×102kN/m,選取輪胎阻尼系數CT=1.5 kN·s/m。

2.5 路面不平度的時域激勵加載

多體動力學模型采用激振臺加載形式,以路面不平度為激勵信號。大量試驗測試結果表明:路面不平度的數據特性是平穩的,是具有零均值與各態歷經的Gauss隨機過程,可用功率譜密度來描述路面的統計特性。基于GB7031-1986《車輛振動輸入——路面平度表示》標準[13]以及路面不平度的Fourier逆變換數值模擬方法[14,15],考慮車速的影響,選取路面等級為B級、C級與D級;車速分別為2.22m/s、3.33 m/s、4.44 m/s與5.55 m/s;時間頻率范圍為0.5Hz≤f≤30 Hz;時間間隔為0.01 s,采樣數據點數為212,仿真時長為40.96 s。運用MATLAP編制程序,計算獲得路面不平度數據,ADAMS中用樣條插值函數AKIPSPL編輯數據成為激振臺的位移加載曲線。

圖2 正弦形凹凸路面的幾何

表2 正弦形凹凸路面的參數

為了描述凹凸不平度的復雜道路,采用正弦形曲線來實現。正弦形幾何參數由正弦形幅值h、軸距L=1.6 m和正弦線凸起部分與前后車輪外公切線來確定,如圖2所示。在此基礎上,設計正弦形幾何參數如表2所示。多體動力學模型中,縱向中間車輪(j=2)由C級路面不平度曲線分別加上正弦形長度2T的S1與S2為激勵信號;縱向左右兩行車輪(j=1, 3)由C級路面不平度為激勵信號。

液壓模塊掛車剛柔多體動力學模型如圖3所示:

圖3 掛車剛柔多體動力學模型

3 掛車道路振動試驗與仿真計算結果的對比

通過對掛車進行道路振動試驗,獲得車架與擺臂的動態垂向加速度頻譜曲線,結合仿真計算結果來驗證剛柔耦合多體系統動力學模型的準確性。

3.1 掛車振動試驗分析

試驗掛車裝貨物載重為282噸如圖4所示,采用六個傳感器單向加速度測量式:測點1到測點5位置分別在車架平面裝載如圖5所示;測點6位置在擺臂S5-2的軸頭處。掛在水泥混凝土路面上運行,穩定車速為2.22 m/s。提取振動測試信號及其變換得到測試點的頻譜曲線,如圖6所示。

圖4 15軸線式液壓模塊掛車

圖5 車架各測點位置

3.2 振動試驗與仿真計算結果的對比分析

掛車在B級路面、車速為2.22 m/s下的運行,仿真計算得到相應于各測點垂向加速度頻譜曲線如圖7所示。

從圖6可看出,各測點都存在2.4 Hz的共振頻率,表明該頻率是實際整車振動頻率,相應仿真計算結果出現在2.6 Hz的共振頻率(加速度振幅為0.02 m/s2左右,如圖7所示)。在ADAMS中運用線性模態分析功能對整車模型進行分析計算,得出在2.6 Hz是整車垂向平動模態振型。

各測點除了在2.4 Hz的共振頻率,還存在比較明顯的共振峰頻率如測點1、2、3與4分別在9.9 Hz、10.45 Hz、10.45 Hz與9.6 Hz;測點5在17.4 Hz;懸架擺臂的測點6在25.2Hz。這些共振峰是掛車各組成部分的固有頻率,與仿真計算結果對比偏差較小。

分析結果表明:仿真計算得到相應于振動試驗測試點的頻率譜變化規律相同,模態頻率值基本一致,所建立的模型具有一定的準確性。

圖6 試驗測點的加速度頻譜

圖7 仿真計算的加速度頻譜

4 掛車多體動力學仿真及其動態響應分析

4.1 掛車在不同路面等級與車速下的運行

掛車滿載為500 t,裝載方式通過三個鞍座支架,在不同路面等級與車速下的運行工況進行仿真計算,結果表明:擺臂S8-2質心垂向加速度均方根(RMS)值與地面承受車輪W8-2最大動載荷值隨路面不平度與車速增加而增大,如圖8與圖9所示。

車輪Wi-j動載荷RMS值如圖10所示,可看出:同屬于一個油缸組的動載荷RMS值相對偏差較小。

柔性車架與懸臂Xi-j連接點的垂向加速度RMS值如圖11所示。由于車架其他部位不受鞍座支架剛度的影響,其垂彎剛度小于車架與鞍座支架接觸處。車輪受來自路面不平度激勵時,振動能量轉換為油液流動動能,其瞬時轉化為油缸內的油液彈性勢能。油缸通過對外做功,將其彈性勢能轉換為車架結構的彈性勢能并使車架變形,車架與懸臂連接處的剛度較小部位會產生垂向變形量較大。因此,車架兩邊與前后車頭部位的加速度向量大于縱向對稱面與接近支架部分(第3、4、8、12與第13軸線)。

圖8 擺臂S8-2質心加速度RMS值

圖11 車架與懸臂連接點的加速度RMS值

圖12 擺臂質心加速度值

圖13 車架與懸臂X8-2連接點的加速度

4.2 掛車在正弦形凹凸路面沖擊下的運行

掛車在正弦形凹凸路面S2類、車速為2.22 m/s下的運行,擺臂S8-2質心垂向加速度如圖12所示,輪胎與凹凸沖擊時加速度值很明顯劇烈增加。

由于液壓懸掛系統能調整車架高度以保持車架水平,因此當任意車輪Wi-2受到正弦形凹凸路面沖擊時,柔性車架的加速度也呈正弦形變化,如圖13所示柔性車架與懸臂X8-2連接點的垂向加速度變化過程。

5 擺臂結構疲勞壽命預測

掛車運行中,擺臂已多次發生斷裂故障,其斷裂斷面位置如圖14所示。將作用于擺臂的動載荷,包括慣性載荷信息,轉換成適用于有限元分析的數據格式,ANSYS中采用準靜態疊加法[16]在離散的時間點上進行分析計算。根據擺臂結構的應力時間歷程結果,選擇合理的疲勞分析方法,對擺臂進行疲勞壽命預測。

圖14 擺臂斷裂故障的斷面位置

5.1 擺臂結構有限元分析

選擇對應車輪有相對較大動載荷的擺臂S8-2為疲勞分析的對象。采用實體單元Solid45進行擺臂劃分網格,材料屬性為鑄鐵QT450-10。提取作用于擺臂S8-2的載荷時間歷程為1 000步數,進行有限元準靜態疊加分析計算,獲得擺臂有限元等效應力(Von-mises)分布如圖15所示。

圖15 擺臂有限元等效應力分布圖

圖16 擺臂最大等效應力值

圖17 應力范圍與平均值分布

擺臂應力集中部位出現在兩邊側翼與底板的交接處,最大等效應力值σmax隨路面不平度與車速增加而增大,其均大于材料屈服極限σS=331 MPa[17],如圖16所示。這表明危險部位應力水平已進入塑性狀態。在D級路面、車速為5.55 km/h下的運行工況,,運用雨流計數法對危險點的應力時間歷程進行應力統計分析,獲得應力范圍與平均值分布如圖17所示。

5.2 疲勞壽命估算方法

根據有限元分析結果,危險點應力水平已進入塑性狀態,其應力與應變不再是線性關系,塑性應變成為影響其疲勞壽命的主要因素。基于局部應力應變法,考慮平均應力的影響,采用Manson-Coffin-Baquin的Morrow平均應力修正公式[18]:

(15)

5.3 Miner線性累積損傷理

運用Miner線性累積損傷理論[19],求解結構疲勞損傷及其疲勞壽命。首先求得在第i應力水平的應力幅范圍Δσi與循環次數ni造成損傷Di,壽命為Ni,總累積損傷D等于1時結構發生疲勞破壞。Miner理論可用式(16)來描述

(16)

5.4 擺臂疲勞壽命預測

圖18 擺臂疲勞壽命分布圖

掛車在不同路面等級與車速下的運行,危險點的疲勞壽命如圖19所示。若掛車使用系數為0.5,最小循環次數Nfmin=4.75×107,即壽命為30.1年,大于掛車使用年限為15年(包括,報廢年限為10年和可延緩報廢為5年)。這表明掛車在B級、C級與D級普通路面下的運行,擺臂不會出現疲勞斷裂。

圖19 掛車在普通路面下的運行,擺臂危險點的疲勞壽命

圖20 掛車在凹凸路面下的運行工況,擺臂危險點的疲勞壽命

掛車在正弦形凹凸路面下的運行,擺臂危險點的疲勞壽命隨正弦線幅值h與車速增加而明顯減短,如圖20所示。計算結果符合物流公司的運行記錄:擺臂斷裂故障常發生在凹凸不平的復雜路段。在凹凸路面S2類,若車速減小為1.67 m/s與1.11 m/s,危險點疲勞壽命可延長分別為18.83年與29.79年。

根據仿真計算結果,各運行工況可分為:適用的運行工況,即擺臂危險點的疲勞壽命大于掛車使用年限的;其余是不宜使用的運行工況,其疲勞壽命小于掛車使用年限。

6 結 論

(1) 基于有限元分析與剛柔多體動力學理論的相結合方法,建立了液壓模塊掛車剛柔耦合多體動力學模型。通過振動試驗與仿真計算結果的對比分析表明:所建立的模型具有一定的準確性。掛車的部件垂向加速度與部件之間動載荷隨路面不平度與車速增加而增大。模型中的液壓回路系統能夠有效地實現液壓懸掛系統的基本功能。

(2) 根據結構有限元分析與疲勞壽命預測的計算結果得出:擺臂疲勞危險部位出現在已發生斷裂的斷面位置。掛車在B級~D級普通路面下的運行,危險點的疲勞壽命均大于掛車使用年限。而在正弦形凹凸路面沖擊下,正弦形幅值與車速增加是引起擺臂危險點的疲勞壽命急劇減短。計算結果符合物流公司的運行記錄:擺臂斷裂故障常發生在凹凸不平的復雜路段。

(3) 根據仿真計算結果,可提出適用的運行工況以提高掛車的運行安全性與運輸效率。研究方法與成果可為類似車輛的結構設計、優化及工程應用提供了參考依據。后續工作為對擺臂進行結構優化研究,改善危險部位的應力水平,并提高結構抗疲勞強度。

參 考 文 獻

[1]張衛東, 莫旭輝, 彭勁松. 液壓平板車懸架系統的有限元分析與結構優化[J]. 機械工程師, 2007, 10: 34-36.

ZHANG Wei-dong, MO Xu-hui, PENG Jin-song. Simulation and structural optimization about suspension of hydraulic conveyance vehicle[J]. Mechanical Engineer,2007,10:34-36.

[2]涂啟養, 董大偉, 閆兵, 等. 某液壓模塊式組合掛車懸架結構件的強度及疲勞損傷分析[J]. 汽車科技, 2011, 3: 54-56.

TU Qi-yang, DONG Da-wei, YAN Bing, et al. Strength and fatigue analysis on the hydraulic modular assembled trailers suspension structure[J]. Auto Mobile Science & Technology, 2011, 03: 54-56.

[3]張宇探,馬力,李冰. 液壓模塊組合掛車整體結構有限元計算分析[J]. 專用汽車, 2010, 1: 53-55.

ZHANG Yu-tan, MA Li, LI Bing. Calculation and analysis by FEM for whole structure of hydraulic module assembled trailer[J]. Special Purpose Vehicle, 2010, 1: 53-55.

[4]陳黝生. 組合掛車液壓懸掛動載荷分析[D]. 上海: 同濟大學工學, 2007.

[5]Guyan R J. Reduction of stiffness and mass matrices[J]. AIAA Journal, 1965, 3(2): 380.

[6]Craig R J, Bampton M. Coupling of substructures for dynamic analyses[J]. AIAA Journal, 1968, 6(7): 1313-1319.

[7]馬星國, 尤小梅, 聞邦椿. 基于虛擬樣機技術的曲軸多體動力學仿真[J]. 振動與沖擊, 2008, 27(9): 155-157.

MA Xing-guo, YOU Xiao-mei, WEN Bang-chun. Multi-body dynamics simulation for a crankshaft system based on vitual prototype technology[J]. Journal of Vibration and Shock, 2008, 27(9): 155-157.

[8]隋立起, 鄭鈺琪, 王三民. 剛柔耦合多體系統的沖擊響應分析方法及應用研究[J]. 振動與沖擊, 2012, 31(15): 26-29.

SUI Li-qi, ZHENG Yu-qi, WANG San-min. Method of analysing impact response of multi-rigid-flexible system and its application[J]. Journal of Vibration and Shock, 2012, 31(15): 26-29.

[9]梁智權. 流體力學[M]. 重慶,重慶大學出版社, 2002.

[10]Herbert E. Merritt. Hydraulic Control Systems[M]. John Wiley & Sons Inc, 1967.

[11]謝基龍, 張燕, 謝云葉. 鐵路凹底平車凹底架動態響應及其疲勞強度[J]. 機械工程學報, 2010, 46(16): 16-22.

XIE Ji-long, ZHANG Yan, XIE Yun-ye. Dynamic response and fatigue strength of dePublisheded center flat frame[J]. Journal of Mechanical Engineering, 2010, 46(16): 16-22.

[12]莊繼德. 計算汽車地面力學[M]. 北京: 機械工業出版社, 2001.

[13]GB/T 7031-1986 車輛振動輸入與路面平度表示方法[S]. 北京: 中國標準出版社, 1987.

[14]段虎明, 石峰, 謝飛, 等. 路面不平度研究綜述[J]. 振動與沖擊, 2009, 28(9): 95-101.

DUAN Hu-ming, SHI Feng, XIE Fei, et al. A survey of road roughness study[J]. Journal of Vibration and Shock, 2009, 28(9): 95-101.

[15]劉獻棟, 鄧志黨, 高峰. 公路路面不平度的數值模擬方法研究[J]. 北京航空航天大學學報, 2003, 29(9), 843-846.

LIU Xian-dong, DENG Zhi-dang, GAO Feng. Research on the method of simulating road roughness numerically[J]. Journal of Beijing University of Aeronautics and Astronautics, 2003, 29(9): 843-846.

[16]Bishop N W M, Sherratt F. Finite element based fatigue calculations[M]. National Agency for Finite Element Methods and Standards NAFEMS Ltd Published, 2000.

[17]陸明炯. 實用機械工程材料手冊[M]. 沈陽: 遼寧科學技術出版社, 2004.

[18]David D W. Limits on morrow mean stress correction of manson-coffin life prediction models[C]. Proceedings of ASME Turbo Expo 2011, Vancouver, Canada, 2011: 85-92.

[19]李春祥, 李薇薇. 斜拉索風致動力疲勞損傷的研究[J]. 振動與沖擊, 2009, 28(11): 61-66.

LI Chun-xiang, LI Wei-wei. Analysis of wind-induced fatigue damage of inclined cables[J]. Journal of Vibration and Shock, 2009, 28(11): 61-66.

[20]Atzori B, Meneghetti G, Ricotta M. A compatible method to summarise the low-and high-cycle fatigue test results of ductile irons and structural steels[J]. Giornata IGF Forni di Sopra (UD), 2012, 153-166.

[21]聶宏, 喬新. 一種新的應變疲勞材質參數估算方法[J]. 南京航空學院學報, 1992, 24(2): 226-231.

NIE Hong, QIAO Xin. A new prediction method of material behavior parameters for strain fatigue[J]. Journal of Nanjing Aeronautical Institute, 1992, 24(2): 226-231.

[22]聶宏, 常龍. 基于局部應力應變法估算高周疲勞壽命[J]. 南京航空航天大學學報, 2000, 32(1): 75-79.

NIE Hong, CHANG Long. High-Cycle fatigue life prediction based on local stress-strain method[J]. Ournal of Nanjing University of Aeronautics &Astronautiics,2000,32(1):75-79.

猜你喜歡
模態有限元模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
國內多模態教學研究回顧與展望
基于HHT和Prony算法的電力系統低頻振蕩模態識別
磨削淬硬殘余應力的有限元分析
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
基于SolidWorks的吸嘴支撐臂有限元分析
箱形孔軋制的有限元模擬
上海金屬(2013年4期)2013-12-20 07:57:18
主站蜘蛛池模板: 国产成+人+综合+亚洲欧美| h视频在线观看网站| 国产成人精品高清不卡在线| 中国一级特黄视频| 91成人在线观看| 国产第一福利影院| 三上悠亚一区二区| 强乱中文字幕在线播放不卡| 国产91丝袜在线播放动漫 | 欧美亚洲欧美| 成·人免费午夜无码视频在线观看| 久久这里只有精品8| 久久99久久无码毛片一区二区| 国产真实乱了在线播放| 国产第一页亚洲| 白丝美女办公室高潮喷水视频| 青青极品在线| 中文字幕在线视频免费| 国产在线视频欧美亚综合| 美女一级毛片无遮挡内谢| 久久国产亚洲欧美日韩精品| 久久久黄色片| 久久综合色播五月男人的天堂| 欧美国产综合视频| 岛国精品一区免费视频在线观看 | 99er精品视频| P尤物久久99国产综合精品| 久久婷婷六月| 欧美日在线观看| 亚洲天堂视频在线观看免费| 999精品色在线观看| 国产精品久线在线观看| 四虎在线观看视频高清无码| 久久超级碰| 国产99免费视频| 欧美综合激情| 午夜精品影院| 极品私人尤物在线精品首页| 在线国产91| 国产午夜精品鲁丝片| 色综合中文综合网| 亚洲AV无码一区二区三区牲色| 国产美女无遮挡免费视频| 亚洲 欧美 偷自乱 图片| 夜精品a一区二区三区| 国产人免费人成免费视频| 在线看片免费人成视久网下载| 欧美精品综合视频一区二区| 亚洲成aⅴ人片在线影院八| 2021最新国产精品网站| 伊人久久精品亚洲午夜| 久久毛片基地| 国产精品久久国产精麻豆99网站| 国产噜噜噜视频在线观看| 亚洲精品色AV无码看| 99热这里只有精品久久免费| 极品性荡少妇一区二区色欲| 老司机精品99在线播放| 国产成人精品一区二区三区| 91毛片网| 91视频免费观看网站| 香港一级毛片免费看| 中文字幕无码制服中字| 国产肉感大码AV无码| 免费国产高清视频| 中文字幕亚洲电影| 精品久久久久久成人AV| 国产香蕉在线视频| 亚洲第一中文字幕| 5388国产亚洲欧美在线观看| 超碰91免费人妻| 国产精品福利尤物youwu| 亚洲区一区| 中文无码精品A∨在线观看不卡| 免费看久久精品99| 欧美国产在线精品17p| 动漫精品中文字幕无码| 国产亚洲精品自在线| 丁香综合在线| 美美女高清毛片视频免费观看| 成人午夜网址| 97人人模人人爽人人喊小说|