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

中國50百分位男性小腿有限元模型的建立與驗證*

2015-04-12 06:27:12曹立波杜現平張冠軍胡躍群
汽車工程 2015年11期
關鍵詞:有限元實驗模型

曹立波,杜現平,張冠軍,胡躍群,張 愷

(1.湖南大學,汽車車身先進設計制造國家重點實驗室,長沙 410082; 2.中南大學湘雅三醫院放射科,長沙 410013)

?

2015219

中國50百分位男性小腿有限元模型的建立與驗證*

曹立波1,杜現平1,張冠軍1,胡躍群2,張 愷1

(1.湖南大學,汽車車身先進設計制造國家重點實驗室,長沙 410082; 2.中南大學湘雅三醫院放射科,長沙 410013)

基于下肢螺旋CT掃描,保留主要的解剖學結構,建立了中國50百分位男性小腿有限元模型。應用中國人體測量學數據進行模型縮放,對皮質骨厚度進行精確的模擬。基于Ls-Dyna材料模型模擬皮質骨拉、壓特性和應變率特性。對模型進行小腿的準靜態軸向壓縮驗證以及脛骨、腓骨與小腿的準靜態三點彎曲驗證和近心端1/3、中部和遠心端1/3的動態三點彎曲驗證。驗證結果顯示:仿真驗證曲線與試驗曲線走勢吻合較好,峰值力大小與出現時刻基本一致,表明所建立的模型生物逼真度較高。

小腿損傷;有限元模型;驗證

前言

美國國家汽車采樣系統(national automotive sampling system, NASS)統計顯示,約1/5的乘員損傷[1]和4/5的行人損傷[2]包含下肢損傷,下肢AIS 2+損傷在身體各部位損傷中位居第二[1]。下肢損傷不僅會給受害者造成巨大的傷痛,還會造成生活不便和很大的社會損失。建立下肢有限元模型,研究下肢損傷機理對優化車輛設計參數,減少下肢損傷具有重要意義。

在下肢損傷的研究方法中,有限元分析因其模型開發周期短,成本低,可以模擬真實的應力應變和骨折等優點,逐漸成為主要的研究方法之一[3-4]。

下肢有限元模型的發展經歷了一個較長的歷程。文獻[5]中將下肢骨定義為剛體,建立了乘員下肢有限元模型。文獻[6]中用可變形的殼單元模擬長骨皮質骨,肌肉和皮膚采用假人有限元模型相應材料參數,建立了行人下肢有限元模型。文獻[7]中用單層實體單元模擬皮質骨,采用非線性、彈塑性和應變率材料模擬長骨材料,建立了50百分位行人下肢有限元模型。文獻[8]中用實體單元模擬骨干皮質骨,用五面體單元實現骨干體單元到骨骺殼單元皮質骨的過渡,建立了男性50百分位行人下肢有限元模型。文獻[9]中用較新型的材料模擬骨骼,建立了行人下肢有限元模型。

這些模型具有一定應用價值,但皮質骨幾何和材料特性模擬存在不足,且多基于歐美人體建立。因此,須要建立中國50百分位男性下肢有限元精確模型[10],用于適合中國人體的汽車參數的優化和相關安全法規的修訂研究。

1 小腿有限元模型的建立

通過一個身高173.1cm、體質量69.7kg的成年男性患者(無下肢病患)螺旋CT掃描(層距1mm),提取幾何數據,對照解剖學圖譜進行光順,保留主要解剖學特征,建立了小腿幾何模型。

1.1 小腿網格模型

采用六面體網格劃分,并對網格質量進行控制。脛骨干皮質骨采用4層實體單元,骨骺皮質骨采用殼單元模擬;為模擬皮質骨厚度實體單元到殼單元的漸變,將骨干實體皮質骨階梯過渡到骨骺端殼單元皮質骨,如圖1所示。

基于CT測量,將脛骨兩骨骺端皮質骨沿軸向分別劃分為6層,如圖2所示;在每層相應的CT圖像中選取6個測量點進行測量并取均值,將均值作為對應殼單元的厚度,如表1所示。腓骨干皮質骨用兩層實體單元模擬,骨骺皮質骨用1.5mm厚的殼單元模擬[11]。肌肉和皮膚分別用實體單元和1mm厚殼單元模擬[12],并采用共節點連接;肌肉和骨骼之間采用節點綁定連接。

表1 骨骺端皮質骨厚度CT測量結果

由于缺少較為全面準確的中國人體統計數據,且GB 10000—88是人體外部尺寸標準[13],不能作為長骨的縮放依據。小腿長度與脛骨長度測量基準相近,故采用脛骨全長與GB 10000—88中小腿長度相近的文獻[14]為基準,如表2所示。按比例進行模型整體縮放,獲得了中國50百分位男性小腿模型。

表2 脛骨參數縮放表

建立的小腿模型包含:44 640個節點,36 115個實體單元和2 634個殼單元,如圖3所示。

1.2 小腿材料模型

基于LS-Dyna建立小腿有限元模型。在交通事故中,長骨主要受拉、壓載荷。由于皮質骨具有拉、壓不同的材料特性,且材料性能受到應變率的影響顯著[8,15]。因此,采用124號彈塑性材料模擬皮質骨的拉、壓特性,并采用Cowper-Symonds模型定義應變率的影響,如圖4所示。選用105號黏彈塑性材料模擬松質骨的黏彈性和其塑性[7-8],松質骨和皮質骨材料均采用單元失效刪除算法模擬骨折。

依據已有模型[8,16-17]的材料參數范圍,選取合適的參數值,骨骼的主要材料參數如表3所示。

表3 骨骼材料參數

肌肉和皮膚在壓縮時,容易因大變形而產生負體積和沙漏能,為降低其出現的幾率,采用單點積分和6號沙漏控制模式進行負體積和沙漏能的控制,并分別選用黏彈性和彈性材料模擬肌肉和皮膚材料。參數依據文獻[11]和文獻[12]中的模型選取,如表4所示。

表4 小腿軟組織模型材料參數

2 基于交通損傷的小腿模型驗證

針對模型應用工況進行合理的仿真驗證,是模型具有較好的生物逼真度和精度的基礎。基于交通損傷研究的小腿有限元模型,須要針對交通事故中小腿的受力和損傷特點進行驗證。

2.1 小腿交通損傷的力學特性

在交通事故中,小腿損傷按照對象不同,主要有行人小腿損傷和乘員小腿損傷。由于載荷形式不同,其損傷機理和防護方法等也不同。

在交通事故中行人小腿側面首先與汽車保險杠接觸,由于上半身慣性和腳部與地面的摩擦作用,在保險杠的沖擊作用下,小腿主要受到三點彎曲和沖擊擠壓的作用。長骨是小腿的主要承力組織,當應變超過極限值時則發生骨折。骨折的發生,對膝關節和踝關節損傷的影響較大。一般認為,在不是特別嚴重的交通事故中,行人膝關節韌帶損傷和骨折只發生一種的幾率較大。小腿受到沖擊的位置,因不同車型的保險杠高度不同,會產生較大的變化[3];因此,有必要對小腿骨和小腿不同位置的動力學彎曲特性進行驗證分析。

在正面碰撞的交通事故中,乘員的小腿在儀表板、前輪、地板和腳踏板入侵的作用下,主要受到軸向力、背屈、內翻和外翻4種載荷工況。背屈和內、外翻主要導致踝關節相關腳骨和軟組織的損傷。軸向力則主要導致小腿中下部的骨折,在沖擊能量較高時,甚至出現粉碎性骨折和踝關節骨的骨折,即pilon骨折[18]。因此,小腿模型遠端在軸向載荷作用下的響應,對于模型準確模擬pilon骨折具有重要影響。

2.2 小腿模型的驗證

基于交通損傷中行人和乘員小腿的載荷特點,依據文獻[19]~文獻[25]中的實驗,對脛骨、腓骨和小腿模型進行了準靜態(QS)和動態(DY)驗證,如表5所示。

表5 小腿有限元模型仿真驗證

注:L-M:側向;A-P:前后;Pro1/3:近心端1/3;Mid:中部;Dist1/3:遠心端1/3。

2.2.1 小腿模型準靜態驗證

脛骨和腓骨是小腿模型的主要部分,準靜態驗證是驗證材料模型有效性的基礎。模型主要采用文獻[19]~文獻[21]、文獻[24]~文獻[25]中的實驗,對脛骨、腓骨和小腿進行準靜態驗證。

依據文獻[19]關于長骨準靜態實驗的描述,建立長骨驗證有限元模型。將長骨兩端支撐在剛性平面,設置支承平面與長骨兩端的接觸,約束剛性平面的6個自由度。用直徑為25mm的剛性圓筒沖擊器的圓柱面,以0.01m/s的速度對長骨中部進行加載,沖擊器與長骨設置接觸。通過設置接觸面之間的摩擦和長骨兩端相應節點的約束,來限制驗證時長骨在剛性平面上的異常運動。

脛骨中部截面為三角形,A-P方向三點彎曲時,其頂角與沖擊器接觸,為防止剛性接觸而產生較大的應力集中,造成不合實際的單元失效,對剛性接觸處個別單元不應用單元失效刪除準則[1,26]。

依據實驗獲得的結果,輸出模型驗證中的沖擊力-位移曲線,與實驗曲線進行對比。脛骨和腓骨準靜態驗證的仿真設置,如圖5所示。

針對汽車正面碰撞的特點,依據文獻[24]和文獻[25]中的實驗,對小腿進行準靜態軸向壓縮驗證,確定小腿的軸向受壓極限和材料準確性。將脛骨近心端與剛性的方盒剛性連接,約束方盒6個方向自由度。小腿遠心端與剛性方盒內的聚氨酯表面之間初始間距為2~3mm,并定義接觸。限制遠心端方盒的自由度,只允許其軸向移動和矢狀面轉動,沖擊器以2mm/s的速度對遠心端方盒進行軸向加載,直至骨折發生,如圖6所示。

2.2.2 小腿模型動態驗證

在道路交通事故中,下肢損傷多是由碰撞引起的動載荷導致。因此,在準靜態驗證的基礎上,根據交通事故中小腿受到動載荷的特點,對脛骨、腓骨和小腿模型的近心端、中部和遠心端分別進行了動態三點彎曲驗證。早期實驗中實驗設置、邊界條件和樣本尺寸缺乏較為完善的描述,不利于仿真驗證的設置。因此,采用描述較完善的文獻[22]中的實驗進行動態驗證。

依據該實驗,對脛骨、腓骨和小腿有限元模型進行驗證設置。將模型近心端和遠心端與兩端剛性盒固接,剛性盒與下端剛性弧面固接。剛性圓弧面可以在剛性平面上自由轉動,約束剛性平面6個自由度。由于聚氨酯密度比剛性盒小很多,故可忽略其質量。剛性沖擊器前端包裹25mm厚的泡沫材料,以1.45m/s的速度對模型進行加載,并在弧面與剛性平面、泡沫與長骨之間定義接觸。

圖7為同一模型不同部位的仿真驗證設置。

3 仿真驗證結果分析

通過對實驗結果與仿真驗證結果的曲線走勢、峰值大小和出現時刻進行對比,評價模型的有效性和生物逼真度。

3.1 準靜態仿真驗證結果對比

模型主要通過輸出仿真的沖擊力-位移曲線與實驗曲線進行對比。主要依據文獻[19]~文獻[21]、文獻[24]和文獻[25]中的實驗,對脛骨和腓骨進行準靜態A-P和L-M方向的仿真驗證及小腿軸向壓縮驗證,并將仿真與實驗結果進行對比分析,如圖8~圖10所示。

由圖8可知:脛骨在A-P方向準靜態驗證沖擊力-位移曲線與實驗結果基本一致,骨折力和位移與實驗吻合較好;L-M方向曲線略有差別,峰值力大小與出現的時刻一致,仿真曲線總體較實驗曲線略小,這與實驗樣本尺寸有關,但仿真仍能夠和實驗大體吻合。說明所建立的脛骨有限元模型能夠較準確地模擬實驗中的脛骨響應。

由圖9可知:腓骨A-P方向仿真驗證結果與實驗曲線吻合較好,骨折最大位移與最大力基本一致;但是,腓骨L-M方向準靜態驗證曲線較實驗要小,這與兩實驗曲線樣本不同有關。而且,腓骨中部A-P方向尺寸較L-M方向大,因此,同一腓骨的三點彎曲,在A-P方向比L-M方向能夠承受更大的力,這與模型仿真結果吻合,而文獻[19]中的實驗在兩個方向沒有顯示出峰值力的明顯差別,這可能是由于兩曲線非同一樣本實驗所致,但曲線走勢仍能與實驗吻合較好,說明所建立的腓骨有限元模型仍具有較好的仿真精度。

由圖10可見:所建立的模型在軸向壓縮時,沖擊力-位移曲線走勢與實驗曲線一致,骨折力為10.2kN,在通道峰值13.6~7.8kN范圍內,相對通道均值偏小,較好反映了中國人體比西方人體平均尺寸小的特點,骨折出現的時刻也較吻合,說明所建立的模型能夠較好地模擬小腿在受到軸向載荷時的響應,仿真精度較好。

3.2 動態仿真驗證結果對比

基于行人載荷特點,對下肢脛骨、腓骨和小腿模型進行了近端1/3、中部、遠端1/3的動態沖擊仿真驗證,依據文獻[22]和文獻[23]中的實驗,對模型加載至骨折,輸出仿真中的力-位移曲線和骨折力矩,與實驗曲線和耐受限度值進行對比。

為補償由于樣本尺寸不同而造成的實驗結果與50百分位中國人體下肢仿真結果的差異,采用力矩、力和位移縮放系數(分別為λM,λF,λD)[7]對尸體實驗的結果進行縮放后,與仿真驗證結果進行對比,所采用的縮放系數計算公式如表6所示。

表6 縮放系數計算公式

表6中λL為尺寸縮放系數,實驗結果除受到尺寸影響外,也會受實驗樣本質量的影響,因此,引入質量縮放系數λma,得到等效縮放系數λLe。

將脛骨和腓骨實驗結果,依據上述方法,縮放至中國50百分位男性(脛骨[14]:353.8mm;體質量[13]:59kg)。至于小腿動態驗證結果,由于進行對比的是由大量實驗統計獲得的響應通道范圍,故未進行縮放。脛骨、腓骨和小腿仿真驗證曲線與實驗曲線對比如圖11~圖13所示。

仿真峰值和縮放后的耐受限度值[3]見表7。

由圖11可知,脛骨沖擊力-位移驗證曲線的走勢與縮放后的實驗曲線[22]吻合較好,且峰值力大小與出現的時刻基本一致,峰值力均在耐受限度范圍內,說明所建立的脛骨模型能夠較好地模擬人體脛骨的動態響應。

由圖12可知,腓骨仿真驗證曲線與縮放后的實驗曲線[22]基本一致,最大峰值力大小與出現的時刻均在實驗曲線范圍之內。由于遠心端動態三點彎曲實驗樣本較少,仿真驗證峰值偏大,但是仿真曲線的走勢與實驗曲線基本吻合,且峰值力0.5kN與耐受限度值0.43kN接近,骨折時刻基本一致,說明所建立的腓骨模型能夠模擬真實實驗,具有較好的精度。

由圖13可知,小腿動態三點彎曲沖擊力-位移曲線在通道范圍[23]內變化,且曲線走勢基本一致,最大峰值力也與實驗通道最大峰值力基本吻合。遠端三點彎曲仿真曲線開始階段走勢與通道邊界曲線不同,可能是由于CT掃描時,小腿遠端支撐在CT臺致使肌肉變形,未提取遠端變形較嚴重的肌肉組織,導致遠心端沖擊點位置相對實驗設置偏向中部,肌肉組織較厚,沖擊器不能很快接觸骨骼,仿真曲線在開始階段上升較慢,但仿真曲線仍在實驗曲線通道范圍內。近心端與中部仿真結果在縮放后的通道內變化,且走勢基本相同,骨折力基本吻合,說明模型能夠模擬真實的小腿動態沖擊響應,可以用于下肢模型的構建。

表7 動態仿真結果

4 結論

通過近似50百分位中國男性下肢螺旋CT掃描,提取下肢幾何模型,劃分網格,控制網格質量,獲得高質量的網格模型。基于中國50百分位人體統計學數據進行縮放,采用更精細的解剖學有限元描述,有效提高模型的幾何精度。模型的皮質骨采用拉、壓不同特性和應變率特性的材料模型模擬,松質骨采用彈黏塑性材料模擬,提高了模型的材料精度。為降低肌肉和皮膚出現負體積的幾率,采用黏彈性和彈性材料模擬肌肉和皮膚,長骨與肌肉之間采用節點固聯接觸,并設置長骨之間的接觸,采用單元失效刪除算法模擬骨折。基于行人和乘員的損傷形式和載荷特點,對脛骨、腓骨和小腿模型進行了準靜態和近心端1/3、中部、遠心端1/3的動態三點彎曲驗證以及小腿軸向準靜態壓縮驗證。仿真結果與實驗曲線吻合較好,峰值力與彎矩出現的時刻基本一致,表明建立的小腿有限元模型能夠較好地模擬人體小腿在實驗中的響應,可用于后續中國男性50百分位人體下肢有限元模型的開發,下肢交通損傷機理的研究和汽車安全參數的優化,對符合國情的汽車安全法規的修訂具有重要的參考價值。

[1] Kim Y S, Choi H H, Cho Y N, et al. Numerical Investigations of Interactions Between the Knee-thigh-hip Complex with Vehicle Interior Structures[J]. Stapp Car Crash Journal,2005,49:85-115.

[2] Saukko P, Knight B. Transportation Injuries[J]. Knight’s Forensic Pathology, 3rd ed. London: Arnold,2004:281-300.

[3] 張冠軍.行人下肢的碰撞損傷特性及相關參數研究[D].長沙:湖南大學,汽車車身先進設計制造國家重點實驗室,2009.

[4] 杜現平,曹立波,張冠軍,等.人體下肢三維有限元模型在交通損傷中的應用研究進展[J].汽車工程學報,2014(4):235-244.

[5] Beaugonin M, Haug E, Cesari D. A Numerical Model of the Human Ankle/Foot under Impact Loading in Inversion and Eversion[C]. SAE Paper 962428.

[6] Schuster P J, Chou C C, Prasad P, et al. Development and Validation of a Pedestrian Lower Limb Non-linear 3-d Finite Element Model[J]. Stapp Car Crash Journal,2000,44:315-334.

[7] Takahashi Y, Kikuchi Y, Mori F, et al. Advanced FE Lower Limb Model for Pedestrians[C]. 18th International ESV Conference, Paper 2003.

[8] Untaroiu C, Darvish K, Crandall J, et al. A Finite Element Model of the Lower Limb for Simulating Pedestrian Impacts[J]. Stapp Car Crash J,2005,49:157-181.

[9] Li Z D, Zou D H, Liu N G, et al. Finite Element Analysis of Pedestrian Lower Limb Fractures by Direct Force: the Result of Being Run Over or Impact?[J]. Forensic Sci Int,2013,229(1-3):43-51.

[10] 杜現平.中國人體50百分位小腿有限元模型開發及骨折特性研究[D].長沙:湖南大學,2015.

[11] Iwamoto M, Tamura A, Furusu K, et al. Development of a Finite Element Model of the Human Lower Extremity for Analyses of Automotive Crash Injuries[C]. SAE Paper 2000-01-0621.

[12] Snedeker J G, Muser M H, Walz F H. Assessment of Pelvis and Upper Leg Injury Risk in Car-pedestrian Collisions: Comparison of Accident Statistics, Impactor Tests and a Human Body Finite Element Model[J]. Stapp Car Crash Journal,2003,47:437-457.

[13] 國家技術監督局.中國成年人人體尺寸GB 10000—88[S].北京:中國標準出版社,1989.

[14] 劉俊先,張興和.中國正常人體測量值[M].北京:中國醫藥科技出版社,1994.

[15] Novitskaya E, Chen P Y, Hamed E, et al. Recent Advances on the Measurement and Calculation of the Elastic Moduli of Cortical and Trabecular Bone: A Review[J]. Theoretical and Applied Mechanics,2011,38(3):209-297.

[16] Takahashi Y, Kikuchi Y, Konosu A, et al. Development and Validation of the Finite Element Model for the Human Lower Limb of Pedestrians[J]. Stapp Car Crash Journal,2000,44:335-355.

[17] Beillas P, Begeman P C, Yang K H, et al. Lower Limb: Advanced FE Model and New Experimental Data[J]. Stapp Car Crash J,2001,45:469-494.

[18] Bedewi P. The Biomechanics of Human Lower Extremity Injury in the Automotive Crash Environment[D]. The School of Engineering and Applied Science, Washington DC: George Washington University,1998:286.

[19] Yamada H, Evans F G. Strength of Biological Materials[M]. Baltimore: The Williams & Wilkins Company,1970.

[20] Asang E. Experimental Biomechanics of the Human Leg:. A Basis for Interpreting Typical Skiing Injury Mechanisms[J]. Orthop Clin North Am,1976,7(1):63-73.

[21] Ehler E, Ouml, Sche H. Human Tibia Under Bending Load[J]. Beitr Orthop Traumatol,1970,17(5):291-304.

[22] Kerrigan J R, Bhalla K S, Funk J R, et al. Experiments for Establishing Pedestrian-impact Lower Limb Injury Criteria[C]. SAE Paper 2003-01-0895.

[23] Ivarsson B J, Kerrigan J R, Lessley D J, et al. Dynamic Response Corridors of the Human Thigh and Leg in Non-midpoint Three-point Bending[C]. SAE Paper 2005-01-0305.

[24] Begeman P, Paravasthu N. Static and Dynamic Compression Loading of the Lower Leg[C]. Proceedings of the Seventh Injury Prevention Through Biomechanics Symposium,1997.

[25] Begeman P, Paravasthu N. Static and Dynamic Compression and Torsion Loading of the Lower Leg[R]. Final Report Submitted to the AAMA,1997.

[26] Untaroiu C D, Yue N, Shin J. A Finite Element Model of the Lower Limb for Simulating Automotive Impacts[J]. Ann Biomed Eng,2013,41(3):513-526.

Development and Validation of the Lower Leg FEModel for the 50th Percentile Chinese Male

Cao Libo1, Du Xianping1, Zhang Guanjun1, Hu Yuequn2& Zhang Kai1

1.HunanUniversity,StateKeyLaboratoryofAdvancedDesignandManufacturingforVehicleBody,Changsha410082;2.DepartmentofRadiology,TheThirdXiangyaHospitalofCentralSouthUniversity,Changsha410013

Based on the spiral CT scan of lower limbs with major anatomical structures reserved, a lower leg FE model for the Chinese 50th percentile male is established. Model scaling is conducted based on the anthropometry for Chinese with cortical bone thickness accurately simulated, and the tension, compression and strain characteristics of cortical bones are simulated by the material model in LS-DYNA. Then the quasi-static axial compression verification for lower leg and the quasi-static three-point-bending verifications (including proximal 1/3, middle and distal 1/3) for tibia, fibula and lower leg are performed on the model. The results indicate that the trends of simulation verification curves agree well with that of test curves, with the magnitudes and appearing moments of peak forces basically coincided, demonstrating the high biofidelity of the model built.

lower leg injury; FE model; validation

*國家自然科學基金(51205118)和湖南大學汽車車身先進設計制造國家重點實驗室自主研究課題(51275001)資助。

原稿收到日期為2013年1月13日,修改稿收到日期為2014年5月27日。

猜你喜歡
有限元實驗模型
一半模型
記一次有趣的實驗
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
做個怪怪長實驗
3D打印中的模型分割與打包
NO與NO2相互轉化實驗的改進
實踐十號上的19項實驗
太空探索(2016年5期)2016-07-12 15:17:55
磨削淬硬殘余應力的有限元分析
基于SolidWorks的吸嘴支撐臂有限元分析
主站蜘蛛池模板: www.99精品视频在线播放| 99视频在线免费| 日韩av在线直播| 午夜毛片免费看| 熟妇人妻无乱码中文字幕真矢织江| 国产剧情国内精品原创| 手机精品福利在线观看| AV天堂资源福利在线观看| 国产鲁鲁视频在线观看| 全裸无码专区| 波多野结衣二区| 久久美女精品国产精品亚洲| 国产一区二区免费播放| 全部免费毛片免费播放| 一级不卡毛片| 97视频免费在线观看| 国产精品爆乳99久久| 亚洲人成网站在线观看播放不卡| 国产激情无码一区二区APP| 婷婷丁香色| 亚洲高清中文字幕在线看不卡| 日本一区高清| 91网在线| 国产永久在线视频| 国产成人精品一区二区| 国产精品毛片在线直播完整版| 亚洲欧洲一区二区三区| 亚洲天堂日韩在线| 蜜臀av性久久久久蜜臀aⅴ麻豆| 特黄日韩免费一区二区三区| 国内精品91| 在线观看免费黄色网址| 一级福利视频| 92精品国产自产在线观看| 试看120秒男女啪啪免费| 国产97视频在线| 91青青在线视频| 国产在线视频福利资源站| 激情无码字幕综合| 91亚洲精选| 无码精品一区二区久久久| 亚洲成a人片| 99视频全部免费| 国产成人8x视频一区二区| 免费毛片视频| 婷婷色婷婷| 亚洲欧美日韩综合二区三区| 激情乱人伦| 一级做a爰片久久毛片毛片| 91福利片| 综合久久五月天| 成人精品午夜福利在线播放| 亚洲国产精品日韩av专区| 欧洲欧美人成免费全部视频| 国产内射在线观看| 伊在人亚洲香蕉精品播放| 一级毛片a女人刺激视频免费| 欧美日韩成人在线观看| 国产欧美日韩综合一区在线播放| 香蕉精品在线| 久久精品人妻中文视频| 亚洲天堂免费在线视频| 日韩AV无码免费一二三区| 国产精品偷伦在线观看| 2021国产精品自拍| 亚洲成肉网| 91视频国产高清| 毛片免费在线视频| 国产人成乱码视频免费观看| 国产高颜值露脸在线观看| 亚洲日本中文字幕乱码中文 | 二级特黄绝大片免费视频大片| 欧美色综合网站| 啪啪国产视频| 亚洲成人77777| 欧洲高清无码在线| 欧美有码在线观看| 在线视频一区二区三区不卡| 伊人久久久大香线蕉综合直播| 久久公开视频| 有专无码视频| 亚洲人成日本在线观看|