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

電力系統電壓穩定鞍結分岔點的快速求取方法

2017-10-09 13:13:21朱永強劉光曄匡一雷康志豪劉媛媛
電力系統及其自動化學報 2017年9期

朱永強,劉光曄,曹 斌,匡一雷,康志豪,劉媛媛

(1.湖南大學電氣與信息工程學院,長沙 410082;2.國網湖南省電力公司岳陽供電分公司,岳陽414000;3.國網湖南省電力公司檢修公司,長沙 410000)

電力系統電壓穩定鞍結分岔點的快速求取方法

朱永強1,2,劉光曄1,曹 斌3,匡一雷1,康志豪1,劉媛媛1

(1.湖南大學電氣與信息工程學院,長沙 410082;2.國網湖南省電力公司岳陽供電分公司,岳陽414000;3.國網湖南省電力公司檢修公司,長沙 410000)

為快速求取電力系統電壓穩定鞍結分岔點,提出一種自動變步長的連續潮流法和局部曲線擬合法相結合的方法。運用廣義戴維南動態等值原理,證明非線性電力網絡傳輸功率達到極值必要條件,定義阻抗模裕度指標,判定系統靜態電壓穩定性最薄弱節點。由當前運行點阻抗模預測下一運行點潮流狀態,得到增長步長,實現自動變步長,快速逼近而不越過鞍結分岔點。當阻抗模裕度值達到滿足局部曲線擬合要求的設定值時,滿足阻抗模裕度判據,對鞍結分岔點精確擬合,得到較準確的鞍結分岔點值。多個IEEE測試系統的仿真表明,該算法簡單、精度高、速度快、計算次數少,且適用于規模較大的電力網絡。

電力系統;電壓穩定;鞍結分岔;廣義戴維南等值;阻抗模裕度判據;局部曲線擬合

Abstract:In order to fast calculate the saddle node bifurcation point of voltage stability of power system,a combination method with automatic variable step is proposed based on continuation power flow method and local curve fitting meth?od.By using generalized Thevenin dynamic equivalence theory,the essential conditions for reaching the extreme values of transmission power in a nonlinear power network are proved.Moreover,impedance modulus margin index is defined to determine the weakest node of static voltage stability of the system.An increasing step length is obtained according to the power flow state at the next operation point predicted by the impedance modulus at the current operation point,thus an automatic variable step is realized,and the saddle node bifurcation point is approached fast but not crossed over.When the impedance modulus margin satisfies the value required by local curve fitting,i.e.,impedance modulus mar?gin criterion,the saddle node bifurcation point will be accurately fitted and a more accurate value will be obtained.Simulations of several IEEE test systems show that the proposed method is simple,accurate,fast and with less computa?tional load,thus it can be applied to large-scale power networks.

Key words:power system;voltage stability;saddle node bifurcation;generalized Thevenin equivalence;impedance modulus margin criterion;local curve fitting

隨著我國電力網絡的日益擴大,電壓穩定問題成為電網安全運行的主要研究方向[1]。電壓穩定研究的重點是反映電壓穩定性的安全指標的確定。裕度指標包含信息量大,且有較好的直觀性和線性,受到電力研究者的高度重視[2]。負荷裕度是指在某種負荷功率增長方式下,基態運行點到電壓極限點的距離,重點是求取鞍結分岔點SNBP(saddle node bifurcation point)[3]。

求取SNBP的方法主要有直接法[4]、非線性規劃法[5]以及連續潮流法 CPF(continuation power flow)[6]。直接法計算量小,且可解決雅可比矩陣奇異問題,但提供有用信息極少,難以滿足運行人員要求。非線性規劃法可避免潮流不收斂,但計算量過大,難以滿足大系統要求。CPF法通過增加參數方程,避免臨界崩潰點潮流發散,但是通過固定步長逐步計算,獲取鼻形曲線,將鼻尖作為SNBP,若步長過小,計算速度太慢,在線應用困難;若步長過大,極易越過極限點,造成迭代不收斂。

為了快速精確地求取SNBP,本文提出一種自動變步長的CPF法,快速逼近卻不越過SNBP,并采用局部曲線擬合法對SNBP進行精確擬合。對多個IEEE測試系統的仿真計算表明,所提方法不僅有效,而且實用。

1 基于阻抗模裕度的系統薄弱節點分析

1.1 廣義戴維南動態等值方法

文獻[7]將一個復雜的線性系統通過戴維南等值原理等值為圖1所示的兩節點簡單系統,并證明電力網絡輸送有功功率達到極值必要條件是負荷靜態阻抗的模值與負荷節點戴維南等值阻抗的模值相等。本文提出了一種命名廣義戴維南的動態等值方法,證明該結論在非線性電力網絡中且負荷的功率因素變化的情況下也正確,且只需一個潮流點就可求解出各種等值參數。

圖1 戴維南等值兩節點系統Fig.1 Two-node system with Thevenin equivalence

式中:?LD為負荷兩端電壓相量;為流過負荷電流相量。將非線性電力網絡的廣義戴維南動態等值阻抗定義為,即

由式(4)可得,電壓實部和虛部分別對電流實部和虛部求偏導均可得到RTHEV和XTHEV。

負荷節點注入的復功率為

式中,PLD、QLD分別為負荷有功和無功功率。

考慮負荷的功率因數是變化的,則限制條件為

式中:k為負荷功率因數;C為常數。當C=0時,負荷的功率因數是恒定的;當k=0時,負荷的無功功率也是恒定的。

通過構造拉格朗日函數來求取PLD的極大值,其中γ為拉格朗日乘數,即

將式(4)和式(5)代入式(7),得到關于Ix、Iy的函數,分別對Ix、Iy求偏導,并令其等于零可得

又因為電壓電流還存在關系

將式(4)和式(9)代入式(8)可得到關于Ix、Iy為變量的二元齊次線性方程組。因為電流不可能恒為零,所以系數矩陣行列式必為零,由此得到負荷的等值阻抗的模值等于系統的廣義戴維南動態等值阻抗的模值,即

這樣證明了負荷功率因素變化的非線性電力網絡輸送有功功率達到極值必要條件。

1.2 阻抗模裕度的定義

通過SNBP判斷電壓穩定臨界崩潰點,重點研究電壓穩定最薄弱節點的SNBP,參考文獻[2]給出了阻抗模裕度的定義為

該指標反映系統正常工作點距離電壓穩定崩潰點有“多遠”,且指標越小電壓穩定性越差。因此由該指標可以判斷電壓穩定性最薄弱節點,然后去單獨求取該薄弱節點的SNBP。當ΔZ>0,節點電壓穩定;當ΔZ=0,電壓臨界穩定,發生鞍結分岔;當ΔZ<0,節點電壓不穩定。

2 變步長的CPF法

2.1 CPF法基本思想

CPF法的基礎是由狀態變量和可變參數構成的潮流方程,即

式中:x為狀態變量,表示電壓幅值及相角;λ為控制功率變化的參數,稱為負荷因子。

將式(12)更具體表述為

式中:PG,i、QG,i為發電機的輸入功率;PL,i、QL,i為負荷消耗的功率;Ps,i、Qs,i為節點的注入功率;g1為PV、PQ節點的并集;g2為PQ節點的集合。若不是發電機節點,則QG,i=0。當發電機節點無功越限,才有無功平衡方程,QG,i=QGmax或QGmin。上下限分別對應最大負荷與基態負荷。

CPF法有4個主要環節:參數化、預估、步長控制、校正[8]。從基態點(x0,λ0)開始,逐步增加負荷因子,求解下一個準確解,最終求得SNBP,具體計算如下。

1)參數化

為了防止在臨近鼻形曲線鼻尖時系數矩陣奇異,以狀態變量的某個分量為參變量,增加一個方程,通常采用局部參數法為

式中:xk為x的第k個分量;η為xk的預測值。

2)預估步

預估步通常采用切線法,即負荷功率控制參變量沿著切線方向增長,逐步預測下一個潮流解[9]。首先確定切線方向為

式中:dx、dλ表示切線方向;±1表示崩潰點前為正,之后為負;為潮流函數對狀態變量和負荷因子的導數;表示與選取狀態變量對應的第k個元素為1,其余為0。再預測下一個潮流點為

式中:λk為當前狀態量;σ為預估的步長。

3)控制步長

常規CPF法的定步長是有極大缺陷的,步長過大能提高計算速度,但極易越過極限點造成潮流發散,步長過小則效率過低[10]。而自動地變步長,既能大大提高運算效率,又不會越限。

4)校正步

2.2 廣義戴維南動態等值預測步長

隨著負荷因子增長,對廣義戴維南動態等值阻抗模值與負荷靜態等值阻抗模值大小及變化趨勢作簡要分析,以圖2所示的簡單兩節點系統為例,線路只計電抗X,超前的相位為δ,δ為功角。

圖2 兩節點系統結構Fig.2 Structure of two-node system

負荷節點的電壓表示為

兩邊同時取負號,并對I求導得

若式(19)為定義的廣義戴維南動態等值阻抗,則可將式(18)改寫為

式中,K為式(19)中廣義戴維南動態等值阻抗的模值,K越來越大,則越來越大。

以上分析表明,隨著負荷因子λ逐漸增大,負荷靜態等值阻抗的模值越來越小,負荷節點的廣義戴維南動態等值阻抗的模值越來越大。可見在基態運行點相差最大,越接近SNBP,差值越小?;鶓B運行時系統穩定,阻抗模裕度應大于0,由式(11)可得,為了使預測到的負荷靜態等值阻抗模能快速地向廣義戴維南動態等值阻抗??拷囱杆傧騍NBP靠近,令預測到的負荷靜態等值阻抗模為上個運行狀態兩者的中值,即

由式(10)知,負荷節點獲取極限功率必要條件是靜態模與動態模相等,則由式(23)知,預測到的運行狀態點既向極限點靠近,而不會越過極限點。

假設采用同步負荷功率擾動方式,詳見式(13),負荷的功率因數也是恒定的,所以預測的負荷視在功率的增長倍數與預測的有功功率的增長倍數相等,那么預測到的負荷因子可表示為

式中:Pi,0、Qi,0分別為節點i基態運行的有功功率與無功功率;Si,k+1為預測到的負荷視在功率。

預測節點i的視在功率平方為

預測到的負荷靜態等值阻抗模與上一個運行狀態的靜態等值阻抗模之差為

以上分析可知,在靠近基態點時,負荷的靜態等值阻抗模值與廣義戴維南動態等值阻抗的模值差值比較大,越靠近SNBP,差值越小。由式(26)知,在靠近基態點時,差值較大,預測到的負荷增長因子距離較大;靠近SNBP時,兩者差值較小,預測到的負荷因子距離也較小。以上方案剛好滿足CPF法自動變步長的要求。

將預測到的負荷因子λk+1與當前運行點的負荷因子λk求差值,并除以預測方向分量dλ,可求出預測步長σ,代入式(16)求出預測運行點,再按照CPF法逐步進行計算。

2.3 阻抗模裕度的計算過程

注入節點的功率可表示為

系統潮流方程表示為

式中:W=[Ps,1,Qs,1,…,Ps,n-1,U2(n-1)]T;Un、U分別為平衡節點與其他節點的電壓,且電壓是在直角坐標表示的,U=[e1,f1,…,en-1,fn-1]T。

注入負荷節點i電流方程為

負荷節點i靜態等值阻抗為

式(29)兩邊同對λ求導得

式(28)兩邊同時對λ求導可得

由式(32)可得

式(31)和式(33)通過求導鏈式法則可得廣義戴維南動態等值阻抗為

將式(30)和式(34)代入式(11)可求出阻抗模裕度值。

3 局部曲線擬合技術

若擬合點距離SNBP較遠,則經過曲線擬合得到SNBP值與準確SNBP值存在很大誤差[12]。一般擬合法至少需要兩個運行點,本文采用一種局部曲線擬合法,通過廣義戴維南動態等值計算的阻抗模來預測CPF的負荷增長因子,實現自動變步長,從而迅速地向SNBP靠近(見第2節)。當預測的運行點潮流收斂,且滿足阻抗模裕度判據ΔZ≤ε(ε為很小的正數,保證此時運行點離SNBP很近),用該運行點就可擬合出SNBP,使擬合精度大大提高。

局部曲線擬合法的數學表達式為

式中:ΔUw、Δλ分別為電壓穩定性薄弱節點滿足要求的擬合點到SNBP的電壓大小變化量和增長的負荷因子(步長);α、β為擬合系數。

對式(35)兩邊同時對Δλ求導可得

因為λ=λ′+Δλ,λ′為當前狀態的負荷因子,Δλ為增長負荷因子,λ為待求的負荷因子,所以dΔUw/dΔλ=dUw/dλ,代入式(36)得

對式(37)繼續對Δλ求導可得

擬合點離SNBP很近,ΔUw=0,dUw/dλ通過式(33)可得,式(32)繼續對λ求導得

dJ/dλ通過J中的電壓對λ求導即可求取,由式(39)就可求出d2Uw/dλ2。

將ΔUw、dUw/dλ、d2Uw/dλ2代入式(36)和式(38)可求得擬合系數α、β。式(35)對ΔUw求導得

在SNBP,dλ/dΔUw=0,所以2αΔUw+β=0,可求出ΔUw。再將α、β、ΔUw代入式(35)可求出擬合點到SNBP的負荷因子增量Δλ,從而求出SN?BP值。

4 算法步驟

步驟1 初始化,輸入原始數據,設置常數ε,并置負荷因子λ0=1,用牛頓法求出基態點潮流數據。

步驟2 潮流若收斂,則按第2.3節步驟計算出各PQ節點的負荷的靜態等值阻抗模和廣義戴維南動態等值阻抗模并計算出阻抗模裕度值ΔZ,將ΔZ最小節點確定為電壓穩定最薄弱節點,假設該節點為j;若潮流不收斂,則返回步驟1。

步驟3 單獨針對節點j,使用局部參數法參數化,并采用切線法得出步長預測的方向dx、dλ。

步驟4 采用廣義戴維南動態等值原理,按照第2.2節方法,預測出下一運行點的負荷因子λ1,并按式(λ1-λ0)/dλ預測步長σ,并計算預測的下一運行點通過垂直校正法求出準確解(x1,λ1)。若不收斂,則采用局部校正法進行校正。

步驟6 采用第3節的局部曲線擬合方法,求出該運行點到SNBP的負荷增長因子Δλcr,繼而求出極限負荷因子λcr=λ+Δλcr。因為擬合點離SN?BP很近,到SNBP的負荷增長因子很小,因此一般不會出現無功越限的情形,該預測點便是SNBP。

5 仿真計算與分析

采用同步功率擾動方式,若節點i為PQ節點,則在直角坐標下的潮流方程為

式中:Pi,0、Qi,0分別為基態下節點i注入的有功與無功;n為潮流計算次數。若節點i為PV節點,則潮流方程表示為

式中:Us,i為節點給出的電壓幅值;ei、fi分別為要求的實際節點電壓實部與虛部。

對多個IEEE系統(IEEE14、IEEE30、IEEE57、IEEE118)進行仿真計算,并考慮無功越限,若無功越限,則將PV節點改成PQ節點。

在基態(λ0=1)運行時,計算4個系統各PQ節點的阻抗模裕度值,并比較大小,確定阻抗模裕度最小的節點為系統靜態電壓穩定性最薄弱節點。4個系統薄弱節點分別為14、30、31、41,阻抗模裕度分別為0.723 1、0.688 0、0.625 4、0.787 0,重點是求這些薄弱節點的SNBP。

先以IEEE30和IEEE57節點系統為例,詳細計算結果見表1和表2。其中λ0為當前運行點的負荷因子,λ1為預測的下一個運行點的負荷因子。

由表1可知,n在1~6時是自動變步長的連續潮流計算過程,可見在基態點附近(1~1.5),預測到的相鄰兩個運行點的負荷因子之差較大,即步長較大。隨著負荷因子逐漸向SNBP(1.5以后)靠近,差值越來越小,即步長越來越小,滿足理想的步長控制要求。阻抗模裕度隨著負荷因子的增大而減小,當n=6時,ΔZ<ε(ε=0.03 p.u.),該運行點離SNBP很近,滿足局部曲線擬合的要求。n=7是擬合的結果,擬合點負荷因子λ=1.546 34,從擬合點到SN?BP增長負荷因子很小Δλcr=0.000 25,擬合后無PV節點無功越界,計算結束。只用了7次潮流計算就求取了SNBP的值。

表1 IEEE30節點系統計算結果Tab.1 Calculation results of IEEE 30-node system

表2 IEEE57節點系統計算結果Tab.2 Calculation results of IEEE 57-node system

由表2可知,對于節點數較多的IEEE57節點系統,前6次潮流計算也是向SNBP快速逼近的過程,也滿足連續潮流計算的步長控制要求。當n=6時,ΔZ<ε(ε=0.03 p.u.),滿足擬合要求,擬合后也無PV節點無功越界,只用了7次潮流計算就求取了SNBP值??梢姡嬎愦螖挡粫S著系統規模的擴大而明顯增加。

將所提方法和文獻[13]提出的改進CPF法以及固定步長(取0.01)的常規CPF法計算結果進行比較,比較結果見表3(含IEEE14、IEEE118的結果)。其中NB、λB分別為改進CPF法的計算次數與SN?BP的值;NG、λG分別為所提方法總共計算次數與SNBP的值;NC、λC分別為常規的定步長連續潮流法計算次數與SNBP的值。

表3 3種方法的比較Tab.3 Comparison among three methods

由表3可知,本文所提方法求取的SNBP值與改進CPF法計算值最大相對誤差為0.028%,與定步長常規CPF法相比,最大相對誤差為0.028%。但本文所提方法的潮流計算次數比另外兩種方法少很多,并且隨著系統節點數的增多,計算次數變化不明顯,計算穩定,而另外兩種方法卻不具備這個特點??梢?,本文所提方法不僅計算精度高,而且計算量小,計算速度快,適用于多節點的電力網絡。

6 結語

運用廣義戴維南動態等值方法,提出了一種自動變步長的CPF法,預測步長在基態點附近的概率大,而在崩潰點附近的概率小,能快速地逼近而不越過SNBP。由于阻抗模裕度反映當前運行點到SNBP距離,提出了阻抗模裕度判據,由阻抗摸裕度判據判定滿足擬合要求的運行點,進行精確局部擬合,可以求取較準確的SNBP。算法簡單,只用到雅可比矩陣保留的因子表,多個IEEE測試系統的仿真表明,該法計算精度高,次數少,速度快,適用于規模較大的電力網絡,具有一定的在線應用價值。

[1]李娟,黃喜旺,關程宇,等(Li Juan,Huang Xiwang,Guan Chengyu,et al).SVC和IPC聯合改善異步機風電場的電壓穩定性(Improvement of voltage stability of asynchro?nous wind farms based on SVC and IPC)[J].電力系統及其自動化學報(Proceedings of the CSU-EPSA),2016,28(8):79-84.

[2]吳政球,李日波,鐘浩,等(Wu Zhengqiu,Li Ribo,Zhong Hao,et al).電力系統靜態電壓穩定極限及裕度計算綜述(Summary of power system's static voltage stability lim?itation and load margin calculation)[J].電力系統及其自動化學報(Proceedings of the CSU-EPSA),2010,22(1):126-132.

[3]牟曉明,李志民(Mou Xiaoming,Li Zhimin).一種計算電壓穩定邊界的兩階段潮流方法(Two-stage power flow method for calculation voltage stability margin)[J].電力自動化設備(Electric Power Automation Equipment),2013,33(3):72-77.

[4]楊小煜,周孝信(Yang Xiaoyu,Zhou Xiaoxin).基于極小擴張系統方法的靜態電壓穩定臨界點計算(Calcula?tion of the critical points of static voltage stability with minimally extended system method)[J].中國電機工程學報(Proceedings of the CSEE),2009,29(25):32-36.

[5]Irisarri G D,Wang X,Tong J,et al.Maximum loadability of power systems using interior point non-linear optimiza?tion method[J].IEEE Trans on Power Systems,1997,12(1):162-172.

[6]趙晉泉,張伯明(Zhao Jinquan,Zhang Boming).改進連續潮流計算魯棒性的策略研究(A study on the strategy for improving robustness of continuation power flow com?putation)[J].中國電機工程學報(Proceedings of the CSEE),2005,25(22):7-11.

[7]湯涌.電力系統電壓穩定性分析[M].北京:科學出版社,2011.

[8]張伯明,陳壽孫,嚴正.高等電力網絡分析[M].北京:清華大學出版社,2007.

[9]Ajjarapu V,Christy C.The continuation power flow:A tool for steady state voltage stability analysis[J].IEEE Trans on Power Systems,1992,7(1):416-423.

[10]張堯,張建設,袁世強(Zhang Yao,Zhang Jianshe,Yuan Shiqiang).求取靜態電壓穩定極限的改進連續潮流法(Improved continuation power flow algorithm for obtain?ing the limit of static voltage stability)[J].電力系統及其自動化學報(Proceedings of the CSU-EPSA),2005,17(2):21-25.

[11]王成山,魏煒(Wang Chengshan,Wei Wei).一種改進的步長控制連續性潮流計算方法(An improved continua?tion method with controlled step size)[J].電工技術學報(Transactions of China Electrotechnical Society),2004,19(2):59-63.

[12]Zhao Jinquan,Chiang H D.Enhanced look-ahead load margin estimation for voltage security assessment[J].IEEE Trans on Power Systems,2003,6(3):2640-2645.

[13]Zarate L A L,Castro C A.Fast computation of security margins to voltage collapse based on sensitivity analysis [J].IEE Proceedings-Generation,Transmission and Distri?bution,2006,153(1):35-43.

Fast Calculation Method for Saddle Node Bifurcation Point of Voltage Stability of Power System

ZHU Yongqiang1,2,LIU Guangye1,CAO Bin3,KUANG Yilei1,KANG Zhihao1,LIU Yuanyuan1
(1.College of Electrical and Information Engineering,Hunan University,Changsha 410082,China;2.Yueyang Power Supply Company,State Grid Hunan Electric Power Company,Yueyang 414000,China;3.Maintenance Company,State Grid Hunan Electric Power Company,Changsha 410000,China)

TM712

A

1003-8930(2017)09-0086-07

10.3969/j.issn.1003-8930.2017.09.015

2015-09-23;

2017-01-11

國家自然科學基金資助項目(51577053)

朱永強(1989—),男,碩士,助理工程師,研究方向為電力系統電壓穩定分析與控制。Email:1090321338@qq.com

劉光曄(1960—),男,博士,教授,博士生導師,研究方向為電力系統分析與控制、輸變電技術、鐵路牽引供電系統、電力系統繼電保護。Email:liuguangye@21cn.com

曹 斌(1991—),男,碩士研究生,研究方向為高電壓與絕緣技術。Email:791589645@qq.com

主站蜘蛛池模板: 无码综合天天久久综合网| 久久这里只有精品8| 国产精品久久精品| 亚洲欧美不卡| 99在线视频精品| 久久综合激情网| 免费人成网站在线观看欧美| 欧美一级黄片一区2区| 亚欧成人无码AV在线播放| 欧洲成人在线观看| 激情六月丁香婷婷四房播| 99re视频在线| 中文字幕乱妇无码AV在线| 亚洲国产黄色| 一区二区三区在线不卡免费| 欧美黄网站免费观看| 精品无码国产自产野外拍在线| 欧美三级日韩三级| 国产在线自在拍91精品黑人| 精品国产女同疯狂摩擦2| 全午夜免费一级毛片| 秋霞国产在线| 中文字幕无码电影| 精品欧美视频| 四虎永久免费地址| 亚洲va在线观看| 国产女人在线| 伊人网址在线| 最新亚洲人成网站在线观看| 国产第一页第二页| 夜夜拍夜夜爽| 中文字幕人成乱码熟女免费| 99r在线精品视频在线播放| 一级毛片视频免费| 狼友视频国产精品首页| 亚洲中文在线视频| 在线高清亚洲精品二区| 国产不卡国语在线| 九九久久精品免费观看| 拍国产真实乱人偷精品| 欧美www在线观看| 欧美成人精品在线| 国产精品嫩草影院av| 欧美色图第一页| 免费无码AV片在线观看中文| 成人福利在线视频| 亚洲精品色AV无码看| 免费 国产 无码久久久| 三上悠亚精品二区在线观看| 日韩精品一区二区三区中文无码 | 免费高清毛片| 黄色网址免费在线| 免费毛片全部不收费的| 国产人成乱码视频免费观看| 亚洲第一天堂无码专区| 丁香五月亚洲综合在线| 真实国产精品vr专区| 在线精品自拍| 91精品伊人久久大香线蕉| 欧美国产视频| 久久综合丝袜日本网| 波多野结衣视频一区二区| 欧美一级黄色影院| 成人国产精品2021| 色窝窝免费一区二区三区| 亚洲视频免| 日韩精品欧美国产在线| 一区二区三区高清视频国产女人| 国产午夜看片| 国产白浆一区二区三区视频在线| 视频一区亚洲| 国产99视频精品免费视频7| 欧美性色综合网| 国产超碰在线观看| 国产精品久久自在自2021| 国产乱码精品一区二区三区中文 | 国产精品手机视频一区二区| 9cao视频精品| 四虎永久在线| 国产成人福利在线视老湿机| 大陆国产精品视频| 国产成人成人一区二区|