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

利用本地頻率線性分量的電力系統(tǒng)功率缺額估計方法

2023-09-21 09:18:42秦穎婕李文博伍雙喜李佳朋譚嫣李宇駿
南方電網(wǎng)技術(shù) 2023年8期
關(guān)鍵詞:發(fā)電機(jī)模型系統(tǒng)

秦穎婕,李文博,伍雙喜,李佳朋,譚嫣,李宇駿

(1. 廣東電網(wǎng)有限責(zé)任公司電力調(diào)度控制中心,廣州 510060;2. 西安交通大學(xué)電氣工程學(xué)院,西安 710049)

0 引言

近年來,新能源發(fā)電以其清潔、可再生等優(yōu)勢得到了大力發(fā)展,隨著電力電子換流器型電源(converter-interfaced generators,CIGs)占比的提升,電網(wǎng)的慣量水平不斷下降,頻率穩(wěn)定問題日漸突出[1-3]。換流器型電源具有控制靈活、響應(yīng)迅速的優(yōu)點,利用這一特點可以快速調(diào)整換流器型電源的出力,減少系統(tǒng)中功率缺額,從而降低系統(tǒng)頻率偏移、提升系統(tǒng)頻率穩(wěn)定性[4-6]。為實現(xiàn)這種控制結(jié)構(gòu),需要開展系統(tǒng)功率缺額的快速估計。系統(tǒng)功率缺額可以通過系統(tǒng)總慣量乘以慣量中心(center of inertia,COI)頻率變化率獲得,在系統(tǒng)總慣量已知的情況下,系統(tǒng)功率缺額估計依賴于COI頻率變化率(rate of change of frequency,ROCOF)的快速計算。

COI 頻率可以通過許多方法獲得,主要分為兩類:單機(jī)等值模型法和多機(jī)模型法。傳統(tǒng)的單機(jī)等值模型將原始多機(jī)系統(tǒng)進(jìn)行單機(jī)等值聚合,僅用平均頻率近似描述系統(tǒng)的動態(tài)過程。文獻(xiàn)[7]對各發(fā)電機(jī)頻率按其慣性常數(shù)做加權(quán)平均,提出了系統(tǒng)平均頻率模型(average system frequency,ASF)。文獻(xiàn)[8]將ASF 模型擴(kuò)展至含頻率控制的風(fēng)力機(jī)組,并計及了負(fù)荷的頻率與電壓特性。為簡化分析,文獻(xiàn)[9]將系統(tǒng)內(nèi)所有發(fā)電機(jī)調(diào)速器等值為單臺調(diào)速器特性,提出了系統(tǒng)頻率響應(yīng)模型(system frequency response,SFR),并導(dǎo)出了系統(tǒng)頻率的解析表達(dá)式?;赟FR 模型,文獻(xiàn)[10]提出了自適應(yīng)低頻減載方法,實現(xiàn)了功率缺額的快速估算與功率缺額在所有切負(fù)荷輪次間的合理分配。文獻(xiàn)[11]進(jìn)一步考慮了系統(tǒng)旋轉(zhuǎn)備用容量、負(fù)荷頻率特性等因素的影響,提出了改進(jìn)的SFR 模型,提高了系統(tǒng)功率缺額的計算精度。文獻(xiàn)[12]考慮風(fēng)機(jī)參與調(diào)頻,推導(dǎo)了風(fēng)機(jī)運(yùn)行在減載和最大功率跟蹤方式下系統(tǒng)的單機(jī)等值頻率模型。采用單機(jī)等值模型法估算系統(tǒng)功率缺額時,認(rèn)為任意母線測量頻率可近似替代系統(tǒng)頻率,這一假設(shè)忽略了系統(tǒng)暫態(tài)過程中多機(jī)間轉(zhuǎn)子的搖擺,可能使功率缺額估算不準(zhǔn)確,甚至?xí)?dǎo)致預(yù)測不平衡功率方向的錯誤。

多機(jī)模型法考慮了發(fā)電機(jī)轉(zhuǎn)子間的相對搖擺和頻率在空間中的分布特性,現(xiàn)在已經(jīng)成為研究這一問題的主流方法。文獻(xiàn)[13]采用慣量中心坐標(biāo)建立模型,通過多機(jī)系統(tǒng)的仿真軌線獲得關(guān)鍵系頻率特性參數(shù),彌補(bǔ)了單機(jī)等值模型法的局限性。文獻(xiàn)[14-15]利用廣域測量測得的發(fā)電機(jī)相關(guān)電氣量數(shù)據(jù),使功率缺額的估算不依賴于數(shù)值仿真。文獻(xiàn)[16-17]通過廣域測量系統(tǒng)測得系統(tǒng)中少部分特定母線的頻率,再通過這些信息估計COI頻率。文獻(xiàn)[18]考慮了風(fēng)機(jī)的虛擬慣量對COI頻率的影響,進(jìn)一步提升了功率缺額估計的精度和魯棒性。但是,此類方法不僅需要通信,而且廣域測量數(shù)據(jù)不一定能保證系統(tǒng)完全能觀,會導(dǎo)致估計的功率缺額不準(zhǔn)確。為了克服這些問題,國內(nèi)外學(xué)者開始探索通過本地頻率快速估計慣量中心頻率,文獻(xiàn)[19]指出本地頻率在其拐點處與COI頻率近似相等,然而這一性質(zhì)只在兩機(jī)系統(tǒng)中得到證明。類似地,文獻(xiàn)[20]提出了基于本地頻率偏移面積的功率缺額估算方法。前一種方法需要等待兩個頻率拐點,后一種方法需要較長的時間窗才能滿足精度要求,這兩種方法耗時較長,因此并不適合系統(tǒng)功率缺額的快速估計。

本文聚焦于基于本地測量信息的COI頻率估算方法。首先,證明了在系統(tǒng)受擾初期,本地頻率變化量的線性分量與COI 頻率變化量相等這一性質(zhì)。然后,基于最小二乘擬合提取本地頻率變化量線性分量實現(xiàn)了快速的功率缺額估計。最后,在新英格蘭39節(jié)點測試系統(tǒng)上驗證了所提方法的有效性。

1 多機(jī)系統(tǒng)線性化模型

1.1 COI頻率與系統(tǒng)不平衡功率

在一個有n臺發(fā)電機(jī)的電力系統(tǒng)中,描述每臺發(fā)電機(jī)動態(tài)行為的轉(zhuǎn)子運(yùn)動方程可以表示為:

式中:Hi為第i臺發(fā)電機(jī)的慣量時間常數(shù);ωi為第i臺發(fā)電機(jī)的角速度;ωB為系統(tǒng)的基準(zhǔn)頻率;δi為第i臺發(fā)電機(jī)的功角;Pmi和Pei分別為第i臺發(fā)電機(jī)的機(jī)械功率和電磁功率標(biāo)幺值;n為發(fā)電機(jī)數(shù)量。在標(biāo)幺制下,角速度與頻率的大小相等,下文中將用發(fā)電機(jī)頻率代指角速度。為了便于描述多機(jī)系統(tǒng)在受擾后的動態(tài)行為,常在慣量中心坐標(biāo)下進(jìn)行建模分析。將系統(tǒng)內(nèi)所有發(fā)電機(jī)轉(zhuǎn)子運(yùn)動方程中的第一式相加,可得COI搖擺方程為:

式中:HCOI為系統(tǒng)總慣量;ωCOI為慣量中心頻率,其定義為:

在系統(tǒng)受擾初期,可以認(rèn)為調(diào)速器還未啟動,發(fā)電機(jī)的機(jī)械功率不發(fā)生變化。因此,將式(2)改寫為如下形式。

式中:ΔωCOI為COI 頻率變化量;Pd為系統(tǒng)功率缺額。將式(3)左右兩邊同時從0~t時刻進(jìn)行積分,可得:

式中:Δωi為第i臺發(fā)電機(jī)的頻率偏差。式(4)表明,COI 頻率變化量與系統(tǒng)不平衡功率及系統(tǒng)總慣量有關(guān),且在系統(tǒng)受擾初期呈線性變化。

1.2 多機(jī)系統(tǒng)線性化模型

本文主要關(guān)注系統(tǒng)受擾初期的動態(tài)行為。首先,將所有負(fù)荷等效為阻抗,并入節(jié)點導(dǎo)納矩陣中,使用Kron 節(jié)點消除技術(shù)消除網(wǎng)絡(luò)中除去發(fā)電機(jī)內(nèi)電勢節(jié)點外的所有節(jié)點[21-22],消除后的節(jié)點導(dǎo)納矩陣記為Y。對于發(fā)電機(jī)內(nèi)電勢節(jié)點,其節(jié)點注入功率與節(jié)點注入電流及節(jié)點電壓之間的關(guān)系為:

各節(jié)點注入電流可由節(jié)點電壓方程計算。

式中Y=。

將式(6)取共軛后代入式(5)可得:

式中:下標(biāo)i、j為發(fā)電機(jī)的序號;Yij為導(dǎo)納矩陣Y第i行j列的元素;Qei為發(fā)電機(jī)i注入電網(wǎng)的無功功率。電壓向量用極坐標(biāo)形式表示為:

式中:Ei和δi分別為發(fā)電機(jī)內(nèi)電勢電壓的有效值與相角,并且該電壓相角等于發(fā)電機(jī)功角。將導(dǎo)納矩陣中的元素用直角坐標(biāo)的形式表示。

則式(5)可以寫為:

對式(10)使用歐拉公式,并進(jìn)行虛實部分離,則在擾動發(fā)生前,系統(tǒng)中各臺發(fā)電機(jī)的電磁功率可以由式(11)給出。

式中:Ei、Ej分別為第i臺、第j臺發(fā)電機(jī)的內(nèi)電勢;Gij和Bij分別為發(fā)電機(jī)i、j之間的轉(zhuǎn)移電導(dǎo)和轉(zhuǎn)移電納;δij為發(fā)電機(jī)i、j的功角差。通過泰勒展開將上式線性化,可得:

式中:δij(0)為發(fā)電機(jī)i、j的功角差初始值;Δδij為發(fā)電機(jī)i、j功角差的變化量。根據(jù)式 (12),可得每臺發(fā)電機(jī)的電磁功率變化量為:

式中:

式中:Pei(0)為電磁功率的初值,在穩(wěn)態(tài)時該值等于發(fā)電機(jī)的機(jī)械功率。

聯(lián)立式(1)和式(13)可得:

記為:

式(15)中矩陣A(0)含有2n個特征根。一般來說,如果忽略發(fā)電機(jī)阻尼,系統(tǒng)正常運(yùn)行時矩陣A(0)有兩個零根、n-1 對共軛復(fù)根[23]。其中一個零根是由式(16)引起的。

另一個零根則是由忽略系統(tǒng)阻尼引起的。通過選取一臺機(jī)作為參考,可以消去這兩個零根,不妨假設(shè)第n臺機(jī)為參考機(jī)。將式(15)的前n-1行與第n行相減,構(gòu)造新的狀態(tài)向量,式(15)可化為2n-1階微分方程,并消去其中一個零根,化簡后的微分方程如式(17)所示。

式中:

將式(17)第n行到第2n-2 行與第2n-1 行相減,構(gòu)造新的狀態(tài)向量,式(17)可化為2n-2 階微分方程,消去另外一個零根,化簡后的微分方程如下。

式中 :A的表達(dá)式如式(18B)所示;x=[Δδ1n…。

2 基于本地頻率的功率缺額快速估計方法

系統(tǒng)發(fā)生線路故障、發(fā)電機(jī)脫網(wǎng)等事故時,會產(chǎn)生功率缺額??焖俟浪愎β嗜鳖~可以評估事故規(guī)模,并指導(dǎo)頻率穩(wěn)定控制。現(xiàn)有功率缺額估算方法大都依賴通信以獲取系統(tǒng)慣量中心頻率,存在通信延時且可靠性較低。針對這一問題,本節(jié)提出了基于本地頻率信息的系統(tǒng)功率缺額估算方法。

2.1 本地頻率線性分量性質(zhì)

在系統(tǒng)受擾后,網(wǎng)絡(luò)結(jié)構(gòu)的改變使得發(fā)電機(jī)的電磁功率產(chǎn)生突變并與機(jī)械功率產(chǎn)生偏差,即Pe(i0)≠Pm(i0)。因此,式(18)可以重寫為:

對式(19)求一階導(dǎo)使得原微分方程齊次化。

由于發(fā)電機(jī)的功角和轉(zhuǎn)速均不能突變,因此x的初值為零?;谑剑?9),x?的初值為:

為了求解微分方程初值問題式(20)、(21),構(gòu)造如下的線性變換。

式中:XR∈C2(n-1)×2(n-1)為由矩陣A的所有右特征向量構(gòu)成的矩陣。基于式(22)和(20),可得:

式中:λ1,λ2, ???,λ2n-2為矩陣A的特征值。由于狀態(tài)變量z中各分量解耦,式(23)中含有2n-2 個獨立的一階線性微分方程,易得式(23)的解為:

基于式(22)和(23),原方程的解為:

矩陣A有n-1 對共軛復(fù)根,對應(yīng)著系統(tǒng)的n-1個機(jī)電模式。根據(jù)式(25),各發(fā)電機(jī)頻率變化量與參考機(jī)頻率變化量之差均可表示為n-1 個機(jī)電模式的線性組合。

式中:c2,k為第k個振蕩模式的振蕩頻率;c1,ik和c3,ik分別為第k個振蕩模式在Δωin中正、余弦分量的幅值。

將式(26)中第i個方程乘以Hi,然后再將這n-1個方程相加可得:

式(27)左右兩端同時除以HCOI,并將式(4)代入,則第n臺機(jī)的頻率變化量為:

類似地,選擇不同的發(fā)電機(jī)作為參考機(jī),每臺發(fā)電機(jī)的頻率變化量都可以被表示為:

式中:c4,ik和c5,ik分別為第k個振蕩模式在Δωi中的幅值和相位。由式(29)和式(4)可知,本地頻率變化量的線性分量與COI頻率變化量相等。

根據(jù)前文所述,在電力系統(tǒng)中,本地頻率變化量的線性分量與擾動后一段時間內(nèi)的COI頻率變化量相等。基于這一性質(zhì),通過提取任意本地頻率的線性分量,可以快速獲取COI頻率變化率并實現(xiàn)系統(tǒng)功率缺額估計。

2.2 基于最小二乘擬合的系統(tǒng)功率缺額估計方法

通過最小二乘擬合,可以提取任意本地測量頻率的線性分量。具體來講,需要構(gòu)造一個擬合函數(shù)p(x),作為定義在點集X={x1,x2, ???,xm}上列表函數(shù)f(x)的近似表達(dá)式,并使得擬合函數(shù)與原列表函數(shù)的誤差在某種意義上達(dá)到最小。最小二乘擬合技術(shù)以誤差二范數(shù)平方為衡量指標(biāo),用這種方式求得的擬合函數(shù)在每個點上與原列表函數(shù)的距離di=p(xi)-f(xi)的平方和最小。

令{(t1,y1),(t2,y2),···,(tm,ym)}為采樣得到的本地頻率變化量的時間序列。由式(29)可知,本地頻率變化量的線性分量可通過求解以下優(yōu)化問題獲得。

式 中 :t=[t1,t2, ???,tm]T,y=[y1,y2, ???,ym]T為采樣序列,k,ai,bi,ci,i∈{1,2...n-1}為決策變量。式(30)是一個非線性優(yōu)化問題,可以使用梯度下降算法進(jìn)行求解。梯度下降是迭代法的一種,可以用于求解非線性最小二乘擬合問題,其主要思想是沿著函數(shù)梯度下降的方向進(jìn)行搜索,經(jīng)過多步迭代后收斂到函數(shù)最小值附近,因其簡單、求解效率高等優(yōu)點,已成為求解無約束優(yōu)化問題時最常采用的方法之一。

根據(jù)式(4),在系統(tǒng)慣量已知的情況下,系統(tǒng)的功率缺額估計由下式給出。

式中?為由式(30)得到的COI頻率變化率。

3 算例分析

為了驗證所發(fā)現(xiàn)特性的正確性和所提系統(tǒng)功率缺額估計方法的有效性,在新英格蘭39 節(jié)點系統(tǒng)上進(jìn)行了測試。該系統(tǒng)額定頻率為60 Hz,容量基準(zhǔn)值為100 MVA,總裝機(jī)容量為7 367 MW,總有功負(fù)荷為6 254.2 MW,負(fù)荷模型分別采用恒阻抗與恒功率模型。系統(tǒng)拓?fù)淙鐖D1 所示。測試系統(tǒng)的詳細(xì)參數(shù)可參見文獻(xiàn)[24-25]。

圖1 新英格蘭39節(jié)點系統(tǒng)拓?fù)涫疽鈭DFig.1 Sketch diagram of New England 39-bus system

首先,負(fù)荷模型采用恒阻抗模型。圖2 為0 時刻時節(jié)點21 處發(fā)生負(fù)荷突增,系統(tǒng)產(chǎn)生3.6 p. u.有功功率缺額時各發(fā)電機(jī)的頻率響應(yīng)。從圖1 可以看出,在擾動發(fā)生后1 s 內(nèi),COI 頻率呈線性下降,每臺發(fā)電機(jī)的頻率均在COI頻率附近振蕩。

圖2 負(fù)荷突增后測試系統(tǒng)頻率響應(yīng)Fig.2 Frequency response of the test system following a sudden load increase

本文所提的快速功率缺額估計中,時間序列的采樣頻率為1 kHz,時間窗長為0.8 s。7 號發(fā)電機(jī)的動態(tài)響應(yīng)如圖3 所示。從圖中可以看出通過擬合得到的本地頻率(黑色虛線)非常接近實際本地頻率(綠色實線),通過求解式(30)提取的本地頻率的線性分量(橙色虛線)幾乎與COI 頻率(紅色實線)相等。

圖3 負(fù)荷突增時7號發(fā)電機(jī)頻率響應(yīng)Fig.3 Frequency response of SG7 following a sudden load increase

除此之外,文獻(xiàn)[19]通過本地頻率每兩個相鄰拐點的連線來近似COI 頻率曲線,其估計結(jié)果(紫色實線)如圖3 所示。采用該方法估算的COI 頻率與實際COI頻率的偏差較大。如果連接之后的相鄰拐點,甚至?xí)a(chǎn)生錯誤的結(jié)果。

表1展示了節(jié)點24處負(fù)荷有功功率突增時,系統(tǒng)功率缺額估計結(jié)果。

表1 節(jié)點24負(fù)荷突增時的系統(tǒng)功率缺額估計結(jié)果Tab.1 System power deficit estimation result following a sudden load increase at bus No.24

如表1 所示,不同發(fā)電機(jī)節(jié)點估算的功率缺額與真實值的誤差基本在10%以內(nèi)。估計誤差的產(chǎn)生原因主要包括兩個方面:1)分析模型中線性化引入的誤差;2)忽略受擾初期機(jī)械功率變化引入的誤差。值得注意的是,進(jìn)行功率缺額估計的主要目的是為快速調(diào)節(jié)新能源出力提供依據(jù),減小功率缺額,從而降低頻率偏移,現(xiàn)有方法均不能精準(zhǔn)補(bǔ)償功率缺額,但可以顯著改善頻率暫態(tài)特性。應(yīng)用于旨在為新能源調(diào)頻提供參考的快速不平衡功率估計時,估計誤差在可接受范圍內(nèi)。

負(fù)荷采用恒功率模型時,所提算法仍然能較為精確地估計有功功率缺額。圖4、圖5 分別為負(fù)荷采用恒功率及恒阻抗建模、節(jié)點15 處負(fù)荷有功突增導(dǎo)致系統(tǒng)中產(chǎn)生 3.0 p.u. 有功缺額時,9 號發(fā)電機(jī)的頻率響應(yīng)和擬合結(jié)果,該場景下所提算法誤差僅為約+6.24%和+5.51%。從圖4 和圖5 的對比中可以看出,所提算法對不同的負(fù)荷建模有較好的適應(yīng)性。

圖4 負(fù)荷突增時9號發(fā)電機(jī)頻率響應(yīng)(恒功率負(fù)荷)Fig.4 Frequency dynamic of SG9 following a sudden load increase (constant power load)

圖5 負(fù)荷突增時9號發(fā)電機(jī)頻率響應(yīng)(恒阻抗負(fù)荷)Fig.5 Frequency response of SG9 following a sudden load increase (constant impedance load)

圖6、圖7 分別為負(fù)荷采用不同模型,節(jié)點15處負(fù)荷突增導(dǎo)致系統(tǒng)中產(chǎn)生 5.0 p.u. 有功功率缺額時9 號發(fā)電機(jī)的頻率響應(yīng)和擬合結(jié)果。兩種情況下,所提方法的不平衡功率估計誤差分別為+1.38%和+4.45%,從擬合結(jié)果可以看出,該方法對不平衡功率大小有較好的適應(yīng)性。

圖6 負(fù)荷突增時9號發(fā)電機(jī)頻率響應(yīng)(恒功率負(fù)荷)Fig.6 Frequency response of SG9 following a sudden load increase (constant power load)

圖7 負(fù)荷突增時9號發(fā)電機(jī)頻率響應(yīng)(恒阻抗負(fù)荷)Fig.7 Frequency dynamic of SG9 following a sudden load increase (constant impedance load)

表2為負(fù)荷采用恒功率模型、節(jié)點15處負(fù)荷突增并產(chǎn)生-3.887 p. u. 實際有功缺額時由各發(fā)電機(jī)頻率估計所得的有功功率缺額估計結(jié)果,此場景下不同發(fā)電機(jī)處估計的功率缺額與真實值的誤差仍然在10%以內(nèi)。如表1—2 所示,由各發(fā)電機(jī)頻率估計所得的功率缺額與實際值的平均誤差為-4.51%,該方法具有很高的應(yīng)用價值。

表2 節(jié)點15負(fù)荷突增時的系統(tǒng)功率缺額估計結(jié)果Tab.2 Estimation result of system power deficit following a sudden load increase at bus No.15

此外,當(dāng)系統(tǒng)功率過剩時,所提方法仍能較好的估計出過剩功率的大小。圖8所示為0時刻節(jié)點15處負(fù)荷突降,系統(tǒng)中產(chǎn)生2.9 p.u. 有功過剩時所有發(fā)電機(jī)的頻率響應(yīng)。從圖8 可以看出,當(dāng)系統(tǒng)中存在功率過剩時,在擾動后1 s內(nèi)COI頻率基本呈線性上升,每臺發(fā)電機(jī)的頻率均在COI頻率附近振蕩。

圖8 負(fù)荷突降后測試系統(tǒng)頻率響應(yīng)Fig.8 Frequency response of the test system following a sudden load decrease

該場景下6號發(fā)電機(jī)的頻率響應(yīng)如圖9所示。

圖9 負(fù)荷突降時6號發(fā)電機(jī)頻率響應(yīng)Fig.9 Frequency response of SG6 following a sudden load decrease

從圖9 可以看出,擬合曲線與本地頻率曲線基本重合;6 號發(fā)電機(jī)頻率曲線中提取出的線性分量斜率為1.782×10-3p.u./s,與慣量中心頻率變化率十分接近。所提方法也適用于系統(tǒng)中存在功率過剩時的不平衡功率估計。

4 結(jié)語

本文首先證明了在系統(tǒng)受擾初期所有發(fā)電機(jī)頻率變化量的線性分量與COI頻率變化量相等。利用這一性質(zhì),提出了一種基于本地頻率的系統(tǒng)功率缺額估計方法。

該算法利用最小二乘擬合提取本地頻率的線性分量,實現(xiàn)了無需實時通信的快速系統(tǒng)功率缺額估計,有利于充分利用CIGs 的快速功率調(diào)節(jié)能力,提高系統(tǒng)頻率穩(wěn)定性。最后在新英格蘭39 節(jié)點測試系統(tǒng)上證明了該方法的有效性。

猜你喜歡
發(fā)電機(jī)模型系統(tǒng)
一半模型
Smartflower POP 一體式光伏系統(tǒng)
WJ-700無人機(jī)系統(tǒng)
ZC系列無人機(jī)遙感系統(tǒng)
北京測繪(2020年12期)2020-12-29 01:33:58
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
連通與提升系統(tǒng)的最后一塊拼圖 Audiolab 傲立 M-DAC mini
大型發(fā)電機(jī)勵磁用旋轉(zhuǎn)變換器的開發(fā)和應(yīng)用
3D打印中的模型分割與打包
隨身攜帶的小發(fā)電機(jī)
軍事文摘(2016年16期)2016-09-13 06:15:49
主站蜘蛛池模板: 茄子视频毛片免费观看| 最新精品久久精品| 99草精品视频| 精品久久久久久成人AV| 欧美精品在线视频观看| 国产精品亚洲五月天高清| 国产在线视频福利资源站| 精品1区2区3区| 色综合天天综合中文网| 激情乱人伦| 999国内精品视频免费| 狠狠躁天天躁夜夜躁婷婷| 欧美成人午夜视频| 免费在线看黄网址| 亚洲黄色视频在线观看一区| 4虎影视国产在线观看精品| 日韩高清欧美| 久久综合婷婷| 精品丝袜美腿国产一区| 亚洲精选高清无码| 99热这里只有精品2| 97精品伊人久久大香线蕉| 亚洲日韩高清在线亚洲专区| AV无码一区二区三区四区| 91精品国产91久无码网站| 亚洲日本一本dvd高清| 国产99久久亚洲综合精品西瓜tv| 国产人人乐人人爱| 国产成人禁片在线观看| 色九九视频| 久久久久久午夜精品| 精品国产成人av免费| 3D动漫精品啪啪一区二区下载| 成人综合在线观看| 成人毛片免费在线观看| 19国产精品麻豆免费观看| 一级黄色片网| 亚洲女同欧美在线| 久久香蕉国产线看观看亚洲片| 久久99蜜桃精品久久久久小说| 亚洲色图另类| 日韩人妻精品一区| 亚洲bt欧美bt精品| 欧美一区二区自偷自拍视频| 999在线免费视频| 区国产精品搜索视频| 国产区免费精品视频| 色亚洲成人| 欧美精品亚洲二区| 国产精品第5页| 欧美成人免费| 亚洲中文精品人人永久免费| 久久精品中文字幕免费| 亚洲男女天堂| 五月婷婷丁香综合| 国产一级妓女av网站| 欧美综合一区二区三区| 97青草最新免费精品视频| 美女高潮全身流白浆福利区| 亚洲香蕉在线| 蜜芽一区二区国产精品| 婷婷色一二三区波多野衣 | 91精品国产情侣高潮露脸| 噜噜噜久久| 少妇极品熟妇人妻专区视频| 国产精品不卡永久免费| 色婷婷电影网| 98超碰在线观看| 国产97色在线| 91精品小视频| 全部无卡免费的毛片在线看| 国产福利一区在线| 全部无卡免费的毛片在线看| www.91中文字幕| 全部无卡免费的毛片在线看| 日本一区二区三区精品AⅤ| 2020久久国产综合精品swag| 国产激情无码一区二区APP | 99精品伊人久久久大香线蕉| 国产另类乱子伦精品免费女| 精品国产一区二区三区在线观看| 99福利视频导航|