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

基于NSGA-II遺傳算法的磁流變懸置多目標優化

2014-02-05 03:50:00段緒偉李以農鄧召學
振動與沖擊 2014年3期
關鍵詞:優化

段緒偉,李以農,鄭 玲,鄧召學

(重慶大學 機械傳動國家重點實驗室,重慶 400044)

基于NSGA-II遺傳算法的磁流變懸置多目標優化

段緒偉,李以農,鄭 玲,鄧召學

(重慶大學 機械傳動國家重點實驗室,重慶 400044)

磁流變懸置集總參數優化是設計高性能發動機懸置的關鍵。為克服以往懸置優化中優化目標單一、優化目標選取不合理、未考慮實際加工可行性等問題,建立單自由度磁流變懸置隔振系統數學模型,提出倍程區間靈敏度分析法,對各集總參數靈敏度進行分析,并以此為依據選取優化變量。以發動機常用轉速激振頻率段的力傳遞率積分為優化目標,采用改進型非支配排序遺傳算法(NSGA-II)進行多目標優化。在一定范圍內將結構尺寸進行離散化處理,計算各組離散尺寸對應的集總參數值,以離散集總參數與集總參數Pareto非劣解之間的綜合距離為準則篩選最優解。

磁流變懸置;集總參數;優化;倍程區間靈敏度;NSGA-II算法

磁流變液是近年來研究較多的新型智能材料,它通常是由非磁性載體液、磁性微粒和穩定劑組成的兩相懸浮液體,其基本特征:在外加磁場作用下,能在瞬間從自由流動的液體轉變為半固體,呈現可控的屈服強度,而且這種變化是可逆的[1]。利用磁流變液的這種流變特性,Ciocanel等[2]設計了混合模式的磁流變懸置,并對其低頻動特性進行了理論和實驗研究;Choi等[3]設計了混合模式磁流變懸置,對其在結構振動控制中的應用進行了相關研究;Farjoud等[4]設計了擠壓模式磁流變懸置,對其磁路等進行了研究;王雪婧[5]、魏付俊[6]、李銳[7]對不同形式磁流變懸置的動特性進行了研究。

在懸置優化設計方面,Ahn等[8]以普通液壓懸置隔振力傳遞率各極值點及其發生頻率的加權組合為優化目標,采用EALA法進行單目標優化,但優化效果受權重系數影響大,且幾個頻率點上傳遞率的疊加不能很好地反映隔振系統的隔振效果。Li等[9]以普通液壓懸置動剛度及阻尼滯后角在特征頻率處取值的加權組合為優化目標,在獲得其實驗數據的基礎上,采用基因神經網絡算法對懸置動特性進行單目標優化,其神經網絡訓練的好壞依賴于實驗數據準確性,對于懸置初期設計不適用。且目前很少有針對磁流變懸置進行多目標優化的研究,而優化對于獲得具有良好隔振性能的磁流變懸置尤為重要,因此對磁流變懸置集總參數進行優化具有重要的工程意義。

1 磁流變懸置隔振

1.1 磁流變懸置結構及工作原理

本文設計了新型流動模式磁流變懸置,其結構示意圖如圖1所示。主要由橡膠主簧8、擾流盤9、磁芯內上1、磁芯外上3、磁芯內下6、磁芯外下5、電磁激勵線圈11、可控阻尼通道2、解耦膜4、上液室10、下液室7以及橡膠底膜12等組成。

低頻大振幅激振時,橡膠主簧發生變形,上液室壓力變化,使磁流變液流經可控阻尼通道。利用磁流變液的流變特性,調節激勵線圈電流大小控制可控阻尼通道有效磁極對應處的磁場強度,即可調節其屈服強度,實現可控阻尼通道液阻調控,從而調控磁流變懸置低頻動態特性。浮動式解耦膜與解耦器隔板之間間隙小,此時因激振頻率低、激振幅值大,解耦膜“阻斷”解耦器,使流經解耦器的流量極少[10]。高頻小振幅激振時,在擾流盤擾流作用下,與擾流盤相鄰的區域紊流度增加,使得解耦膜在其平衡位置附近振動,部分磁流變液流經解耦器,緩解上下液室壓力波動,有效抑制懸置高頻動態硬化,可改善磁流變懸置高頻動態特性。

圖1 磁流變懸置結構示意圖Fig.1 Structural sketch of MR mount

圖2 單自由度隔振力學模型Fig.2 Mechanical model of vibration isolation

1.2 磁流變懸置集總參數模型

為研究磁流變懸置隔振性能,建立含磁流變懸置集總參數的單自由度隔振系統力學模型,如圖2所示。

Me為單個懸置支撐的等效發動機質量,Kr、Br分別為橡膠主簧的靜剛度和阻尼,C1、C2分別為上、下液室的體積柔度。Ap、Ad分別為橡膠主簧等效活塞面積、解耦膜面積,Im、Id分別為可控阻尼通道的液感、解耦器液感,Rm(I)、Rd分別表示可控阻尼通道的液阻、解耦器液阻。通過可控阻尼通道、解耦器的流量分別為Qm(t)、Qd(t)。上、下液室的壓力分別為P1(t)、P2(t)。Xe(t)為發動機振動位移激勵,F(t)為傳遞給基體的力。

低頻大振幅激振時,通過解耦器的磁流變液極少,高頻小振幅激振時,通過可控阻尼通道的磁流變液也極少,兩者均可忽略。利用流體力學理論等,可推導出懸置低頻復剛度式(1)及力傳遞率式(2),懸置高頻復剛度式(3)及力傳遞率式(4)。

式中:Rm(I)=Rm0+RIm,Rm0為可控阻尼通道零場粘性液阻=6lcζ(NI/2g)τ/為屈服液阻,lc為有效磁極長度,N、I分別為激勵線圈的匝數和電流,g為通道間隙,Qm為通道流量,ζ、τ為磁流變液屬性參數[5],KL1、KH1、KL2、KH2分別為低、高頻復剛度的實部和虛部,TLF、THF分別為低、高頻懸置力傳遞率。

2 隔振優化問題

為有效隔離發動機的振動向車架傳遞,磁流變懸置隔振系統力傳遞率峰值、怠速工況激振以及常用高轉速激振時力傳遞率都應盡可能小,力傳遞率峰值發生頻率應避開發動機怠速工況,同時各集總參數的取值范圍受實際結構等因素制約。因此,本文以低頻段(f1,f2)=(6 Hz,30 Hz)和高頻段(f3,f4)=(100 Hz,200 Hz)的力傳遞率積分值作為目標,以峰值、峰值發生頻率和集總參數的取值為約束條件進行優化,則優化問題可描述如下:

其中:TL(f,X)、TH(f,X)分別為低、高頻傳遞率,TL_max為低頻力傳遞率最大值,fL_max為低頻力傳遞率最大值發生頻率,X為含各集總參數的優化變量,X0為各集總參數初始值,K1、K2為集總參數取值范圍的系數矩陣。

3 選取優化變量

磁流變懸置集總參數的調節,可通過調整結構尺寸、橡膠硬度、激勵線圈匝數、激勵電流大小等來實現,但受尺寸、加工可行性、密封要求等限制,以最小的調整獲得最佳的參數匹配尤為重要,為此通過集總參數的靈敏度分析各參數對隔振性能的影響程度,并以此為依據選取合適的優化變量。

3.1 倍程區間靈敏度法

集總參數的物理意義各異,參數值數量級相差大,參數變化范圍大,且具有較強的非線性。利用差分法、攝動法等常規靈敏度分析法計算的靈敏度值可比性差、可信度低。針對集總參數的特點,采用倍程區間靈敏度分析法對集總參數的靈敏度進行分析,此法是對區間靈敏度法[12]的改進,即將集總參數攝動范圍進行倍程化處理,利用區間數學計算其相對靈敏度。該法對參數變化量Δx的選擇比差分法要求低,且求解非單調非線性問題靈敏度的結果比差分法更合理。倍程區間靈敏度分析法簡述如下:

3.2 集總參數靈敏度分析

為分析磁流變懸置的各集總參數對隔振性能的影響,可等價轉換為各集總參數對優化目標F1和F2的倍程區間靈敏度,由此便可選出合適的優化變量。

針對某型號發動機,參考原液壓懸置設計流動模式磁流變懸置,其集總參數初始值如表1所示。根據集總參數初始值,采用倍程區間靈敏度法,利用MATLAB編程可計算各集總參數對優化目標F1和F2的倍程區間靈敏度值。為具有可比性,各參數的區間倍程化系數均為β=0.1,結果如表2和表3所示。

表1 磁流變懸置集總參數初始值Tab.1 Lumped parameters of original MR mount

表2 集總參數對F1的倍程區間靈敏度值Tab.2 Multiple interval sensitivity of lumped parameters

表3 集總參數對F2的倍程區間靈敏度值Tab.3 Multiple interval sensitivity of lumped parameters

從計算的倍程區間靈敏度值可知,各集總參數對目標F1的靈敏度絕對值由大到小為Ap>C1>Rm>Kr>Im>C2>Br>Rd=Id=0,各集總參數對目標 F2的靈敏度絕對值由大到小為Kr>C1>Rd>Ap>Id>C2>Br>Rm=Im=0。參數C2、Br的靈敏度絕對值很小,即表明其對隔振性能的影響很小,可不予考慮。

4 集總參數優化

為實現集總參數的多目標約束優化,本文借鑒具有較高計算性能和較少人工經驗參與的非支配排序遺傳算法NSGA-II,并針對磁流變懸置的集總參數取值等特點,改進NSGA-II遺傳算法,使之適用于該問題。

4.1 NSGA-II遺傳算法及其改進

非支配排序遺傳算法(NSGA-II)是Deb對Srinivas所提出的NSGA算法的改進[14],相對于NSGA算法具有以下優越性:① 計算復雜度由O(mN3)降低為O(mN2)。其中,m為目標函數個數;N為種群規模。② 引入精英策略,可保留父代精英個體,且無需人為確定共享半徑參數。③ 采用擁擠度比較算子使種群具有多樣性。

采用非支配排序遺傳算法(NSGA-II)進行優化時,由于磁流變懸置集總參數的數量級相差大,若采用二進制編碼,染色體長度驟增,計算量大;若采用實數編碼,計算精度受到限制,且采用標準NSGA-II遺傳算法進行優化時未考慮約束條件。因此,本文對標準的非支配多目標遺傳算法進行以下改進:① 在進行遺傳算子之前對變量X進行同級化處理XR=RX,其中R=diag(r1,r2,r3...r7)=diag(103,10-5,1011,10-5,10-8,10-4,10-7)為同級化轉換對角矩陣,ri為第 i個參數同級化因子。在完成遺傳算子操作后,還原參數值進行適應度計算。②在選擇個體時,增加判斷個體是否滿足約束條件的環節,若滿足則保留,反之則直接淘汰。改進后的算法流程如圖3所示。

圖3 改進型NSGA-II算法流程圖Fig.3 Process of improved NSGA-II algorithm

4.2 優化計算

在磁流變懸置的集總參數中,Ap、Kr、C1對磁流變懸置隔振性能影響較大,且通過調節結構尺寸可使其在較大范圍內變化,其倍程系數β可設置較大。可控阻尼通道的液阻Rm與液感Im、解耦器的液阻Rd與液感Id分別為同一結構的不同集總參數,為兩對非獨立參數。

此外,因可控阻尼通道液阻隨激勵電流的增大而增大,而液感基本不變,且Lr(Rm)≈Lr(Im),為實現兩者的匹配,則倍程系數β(Rm)應略大于β(Im);解耦器液阻Rd與液感Id不受激勵電流影響,Lr(Rd)>Lr(Id),為實現兩者的匹配,則倍程系數 β(Rd)應略大于β(Id)。因此,集總參數取值范圍的系數對角化矩陣可設置為:

K1=diag(0.7,0.7,0.7,0.7,0.65,0.75,0.7)

K2=diag(1.3,1.3,1.3,1.3,1.35,1.25,1.3)

根據磁流變懸置初始結構的集總參數,利用MATLAB編寫改進型NSGA-II算法進行優化計算,在頻域區間內求解力傳遞率積分時,采用4節點Newton-Cotes公式進行數值積分,設置種群大小為100,遺傳代數為300代,優化后所得 Pareto非劣解集前沿面如圖4所示。

圖4 Pareto非劣解集前沿面Fig.4 The Pareto solution frontier

4.3 篩選Pareto非劣解集中的最優個體

考慮結構尺寸限制、加工精度、結構布局合理性以及密封要求等實際問題,在篩選經優化所得的非劣解集中最優個體時,由各集總參數關于結構尺寸以及流體屬性的計算公式[5,15],根據實際結構可取的離散尺寸值,可獲得各集總參數與實際結構尺寸相對應的取值離散集,然后計算所有Pareto非劣解個體的個體綜合距離DR,以個體綜合距離DR最小值個體作為最優個體。考慮到各集總參數的倍程區間靈敏度以及各集總參數量級不一致等因素,定義同級化個體綜合距離DR如下:

其中:Y表示獨立集總參數 Ap、Kr、C1,rY表示獨立集總參數Y的同級化因子。Yj為離散集中與第j組尺寸對應的集總參數值,YPi為Pareto非劣解集中第i個個體對應的集總參數值,j=1,2,3…n,i=1,2,3…k,n 為離散尺寸組數。

圖5 液阻Rm液感Im組離散集和Pareto集Fig.5 Fluid Rm|Imin discretization and Pareto set

圖6 Pareto集個體綜合距離Fig.6 Synthesized distance of Pareto set’s individuals

根據個體綜合距離公式(7),可計算改進型NSGA-II算法多目標優化所得的Pareto非劣解集中每個個體的個體綜合距離,選取個體綜合距離最小個體作為最符合實際結構要求的最優個體。Pareto非劣解集中每個個體的個體綜合距離DiR如圖6所示,與最優個體對應的同級化集總參數值如圖7所示。

圖7 最優個體同級化值Fig.7 Peer parameter value of the best individual

4.4 最優隔振性能及懸置最優動特性

根據綜合距離準則選取的最優集總參數,由式(2)編寫MATLAB程序可得磁流變懸置隔振系統低頻最優力傳遞率,如圖8所示。

圖8 低頻力傳遞率對比Fig.8 Contrast the force transmissibility in low frequency range

由圖8可知,優化后低頻大振幅激振時隔振總體效果優于原結構。優化后的力傳遞率峰值TLmax由3.07降為2.564,峰值發生頻率fLmax由 8.8 Hz 降為 7.2 Hz。優化后的力傳遞率在0~7.36 Hz內略大于原結構,而在7.36~50 Hz內明顯小于原結構。

最符合實際結構要求的最優集總參數組中,可控阻尼通道的液阻Rm0是零磁場粘性液阻,而Rm(I)=Rm0+RIm隨激勵電流增大而呈非線性增大。因此,由式(1)和(2)進行MATLAB編程可得磁流變懸置集總參數最優化后的低頻力傳遞率曲面以及此時磁流變懸置的動態特性曲面,分別如圖9、圖10和圖11所示。

從圖9低頻傳遞率曲面中可知,隨著激勵電流的增大,力傳遞率峰值呈先減小后略增大變化趨勢,峰值發生頻率基本不變。在實際運用中,為獲得低頻段最優的隔振效果,可由轉速傳感器采集發動機實時轉速,計算其二階激振頻率,以該激振頻率下傳遞率最小的激勵電流作為最優控制電流Iopt(f),從而實現磁流變懸置低頻大振幅激振時的最優隔振控制。

從圖10和圖11可知,在最優激勵電流Iopt(f)下可得的磁流變懸置最優動態特性,這為設計具有優良隔振性能的磁流變懸置提供了動特性最優目標,具有一定的參考價值。

諸多研究表明,發動機懸置隔振系統對懸置的理想特性要求是:低頻大振幅激振時懸置應具有高剛度大阻尼特性,高頻小振幅激振時應具有低剛度小阻尼特性[16]。因此,在高頻小振幅激振時,磁流變懸置不施加激勵電流,而通過擾流盤的擾流作用使更多的磁流變液流經解耦器,緩解高頻動態硬化,從而獲得較小的動剛度和阻尼滯后角,便可使磁流變懸置隔振系統具有更理想的隔振效果。由優化后所得的解耦器最優液阻Rd與液感Id參數,根據式(4)編程可得磁流變懸置隔振系統高頻最優力傳遞率曲線,如圖12所示。

圖9 低頻力傳遞率曲面Fig.9 Curved surface of force transmissibility in low frequency range

圖10 低頻懸置動剛度曲面Fig.10 Curved surface of dynamic stiffness in low frequency range

圖11 低頻懸置阻尼滯后角曲面Fig.11 Curved surface of loss angle in low frequency range

圖12 高頻力傳遞率對比Fig.12 Contrast the force transmissibility in high frequency range

由圖12可知,優化后高頻激振時隔振效果優于原結構,優化后力傳遞率峰值THmax由0.030 1降為0.010 3,峰值發生頻率fHmax由107.2 Hz增為128.8 Hz。優化后的力傳遞率在50~75 Hz內略小于原結構,在75~139 Hz內明顯小于原結構,在139~250 Hz內略小于原結構。

5 結論

本文建立了磁流變懸置隔振系統的等效集總參數模型,推導了磁流變懸置復剛度和力傳遞率的解析式,提出倍程區間靈敏度分析法,通過集總參數靈敏度分析,確定了磁流變懸置隔振系統的優化變量,以常用轉速激振頻率力傳遞率積分為目標函數,采用改進型非支配排序遺傳算法(NSGA-II),對磁流變懸置隔振系統參數進行了多目標優化,提出了以個體綜合距離來衡量Pareto非劣解集中個體優越性的方法。研究結果表明:本文提出的倍程區間靈敏度分析法和多目標優化方法是正確、可行的,優化后的懸置隔振系統,無論在低頻段還是高頻段,隔振效果均明顯優于原系統,為磁流變懸置隔振系統的優化設計提供了科學方法與手段。在實際中,應提高加工和裝配的精度、控制橡膠主簧硬度、選用性能穩定的磁流變液等,從而使懸置的各集總參數與最優值盡可能一致。

[1] Olabi A G,Grunwald A.Design and application of magnetorheological fluid[J].Materials and Design,2007,28:2658-2664.

[2] Ciocanel C,Nguyen,Elahinia M.Design and modeling of a mixed mode magnetorheological(MR)fluid mount[C].Active and Passive Smart Structures and Integrated Systems,2008,6928:69281C1-69281C10.

[3]Choi S B,Hong S R,Sung K G,et al.Optimal control of structural vibrations using a mixed-mode magnetorheological fluid mount [J]. InternationalJournalofMechanical Sciences,2008,50(3):559-568.

[4]Farjoud A, CraftM, BurkeW, etal. Experimental investigation ofMR squeeze mounts [J]. Journalof Intelligent Material Systems and Structures,2011,22:1645-1652.

[5]王雪婧.磁流變半主動發動機懸置隔振性能與控制方法研究[D],長春:吉林大學,2011.

[6]魏付俊.汽車動力總成磁流變懸置的設計和半主動控制研究[D],南京:南京航空航天大學,2007.

[7]李 銳.基于磁流變液懸置的發動機隔振方法與試驗研究[D],重慶:重慶大學,2009.

[8] Ahn Y K,Song J D,Yang B S.Optimal design of engine mount using an artificial life algorithm[J].Journal of Sound and Vibration,2003,261:309 328.

[9]Li Q,Zhao J C,Zhao B,et al.Parameter optimization of a hydraulic engine mount based on a genetic neural network[J].Proceedings of the Institution of Mechanical Engineers,Part D:Journal of Automobile Engineering,2009,223(9):1109-1117.

[10] Singh R, Kim G, Ravindra P V.Linearanalysisof automotive hydro-mechanical mount with emphasis on decoupler characteristics [J]. JournalofSound and Vibration,1992,158(2):219-243.

[11] Ahn Y K,Song J D,Yang B S,et al.Optimal design of nonlinear hydraulic engine mount[J].Journal of Mechanical Science and Technology,2005,19(3):768-777.

[12]邱志平,王曉軍.結構靈敏度分析的區間方法[J].兵工學報,2005,26(6):798 -802.

QIU Zhi-ping,WANG Xiao-jun.An interval method for sensitivity analysis of the structures[J].Acta Armamentarii,2005,26(6):798-802.

[13]沈祖和.區間分析方法及其應用[J].應用數學與計算數學,1983,2:1 -29.

SHEN Zhu-he. The method ofintervalanalysis and application [J].Communication on Applied Mathematics and Computation,1983,2:1-29.

[14]Kalyanmoy D,Amrit P,Sameer A,et al.A fast and elitist multi-objective genetic algorithm:NSGA - II [J].IEEE Transactions on Evolutionary Computation,2002,6(2):182-197.

[15] Nguyen,Ciocanel C,Elahinia M.Analytical modeling and experimental validation of a magnetorheological mount[C].Active and Passive Smart Structures and Integrated Systems,2009,7288:72881D1-72881D7.

[16] Kim G.Study of passive and adaptive hydraulic engine mount[D].Columbus:The Ohio State University,1992.

Multi-objective optimization for a MR engine mount based on NSGA-II algorithm

DUAN Xu-wei,LI Yi-nong,ZHENG Ling,DENG Zhao-xue

(State Key Laboratory of Mechanical Transmission,Chongqing University,Chongqing 400044,China)

To achieve a high performance,the design optimization of lumped parameters for a magneto-rheological(MR)engine mount is essential.Mathematical model of a single DOF vibration isolation system was established and the multiple interval sensitivity method was proposed to overcome drawbacks of conventional optimization designs,such as,single objective optimization,improper optimization objective,unfeasible machining and so on.Optimization variables in a MR engine mount were selected with a lumped parameter multiple interval sensitivity analysis.The integral of force transmissibility within normal frequency ranges of an engine was assigned as an objective function,the non-dominated sorting genetic algorithm(NSGA-II)was improved and used to optimize design variables.The synthesized distances between Pareto lumped parameters and discontinuous lumped parameters matching along with physical discretization dimensions were calculated to select the most appropriate solution from Pareto lumped parameters.

MR mount;lumped parameters;optimization;multiple interval sensitivity;NSGA-II algorithm

U463.33+

A

10.13465/j.cnki.jvs.2014.03.035

機械傳動國家重點實驗室項目(0301002109165)資助

2013-01-04 修改稿收到日期:2013-03-11

段緒偉 男,碩士生,1987年9月生

鄭 玲 女,博士,教授,博士生導師,1963年生

猜你喜歡
優化
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
PEMFC流道的多目標優化
能源工程(2022年1期)2022-03-29 01:06:28
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
圍繞“地、業、人”優化產業扶貧
今日農業(2020年16期)2020-12-14 15:04:59
事業單位中固定資產會計處理的優化
消費導刊(2018年8期)2018-05-25 13:20:08
4K HDR性能大幅度優化 JVC DLA-X8 18 BC
幾種常見的負載均衡算法的優化
電子制作(2017年20期)2017-04-26 06:57:45
主站蜘蛛池模板: 国产成人高精品免费视频| 五月婷婷中文字幕| 欧美国产日产一区二区| 国产精品美乳| 91精品aⅴ无码中文字字幕蜜桃 | 免费又黄又爽又猛大片午夜| 亚洲免费人成影院| 午夜一区二区三区| 成人日韩欧美| 日韩高清在线观看不卡一区二区| 中文字幕无码电影| 国产综合日韩另类一区二区| 精品色综合| 成人免费网站久久久| 亚洲色图另类| 亚洲无码免费黄色网址| 国产精品福利导航| 国产精品私拍在线爆乳| 欧美成人综合在线| 国产爽爽视频| 日韩a级片视频| 亚洲欧美另类日本| 在线视频97| 欧美成人二区| 国产欧美性爱网| 欧洲一区二区三区无码| 久久人体视频| 午夜国产精品视频黄 | 第一区免费在线观看| 亚洲久悠悠色悠在线播放| 免费一级毛片不卡在线播放| 国产真实乱了在线播放| 亚洲欧州色色免费AV| 操国产美女| 99re免费视频| 国产乱人视频免费观看| 国产精品美乳| 国产精品伦视频观看免费| 亚洲区欧美区| 亚洲狼网站狼狼鲁亚洲下载| 最新精品国偷自产在线| 狠狠干综合| 亚洲成人福利网站| 无码视频国产精品一区二区 | 亚洲自拍另类| 精品欧美日韩国产日漫一区不卡| 亚洲国产亚综合在线区| 国产第一页屁屁影院| 免费观看国产小粉嫩喷水| 久久96热在精品国产高清| 国产精品美女网站| 午夜天堂视频| 亚洲欧洲国产成人综合不卡| 亚洲高清在线播放| 欧美无专区| 欧美色伊人| 亚洲一区第一页| 亚洲午夜18| 国产微拍精品| AV不卡无码免费一区二区三区| 亚洲 欧美 偷自乱 图片| 国产成人免费手机在线观看视频| 日本成人不卡视频| a毛片在线播放| 欧美激情一区二区三区成人| 国产成人免费观看在线视频| 思思99热精品在线| 欧美一级片在线| 亚洲成aⅴ人片在线影院八| 丁香婷婷激情综合激情| 福利小视频在线播放| 国产二级毛片| 怡春院欧美一区二区三区免费| 国产成人区在线观看视频| 一本色道久久88综合日韩精品| 国产美女一级毛片| 亚洲一区无码在线| 成人永久免费A∨一级在线播放| 国内精品小视频福利网址| 亚洲一区无码在线| 国产91色在线| 日本a级免费|