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

電動汽車車身的回傳射線矩陣剛度鏈分析方法

2016-12-23 02:58:13陳毅強劉子建
中國機械工程 2016年23期
關鍵詞:優化結構設計

陳毅強 劉子建

湖南大學汽車車身先進設計制造國家重點實驗室,長沙,410082

?

電動汽車車身的回傳射線矩陣剛度鏈分析方法

陳毅強 劉子建

湖南大學汽車車身先進設計制造國家重點實驗室,長沙,410082

基于LifeDrive純電動汽車模塊化結構形式,通過拓撲優化得到一款純電動汽車車架結構;建立了基于梁單元的車架結構回傳射線矩陣計算模型;由模型得到車架梁單元尺寸以及車身與車架之間的多個耦合力,以更合理的剛度鏈方法優化了車身簡化模型的主斷面參數;利用有限元方法計算了車身與車架的承載度。結果證明設計的車架與車身結構滿足目標承載度和剛度性能的要求,該方法為非承載式電動汽車車身及車架的概念設計提供了參考。

純電動汽車;車架設計;承載度;回傳射線矩陣法

0 引言

能源危機和環境保護要求汽車產品向微型化、低排放方向發展,作為解決能源和環保問題方案之一的電動汽車越來越受到人們的關注。純電動汽車車身的正向設計技術已成為各國汽車研發的一個重點。

新型材料的應用、客戶需求多樣化正在影響純電動汽車車身結構朝著模塊化、平臺化方向發展。典型的如寶馬i系列車采用了LifeDrive架構,它由Life和Drive兩個獨立的模塊組成,Life模塊是乘員艙部分,采用超輕量化且高強度CFRP碳纖維復合材料構成,Drive模塊包括懸架、蓄電池組、驅動系統和碰撞防護結構等。

這種模塊化架構車輛采取非承載式車身時,通常認為載荷均由車架承擔,并沒有考慮車身的承載作用。此外,如何合理設計純電動汽車的車架,目前還沒有成熟的理論和方法[1]。

國內外學者研究車身及車架設計優化方法取得了多項進展。如田海豹[2]提出了用于車身結構概念設計階段的剛度鏈方法,通過建立以主斷面、接頭等為節點的車身剛度鏈模型,進行車身載體剛度優化和輕量化設計,并給出了承載式車身設計的實例。文獻[2]使用傳遞矩陣法計算車身剛度,用于梁結構分支較多的非承載式電動車車架分析時會產生子剛度鏈分解復雜、計算效率低的問題。扶原放等[3]依據微型電動汽車車架結構的受力特性及材料性能要求,考慮多種行駛沖擊載荷對車架的作用,建立了設計優化數學模型。高云凱等[4]基于拓撲優化方法研究了某非承載式電動汽車車身。周姍姍[5]利用客車有限元模型,基于車身承載度分析對車身結構進行了改進研究。Cavazzuti等[6]基于有限元方法研究了以車架結構性能為約束、以質量為優化目標的汽車底盤框架結構的設計方法。Hodkinson等[7]研究了電動汽車的輕量化設計流程、車身結構分析方法以及有限元仿真方法。梁晨等[8]以某非承載式車身及車架為例,利用有限元方法研究了彎曲和扭轉工況下的車身及車架對整車的剛度貢獻率。

如何在非承載式車身的概念設計階段就充分考慮車架和車身對于載荷的貢獻,形成合理的設計方法,并獲得車身整體輕量化效果,這一問題仍然沒有較好的解決方法。因此,研究非承載式車身與車架承載度優化分配的方法,以及先進的車架設計方法,對于提升純電動汽車整車輕量化水平和綜合性能都具有重要意義。

1 車架結構設計

1.1 承載度分配與剛度設計目標值的確定

本文研究的純電動汽車結構形式如圖1所示,車架承載度設計流程如圖2所示。

圖1 純電動汽車結構形式

圖2 承載度設計流程

在總體結構方案基本確定的情況下,選擇不同的車身承載度目標值會得到不同的車身與車架的質量比。本文研究A00級電動汽車,參考文獻[8]中的非承載式SUV車身對整車剛度的貢獻,設定車身與車架承載度的目標比值為3∶7,車身與車架通過接口連接,前座椅中心位置下方為撓度測量點,根據文獻[9]中電動汽車設計過程中的性能指標確定方法,此處最大撓度設定為0.497mm。考慮彎曲工況下的載荷為多點集中載荷,分析過程中可將車身、車架都按照簡支梁計算。先繪制車架簡化受力圖(圖3),此時可得車架的彎曲剛度計算公式[9]:

圖3 車架受多點集中載荷受力圖

(1)

式中,EI為彎曲剛度值,N·m2;yx為測定點撓度值,m;x為前軸中心位置到測定點的距離,m;l為軸距,m;Pi為各點加載載荷,N;ai為前軸中心到加載點的距離,m;bi=l-ai,m;j為前軸前載荷項數;k為至撓度測定點的載荷項數;n為后軸前載荷項數;m為總載荷項數。

由式(1)計算出EI=2.76 MN·m2。車架承載度目標值為0.7時,測量點最大撓度約束值為1.43yx=0.7115mm。

根據普通車型扭轉剛度結果統計數據[10],此A00級車扭轉角目標值θ不應超過0.5°,前懸架支座處位移可表示為

ΔL=(b/2)tanθ

(2)

其中,ΔL為前懸支座位移;b為前懸支座間距。本文中b=610mm,由式(2)可計算出前懸支座處最大位移約束為2.71mm。

1.2 拓撲優化設計

車架所受靜載荷F1包括乘客、電池組、電機的載荷,大小約為4500N,如表1所示。

表1 主要載荷部件

部件名稱數量質量(kg)位置電池包165車架中部電機驅動系統140前軸中部電機控制器15前輪內側乘員及座椅4340車架中部

載荷均施加于設計域的上表面,并將電池對車架的載荷考慮為多點集中載荷[11],加載方向均垂直于地面。在靜載荷基礎上對模型施加彎矩載荷F=1.8F1與扭轉載荷M=0.5F2l,其中F2為前軸載荷。根據剛度和頻率要求,以結構整體的體積約束作為優化的目標函數,以結構的剛度最大化作為優化的約束條件,以前座椅測量點處撓度及前懸支座處最大位移為邊界條件[12],建立拓撲優化的設計域,如圖4所示。

圖4 車架拓撲優化設計域

選擇材料為45鋼,彈性模量為200GPa,泊松比為0.3,屈服強度為355MPa。剛度約束條件下以體積為目標函數,剛度要求為約束條件(文中反映為測量點撓度及最大扭轉角度),單元密度為設計變量;頻率約束條件下以一階模態頻率為目標函數,體積分數為約束條件,單元密度為設計變量,設定收斂容差為0.001,若兩次迭代之差小于0.001,則認為優化收斂并停止計算[12]。經迭代計算后單元密度取0.17,優化結果如圖5所示。

圖5 車架結構材料分布

考慮各總成的布置以及工藝要求,采用截面為矩形的型鋼制作車架,將拓撲結構簡化成圖6所示的形式,但此時只確定了車架的結構形式,還沒有確定矩形鋼管的厚度,因此還需要進一步的優化設計。

圖6 車架結構圖

2 基于剛度鏈的車架建模與分析

圖6所示車架模型的分支梁較多,當一個節點與多個單元相連時,采用基于傳遞矩陣的剛度鏈方法建模時不易確定傳遞路徑和節點的傳遞矩陣[13],且計算量大,而利用回傳射線矩陣法(methodofreverberaton-raymatrix,MRRM)則可以避免這些缺陷。MRRM方法對所有單元列相位關系的列式方式與對所有節點列散射關系的方式相同,總體相位矩陣和散射矩陣分別是將各節點的局部相位矩陣和散射矩陣放置在矩陣的對角線上得到的分塊對角矩陣,列式非常統一。MRRM基本思路是:依據兩組關系建立整體結構的回傳矩陣,第一組關系為節點的近端位移和遠端位移的關系,由節點的力和力矩的平衡關系建立;第二組關系為桿件的近端位移與遠端位移的關系,由局部坐標系的設定得到。由回傳矩陣即可求得整體結構各節點的位移和內力的精確解[14]。

2.1 結構描述

將梁的中心線相連,得到車架結構的簡化線框模型,節點以大寫字母I、J、K…表示,則桿件表示為〈I,J〉、〈J,K〉等。圖7中,各桿件的編號表示為〈1,2〉,〈2,3〉,…,〈27,28〉。

圖7 車架簡化線框模型

車架結構節點總數為28,梁單元總數為48,f1~f5為車身與車架之間的耦合力,利用對稱關系,先考慮1/2模型的分析。

2.2 建模方法

2.2.1 定義坐標系

首先,對車架結構建立整體坐標系(X,Y,Z),如圖7所示。以與節點1相連的三根梁為例,對梁單元建立對偶坐標系,如圖8所示,各梁單元的長度記為l〈1,2〉,l〈2,3〉,…,l〈27,28〉。

圖8 梁單元局部坐標系

2.2.2 單元物理量的表示方法

在整體坐標系(X,Y,Z)下,定義節點J(J=1,2,…,28)的廣義力和廣義位移向量分別為

(3)

其中,pJ中包含沿X軸、Y軸、Z軸的集中力和繞三個坐標軸的力矩,uJ中包含沿X軸、Y軸、Z軸三個方向的位移和轉角。

在局部坐標系(x,y,z)〈J,K〉下,梁單元任意截面的廣義力和廣義位移向量分別為

(4)

其中,f〈J,K〉中包含z方向的軸力、x和y方向上的剪力、繞z軸的扭矩、繞x和y軸的彎矩;δ〈J,K〉中包含三個沿x、y和z軸的線位移和三個繞x、y、z的轉角。

2.2.3 力平衡關系與位移協調關系

車架結構的梁與梁之間的連接點簡化為無集中質量的剛性節點,在節點J處有平衡關系:

(5)

式中,T〈J,K〉為66階轉換矩陣;F〈J,K〉為節點J處61階內力向量,包含軸剪力和扭矩彎矩和;F為外力向量;fi(i=1,2,3,4,5)為車架與車身的耦合力向量。

(6)

A〈J,K〉和D〈J,K〉由空間桿系結構的剛度矩陣得到:

A〈J,K〉=

(7)

D〈J,K〉=

(8)

式中,EA、GJ分別為車架梁單元的抗壓剛度和剪切剛度。

根據式(5)和式(6)可得第一組關系,即包含節點近端位移、遠端位移的力和力矩平衡關系:

(9)

(10)

式中,mJ為與節點J相連的桿件數量。

在節點3與節點17處,F≠0,fi≠0;在節點1、4、5、7處,F=0,fi≠0;其他節點處,F=0,fi=0。

除了上述力和力矩平衡關系外,車架的梁結構在連接節點處還滿足位移協調條件:

(11)

(12)

2.2.4 車架結構的回傳矩陣

通過力平衡方程式(9)和位移協調方程式(11)可得到節點的波源向量s和傳遞分配矩陣S。則有如下關系:

a=Sd+s

(13)

其中,a和d是所有節點的近端位移向量aJ與遠端位移向量dJ組成的向量:

(14)

式(13)中有12mJ個方程,24mJ+5個未知數,分別為節點近端和遠端沿三個坐標軸方向上的位移、轉角、力和力矩,如果考慮5個未知耦合力,則還需要建立另外一組關系。

在局部坐標系中,桿件的遠端位移可以用近端位移表示:

d=PUa

(15)

對于某一桿件〈J,K〉,式(15)中的d包含節點J的遠端位移d〈J,K〉實際上是節點K的近端位移a〈K,J〉,但兩者可能存在正負號的差異。利用相位矩陣p可以將兩者關系表示為

d〈J,K〉=pa〈J,K〉

(16)

J=1,2,…,mJK=1,2,…,n

式(15)中的U為轉列矩陣,用于改變向量a中元素的排列次序,使之與向量d中的元素相對應,例如d〈J,K〉對應a〈K,J〉;P為整體相位矩陣,表示為

P=diag(p,p,…,p)

由式(12)和式(15)可得:

a=(I-R)-1s

(17)

R=SPU

式中,I為6×6階單位矩陣。

由于耦合力未知,為了避免未知數個數多于方程個數,需要確定彎曲工況下部分節點的位移。根據已確定的車架最大撓度以及邊界條件進行三次擬合可以近似得到車架彎曲工況下的變形曲線,如圖9所示。

圖9 車架結構變形擬合曲線

由圖9可以得到yi(i=1,3,4,5,7),即節點1、節點3、節點4、節點5、節點7在整體坐標系中沿Z軸負方向的位移:

(18)

由于與節點1、3、4、5和7連接的單元的近端或遠端狀態向量a〈J,K〉中的元素v〈J,K〉可以由擬合曲線的坐標得到,由式(13)和式(15)可以得到的12mJ個方程中,含12mJ個未知數,其中包含5個已知的節點位移v1、v3、v4、v5、v7,其大小對應式(18)中的y值和5個未知的耦合力。車架梁截面采用30 mm×30 mm的矩形,厚度設為t,因此由不受外力作用的節點的關系可解出梁單元厚度t。由式(17)可以得到a, 由式(13)可以得到d,將a和d代入式(9)可得到內力向量F〈J,K〉,由此可以求得車架結構與車身結構之間的耦合力fi。

3 剛度分析與承載度驗證

參考白車身彎曲剛度試驗加載方法[15],加載力F=3000N,加載點位于前座椅安裝位置左右對稱處,方向垂直向下(整體坐標系Y軸負方向)。根據式(9)、式(15)、式(16)和式(18)解得梁單元厚度ti以及車身與車架之間耦合力fi如表2所示,fi為正值表示其為拉力,為負值表示其為壓力。車架最終梁單元厚度t取3.686mm。

表2 車架梁單元厚度及車身與車架耦合力

耦合點撓度yi梁厚度ti(mm)耦合力fi(N)y13.560-169.680y24.110233.480y33.4905.116y43.74029.490y53.530-101.280

為了達到設定的承載度目標值,需要優化車身在耦合力fi(i=1,2,…,5)作用下的主斷面參數,車身結構采用基于線框的簡化幾何模型,簡化模型梁單元的中心線由車身A級曲面(圖10)上采集的特征點連接而成,并且遵照以下簡化原則:①先考慮車身主要承載梁結構,暫不考慮地板、頂蓋等覆蓋件;②暫不考慮焊點特性;③曲梁用多段直梁近似逼近。車身側圍及頂板的梁單元分布參考傳統車型,底板梁單元分布根據乘員座椅位置及與車架的接合位置確定,得到車身梁單元簡化模型如圖11所示。

圖10 車身A級曲面

圖11 車身梁單元分布及主斷面位置

在車身剛度鏈分析方法[16]中,只在前座椅中心處施加垂直于地面向下的載荷。對于本文研究的非承載式電動汽車車身,在進行車身的截面優化時,需要用f1、f2、f3、f4、f5代替彎曲工況下座椅安裝點處左右對稱的加載力F,將耦合力加載到車身與車架的接合點后的力學模型上,如圖12所示(圖12a為文獻[2]和文獻[17]中模型加載方式,圖12b為本文對應的加載方式)。車身主斷面簡化為薄壁矩形,如圖13所示。

[17]的車身剛度鏈計算方法,以輕量化為目標建立約束函數,優化設計模型如下:

(a)車身剛度鏈單元劃分及加載方式[17]

(b)非承載式車身剛度鏈單元劃分及加載方式圖12 車身剛度鏈單元劃分及加載方式

圖13 簡化梁截面

(19)式中,Δz為位于前座椅中心處的測量點最大撓度,其大小不應超過車身與車架整體變形時的約束值4.97×10-4m;ΔL為前懸支座處位移量。

利用MATALB遺傳算法工具箱計算各主斷面的最優截面尺寸,計算結果如表3所示。根據表3中的計算結果可以得到主斷面面積A、慣性矩Iy=∫Ay2dA、慣性矩Iz=∫Az2dA和極慣性矩Ip=∫A(y2+z2)dA在車身主斷面節點上的分布情況,如圖14所示。

根據圖14可以看出,編號為3、8、10、12的梁對車身剛度影響較為顯著。這與此類車身梁單元剛度靈敏度分布情況[17]吻合,說明優化結果合理。

為驗證上述方法的準確性,利用ANSYS分別建立車身、車架以及車身與車架耦合的beam188梁單元模型,彎曲工況下變形云圖見圖15~圖17,車身與車架耦合整體扭轉工況下的變形云圖見圖18。

根據位移云圖得到前座椅中心測量點處變形量如表4所示,測量點處的撓度在容許的變形范圍內,扭轉工況下前懸支座處位移ΔL=2.627mm<2.71mm,滿足扭轉性能要求。

表3 車身主斷面截面參數 mm

圖14 主斷面面積、慣性矩和極慣性矩分布情況

圖15 彎曲工況下車身變形

圖16 彎曲工況下車架變形

圖17 彎曲工況下車身與車架耦合變形

圖18 扭轉工況下車身與車架耦合變形

根據基于剛度的承載度評價方法,車架的剛度與車身及車架總體剛度的比值即為車架的承載度[18],根據式(1)分別得到車架的彎曲剛度K1=2.803MN·m2、車身和車架耦合整體的彎曲剛度K=3.84MN·m2。可以計算出車架彎曲剛度與整體彎曲剛之比即承載度為0.73,與設定的目標承載度非常接近,滿足設計要求。

表4 測量點撓度值 mm

4 結語

本文根據LifeDrive模塊化車架的載荷情況,拓撲優化后得到了電動汽車車架結構。建立了車架的回傳射線矩陣計算模型,當車架梁結構變動時,只需增加或減少相應的方程即可,避免了基于傳遞矩陣法的剛度鏈方法難以確定車架載荷的傳遞路徑和計算量大的缺陷,便于建立車架與車身之間的耦合計算模型。本方法用于非承載式車身整體的優化設計,明顯拓展了車身剛度鏈設計方法的應用范圍,提高了計算效率。最后,利用有限元方法計算驗證了實例車身、車架簡化模型剛度與承載度的分配滿足預先設定的目標值要求。

參考文獻:

[1]SchelkleE,ElsenhansH.VirtualVehicleDevelopmentintheConceptStageCurrentStatusofCAEandOutlookontheFuture[C]//3rdMSCWorldwideAerospaceConference&TechnologyShowcase.Toulouse,2001:24-26.

[2] 田海豹. 基于剛度鏈方法的車身概念設計研究[D]. 長沙:湖南大學, 2013.

[3] 扶原放, 金達鋒, 喬蔚煒. 微型電動車車架優化設計研究[J]. 機械制造, 2009, 47(1):12-15.FuYuanfang,JinDafeng,QiaoWeiwei.MiniResearchElectricVehicleFrameOptimizationDesign[J].MachineManufacturing, 2009, 47(1):12-15.

[4] 高云凱, 邵力行, 張海華. 微型電動車非承載式車身輕量化研究[J]. 汽車工程, 2008, 30(9):808-810.GaoYunkai,ShaoLixing,ZhangHaihua.DevelopmentofSkeletonStructurefortheFrame/bodyConstructionofaMiniElectricVehicle[J].AutomotiveEngineering, 2008, 30(9):808-810.

[5] 周姍姍. 基于車身承載度的半承載式客車結構改進研究[D]. 長春:吉林大學,2011.

[6]CavazzutiM,MerullaA,BertocchiE,etal.AdvancedHighPerformanceVehicleFrameDesignbyMeansofTopologyOptimization[C]//Internati-onalConferenceonStructuralEngineering,MechanicsandComputation.CapeTown,2010:279.

[7]HodkinsonR,FentonJ.LightweightElectric/HybridVehicleDesign[M].Atlanta:Elsevier, 2001:36-70.

[8] 梁晨, 何寧寧, 姚俊賢. 非承載式車身對整車剛度貢獻率研究[J]. 機械強度, 2009, 31(6):887-891.LiangChen,HeNingning,YaoJunxian.StudyonStiffnessContributionofSeparateTypeBodytoWholeFrame-bodyStructure[J].MechanicalStrength, 2009, 31(6):887-891.

[9] 喬蔚煒, 金達鋒, 于興林. 微型電動車車身結構優化設計中性能指標的確定方法[J]. 天津汽車, 2008(8):36-38.QiaoWeiwei,JinDafeng,YuXinglin.ResearchonPerformanceDeterminationMethodofaMiniElectricCarBody[J].TianjinAuto, 2008(8):36-38.

[10]PorscheEngineeringServices.Ulsab-Avc-PesEngineeringReport[R]. 2001:731-800.

[11]HifnieMAI.DesignandAnalyzetheChassisofanElectricVehicleforUseinCampusCondition[J].Utem, 2009(1):46-54.

[12] 扶原放, 金達鋒, 喬蔚煒. 微型電動車車架結構優化設計方法[J]. 機械工程學報, 2009, 45(9):210-213.FuYuanfang,JinDafeng,QiaoWeiwei.OptimazationMethodforaMiniElectricVehicleFrameStructure[J].JournalofMechanicalEngineering, 2009, 45(9):210-213.

[13] 郭永強. 回傳射線矩陣法的理論及其應用[D]. 杭州:浙江大學, 2008.

[14]McGuireW,GallagherRH,ZiemianRD.MatrixStructuralAnalysis[M]. 2nded.Amazon:CreateSpace, 2000.

[15] 林智平.JL030轎車白車身結構優化研究 [D]. 長春: 吉林大學, 2011.

[16] 劉子建, 周小龍, 田海豹, 等. 基于主斷面剛度優化分配的車身正向概念設計[J]. 中國機械工程,2015,26(6):837-843.LiuZijian,ZhouXiaolong,TianHaibao,etal.ForwardConceptualDesignofCarBodyUsingStiffnessOptimalAllocationofMainSections[J].ChinaMechanicalEngineering, 2015,26(6):837-843.

[17]EuropeanAluminiumAssociationStiffnessRelevanceandStrengthRelevanceinCrashofCarBodyComponents[R].Nordrhein-Westfalen:UniversityofAachen, 2010.

[18] 林標華, 姚成. 基于剛度的客車車身承載度評價[J]. 客車技術與研究, 2014, 36(4):7-9.LinBiaohua,YaoCheng.EvaluationofCoachBodyworkLoadingDegreeBasedonStiffness[J].Bus&CoachTechnologyandResearch, 2014, 36(4):7-9.

(編輯 王旻玥)

AnalysisMethodofReverberation-rayMatrixStiffnessChainforElectricVehicleBodies

ChenYiqiangLiuZijian

StateKeyLaboratoryofAdvancedDesignandManufactureforVehicleBody,HunanUniversity,Changsha, 410082

BasedonLifeDrivepureelectricvehiclemodularstructures,apureelectricvehicleframestructurewasobtainedbytopologyoptimizationmethod.Thereverberation-raymatrixmodelofframestructurewasestablishedbasedonthebeamelements.Accordingtothereverberation-raymatrixcalculationmodel,thethicknessesoftheframebeamelementsandtheforcesbetweenthebodyandframewereobtained,themainsectionalparametersofsimplifiedmodelofvehiclebodywasoptimizedwithamorereasonablemethod.Theloadingdegreebetweenbodyandframewascalculatedbyfiniteelementmethod.Resultsshowthattheframeandthebodystructuresmeettherequirementsofthetargetloadingdegreeandstiffnessperformance,andthismethodprovidesreferencetothebodyandframeconceptualdesignofelectricvehicle.

pureelectricvehicle;framedesign;loadingdegree;methodofreverberation-raymatrix

2016-01-13

國家重點基礎研究發展計劃(973計劃)資助項目(2010CB328002);國家自然科學基金資助項目(51475152)

U463.82

10.3969/j.issn.1004-132X.2016.23.021

陳毅強,男,1991年生。湖南大學機械與運載工程學院碩士研究生。主要研究方向為電動汽車車身結構正向設計方法與結構優化。劉子建(通信作者),男,1953年生。湖南大學機械與運載工程學院教授、博士研究生導師。

猜你喜歡
優化結構設計
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
瞞天過海——仿生設計萌到家
藝術啟蒙(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
論《日出》的結構
主站蜘蛛池模板: 无码人妻免费| 播五月综合| 色偷偷男人的天堂亚洲av| 日韩免费中文字幕| 亚洲Av激情网五月天| 亚洲天堂成人| 亚洲香蕉久久| 91小视频在线播放| 免费看a级毛片| 青青草原国产| 成人毛片免费观看| 青青草原偷拍视频| 国产在线小视频| 欧美亚洲激情| 亚洲日本在线免费观看| 伊人久综合| 国产人成网线在线播放va| 国产视频大全| 免费在线观看av| аⅴ资源中文在线天堂| 波多野结衣在线一区二区| www.精品国产| 国产网友愉拍精品视频| 欧美一道本| 国产网友愉拍精品视频| 在线观看国产一区二区三区99| 亚洲AV无码乱码在线观看裸奔| 一级黄色片网| 日韩国产精品无码一区二区三区| h视频在线观看网站| 精品国产成人a在线观看| 国产精品久久自在自线观看| 国产精品不卡片视频免费观看| 国产激情无码一区二区APP| 久久精品中文字幕少妇| 婷婷亚洲视频| 国产手机在线观看| 91久久精品日日躁夜夜躁欧美| 毛片免费高清免费| 老司国产精品视频91| 久久99精品国产麻豆宅宅| 亚洲精品日产AⅤ| 国产青青草视频| 91精品最新国内在线播放| 在线色国产| 色视频国产| 国产女同自拍视频| 免费在线观看av| 精品国产香蕉伊思人在线| 亚洲国产成人超福利久久精品| 国产福利在线免费观看| 国产一区自拍视频| 在线精品视频成人网| 欧美精品色视频| 久久永久免费人妻精品| 国产91小视频| 国产女人在线| 精品国产免费观看| 最新亚洲人成网站在线观看| 日韩精品欧美国产在线| 色综合中文综合网| 欧美色伊人| 国产精品永久免费嫩草研究院| 亚洲综合第一区| 亚洲AV无码一区二区三区牲色| 国产精品免费p区| 在线无码九区| 国产另类视频| 伊人激情久久综合中文字幕| 无码在线激情片| 国产另类视频| 波多野结衣中文字幕久久| 午夜福利免费视频| 色偷偷一区| 国产成人高清精品免费| 色综合婷婷| 午夜毛片免费观看视频 | 97久久人人超碰国产精品| 精品国产污污免费网站| 在线中文字幕日韩| 亚洲swag精品自拍一区| 国产视频入口|