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

基于多工況并行任務的摩托車懸架參數多目標優化

2013-09-08 07:22:22徐中明楊建國張志飛李仕生李曉靈
振動與沖擊 2013年4期
關鍵詞:振動優化模型

徐中明,楊建國,張志飛,李仕生,李曉靈

(1.重慶大學 機械傳動國家重點實驗室,重慶 400030;2.重慶大學 機械工程學院,重慶 400030;3.重慶建設摩托車股份有限公司,重慶 400054)

摩托車懸架系統是整車的重要組成部分,對整車的行駛平順性、操縱穩定性、騎乘的舒適性等有著很大影響[1]。用戶對摩托車綜合性能要求的不斷提高,國內外對摩托車懸架系統的研究越來越受到重視。當前國內外對摩托車懸架系統的研究較多的集中在懸架系統零部件性能改善[2-3]、單一工況下的性能提升等方面,從人-車系統和整車綜合性能角度對懸架系統進行設計研究的還較少。在平順性方面,當前大部分研究都是以隨機路面工況下降低垂向振動為目標[4-5],并沒有考慮俯仰振動的影響,由于摩托車的懸掛質量分配系數ε遠小于1,前后懸架系統具有高度的耦合性,因此,車身局部垂向振動的改善可能導致俯仰振動的加劇,從而使整車舒適性變差;此外,隨機路面工況下平順性的改善可能導致其它工況下舒適性的降低,由于摩托車前懸架剛度較小,加速制動時俯仰角太大會導致前叉管撞擊限位,從而影響整車舒適性。半主動懸架能夠適應復雜多變的行駛工況,對于改善摩托車的垂向和俯仰振動能夠取得良好的效果[6-8],并且能夠保證其它工況下的整車舒適性,但由于成本昂貴,在摩托車上的應用受到限制,目前在國內的摩托車上未見使用。俯仰和垂向振動是影響車輛乘坐舒適性的關鍵問題[9],國內外從考慮俯仰振動角度對摩托車平順性研究的文獻很少,而從多工況角度對摩托車懸架系統進行匹配研究的文獻還未有報道。

本文針對某125摩托車存在的制動俯仰角過大的問題對懸架系統進行優化匹配,通過懸架特性參數的合理匹配實現制動俯仰角的降低,避免發生前懸架撞擊限位提高騎乘舒適性,同時保證隨機路面工況下的平順性能。采用摩托車動力學專用軟件BikeSim建立包含駕駛員在內的人-車系統仿真模型,以加速/制動工況下的車身俯仰角最大值和B級路面工況下的俯仰、垂向振動加速度均方根值為優化目標,在iSIGHT中建立兩種工況的并行優化任務,集成BikeSim、MATLAB進行多目標優化,取得了明顯效果。

1 加速/制動與平順性仿真

1.1 BikeSim中人-車系統模型的建立

整車模型的建立。BikeSim是一種基于數學模型的面向特性的參數化建模的動力學軟件,該模型共包含10個剛性體29個自由度,將整車分為車體、轉向、輪胎、懸架、制動、傳動等子系統,通過輸入各子系統的特性參數或曲線進行建模,其整車模型如圖1所示。

圖1 摩托車整車模型Fig.1 Vehicle model of motorcycle

由于摩托車后懸架結構形式多樣,BikeSim將后懸架統一簡化為搖臂機構,通過搖臂杠桿比來表達不同的懸架結構。對于制動系,BikeSim用杠桿比來表達制動操縱機構,將操縱機構對制動系輸入的制動力簡化為制動杠桿力,制動杠桿力經過增壓器、主缸、比例閥、輪缸等一系列部件的傳遞作用最終形成車輪上的制動力矩,由于本文只是研究相同制動力條件下不同懸架特性參數的車身俯仰角,只需制動系產生相同的減速度,因此制動系可采用BikeSim自帶的制動系模型。BikeSim提供的輪胎模型包括:“魔術公式”、查表輪胎模型和第三方輪胎模型,本文采用的是“魔術公式”。整車模型中制動系、傳動系和輪胎對本文的研究不產生影響,考慮到單獨建模的復雜性,這里采用了BikeS-im自帶的模型,其他子系統均是根據對象車的實際結構提取的參數建立。建立的某125摩托車部分參數如表1所示。

表1 某125摩托車部分參數Tab.1 Some parameters of a 125 motorcycle

駕駛員模型的建立。由于摩托車駕駛員的質量與摩托車整備質量相當,因此在整車建模時必須考慮駕駛員的影響。當前大部分研究中對于駕駛員模型過于簡化,并不能全面的反映整車動力學特性;有的甚至忽略了駕駛員的影響,不符合真實情況。BikeSim中將駕駛員質量分為上下半身兩部分,考慮了上下半身的質量、轉動慣量、尺寸、姿勢等關鍵因素,上下半身通過轉動副連接,上半身可以相對于車架做側傾和側向運動。按照NASA的建議駕駛員上半身質量占62%,下半身質量占38%進行質量分配。這里取駕駛員質量為70 kg,建立的駕駛員模型如圖2所示。

圖2 駕駛員模型Fig.2 Driver model

1.2 路面模型的建立

加速制動仿真路面的建立。在BikeSim的道路建模部分根據GB20073的要求建立平直路面,設置道路附著系數為0.85,建立的平直路面如圖3所示。

圖3 平直路面Fig.3 Flat straight road model

平順性仿真路面的建立。摩托車是單軌車輛,平順性仿真時只需在BikeSim中建立道路高程隨路長變化的二維路面模型。按照GB/T7031對于道路分級的要求以及道路平度的表示方法,本文采用諧波疊加法[10-11]對 B級路面的平度進行模擬。

式中:θ為[0,2π]上的隨機數;x為道路的x方向。通過離散道路x方向的值,就可以得到隨機路面z方向的值。

根據以上原理,采用MATLAB編制B級路面程序,得到的B級路面的平度數據導入BikeSim中,如圖4。

圖4 B級路面不平度Fig.4 Grade B road stochastic excitation

1.3 加速/制動工況的設定與仿真計算

圖5 速度變化曲線Fig.5 Target speed curve

在平直路面上進行加速/制動工況的仿真,加速工況設定為0~5 s速度從0均勻增加到60 km/h,之后進入勻速行駛;從第15 s開始制動,為了獲得更好的制動效果,采用聯動制動的方式給前后輪同時添加制動力。作用在摩托車上的制動力的大小與眾多因素有關,包括操縱機構的杠桿比、液壓傳動系統、地面附著情況等,但制動力與制動系統并非本文研究重點,這里只是為了在附著良好的情況下獲得滿足要求的制動加速度,因此只要添加合適的制動力即可。在BikeSim中需要添加的是制動杠桿力(即制動踏板力與杠桿比的乘積),不需要考慮復雜的制動操縱機構,這里給前后輪添加的制動杠桿力分別為120 N和80 N。速度和制動力曲線如圖5、圖6。

設定仿真時間為25 s,運行仿真,得到的速度、俯仰角和前后懸架變形曲線如下圖所示。雖然制動時俯仰角加速度也會對行駛舒適性產生一定影響,但是由于俯仰角加速度主要受制動力大小的影響,懸架特性參數的改變對其降低效果并不明顯,因此,這里沒有給出俯仰角加速度的曲線。

從圖7可以發現,在制動力作用下17.5 s時速度降低為0,前進已經停止,計算可知制動減速度大小為6.67 m/s2,完全滿足GB20073對于L3類摩托車制動性能的要求,說明添加的制動杠桿力大小是合適的。

從圖8可以看出加速階段俯仰角較小,俯仰角最大值出現在制動即將停止時,最大值為6.4°,制動階段的車身俯仰運動近似于整車繞后軸的旋轉,會導致前叉管發生壓縮運動,由圖9可以看出前叉管壓縮最大值為199.7 mm,而前叉管的行程僅為105 mm,因此必然會發生撞擊限位而影響整車的騎乘舒適性。

圖6 前后制動杠桿力曲線Fig.6 The front and rear break lever forces

圖7 加速/制動時速度曲線Fig.7 Speed curve during acceleration-braking

圖8 加速/制動時俯仰角曲線Fig.8 Pitch angle curve during acceleration-braking

圖9 加速/制動時前后懸架變形Fig.9 Suspension deformation during acceleration-braking

圖10 車身俯仰振動角加速度功率譜密度Fig.10 PSD of body pitching vibration acceleration

1.4 平順性工況的設定與仿真計算

B級路面下的平順性仿真。設定摩托車以車速60 km/h等速行駛于B級路面上,仿真時間90s,將計算結果導入MATLAB進行后處理得到俯仰和垂向振動加速度均方根值分別為 5.12 rad/s2、1.69 m/s2,圖 10、圖11分別是車身俯仰、垂向振動加速度功率譜密度。

圖11 車身質心垂向振動加速度功率譜密度Fig.11 PSD of sprung mass vertical vibration acceleration

2 多工況并行任務多目標優化的實現

2.1 問題描述

為了實現該車在加速/制動時行駛舒適性的改善,避免前叉管撞擊限位,同時保證隨機路面的平順性,防止俯仰或垂向振動的惡化;以前、后懸架的剛度、阻尼為設計變量,以加速/制動時俯仰角最大值和隨機路面工況下車身質心俯仰、垂向振動加速度均方根值為優化目標,采用NSGA-Ⅱ算法進行多目標優化。通過前后懸架的剛度、阻尼的優化匹配實現該車在這兩種工況下綜合性能的提高。

該優化問題的數學模型可描述為:式中:fyjd是加速/制動工況下俯仰角最大值,rmsp、rmsv分別是俯仰和垂向振動加速度均方根值,Front_K、Front_C、Rear_K、Rear_C分別是前、后懸架的剛度、阻尼。

2.2 多工況并行任務優化平臺的搭建

BikeSim具有開放的接口,可以與多種軟件連接實現聯合仿真優化。本文的多目標優化平臺的搭建采用專業的優化軟件iSIGHT通過接口技術集成BikeSim來實現,為了處理仿真計算過程中產生的大量數據,這里采用MATLAB編寫數據處理程序,并集成到iSIGHT中實現數據的自動化處理。在這三個軟件的集成中,BikeSim負責不同工況的仿真,MATLAB負責數據處理,iSIGHT的優化控件提供優化算法并驅動BikeSim和MATLAB實現迭代計算。每一次計算過程中,優化控件給出的一組懸架特性參數將同時傳遞給加速/制動工況和B級路面工況,兩種工況將同時開始仿真,形成并行計算任務。建立的并行任務優化工作流程如圖12,并行任務的計算結果將同時傳遞給優化控件。

圖12 并行任務優化工作流程Fig.12 The workflow of parallel optimized tasks

iSIGHT驅動BikeSim和MATLAB進行聯合仿真優化,能夠使各個軟件發揮專長使優化計算的精度和效率都極大提高。

2.3 參數設置與計算

本文以摩托車前后懸架的剛度和阻尼作為設計變量。由于摩托車的懸掛質量分配系數ε遠小于1,與轎車的差別較大,前后懸架系統耦合緊密,因此不宜采用偏頻法來確定前后懸架剛度的取值范圍。這里是根據懸架的靜撓度來確定剛度的上下限的,在BikeSim中通過在靜平衡時懸架變形量在一定范圍內來確定剛度的上下限,阻尼的上下限按照與剛度初始值擴大相同的比例來確定。這里將表1中懸架原始特性參數值作為設計變量初始值,初始值和上下限如表2。

表2 設計變量參數表Tab.2 The list of design variables

在平順性工況中,由于懸架動撓度的增大會使撞擊限位的概率增大,使平順性變差,因此在優化的同時必須對前后懸架的動撓度進行約束。一般要求為:

式中:fd1,fd2分別為前后懸的動撓度;[fd1]=50 mm,[fd2]=25 mm分別是前后懸的限位行程。

遺傳算法具有魯棒性、全局最優性、高效并行性、不要求函數連續可導等特點,近年來在車輛懸架方面的應用越來越多[12],然而基本遺傳算法只適用于傳統的單目標優化問題。對于本文的多目標優化問題,采用iSIGHT提供的NSGA-Ⅱ算法來實現,設定種群數量50、遺傳代數50、交叉率0.9。

經過2 500次計算,得到的計算結果如圖13所示,將計算結果向坐標平面投影可以得到圖14、圖15,從圖中可以看出計算結果已收斂,在兩個坐標面上Pareto前沿已經形成,由多目標優化的基礎理論可知Pareto前沿是各種權重下最優值點的集合,從計算結果可以發現兩個Pareto前沿分別近似為直線段。

3 優化計算結果分析

3.1 目標權重與最優值的確定

多目標優化的最優值是針對某一權重下的最優值,為了找出最優值首先要確定目標權重。但是先決策后尋優的模式無法準確把握目標權重,該模式在應用中存在明顯的缺陷,本文采用的是先尋優后決策的模式[13-14]。由多目標優化的基礎理論可知:多目標優化的最優值通常存在于Pareto前沿上。因此本文根據計算得到的Pareto前沿的特點來確定權重關系。從圖14、圖15可以看出兩個圖中的Pareto前沿近似為直線,通過數據擬合可以得到Pareto前沿的斜率,Pareto前沿的斜率反映了目標之間的權重關系。

設加速/制動工況下俯仰角最大值fyjd的權重為P;等速B級路面工況下俯仰振動加速度均方根值rmsp的權重為M,垂向振動加速度均方根值rmsv的權重為N。在rmsp-fyjd平面上,通過數據擬合可知Pareto前沿的斜率為 -0.39,為了計算方便這里取 -0.4。同樣,在rmsv-fyjd平面上得到的 Pareto前沿斜率為-0.83,這里取 -0.8,即權重關系為:

設X、Y、Z分別代表俯仰振動、垂向振動和俯仰角最大值,由于本文的優化目標是希望這三個值均能降低,即求函數T的最小值,其表達式為:

將iSIGHT中2500次優化計算的全部結果導入MATLAB中進行數據處理,第2195次計算結果為最優值,從圖14、圖15可以看出2 195次與初始點的對比。

圖13 計算散點圖Fig.13 Optimization 3D scatter plot

圖14 rmsp-fyjd平面投影Fig.14 The rmsp-fyjd plane projection

將圖13的計算散點圖向兩個坐標平面投影形成的Pareto前沿是兩條近似直線,實際上在空間坐標系中Pareto前沿應該是一個前沿面,在MATLAB中的數據處理即是找出平面T=0.4X+0.8Y+Z與Pareto前沿面的切點,前沿面上的切點即為最優值點。

3.2 優化前后懸架性能比較

優化前后懸架特性參數如表3所示;優化前后的動撓度均方根值對比如表4所示,優化后的前后懸動撓度略有增大;優化前后的制動俯仰角最大值fyjd、俯仰振動與垂向振動加速度均方根值rmsp、rmsv如表5所示,從表中可以看出,加速制動工況下的俯仰角最大值優化效果最為明顯,降低幅度高達57.7%,同時平順性工況下的俯仰和垂向振動均有不同程度下降,完全實現了預期的優化目標。

表3 優化前后懸架特性參數對比Tab.3 Suspension parameters before and after optimization

圖15 rmsv-fyjd平面投影Fig.15 The rmsv-fyjd plane projection

表4優化前后動撓度均方根值對比Tab.4 Suspension dynamic deflection RMS values before and after optimization

表5 優化前后懸架性能對比Tab.5 Comparison of optimization objectives

圖16 優化前后的俯仰角變化Fig.16 Pitch angle before and after optimization

圖17 優化前后前懸架變形Fig.17 Front suspension deformation before and after optimization

圖18 優化前后的后懸架變形Fig.18 Rear suspension deformation before and after optimization

圖19 優化前后俯仰振動對比Fig.19 Comparison of pitching vibration

圖20 優化前后車身質心垂向振動對比Fig.20 Comparison of vertical vibration

圖16、17是優化以后加速/制動工況下的車身俯仰角的變化以及前懸架變形關系圖,從圖中可以看出俯仰角和懸架變形量明顯降低了,優化后前懸架變形最大值降低為88 mm,叉管行程為105 mm,因此前懸架撞擊限位的現象消除,騎乘舒適性得到了提高。圖18是優化前后的后懸架變形圖,從圖中可以看出優化后的后懸架變形比優化前有相應的增大,優化后的懸架變形在前5 s的加速階段和勻速階段無明顯突變,最大值41 mm出現在制動停止后的回彈階段,小于后懸架的限位行程50 mm,因此在加速階段和停止后的回彈階段都不會引起后懸架的撞擊。圖19、20是優化前后俯仰振動和垂向振動的功率譜對比,從圖中可以看出采用匹配后的懸架參數整車的俯仰和垂向振動均有降低,在整個頻率范圍內都沒有出現惡化的現象,并且在低頻范圍內垂向振動降低的效果更加明顯,優化后的整車性能滿足了預期的要求。

通過以上優化前后兩種工況下的懸架各項性能的對比分析,加速/制動工況的俯仰角峰值明顯降低,騎乘舒適性得到改善,隨機路面工況下的平順性得到了提高,驗證了多工況并行任務優化方法的有效性。

4 結論

(1)采用摩托車動力學專用軟件BikeSim建立了路面和人-車系統模型,進行了加速制動工況和平順性工況的仿真;以降低加速/制動工況下的俯仰角最大值、等速B級路面工況下的俯仰和垂向振動為目標,采用iSIGHT集成BikeSim和MATLAB的聯合優化仿真的方法建立了兩種工況的并行優化任務,對摩托車懸架系統進行了多目標優化,得到了Pareto前沿,實現了預期的優化目標。

(2)采用先優化后決策的模式,根據優化計算結果對Pareto前沿進行了分析,提出了一種確定目標權重的方法,找出了該權重下的最優值;通過對優化前后的俯仰角、前懸架變形、俯仰和垂向振動加速度均方根值、功率譜密度的結果分析表明,消除了制動時的前懸架撞擊現象,騎乘舒適性得到了提高,同時隨機路面平順性得到了改善,驗證了該方法的有效性,為車輛綜合性能的優化提供了重要參考。

[1]李德寬.汽車工程手冊摩托車篇[M].北京:人民交通出版社,2001.

[2]Cossalter V,Doria A,Pegoraro R,et al.On the non-linear behavior of motorcycle shock absorbers[J].Journal of Automobile Engineering,2010,224(1):15 -27.

[3]Audenino A L,Belingardi G.Modelling the dynamic behavior ofa motorcycle damper[J]. JournalofAutomobile Engineering,1995,209(4):249 -262.

[4]徐中明,盧少波,文 俐,等.基于MATLAB的摩托車舒適性仿真分析[J].計算機仿真,2006,23(7):249-252.XU Zhong-ming,LU Shao-bo,WEN Li,et al.Simulation of ride comfort for motorcycle based on MATLAB[J].Computer Simulation,2006,23(7):249 -252.

[5]陳 松,蹇開林,黃澤好,等.摩托車平順性的仿真研究[J].機械強度,2006,28(3):465 -469.

CHEN Song, JIAN Kai-Lin, HUANG Ze-hao, etal.Simulation and research of motorcycle's rider comfort[J].Journal of Mechanical Strength,2006,28(3):465 -469.

[6]Spelta C,Savaresi S M,Fabbri L.Experiment analysis of a motorcycle semi-active rear suspension[J]. Control Engineering Pratice,2010,18(11):1239 -1250.

[7] Wu L,Zhang W J.Hierarchical modeling of semi-active control of a full motorcycle suspension with six degrees of freedoms[J]. International Journal of Automotive Technology,2010,11(1):27 -32.

[8]吳 龍,陳花玲,陳玲莉.摩托車半主動懸架分層預測控制及仿真[J].系統仿真學報,2006,18(8):2239 -2243.

WU Long,CHEN Hua-ling,CHEN Ling-li.Hierarchical preview control and simulation of semi-cative motorcycle suspension[J].Journal of System Simulation,2006,18(8):2239-2243.

[9]于學華,丁康.汽車俯仰與跳動復合運動的分析[J].噪聲與振動控制,2005,4:26 -29.

YU Xue-hua,DING Kang.Analysis of the pitching and bouncing motion of the car[J].Noise and Vibration Control,2005,4:26 -29.

[10]張永林.用諧波疊加法重構隨機道路不平順高程的時域模型[J].農業工程學報,2003,19(6):32 -34.

ZHANG Yong-lin.Time domain model of road irregularities simulated using the harmony superposition method[J].Transactions of The Chinese Society of Agricultural Engineering,2003,19(6):32 -34.

[11]張志飛,徐中明,賀巖松.基于多體動力學的摩托車平順性分析[J].中國機械工程,2010,21(8):988 -992.

ZHANG Zhi-fei,XU Zhong-ming,HE Yan-song.Analysis of motorcycle riding comfort based on multi-body dynamics[J].China Mechanical Engineering,2010,21(8):988 -992.

[12]奉銅明,鐘志華,閆曉磊,等.基于NSGA-Ⅱ算法的多連桿懸架多目標優化[J].汽車工程,2010,32(12):1063-1066.

FENG Tong-ming,ZHONG Zhi-hua,YAN Xiao-lei,et al.Multi-objective optimization for multi-link suspension based on NSGA-Algorithm[J].Automotive Engineering,2010,32(12):1063-1066.

[13]王 濤.汽車懸架參數的多目標多標準決策優化[J].農業機械學報,2009,40(4):27 -32.

WANG Tao.Multi-objective and multi-criteriadecision optimization ofautomobile suspension parameters[J].Transactions of the Chinese Society for Agricultural Machinery,2009,40(4):27 -32.

[14]王 濤,陶 薇.考慮隨機因素的汽車懸架參數多目標穩健優化[J].振動與沖擊,2009,28(11):146 -149.

WANG Tao,TAO Wei.Multi-objective robust optimization of automobile suspension parameters considering random factors[J].Journal of Vibration and Shock,2009,28(11):146-149.

猜你喜歡
振動優化模型
一半模型
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
重要模型『一線三等角』
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
重尾非線性自回歸模型自加權M-估計的漸近分布
中立型Emden-Fowler微分方程的振動性
主站蜘蛛池模板: 国产精品一老牛影视频| 亚洲成人高清无码| 在线欧美日韩国产| 黄色在线不卡| 狠狠干综合| 亚洲欧美日本国产综合在线| 97国产一区二区精品久久呦| 国产精品吹潮在线观看中文| 在线观看亚洲精品福利片| 国产精品综合久久久| 激情视频综合网| 狠狠干欧美| 亚洲色精品国产一区二区三区| 国产精品观看视频免费完整版| 国产成人精品高清不卡在线| 在线国产毛片| 久久男人资源站| 国产成人亚洲精品无码电影| 91国内外精品自在线播放| 国产一区二区三区精品欧美日韩| 一级爱做片免费观看久久| 99er这里只有精品| 精品国产美女福到在线不卡f| 区国产精品搜索视频| 国产精品手机在线观看你懂的| 亚洲男人天堂2020| 激情六月丁香婷婷| 国产女同自拍视频| 在线五月婷婷| 2020极品精品国产| 久久精品日日躁夜夜躁欧美| 欧洲日本亚洲中文字幕| 国产尤物在线播放| 伊人91视频| 国产欧美日韩综合在线第一| 在线观看国产精美视频| 久久亚洲国产一区二区| 在线观看亚洲精品福利片| 风韵丰满熟妇啪啪区老熟熟女| 亚洲精品福利视频| 国产美女91视频| 毛片视频网址| 国产日韩精品欧美一区灰| 国产精品无码影视久久久久久久 | 99国产精品一区二区| 亚洲人成在线免费观看| 国产精品开放后亚洲| 一级毛片网| 亚洲黄色网站视频| 伊人久久大香线蕉aⅴ色| 亚洲妓女综合网995久久 | 免费jizz在线播放| 色偷偷男人的天堂亚洲av| 亚洲av无码久久无遮挡| 欧美日韩中文字幕在线| 久久中文字幕2021精品| 日韩区欧美区| 99偷拍视频精品一区二区| 综合色区亚洲熟妇在线| 色成人亚洲| 亚洲成a人片在线观看88| 欧美视频免费一区二区三区| 欧美日韩激情| 99视频精品全国免费品| 日韩精品无码免费一区二区三区| 国产成人啪视频一区二区三区| 国产在线一二三区| 91精品情国产情侣高潮对白蜜| 一区二区三区四区日韩| 国产麻豆va精品视频| 欧美另类一区| 亚洲IV视频免费在线光看| 91国内在线观看| 98精品全国免费观看视频| 91在线精品免费免费播放| 亚洲中文无码av永久伊人| 一级毛片网| 国产亚洲第一页| AV片亚洲国产男人的天堂| 中国黄色一级视频| 日本成人不卡视频| 国产高清国内精品福利|