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

基于缸蓋振動信號的柴油機多工況氣缸壓力識別方法研究*

2015-04-12 07:11:45劉建敏李曉磊喬新勇
汽車工程 2015年8期
關(guān)鍵詞:振動信號模型

劉建敏,李曉磊,2,喬新勇,張 杰

(1.裝甲兵工程學院機械工程系,北京 100072; 2.解放軍77136部隊,璧山 402760;3.中國人民解放軍駐318廠軍事代表室,北京 100053)

?

2015151

基于缸蓋振動信號的柴油機多工況氣缸壓力識別方法研究*

劉建敏1,李曉磊1,2,喬新勇1,張 杰3

(1.裝甲兵工程學院機械工程系,北京 100072; 2.解放軍77136部隊,璧山 402760;3.中國人民解放軍駐318廠軍事代表室,北京 100053)

對某12150型柴油機進行了缸內(nèi)燃燒激勵的瞬態(tài)動力學計算,分析了其缸蓋振動的位移、速度和加速度與缸內(nèi)燃燒特征參數(shù)的對應(yīng)關(guān)系。接著在此基礎(chǔ)上,對實測振動加速度進行數(shù)字積分和平均濾波得到振動位移信號,并利用希爾伯特包絡(luò)和滑動平均法提取了振動位移的趨勢項。再以該趨勢項為輸入?yún)?shù)構(gòu)建了Adaboost_BP集成神經(jīng)網(wǎng)絡(luò)模型,最后利用此模型對不同工況下的缸內(nèi)壓力進行識別。結(jié)果表明:振動位移趨勢項與缸內(nèi)壓力的良好對應(yīng)關(guān)系和參數(shù)本身的簡潔性有效降低了神經(jīng)網(wǎng)絡(luò)輸入的復(fù)雜度,提高了神經(jīng)網(wǎng)絡(luò)的訓練效率;集成神經(jīng)網(wǎng)絡(luò)模型能夠準確識別不同工況下的缸內(nèi)壓力,其泛化性和精度均有大幅度提高。

柴油機;振動;氣缸壓力;希爾伯特包絡(luò);集成神經(jīng)網(wǎng)絡(luò)

前言

缸內(nèi)壓力直接反映柴油機燃燒過程,是表征柴油機技術(shù)狀況的重要參數(shù)之一。對于臺架試驗和某些設(shè)計完善的車輛通常采用缸壓傳感器直接測量缸內(nèi)壓力,可實現(xiàn)發(fā)動機工作狀況的實時監(jiān)測。軍用車輛沒有預(yù)設(shè)缸內(nèi)壓力傳感器,同時由于空間的限制,實車安裝缸內(nèi)壓力傳感器非常繁瑣,嚴重影響了檢測效率。相比之下,柴油機表面振動信號獲取較為方便快捷,且其中蘊含了豐富的缸內(nèi)燃燒信息,可以用于間接檢測缸內(nèi)壓力。

目前利用振動信號檢測缸內(nèi)壓力的研究比較多,方法也是多種多樣[1-4],但通常不進行振動信號產(chǎn)生機理的分析,忽略了振動信號與缸內(nèi)壓力的本質(zhì)聯(lián)系。本文中通過瞬態(tài)動力學仿真,深入分析了不同類型振動信號與缸內(nèi)壓力的對應(yīng)關(guān)系,發(fā)現(xiàn)振動位移信號與缸內(nèi)壓力的變化趨勢最為接近,適合用于缸內(nèi)壓力的精確識別。在此基礎(chǔ)上,利用希爾伯特包絡(luò)和滑動平均,提取了振動位移信號的變化趨勢,并建立了Adaboost_BP神經(jīng)網(wǎng)絡(luò)模型,通過改進權(quán)重系數(shù)和誤差函數(shù),實現(xiàn)了多工況下缸內(nèi)壓力的識別,并對缸內(nèi)壓力異常現(xiàn)象進行了準確檢測。

1 缸內(nèi)燃燒激勵振動響應(yīng)機理分析

針對某12150型柴油機建立了缸蓋裝配體的有限元模型,如圖1所示。由于主要研究對象為缸蓋,為減小計算量,將機體設(shè)置為剛體。為了便于網(wǎng)格劃分,對缸蓋和螺栓也進行了簡化,忽略了部分倒角、凹槽和小的螺栓孔。經(jīng)過檢驗網(wǎng)格無關(guān)性,最終確定計算中采用四面體網(wǎng)格單元,裝配體單元數(shù)為226 115,節(jié)點數(shù)119 829,具體如圖1所示。缸蓋材料為鋁合金ZL702,缸蓋螺栓的材料為42CrMo。

通過臺架試驗得到右1缸缸內(nèi)壓力曲線,將兩個工作循環(huán)的氣缸壓力施加在缸蓋火力面上,壓力載荷方向垂直于火力面向上。施加的氣缸壓力曲線如圖2所示。

設(shè)置邊界條件和螺栓預(yù)緊力,進行瞬態(tài)動力學計算,得到振動加速度響應(yīng),如圖3所示,與實測振動加速度變化趨勢較為接近,表明建立的仿真模型基本合理,可以用來分析缸內(nèi)燃燒特征參數(shù)及其響應(yīng)信號之間的關(guān)系。

1.1 振動位移信號產(chǎn)生機理分析

實測缸內(nèi)壓力信號及其計算振動位移響應(yīng)如圖4所示。由圖可見:二者在變化趨勢上非常相似,尤其是在最大爆發(fā)壓力之前;最大爆發(fā)壓力之后,位移響應(yīng)的波動開始變的劇烈,但趨勢仍舊比較相似。

對缸內(nèi)壓力信號與振動位移信號的變化關(guān)系做如下解釋:缸蓋通過聯(lián)接螺栓與機體固定,隨著缸內(nèi)壓力的升高,缸蓋以及螺栓在活塞軸線方向產(chǎn)生變形,峰值壓力前,該變形量出現(xiàn)一致性增大,但由于作用力持續(xù)增加,不會產(chǎn)生振動;峰值壓力后,缸內(nèi)壓力逐漸降低,各部件的彈性回復(fù)力開始發(fā)揮作用,導(dǎo)致缸蓋系統(tǒng)開始產(chǎn)生振動。由此缸蓋振動位移響應(yīng)以峰值壓力為界分為兩個階段,峰值壓力之前的位移信號與缸內(nèi)壓力相關(guān)性較強,峰值壓力后的位移信號由于回復(fù)力的作用存在一定的振蕩,但整體趨勢仍與缸內(nèi)壓力有較高的相似性。

1.2 振動速度信號產(chǎn)生機理分析

由以上分析,不難推斷振動速度信號應(yīng)該與缸內(nèi)壓力信號的導(dǎo)數(shù)即壓力升高率存在一定的關(guān)系。

計算實測缸內(nèi)壓力信號的壓力升高率,其計算公式[5]為

(1)

式中:PIRi為第i點壓力升高率;pi為第i點壓力值;Δφ為兩采樣點之間的角度,假定發(fā)動機一個循環(huán)內(nèi)的轉(zhuǎn)速恒定,則Δφ=720/一個循環(huán)的采樣點數(shù)。

圖5為由實測缸內(nèi)壓力計算得到的壓力升高率與仿真振動速度的對比。由圖可見:二者在幅值和相位上有一定偏差,但大體趨勢基本一致,這表明振動速度與缸內(nèi)壓力升高率存在較強的對應(yīng)關(guān)系。

1.3 振動加速度信號產(chǎn)生機理分析

振動加速度是振動速度的導(dǎo)數(shù),由此可推斷振動加速度應(yīng)該與壓力升高率的導(dǎo)數(shù)即壓力升高加速度存在對應(yīng)關(guān)系。圖6為計算振動加速度與缸內(nèi)壓力2階導(dǎo)數(shù)的對比。由圖可見,二者在變化趨勢上較為相似,但由于缸蓋、機體系統(tǒng)的響應(yīng)延遲以及其他激勵的干擾,二者在幅值及相位上存在一定差異。

圖7為實測振動加速度與缸內(nèi)壓力2階導(dǎo)數(shù)的對比。由圖可見,二者變化趨勢也較為相似,但由于針閥落座激勵和燃燒激勵的干擾,信號對最大壓力升高率附近的變化反映不夠準確,另外振動加速度對高頻信號較為敏感,導(dǎo)致信號的局部振蕩較劇烈,不利于識別缸內(nèi)燃燒的特征信息。

2 信號分析及預(yù)處理

通過機理分析,不難看出振動位移信號與缸內(nèi)壓力的變化趨勢最為接近,因此應(yīng)采用振動位移識別缸內(nèi)壓力。由于車輛空間限制,不便于采集柴油機缸蓋的振動位移信號,試驗中通常采集缸蓋振動加速度信號。由3種信號的關(guān)系可知,將缸蓋振動加速度進行二次積分即可得到振動位移信號。數(shù)字積分的常用算法有梯形積分公式和Simpson積分公式等,其中梯形公式計算簡單,但誤差稍大,而Simpson公式計算量稍大但積分精度較高,因此本文中采用Simpson積分法,其表達式[6]為

(2)

式中:v(i)為積分后的速度信號;a(i)為采集到的加速度信號;Δt為采樣間隔時間。

振動信號積分時,由于積分初值無法確定,導(dǎo)致積分后的信號包含直流分量,同時由于缸蓋表面溫度多變,應(yīng)力環(huán)境復(fù)雜,會在信號中存在不規(guī)則的趨勢項,趨勢項對變換結(jié)果影響比較突出,可能導(dǎo)致積分結(jié)果完全失真。針對此問題,本文中采用滑動平均法消除趨勢項,該方法具有較高的噪聲減少比,且實現(xiàn)簡單[7]。圖8為振動加速度信號二次積分得到的振動位移信號。由圖可見,受針閥落座和活塞敲擊等激勵的影響,振動位移信號對缸內(nèi)壓力變化的反應(yīng)不太明顯,僅在燃燒始點至最大爆發(fā)壓力段保持了較強的相似度。

為弱化其他激勵引起的局部波動,提取振動位移信號的變化趨勢,首先采用希爾伯特(Hilbert)變換求取振動位移信號的包絡(luò)信號。Hilbert變換是一個定義信號幅值和相位的特殊函數(shù),它能充分體現(xiàn)原始信號的包絡(luò)信息。設(shè)振動信號序列為X(t),則其解析信號表達式Y(jié)(t)為

Y(t)=X(t)+jH(t)=A(t)expjφ(t)

(3)

(4)

式中H(t)為X(t)的希爾波特變換。

從式(3)中可知振動信號的幅值信息全部蘊含于幅值函A(t)中,該函數(shù)直觀反映了激勵(氣缸壓力)的變化信息[8]。圖9為上止點前后振動位移包絡(luò)信號與缸內(nèi)壓力的對比。由圖可見:由于幅值負值的不存在,希爾波特變換增強了信號幅值的變化效應(yīng),較為明顯地反映了氣缸壓力的變化;但由于其他激勵的低頻成分無法徹底清除,包絡(luò)信號局部仍存在較大波動。為進一步突出信號的整體變化趨勢,利用滑動平均求取包絡(luò)信號的趨勢項,具體如圖10所示。由圖可見:趨勢信號的變化更加平滑,與缸內(nèi)壓力的相似度大大增加。

3 缸內(nèi)壓力識別

利用缸蓋振動信號識別氣缸壓力的方法可以歸納為以下3種:傳遞函數(shù)法、倒譜法和神經(jīng)網(wǎng)絡(luò)法。前兩種方法均須假設(shè)研究對象為線性系統(tǒng),且需要獲取系統(tǒng)的傳遞函數(shù)。柴油機缸蓋系統(tǒng)結(jié)構(gòu)復(fù)雜,振動信號在傳遞過程中受到諸多非線性因素的影響,將其假設(shè)為線性系統(tǒng)顯然不合理。BP神經(jīng)網(wǎng)絡(luò)能夠以任意精度逼近任意非線性映射,且結(jié)構(gòu)簡單,在機械故障診斷中得到了廣泛的應(yīng)用。但也存在一些缺陷,如學習收斂速度慢、易陷入局部最小等,同時利用單一神經(jīng)網(wǎng)絡(luò)往往只能構(gòu)建針對某一特定工況的缸內(nèi)壓力識別模型,將其應(yīng)用于多個工況缸內(nèi)壓力的識別,會導(dǎo)致網(wǎng)絡(luò)結(jié)構(gòu)過于復(fù)雜,識別精度低,泛化性差等一系列問題。為此,本文中利用Adaboost算法建立了組合BP神經(jīng)網(wǎng)絡(luò)模型(Adaboost_BP模型),針對多個工況下的缸蓋振動位移信號進行缸內(nèi)壓力識別,取得了較好的效果。

3.1 Adaboost_BP模型的建立方法

Adaboost算法的雛形為Boosting算法,其基本思想就是試圖產(chǎn)生數(shù)個簡單的、精度比隨機猜測略好的弱規(guī)則,再將這些規(guī)則集成構(gòu)造一個高精度的強規(guī)則,這種思想起源于Valiant提出的PAC學習模型[9]。

針對缸內(nèi)壓力識別問題,AdaBoost集成算法將BP神經(jīng)網(wǎng)絡(luò)作為弱預(yù)測器,反復(fù)訓練BP神經(jīng)網(wǎng)絡(luò),進而加權(quán)得到BP神經(jīng)網(wǎng)絡(luò)強預(yù)測器。其算法具體步驟如下。

步驟1:數(shù)據(jù)選擇和網(wǎng)絡(luò)初始化。從樣本空間中隨機抽取m組數(shù)據(jù)作為訓練樣本,并初始化分布權(quán)重Dt(i)=1/m,根據(jù)樣本輸入輸出維數(shù)確定BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu),初始神經(jīng)網(wǎng)絡(luò)的權(quán)值和閾值。

步驟2:弱預(yù)測器預(yù)測。以第t個預(yù)測器為例,用訓練樣本數(shù)組g(t)訓練BP神經(jīng)網(wǎng)絡(luò)并預(yù)測訓練樣本的輸出結(jié)果,得到訓練樣本數(shù)組g(t)的預(yù)測誤差和et,同時計算下一個預(yù)測器中樣本的分布權(quán)重Dt+1(i)。

(5)

(6)

步驟3:計算預(yù)測序列權(quán)重。根據(jù)訓練樣本組g(t)的預(yù)測誤差et計算序列的權(quán)重at:

(7)

(8)

式中:Bt為歸一化因子,目的是在權(quán)重比例不變的情況下視分布權(quán)重和為1。

步驟5:訓練停止條件。當所有訓練樣本的誤差均方值小于指定值時,停止訓練下一個弱預(yù)測器,即

errt(i)

(9)

此時,et=0,at達到最大值0.5,即該預(yù)測器的權(quán)重最大。

步驟6:強預(yù)測函數(shù)。訓練T輪后得到T組弱預(yù)測函數(shù)f(g(t),at),將T組弱預(yù)測器組合得到強預(yù)測函數(shù)h(x)為

(10)

3.2 基于Adaboost_BP模型的缸壓識別

以上止點前20~上止點后30°CA的缸內(nèi)壓力為識別對象,選取振動位移包絡(luò)信號的趨勢項作為神經(jīng)網(wǎng)絡(luò)輸入進行識別。為了統(tǒng)一輸入和輸出向量的維數(shù),以1°CA為單位抽取50個點的趨勢項數(shù)據(jù)及缸壓數(shù)據(jù),并利用式(11)對樣本進行歸一化。

(11)

對于BP網(wǎng)絡(luò),單隱層的網(wǎng)絡(luò)可以完成任意的非線性映射。因此采用3層BP神經(jīng)網(wǎng)絡(luò)進行訓練。隱層神經(jīng)元數(shù)確定是個十分復(fù)雜的問題,通常需要根據(jù)設(shè)計者的經(jīng)驗和多次實驗來確定,不存在理想的解析表達式。隱層神經(jīng)元數(shù)目與問題的要求和輸入/輸出的單元數(shù)均有密切聯(lián)系。隱層神經(jīng)元數(shù)目太多會導(dǎo)致訓練時間過長,容錯性變差。參考式(12)經(jīng)驗公式可知,BP神經(jīng)網(wǎng)絡(luò)輸入層有50個神經(jīng)元,則隱層神經(jīng)元數(shù)目約為6,可初步確定網(wǎng)絡(luò)結(jié)構(gòu)為50-6-50,經(jīng)后期實驗對比最佳網(wǎng)絡(luò)結(jié)構(gòu)為50-8-50。其他神經(jīng)網(wǎng)絡(luò)參數(shù)設(shè)置如表1所示。

n1=log2n2

(12)

式中n1和n2分別為隱層和輸入層神經(jīng)元數(shù)目。

表1 參數(shù)設(shè)置

為了提高網(wǎng)絡(luò)的工況適應(yīng)性,選擇1 200,1 600和2 000r/min空載工況建立神經(jīng)網(wǎng)絡(luò)模型。每個工況選擇10個樣本共30個樣本組成訓練樣本。另外各抽取10個樣本作為測試樣本。首先單獨采用BP神經(jīng)網(wǎng)絡(luò)進行訓練,其誤差曲線如圖11所示,由于輸入輸出數(shù)據(jù)對應(yīng)關(guān)系較好,網(wǎng)絡(luò)經(jīng)過8步訓練即達到預(yù)定誤差并收斂。利用該網(wǎng)絡(luò)對不同工況下未參與訓練的振動位移數(shù)據(jù)進行測試。圖12~圖14為識別缸內(nèi)壓力與實測缸內(nèi)壓力的對比。由圖可見,模型的識別精度較低,泛化性較差,無法滿足多工況缸內(nèi)壓力的檢測要求,分析認為主要是由于不同工況下缸內(nèi)壓力與振動位移的非線性關(guān)系有所差別,而單一神經(jīng)網(wǎng)絡(luò)的泛化能力較弱,無法同時建立多個不同的映射關(guān)系且保持較高的識別精度。

采用Adaboost集成BP神經(jīng)網(wǎng)絡(luò)以提高模型的泛化性和準確度。將BP神經(jīng)網(wǎng)絡(luò)作為弱預(yù)測器,設(shè)定循環(huán)訓練次數(shù)T=10,誤差閾值為c=0.001。

相比單一神經(jīng)網(wǎng)絡(luò)模型(弱預(yù)測器),集成BP神經(jīng)網(wǎng)絡(luò)的識別誤差明顯下降,如圖15所示。模型對不同工況均保持了較高的識別精度,具體識別情況如圖16~圖18所示。由圖可見識別缸內(nèi)壓力與實測缸內(nèi)壓力幾乎重合,泛化性及精確度均有了大幅度提升。

4 結(jié)論

(1) 通過瞬態(tài)動力學仿真,分析了缸內(nèi)壓力作用下缸蓋振動位移、振動速度和振動加速度的產(chǎn)生機理,發(fā)現(xiàn)三者分別與缸內(nèi)壓力、缸內(nèi)壓力升高率和缸內(nèi)壓力升高加速度存在對應(yīng)關(guān)系。

(2) 利用數(shù)字積分和平均濾波處理振動加速度

信號,能夠得到可用的振動位移信號;希爾伯特包絡(luò)和滑動平均處理能夠有效削弱振動位移信號中其他激勵引起的局部波動,準確提取振動位移的趨勢信息。

(3) 振動位移趨勢項的簡潔性大幅度降低了神經(jīng)網(wǎng)絡(luò)輸入的復(fù)雜度,利用它與缸內(nèi)壓力良好的對應(yīng)關(guān)系,建立了Adaboost_BP集成的神經(jīng)網(wǎng)絡(luò)模型,經(jīng)驗證,該模型能夠準確識別不同工況下的缸內(nèi)壓力,相比單一神經(jīng)網(wǎng)絡(luò),泛化性和精確度均有大幅度提高。

[1] 喬新勇,劉建敏,劉瑋.基于振動測量的發(fā)動機氣缸壓縮壓力檢測方法[J].內(nèi)燃機工程,2008,29(4):63-67.

[2] 朱繼軍.沖擊響應(yīng)信號反演柴油機氣缸壓力的研究[D].武漢:武漢理工大學,2007.

[3] 紀少波,程勇,唐娟,等.基于缸蓋振動信號時域特征識別氣缸壓力的研究[J].內(nèi)燃機工程,2008,29(2):76-80.

[4] 奚銀華,林瑞霖,劉伯運.基于獨立分量分析與傳遞函數(shù)的氣缸壓力重構(gòu)[J].船海工程,2011,40(5):82-85.

[5] PEI Operation Manual: Engine Cycle Analysis[G]. Version 892.

[6] Ribeiro J G T, Castro J T P, Freire J L F. New Improvements in the Digital Double Integration Filtering Method to Measure Displacements Using[C]. Proceedings of the 19th International Modal Analysis Conference, CR-ROM, Orlando, Florida,2001.

[7] Sophocles J O. Introduction to Signal Processing[M]. Prentice Hall,1996.

[8] 王珍,李吉,丁子佳.基于局域波和Hilbert變換的柴油機氣缸壓力識別方法[J].內(nèi)燃機學報,2005,23(4),380-383.

[9] Freund Y, Schapire R E. A Short Introduction to Boosting[J]. Journal of Japanese Society for Artificial Intelligence,1999,14(5):771-780.

A Study on the Cylinder Pressure Identification Method for Diesel EngineUnder Multiple Working Conditions Based on Cylinder Head Vibration Signals

Liu Jianmin1, Li Xiaolei1,2, Qiao Xinyong1& Zhang Jie3

1.DepartmentofMechanicalEngineering,AcademyofArmoredForcesEngineering,Beijing100072; 2.77136Troops,PLA,Bishan402760; 3.TheMilitaryDelegacyofPLAtoNo.318Factory,Beijing100053

The transient dynamics calculation of in-cylinder combustion excitation in a 12150 diesel engine is carried out and the correlations between the displacement, velocity and acceleration of cylinder head vibration and the characteristic parameters of in-cylinder combustion are analyzed. Then on these bases, with the digital integration and average filtering of the vibration acceleration measured, vibration displacement signals are obtained, and by using Hilbert envelope and moving average method the trend of vibration displacement is extracted, with which as input parameter, Adaboost_BP integrated neural network model is built. Finally based on the model the cylinder pressures in different working conditions are identified. The results show that the good correlation between the trend of vibration displacement and cylinder pressure as well as the brevity of parameters themselves effectively reduce the complexity of neural networks input, and hence improve the training efficiency of neural network, while the integrated neural networks model can accurately identify cylinder pressures under different working conditions, with its generalization and accuracy greatly increased.

diesel engine; vibration; cylinder pressure; Hilbert envelope; integrated neural network

*裝備預(yù)先研究項目(40402020101)資助。

原稿收到日期為2013年12月31日,修改稿收到日期為2014年3月7日。

猜你喜歡
振動信號模型
一半模型
振動的思考
科學大眾(2023年17期)2023-10-26 07:39:14
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
重要模型『一線三等角』
完形填空二則
振動與頻率
天天愛科學(2020年6期)2020-09-10 07:22:44
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
基于FPGA的多功能信號發(fā)生器的設(shè)計
電子制作(2018年11期)2018-08-04 03:25:42
中立型Emden-Fowler微分方程的振動性
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产日本欧美亚洲精品视| 91在线播放国产| 亚洲男人天堂网址| 久久久精品国产SM调教网站| 亚洲欧美日本国产综合在线 | 日韩精品成人网页视频在线| a色毛片免费视频| 福利国产微拍广场一区视频在线| 亚洲三级片在线看| 五月激情综合网| 99视频在线观看免费| 国产福利小视频高清在线观看| 久久伊人久久亚洲综合| 91麻豆国产视频| 9cao视频精品| 亚洲中文字幕国产av| 美女无遮挡拍拍拍免费视频| 免费一级毛片| 欧洲亚洲一区| 久久五月视频| 99热这里只有精品在线观看| 亚洲国产AV无码综合原创| 秘书高跟黑色丝袜国产91在线 | 精品视频一区在线观看| 露脸一二三区国语对白| 欧美午夜在线视频| 99一级毛片| 色窝窝免费一区二区三区| 美女一区二区在线观看| 麻豆AV网站免费进入| 青青草一区| 国产成人欧美| 在线播放91| 69综合网| 麻豆国产原创视频在线播放 | 欧美日韩成人在线观看| 日韩精品无码免费专网站| 爱色欧美亚洲综合图区| 日韩精品免费一线在线观看| 青青青伊人色综合久久| 毛片一级在线| 狠狠色婷婷丁香综合久久韩国| 亚洲国产精品美女| 日本成人一区| 欧美精品v欧洲精品| 亚洲一区波多野结衣二区三区| 波多野结衣国产精品| 99爱视频精品免视看| 亚洲国产黄色| 国内精品视频| 久久香蕉国产线看观| 欧美成人aⅴ| 67194在线午夜亚洲| 欧美日韩资源| 免费毛片全部不收费的| 国产嫩草在线观看| 色国产视频| 欧美日韩免费观看| 在线不卡免费视频| 亚洲国产成人超福利久久精品| 国产女人在线观看| 激情视频综合网| 日韩123欧美字幕| 日韩精品毛片| AV网站中文| 亚洲国产精品久久久久秋霞影院 | 婷婷色狠狠干| 国产成人精品一区二区三在线观看| 国产视频你懂得| 热九九精品| 国产欧美日韩va另类在线播放| 国产自视频| 国产成年女人特黄特色毛片免| 国产黑丝视频在线观看| 日韩在线永久免费播放| 国产99视频免费精品是看6| 激情爆乳一区二区| 无码国产伊人| 婷婷99视频精品全部在线观看| 性做久久久久久久免费看| 无码国产伊人| 亚洲一级毛片在线观|