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

基于定點諧波平衡法的鐵心磁滯與損耗特性分析

2014-11-25 09:34:04趙小軍趙志剛李慧奇
電工技術學報 2014年7期

趙小軍 崔 燦 李 琳 趙志剛 劉 剛 李慧奇

(1.華北電力大學電力工程學院 保定 071003 2.華北電力大學電氣與電子工程學院 北京 102206 3.河北工業大學電磁場與電器可靠性省部共建重點實驗室 天津 300130)

1 引言

不同勵磁條件下電工材料的磁性能并不相同。與交流勵磁相比,當電工鋼片或變壓器鐵心承受直流偏磁時,其磁化特性及損耗特性均會發生明顯的改變[1,2]。在實際工程中,電力變壓器經常承受直流偏磁或諧波形式的勵磁,這會給變壓器本身及電網的運行帶來一系列危害[3,4]。各種不同勵磁條件下電工材料磁特性的準確模擬對于電磁場理論研究和輸變電設備的安全運行具有重要的意義。目前有以下相關研究工作:一方面,基于磁化特性、損耗特性的測量和模擬對電工鋼片及鐵心的磁性能進行研究,如等效磁路長度的測定[5],鐵心中不同區域的有功及無功分離[6],各種磁滯模型的提出和改進[7-9],開發不同的測量設備及設計相關的測量方案等[10,11],這些工作取得的研究成果為面向工程應用的數值仿真奠定了良好的數據基礎;另一方面,時步有限元,諧波平衡有限元和時間周期有限元等方法[12-14]被用于非線性磁場的仿真和分析,不同的數值計算方法和多種仿真技術可以為材料性能模擬的深入研究提供準確性及有效性方面的支持。

為了準確計算變壓器鐵心中的非線性磁場,深入分析不同勵磁條件下鐵心的磁化特性和損耗特性,需要在場路耦合計算中考慮鐵心的磁滯效應[15]。解決該問題的關鍵在于磁滯模型的選擇和如何對磁場中的非線性本構關系進行處理。選擇磁滯模型時,應該遵循簡單有效、數值計算時易實現的原則。對于含磁滯效應的非線性磁場本構關系的處理,目前主要有兩種方法:①引入基于空氣磁導率的μ0-B-H-M 關系[16],②是引入基于定點磁阻率的νFP-B-H-M[17]關系。二者都能解決基于B-H 磁場本構關系中磁阻率ν 的不連續性問題。

本文利用疊片鐵心模型分別進行了正弦勵磁和直流偏磁勵磁下的實驗,得到相應的無偏磁磁滯回線和直流偏磁磁滯回線。利用基于損耗函數的磁滯模型對磁滯回線進行擬合和預測,并將該磁滯模型與定點諧波平衡有限元算法相結合,計算無偏磁和有偏磁條件下繞組勵磁電流及鐵心內的非線性磁場。比較了不同區域內的鐵心磁滯特性,分析了不同勵磁條件下鐵心損耗的分布特征。

2 鐵心磁特性實驗

圖1 疊片鐵心模型Fig.1 Laminated core model

圖1 所示為疊片鐵心模型,鐵心上繞有匝數相同的勵磁線圈和測量線圈?;诏B片鐵心進行空載實驗,利用功率分析儀分別測量鐵心損耗,勵磁線圈的勵磁電流i 和測量線圈中的感應電壓u。疊片鐵心所用硅鋼片為30Q140 型取向硅鋼片。

在無偏磁條件下,在勵磁線圈端口施加不同的交流電壓勵磁,由下式可以得到相應的無偏磁磁滯回線。

式中,S 和L 分別為鐵心截面積和等效磁路長度;φ為鐵心磁通;Ncoil為線圈匝數。

本文中等效磁路長度L 取疊片鐵心的幾何平均磁路長度。

由無偏磁條件下的測量結果,得到不同勵磁條件下的磁滯回線,如圖2 所示。

圖2 無偏磁磁滯回線Fig.2 Hysteresis loops under sinusoidal excitation

將一直流電流源與勵磁線圈所在回路相串聯,模擬變壓器鐵心的直流偏磁情況。此時鐵心磁通φ包含兩部分,即直流磁通φdc和交流磁通φac。φac可按照式(1)進行計算,φdc則可以通過迭代法進行計算[18,19]。

由直流偏磁條件下的測量結果,得到不同勵磁條件下的磁滯回線,如圖3 所示。

圖3 直流偏磁磁滯回線(Idc=0.426A)Fig.3 Hysteresis loops under DC-biased excitation

3 磁滯模型與磁滯回線的擬合

3.1 基于損耗函數的磁滯模型

在正弦勵磁條件下,如圖4 所示,可以將勵磁電流i 分為兩部分,一部分im與磁通φ 同相位,一部分ih與感應電動勢e 同相位[20]。于是,在i-φ 關系中可以引入以下函數

式中,i1對應于圖4 中im;i2對應于圖4 中ih。

圖4 無偏磁條件下的勵磁電流Fig.4 Magnetizing current under sinusoidal excitation

基于以上分析,對于如圖5 所示的i-φ 磁滯回線,曲線“aOc”為中間磁化曲線(即i-φ 磁滯回線的中點軌跡),其所代表的勵磁特性與式(4)中的i1相對應,i-φ 磁滯回線外圍上任意一點與中間磁化曲線“aOc”的間距“ef”反映了鐵心的磁滯效應且與損耗相關,與式(4)中的i2相對應,稱之為損耗函數[21,22]。損耗函數的變化趨勢與中間磁化曲線相反,當磁化曲線上升至最高點時,損耗函數值為零;當磁化曲線下降至零點時,損耗函數值則為最大,與“Ob”相對應。根據損耗函數的特點,可作如下定義

圖5 基于損耗函數的磁滯模型Fig.5 Hysteresis model based on consuming function

式中,φm為磁通的幅值;D 是與損耗函數相關的系數。

損耗函數i2可進一步寫為以下形式

式中,Iob為損耗系數,與磁通幅值φm及頻率f 均相關,可以通過實驗中的測量結果計算得到。

由式(2)、式(3)可知,對于疊片鐵心,可以提出基于損耗函數的B-H 磁滯模型為

式中,Bm為無偏磁條件下交流磁通密度的幅值;H為鐵心中總的磁場強度。

同理,在直流偏磁條件下,基于損耗函數的B-H磁滯模型如下

式中,Bdc為直流磁通密度,Bdc=φdc/S。

由于各中間磁化曲線均通過磁滯回線的頂點,其軌跡為基本磁化曲線,在有限元計算中可以分別選擇無偏磁和直流偏磁下的基本磁化曲線作為式(8)中的H1(B) 和H1d(B)[19]。

3.2 磁滯模型的驗證

通過實驗,可以對測量結果進行處理,進而得到不同勵磁條件下磁滯回線所對應的損耗系數Hob,結果見表1~表3。

表1 無偏磁條件下的損耗系數(Idc=0A)Tab.1 Consuming coefficients under sinusoidal excitation(Idc=0A)

表2 直流偏磁條件下的損耗系數(Idc=0.426A)Tab.2 Consuming coefficients under DC-biased excitation(Idc=0.426A)

表3 直流偏磁條件下的損耗系數(Idc=0.847A)Tab.3 Consuming coefficients under DC-biased excitation(Idc=0.847A)

由此可利用基于損耗函數的磁滯模型對測量得到磁滯回線進行仿真,無偏磁磁滯回線的結果如圖6 所示,直流偏磁磁滯回線的結果如圖7 所示。通過比較可以看出,基于損耗函數的磁滯模型能夠較好的模擬無偏磁和直流偏磁條件下疊片鐵心的磁滯效應。

圖6 無偏磁磁滯回線的測量與仿真結果Fig.6 Simulated and measured hysteresis loop under sinusoidal excitation

圖7 直流偏磁磁滯回線的測量與仿真結果Fig.7 Simulated and measured hysteresis loop under DC-biased excitation

4 計算結果及其分析

4.1 諧波平衡法

在無偏磁條件下,變壓器端口承受穩態電壓勵磁,勵磁電流中只含有奇次諧波,因此勵磁電流i可以表達成如下形式

式中,I 為勵磁電流密度i 的諧波矢量表達式。

同理可知在直流偏磁條件下,i 的諧波矢量表達式為

在穩態勵磁下,電磁場中的各物理量均具有周期性,因此各變量均可表達為式(10)、式(11)所示的形式。

4.2 基于B-H-M 的定點諧波平衡方程

基于巴拿赫不動點定理,當不考慮各向異性時,可在二維非線性磁場中引入如下關系[17]

式中,νFP為定點磁阻率,諧波矢量Bx、By、Hx、Hy、Mx、My的表達式同式(10)、式(11)。

由式(12)可以看出,定點磁阻率νFP將磁場強度H 分為了線性和非線性兩部分,線性部分與定點磁阻率相關,非線性部分與類磁化強度矢量M相關。

基于定點磁阻率的二維非線性磁場方程為

利用伽遼金法可以得到關于式(13)的定點諧波平衡有限元方程

式中,Q 為系數矩陣;Uk為線圈k 的端口電壓諧波矢量;Zk為阻抗矩陣;Ik為流過線圈k 的勵磁電流密度的諧波矢量;Ck為場路耦合矩陣;P 是與類磁化強度矢量M 相關的諧波矢量[23]。

4.3 νFP的取值與非線性場的求解

依據式(14)可對勵磁電流和磁矢量位的諧波矢量同時進行求解,在求解過程中,諧波解的收斂性取決于定點磁阻率νFP的取值。在非線性靜態場的計算中,定點磁阻率可按照下式確定[24]

式中,νdmax和νdmin分別為磁化曲線H=F(B)上微分磁阻率的最大值和最小值。

在時域有限元計算中,Dlala 等也給出了時步迭代中定點磁阻率的確定方案[25]。

在定點諧波平衡有限元計算中,初始迭代過程中采用較大的定點磁阻率(ν1=3ν0~ν0)可以保證諧波解的穩定收斂,之后采用較小的定點磁阻率(ν2=ν0/10~ν0/40)則能夠實現諧波解的快速收斂。

圖8 所示為考慮磁滯效應時的計算流程圖,其中H=F(B) 與不同勵磁條件下的磁滯回線相對應。忽略疊片鐵心的各向異性,Bx、By、Hx、Hy滿足以下關系

圖8 計算流程Fig.8 Flow chart of the computational process

通過式(16),可以實現迭代計算中由B 求解H的過程。

5 計算結果及分析

5.1 勵磁電流

對無偏磁和有偏磁條件下的勵磁電流和磁場進行計算,勵磁電流的計算結果如圖9~圖11 所示。從比較結果可以看出,計算結果與測量結果吻合較好。無偏磁條件下,由于磁滯效應的影響,勵磁電流的前半周與后半周沿橫軸(時間)不對稱,但是沿縱軸(電流)呈反對稱的特點。在直流偏磁條件下,由于勵磁電流中的各次諧波分量迅速增大,勵磁電流波形呈尖頂波形狀,此時磁滯效應對于電流波形的影響不再明顯,且電流波形也不再具有對稱性。由圖9 可以看出,在交流勵磁較小,磁滯效應的影響相對較大時,勵磁電流的計算結果與測量結果仍然存在一定的誤差。這是因為在無偏磁條件下,交流勵磁較小時疊片鐵心的磁滯回線呈橢圓形,而本文提出的基于損耗函數的磁滯模型不能準確模擬橢圓形磁滯回線,因此數值計算結果的誤差相對較大。要得到更加準確的結果,需要對該磁滯模型作進一步的改進和完善。

圖9 無偏磁條件下的勵磁電流Fig.9 Magnetizing current under sinusoidal excitation

圖10 直流偏磁條件下的勵磁電流(Idc=0.426A)Fig.10 Magnetizing current under DC-biased excitation(Idc=0.426A)

圖11 直流偏磁條件下的勵磁電流(Idc=0.847A)Fig.11 Magnetizing current under DC-biased excitation(Idc=0.847A)

5.2 鐵心內不同區域磁滯特性

由計算得到的磁場,可以得到鐵心內不同區域的磁滯特性。如圖12 所示,為用于計算的疊片鐵心幾何模型,其中陰影部分代表線圈所在區域。根據鐵心結構及磁場在鐵心中的分布特點,可以將鐵心分為兩個區域,柱-軛區和接縫區。在兩個區域中分別選擇A 點和B 點,基于計算結果可以得到各點的磁場及對應的磁滯特性。

圖12 疊片鐵心的幾何計算模型Fig.12 Geometric model of the laminated core for computation

由圖13~圖18 可以看出,在無偏磁條件下,疊片鐵心內各點的磁滯回線是對稱的,而直流偏磁條件下,疊片鐵心內各點的磁滯回線呈現出明顯不對稱的特征。這是導致無偏磁和有偏磁條件下疊片鐵心勵磁特性和損耗特性存在差異的主要原因。由圖14 和圖15 的比較還可以看出,A 點處不同方向的磁滯特性存在差異,主要原因是在疊片鐵心接縫區內的同一點處,磁通密度在不同方向上的分量并不相等。圖17 和圖18 的比較結果也反映了直流偏磁條件下接縫區內同一點處不同方向上磁滯特性的差異。

圖13 無偏磁條件下的B 點y 方向磁滯回線Fig.13 Hysteresis loop on point B under sinusoidal excitation

圖14 無偏磁條件下的A 點y 方向磁滯回線(Hy-By)Fig.14 Hysteresis loop on point A under sinusoidal excitation(Hy-By)

圖15 無偏磁條件下的A 點x 方向磁滯回線(Hx-Bx)Fig.15 Hysteresis loop on point A under sinusoidal excitation(Hx-Bx)

圖16 直流偏磁條件下的B 點y 方向磁滯回線(Idc=0.426A)Fig.16 Hysteresis loop on point B under DC-biased excitation(Idc=0.426A)

圖17 直流偏磁條件下的A 點x 方向磁滯回線(Hx-Bx,Idc=0.426A)Fig.17 Hysteresis loop on point A under DC-biased excitation(Hx-Bx,Idc=0.426A)

圖18 直流偏磁條件下的A 點y 方向磁滯回線(Hy-By,Idc=0.426A)Fig.18 Hysteresis loop on point A under DC-biased excitation(Hy-By,Idc=0.426A)

5.3 鐵心的損耗特性

由計算得到每個單元的B 和H,可以按照下式計算鐵心的損耗[26,27]

式中,Pi為第i 個單元中的比損耗;ρ 為鐵心疊片的密度;T 為周期;Ns和ds分別為鐵心中的疊片數量和單片厚度;Ne為單元總數;Si為單元面積;Pt為鐵心損耗。

表4 中給出了無偏磁條件下鐵心損耗的計算結果與測量結果,其中Pc為鐵心比損耗的計算值,Pm為鐵心比損耗的測量值,Bacm為鐵心中的交流磁通密度的幅值,可由式(2)計算得到。由比較可以看出,計算結果與測量結果相吻合。

表4 無偏磁條件下的鐵心損耗(Idc=0)Tab.4 The iron loss under sinusoidal excitation

由圖12 中A、B 兩點的磁場計算結果可以看出,鐵心中的磁場分布并不均勻,尤其是接縫區和柱-軛區的磁場分布明顯不同,由此可知,鐵心損耗在不同區域的分布也將不同。分別計算鐵心接縫區和柱-軛區的比損耗,觀察不同勵磁條件和勵磁形式對二者的影響,結果見表5 和表6,其中Pj為接縫區比損耗,PL為柱-軛區比損耗。

表5 無偏磁下鐵心不同區域比損耗(Idc=0)Tab.5 Specific loss in different areas in the laminated core under sinusoidal excitation

表6 偏磁下鐵心不同區域比損耗(Idc=0.426A)Tab.6 Specific loss in different areas in the laminated core under dc-biased excitation

由計算結果可以看出,當鐵心中交流磁通密度相同時,直流偏磁條件下的鐵心損耗大于無偏磁條件下鐵心損耗。同時,鐵心內的損耗并不均勻,主要體現在接縫區與柱-軛區的差異。在相同勵磁條件下,柱-軛區比損耗大于接縫區的比損耗。

6 結論

(1)基于疊片鐵心的實驗和計算表明,在無偏磁條件下和直流偏磁條件下,變壓器鐵心的磁滯特性和損耗特性明顯不同,磁滯特性表現為偏磁后的磁滯回線的不規則與不對稱性,損耗特性表現為偏磁后鐵心損耗的增大。

(2)基于損耗函數的磁滯模型簡單有效,能夠模擬鐵心在無偏磁和有偏磁條件下的磁滯效應。將其與定點技術和諧波平衡法相結合時,在數值計算中易實現,適用于頻域有限元計算。

(3)在無偏磁和有偏磁條件下,鐵心接縫區與柱-軛區磁場分布并不均勻,不同位置處的磁滯特性不同,并由此導致鐵心損耗分布的不均勻性。在鐵心接縫區內,同一位置處不同方向上的磁滯特性存在差異。

[1]趙志剛,劉福貴,張俊杰,等.直流偏磁條件下變壓器勵磁電流的實驗與分析[J].電工技術學報,2010,25(4):71-76.Zhao Zhigang,Liu Fugui,Zhang Junjie,et al.Measurement and analysis of magnetizing current in DC-biased transformers[J].Transactions of China Electrotechnical Society,2010,25(4):71-76.

[2]郭滿生,梅桂華,劉東升,等.直流偏磁條件下電力變壓器鐵心B-H 曲線及非對稱勵磁電流[J].電工技術學報,2009,24(5):46-52.Guo Mansheng,Mei Guihua,Liu Dongsheng,et al.B-H curve based on core and asymmetric magnetizing current in DC-Biased transformers[J].Transactions of China Electrotechnical Society,2009,24(5):46-52.

[3]段旭.直流偏磁下變壓器振動機理與振動信號分析研究[D].重慶:重慶大學,2013.

[4]李泓志.直流接地極接地性能及入地電流對變壓器影響的研究[D].北京:華北電力大學,2010.

[5]Cheng Z,Takahashi N,Forghani B,et al.Modeling of magnetic properties of GO electrical steel based on Epstein combination and loss data weighted processing[J].IEEE Transactions on Magnetics,2014,50(1):6300209.

[6]杜永,程志光,顏威利,等.電力變壓器全斜接縫疊片鐵心工作條件下的磁性能模擬[J].電工技術學報,2010,25(3):14-19.Du Yong,Cheng Zhiguang,Yan Weili,et al.Working magnetic property modelling of the power transformer laminated core with mitred joints[J].Transactions of China Electrotechnical Society,2010,25(3):14-19.

[7]曹林,何金良,張波,等.直流偏磁狀態下電力變壓器動態磁滯損耗模型及驗證[J].中國電機工程學報,2008,28(24):141-146.Cao Lin,He Jinliang,Zhang Bo,et al.Dynamic hysteresis loss model of power transformer under DC current biasing and its verification[J].Proceedings of the CSEE,2008,28(24):141-146.

[8]張艷麗,劉洋,謝德馨,等.耦合改進型矢量磁滯模型的變壓器磁場分析及實驗研究[J].中國電機工程學報,2010,30(21):109-113.Zhang Yanli,Liu Yang,Xie Dexin,et al.Finite element analysis of magnetic field in transformer core coupled with improved vector hysteresis model and its experimental verification[J].Proceedings of the CSEE,2010,30(21):109-113.

[9]劉洋,張艷麗,謝德馨,等.考慮硅鋼片矢量磁特性的復數 E&S 模型[J].中國電機工程學報,2012,32(3):144-149.Liu Yang,Zhang Yanli,Xie Dexin,et al.A complex E&S model considering 2D vector magnetic properties of silicon steel sheet[J].Proceedings of the CSEE,2012,32(3):144-149.

[10]Miyagi D,Yoshida T,Nakano M,et al.Development of measuring equipment of dc-biased magnetic properties using open-type single-sheet tester[J].IEEE Transactions on Magnetics,2006,42(10):2846-2848.

[11]Ishikawa S,Yanase S,Okazaki Y.AC magnetic properties of electrical steel sheet under twodimensional dc-biased magnetization[J].IEEE Transactions on Magnetics,2012,48(4):1413-1416.

[12]Yao Y,Koh C,Ni G.3-D nonlinear transient eddy current calculation of online power transformer under DC bias[J].IEEE Transactions on Magnetics,2005,41(5):1840-1843.

[13]Zhao X,Li L,Lu J,et al.Characteristics analysis of the square laminated core under dc-biased magnetization by the fixed-point harmonic-balanced FEM[J].IEEE Transactions on Magnetics,2012,48(2):747-750.

[14]Takahashi Y,Iwashita T,Nakashima H,et al.Parallel time-periodic finite-element method for steady-state analysis of rotating machines[J].IEEE Transactions on Magnetics,2012,48(2):1019-1022.

[15]劉福貴,楊慶新,顏威利,等.一種考慮磁滯特性的磁場數值計算方法[J].華北電力大學學報,2005,32(6):59-61.Liu Fugui,Yang Qingxin,Yan Weili,et al.A numerical method of magnetic field computation considering magnetic hysteresis property[J].Journal of North China Electric Power University,2005,32(6):59-61.

[16]Matsuo T,Shimasaki M.Time-periodic finite-element method for hysteretic eddy-current analysis[J].IEEE Transactions on Magnetics,2012,38(2):549-552.

[17]Mathekga M,Mcmahon R.Application of the fixed point method for solution in time stepping finite element analysis using the inverse vector Jiles-Atherton model[J].IEEE Transactions on Magnetics,2011,47(10):3048-3051.

[18]王祥珩,徐伯雄.變壓器的偏磁問題[J].變壓器,1992,29(8):11-14.Wang Xiangheng,Xu Boxiong.The problems on biased transformer[J].Transformer,1992,29(8):11-14.

[19]趙小軍,李琳,程志光,等.基于直流偏磁實驗的疊片鐵心磁化特性分析[J].電工技術學報,2011,26(1):7-13.Zhao Xiaojun,Li Lin,Cheng Zhiguang,et al.Analysis of magnetizing characteristic of laminated core based on the dc-biasing experiment[J].Transactions of China Electrotechnical Society,2011,26(1):7-13.

[20]Kulkarni S V,Khaparde S A.Transformer engineering:design and practice[M].New York:Marcel Dekker,2004.

[21]Lin C E,Wei J B,Huang C L,et al.A new method for representation of hysteresis loops[J].IEEE Transactions on Power Delivery,1989,4(1):413-420.

[22]Faiz J,Sharifian M.Hysteresis loop modeling techniques and hysteresis loss estimation of soft magnetic materials[J].COMPEL,2001,20(4):981-1001.

[23]趙小軍,李琳,程志光,等.定點諧波平衡有限元法與疊片鐵心直流偏磁磁化特性研究[J].中國電機工程學報,2011,31(9):126-132.Zhao Xiaojun,Li Lin,Cheng Zhiguang,et al.Fixed-point harmonic-balanced finite element methodand dc-biasing magnetizing characteristics of laminated core[J].Proceedings of the CSEE,2011,31(9):126-132.

[24]Hantila F,Preda G,Vasiliu M.Polarization method for static field[J].IEEE Transactions on Magnetics,2000,36(4):672-675.

[25]Dlala E,Arkkio A.Analysis of the convergence of the fixed-point method for solving time-stepping nonlinear field problems[J].IEEE Transactions on Magnetics,2008,44(4):473-478.

[26]張艷麗,李玉梅,劉洋,等.考慮不同磁特性模型的感應電機鐵心損耗分析[J].中國電機工程學報,2013,33(27):120-127.Zhang Yanli,Li Yumei,Liu Yang,et al.Analysis of core losses in induction motors considering different magetic property models[J].Proceedings of the CSEE,2013,33(27):120-127.

[27]Enokizono M,Shimoji H,Horibe T.Loss evaluation of induction motor by using magnetic hysteresis E&S2model[J].IEEE Transactions on Magnetics,2002,38(5):2379-2382.

主站蜘蛛池模板: 国产福利免费在线观看| 99热这里只有精品国产99| 久久精品人人做人人综合试看| 全部毛片免费看| 欧美色视频网站| 午夜国产理论| 国产夜色视频| 亚洲欧美不卡| 色偷偷综合网| 精品一区二区三区水蜜桃| 久久这里只有精品8| 伊人成人在线视频| 亚洲男人天堂2020| 国产精品久久久久久久久| 日韩福利在线观看| 日本91在线| 国产簧片免费在线播放| 美女无遮挡被啪啪到高潮免费| 亚洲av无码人妻| 亚洲色偷偷偷鲁综合| 日本在线亚洲| 综合色亚洲| 亚洲精品视频免费观看| 色婷婷狠狠干| 亚洲日韩第九十九页| 国产精品亚洲va在线观看| 亚洲精品第五页| 国产亚洲高清视频| 国产精品亚洲天堂| 国产精品久线在线观看| 欧美国产视频| 国内黄色精品| 好吊妞欧美视频免费| 欧美在线观看不卡| 久久黄色毛片| 激情综合网址| 日韩欧美高清视频| 久久亚洲天堂| 漂亮人妻被中出中文字幕久久| 精品国产Ⅴ无码大片在线观看81| 成人在线观看不卡| 国产打屁股免费区网站| 福利国产微拍广场一区视频在线| 国产另类视频| 在线人成精品免费视频| 中文字幕2区| 国产成人免费观看在线视频| 内射人妻无码色AV天堂| 中文字幕在线视频免费| 欧美三级视频网站| 91免费观看视频| 亚洲第一页在线观看| 97se亚洲综合不卡 | 国产亚洲欧美另类一区二区| 第一页亚洲| 国产va在线观看免费| 亚洲黄色激情网站| 亚洲成A人V欧美综合| 亚洲伊人天堂| 欧美成在线视频| 亚洲无码精彩视频在线观看| 亚洲美女久久| 日韩高清在线观看不卡一区二区| 国产乱子伦视频在线播放 | 亚洲IV视频免费在线光看| 亚洲大尺码专区影院| AV在线麻免费观看网站| 国产成人精品高清不卡在线| 一级毛片高清| 国产亚洲第一页| 国产成人夜色91| 欧美一区日韩一区中文字幕页| 亚洲无码高清视频在线观看| 亚洲精品午夜天堂网页| 在线国产毛片| 久久综合色天堂av| 国产91九色在线播放| 亚洲最新在线| 欧美日韩在线国产| 青青青亚洲精品国产| 国产91色| 亚洲人成人无码www|