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

基于gPC理論的不確定參數電動汽車脈沖響應研究

2021-09-27 07:05:08田國英張大偉易興利鄧鵬毅孫樹磊
振動與沖擊 2021年16期

田國英,張大偉,易興利,鄧鵬毅,孫樹磊

(1.西華大學 汽車與交通學院,成都 610039;2.西華大學 汽車測控與安全四川省重點實驗室,成都 610039;3.西華大學 四川省新能源汽車智能控制與仿真測試技術工程研究中心,成都 610039;4.長安大學 汽車學院,西安 710064;5.中車株洲電力機車有限公司,湖南 株洲 412001)

平順性是評價汽車性能的主要指標之一,其優劣對乘坐舒適性、零部件及路面的疲勞壽命和行駛安全性均有重要影響,如何提高平順性一直受到整車廠、各國研究機構和學者的重視[1-2]。

目前,汽車平順性動力學方面的相關研究已經取得了很大發展,但是,多數研究仍基于車輛結構參數為確定性的動力學模型[3-6]。汽車平順性隨機振動研究中主要以平穩高斯隨機路面不平作為輸入,較少考慮參數不確定的影響。實際上,受生產制造、長期運營磨損老化及運載條件變化等因素影響,車輛自身結構參數相對于設計值是隨機變化的,在研究汽車平順性時有必要考慮參數不確定的影響。考慮參數不確定的汽車動力學計算方法主要有蒙特卡洛(Monte Carlo,MC)法[7]、協方差法[8]、區間法[9]及廣義多項式混沌(generalized polynomial chaos,gPC)法[10]等。其中,MC法計算效率較低,協方差法僅能給出響應的統計值,區間法則需確定參數的具體范圍,但汽車大多數結構參數的分布范圍并不明確。

gPC法能有效緩解上述問題,近年來發展迅速。該方法源于Wiener[11]的均勻混沌理論,他最早給出用埃爾米特正交多項式組合描述高斯隨機過程的表達式。此后,Ghanem等[12]將該方法與有限元法結合,并在結構動力學領域開展了大量研究。在此基礎上,Xiu等[13-14]基于Askey框架,將埃爾米特多項式表示高斯隨機過程的方式拓展到其他隨機過程的正交多項式混沌展開,建立了正交多項式和多數概率分布間的對應關系,擴展為廣義多項式混沌理論。Sandu等[15-16]最早將gPC理論引入多體系統動力學和汽車動力學中,其研究團隊開展了包括汽車垂向動力學[17]、胎地相互作用[18]、操縱穩定性[19]等方面的不確定參數隨機振動分析。之后,學者們利用gPC理論在汽車及其零部件參激隨機振動方面開展了較多工作[20-24]。與汽車平順性相關的研究主要有:Sandu等建立了懸架剛度和阻尼非線性的1/4汽車垂向模型,通過計算不確定懸架阻尼、輪胎剛度和路面不平下的振動響應,與MC法對比,說明了gPC方法計算的準確性和高效性。Li等建立7自由度越野車輛動力學模型,分析了輕/重車情況下不確定結構參數和土壤濕度對車輛牽引能力的影響。Kewlani等借助1/4車輛垂向模型對比了gPC法、多元gPC法和MC法的差異,從理論和實驗兩方面驗證了gPC方法的正確性。凌鋒利用Adams/Car軟件建立某型汽車整車模型,與MATLAB軟件聯合仿真懸架剛度和阻尼不確定時的動力學響應統計值,與MC法對比驗證了PC法的正確性。Wu等提出了一種多項式混沌-切比雪夫區間的混合方法,以4自由度汽車垂向模型為例,驗證了方法在計算不確定參數影響中的準確性和高效性。上述工作已經證明了gPC理論在汽車動力學分析中的可行性和計算精度,但這些研究側重方法準確性的驗證,且其汽車平順性模型做了大量簡化,缺乏參數不確定對汽車平順性的影響程度及規律的研究。因此,基于gPC理論,在以往理論方法可行性及準確性驗證基礎上,利用完整汽車平順性模型,從工程應用角度開展關鍵參數不確定的定量影響研究顯得十分必要,可為合理預測汽車平順性響應范圍,并對關鍵零部件的設計和可靠性評估提供幫助。

本文將gPC理論運用于汽車平順性分析中,建立考慮電機參振的含不確定參數電動汽車的平順性模型。在與MC法對比驗證基礎上,分別在整車和電機總成兩類參數不確定條件下,對某型電動乘用車以不同車速通過脈沖型路面時的平順性和電機振動響應的統計指標開展細致的分析工作,獲得參數不確定對各個振動響應的影響規律和定量化程度,以期為汽車平順性方面的參激隨機振動分析提供參考。

1 基于gPC理論的不確定參數電動汽車平順性模型

1.1 考慮電機參振的電動汽車平順性模型

考慮電機參振的電動汽車平順性模型,如圖1所示。模型中,假設車體、人體、電機和車輪均為剛體。車體由4個懸架共同支撐,考慮垂向(Zb)、俯仰(βb)和側傾(Φb)運動3個自由度;電機通過3個懸置彈性連接至車體前部,考慮垂向(Ze)、俯仰(βe)和側傾(Φe)運動3個自由度;人體和車輪僅考慮垂向運動(Zp和Zw),共計11個自由度。圖1中:Mb,Mp,Me和Mw分別為車體、人體、電機和車輪質量;Ibx,Iby分別為車體側傾和俯仰慣量;Iex和Iey分別為電機側傾和俯仰慣量;Ks和Cs分別為懸架彈簧剛度和減振器阻尼;Kc和Cc分別為座椅剛度和阻尼;Kt為輪胎垂向剛度;Km和Cm分別為電機懸置垂向剛度和阻尼;Lfc和Lrc分別為質心距前軸和后軸的距離;Wf和Wr分別為前軸和后軸的輪距;其他尺寸參數如附錄A所示;Z0為路面不平度輸入。

圖1 考慮電機參振的電動汽車平順性模型Fig.1 Vertical dynamics model of an electric vehicle considering electric motor vibration

根據牛頓第二定律和車輛各部件受力分析,可得各自由度的動力學方程

(1)車體垂直運動

(1)

(2)車體俯仰運動

(2)

(3)車體側傾運動

(3)

(4)電機垂向運動

(4)

(5)電機俯仰運動

(5)

式中,Te為電機驅動扭矩,根據汽車驅動理論[25],可得

式中:G為整車重力;f為滾動阻力系數;α為坡度角;CD為空氣阻力系數;A為迎風面積;v為車速;δ為旋轉質量換算系數;r為車輪半徑;iz為傳動系總傳動比;ηt為傳動系傳動效率。根據上式可計算不同車速、坡度和加速條件下的電機驅動扭矩,進而影響電機俯仰振動。

(6)電機側傾運動

(6)

(7)人體垂向運動

(7)

(8)左前車輪垂向運動

(8)

(9)右前車輪垂向運動

(9)

(10)左后車輪垂向運動

(10)

(11)右后車輪垂向運動

(11)

式中:Ffl,Ffr,Frl,Frr為4個懸架垂向作用力;Fbm為車體對電機的垂向作用力;Fpb為座椅對車體的垂向作用力;Ft為車輪與路面作用力;以上各作用力的具體表達式可通過受力分析獲得,此處不再贅述。為真實模擬車輪與路面相互作用,考慮了當輪胎上跳過大時可能出現車輪跳離地面的情況,因此,車輪與路面作用力可表達為

(12)

1.2 廣義多項式混沌理論的應用

目前,基于gPC理論的計算方法主要有兩種,隨機伽遼金方法(stochastic Galerkin method,SGM)和隨機配點法(stochastic collocation method,SCM)。由于隨機伽遼金方法需要對原確定性模型進行不確定變量的混沌展開推導,工作量大,且不同的不確定變量需重新推導,不利于快速建模和計算,故隨機伽遼金方法一般用于自由度不大的簡單模型。相比之下,隨機配點法無需對原確定性模型進行修改,只需將原模型封裝成子程序開展計算,大大簡化了復雜模型的隨機分析工作。因此,本文利用隨機配點法開展不確定參數的影響計算。

對于式(1)~式(11)的各部件動力學方程,其系統動力學方程可總體表達為

(13)

當考慮系統存在參數不確定時,式(13)可變形為

(14)

式中:Z為總不確定參數向量;Z1,Z2,Z3和Z4分別為對應質量、阻尼、剛度和外部激勵的不確定參數向量,可展開表達為

(15)

式中,D1,D2,D3和D4分別為對應質量、阻尼、剛度和外部激勵的不確定參數數量,記不確定參數總數D=D1+D2+D3+D4。

根據gPC理論,各類不確定參數可表達為多項式混沌展開形式,即

(16)

(17)

(18)

(19)

類似的,系統輸出響應也可表達為多項式混沌展開形式

(20)

(21)

(22)

式中:i=(i1,…,ij,…,iD),j=1,…,D,為按正交多項式階數的分級詞典順序排列的多重指標,ij為第j個隨機變量的多項式階數,|i|=i1+…+iD為其指標值;Φi(Z)為不同組合下各不確定參數正交多項式的乘積,即

Φi(Z)=φi1(Z11)…φiD(Z4D4)

(23)

H為D個不確定參數的多項式組合數,記Q=max(Q1,Q2,Q3,Q4),則

(24)

對于隨機伽遼金法來說,其通過將式(16)~式(22)代入式(14),把原確定性方程組在隨機域擴維后獲得一組新的包含不確定參數的方程組,從而開展計算,其對復雜動力學系統的推導過程繁瑣耗時。隨機配點法則首先在式(15)構成的隨機空間中選取節點,節點的選取依據高斯求積法則確定,設選取的節點集為ΘN={Z(n)},N=N1×…×Nj×…×ND為節點集總數,Nj為第j個不確定參數的節點數,n=1,…,N。將ΘN代入式(16)~式(19),進一步代入式(14)可得

(25)

(26)

(27)

(28)

此外,為獲得較準確的響應概率密度及分布,可通過對式(20)~式(22)中的不確定參數向量Z進行抽樣得到響應的大樣本數據,然后擬合獲得其概率密度及相關概率統計;響應的變異系數εX可利用μX和σX的比值計算

(29)

2 不確定參數對電動汽車平順性的影響分析

2.1 不確定參數的選取

針對所建立的電動汽車平順性模型,以某型電動乘用車為例,將待分析的不確定參數分為兩種工況:一種是對汽車平順性有重要影響的整車參數;另一種是對電機振動有重要影響的電機總成參數,不確定參數的選取情況,如表1所示。

表1 不確定參數的選取Tab.1 Selection of uncertain parameters

從表1可知:在整車平順性的主要影響參數中,考慮車體質量、懸架剛度和阻尼、輪胎剛度等4種參數隨機變化。其中:考慮車輛運行中乘客Mp及行李質量Ml可能發生變化,且車體處于各種質量狀態的概率并無差別,因此選取車體質量服從均勻分布,其上、下限分別為Mbu和Mbl,質量上限Mbu=Mb+4Mp+5Ml,Ml為每位乘客的行李質量,質量下限Mbl=Mb;其余不確定參數均服從正態分布,其均值取設計值,方差在均值的10%~30%。

在電機總成參數中,考慮電機質量、懸置剛度和阻尼等3種參數隨機變化。其中,考慮到設計過程中電機選型的不同,且電機選型時處于各種質量狀態的概率并無差別,因此選取電機質量服從均勻分布,其上、下限分別為Meu和Mel,分別取Meu=150%Me和Mel=50%Me;其余不確定參數服從正態分布,其均值取設計值,方差在均值的10%~30%。

此外,根據GB/T 4970—2009《汽車平順性試驗方法》規定,路面不平度選取典型的脈沖路面輸入,其波形如圖2所示。車速以10 km/h為間隔在10~60 km/h變化。

圖2 脈沖型路面不平度輸入(mm)Fig.2 Road impulse input (mm)

為驗證基于SCM建立的不確定模型的正確性,在以上兩類不確定參數工況下,將其與MC法的計算結果進行了對比,如圖3和圖4所示。計算中,車速取30 km/h,正態分布的不確定參數變異系數均取20%,MC法的樣本數取10 000組。

圖3 車體加速度響應的對比Fig.3 Comparison of body acceleration responses

圖4 電機加速度響應的對比Fig.4 Comparison of motor acceleration responses

由圖3和圖4可知:兩種工況下,車體和電機加速度均值響應均與MC法計算結果完全重合,方差響應十分接近,響應最大值的概率密度曲線也幾乎重合,從而說明了本文所建模型的正確性。此外,圖中還給出了相應的確定性系統響應和標準正態分布的概率密度。可以看出,確定性系統響應與不確定系統響應在峰值位置處存在較明顯差異,車體和電機加速度最大值的概率密度與標準正態分布較為接近。兩種方法計算效率的對比數據,如表2所示。由表2可知:相比MC法,SCM利用較低的多項式階數和較少組合數便可快速計算出準確結果,當不確定參數較少時,這一優勢更為明顯。

表2 SCM和MC法計算效率對比Tab.2 Comparison of computational efficiency between SCM and MC method

基于以上不確定參數選取和模型驗證,下面分別在兩種工況下對不確定參數對汽車平順性的影響進行詳細計算和分析,SCM的計算參數與表2相同。

2.2 整車參數對平順性的影響(工況1)

根據表1,計算了工況1條件下的汽車平順性響應,并進行了相關統計分析。60 km/h條件下,當正態分布的不確定參數變異系數均為30%時的車體和人體加速度、懸架壓縮量和輪胎載荷的均值和方差響應,以及其概率密度隨時間的演化特征,如圖5所示。

圖5 整車參數不確定條件下的汽車平順性響應Fig.5 Vertical dynamics responses under whole vehicle parametric uncertainties

由圖5可知:不確定參數影響下各響應的動態分布范圍清晰可見,在峰值處的方差較大。其中,懸架壓縮量受參數不確定的影響較大,各時刻的響應方差均較大。由概率密度演化結果可明顯看出大概率響應值的范圍,車體加速度和輪胎載荷的大概率響應值基本分布在均值附近,而懸架壓縮量的大概率響應值則與其均值有明顯差異。也就是說,懸架壓縮量的均值和方差不能完全表征參數不確定情況下的響應特征,還應結合概率密度演化來補充分析。

為細致考察輸入變異系數和車速對各動力學響應的影響,各動力學均值響應最大值的變異系數隨輸入變異系數和車速的變化規律,如圖6所示。一般來說,響應變異系數越大,說明不確定參數對該響應的影響越明顯。

圖6 平順性響應最大值變異系數的變化規律Fig.6 Changing rules of variation coefficients for maximums of vertical dynamics responses

由圖6可知:總體來看,隨輸入變異系數的增大,各響應的變異系數不斷增大,且隨車速提高,各響應的變異系數總體呈現增大趨勢;懸架壓縮量的變異系數最大,車體和人體加速度次之,輪胎載荷變異系數最小。具體來看,車體加速度變異系數較人體加速度大,二者最大可達34%和26%。前懸架壓縮量的變異系數略大于后懸架,隨輸入變異系數增大,差別越明顯;在不同輸入變異系數下,懸架壓縮量的最大變異系數分別為10%,20%和47%,響應變異系數甚至大于輸入的變異系數;此外,車速對懸架壓縮量響應的變異系數影響不顯著。與懸架壓縮量類似,前輪載荷的變異系數略大于后輪,但車速對輪胎載荷的變異系數影響較明顯,尤其在大的輸入變異系數情況下,幾乎隨車速提高線性增大;在不同輸入變異系數下,輪胎載荷的最大變異系數分別為6%,12%和17%。

基于各響應最大值的概率密度,其分布范圍、均值和方差隨輸入變異系數和車速的變化規律,如圖7所示。

圖7 平順性響應最大值的均值、方差和分布范圍變化規律Fig.7 Changing rules of means,variances and distribution ranges for maximums of vertical dynamics responses

由圖7可知:總體來看,隨輸入變異系數增大,各響應最大值的均值變化不大,但其方差和分布范圍均有明顯提高,說明輸入變異系數主要影響方差和分布范圍,對均值影響不大。隨車速提高,車體加速度和懸架壓縮量最大值的均值總體略有下降,輪胎載荷則明顯增大;車體加速度和輪胎載荷方差和分布范圍總體不斷增大,懸架壓縮量則略有下降。

根據各動力學響應最大值的概率密度,圖8給出其概率分布,并基于此計算了各動力學響應最大值超過設定閾值的概率。其中,車體加速度、懸架壓縮量、輪胎最大載荷和最小動載的閾值分別設置為0.5g,0.14 m(懸架設計行程)、15.5 kN(3倍靜輪載)和0(離地)。

由圖8可知:總體來看,隨輸入變異系數的降低,各響應概率分布曲線越來越陡峭,這與前述的響應最大值的分布范圍變窄一致。從超限概率來看,不同輸入變異系數下車體加速度超過0.5g的最大概率分別為53%,47%和45%,隨輸入變異系數增大超限概率最大值略有下降;除10 km/h車速外,車體加速度超限概率總體隨車速提高而降低,反映出在小變異系數和低速時的車體加速度最大值略微偏大;懸架壓縮量超過設計值的最大概率分別為14%,29%和36%,隨變異系數提高明顯增大,相反隨車速提高而降低,說明在大變異系數和低速時懸架壓縮量較大;輪胎載荷最大值超過3倍靜載的最大概率分別為8%,23%和30%,其超限概率隨變異系數和車速提高均明顯增大,說明在大變異系數和高速時的輪胎動載較高;輪胎載荷最小值離地最大概率分別為64%,54%和42%,除了車速10 km/h車輪不發生離地外,離地概率總體上隨變異系數和車速提高而下降,說明在小變異系數和低速時車輪更容易與地面脫離。

圖8 平順性響應最大值的概率分布及超限概率Fig.8 Probability distributions and overrun probabilities for maximums of vertical dynamics responses

2.3 電機總成參數對電機振動的影響(工況2)

根據表1,計算了工況2條件下的電機振動響應,并給出了相關統計分析。60 km/h條件下,當正態分布的不確定參數變異系數均為30%時的電機加速度和懸置壓縮量的均值和方差響應,以及其概率密度隨時間的演化特征,如圖9所示。

圖9 電機總成參數不確定條件下的電機振動響應Fig.9 Motor dynamics responses under motor assembly parametric uncertainties

由圖9可知:不確定參數影響下各響應的動態分布范圍清晰可見,在峰值處的方差較大。其中,電機懸置壓縮量受參數不確定的影響較大,各時刻的響應方差均較大。由概率密度演化結果可明顯看出大概率響應值的范圍,電機加速度的大概率響應值基本分布在均值響應附近,而懸置壓縮量的大概率響應值則與其均值有明顯差異。同樣,對于懸置壓縮量來說,還應結合概率密度演化來補充分析參數不確定情況下的響應特征。

為細致考察輸入變異系數和車速對電機振動響應的影響,電機振動均值響應最大值的變異系數隨輸入變異系數和車速的變化規律,如圖10所示。

圖10 電機振動響應最大值變異系數的變化規律Fig.10 Changing rules of variation coefficients for maximums of motor dynamics responses

由圖10可知:總體來看,隨輸入變異系數的增大,各響應的變異系數不斷增大,且隨車速提高,各響應的變異系數總體呈現增大趨勢;懸置壓縮量響應的變異系數明顯大于電機加速度,說明輸入參數不確定對懸置壓縮量影響更大。具體來看,在不同輸入變異系數情況下,電機加速度和懸置壓縮量最大值的變異系數可分別達到5%,9%,15%和14%,24%,51%,懸置壓縮量最大值的變異系數甚至大于輸入的變異系數。

基于各均值響應最大值的概率密度,其分布范圍、均值和方差隨輸入變異系數和車速的變化規律,如圖11所示。

圖11 電機振動響應最大值的均值、方差和分布范圍變化規律Fig.11 Changing rules of means,variances and distribution ranges for maximums of motor dynamics responses

由圖11可知:總體來看,隨輸入變異系數增大,各響應均值最大值變化不大,但其方差和分布范圍增大,同樣說明電機參數不確定主要影響響應方差和分布范圍,對均值影響不大。隨車速提高,電機加速度最大值均值總體下降,其方差和分布范圍總體增大;懸置壓縮量最大值均值、方差和分布范圍隨車速變化不明顯。

根據電機各響應最大值的概率密度,圖12給出其概率分布,并基于此計算了電機振動響應最大值超過設定閾值的概率。其中,電機加速度、懸置壓縮量的閾值分別設置為0.9g和15 mm。

圖12 電機振動響應最大值的概率分布及超限概率Fig.12 Probability distributions and overrun probabilities for maximums of motor dynamics responses

由圖12可知:總體來看,隨輸入變異系數的降低,各響應概率分布曲線越來越陡峭,這與前述的響應最大值的分布范圍變窄一致。從超限概率來看,不同輸入變異系數下電機加速度超過0.9g的最大概率分別為75%,71%和65%,隨輸入變異系數增大超限概率最大值略有下降,除10 km/h車速外,超限概率總體隨車速提高而降低,說明在小變異系數和低速時的電機加速度偏大;懸置壓縮量超過設定值的最大概率分別為68%,62%和55%,隨變異系數提高也略有降低,隨車速提高超限概率也總體下降,說明在小變異系數和低速時懸置壓縮量較大,其與電機加速度的規律類似。

3 結 論

本文基于廣義多項式混沌理論和汽車動力學理論,建立了考慮電機參振的含不確定參數電動汽車的平順性分析模型。以某型電動乘用車為例,計算分析了不同輸入變異系數和車速條件下,車輛通過脈沖路面時兩類參數不確定對各振動響應的影響程度和規律,得到如下結論:

(1)在電動汽車整車平順性模型中,相比傳統MC法,運用gPC理論能更快速有效度量參數不確定對各統計指標的影響程度。考慮參數的不確定后,對動力學響應的量值和規律分析更加全面,可為合理預測響應范圍并對關鍵零部件的設計和可靠性評估提供建議。

(2)在所分析速度域內,隨參數不確定程度的增加,各響應最大值的均值基本不變,但其方差、變異系數和分布范圍均不斷增大。當不確定參數變異系數為30%,車速為60 km/h時,各響應最大值的變異系數達到最大,部分響應的變異系數甚至超過輸入變異系數。

(3)不同參數變異系數和車速條件下,超限概率的變化規律具有多樣性,與所設置響應限值的大小及該響應隨車速的變化規律有關。在小變異系數和低速時,車體加速度、輪胎載荷最小值、電機加速度和懸置壓縮量超限概率偏大,在大變異系數和低速時懸架壓縮量超限概率較大,而在大變異系數和高速時的輪胎載荷最大值超限概率較高。

附錄A

表A.1 考慮電機振動的電動汽車平順性確定性模型質量及慣量參數Tab.A.1 Mass and inertia of the determined vertical dynamics model of electric vehicle considering electric motor vibration

表A.2 考慮電機振動的電動汽車平順性確定性模型剛度及阻尼參數Tab.A.2 Stiffness and damping of the determined vertical dynamics model of electric vehicle considering electric motor vibration

表A.2(續)

表A.3 考慮電機振動的電動汽車平順性確定性模型幾何尺寸及其他參數Tab.A.3 Geometric and other parameters of the determined vertical dynamics model of electric vehicle considering electric motor vibration

主站蜘蛛池模板: 亚洲一级毛片在线播放| 亚洲黄色激情网站| 欧美日韩午夜| 一区二区三区国产精品视频| 风韵丰满熟妇啪啪区老熟熟女| 精品成人一区二区| 免费人成在线观看成人片| 亚洲成人77777| 精品一区二区三区波多野结衣| 人妻91无码色偷偷色噜噜噜| 91在线国内在线播放老师| 一区二区理伦视频| 91黄色在线观看| 国产人成在线视频| 看你懂的巨臀中文字幕一区二区| 日本91在线| 国产成人综合日韩精品无码不卡| 扒开粉嫩的小缝隙喷白浆视频| 日本一区高清| 九九九精品成人免费视频7| 9丨情侣偷在线精品国产| 亚洲无码熟妇人妻AV在线| 青青草原偷拍视频| 国产精品夜夜嗨视频免费视频| 国产精品尤物铁牛tv| 成人在线视频一区| 亚洲成A人V欧美综合| 制服丝袜在线视频香蕉| 91啪在线| 丁香亚洲综合五月天婷婷| 国产精品久线在线观看| 毛片网站在线播放| 广东一级毛片| 亚洲香蕉伊综合在人在线| 国产精品美女网站| 国产高清精品在线91| 日韩毛片免费| 一级一级特黄女人精品毛片| 午夜丁香婷婷| WWW丫丫国产成人精品| 国产在线精品99一区不卡| 国产精品熟女亚洲AV麻豆| 欧美日韩免费| 亚洲欧美国产五月天综合| 欧美专区日韩专区| 久久久久亚洲av成人网人人软件| 四虎永久在线精品影院| 欧美区一区| 91精品aⅴ无码中文字字幕蜜桃| 亚洲侵犯无码网址在线观看| 萌白酱国产一区二区| 国产精品原创不卡在线| 亚洲不卡av中文在线| 免费人成黄页在线观看国产| 在线中文字幕日韩| 综合色亚洲| 日韩国产黄色网站| 久草网视频在线| 美女毛片在线| 久久亚洲黄色视频| 成年免费在线观看| 欧美一级视频免费| 国产成人免费| 亚洲伊人天堂| www中文字幕在线观看| 色噜噜狠狠狠综合曰曰曰| 国产丝袜一区二区三区视频免下载| 四虎永久免费地址在线网站| 国产精品片在线观看手机版| 在线观看国产一区二区三区99| 亚洲人网站| 久久77777| 日本爱爱精品一区二区| 人妻无码中文字幕第一区| 国产成年女人特黄特色毛片免| 国产伦精品一区二区三区视频优播| 99re66精品视频在线观看| 看你懂的巨臀中文字幕一区二区| 亚洲三级电影在线播放| 久久五月天综合| 91香蕉视频下载网站| 午夜在线不卡|