亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        一種帶有反饋濾波器的自整定Smith預(yù)估器設(shè)計及其在紙張橫幅定量控制中的應(yīng)用

        2023-03-31 04:30:06劉文波徐大海
        中國造紙 2023年2期
        關(guān)鍵詞:模型系統(tǒng)

        湯 偉 張 旭 劉文波 徐大海 王 杰

        (陜西科技大學(xué)電氣與控制工程學(xué)院,陜西西安,710021)

        紙張定量控制技術(shù)包含縱向定量控制技術(shù)和橫幅定量控制技術(shù),其中縱向是指紙機運行方向,橫幅即為垂直紙機運行方向[1-2]??v向定量控制技術(shù)已逐漸成熟,但是不能保證紙張定量的均勻分布,而實施橫幅定量控制技術(shù)可以使紙張的勻度更高。但由于稀釋水閥執(zhí)行器之間存在強耦合、系統(tǒng)存在大時滯等特性使得橫幅定量控制變的更加困難。工程實踐中傳統(tǒng)的比例積分微分(PID)控制器參數(shù)一般是人為設(shè)定,但在面對這種具有不確定性的復(fù)雜對象時,難以獲得理想的控制效果。有學(xué)者采用模糊控制器整定PID參數(shù),雖然提高了系統(tǒng)的抗擾能力,但在面對這種強耦合、大時滯控制對象時,系統(tǒng)動態(tài)性能依然不佳[3]。鑒于Smith預(yù)估器具有較好的控制效果,同時受模糊控制不需要對控制對象進(jìn)行精確建模的啟發(fā),本課題提出將Smith預(yù)估器與模糊自整定PID相結(jié)合的控制方案來解決紙張橫幅定量控制所存在的控制難題。

        本課題首先對帶有稀釋水水力式流漿箱的長網(wǎng)造紙工藝流程進(jìn)行闡述,分析紙張橫幅定量控制所存在的控制難點。其次建立橫幅定量控制過程的數(shù)學(xué)模型,針對其關(guān)聯(lián)矩陣為非方高維矩陣的特點,使用“區(qū)域化點”的方法對系統(tǒng)進(jìn)行降維,并設(shè)計解耦逆矩陣以將系統(tǒng)轉(zhuǎn)換為與稀釋水閥數(shù)量一致的n個單回路。最后針對解耦后的多個單回路存在大時滯的特性,設(shè)計帶有反饋濾波器的Smith預(yù)估器加模糊自整定PID復(fù)合控制方案,并取得良好的控制效果。

        1 帶有稀釋水水力式流漿箱的長網(wǎng)造紙工藝流程

        圖1為稀釋水流漿箱定量控制工藝圖。如圖1所示,在一般的橫幅定量控制工藝中,稀釋水閥位于紙機頭部,沿流漿箱橫幅等距安裝,定量檢測裝置(掃描架)則安裝在紙機尾部,掃描架上的掃描探頭沿紙機橫幅往復(fù)運動,所采集到的瞬時采樣值便是橫幅定量數(shù)據(jù)。稀釋水流漿箱的定量調(diào)節(jié)過程為:當(dāng)尾部的掃描架檢測到某處定量與設(shè)定的定量值出現(xiàn)差距時,改變紙機頭部的稀釋水閥開度,從而相對應(yīng)的階梯擴散管的稀釋水注入量得到改變,最終產(chǎn)生的效果即為該處的漿料流量與稀釋水流量的比率產(chǎn)生變化,即通過調(diào)節(jié)漿料濃度來達(dá)到使該處定量達(dá)到設(shè)定值的目的。因此采用稀釋水流漿箱的橫幅定量控制系統(tǒng)存在以下控制難點。

        圖1 稀釋水流漿箱定量控制工藝圖[4]Fig.1 Basis weight control process diagram of headbox with dilution water[4]

        1.1 非方高維耦合特性

        在橫幅定量控制系統(tǒng)中,一般沿流漿箱橫幅每60 mm安裝一支稀釋水閥,水閥數(shù)量隨著紙幅寬度增大而增多,一般在幾十到幾百只,而紙機尾部的橫幅定量數(shù)據(jù)一般在幾百個,甚至數(shù)千個,這樣就構(gòu)成了一個非方的高維系統(tǒng),對非方高維系統(tǒng)進(jìn)行求解計算量巨大。當(dāng)改變稀釋水閥開度使稀釋水閥動作時,由于流體在流動時具有沖擊特性,會對該稀釋水閥兩側(cè)的鄰近區(qū)域產(chǎn)生影響,即影響鄰近區(qū)域的定量值,在控制上表現(xiàn)為耦合特性[4]。

        1.2 時滯特性

        定量檢測點位于紙機尾部,稀釋水閥位于紙機頭部的流漿箱上,數(shù)據(jù)傳輸與稀釋水閥動作輸出之間橫跨整個紙機,系統(tǒng)存在較大滯后時間,而滯后時間與紙機長度和定量數(shù)值成正比,與紙機車速成反比,在控制上表現(xiàn)出大時滯特性。

        傳統(tǒng)的控制方式難以解決上述提及的控制難點,因此本課題提出設(shè)計一種帶有反饋濾波器的Smith預(yù)估器的模糊自整定PID控制方案來解決上述存在的控制難題。

        2 橫幅定量控制過程數(shù)學(xué)模型

        對于一套含有m個測量點和n個稀釋水閥的橫幅定量控制系統(tǒng),其數(shù)學(xué)模型如式(1)所示[5]。

        式中,Ycd(s)為橫幅檢測點的輸出值;yi(s)為沿紙機橫幅第i(i=1,2,…,n)個檢測點的測量定量值;U(s)為稀釋水閥執(zhí)行輸出的控制量;uj(s)為第j(j=1,2,…,m)個稀釋水閥的輸出控制量;Gt(s)為橫幅定量控制系統(tǒng)中的關(guān)聯(lián)耦合矩陣;gij為第j個稀釋水閥動作時對下游第i個測量點的影響系數(shù);g0(s)為稀釋水閥輸出到紙機尾部檢測點之間的傳遞函數(shù);T為系統(tǒng)慣性時間常數(shù);l為系統(tǒng)的滯后時間。

        橫幅定量系統(tǒng)的關(guān)聯(lián)矩陣是一個m×n的非方矩陣,若想對其求解則需要進(jìn)行降維化方,才能降低解耦運算量。采用“區(qū)域化點”的思想來對系統(tǒng)進(jìn)行降維,即1只稀釋水閥控制1個測量點區(qū)域,假設(shè)1只稀釋水閥影響相鄰兩側(cè)各d個測量點,則系統(tǒng)的耦合寬度為p=2d+1。若通過引入1個前置矩陣對系統(tǒng)進(jìn)行降維,降維計算過程如式(2)所示。

        式中,gp+1+n/m,i為第i個稀釋水閥對應(yīng)的正下游測量點的位置。則降維后的關(guān)聯(lián)矩陣形式即為式(3)所示。

        根據(jù)研究發(fā)現(xiàn),某個稀釋水閥執(zhí)行器動作,其對紙機尾部的測量點產(chǎn)生的影響是中心對稱的,則有g(shù)'ij=g'ji,因此關(guān)聯(lián)矩陣可進(jìn)一步簡化為式(4)。

        式中,g'0代表該稀釋水閥動作對正下游測量區(qū)域的影響系數(shù);g'q代表該稀釋水閥動作對鄰側(cè)第q個稀釋水閥對應(yīng)的測量區(qū)域的影響系數(shù)。

        這樣,新的關(guān)聯(lián)矩陣維數(shù)就變成了與稀釋水閥數(shù)量一致的方形矩陣,降維后的系統(tǒng)數(shù)學(xué)模型如式(5)所示。

        而降維后的關(guān)聯(lián)矩陣需要對其設(shè)計解耦補償矩陣,將降維后的系統(tǒng)轉(zhuǎn)變成對角系統(tǒng),這樣便可分解成與稀釋水閥數(shù)量一致的n個單回路,然后對n個單回路分別設(shè)置控制器進(jìn)行控制,使得控制難度大大降低。

        設(shè)解耦器為N(s),則解耦后系統(tǒng)的關(guān)聯(lián)矩陣形式如式(6)所示。

        根據(jù)式(6)計算可得N(s)G'(s)=G0(s)。若G'(s)為非奇異矩陣,則解耦器可設(shè)置為N(s)=G'-1(s),將其設(shè)置在被控對象與控制器之間,通過解耦器作用于被控對象,將其轉(zhuǎn)變?yōu)榫哂衝個單回路的主對角系統(tǒng),然后通過控制器完成對n個單回路的單獨控制。解耦器的解耦結(jié)構(gòu)如圖2所示。

        圖2 橫幅定量系統(tǒng)解耦框圖Fig.2 Block diagram of decoupling structure of the cross direction basis weight system

        3 模糊自整定PID-Smith控制器設(shè)計

        解耦后的系統(tǒng)為與稀釋水閥數(shù)量一致的n個單回路,由于稀釋水閥之間具有相同的物理特性,因此只需完成1個回路的設(shè)計,同理便可完成其余回路的設(shè)計。

        本研究將模糊控制與PID控制相結(jié)合,同時對Smith預(yù)估器進(jìn)行改進(jìn),即采用改進(jìn)后的Smith預(yù)估器與模糊自整定PID相結(jié)合來加強對橫幅定量控制系統(tǒng)的大時滯過程的控制效果。圖3為系統(tǒng)結(jié)構(gòu)框圖。

        圖3 系統(tǒng)控制框圖Fig.3 Block diagram of system control

        3.1 帶有反饋濾波器的Smith預(yù)估器

        常規(guī)Smith預(yù)估器在理想情況下(即Gm(s)e-lms=Gp(s)e-ls)可以完全消除時滯帶來的不利影響,但Smith預(yù)估器存在的缺點是要求精確建立控制對象的預(yù)估數(shù)學(xué)模型,否則控制效果會大打折扣[6-8]。而在實際應(yīng)用時,大多控制對象都難以進(jìn)行精確建模,因此本研究采用帶有反饋濾波器的Smith預(yù)估器(FuzzyPID-FSmith),其結(jié)構(gòu)圖如圖4所示。

        圖4 帶有反饋濾波器的Smith預(yù)估器結(jié)構(gòu)圖Fig.4 Structure of Smith predictor with feedback filter

        由結(jié)構(gòu)框圖可以推導(dǎo)出輸出y(s)到輸入r(s)的系統(tǒng)傳遞函數(shù)如式(7)所示。

        式中,Gc(s)代表控制器;e-ls代表被控對象的時滯環(huán)節(jié);Gm(s)代表對象模型不含時滯環(huán)節(jié)的模型;Gm(s)e-lms表示預(yù)估模型,在理想狀況下,它等同于實際過程Gp(s)e-ls。

        根據(jù)傳遞函數(shù)可以看出,當(dāng)模型精確時,帶有反饋濾波器的系統(tǒng)傳遞函數(shù)相比傳統(tǒng)Smith預(yù)估器并無區(qū)別[9];而當(dāng)模型失配時,相當(dāng)于在主反饋通道引入一階濾波器,使得模型失配的誤差信號pm(s)經(jīng)過一階濾波器的濾波處理后才能反饋到控制器Gc(s),從而削弱了模型失配產(chǎn)生的模型誤差帶來對系統(tǒng)的影響。根據(jù)研究發(fā)現(xiàn),當(dāng)模型失配時,若Tm取值過小,雖然調(diào)節(jié)時間會變快,但同時也產(chǎn)生了一定的超調(diào),系統(tǒng)穩(wěn)定性變差;若Tm取值過大,雖然會減小系統(tǒng)的超調(diào)量,但反應(yīng)速度也會變慢。根據(jù)多次實驗發(fā)現(xiàn),當(dāng)Tm≥lm時,濾波器才會起到較好的調(diào)節(jié)效果。因此為了兼顧快速性和穩(wěn)定性,本研究濾波時間常數(shù)取值與滯后時間常數(shù)相同,即Tm=lm。

        3.2 模糊自整定PID橫幅定量控制器設(shè)計

        本課題所設(shè)計模糊控制器的2個輸入分別為經(jīng)過量化因子量化后的橫幅定量偏差e(t)=r(t)-y(t)和偏差變化率ec=de(t)dt,輸出分別為kP、kI、kD3個參數(shù)的調(diào)整值,記作ΔkP、ΔkI、ΔkD,紙機尾部所檢測到的橫幅定量值經(jīng)過與設(shè)定的定量值比較,若存在誤差,則計算出該時刻的定量誤差以及其變化率,然后輸送到模糊控制器中并根據(jù)所設(shè)計的模糊規(guī)則來計算出PID參數(shù)的修正量,從而在線調(diào)整PID參數(shù)來滿足系統(tǒng)不同時刻的需求[10-11]。

        3.2.1 模糊集、模糊論域及量化因子選取

        橫幅定量偏差e和ec的基本論域為[-0.5,0.5],模糊論域可選取為[-3,3],通過設(shè)定量化因子可將基本論域映射到模糊論域范圍內(nèi),量化因子為:3/0.5=6。

        圖5為e和ec隸屬度函數(shù)。由圖5可知,e和ec的模糊集分為7個等級,其中NB、NM、NS的含義為負(fù)向大、負(fù)向中和負(fù)向小,分別表示實際定量測量值相比設(shè)定的定量值高得多、較高和高一點,ZE表示實際定量測量值與設(shè)定的定量值相等,PS、PM、PB的含義為正向大、正向中和正向小,分別表示實際定量測量值相比設(shè)定值低一點、較低和低的多。

        圖5 e和ec隸屬度函數(shù)Fig.5 Degree of membership function for e and ec

        圖6為ΔkP、ΔkI、ΔkD的隸屬度函數(shù)。根據(jù)經(jīng)驗輸出變量 ΔkP、ΔkI、ΔkD的基本論域選取為[-0.5,0.5],其模糊論域可取為[-1,1],那么量化因子均為:0.5/1=0.5。輸出變量模糊集同樣分為7個等級,其中NB代表參數(shù)調(diào)整為負(fù)向極大,NM代表參數(shù)調(diào)整為負(fù)向適中,NS代表參數(shù)調(diào)整為負(fù)向極小,ZE則代表參數(shù)調(diào)整為零,PS代表參數(shù)調(diào)整為正向極小,PM代表參數(shù)調(diào)整為正向適中,PB代表參數(shù)調(diào)整為正向極大。

        圖6 ΔkP、ΔkI、ΔkD隸屬度函數(shù)Fig.6 Degree of membership function for ΔkP, ΔkI, ΔkD

        3.2.2 隸屬函數(shù)選取

        考慮到三角隸屬度函數(shù)運算簡單、內(nèi)存空間占用小的特點,以及三角函數(shù)對系統(tǒng)誤差變化反應(yīng)靈敏、分布均勻[12],因此本研究的輸入輸出變量均選取論域均勻劃分的三角函數(shù)作為隸屬函數(shù)。

        3.2.3 紙張橫幅定量控制器模糊規(guī)則選取

        結(jié)合現(xiàn)場參數(shù)整定經(jīng)驗以及前人研究,橫幅定量系統(tǒng)模糊規(guī)則的設(shè)定遵循原則如下。

        當(dāng)定量誤差較大時,此時以盡快消除定量偏差為主要目標(biāo),應(yīng)該選取較大的ΔkP來盡快到達(dá)設(shè)定的定量值,為防止出現(xiàn)超調(diào),ΔkI一般取零即可。

        當(dāng)定量誤差適中時,此時系統(tǒng)到達(dá)響應(yīng)中期,定量偏差變化率ec為中等大小,為了保證系統(tǒng)的響應(yīng)速度同時也要減小超調(diào),應(yīng)該選用較小的ΔkP、ΔkI,而ΔkD適中即可。

        當(dāng)系統(tǒng)到達(dá)響應(yīng)后期時,定量偏差|e|較小,實際定量值已經(jīng)非常接近設(shè)定值,為降低穩(wěn)態(tài)誤差,應(yīng)選取較大的ΔkP以及較小的ΔkI,同時為防止在穩(wěn)態(tài)值附近震蕩,ΔkD適中即可[13-14]。

        綜上所述整定原則,可列模糊規(guī)則表,具體如表1所示。

        表1 ΔkP、ΔkI、ΔkD模糊規(guī)則表Table 1 Table of fuzzy rules for ΔkP, ΔkI, and ΔkD

        根據(jù)所設(shè)計的模糊規(guī)則,模糊控制器進(jìn)行推理運算,實時輸出PID控制器的校正參數(shù),最終得到PID的3個參數(shù)的改變量分別為ΔkP、ΔkI和ΔkD,參數(shù)整定公式如式(8)所示。

        式中,kP0、kI0和kD0為PID參數(shù)的初始值。

        4 仿真研究

        采用Simulink搭建仿真模型,被控對象數(shù)學(xué)模型采用關(guān)聯(lián)矩陣為10×10的托普利茲矩陣形式[15],如式(9)所示。

        對所建立的定量系統(tǒng)模型加單位階躍信號,定量系統(tǒng)的開環(huán)響應(yīng)結(jié)果如圖7所示。由圖7可以看出,橫幅各稀釋水閥之間存在較強耦合關(guān)系,使得橫幅定量嚴(yán)重不均勻,與設(shè)定的定量值之間存在較大誤差。

        圖7 橫幅定量系統(tǒng)開環(huán)響應(yīng)平面圖Fig.7 Open-loop response of cross direction basis weight system

        為驗證本研究所設(shè)計的控制算法的控制效果,分別采用PID控制器、傳統(tǒng)的Smith預(yù)估器(PID-Smith)以及本研究的控制算法(FuzzyPID-FSmith)對定量系統(tǒng)進(jìn)行仿真。用Z-N法得到PID控制器的參數(shù)值kP=2.14、kI=0.06、kD=0.25,采用大林算法得到PIDSmith控制器的參數(shù)值kP=2.14、kI=0.41、kD=0,并將其作為FuzzyPID-FSmith的初始值,即kP0=2.14、kI0=0.41、kD0=0,Smith預(yù)估器的反饋濾波時間常數(shù)一般取值與預(yù)估模型的滯后時間常數(shù)一致tm=8。同時為驗證抗擾動能力分別在t=0 s和t=100 s對系統(tǒng)加入單位階躍信號和幅值為-0.3的階躍負(fù)載擾動。

        4.1 標(biāo)稱情況下系統(tǒng)性能測試

        圖8為系統(tǒng)加入解耦矩陣后,并使用PID控制器的仿真結(jié)果。由圖8可以看出,系統(tǒng)的耦合已經(jīng)被抵消,但由于時滯的存在,導(dǎo)致系統(tǒng)響應(yīng)速度遲鈍,存在著32%的超調(diào)量和長達(dá)50 s的調(diào)節(jié)時間,且在受到擾動后需要較長時間才能恢復(fù)到穩(wěn)態(tài)值,抗擾性能不佳,不適用于橫幅定量的大時滯控制。

        圖8 PID控制器下紙張橫幅平面圖Fig.8 Plan view of cross direction of paper with PID controller

        圖9和圖10分別為在PID-Smith控制器和本研究所設(shè)計的FuzzyPID-FSmith控制器下的仿真結(jié)果。由圖9和圖10可以看出,系統(tǒng)在FuzzyPID-FSmith控制器和傳統(tǒng)的PID-Smith控制器都具有較好的動態(tài)性能,本研究所設(shè)計的控制算法無論是在追蹤設(shè)定定量值方面還是抗擾動方面都略勝一籌,各控制方案的系統(tǒng)動態(tài)性能對比結(jié)果見表2。

        圖9 PID-Smith控制器下紙張橫幅平面圖Fig.9 Plan view of cross direction for paper with PID-Smith controller

        圖10 FuzzyPID-FSmith控制器下紙張橫幅平面圖Fig.10 Plan view of cross direction for paper with FuzzyPIDFSmith controller

        表2 標(biāo)稱情況下性能指標(biāo)Table 2 Performance indicators at nominal conditions

        4.2 模型失配情況下系統(tǒng)性能測試

        為檢驗在模型失配下各控制算法的有效性,將系統(tǒng)模型的增益K、時間常數(shù)L均增大20%,時間常數(shù)T減小20%,然后進(jìn)行仿真。圖11為在傳統(tǒng)的PID控制器下系統(tǒng)的仿真結(jié)果圖。由圖11可知,系統(tǒng)調(diào)節(jié)時間變得更長,抗擾性能更差,系統(tǒng)動態(tài)性能全面下降。圖12為在PID-Smith控制器控制下系統(tǒng)仿真結(jié)果。由于模型失配原因系統(tǒng)已經(jīng)發(fā)散,紙張的橫幅平面起伏嚴(yán)重,控制效果大大降低。在本研究所設(shè)計的FuzzyPID-FSmith控制器下,紙幅平面圖的結(jié)果如圖13所示。由圖13可知,系統(tǒng)雖然產(chǎn)生了超調(diào),但是依然能快速達(dá)到設(shè)定值并進(jìn)入穩(wěn)態(tài),受到干擾后也能在較短時間內(nèi)恢復(fù)到穩(wěn)定狀態(tài),系統(tǒng)的魯棒性得到較大的提升。模型失配下各個控制方案系統(tǒng)的動態(tài)性能指標(biāo)具體如表3所示。

        圖11 PID控制器下紙張橫幅平面圖(模型失配)Fig.11 Plan view of cross direction of paper with PID controller(model mismatch)

        圖12 PID-Smith控制器下紙張橫幅平面圖(模型失配)Fig.12 Plan view of cross direction of paper with PID-Smith controller (model mismatch)

        圖13 FuzzyPID-FSmith控制器下紙張橫幅平面圖(模型失配)Fig.13 Plan view of cross direction of paper with FuzzyPID-Fsmith controller (model mismatch)

        表3 模型失配情況下性能指標(biāo)Table 3 Performance indicators in case of model mismatch

        5 實際應(yīng)用

        以陜西某造紙廠為例,紙機車速600 m/min、幅寬3.5 m、涂布白紙板定量306 g/m2,并由55個稀釋水閥共同控制。圖14為系統(tǒng)運行穩(wěn)定后某個時間段的橫幅定量偏差監(jiān)控曲線。由圖14可知,橫軸坐標(biāo)為55個稀釋水閥的編號,縱軸坐標(biāo)由上到下分別表示橫幅定量差和相對應(yīng)稀釋水閥開度。最大正偏為該處橫幅定量值偏離平均值的正向最大值,最大負(fù)偏為該處橫幅定量值偏離平均值的負(fù)向最大值,2σ的含義為橫幅定量的均方差指標(biāo),其值越小,定量越均勻。

        圖14 現(xiàn)場實際應(yīng)用橫幅定量差曲線Fig.14 Deviation curves of cross direction basis weight for practical application on site

        從圖14可以看出,在應(yīng)用本課題的研究算法后,系統(tǒng)的定量偏差整體波動較小,其橫幅定量差指標(biāo)為6.6/306.1=2%(一般采用2σ/定量平均值的方式來作為橫幅定量差的評價指標(biāo),并要求在5%以內(nèi),同時其指標(biāo)在3%以內(nèi)的為優(yōu)等品),滿足造紙企業(yè)對橫向定量控制的需求,取得了較好的控制效果。

        6 結(jié)論

        本課題通過建立紙張橫幅定量系統(tǒng)數(shù)學(xué)模型,采用“區(qū)域化點”法對紙張橫幅定量系統(tǒng)進(jìn)行降維,并對降維后的系統(tǒng)設(shè)計解耦逆矩陣作為解耦器。設(shè)計了一種帶有反饋濾波器的自整定Smith預(yù)估器(FuzzyPID-FSmith)。仿真結(jié)果表明,本課題所設(shè)計的控制算法調(diào)節(jié)時間快、超調(diào)量小、系統(tǒng)魯棒性強。將本課題所設(shè)計算法投入現(xiàn)場應(yīng)用,將橫幅定量差指標(biāo)控制在2%以內(nèi),取得了較好的控制效果,能夠滿足造紙企業(yè)的需求。

        猜你喜歡
        模型系統(tǒng)
        一半模型
        Smartflower POP 一體式光伏系統(tǒng)
        WJ-700無人機系統(tǒng)
        ZC系列無人機遙感系統(tǒng)
        北京測繪(2020年12期)2020-12-29 01:33:58
        重要模型『一線三等角』
        重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
        基于PowerPC+FPGA顯示系統(tǒng)
        半沸制皂系統(tǒng)(下)
        連通與提升系統(tǒng)的最后一塊拼圖 Audiolab 傲立 M-DAC mini
        3D打印中的模型分割與打包
        中文字幕天天躁日日躁狠狠躁免费| 中文字幕亚洲综合久久久| 国产一区二区三区啊啊| 无码av不卡一区二区三区| av无码免费永久在线观看| 国产日韩精品一区二区在线观看播放 | 无码专区中文字幕DVD| 亚洲国产日韩综一区二区在性色| 国产婷婷成人久久av免费| 品色堂永远免费| 又大又粗弄得我出好多水| 在线你懂| 色和尚色视频在线看网站| 日本一本之道高清不卡免费| 男女超爽视频免费播放| 精品视频在线观看一区二区有| 亚洲精品偷拍自综合网| 啦啦啦www在线观看免费视频| 欧美成人三级一区二区在线观看| 黄色三级视频中文字幕| 日本中文一区二区在线| 国产午夜福利不卡在线观看| 久久国产精品国产精品日韩区| av免费看网站在线观看| 久久久亚洲熟妇熟女av| 国产午夜无码片在线观看影院| 亚洲午夜福利精品久久| 三个黑人插一个女的视频| 欧美激情一区二区三区成人 | 国产特级毛片aaaaaa高潮流水| 国产亚洲精品久久久久婷婷瑜伽| 亚洲中文字幕巨乳人妻| 日本韩国亚洲三级在线| 国产精品嫩草99av在线| 99久久久国产精品免费蜜臀| 国产极品嫩模大尺度在线播放| 国产毛片av最新视频| 久久综合九色综合网站| 99在线国产视频| 国产自拍偷拍视频免费在线观看| 国产精品一区二区久久乐下载 |