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

雙繞組電力變壓器梯形網絡頻率響應的矩陣算法

2018-02-27 10:58:02任富強汲勝昌祝令瑜劉勇楊帆陸偉鋒李熙寧
西安交通大學學報 2018年2期
關鍵詞:變壓器

任富強,汲勝昌,祝令瑜,劉勇,楊帆,陸偉鋒,李熙寧

(西安交通大學電力設備電氣絕緣國家重點實驗室,710049,西安)

電力變壓器安全可靠運行對于保障電力系統的持續運行意義重大[1-2]。變壓器繞組是發生故障損壞的主要部位,并且故障程度隨著短路電流水平的提高而愈發嚴重[3-4]。作為一種離線無損檢測變壓器繞組變形的方法,頻率響應法由于具有良好的測試重復性和靈敏度得到了廣泛應用[5]。頻率響應法一般通過計算不同曲線之間的相關系數來大致判斷繞組變形程度[6]。為了進一步挖掘頻響曲線蘊含的繞組狀態信息,更加準確診斷繞組變形[7-8],需要研究不同變形形式下繞組高頻等效電路中元件參數的變化,進而得到其頻響曲線的變化規律[9]。上述方法的關鍵點之一是電力變壓器高頻等效網絡的仿真計算[10],這也是建立變壓器繞組變形與頻響曲線內在聯系的重要步驟。因而,研究電力變壓器高頻等效模型頻率響應的計算方法具有重要意義。

電力變壓器繞組高頻等效模型包括黑箱模型、物理模型以及混合模型3類,其中物理模型能準確描述高頻激勵下繞組的電磁特性,應用最為廣泛[11-12]。等效梯形網絡模型為國內外學者普遍使用的物理模型,考慮了繞組單元的自感互感關系以及分布電容等參量[13-14],能夠對電力變壓器繞組頻響特性進行精確的仿真研究。計算梯形網絡的頻響特性一般基于電路仿真軟件。文獻[15]忽略繞組二次側影響,建立了單繞組電路模型,仿真獲得變壓器繞組的頻響特性,對等效電路電感、分布電容等參數改變時頻響曲線的變化規律進行了研究。文獻[16]通過忽略單元間互感來簡化計算,提取了繞組發生軸心偏移、輻向變形和軸向位移等故障時等效電路的參數,進而對不同故障狀態下的繞組頻響特性進行比較。文獻[8]對等效電路中不同元件在不同頻段的靈敏度進行了分析,仿真獲得了繞組發生預緊力松動、軸向及輻向位移等故障時頻響曲線在不同頻段的變化規律。上述文獻均基于電路仿真軟件,對變壓器繞組等效梯形網絡的頻響特性進行了研究。但是,通過仿真軟件建立電路模型比較復雜,并且由于忽略了繞組二次側以及單元間互感以簡化計算,等效電路模型并不完整,有一定的計算誤差。

利用基于電路原理的算法也可以實現等效電路頻響特性的計算。文獻[10]將等效梯形網絡的阻抗矩陣(電阻和電感串聯支路)進行轉化,結合電路導納矩陣(分布電容支路),構建整體網絡的矩陣形式,得到關于網絡節點電壓的方程組并對其求解,但是算法多次對矩陣求逆,導致計算誤差較大,且求解效率較低。文獻[17]將等效電路的阻抗支路電流和節點電壓作為未知量,通過阻抗或導納矩陣來聯系不同支路或節點,實現支路電流和節點電壓的求解。由于計算頻響特性僅利用了節點電壓數據,故上述算法增加了未知量,使得計算效率下降。文獻[18]中的矩陣求解算法和文獻[10]類似,通過網絡導納矩陣的轉化,實現節點電壓及部分支路電流的求解,但是算法也存在多變量的求解問題,而且沒有考慮繞組二次側影響,此外對方程組常數項求取的關鍵步驟也缺乏必要推導。綜合上述文獻,基于電路原理算法構建的等效梯形網絡模型更加完整,頻率響應的求解過程更加系統化且求解效率很高。但是,上述文獻對算法的推導均不夠系統,沒有解決網絡矩陣構建以及多變量求解或者矩陣多次求逆導致的低計算效率的問題。

鑒于上述問題,針對考慮變壓器繞組二次側的等效梯形網絡,本文首先構建用于描述網絡的阻抗支路矩陣及導納支路矩陣;其次,基于電路原理,獲得阻抗支路電流、導納支路電流與節點電壓的關系;最后,通過節點電壓列向量的求取,獲得電路的頻率響應。在此基礎上,提出了一種不同于前述文獻的求取頻率響應的矩陣算法。通過與相關文獻的計算結果對比,證明本文提出的算法在保證模型完整性和精確性的同時,具有良好的計算精度及計算效率。

1 矩陣算法推導

1.1 雙繞組電力變壓器等效梯形網絡模型

圖1為雙繞組電力變壓器的等效梯形網絡模型,由電阻、電感、電容及電導等元件組成,分別表征繞組及鐵心的有功損耗、電磁效應、介質損耗及電容效應。元件具體數值可以通過解析法或有限元仿真計算,也可以通過智能算法結合相關的頻響數據確定[11,17,19-20]。

在圖1所示的電路中,基本梯形單元由電阻電感串聯支路(阻抗支路)及電導電容并聯支路(導納支路)組成。基本梯形單元用數匝層式或餅式繞組線圈來等效,線匝數根據最高激勵頻率以及繞組的尺寸和結構等參數確定,并保證流過此單元線匝的電流近似恒定[15]。在頻響法的掃頻范圍(1~1 000 kHz)內,對于尺寸不同或結構相異的繞組,其梯形等效電路拓撲結構一致,區別僅在于基本梯形單元數的不同或單元參數計算方法的差異,故梯形等效電路適用于描述多數雙繞組電力變壓器的高頻電磁特性,因而基于此等效電路建立的矩陣算法也普遍適用于一般結構的雙繞組變壓器。另外,對于同一繞組,健康或故障狀態僅對應為等效網絡相關單元元件參數的不同,而其拓撲結構并不改變,因而本文得到的矩陣算法對健康或故障狀態繞組均適用。

圖1 雙繞組電力變壓器等效梯形網絡

設圖1中高、低壓繞組的梯形單元數均為(n-1),則等效電路是由2(n-1)個單元組成的網絡。網絡中需要求解的節點電壓數為2n,用未知量U1,U2,…,U2n表示。其次,阻抗支路電流是未知的,其數量為2(n-1),用未知量I1,I2,…,I2(n-1)表示,電流正方向規定為由較小編號節點指向較大編號節點的方向。按照關聯的節點數目統計未知的導納支路電流,其數量為2n,用未知量i1,i2,…,i2n表示,電流正方向規定為流出其關聯節點的方向。

1.2 阻抗支路與導納支路的矩陣描述

求解等效梯形網絡模型的頻率響應,首先需要構建用于描述網絡單元參數的矩陣。通過將圖1中復雜網絡分解成阻抗支路與導納支路并求取其矩陣,即可完整描述圖1所示網絡。

阻抗支路由電阻及電感串聯構成,支路間通過阻抗支路電流相聯系,故定義(2n-2)階方陣MRL來描述頻域中的阻抗支路,其形式為

MRL=

(1)

式中:s為拉普拉斯算子;Li,i(i=1,2,…,2n-2)為第i條阻抗支路的電感值;Mi,j(i,j=1,2,…,2n-2,i≠j)代表編號為i、j的阻抗支路互感;Ri,j(i,j=1,2,…,2n-2)為第i條阻抗支路的電阻(當i=j時)或第i、j條支路的互阻(當i≠j時)[20]。電感及電阻在網絡中的具體位置參見圖1。顯然,在矩陣MRL中,對角元素MRL(i,i)(i=1,2,…,2n-2)為第i條阻抗支路的自阻抗,元素MRL(i,j)(i,j=1,2,…,2n-2,i≠j)則代表第i、j條阻抗支路通過互阻抗相互關聯。MRL即為阻抗支路的矩陣描述。

導納支路為分布電導與分布電容并聯支路的等效,支路間通過節點電壓相聯系。為了在頻域中描述導納支路,先定義下式

(2)

式中:i、j為節點編號;Ggi與Cgi分別代表第i個節點的對地分布電導與分布電容;Gi,j與Ci,j分別代表編號為i、j(i,j=1,2,…,2n)的節點間分布電導與分布電容。對于同一繞組上的節點,Gi,j與Ci,j分別為繞組縱向電導與縱向電容;對于相異繞組間的節點,Gi,j與Ci,j分別代表高低壓繞組間的分布電導與分布電容。繞組對地電導、電容,縱向電導、電容及高低壓繞組間分布電導、電容的計算可以參見文獻[17]。電導及電容在網絡中的具體位置參見圖1。

類似地,定義2n階方陣MCG,其形式如下

MCG=

(3)

式中:元素MCG(i,j)代表節點i通過導納與節點j的電氣聯系。矩陣MCG的對角元素為自導納,其值為正;非對角元素為互導納,其值為負。MCG即為導納支路的矩陣描述。式(1)與式(3)即完整描述了圖1所示的等效電路網絡。

1.3 阻抗支路與導納支路電流的關系

根據基爾霍夫電流定律(KCL),可得節點的KCL方程為

(4)

定義n×(n-1)階矩陣A,其元素aij(i=1,2,…,n且j=1,2,…,n-1)的具體值為

(5)

根據式(4)和式(5),可以得到如下關系

(6)

式中:O為n×(n-1)階零矩陣;B為2n×(2n-2)階矩陣。式(6)即為阻抗支路電流與導納支路電流的關系。

1.4 節點電壓與阻抗支路電流的關系

基于阻抗支路的電壓電流關系(VCR),可得

(7)

結合式(1)、式(5)與式(7),通過觀察,可以得到節點電壓與阻抗支路電流關系的最終形式

(8)

1.5 節點電壓與導納支路電流的關系

基于導納支路的VCR,即可得到節點電壓與導納支路電流的關系,其形式如下

(9)

結合式(2)、式(3)與式(9),即可得到節點電壓與導納支路電流的最終形式

(10)

式(10)即為節點電壓與導納支路電流的聯系。

1.6 節點電壓列向量的求解

式(6)、式(8)與式(10)建立了阻抗支路電流、導納支路電流與節點電壓3者之間的聯系。對于頻率響應的求解僅需節點電壓值,因而在上述3式中,消去阻抗支路及導納支路電流列向量,即得到關于節點電壓的方程組

(11)

式中:0為2n維0值列向量;M為2n階方陣,其元素M(i,j)代表編號為i的節點電壓與編號為j的節點電壓間的電氣聯系。

直接求解式(11)是不可行的。由于U(1,1)為進行頻率響應測試時施加的激勵,即向量U的第1個元素為已知量,故需要將式(11)進行進一步變換。

取列向量U刪除第1個元素U(1,1)得到的(2n-1)維未知列向量為U′。將矩陣M中表征其他節點電壓與U(1,1)關聯的第1行元素刪除得到矩陣M1。取M1的第1列構成(2n-1)維列向量N1,取M1的剩余部分作為(2n-1)階方陣N。經過上述變換后即可得到下式

NU′=-U(1,1)N1

(12)

通過求解非齊次線性方程組(12),即可解得U′,即(2n-1)個未知節點電壓。最后,通過下式得到等效電路的頻率響應

(13)

式中:U′(n,1)表示列向量U′的第n個元素,對應圖1中編號為(n+1)節點的電壓值,即采樣電阻的電壓。

2 繞組頻率響應的計算

本節基于電路仿真軟件PSPICE(Version 16.5)與Simulink(MATLAB Version 8.5),利用提出的矩陣算法來計算變壓器繞組等效梯形網絡的頻率響應,驗證矩陣算法的準確性以及計算的高效性,包含變壓器單繞組模型、考慮單元互感的單繞組模型以及雙繞組模型共3個計算案例。

2.1 算例模型及參數

算例1繞組梯形等效電路由7個單元組成,參數來源于文獻[21],各單元的參數相等,其自感、縱向電容、對地電容以及采樣電阻值分別為41.426 mH、19.8 pF、1 214.286 pF以及50 Ω。文獻忽略了對頻率響應測試結果影響較小的開路二次側及數值較小的繞組電阻以簡化計算。綜上,繞組梯形網絡簡化為如圖2所示的等效電路。對于簡化后的單繞組模型,只需要將雙繞組梯形網絡矩陣算法中表征二次側梯形單元的參數刪除即可。

圖2 算例1單繞組7梯形單元模型

算例2來源于文獻[18],為包含10個梯形單元的單繞組變壓器模型,對比算例1,其仿真單元數目增加并計及了單元間互感,等效電路與圖2類似。單元的電阻、自感、縱向電容、對地電容以及采樣電阻值分別為0.137 Ω、10.32 mH、22.54 pF、200 pF以及50 Ω。單元間互感有9個,其值分別為5.815、2.978、1.615、0.912、0.551、0.231、0.056、0.024及0.011 mH。

算例3為雙繞組等效電路的計算,來源于文獻[9],共包含20個梯形單元,文獻忽略了單元間互感耦合作用。其等效電路圖與圖1類似,相關參數如表1所示。表中的物理量Ri、Li、Ki、Cgj、Chlj與Gj(i=1,2,…,10且j=1,2,…,11)分別表示單元的電阻、電感、縱向電容、對地電容、高低壓繞組間的分布電容與電導。

表1 算例3雙繞組20梯形單元模型相關參數

2.2 算例計算結果

采用電路仿真軟件建立的仿真電路與圖2類似,在模型的節點1施加幅值為1的掃頻電源。對于掃頻模式的實現,在PSPICE中可以將仿真模式設為AC Sweep/Noise,在Simulink中則通過編程來實現。通過繞組末端的采樣電阻得到響應端的正弦波幅值U0,最后通過下式計算得到電路的頻率響應

H=20logU0

(14)

利用3種方法計算得到的3個算例中繞組等效梯形網絡的頻率響應匯總分別如圖3~圖5所示。通過矩陣算法計算得到考慮單元間互感的頻率響應曲線見圖4。

圖3 算例1頻率響應的計算結果

圖4 算例2頻率響應的計算結果

2.3 計算結果分析

由2.2節中的計算結果可知,對于上述3個算例,除極少數頻段(高幅值諧振峰對應頻段附近)計算結果有差別外,通過矩陣算法以及軟件仿真得到的計算結果完全一致。為了量化不同計算方法的結果差異,分別計算矩陣算法與兩種仿真得到的頻響曲線的相關系數,如表2所示,表中全頻、低頻、中頻及高頻表示的頻率范圍為1~1 000 kHz,1~100 kHz,100~600 kHz,600~1 000 kHz。組1為矩陣算法結果和PSPICE仿真結果的相關系數;組2為矩陣算法結果和Simulink仿真結果的相關系數。

表2 3個算例頻響曲線的相關系數

圖5 算例3頻率響應的計算結果

通過相關系數判斷,對于算例1,3條曲線的相似度很高,由圖3可知,其差異在中頻段的最大諧振峰,即176 kHz處幅值約為-600 dB的諧振峰附近,而其他頻率點完全重合;算例2與算例1類似,差異集中于330 kHz處的最大諧振峰附近;對于算例3,3種計算方法的結果完全一致。算例1與算例2中的差異是由仿真軟件的設定以及數據在不同軟件之間傳遞產生的截斷誤差所致。在最大諧振峰處,等效電路會發生電流諧振,此時繞組中的電流最小,采樣電阻得到的電壓最低。在算例1中,-600 dB諧振峰對應的采樣電阻電壓僅為10-30,算例2中對應電壓更低。這些極低數值小于電路仿真軟件的最小分辨值,故計算結果出現誤差。在算例3中,由于頻率響應的增益較小(最大峰值僅約為-250 dB),計算結果完全吻合。由此可知,矩陣算法的計算結果準確,即使在頻率響應曲線的幅值差異較大時,依然保持很高的精度。

在同一臺計算機上,利用相同的MATLAB軟件執行矩陣算法和Simulink仿真以及運行PSPICE程序的計算耗時匯總于表3。3者的計算耗時可以通過編程或者訪問軟件日志文件獲得。由表3可知,矩陣算法的計算效率與PSPICE軟件相當。由于M文件與Simulink程序之間頻繁的數據傳遞等原因,Simulink仿真計算耗時較長。綜合比較3種方法,認為矩陣算法的計算效率很高。

表3 3個算例的計算耗時比較

對于算例2,在電路仿真軟件中,當梯形單元數較多時,實現在所有單元間引入互感的難度較大。但是,基于矩陣算法,可以通過式(1)方便引入單元間互感,此時的頻響曲線如圖4所示,可見考慮互感與否得到的頻響曲線差別較大。由于互感是繞組單元電磁耦合的重要參數,引入互感使等效電路更加符合繞組實際狀態,可保證建模的完整性,故在仿真時必須將單元間互感考慮在內[18]。由圖4可知,在考慮單元間互感時,頻率響應曲線諧振點有向高頻偏移的趨勢,且曲線幅值在下降。頻率響應曲線幅值下降意味著采樣電阻電壓升高,此時的計算結果會更加準確。因而,通過矩陣算法可以引入互感,可保證仿真模型的完整性,使等效模型更加符合繞組實際狀態,得到的仿真結果也更加精確。

綜上,基于矩陣算法計算變壓器繞組等效梯形網絡模型的頻率響應,具有很高的計算精度以及計算效率,并且能夠保證模型的完整性,使等效模型更加符合繞組的實際狀態。

3 結 論

本文針對雙繞組電力變壓器等效梯形網絡模型,首先構建描述網絡的阻抗支路和導納支路矩陣,其次基于電路原理,獲得支路電流和節點電壓間的關系,最后通過矩陣運算求解節點電壓列向量,從而得到等效電路的頻率響應。通過矩陣算法與電路仿真軟件計算結果相對比,證明了本文矩陣算法的計算精確性和高效性。另外,在矩陣算法中,可以通過更改相應的矩陣將單元間互感考慮在內,使計算模型更加精確完整。綜上所述,通過矩陣算法得到變壓器繞組的頻率響應是兼顧計算精度、計算效率以及保證計算模型完整性的最優方法。

[1] WANG M, VANDERMAAR A J, SRIVASTAVA K D. Review of condition assessment of power transformers in service [J]. IEEE Electrical Insulation Magazine, 2002, 18(6): 12-25.

[2] 劉勇, 楊帆, 張凡, 等. 檢測電力變壓器繞組變形的掃頻阻抗法研究 [J]. 中國電機工程學報, 2015, 35(17): 4505-4516. LIU Yong, YANG Fan, ZHANG Fan, et al. Study on sweep frequency impedance to detect winding deformation within power transformer [J]. Proceedings of the CSEE, 2015, 35(17): 4505-4516.

[3] 朱葉葉, 汲勝昌, 張凡, 等. 電力變壓器振動產生機理及影響因素研究 [J]. 西安交通大學學報, 2015, 49(6): 115-125. ZHU Yeye, JI Shengchang, ZHANG Fan, et al. Vibration mechanism and influence factors in power transformers [J]. Journal of Xi’an Jiaotong University, 2015, 49(6): 115-125.

[4] 張凡, 汲勝昌, 祝令瑜, 等. 短路沖擊下變壓器振動頻響函數研究 [J]. 西安交通大學學報, 2017, 51(2): 97-103, 154. ZHANG Fan, JI Shengchang, ZHU Lingyu, et al. Frequency response function of short circuit viberation for power transformer [J]. Journal of Xi’an Jiaotong University, 2017, 51(2): 97-103, 154.

[5] REYKHERDT A A, DAVYDOV V. Case studies of factors influencing frequency response analysis measurements and power transformer diagnostics [J]. IEEE Electrical Insulation Magazine, 2011, 27(1): 22-30.

[6] RYDER S A. Diagnosing transformer faults using frequency response analysis [J]. IEEE Electrical Insulation Magazine, 2003, 19(2): 16-22.

[7] KIM J W, PARK B, JEONG S, et al. Fault diagnosis of a power transformer using an improved frequency-response analysis [J]. IEEE Transactions on Power Delivery, 2005, 20(1): 169-178.

[8] 王曙鴻, 張海軍. 一種基于模型修正的變壓器繞組變形定量診斷方法: CN105180792A[P]. 2015-07-01.

[9] ABU-SIADA A, HASHEMNIA N, ISLAM S, et al. Understanding power transformer frequency response analysis signatures [J]. IEEE Electrical Insulation Magazine, 2013, 29(3): 48-56.

[10]BJERKAN E. High frequency modeling of power transformers [D]. Trondheim, Norwegian: Norwegian University of Science and Technology, 2005: 192-195.

[11]SHABESTARY M M, GHANIZADEH A J, GHAREHPETIAN G B, et al. Ladder network parameters determination considering nondominant resonances of the transformer winding [J]. IEEE Transactions on Power Delivery, 2014, 29(1): 108-117.

[12]RASHTCHI V, RAHIMPOUR E, REZAPOUR E M. Using a genetic algorithm for parameter identification of transformer RLCM model [J]. Electrical Engineering, 2006, 88(5): 417-422.

[13]ABEYWICKRAMA N, SERDYUK Y V, GUBANSKI S M. High-frequency modeling of power transformers for use in frequency response analysis (FRA) [J]. IEEE Transactions on Power Delivery, 2008, 23(4): 2042-2049.

[14]SHINTEMIROV A, TANG W H, WU Q H. A hybrid winding model of disc-type power transformers for frequency response analysis [J]. IEEE Transactions on Power Delivery, 2009, 24(2): 730-739.

[15]張震. 變壓器繞組變形檢測裝置的研究 [D]. 沈陽: 沈陽工業大學, 2013: 18-27.

[16]陳曉晗. 基于有限元法的電力變壓器繞組變形檢測與識別的仿真研究 [D]. 重慶: 重慶大學, 2015: 42-55.

[17]ABEYWICKRAMA K B, SERDYUK Y V, GUBANSKI S M. Exploring possibilities for characterization of power transformer insulation by frequency response analysis (FRA) [J]. IEEE Transactions on Power Delivery, 2006, 21(3): 1375-1382.

[18]朱明林. 頻響分析法及其模型辨識檢測和診斷變壓器繞組變形的研究 [D]. 上海: 上海交通大學, 2001: 12-20.

[19]PODOLTSEV A D, NILANGA K G, ABEYWICKRAMA K B, et al. Multiscale computations of parameters of power transformer windings at high frequencies: part I Small-scale level [J]. IEEE Transactions on Magnetics, 2007, 43(11): 3991-3998.

[20]PODOLTSEV A D, ABEYWICKRAMA K B, SERDYUK Y V, et al. Multiscale computations of parameters of power transformer windings at high frequencies: part II Large-scale level [J]. IEEE Transactions on Magnetics, 2007, 43(12): 4076-4082.

[21]武劍利. 頻響分析法檢測變壓器繞組變形的理論研究 [D]. 武漢: 武漢大學, 2004: 21-22.

猜你喜歡
變壓器
10KV變壓器配電安裝柜中常見問題分析
理想變壓器的“三個不變”與“三個變”
超導變壓器技術展望
開關電源中高頻變壓器的設計
變壓器光纖測溫探頭的安裝固定
電子制作(2017年8期)2017-06-05 09:36:15
探究供配電變壓器的推廣運行
電子制作(2017年8期)2017-06-05 09:36:15
變壓器運行中的異常現象與故障處理
電子制作(2017年7期)2017-06-05 09:36:14
變壓器差動保護誤動原因探討
電子制作(2017年10期)2017-04-18 07:23:16
一種不停電更換變壓器的帶電作業法
變壓器免維護吸濕器的開發與應用
主站蜘蛛池模板: 91久久国产综合精品| 区国产精品搜索视频| 国产精品原创不卡在线| 亚洲一区网站| 日本少妇又色又爽又高潮| 日韩高清一区 | 中文毛片无遮挡播放免费| 日韩免费成人| 日韩精品无码免费一区二区三区| 国产va视频| 青青草国产一区二区三区| 国产SUV精品一区二区6| 黄色不卡视频| 亚洲手机在线| 国产精品亚洲片在线va| 重口调教一区二区视频| 中文精品久久久久国产网址| 夜夜拍夜夜爽| 中文无码毛片又爽又刺激| 亚洲精品午夜无码电影网| 色哟哟色院91精品网站| 久久婷婷五月综合97色| 久久亚洲中文字幕精品一区| 制服丝袜亚洲| 国产菊爆视频在线观看| 日韩欧美国产成人| 无码国产偷倩在线播放老年人| 亚洲色图欧美在线| 高清无码不卡视频| 最新无码专区超级碰碰碰| 成年看免费观看视频拍拍| 91丝袜在线观看| 国产精品999在线| 午夜无码一区二区三区| 一区二区在线视频免费观看| 色婷婷在线影院| 国产美女无遮挡免费视频网站 | 午夜精品久久久久久久99热下载| 亚洲Av综合日韩精品久久久| 色综合中文综合网| 欧美日韩一区二区三| 天天做天天爱夜夜爽毛片毛片| 久久免费看片| 国产成人综合亚洲欧美在| 在线观看视频一区二区| 国产激爽爽爽大片在线观看| 黄色三级网站免费| 99久视频| 国产91导航| 久久99国产乱子伦精品免| 看你懂的巨臀中文字幕一区二区 | 91欧美亚洲国产五月天| 国产精品区网红主播在线观看| 在线不卡免费视频| 亚洲欧美一区二区三区麻豆| 鲁鲁鲁爽爽爽在线视频观看| 中文字幕在线不卡视频| 亚洲欧美自拍一区| 国产免费福利网站| 亚洲欧洲国产成人综合不卡| 欧美亚洲第一页| 国产网友愉拍精品视频| 亚洲欧美日韩成人高清在线一区| 91热爆在线| 孕妇高潮太爽了在线观看免费| 99国产精品国产高清一区二区| 99这里只有精品在线| 亚洲Av激情网五月天| 99在线观看国产| 国产免费人成视频网| 国产精品一区二区无码免费看片| 午夜电影在线观看国产1区| 欧美人人干| 国产精品无码AV片在线观看播放| 波多野结衣一区二区三区88| 久久这里只精品国产99热8| 欧美视频免费一区二区三区| 2020精品极品国产色在线观看| 黄色网址免费在线| 免费看一级毛片波多结衣| 国产精品网址你懂的| 亚洲IV视频免费在线光看|