• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    一種考慮動態(tài)磁滯效應(yīng)的高效穩(wěn)定時域有限元計(jì)算方法

    2023-11-11 06:12:42井立兵
    電工技術(shù)學(xué)報 2023年21期
    關(guān)鍵詞:磁滯回線磁阻磁通

    魏 鵬 陳 龍 賁 彤 井立兵 張 獻(xiàn)

    (1.三峽大學(xué)電氣與新能源學(xué)院 宜昌 443002 2.湖北省微電網(wǎng)工程技術(shù)研究中心(三峽大學(xué)) 宜昌 443002 3.省部共建電工裝備可靠性與智能化國家重點(diǎn)實(shí)驗(yàn)室(河北工業(yè)大學(xué)) 天津 300130)

    0 引言

    時域有限元方法廣泛應(yīng)用于電工裝備的運(yùn)行狀態(tài)預(yù)測與損耗評估之中[1],準(zhǔn)確預(yù)測電機(jī)、變壓器等電工裝備的時域暫態(tài)損耗分布特性不僅對降低設(shè)備能耗、改善局部過熱具有重要意義,同時考慮損耗對暫態(tài)磁場求解的反饋?zhàn)饔茫瑢M(jìn)一步分析電工裝備的瞬時功率平衡、暫態(tài)勵磁特征等電磁特性至關(guān)重要[2-4]。

    目前,對于電工裝備鐵心損耗的有限元計(jì)算方法可大致分為兩類:一是基于場計(jì)算結(jié)果的后處理方法,該類方法利用Steinmetz 損耗模型或Bertotti統(tǒng)計(jì)學(xué)損耗分離理論對僅單值非線性磁場的計(jì)算結(jié)果進(jìn)行后處理,得到每一個單元的損耗,最終進(jìn)行積分得到電工裝備整體的損耗[5-7]。這一類方法主要用于對電機(jī)或變壓器的損耗進(jìn)行頻域分析,并通過傅里葉變換來分析含有諧波激勵時的鐵心損耗。Ansys 公司D.Lin 等基于等效橢圓環(huán)思想,將上述方法進(jìn)行了時域瞬時功率損耗計(jì)算擴(kuò)展[8]。該方法的優(yōu)勢在于在損耗計(jì)算過程中無需迭代,在磁路簡單的結(jié)構(gòu)中具有較好的精度,但由于沒有考慮鐵心磁滯、渦流損耗、剩余損耗等效動態(tài)磁場對鐵心中磁場分布的反饋?zhàn)饔茫陔姍C(jī)等復(fù)雜磁路電工裝備的損耗計(jì)算中將帶來較大誤差。另一類鐵心損耗計(jì)算方法則是在時域有限元分析中,考慮動態(tài)磁滯效應(yīng)對每一瞬時的功率損耗進(jìn)行精確求解[9-10]。這一類方法可最大限度地模擬鐵心的真實(shí)磁化過程,在處理含有磁通密度波形畸變工況下的損耗計(jì)算問題時具有較高精度。同時,在考慮鐵心動態(tài)磁滯效應(yīng)后,鐵心中動態(tài)磁滯回線波形的精確模擬可在時域中對勵磁電流的波形進(jìn)行實(shí)時反饋,進(jìn)而可對電機(jī)或變壓器中暫態(tài)運(yùn)行特性進(jìn)行有效預(yù)測。因此,電工裝備損耗的精確預(yù)測與運(yùn)行狀態(tài)的有效評估都離不開在時域有限元分析時對動態(tài)磁滯特性的影響的考慮。

    然而,目前在有限元分析中常用的磁滯模型,如J-A 模型[11-13]、Preisach 模型[14-17]等均是與磁化速率無關(guān)的磁滯模型,并不能直接考慮電機(jī)或變壓器疊片宏觀渦流損耗與剩余損耗等動態(tài)因素對鐵心磁場分布的影響。為此,需要進(jìn)一步將二者的等效動態(tài)磁場考慮到有限元的計(jì)算方程之中。傳統(tǒng)方法是將動態(tài)效果等效為電導(dǎo)率乘以磁矢位A對時間的一階導(dǎo)數(shù),但該種等效僅適合均勻?qū)嵭膮^(qū)域渦流路徑的求解,并不適合疊片鐵心結(jié)構(gòu)[18-19]。對于鐵心疊片結(jié)構(gòu),進(jìn)行實(shí)際三維建模將會帶來無法承受的計(jì)算成本。為此,E.Dlala 等提出了一種耦合一維疊片有限元方法的二維時域有限元計(jì)算方法[20-21]。在疊片維度,忽略邊緣效應(yīng),認(rèn)為磁通密度僅是疊片厚度方向上的函數(shù),進(jìn)而僅需要求解一維有限元方程,即可得到疊片的渦流損耗等效磁場。然而,該方法實(shí)質(zhì)上是在二維有限元計(jì)算中嵌套了一維有限元程序,仍具有較高的求解難度。為了對上述求解方法進(jìn)行簡化,目前常采用Bertotti 模型的變形形式,即在耦合靜態(tài)磁滯模型的過程中,通過耦合等效微分磁阻率來考慮渦流損耗磁場與剩余損耗磁場的效果,該方法避免了一維有限元的求解,在一定程度上大大簡化了計(jì)算,但在每次迭代過程中仍需更新磁阻率和剛度矩陣,導(dǎo)致計(jì)算效率低下[22-23]。因此,在進(jìn)一步求解含有等效磁阻率的方程時,需要進(jìn)一步耦合非線性磁場計(jì)算中的固定點(diǎn)迭代技術(shù)。在僅考慮靜態(tài)磁滯模型的有限元計(jì)算中,固定點(diǎn)迭代方法已經(jīng)取得了較好的收斂效果與計(jì)算速度[24-26],并已在目前的商業(yè)軟件中取得了較好的計(jì)算效果。然而,在引入渦流損耗磁場與剩余損耗磁場后,收斂將難以得到保證,特別是在計(jì)算至磁通密度零點(diǎn)附近時間步時算法極易發(fā)散,分析其原因主要是渦流損耗、剩余損耗等效動態(tài)磁場的介入改變了固定點(diǎn)法中磁通密度與磁場強(qiáng)度的收斂過程,使得傳統(tǒng)微分磁阻率的定義無法準(zhǔn)確描述收斂路徑的斜率,計(jì)算得到微分磁阻率過低導(dǎo)致考慮動態(tài)磁滯效應(yīng)的有限元方程求解不收斂。

    為解決在考慮動態(tài)磁滯特性的時域有限元計(jì)算難以收斂的問題,本文提出一種基于等效磁阻率的固定點(diǎn)迭代求解策略,進(jìn)而提出一種耦合動態(tài)磁滯特性的時步有限元分析方法。首先,為了獲得電工鋼片的動態(tài)、靜態(tài)磁滯特性,基于愛潑斯坦方圈法搭建了一維磁特性測量系統(tǒng);其次,考慮到逆Preisach模型可以考慮磁化歷史具有較高的計(jì)算精度,同時在耦合有限元計(jì)算時易于數(shù)值實(shí)現(xiàn)的特點(diǎn),本文構(gòu)建了基于參數(shù)化解析一階回轉(zhuǎn)曲線的逆 Preisach模型,并進(jìn)行了參數(shù)辨識;再次,分析了渦流損耗和剩余損耗對算法收斂特性的影響,提出基于等效磁阻率的固定點(diǎn)迭代策略,并在有限元迭代過程中加入松弛因子保證穩(wěn)定性;最后,以愛潑斯坦方圈為例,進(jìn)行了考慮動態(tài)磁滯效應(yīng)的時域有限元分析,求解了瞬時功率損耗特性,并進(jìn)行了實(shí)驗(yàn)對比驗(yàn)證。結(jié)果驗(yàn)證了本文所提方法的準(zhǔn)確性、高效性與穩(wěn)定性。

    1 準(zhǔn)靜態(tài)與動態(tài)磁滯特性測量

    為建立磁滯模型及驗(yàn)證有限元計(jì)算準(zhǔn)確性,基于 IEC 60404—2 標(biāo)準(zhǔn)建立了一維磁性能測量平臺,對牌號為B30P105 的取向硅鋼片進(jìn)行準(zhǔn)靜態(tài)與動態(tài)磁滯特性測量。電工鋼片一維磁特性測量平臺如圖1 所示,該平臺包括寬頻功率放大器、電壓前置放大器、愛潑斯坦方圈、多功能數(shù)據(jù)采集卡和 LabVIEW 控制系統(tǒng)。本實(shí)驗(yàn)采用的硅鋼片規(guī)格及測試條件見表1。分別設(shè)定正弦交流電壓激勵的頻率為5 Hz 與50 Hz,測量得到準(zhǔn)靜態(tài)磁滯回線和50 Hz 動態(tài)磁滯回線,結(jié)果如圖2 和圖3 所示。

    表1 硅鋼片規(guī)格及測試條件Tab.1 Specifications and test conditions of silicon steel sheets

    圖1 電工鋼片一維磁特性測量平臺Fig.1 Measuring platform for 1-D magnetic propertiesof electrical steel sheets

    圖2 5 Hz 準(zhǔn)靜態(tài)磁滯回線測量結(jié)果Fig.2 Measured quasi-static hysteresis loop (5 Hz)

    圖3 50 Hz 動態(tài)磁滯回線測量結(jié)果Fig.3 Measured dynamic hysteresis loop (50 Hz)

    2 靜態(tài)逆Preisach 磁滯模型的構(gòu)建

    本文采用離散形式的靜態(tài)逆Preisach 磁滯模型描述硅鋼片的磁滯特性,該模型采用一種一階回轉(zhuǎn)曲線(First Order Reversal Curves, FORCs)的解析計(jì)算方法,并基于少量準(zhǔn)靜態(tài)磁滯回線測量數(shù)據(jù)實(shí)現(xiàn)解析式的參數(shù)辨識,從而構(gòu)造出辨識逆Everett 函數(shù)所需的一階回轉(zhuǎn)曲線,并進(jìn)一步建立靜態(tài)逆Preisach 磁滯模型。解析一階回轉(zhuǎn)曲線構(gòu)造示意圖如圖4 所示:一條起始于準(zhǔn)靜態(tài)極限磁滯回線下降支的一階回轉(zhuǎn)曲線R-P-N-T,該曲線起點(diǎn)為回轉(zhuǎn)點(diǎn)R,極限磁滯回線上升支Ha和下降支Hd交匯于頂點(diǎn)T。該一階回轉(zhuǎn)曲線構(gòu)造時需首先確定極限磁滯回線在回轉(zhuǎn)點(diǎn)R 處寬度ΔHrev和R 點(diǎn)與T 點(diǎn)的垂直高度ΔBrev為

    圖4 解析一階回轉(zhuǎn)曲線構(gòu)造示意圖Fig.4 Construction ofanalytical first order reversal curve

    圖4 中一階回轉(zhuǎn)曲線上點(diǎn)P 所在磁通密度為BP,磁場強(qiáng)度HP的計(jì)算式為

    式中,ΔHout為極限磁滯回線在高度為BP處的寬度;x表征不同的P 點(diǎn)在一階回轉(zhuǎn)曲線上的相對位置,定義為

    a、b、c為待擬合的模型參數(shù)。當(dāng)滿足a>0,0<b<1,c>0 等條件時,可使得一階回轉(zhuǎn)曲線斜率始終為正且不會超出最大極限磁滯回線。同時,以上模型參數(shù)與回轉(zhuǎn)點(diǎn)R 在下降支上的位置有關(guān),此處定義位置系數(shù)β為

    式中,ΔBout為最大環(huán)磁滯回線的高度。

    參數(shù)a、b、c可表示為β的多項(xiàng)式,即

    式中,y為多項(xiàng)式系數(shù)向量。

    以上多項(xiàng)式系數(shù)通過準(zhǔn)靜態(tài)磁滯回線測量數(shù)據(jù)進(jìn)行辨識,由于對稱性,辨識過程僅需要磁滯回線的上升支。逆Everett 函數(shù)值可利用上升支一階回轉(zhuǎn)曲線Hforc通過式(7)計(jì)算。

    磁通密度幅值為Bm的同心磁滯回線上升支可通過逆Everett 函數(shù)進(jìn)行插值運(yùn)算得到。

    辨識程序中一共使用9 條磁滯回線,將每條磁滯回線上升支根據(jù)磁通密度均勻等分為50 個點(diǎn)?;谑剑?)、式(8)可得到由一階回轉(zhuǎn)曲線計(jì)算磁滯回線的數(shù)值關(guān)系,進(jìn)而建立適應(yīng)度函數(shù)F為

    式中,Hmeas、Hcal分別為磁滯回線磁場強(qiáng)度的測量值和計(jì)算值;Bki為辨識程序中第k條磁滯回線上第i個點(diǎn)的磁通密度值。采用遺傳算法基于適應(yīng)度函數(shù)進(jìn)行參數(shù)尋優(yōu),結(jié)果見表2。獲得全部參數(shù)后,采用該解析式模擬所需一階回轉(zhuǎn)曲線數(shù)據(jù),進(jìn)一步可計(jì)算逆Everett 函數(shù)E(bu,bv),通過解析一階回轉(zhuǎn)曲線辨識的逆Everett 函數(shù)如圖5 所示。

    表2 B30P105 型電工鋼片多項(xiàng)式參數(shù)辨識結(jié)果Tab.2 The identified polynomial parameters of B30P105 elctrical steel sheet

    圖5 通過解析一階回轉(zhuǎn)曲線辨識的逆Everett 函數(shù)Fig.5 Inverted Everett functionidentified by analytical first order reversal curves

    通過上述逆Everett 函數(shù)辨識結(jié)果,可計(jì)算出一個周期內(nèi)不同磁通密度下電工鋼片的靜態(tài)磁滯回線。B30P105 型取向硅鋼的準(zhǔn)靜態(tài)磁滯回線實(shí)驗(yàn)結(jié)果和模型計(jì)算結(jié)果對比如圖6 所示,可以看出,本文所構(gòu)建的逆Preisach 磁滯模型具有一定的全局意義,可準(zhǔn)確模擬不同磁通密度幅值下的準(zhǔn)靜態(tài)磁滯回線,為進(jìn)一步考慮動態(tài)磁滯特性的有限元計(jì)算提供了可靠的模型基礎(chǔ)。

    圖6 準(zhǔn)靜態(tài)磁滯回線測量值與計(jì)算值對比Fig.6 Comparison between measured and calculated quasi-static hysteresis loops

    3 考慮動態(tài)磁滯特性的時域有限元方法

    3.1 動態(tài)磁場的表征方法

    Bertotti 基于磁疇理論和統(tǒng)計(jì)學(xué)原理提出損耗分離理論,將鐵磁材料單位體積的總損耗Ptot分為磁滯損耗Phys、渦流損耗Peddy和異常損耗Pexc三項(xiàng)。即

    式中,σ為材料的電導(dǎo)率;d為硅鋼片厚度;S為橫截面積;G為無量綱耦合常數(shù),G=0.135 6;V0為與材料內(nèi)部磁性單元有關(guān)的統(tǒng)計(jì)參數(shù),同磁通密度幅值Bm相關(guān),可通過不同頻率下的磁滯回線實(shí)驗(yàn)數(shù)據(jù)擬合得到。

    磁場損耗W與B、H之間的關(guān)系為

    將式(11)代入式(10)可得動態(tài)磁場表達(dá)式為

    式中,Htot、Hhys、Heddy、Hexc分別為動態(tài)磁場總磁場強(qiáng)度、靜態(tài)磁滯磁場強(qiáng)度、渦流磁場強(qiáng)度、異常磁場強(qiáng)度;δB為符號函數(shù),δB=sign(dB/dt)=±1。然而,由于較強(qiáng)的非線性特征,上述磁場難以直接在有限元分析中進(jìn)行迭代耦合計(jì)算,需進(jìn)一步研究其有限元迭代形式。

    3.2 動態(tài)磁滯特性在有限元中的直接耦合形式

    基于后向差分法,總磁場強(qiáng)度的三個分量的時域離散形式為

    式中,t為時間步;Δt為步長間隔;vh為靜態(tài)磁滯微分磁阻率,定義為

    相關(guān)麥克斯韋磁場方程為

    結(jié)合方程式(13)~式(19)可得基于矢量磁位A的控制方程為

    式中,ve為等效磁阻率,表達(dá)式為

    3.3 基于等效磁阻率的固定點(diǎn)迭代策略

    固定點(diǎn)法通過將非線性磁滯關(guān)系分為線性和非線性兩部分,通過在迭代過程內(nèi)不斷修改非線性部分以達(dá)到收斂,在解決磁滯非線性問題中具有較好的收斂性。為此,本文考慮將固定點(diǎn)迭代技術(shù)引入考慮動態(tài)磁滯特性的有限元迭代過程之中,在考慮靜態(tài)磁滯特性的基礎(chǔ)上,引入考慮動態(tài)磁場分量影響的等效磁阻率,采用式(21)中等效磁阻率代替?zhèn)鹘y(tǒng)固定點(diǎn)法的微分磁阻率vFP。應(yīng)用固定點(diǎn)技術(shù),總磁場強(qiáng)度Htot和磁通密度B的非線性關(guān)系可表示為

    式中,R為類磁化余量。

    結(jié)合麥克斯韋方程,可得控制方程為

    在二維求解域Ω中,將方程式(23)展開為離散形式為

    式中,N為離散單元的基函數(shù)。將方程式(24)改寫為矩陣形式,即

    式中,S為剛度矩陣;A為磁矢位向量;Y為電流區(qū)域系數(shù)向量;I為繞組電流向量;G為余量R的旋度矩陣。該方法中,渦流損耗等效磁場和剩余損耗等效磁場表達(dá)在類磁化余量R的計(jì)算中。通過這種方式,不需要大幅改動有限元方程,使得算法流程更簡潔,同時,將動態(tài)特性考慮在余量中,便于在以后的工作中對渦流和剩余損耗計(jì)算模型的改進(jìn)與修正。而且,在每一時刻的迭代過程,僅需計(jì)算一次等效磁阻率ve和剛度矩陣S,顯著提升了計(jì)算效率。

    為保證固定點(diǎn)算法收斂,磁阻率v應(yīng)滿足式(26)的必要條件[25,27]。

    在傳統(tǒng)固定點(diǎn)法中,所有網(wǎng)格單元采用統(tǒng)一的磁阻率,即

    式中,vmin、vmax分別為磁滯回線上斜率最小值與最大值。該方法稱為GCM(global-coefficient method)。顯然,采用式(27)計(jì)算的磁阻率總能滿足算法收斂的必要條件。然而,GCM 采用統(tǒng)一的磁阻率導(dǎo)致計(jì)算時間過長,需采用更具效率的微分磁阻率進(jìn)行迭代,微分磁阻率vFP定義為

    式中,dHtot、dB分別為磁場強(qiáng)度和磁通密度的差分。

    對于標(biāo)量磁滯模型,二維電磁場量間的關(guān)系可通過取模值降階為一維問題。如圖7 所示,藍(lán)色實(shí)線表示靜態(tài)磁滯回線,藍(lán)色虛線表示動態(tài)磁滯曲線,曲線L 為固定點(diǎn)迭代過程中動態(tài)磁滯磁場Htot與磁通密度B的軌跡。由于磁場分量Heddy、Hexc的計(jì)算是基于磁通密度B的后向差分對時間的偏導(dǎo)數(shù),因此,該曲線L 起始于t時刻靜態(tài)磁滯回線上的點(diǎn)P1,止于t+1 時刻動態(tài)磁滯回線上點(diǎn)P3。顯然,由于渦流損耗和剩余損耗的影響,曲線L 的斜率明顯不同于動態(tài)磁滯回線的斜率,在磁通密度的零點(diǎn)附近,兩者相差很大。傳統(tǒng)的固定點(diǎn)迭代法采用式(28)來計(jì)算微分磁阻率vFP,由于只考慮到動態(tài)磁滯磁場的軌跡,其計(jì)算結(jié)果偏小,往往不滿足收斂的必要條件,在迭代計(jì)算中算法極易發(fā)散。相反,在本文所提式(21)中等效磁阻率ve的計(jì)算考慮了渦流損耗和剩余損耗的等效磁場,反映了鐵心動態(tài)磁滯損耗對迭代過程的影響大小,可準(zhǔn)確估計(jì)動態(tài)B-Htot軌跡L 的斜率,從而保證迭代過程穩(wěn)定收斂。在有限元程序中ve根據(jù)前兩個時間點(diǎn)的結(jié)果進(jìn)行估算,即

    圖7 動態(tài)磁滯收斂路徑示意圖Fig.7 Schematic diagramof dynamic hysteresis convergence path

    3.4 考慮動態(tài)磁滯特性的場-路耦合有限元分析

    鑒于電機(jī)、變壓器等大多數(shù)電氣設(shè)備工作于50 Hz正弦交流電壓激勵下,為反映動態(tài)磁場和電壓源的相互作用,需在時域有限元分析中進(jìn)行“場-路”耦合。由于愛潑斯坦方圈實(shí)質(zhì)可看作一臺等比的空載變壓器,可以滿足對變壓器鐵心損耗分析的需要。為驗(yàn)證本文所提方法的有效性,本文以愛潑斯坦方圈為例構(gòu)建了二維有限元計(jì)算模型。

    在模型的計(jì)算過程中,以50 Hz 正弦交流電壓作為激勵,從電路角度分析,不考慮勵磁線圈漏磁和匝間電容,則端口電路方程為

    式中,U為激勵電壓;Uin為繞組上的感應(yīng)電動勢;r為線路電阻;I為勵磁電流。對于一階三角形單元,Uin的計(jì)算式為

    式中,Φ為繞組磁通;Nc為線圈匝數(shù);Sc為繞組截面積;Ωc為電流區(qū)域;se為單元面積;lz為二維求解系統(tǒng)縱向深度;Ae為單元節(jié)點(diǎn)磁矢位。將式(31)代入式(30)得到

    式中,U為激勵電壓向量;r為支路電阻矩陣。采用Crank-Nicolson 差分格式,將式(25)、式(32)聯(lián)立求解。

    通過求解方程式(33)得到磁矢勢向量,然后根據(jù)B=?×A計(jì)算出每個單元的磁通密度大小。如果計(jì)算結(jié)果未達(dá)到收斂條件eB=∣Bi+1-Bi∣<Bε(相對殘差閾值一般取εB= 1× 1 0-5),則通過磁滯模型計(jì)算靜態(tài)磁滯磁場Hhys、渦流磁場Heddy、剩余損耗磁場Hexc,進(jìn)而計(jì)算總磁場強(qiáng)度Htot,并通過松弛因子λ對Htot進(jìn)行修正。

    式中,i為迭代次數(shù)。若磁通密度B相鄰迭代誤差滿足收斂條件,則判斷是否達(dá)到最大時間步,若t=tmax,則停止迭代,輸出全部計(jì)算結(jié)果,否則跳轉(zhuǎn)至下一時間步。整個有限元迭代過程如圖8 所示。

    圖8 考慮動態(tài)磁滯特性的有限元分析流程Fig.8 Finite element iterative calculation process considering dynamic hysteresis characteristics

    需要說明的是:不同于靜態(tài)磁滯磁場,動態(tài)磁場強(qiáng)度和磁通密度間的非線性關(guān)系不是一開始就確定的,由于渦流磁場和剩余損耗磁場不是直接通過有限元磁場方程表達(dá),而是在固定點(diǎn)迭代過程加入,這意味著在算法迭代過程中,非線性磁場關(guān)系不斷改變,由動態(tài)磁場分量帶來的數(shù)值擾動使得磁矢位A和磁通密度B計(jì)算結(jié)果變化劇烈,造成收斂困難,尤其在磁通密度的零點(diǎn)附近,dB/dHtot較大,磁場強(qiáng)度的計(jì)算誤差引起的數(shù)值振蕩更大。因此松弛因子λ的引入是十分必要的。

    4 結(jié)果與討論

    4.1 愛潑斯坦方圈動態(tài)磁滯回線有限元模擬

    為說明本文所提方法的有效性,在計(jì)算模型中,選取了愛潑斯坦方圈臂上的三角形單元作為磁滯回線與損耗的觀測點(diǎn),其二維網(wǎng)格剖分情況如圖9 所示,三角網(wǎng)格M 為進(jìn)行實(shí)驗(yàn)數(shù)據(jù)與模型計(jì)算數(shù)據(jù)的對比的參考單元。在50 Hz 正弦條件下,圖10a、圖11a 分別給出了激勵電壓幅值為Umax=9 V 和Umax=13 V 時M 單元處磁滯回線測量值與計(jì)算值的對比??梢钥闯?,再考慮動態(tài)磁滯特性后,本文所提方法具有較高的磁滯回線模擬精度。進(jìn)一步地,圖10b、圖11b 給出了相應(yīng)電壓激勵條件下愛潑斯坦方圈勵磁電流波形的模擬結(jié)果。由于考慮了動態(tài)磁場與磁滯效應(yīng)對電流波形的反饋?zhàn)饔?,本文所提方法可以較好地模擬勵磁電流時域波形。

    圖9 愛潑斯坦方圈有限元剖分Fig.9 Finite element mesh of Epstein frame

    圖10 Umax=9 V 電壓激勵下測量值與計(jì)算值對比Fig.10 Comparison between measured and calculated values under voltage excitation Umax=9 V

    圖11 Umax=13 V 電壓激勵下測量值與計(jì)算值對比Fig.11 Comparison between measured and calculated values under voltage excitation Umax=13 V

    為考察松弛因子λ對程序收斂性的影響,控制電壓激勵分別為Umax=9 V,調(diào)節(jié)松弛因子大小,記錄程序是否收斂,以及程序收斂情況下一周期內(nèi)的平均迭代次數(shù)和程序總計(jì)算時間。結(jié)果見表3,當(dāng)λ=0.8時,有限元迭代過程不收斂,在確保收斂的情況下,采用較大的松弛因子可減少每時步的平均迭代次數(shù)和總計(jì)算時間,提升計(jì)算效率,最短總計(jì)算時間為228 s。通常,為同時滿足穩(wěn)定性和計(jì)算速度的要求,設(shè)置λ=0.5~0.7 即可。

    表3 9 V 電壓激勵時不同松弛因子對迭代次數(shù)和總迭代時間影響Tab.3 Effects of different relaxations on the iterations and computation time with voltage excitation amplitude of 9 V

    4.2 磁通密度分布與瞬時損耗特性

    為進(jìn)一步檢驗(yàn)算法的有效性,本文對愛潑斯坦方圈的磁通密度分布特性與瞬時功率損耗特性進(jìn)行了分析與計(jì)算。在t=0.005 s 時的磁通密度分布計(jì)算結(jié)果如圖12 所示??梢钥闯?,主磁路上磁通密度分布較為均勻,側(cè)壁上磁通密度幅值計(jì)算結(jié)果與實(shí)驗(yàn)測量值基本一致;而在疊片的內(nèi)角處,磁通密度較大,在外角部分磁通密度較小。

    圖12 磁通密度分布(Umax=13 V、t=0.005 s)Fig.12 Distribution of magnetic flux density (Umax=13 V、t=0.005 s)

    根據(jù)4.1 節(jié)中磁通密度B和動態(tài)磁場Htot計(jì)算結(jié)果,本文進(jìn)一步計(jì)算了Umax=9 V 和Umax=13 V 兩種激勵大小在參考單元M 處一周期內(nèi)的瞬時損耗功率tp并與測量值進(jìn)行對比,結(jié)果如圖13 所示。在t=0.00 2 s 和t=0.012 s 附近,測量與計(jì)算的瞬時損耗達(dá)到極大值;在t=0.005 s 和t=0.015 s 附近,瞬時損耗達(dá)到極小值。損耗負(fù)值來源于可逆磁化的貢獻(xiàn)??梢钥闯觯疚乃岱椒ㄔ谀M瞬時功率方面具有較高的精度,瞬時損耗的最大誤差不超過6%。

    圖13 動態(tài)磁滯瞬時損耗測量值與計(jì)算值對比Fig.13 Comparison between measured and calculated instantaneous loss when considering dynamic hysteresis

    完成一個周期有限元計(jì)算后,積分得到每一個單元磁滯回線的面積并求和即可得到鐵心區(qū)域整體鐵心總損耗PFe為

    式中,f為激勵的頻率;ne為鐵心離散單元總數(shù)。圖14 所示為不同磁通密度幅值時的鐵心總損耗計(jì)算值和測量值大小??梢钥闯觯ㄟ^提出的有限元算法計(jì)算得到的一個周期內(nèi)的鐵心損耗與實(shí)驗(yàn)測量結(jié)果基本一致。

    圖14 鐵心損耗測量值與計(jì)算值對比Fig.14 Comparison between measured and calculated power losses

    5 結(jié)論

    為解決時域有限元分析中考慮鐵心材料的動態(tài)磁滯特性出現(xiàn)難以收斂的問題,本文提出了一種基于等效磁阻率的高效穩(wěn)定改進(jìn)固定點(diǎn)迭代策略,并編寫了相應(yīng)的有限元程序,實(shí)現(xiàn)了硅鋼片的動態(tài)磁滯特性進(jìn)行模擬計(jì)算,并通過實(shí)驗(yàn)測量結(jié)果對比,驗(yàn)證了該有限元算法的有效性與準(zhǔn)確性,具體貢獻(xiàn)如下:

    1)采用解析方法生成一階回轉(zhuǎn)曲線,辨識逆Everett 函數(shù),建立逆Preisach 磁滯模型,可準(zhǔn)確模擬電工鋼片靜態(tài)磁滯特性。

    2)采用固定點(diǎn)迭代技術(shù)耦合動態(tài)磁滯特性,通過分析渦流損耗和剩余損耗等效磁場對有限元迭代過程的影響,提出基于等效磁阻率的固定點(diǎn)策略,并在有限元迭代過程引入松弛因子,通過實(shí)驗(yàn)數(shù)據(jù)與計(jì)算結(jié)果的對比分析,驗(yàn)證了本文提出的方法可實(shí)現(xiàn)低頻動態(tài)磁滯磁場的準(zhǔn)確計(jì)算,并具有良好的計(jì)算速度和穩(wěn)定性。

    3)通過場-路耦合有限元算法與動態(tài)磁滯模型結(jié)合,可充分考慮鐵心總損耗對電路和磁場的反饋效應(yīng),實(shí)現(xiàn)瞬時損耗的準(zhǔn)確計(jì)算,最大計(jì)算誤差不超過6%,可滿足工程計(jì)算要求。

    猜你喜歡
    磁滯回線磁阻磁通
    軸向磁通電勵磁雙凸極電機(jī)及容錯運(yùn)行控制策略
    基于MATLAB處理大學(xué)物理實(shí)驗(yàn)數(shù)據(jù)探究
    永磁磁阻電動機(jī)的研究
    磁場強(qiáng)度波形畸變對交流磁滯回線形狀的影響
    基于LabVIEW的微型磁通門磁強(qiáng)計(jì)測試系統(tǒng)搭建
    基于磁通門原理的零磁通交直流電流傳感器
    高頻脈沖激勵下磁滯回線動態(tài)測量裝置的設(shè)計(jì)及分析
    巨磁阻電渦流傳感器設(shè)計(jì)
    基于FPGA的數(shù)字磁通計(jì)設(shè)計(jì)
    電測與儀表(2015年3期)2015-04-09 11:37:52
    四相開關(guān)磁阻電機(jī)的四電平DITC調(diào)速系統(tǒng)
    国产视频内射| 亚洲精品在线观看二区| 国产精品一及| а√天堂www在线а√下载| 国产中年淑女户外野战色| 国产高清三级在线| 国产毛片a区久久久久| 在线观看av片永久免费下载| 变态另类丝袜制服| 欧美日韩福利视频一区二区| 美女被艹到高潮喷水动态| 婷婷精品国产亚洲av| 亚洲自拍偷在线| 久久久久久久午夜电影| 精品一区二区三区视频在线| 免费人成在线观看视频色| 午夜两性在线视频| 五月玫瑰六月丁香| 午夜激情欧美在线| 99精品在免费线老司机午夜| 久久热精品热| 国产精品爽爽va在线观看网站| 99在线视频只有这里精品首页| 国产在视频线在精品| 少妇高潮的动态图| 婷婷色综合大香蕉| 亚洲国产精品成人综合色| 性欧美人与动物交配| 国产亚洲精品综合一区在线观看| 亚洲真实伦在线观看| 亚洲熟妇熟女久久| 在线免费观看不下载黄p国产 | 亚洲欧美精品综合久久99| 人人妻人人看人人澡| 国产高清视频在线播放一区| 久久精品综合一区二区三区| 精品无人区乱码1区二区| 人妻丰满熟妇av一区二区三区| 午夜两性在线视频| 国产av不卡久久| 欧美成狂野欧美在线观看| 最近最新中文字幕大全电影3| 国语自产精品视频在线第100页| 欧美在线黄色| 色精品久久人妻99蜜桃| 最近最新中文字幕大全电影3| 99视频精品全部免费 在线| 亚洲片人在线观看| 一边摸一边抽搐一进一小说| 国产又黄又爽又无遮挡在线| 又紧又爽又黄一区二区| 美女黄网站色视频| 一本一本综合久久| 看黄色毛片网站| 91午夜精品亚洲一区二区三区 | 久久久国产成人免费| 国产在线精品亚洲第一网站| 日韩免费av在线播放| 床上黄色一级片| 99在线视频只有这里精品首页| 啦啦啦观看免费观看视频高清| 亚洲精品在线美女| 免费高清视频大片| 男人和女人高潮做爰伦理| 淫妇啪啪啪对白视频| 亚洲电影在线观看av| 黄片小视频在线播放| a在线观看视频网站| 久久久久精品国产欧美久久久| 一个人看的www免费观看视频| 男人舔女人下体高潮全视频| 18禁裸乳无遮挡免费网站照片| 精品人妻一区二区三区麻豆 | 国产久久久一区二区三区| 亚洲av五月六月丁香网| 亚洲av美国av| 精品人妻熟女av久视频| 好看av亚洲va欧美ⅴa在| 精品人妻1区二区| 我要看日韩黄色一级片| 欧美丝袜亚洲另类 | 国产男靠女视频免费网站| av在线天堂中文字幕| 国产69精品久久久久777片| 国产伦在线观看视频一区| 国产男靠女视频免费网站| 免费观看人在逋| 婷婷丁香在线五月| 婷婷精品国产亚洲av| 亚洲激情在线av| 久久久久九九精品影院| 欧美一级a爱片免费观看看| 麻豆久久精品国产亚洲av| 久久人人爽人人爽人人片va | 亚洲 欧美 日韩 在线 免费| 色综合站精品国产| 久久草成人影院| 看十八女毛片水多多多| www.999成人在线观看| 久9热在线精品视频| 99精品在免费线老司机午夜| 精华霜和精华液先用哪个| 亚洲av熟女| 国产精品久久电影中文字幕| 小说图片视频综合网站| 免费高清视频大片| 男人狂女人下面高潮的视频| 国产老妇女一区| 深夜精品福利| 97超级碰碰碰精品色视频在线观看| 欧美乱色亚洲激情| eeuss影院久久| 久久亚洲真实| 亚洲美女搞黄在线观看 | 成人一区二区视频在线观看| 黄色日韩在线| 此物有八面人人有两片| 欧美精品国产亚洲| 亚洲成av人片免费观看| av在线老鸭窝| h日本视频在线播放| 久久精品国产自在天天线| 丰满人妻一区二区三区视频av| 亚洲欧美日韩卡通动漫| 三级毛片av免费| 国产精品一区二区免费欧美| 好看av亚洲va欧美ⅴa在| 在线观看av片永久免费下载| 亚洲,欧美精品.| 国产精品久久视频播放| 精品日产1卡2卡| 亚洲,欧美,日韩| 老熟妇仑乱视频hdxx| 99在线人妻在线中文字幕| 一边摸一边抽搐一进一小说| 欧美极品一区二区三区四区| 一本综合久久免费| 日本黄大片高清| 五月伊人婷婷丁香| 国产国拍精品亚洲av在线观看| 成人午夜高清在线视频| av女优亚洲男人天堂| 国产成人啪精品午夜网站| 欧美黄色片欧美黄色片| 超碰av人人做人人爽久久| 禁无遮挡网站| 特大巨黑吊av在线直播| 精品国内亚洲2022精品成人| 国产乱人伦免费视频| 中文字幕av成人在线电影| 国产91精品成人一区二区三区| 欧美在线一区亚洲| 最近最新中文字幕大全电影3| 国产爱豆传媒在线观看| 久久亚洲真实| 成人一区二区视频在线观看| 亚洲精品日韩av片在线观看| 99在线视频只有这里精品首页| 91九色精品人成在线观看| 国内揄拍国产精品人妻在线| 亚洲人与动物交配视频| 国产精品1区2区在线观看.| 日本精品一区二区三区蜜桃| 成年女人毛片免费观看观看9| 国内精品一区二区在线观看| 国产精品三级大全| 黄色丝袜av网址大全| 两个人的视频大全免费| 久久久久久久精品吃奶| a级一级毛片免费在线观看| 色在线成人网| 亚洲精品影视一区二区三区av| 午夜精品一区二区三区免费看| 精品人妻1区二区| 亚洲精华国产精华精| 成人一区二区视频在线观看| 精品一区二区免费观看| 夜夜夜夜夜久久久久| 一个人看的www免费观看视频| 欧美精品国产亚洲| 久久精品国产亚洲av天美| 国产高潮美女av| 超碰av人人做人人爽久久| 天堂影院成人在线观看| 成人av一区二区三区在线看| 久久久久久久精品吃奶| 中文字幕高清在线视频| 草草在线视频免费看| 别揉我奶头~嗯~啊~动态视频| 国产精品99久久久久久久久| 精品国产三级普通话版| 午夜激情欧美在线| 久久久久国产精品人妻aⅴ院| 日韩欧美三级三区| 日本黄色片子视频| 成人av一区二区三区在线看| 精品久久久久久久末码| av专区在线播放| 一级a爱片免费观看的视频| 男人的好看免费观看在线视频| 国产激情偷乱视频一区二区| 他把我摸到了高潮在线观看| 国产69精品久久久久777片| 精品不卡国产一区二区三区| 国产伦人伦偷精品视频| av在线老鸭窝| 久久亚洲精品不卡| 亚洲成人久久爱视频| 欧美不卡视频在线免费观看| 国产一区二区三区在线臀色熟女| av在线老鸭窝| 深夜a级毛片| 午夜福利视频1000在线观看| 激情在线观看视频在线高清| 欧美性猛交╳xxx乱大交人| 丰满人妻熟妇乱又伦精品不卡| 国内少妇人妻偷人精品xxx网站| 两个人视频免费观看高清| 亚洲欧美日韩高清在线视频| 黄色女人牲交| 亚洲一区高清亚洲精品| 国内少妇人妻偷人精品xxx网站| 69av精品久久久久久| 成人欧美大片| 免费av毛片视频| x7x7x7水蜜桃| 亚洲第一区二区三区不卡| 99国产综合亚洲精品| 国产又黄又爽又无遮挡在线| 少妇丰满av| 精品欧美国产一区二区三| 成人高潮视频无遮挡免费网站| 一二三四社区在线视频社区8| 一卡2卡三卡四卡精品乱码亚洲| 国产精品综合久久久久久久免费| 欧美成人a在线观看| 午夜日韩欧美国产| 丰满人妻一区二区三区视频av| 美女高潮的动态| 国产亚洲精品久久久久久毛片| 欧美性猛交╳xxx乱大交人| 九色国产91popny在线| 深夜a级毛片| 欧美色视频一区免费| 成人高潮视频无遮挡免费网站| 有码 亚洲区| 久久伊人香网站| 男人舔奶头视频| 高清毛片免费观看视频网站| 久久久久久久久久成人| 久久久精品欧美日韩精品| 九色国产91popny在线| 亚洲国产精品sss在线观看| 久久久国产成人免费| 白带黄色成豆腐渣| 噜噜噜噜噜久久久久久91| 中文字幕免费在线视频6| 久久热精品热| 欧美成狂野欧美在线观看| 国内精品一区二区在线观看| 深爱激情五月婷婷| 中文字幕av在线有码专区| 亚洲av.av天堂| 最新在线观看一区二区三区| av在线蜜桃| 久久亚洲真实| 亚洲欧美日韩东京热| 亚洲国产欧美人成| 人妻丰满熟妇av一区二区三区| avwww免费| 在线看三级毛片| 精品国产亚洲在线| 午夜福利在线观看吧| 91狼人影院| 亚洲精品色激情综合| 老司机福利观看| 91午夜精品亚洲一区二区三区 | 欧美另类亚洲清纯唯美| 最好的美女福利视频网| 亚洲成人久久爱视频| 国产精品久久久久久久电影| 日本 欧美在线| 国产伦在线观看视频一区| 少妇的逼好多水| 国产av一区在线观看免费| 91av网一区二区| 无人区码免费观看不卡| 国产高清视频在线播放一区| 免费观看人在逋| 精品久久久久久久久久久久久| 一边摸一边抽搐一进一小说| 中文字幕免费在线视频6| 国语自产精品视频在线第100页| 亚洲综合色惰| 欧美乱妇无乱码| 窝窝影院91人妻| 午夜福利欧美成人| 国产69精品久久久久777片| 久久久久久九九精品二区国产| www.www免费av| 欧美日韩瑟瑟在线播放| 日本免费a在线| 午夜免费激情av| 成人国产综合亚洲| 国产伦人伦偷精品视频| 国产成+人综合+亚洲专区| 少妇的逼好多水| 十八禁国产超污无遮挡网站| 他把我摸到了高潮在线观看| 又爽又黄a免费视频| 国产一区二区在线观看日韩| 特大巨黑吊av在线直播| 97人妻精品一区二区三区麻豆| 欧美xxxx性猛交bbbb| 国产一区二区三区视频了| 嫩草影院入口| 欧美不卡视频在线免费观看| 偷拍熟女少妇极品色| 久久国产精品人妻蜜桃| 香蕉av资源在线| aaaaa片日本免费| 国产爱豆传媒在线观看| 欧美成狂野欧美在线观看| 脱女人内裤的视频| 国产高清视频在线播放一区| 亚洲中文字幕一区二区三区有码在线看| 欧美日韩综合久久久久久 | 看片在线看免费视频| 免费看a级黄色片| 色播亚洲综合网| 亚洲人成网站高清观看| 国产av在哪里看| 成人毛片a级毛片在线播放| 琪琪午夜伦伦电影理论片6080| 成熟少妇高潮喷水视频| 免费在线观看影片大全网站| 丝袜美腿在线中文| 日本一本二区三区精品| 少妇人妻一区二区三区视频| 国产亚洲精品久久久com| 特级一级黄色大片| 精品午夜福利在线看| 久久精品国产99精品国产亚洲性色| 两个人的视频大全免费| 欧美激情在线99| 男女那种视频在线观看| 丁香欧美五月| 别揉我奶头~嗯~啊~动态视频| 在线十欧美十亚洲十日本专区| 欧美午夜高清在线| 无人区码免费观看不卡| 日韩国内少妇激情av| 欧美黑人欧美精品刺激| 很黄的视频免费| 国语自产精品视频在线第100页| 免费高清视频大片| 中出人妻视频一区二区| av天堂中文字幕网| 色哟哟·www| 欧美高清性xxxxhd video| 91字幕亚洲| 国产精华一区二区三区| 白带黄色成豆腐渣| 观看美女的网站| 看黄色毛片网站| 十八禁网站免费在线| 亚洲美女搞黄在线观看 | 午夜久久久久精精品| 国产三级在线视频| 日韩 亚洲 欧美在线| 亚洲电影在线观看av| 久久亚洲真实| 日本熟妇午夜| 最近视频中文字幕2019在线8| 一本精品99久久精品77| 色在线成人网| 男女视频在线观看网站免费| 免费在线观看日本一区| 亚洲国产精品999在线| 国产亚洲av嫩草精品影院| 91九色精品人成在线观看| 久久精品国产清高在天天线| 三级毛片av免费| 亚洲中文字幕日韩| 午夜福利视频1000在线观看| 1000部很黄的大片| 最近在线观看免费完整版| 中文亚洲av片在线观看爽| 日本a在线网址| 国产成+人综合+亚洲专区| 午夜亚洲福利在线播放| 老熟妇乱子伦视频在线观看| 91久久精品国产一区二区成人| 国产精品一区二区三区四区免费观看 | 9191精品国产免费久久| 乱人视频在线观看| 亚洲va日本ⅴa欧美va伊人久久| 亚洲av第一区精品v没综合| 亚洲va日本ⅴa欧美va伊人久久| 嫩草影院新地址| 黄色一级大片看看| 国模一区二区三区四区视频| 国产一区二区三区视频了| 麻豆国产av国片精品| 亚洲专区中文字幕在线| 久久99热6这里只有精品| 国产午夜精品久久久久久一区二区三区 | 美女大奶头视频| 国产精品久久久久久久久免 | 欧美日韩福利视频一区二区| 久久久久亚洲av毛片大全| 免费在线观看成人毛片| 国产三级在线视频| 欧美极品一区二区三区四区| 免费av毛片视频| 国产在线精品亚洲第一网站| 国产三级黄色录像| 好看av亚洲va欧美ⅴa在| 精品免费久久久久久久清纯| 国产精品一区二区三区四区免费观看 | 国产午夜精品久久久久久一区二区三区 | 少妇丰满av| 久久伊人香网站| 亚洲激情在线av| 天堂网av新在线| 麻豆国产97在线/欧美| 国产激情偷乱视频一区二区| 乱人视频在线观看| 日韩精品中文字幕看吧| 一边摸一边抽搐一进一小说| 欧美日韩黄片免| 熟女人妻精品中文字幕| 色在线成人网| 国产久久久一区二区三区| 美女黄网站色视频| 99久久精品一区二区三区| 一个人免费在线观看的高清视频| 国产探花极品一区二区| 精品久久久久久久末码| 人妻夜夜爽99麻豆av| 在线免费观看不下载黄p国产 | 亚洲性夜色夜夜综合| 亚洲,欧美,日韩| 亚洲av熟女| 国产精品1区2区在线观看.| 免费观看的影片在线观看| 麻豆久久精品国产亚洲av| 午夜激情欧美在线| 国产野战对白在线观看| 亚洲avbb在线观看| eeuss影院久久| 美女高潮的动态| av天堂中文字幕网| h日本视频在线播放| 国产一区二区在线av高清观看| 亚洲第一电影网av| 国产一区二区在线观看日韩| 欧美色欧美亚洲另类二区| 999久久久精品免费观看国产| 欧美zozozo另类| 亚洲一区二区三区不卡视频| 亚洲国产精品久久男人天堂| 久久精品人妻少妇| 久久午夜亚洲精品久久| 国产精品亚洲一级av第二区| 夜夜爽天天搞| 国产探花极品一区二区| 午夜福利免费观看在线| 99精品在免费线老司机午夜| avwww免费| 国产单亲对白刺激| 搡女人真爽免费视频火全软件 | av在线天堂中文字幕| 国内精品一区二区在线观看| 国模一区二区三区四区视频| 91久久精品电影网| 又爽又黄无遮挡网站| 男女做爰动态图高潮gif福利片| 成人美女网站在线观看视频| av在线天堂中文字幕| h日本视频在线播放| 男人和女人高潮做爰伦理| 99国产极品粉嫩在线观看| 成熟少妇高潮喷水视频| 波多野结衣高清作品| 国产亚洲精品av在线| 男女视频在线观看网站免费| 不卡一级毛片| 色视频www国产| 久久婷婷人人爽人人干人人爱| 免费在线观看日本一区| 我要看日韩黄色一级片| 亚洲av二区三区四区| 国产欧美日韩一区二区精品| 精品久久久久久,| 级片在线观看| 午夜福利免费观看在线| 免费无遮挡裸体视频| 好男人在线观看高清免费视频| 亚洲精品一区av在线观看| 校园春色视频在线观看| 99在线人妻在线中文字幕| 国产精品99久久久久久久久| 亚洲欧美清纯卡通| 最好的美女福利视频网| 国产在视频线在精品| 精品久久久久久久久av| 日本熟妇午夜| 久久久色成人| 深夜a级毛片| 国产又黄又爽又无遮挡在线| 我要看日韩黄色一级片| 小说图片视频综合网站| 非洲黑人性xxxx精品又粗又长| 男人舔女人下体高潮全视频| 欧美潮喷喷水| 亚洲欧美日韩卡通动漫| 12—13女人毛片做爰片一| 亚洲精品456在线播放app | 欧美成狂野欧美在线观看| 午夜免费成人在线视频| 精品人妻一区二区三区麻豆 | 神马国产精品三级电影在线观看| 国产精品久久电影中文字幕| 精品久久国产蜜桃| 成人高潮视频无遮挡免费网站| 在现免费观看毛片| 丁香欧美五月| 亚洲三级黄色毛片| 国产高清有码在线观看视频| av在线蜜桃| 久久久国产成人精品二区| 久久中文看片网| 韩国av一区二区三区四区| 性色avwww在线观看| 亚洲av成人不卡在线观看播放网| 日本在线视频免费播放| 中文字幕熟女人妻在线| 小说图片视频综合网站| 免费在线观看亚洲国产| 我的女老师完整版在线观看| 国产亚洲精品av在线| 18禁裸乳无遮挡免费网站照片| 婷婷精品国产亚洲av在线| 日韩欧美在线二视频| a级毛片免费高清观看在线播放| 国产乱人伦免费视频| 日韩 亚洲 欧美在线| 男女做爰动态图高潮gif福利片| a级一级毛片免费在线观看| 日韩av在线大香蕉| 亚洲 国产 在线| 欧美国产日韩亚洲一区| 午夜福利免费观看在线| 久久香蕉精品热| 两人在一起打扑克的视频| 在线免费观看的www视频| 亚洲,欧美,日韩| 亚洲欧美清纯卡通| 深夜a级毛片| 中文字幕久久专区| 91久久精品国产一区二区成人| 老熟妇仑乱视频hdxx| 亚洲精品乱码久久久v下载方式| 国产亚洲欧美在线一区二区| 波野结衣二区三区在线| 欧美色欧美亚洲另类二区| 免费电影在线观看免费观看| 熟妇人妻久久中文字幕3abv| 好男人在线观看高清免费视频| 婷婷色综合大香蕉| a级一级毛片免费在线观看| 最近最新免费中文字幕在线| 亚洲精品一区av在线观看| 国产色爽女视频免费观看| 亚洲av第一区精品v没综合| 精品一区二区三区人妻视频| 精品久久久久久成人av| 露出奶头的视频| 全区人妻精品视频| 国产精品久久久久久久久免 | 一本精品99久久精品77| 我要搜黄色片| 国产高清视频在线播放一区| 欧美成人免费av一区二区三区| 国产一区二区激情短视频| 在线观看66精品国产| 夜夜夜夜夜久久久久| 亚洲熟妇中文字幕五十中出| 伦理电影大哥的女人| 免费在线观看亚洲国产| 免费电影在线观看免费观看| 欧美绝顶高潮抽搐喷水| 日韩有码中文字幕| 色综合亚洲欧美另类图片| 欧美乱色亚洲激情| 久久久国产成人免费| 久久久成人免费电影| 欧美+日韩+精品| 麻豆国产av国片精品| 99久久精品热视频| 欧美+日韩+精品| 麻豆国产av国片精品| 久久精品国产亚洲av香蕉五月| 欧美日韩国产亚洲二区| 精品欧美国产一区二区三| 久久精品国产亚洲av香蕉五月| 在线播放无遮挡| 欧美色视频一区免费| 久久久成人免费电影| 熟妇人妻久久中文字幕3abv| 熟女人妻精品中文字幕| 高清日韩中文字幕在线| 伊人久久精品亚洲午夜| 在线免费观看不下载黄p国产 |