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

基于時步有限元的抽水蓄能電機瞬態參數計算方法的對比

2015-11-25 09:30:28王偉華王紅宇許國瑞康錦萍孫玉田
電工技術學報 2015年1期

王偉華 王紅宇 許國瑞 康錦萍 孫玉田

(1.華北電力大學電氣與電子工程學院 北京 102206 2.哈爾濱大電機研究所 哈爾濱 150040)

1 引言

抽水蓄能電站在電力系統中擔負著調峰、填谷、事故備用等重要任務,可以快速跟蹤負荷,且能夠減少煤的損耗,對環境無污染,在現代電網中具有不可替代的位置[1]。抽水蓄能電機作為抽水蓄能電站的主要設備發揮著至關重要的作用。抽水蓄能電機的瞬態電抗和有關時間常數對電機的瞬態和動態行為有著重要影響[2]。因此,研究抽水蓄能電機的瞬態參數具有重要的現實意義。

對電機瞬態參數的計算最早采用傳統的路的設計公式[3],該公式沒有考慮飽和以及渦流對瞬態參數的影響,已不能滿足現代設計和運行的要求,因此,有必要用場的方法對瞬態參數進行求取。用場的方法確定電機的瞬態參數總體上講有兩種方法:一種是頻域法,另一種是包絡線法。文獻[2,4-6]均采用頻域法計算電機的瞬態參數,該方法由于結合頻率特性,所以要求繞組磁鏈為正弦,而對于轉子為凸極的抽水蓄能電機來說,繞組磁鏈分布會偏離正弦波形。文獻[2,7]采用包絡線法對電機瞬態參數進行計算,該方法采用瞬態電磁場模型,其數學模型直接反映電機的實際工況,計算結果較為準確,但計算量大,且在后處理中涉及包絡線的繪制及曲線擬合等操作,使結果容易產生偏差。

本文基于抽水蓄能電機單元電機的場路耦合時步有限元模型,采用包絡線法得到響水澗抽水蓄能電機在發電及電動工況時的穩態、瞬態、超瞬態電抗及時間常數;同時采用凍結磁導率法[8]計算出每一時刻空載突然三相短路時定、轉子繞組漏抗[9-12]及直、交軸電樞反應電抗,將得到的各繞組漏抗值分別取平均值并代入瞬態參數定義式中,得到兩種工況下的瞬態參數。將凍結磁導率法計算值與設計值分別代入短路電流解析式中,計算得到短路電流波形。與時步有限元仿真波形相比,凍結磁導率法計算波形與之更為接近。進而,將包絡線法和凍結磁導率法應用于實驗室6 極反裝水輪發電機,與實驗結果對比驗證了兩種計算方法的有效性和可行性。

2 場路耦合時步有限元模型

2.1 抽水蓄能電機基本參數

響水澗抽水蓄能電機額定數據見表1。對抽水蓄能電機數學模型做如下假設[13]:

(1)將電機分為直線和端部兩部分。直線部分用二維電磁場有限元計算,端部效應用集中參數的方法計入,體現在電路方程中。

(2)將問題作為二維恒定磁場來處理。

(3)渦流損耗及趨膚效應忽略不計,忽略鐵磁材料的磁滯效應,不計定轉子鐵心外緣散磁。

(4)在整個計算過程中,電機轉速不變。

(5)定子繞組為Y 聯結,無中線。

2.2 抽水蓄能電機磁場方程

抽水蓄能電機單元電機的二維電磁場求解區域如圖1 所示。根據上述假設條件,用矢量磁位A 表述的瞬態電磁場邊值問題為[13]

圖1 求解區域Fig.1 Solving region

式中 Js——電流密度的軸向分量;

μ(x,y)——有效磁導率;

σ——電導率。

2.3 抽水蓄能電機回路方程

將抽水蓄能電機二維電磁場方程(1)和繞組電路方程相耦合,即可求得其場路耦合時步有限元模型。發電機工況時定子繞組、勵磁繞組及阻尼繞組的等效電路[14]如圖2 所示。由基爾霍夫定律得定子繞組、勵磁繞組回路方程[13]為

表1 響水澗抽水蓄能電機額定數據Tab.1 Rated data of Xiang Shui Jian pumped storage motor

式中,U=(uA,uB,uC, uf)T;I=( iA,iB,iC, if)T;R=diag(rA,rB,rC,-rf)T;L=diag(lA,lB,lC,-lf)T;E=diag(eA,eB,eC,ef)T分別為定子及勵磁繞組端電壓矩陣、電流矩陣、電阻矩陣、電感矩陣及感應電動勢矩陣。

圖2 定、轉子繞組等效電路Fig.2 Equivalent circuits of stator and rotor windings

電樞繞組的感應電動勢可通過與繞組鉸鏈的磁通的變化得到,每相繞組的感應電動勢[13]為

式中 Ns——一相繞組串聯總匝數;

lef——電機軸向有效長度;

S——一相定子繞組電流流入端(或流出端)槽面積之和,S+為電流流出端,S-為電流的流入端。

對轉子阻尼回路的處理采用文獻[15]中的方法,對阻尼條和端環構成的回路單獨列寫方程表示。阻尼條中電流可表示為[14]

由圖2可得阻尼條支路電流與回路電壓方程[14]

綜合式(4)、式(5),經離散可得阻尼回路方程[14]為

式中,Ud=(ud1…udi…udk)T;Id=(id1…idi…idk)T,各系數矩陣具體元素見文獻[15]。

根據方程(1)~式(3)以及式(6),建立抽水蓄能電機單元電機的場路耦合模型如圖3 所示。

圖3 單元電機場路耦合模型Fig.3 Field-circuit coupled model of unit machine

圖3a 表示場的部分,圖3b 表示外電路部分。電動機工況時定子繞組電流方向與發電機工況時相反。對轉子運動的處理采用運動邊界法[16]。

3 瞬態參數計算

基于上述建立的抽水蓄能電機單元電機的場路耦合模型,分別采用包絡線法和凍結磁導率法對發電及電動工況的瞬態參數進行計算。

3.1 定、轉子繞組漏抗的計算

當抽水蓄能電機分別運行于發電及電動工況時,令電機運行于空載突然三相短路狀態,保存短路狀態下定轉子鐵心的磁導率不變,采用凍結磁導率法[8]計算定子繞組、勵磁繞組、阻尼繞組的漏感及直、交軸電樞反應電感,再用定義式計算電抗值。

3.1.1 電樞繞組漏抗及電樞反應電抗的計算

首先將電機運行狀態設置為空載突然三相短路狀態,當轉子旋轉到+d 軸與+A 軸重合時,保存此時定轉子鐵心的磁導率,令勵磁繞組和阻尼繞組開路,電樞繞組加直軸電樞電流,求得此時電樞繞組的電抗加上端部漏抗即為直軸同步電抗。再將轉子鐵心磁導率賦為空氣磁導率,定子鐵心磁導率仍為該短路狀態下保存的磁導率,令勵磁繞組和阻尼繞組開路,電樞繞組加直軸電樞電流,求解得到電樞繞組漏抗x1。求解完成后,將轉子鐵心磁導率賦為該短路時刻保存的磁導率,令勵磁繞組和阻尼繞組均閉合,電樞繞組所加電流賦為0,將電機狀態恢復到短路狀態。

轉子旋轉到+q 軸與+A 軸重合時,采用與計算直軸同步電抗相同的方法得到交軸同步電抗。

空載突然三相短路時,以計算A相飽和同步電抗為例,通過二維非線性磁場計算,可得繞組橫截面上各單元節點的矢量磁位A,設電機每極每相槽數為q,則A相一個相帶的磁鏈為[6]

式中,AUk、ALk分別為A相一個相帶中第k個線圈上層和下層載流導體中心處的矢量磁位。

計及定子端部漏抗時,A相的同步電抗標幺值為[6]

式中 β——一相串聯的相帶數;

a——并聯支路數;

IAm——A相電流幅值;

UΦN,IΦN——額定相電壓和相電流;

xσE——定子一相繞組的端部漏抗。

同理,計算B、C相繞組電抗標幺值,取三相平均值作為同步電抗或漏抗的飽和值。將同步電抗減去漏抗即得電樞反應電抗。轉子每旋轉一次即求一次同步電抗和電樞繞組漏抗,取平均值作為電樞繞組的電抗值。

3.1.2 勵磁繞組漏抗的計算

當電機運行在空載突然三相短路狀態時,在某一時刻保存短路狀態下定、轉子鐵心的磁導率,保持轉子鐵心磁導率不變,將定子鐵心磁導率賦為1,令電樞繞組和阻尼繞組開路,勵磁繞組加勵磁電流,得到勵磁繞組電流流入端的平均矢量磁位AR+和勵磁電流流出端的平均矢量磁位AR-,由式(9)計算鉸鏈于勵磁繞組的總磁鏈[10]為

式中 p——電機極對數;

Nf——每極下勵磁繞組匝數。

則勵磁繞組漏抗標幺值為

式中 if——所加的勵磁電流;

xl1——勵磁繞組端部漏抗;

Zf——勵磁繞組阻抗基值。

求解完成后,將定子鐵心磁導率賦為該短路狀態下保存的磁導率,令電樞繞組和阻尼繞組均閉合,將電機恢復到短路狀態。

轉子每旋轉一次即求解一次勵磁繞組漏抗值,取平均值作為勵磁繞組漏抗。

3.1.3 阻尼繞組漏抗計算

求解阻尼繞組漏抗時,將阻尼繞組看成同心式繞組[17]。以計算直軸阻尼繞組漏抗為例,當電機運行于空載突然三相短路狀態時,在某一時刻保存短路狀態下定轉子鐵心的磁導率,令電樞繞組和勵磁繞組均開路,給+d 軸兩側的阻尼條分別加大小相等方向相反的阻尼電流,求出每一個阻尼回路的磁鏈,求和,即為整個電機阻尼繞組的總磁鏈。磁鏈計算公式如下所示:

式中 ADi+,ADi-——一個阻尼回路電流流入端和電流流出端的平均矢量磁位;

n——阻尼回路數。

阻尼繞組漏抗標幺值為

式中 iD——阻尼條中的電流;

x1D——阻尼繞組端部漏抗;

Z1D——直軸阻尼繞組阻抗基值。

求解完成后,將定子鐵心磁導率賦為該短路狀態下保存的值,令電樞和勵磁繞組均閉合,將阻尼條的電流恢復為0,將電機恢復到短路狀態。

交軸阻尼繞組漏抗計算方法同直軸。

轉子每旋轉一次即求一次阻尼繞組漏抗值,而在短路時,阻尼繞組中的電流很快衰減到0,此時阻尼繞組便不再起作用,因此取電流衰減為0 前的漏抗的平均值作為阻尼繞組漏抗。

將定、轉子鐵心的磁阻率賦予磁化曲線直線部分的值,采用上述方法得到的即為各繞組不飽和漏抗及直、交軸電樞反應電抗。

3.2 瞬態參數計算

將上述用凍結磁導率法計算得到的各繞組漏抗及電樞反應電抗代入到瞬態參數定義式中,即可得到兩種工況下電機的瞬態參數。

用包絡線法計算瞬態參數采用二階模型,其具體計算方法見文獻[7]。

4 瞬態參數計算結果

4.1 凍結磁導率法計算結果

采用上述方法計算可得空載突然三相短路時每一時刻發電及電動工況的電樞繞組漏抗x1、直軸電樞反應電抗xad、交軸電樞反應電抗xaq、勵磁繞組漏抗xf1、直軸阻尼繞組漏抗xD1、交軸阻尼繞組漏抗xQ1的值(飽和值,標幺值)見表2。

表2 幾個典型時刻發電及電動工況各繞組漏抗(標幺值)Tab.2 Winding leakage reactance in power generation and electric working conditions of several typical moment

由表2 可知,隨著時間的增加,兩種工況下xad、x1、xf1、xD1、xaq、xQ1均呈增大趨勢。表明空載突然三相短路(非線性)時,電機主磁路及各繞組漏磁路的磁導率均逐漸增大,磁場的飽和程度逐漸減小,電樞反應磁場的去磁作用逐漸增大,穩態短路時電樞反應磁場的去磁作用最大、磁場的飽和程度最小。所以傳統理論認為短路時磁路不存在飽和是不完善的。

各繞組漏抗平均值見表3。

表3 發電及電動工況各繞組漏抗平均值(標幺值)Tab.3 The average winding leakage reactances in power generation and electric working conditions

由表3 可知,飽和對抽水蓄能機組運行于發電和電動工況的參數計算結果影響均較大;發電和電動工況的各參數對應值均比較接近。

根據瞬態參數的定義式[17]得到發電及電動工況時瞬態、超瞬態電抗及時間常數的計算值見表4。

表4 瞬態參數計算值與設計值對比(標幺值)Tab.4 Contrast between the calculated and the design values of transient parameters

4.2 包絡線法計算結果

抽水蓄能電機定子由空載變為突然三相短路時,短路初相角為42o,勵磁繞組電壓為空載勵磁電壓,通過對場路耦合模型中與電樞電壓源串聯的電阻Rs的控制,來模擬抽水蓄能電機在空載情況下發生突然三相短路。首先令Rs=1 000MΩ 來模擬抽水蓄能電機的空載穩態運行,然后令Rs=0Ω 模擬抽水蓄能電機的短路運行[13]。

6.5s 內非線性空載突然三相短路電流波形如圖4 所示。圖5為A相短路電流周期、非周期分量的幅值曲線的擬合結果。

圖4 發電機工況定、轉子短路電流Fig.4 Stator and rotor short-circuit currents in generator working conditions

圖5 A相電流擬合曲線Fig.5 Fitting curves of phase A current

由圖5 可見,發電機工況周期分量和非周期分量幅值曲線的擬合效果都非常好。

電動機工況時短路電流波形以及擬合曲線與發電機工況時類似,但數值略有差別。包絡線法計算結果見表4。

由表4 可知:①兩種方法的計算結果及設計值的偏差均在合理的范圍內;②不飽和以及飽和情況下電動機包絡線法及凍結磁導率法計算得到的瞬態參數均比發電機的大;③兩種工況下包絡線法計算得到的時間常數基本相等,凍結磁導率法計算得到的電動機的時間常數比發電機的偏大。

將凍結磁導率法計算值、設計值分別代入到短路電流解析式[18]中得到短路電流的包絡線與時步有限元仿真波形(非線性)的包絡線的對比如圖6 所示。

圖6 發電機工況A相短路電流包絡線的對比Fig.6 Contrast of phase A short circuit current envelope in generator working conditions

由圖6 可見,與設計值計算結果相比,凍結磁導率法計算值代入解析式中得到的短路電流波形與時步有限元仿真波形更接近。

5 模型機實驗對比

為驗證包絡線法及凍結磁導率法計算方法和結果的正確性,據上述假設條件建立了6 極反裝水輪發電機的場路耦合時步有限元模型,得到的兩種計算結果與實測波形包絡線法計算結果非常接近,表明所建模型正確且此兩種計算方法及結果合理。

5.1 反裝水輪發電機模型

反裝水輪發電機模型機基本數據見表5。反裝水輪發電機求解區域如圖7 所示。

表5 反裝水輪發電機基本數據Tab.5 The basic data of anti-mounted hydro-generator

圖7 反裝水輪發電機求解區域Fig.7 Solving region of anti-mounted hydro-generator

5.2 模型機瞬態特性仿真與實驗驗證

反裝水輪發電機在端電壓為220V 時發生突然三相短路,其A相短路電流波形及勵磁電流波形與實測波形對比如圖8 所示。B、C相仿真與實測短路電流波形的對比與A相類似。

由圖8 可見,模型機短路電流仿真波形與實測波形吻合較好,表明建模方法是正確的。

圖8 模型機短路電流仿真與實驗對比Fig.8 Contrast between simulation and experiment of the model machine short circuit currents

5.3 瞬態參數計算

5.3.1 凍結磁導率法計算結果

采用凍結磁導率法計算得到的反裝水輪發電機電樞繞組、勵磁繞組和阻尼繞組漏抗及直、交軸電樞反應電抗(飽和值)見表6。

表6 反裝水輪發電機各繞組漏抗(標幺值)Tab.6 Winding leakage reactance of hydro-generator

將表6 中的計算值代入各瞬態參數定義式中得到的瞬態參數及時間常數見表7。

表7 瞬態參數計算值與實驗值對比(標幺值)Tab.7 Contrast between the calculated and the experimental values of transient parameters

5.3.2 包絡線法計算結果

將A相短路電流仿真波形與實測波形的周期分量和非周期分量的幅值曲線分別進行雙指數和單指數擬合,得到的瞬態參數值(飽和值)見表7。時步有限元仿真波形的曲線擬合圖如圖9 所示。

圖9 反裝水輪發電機A相電流擬合曲線Fig.9 Phase A current fitting curves of anti-mounted hydro-generator

由圖9 可見,反裝水輪發電機短路電流的周期分量和非周期分量擬合效果都非常好。

由表7 知,反裝水輪發電機分別由凍結磁導率法、仿真包絡線法及實驗包絡線法計算得到的xd、x′d、x″d非常接近,T′d相差均小于0.124,T″d、Ta對應值相差均小于0.017。

將圖8 中短路電流波形與表7 中包絡線計算結果相比可知:

(1)圖8 中反裝水輪發電機在突然短路時超瞬態過程并不明顯,短路電流幅值與穩態時相差并不大,表明模型機的漏抗較大;且時步有限元仿真波形的幅值比實測波形的幅值略小,表明仿真模型的漏抗比實測波形的略大;而表7 包絡線的計算結果中,x′d、x″d較大且仿真包絡線計算值比實驗包絡線計算值偏大,表明模型機漏抗較大且仿真模型的漏抗比實測波形的大,兩者結果相符。

(2)圖8 中模型機實測波形在0.5s 時已達到穩定,而不計轉速變化時仿真波形在0.8s 時才穩定,表明仿真波形的T′d比實測波形的大;且勵磁電流仿真波形振蕩約6個周波后開始平緩下降,實測波形振蕩約8個周波后開始平緩下降,表明仿真波形的Ta比實測波形的小,與表7 辨識結果也相符。

將凍結磁導率法計算值及實驗包絡線法計算值分別代入短路電流解析式[18]中得到的電流波形與實測波形對比如圖10 所示。

圖10 反裝水輪發電機A相短路電流對比Fig.10 Contrast of phase A short circuit current of anti-mounted hydro-generator

由圖10 可見:

(1)在短路0.1s 內,與將實驗包絡線法計算值代入短路電流解析式所得波形相比,將凍結磁導率法計算值代入解析式所得電流波形與實測波形更接近,表明凍結磁導率法計算值與實際值更為接近。

(2)0.1s 以后,將實驗包絡線法計算值及凍結磁導率法計算值分別代入解析式所得電流波形幾乎完全重合,兩者的電流幅值與實測波形非常接近,但相位有一定偏差。由于空載突然三相短路電流的解析式是在疊加原理且假設瞬態過程中轉子轉速恒為同步轉速的條件下得到的,所以解析式計算得到的電流波形與實測波形有一定的相位差。

表7、圖8 及圖10 表明凍結磁導率法及包絡線法的計算方法及結果正確。

6 結論

本文基于單元電機建立了抽水蓄能電機的場路耦合時步有限元模型。分別采用包絡線法和凍結磁導率法計算了空載突然三相短路時發電及電動工況的瞬態參數。

(1)包絡線法計算值、凍結磁導率法計算值及設計值非常接近,表明瞬態參數計算方法有效。

(2)空載突然三相短路狀態下,電動機工況飽和及不飽和時電樞繞組、勵磁繞組、阻尼繞組漏抗及直、交軸電樞反應電抗與發電機工況的對應值均比較接近。

(3)與設計值相比,將凍結磁導率法計算值代入短路電流解析式中得到的波形與時步有限元仿真波形吻合更好,表明凍結磁導率法計算結果更加準確。

(4)將這兩種瞬態參數計算方法應用于實驗室6 極反裝水輪發電機,與實驗結果對比驗證了兩種計算方法的有效性和可行性。

[1]王學良,于繼來.分布式抽水蓄能系統的運營策略及其效益評估[J].電力系統保護與控制,2012,40(7):130-142.Wang Xueliang,Yu Jilai.The operation strategy and its benefit assessment of the distributed pumped storage system[J].Power System Protection and Control,2012,40(7):130-142.

[2]梁艷萍,湯蘊璆,周封.大型汽輪發電機瞬態參數的數值計算[J].電工技術學報,1999,14(2):5-10.Lang Yanping,Tang Yunqiu,Zhou Feng.Numerical calculation of transient parameters of large turbogenerator[J].Transactions of China Electrotechnical Society,1999,14(2):5-10.

[3]陳世坤.電機設計[M].2 版.北京:機械工業出版社,1990.

[4]梁艷萍,湯蘊璆.汽輪發電機的交軸瞬態和超瞬態電抗[J].電機與控制學報,1997,1(2):90-93.Lang Yanping,Tang Yunqiu.Q axis transient and subtransient reactances of turbogenerator[J].Electric Machines and Control,1997,1(2):90-93.

[5]梁艷萍,湯蘊璆.水輪發電機瞬態電抗的數值計算[J].電工技術學報,1996,11(6):22-26.Lang Yanping,Tang Yunqiu.Numerical calculation of transient reactance of waterwheel generator[J].Transactions of China Electrotechnical Society,1996,11(6):22-26.

[6]湯蘊璆,梁艷萍.電機電磁場的分析與計算[M].北京:機械工業出版社,2010.

[7]梁艷萍,周封.用時步有限元法計算汽輪發電機直軸瞬態參數[J].電機與控制學報,1998,2(2):69-74.Lang Yanping,Zhou Feng.Calculation of D axis transient parameters of turbogenerator by timestepping finite element method[J].Electric Machines and Control,1998,2(2):69-74.

[8]李和明,張健,羅應立,等.考慮交叉飽和影響的永磁同步電機穩態參數有限元分析[J].中國電機工程學報,2012,32(12):104-110.Li Heming,Zhang Jian,Luo Yingli,et al.Finite element analysis of PMSM steady state parameters considering cross-saturation effect[J].Proceedings of the CSEE,2012,32(12):104-110.

[9]Fuchs E F,Erdelyi E A.Determination of waterwheel alternator transient reactances from flux plots[J].IEEE Transactions on Power Apparatus and Systems,1972:1795-1802.

[10]Shima K,Ide K,Takahashi M.Calculation of leakage inductances of a salient-pole synchronous machine using finite elements[J].IEEE Transactions on Energy Conversion,1999,14(4):1156-1161.

[11]Kazuo Shima,Kazumasa Ide,Miyoshi Takahashi.Finite-element calculation of leakage inductances of a saturated salient-pole synchronous machine with damper circuits[J].IEEE Transactions on Energy Conversion,2002,17(4):463-470.

[12]Hiramatsu D,Koyanagi K,Hirayama K.Estimation of field mutual leakage reactance in synchronous machine by line-to-line sudden short circuit test[C].IEEE Transaction on Power Engineering Society General Meeting,2004:1359-1366.

[13]張燕燕.基于時步有限元的巨型水輪發電機動態特性及參數研究[D].北京:華北電力大學,2012.

[14]胡笳,羅應立,劉曉芳,等.汽輪發電機小擾動特性及靜態穩定性的時步有限元分析[J].中國電機工程學報,2010,30(6):76-82.Hu Jia,Luo Yingli,Liu Xiaofang,et al.Time-step finite element analysis for small disturbance characteristics and static stability of turbo generators[J].Proceedings of the CSEE,2010,30(6):76-82.

[15]胡敏強,黃學良.電機運行性能數值計算方法及其應用[M].南京:東南大學出版社,2003.

[16]孫玉田,楊明,李北芳.電機動態有限元法中的運動問題[J].大電機技術,1997,26(6):35-39.Sun Yutian,Yang Ming,Li Beifang.The moving problem in the dynamic FEM of electric machines[J].Large Electric Machine and Hydraulic Turbine,1997,26(6):35-39.

[17]馬志云.電機瞬態分析[M].北京:中國電力出版社,1998.

[18]李發海,朱東起.電機學[M].北京:科學出版社,2007.

主站蜘蛛池模板: 91青青草视频| 国产亚洲精品va在线| 久久国产黑丝袜视频| 91福利国产成人精品导航| 国产精品无码翘臀在线看纯欲| 国产精品人莉莉成在线播放| 女人天堂av免费| 九色在线观看视频| 无码精品国产dvd在线观看9久| 亚洲—日韩aV在线| 韩日无码在线不卡| 视频一区亚洲| 国产99视频精品免费视频7| 日韩黄色精品| 在线高清亚洲精品二区| 成人免费一级片| 欧美区一区| 国产粉嫩粉嫩的18在线播放91| 免费无码又爽又黄又刺激网站| 人人澡人人爽欧美一区| 欧美激情综合| 国产麻豆91网在线看| 国内精品一区二区在线观看| 伊人久久综在合线亚洲2019| 国产成人一区二区| 国产免费久久精品99re丫丫一| 她的性爱视频| 九九热精品视频在线| 22sihu国产精品视频影视资讯| 亚洲欧美一区二区三区蜜芽| 欧美国产日韩在线| 67194在线午夜亚洲 | 亚洲视频欧美不卡| 免费A级毛片无码免费视频| 国产精品第一区| 国产精品爽爽va在线无码观看 | 国产一区二区福利| 色网在线视频| 一级不卡毛片| 欧美日本中文| 亚洲啪啪网| 国产成人免费观看在线视频| 91久久偷偷做嫩草影院精品| 国产丝袜无码精品| 精品99在线观看| 欧美中文字幕一区二区三区| 亚洲码一区二区三区| 日韩少妇激情一区二区| 人妻中文字幕无码久久一区| 亚洲日本中文字幕天堂网| 91丝袜在线观看| 精品人妻一区二区三区蜜桃AⅤ | 国产又色又爽又黄| 亚洲日本www| 亚洲精品国产自在现线最新| 亚洲精品你懂的| 一级毛片免费不卡在线视频| 久久人妻xunleige无码| 色妺妺在线视频喷水| 一区二区三区国产精品视频| 日韩欧美网址| 国产h视频在线观看视频| 九色91在线视频| 成人无码一区二区三区视频在线观看 | 国产高清又黄又嫩的免费视频网站| 亚洲中文字幕在线观看| 国产人人射| 色综合日本| 中文字幕在线看| 伊人色在线视频| 国产二级毛片| 国产福利一区在线| 国产成年无码AⅤ片在线| 国产理论最新国产精品视频| 亚洲人成影视在线观看| 久久精品人妻中文视频| 亚洲浓毛av| 精品国产一二三区| 亚洲一欧洲中文字幕在线| 欧美日韩国产精品va| 2018日日摸夜夜添狠狠躁| 日本午夜精品一本在线观看|