一种考虑动态磁滞效应的高效稳定时域有限元计算方法

魏 鹏1 陈 龙1 贲 彤1 井立兵2 张 献3

(1.三峡大学电气与新能源学院 宜昌 443002 2.湖北省微电网工程技术研究中心(三峡大学) 宜昌 443002 3.省部共建电工装备可靠性与智能化国家重点实验室(河北工业大学) 天津 300130)

摘要 在时域有限元分析中,考虑铁心磁滞特性与损耗等效磁场对磁场分布的反馈作用,对精细化模拟电工装备的暂态电磁特性具有重要意义。针对耦合动态磁滞特性有限元分析中,存在计算时间长、收敛不稳定等问题,该文提出一种基于等效磁阻率的改进固定点迭代策略,并应用于耦合Preisach 磁滞模型的时域有限元分析中。首先,鉴于逆Preisach 模型易于数值实现且具有较高的计算精度,铁心材料的磁滞行为由逆Preisach 模型进行描述,并基于解析一阶回转曲线对模型进行了参数辨识;其次,通过分析动态磁场分量对收敛性的影响,提出基于等效磁阻率的固定点迭代策略,并在有限元迭代过程中引入松弛因子以加强算法的收敛速度与稳定性;最后,以爱泼斯坦方圈为例,计算铁心中的磁场分布特性与瞬时损耗分布特性并进行实验验证。结果表明,该文所提动态磁滞有限元计算方法在具有较高精度的同时具有较好的收敛速度及稳定性。

关键词:时域有限元 固定点技术 磁滞建模 铁心损耗

0 引言

时域有限元方法广泛应用于电工装备的运行状态预测与损耗评估之中[1],准确预测电机、变压器等电工装备的时域暂态损耗分布特性不仅对降低设备能耗、改善局部过热具有重要意义,同时考虑损耗对暂态磁场求解的反馈作用,对进一步分析电工装备的瞬时功率平衡、暂态励磁特征等电磁特性至关重要[2-4]

目前,对于电工装备铁心损耗的有限元计算方法可大致分为两类:一是基于场计算结果的后处理方法,该类方法利用Steinmetz 损耗模型或Bertotti统计学损耗分离理论对仅单值非线性磁场的计算结果进行后处理,得到每一个单元的损耗,最终进行积分得到电工装备整体的损耗[5-7]。这一类方法主要用于对电机或变压器的损耗进行频域分析,并通过傅里叶变换来分析含有谐波激励时的铁心损耗。Ansys 公司D.Lin 等基于等效椭圆环思想,将上述方法进行了时域瞬时功率损耗计算扩展[8]。该方法的优势在于在损耗计算过程中无需迭代,在磁路简单的结构中具有较好的精度,但由于没有考虑铁心磁滞、涡流损耗、剩余损耗等效动态磁场对铁心中磁场分布的反馈作用,在电机等复杂磁路电工装备的损耗计算中将带来较大误差。另一类铁心损耗计算方法则是在时域有限元分析中,考虑动态磁滞效应对每一瞬时的功率损耗进行精确求解[9-10]。这一类方法可最大限度地模拟铁心的真实磁化过程,在处理含有磁通密度波形畸变工况下的损耗计算问题时具有较高精度。同时,在考虑铁心动态磁滞效应后,铁心中动态磁滞回线波形的精确模拟可在时域中对励磁电流的波形进行实时反馈,进而可对电机或变压器中暂态运行特性进行有效预测。因此,电工装备损耗的精确预测与运行状态的有效评估都离不开在时域有限元分析时对动态磁滞特性的影响的考虑。

然而,目前在有限元分析中常用的磁滞模型,如J-A 模型[11-13]、Preisach 模型[14-17]等均是与磁化速率无关的磁滞模型,并不能直接考虑电机或变压器叠片宏观涡流损耗与剩余损耗等动态因素对铁心磁场分布的影响。为此,需要进一步将二者的等效动态磁场考虑到有限元的计算方程之中。传统方法是将动态效果等效为电导率乘以磁矢位A 对时间的一阶导数,但该种等效仅适合均匀实心区域涡流路径的求解,并不适合叠片铁心结构[18-19]。对于铁心叠片结构,进行实际三维建模将会带来无法承受的计算成本。为此,E.Dlala 等提出了一种耦合一维叠片有限元方法的二维时域有限元计算方法[20-21]。在叠片维度,忽略边缘效应,认为磁通密度仅是叠片厚度方向上的函数,进而仅需要求解一维有限元方程,即可得到叠片的涡流损耗等效磁场。然而,该方法实质上是在二维有限元计算中嵌套了一维有限元程序,仍具有较高的求解难度。为了对上述求解方法进行简化,目前常采用Bertotti 模型的变形形式,即在耦合静态磁滞模型的过程中,通过耦合等效微分磁阻率来考虑涡流损耗磁场与剩余损耗磁场的效果,该方法避免了一维有限元的求解,在一定程度上大大简化了计算,但在每次迭代过程中仍需更新磁阻率和刚度矩阵,导致计算效率低下[22-23]。因此,在进一步求解含有等效磁阻率的方程时,需要进一步耦合非线性磁场计算中的固定点迭代技术。在仅考虑静态磁滞模型的有限元计算中,固定点迭代方法已经取得了较好的收敛效果与计算速度[24-26],并已在目前的商业软件中取得了较好的计算效果。然而,在引入涡流损耗磁场与剩余损耗磁场后,收敛将难以得到保证,特别是在计算至磁通密度零点附近时间步时算法极易发散,分析其原因主要是涡流损耗、剩余损耗等效动态磁场的介入改变了固定点法中磁通密度与磁场强度的收敛过程,使得传统微分磁阻率的定义无法准确描述收敛路径的斜率,计算得到微分磁阻率过低导致考虑动态磁滞效应的有限元方程求解不收敛。

为解决在考虑动态磁滞特性的时域有限元计算难以收敛的问题,本文提出一种基于等效磁阻率的固定点迭代求解策略,进而提出一种耦合动态磁滞特性的时步有限元分析方法。首先,为了获得电工钢片的动态、静态磁滞特性,基于爱泼斯坦方圈法搭建了一维磁特性测量系统;其次,考虑到逆Preisach模型可以考虑磁化历史具有较高的计算精度,同时在耦合有限元计算时易于数值实现的特点,本文构建了基于参数化解析一阶回转曲线的逆 Preisach模型,并进行了参数辨识;再次,分析了涡流损耗和剩余损耗对算法收敛特性的影响,提出基于等效磁阻率的固定点迭代策略,并在有限元迭代过程中加入松弛因子保证稳定性;最后,以爱泼斯坦方圈为例,进行了考虑动态磁滞效应的时域有限元分析,求解了瞬时功率损耗特性,并进行了实验对比验证。结果验证了本文所提方法的准确性、高效性与稳定性。

1 准静态与动态磁滞特性测量

为建立磁滞模型及验证有限元计算准确性,基于 IEC 60404—2 标准建立了一维磁性能测量平台,对牌号为B30P105 的取向硅钢片进行准静态与动态磁滞特性测量。电工钢片一维磁特性测量平台如图1 所示,该平台包括宽频功率放大器、电压前置放大器、爱泼斯坦方圈、多功能数据采集卡和 LabVIEW 控制系统。本实验采用的硅钢片规格及测试条件见表1。分别设定正弦交流电压激励的频率为5 Hz 与50 Hz,测量得到准静态磁滞回线和50 Hz 动态磁滞回线,结果如图2 和图3 所示。

表1 硅钢片规格及测试条件
Tab.1 Specifications and test conditions of silicon steel sheets

硅钢牌号 长度/mm 宽度/mm 厚度/mm 片数B30P105 300 30 0.275 6 20

图1 电工钢片一维磁特性测量平台
Fig.1 Measuring platform for 1-D magnetic propertiesof electrical steel sheets

图2 5 Hz 准静态磁滞回线测量结果
Fig.2 Measured quasi-static hysteresis loop (5 Hz)

图3 50 Hz 动态磁滞回线测量结果
Fig.3 Measured dynamic hysteresis loop (50 Hz)

2 静态逆Preisach 磁滞模型的构建

本文采用离散形式的静态逆Preisach 磁滞模型描述硅钢片的磁滞特性,该模型采用一种一阶回转曲线(First Order Reversal Curves, FORCs)的解析计算方法,并基于少量准静态磁滞回线测量数据实现解析式的参数辨识,从而构造出辨识逆Everett 函数所需的一阶回转曲线,并进一步建立静态逆Preisach 磁滞模型。解析一阶回转曲线构造示意图如图4 所示:一条起始于准静态极限磁滞回线下降支的一阶回转曲线R-P-N-T,该曲线起点为回转点R,极限磁滞回线上升支Ha 和下降支Hd 交汇于顶点T。该一阶回转曲线构造时需首先确定极限磁滞回线在回转点R 处宽度ΔHrev 和R 点与T 点的垂直高度ΔBrev

图4 解析一阶回转曲线构造示意图
Fig.4 Construction ofanalytical first order reversal curve

图4 中一阶回转曲线上点P 所在磁通密度为BP,磁场强度HP 的计算式为

式中,ΔHout为极限磁滞回线在高度为BP处的宽度;x 表征不同的P 点在一阶回转曲线上的相对位置,定义为

abc 为待拟合的模型参数。当满足a>0,0<b<1,c>0 等条件时,可使得一阶回转曲线斜率始终为正且不会超出最大极限磁滞回线。同时,以上模型参数与回转点R 在下降支上的位置有关,此处定义位置系数β

式中,ΔBout为最大环磁滞回线的高度。

参数abc 可表示为β 的多项式,即

式中,y 为多项式系数向量。

以上多项式系数通过准静态磁滞回线测量数据进行辨识,由于对称性,辨识过程仅需要磁滞回线的上升支。逆Everett 函数值可利用上升支一阶回转曲线Hforc 通过式(7)计算。

磁通密度幅值为Bm 的同心磁滞回线上升支可通过逆Everett 函数进行插值运算得到。

辨识程序中一共使用9 条磁滞回线,将每条磁滞回线上升支根据磁通密度均匀等分为50 个点。基于式(7)、式(8)可得到由一阶回转曲线计算磁滞回线的数值关系,进而建立适应度函数F

式中,HmeasHcal 分别为磁滞回线磁场强度的测量值和计算值; Bki 为辨识程序中第k 条磁滞回线上第i 个点的磁通密度值。采用遗传算法基于适应度函数进行参数寻优,结果见表2。获得全部参数后,采用该解析式模拟所需一阶回转曲线数据,进一步可计算逆Everett 函数E(bu,bv),通过解析一阶回转曲线辨识的逆Everett 函数如图5 所示。

表2 B30P105 型电工钢片多项式参数辨识结果
Tab.2 The identified polynomial parameters of B30P105 elctrical steel sheet

模型多项式参数 计算值y1 9.97 y2 7.58 y3 -5.78 y4 0.95 y5 1.91 y6 -4.91 y7 3.01 y8 0.01 y9 0.54

图5 通过解析一阶回转曲线辨识的逆Everett 函数
Fig.5 Inverted Everett functionidentified by analytical first order reversal curves

通过上述逆Everett 函数辨识结果,可计算出一个周期内不同磁通密度下电工钢片的静态磁滞回线。B30P105 型取向硅钢的准静态磁滞回线实验结果和模型计算结果对比如图6 所示,可以看出,本文所构建的逆Preisach 磁滞模型具有一定的全局意义,可准确模拟不同磁通密度幅值下的准静态磁滞回线,为进一步考虑动态磁滞特性的有限元计算提供了可靠的模型基础。

图6 准静态磁滞回线测量值与计算值对比
Fig.6 Comparison between measured and calculated quasi-static hysteresis loops

3 考虑动态磁滞特性的时域有限元方法

3.1 动态磁场的表征方法

Bertotti 基于磁畴理论和统计学原理提出损耗分离理论,将铁磁材料单位体积的总损耗Ptot 分为磁滞损耗Phys、涡流损耗Peddy 和异常损耗Pexc 三项。即

式中,σ 为材料的电导率;d 为硅钢片厚度;S 为横截面积;G 为无量纲耦合常数,G=0.135 6;V0 为与材料内部磁性单元有关的统计参数,同磁通密度幅值Bm 相关,可通过不同频率下的磁滞回线实验数据拟合得到。

磁场损耗WBH 之间的关系为

将式(11)代入式(10)可得动态磁场表达式为

式中,HtotHhysHeddyHexc 分别为动态磁场总磁场强度、静态磁滞磁场强度、涡流磁场强度、异常磁场强度;δB 为符号函数,δB=sign(dB/dt)=±1。然而,由于较强的非线性特征,上述磁场难以直接在有限元分析中进行迭代耦合计算,需进一步研究其有限元迭代形式。

3.2 动态磁滞特性在有限元中的直接耦合形式

基于后向差分法,总磁场强度的三个分量的时域离散形式为

式中,t 为时间步;Δt 为步长间隔;vh 为静态磁滞微分磁阻率,定义为

相关麦克斯韦磁场方程为

结合方程式(13)~式(19)可得基于矢量磁位A 的控制方程为

式中,ve 为等效磁阻率,表达式为

3.3 基于等效磁阻率的固定点迭代策略

固定点法通过将非线性磁滞关系分为线性和非线性两部分,通过在迭代过程内不断修改非线性部分以达到收敛,在解决磁滞非线性问题中具有较好的收敛性。为此,本文考虑将固定点迭代技术引入考虑动态磁滞特性的有限元迭代过程之中,在考虑静态磁滞特性的基础上,引入考虑动态磁场分量影响的等效磁阻率,采用式(21)中等效磁阻率代替传统固定点法的微分磁阻率 vFP。应用固定点技术,总磁场强度Htot 和磁通密度B 的非线性关系可表示为

式中,R 为类磁化余量。

结合麦克斯韦方程,可得控制方程为

在二维求解域Ω 中,将方程式(23)展开为离散形式为

式中,N 为离散单元的基函数。将方程式(24)改写为矩阵形式,即

式中,S 为刚度矩阵;A 为磁矢位向量;Y 为电流区域系数向量;I 为绕组电流向量;G 为余量R 的旋度矩阵。该方法中,涡流损耗等效磁场和剩余损耗等效磁场表达在类磁化余量R 的计算中。通过这种方式,不需要大幅改动有限元方程,使得算法流程更简洁,同时,将动态特性考虑在余量中,便于在以后的工作中对涡流和剩余损耗计算模型的改进与修正。而且,在每一时刻的迭代过程,仅需计算一次等效磁阻率ve 和刚度矩阵S,显著提升了计算效率。

为保证固定点算法收敛,磁阻率v 应满足式(26)的必要条件[25,27]

在传统固定点法中,所有网格单元采用统一的磁阻率,即

式中,vminvmax 分别为磁滞回线上斜率最小值与最大值。该方法称为GCM(global-coefficient method)。显然,采用式(27)计算的磁阻率总能满足算法收敛的必要条件。然而,GCM 采用统一的磁阻率导致计算时间过长,需采用更具效率的微分磁阻率进行迭代,微分磁阻率vFP 定义为

式中,dHtot、dB 分别为磁场强度和磁通密度的差分。

对于标量磁滞模型,二维电磁场量间的关系可通过取模值降阶为一维问题。如图7 所示,蓝色实线表示静态磁滞回线,蓝色虚线表示动态磁滞曲线,曲线L 为固定点迭代过程中动态磁滞磁场Htot 与磁通密度B 的轨迹。由于磁场分量HeddyHexc 的计算是基于磁通密度B 的后向差分对时间的偏导数,因此,该曲线L 起始于t 时刻静态磁滞回线上的点P1,止于t+1 时刻动态磁滞回线上点P3。显然,由于涡流损耗和剩余损耗的影响,曲线L 的斜率明显不同于动态磁滞回线的斜率,在磁通密度的零点附近,两者相差很大。传统的固定点迭代法采用式(28)来计算微分磁阻率 vFP,由于只考虑到动态磁滞磁场的轨迹,其计算结果偏小,往往不满足收敛的必要条件,在迭代计算中算法极易发散。相反,在本文所提式(21)中等效磁阻率ve 的计算考虑了涡流损耗和剩余损耗的等效磁场,反映了铁心动态磁滞损耗对迭代过程的影响大小,可准确估计动态B-Htot 轨迹L 的斜率,从而保证迭代过程稳定收敛。在有限元程序中ve根据前两个时间点的结果进行估算,即

图7 动态磁滞收敛路径示意图
Fig.7 Schematic diagramof dynamic hysteresis convergence path

3.4 考虑动态磁滞特性的场-路耦合有限元分析

鉴于电机、变压器等大多数电气设备工作于50 Hz正弦交流电压激励下,为反映动态磁场和电压源的相互作用,需在时域有限元分析中进行“场-路”耦合。由于爱泼斯坦方圈实质可看作一台等比的空载变压器,可以满足对变压器铁心损耗分析的需要。为验证本文所提方法的有效性,本文以爱泼斯坦方圈为例构建了二维有限元计算模型。

在模型的计算过程中,以50 Hz 正弦交流电压作为激励,从电路角度分析,不考虑励磁线圈漏磁和匝间电容,则端口电路方程为

式中,U 为激励电压;Uin 为绕组上的感应电动势;r 为线路电阻;I 为励磁电流。对于一阶三角形单元,Uin 的计算式为

式中,Φ 为绕组磁通;Nc 为线圈匝数;Sc 为绕组截面积;Ωc 为电流区域;se 为单元面积;lz 为二维求解系统纵向深度;Ae 为单元节点磁矢位。将式(31)代入式(30)得到

式中,U 为激励电压向量;r 为支路电阻矩阵。采用Crank-Nicolson 差分格式,将式(25)、式(32)联立求解。

通过求解方程式(33)得到磁矢势向量,然后根据B=∇×A 计算出每个单元的磁通密度大小。如果计算结果未达到收敛条件 eB=∣Bi+1-Bi∣< Bε(相对残差阈值一般取 εB = 1× 1 0-5),则通过磁滞模型计算静态磁滞磁场Hhys、涡流磁场Heddy、剩余损耗磁场Hexc,进而计算总磁场强度Htot,并通过松弛因子λHtot 进行修正。

式中,i 为迭代次数。若磁通密度B 相邻迭代误差满足收敛条件,则判断是否达到最大时间步,若t=tmax,则停止迭代,输出全部计算结果,否则跳转至下一时间步。整个有限元迭代过程如图8 所示。

图8 考虑动态磁滞特性的有限元分析流程
Fig.8 Finite element iterative calculation process considering dynamic hysteresis characteristics

需要说明的是:不同于静态磁滞磁场,动态磁场强度和磁通密度间的非线性关系不是一开始就确定的,由于涡流磁场和剩余损耗磁场不是直接通过有限元磁场方程表达,而是在固定点迭代过程加入,这意味着在算法迭代过程中,非线性磁场关系不断改变,由动态磁场分量带来的数值扰动使得磁矢位A 和磁通密度B 计算结果变化剧烈,造成收敛困难,尤其在磁通密度的零点附近,dB/dHtot 较大,磁场强度的计算误差引起的数值振荡更大。因此松弛因子λ 的引入是十分必要的。

4 结果与讨论

4.1 爱泼斯坦方圈动态磁滞回线有限元模拟

为说明本文所提方法的有效性,在计算模型中,选取了爱泼斯坦方圈臂上的三角形单元作为磁滞回线与损耗的观测点,其二维网格剖分情况如图9 所示,三角网格M 为进行实验数据与模型计算数据的对比的参考单元。在50 Hz 正弦条件下,图10a、图11a 分别给出了激励电压幅值为 Umax=9 V 和Umax=13 V 时M 单元处磁滞回线测量值与计算值的对比。可以看出,再考虑动态磁滞特性后,本文所提方法具有较高的磁滞回线模拟精度。进一步地,图10b、图11b 给出了相应电压激励条件下爱泼斯坦方圈励磁电流波形的模拟结果。由于考虑了动态磁场与磁滞效应对电流波形的反馈作用,本文所提方法可以较好地模拟励磁电流时域波形。

图9 爱泼斯坦方圈有限元剖分
Fig.9 Finite element mesh of Epstein frame

图10 Umax=9 V 电压激励下测量值与计算值对比
Fig.10 Comparison between measured and calculated values under voltage excitation Umax=9 V

图11 Umax=13 V 电压激励下测量值与计算值对比
Fig.11 Comparison between measured and calculated values under voltage excitation Umax=13 V

为考察松弛因子λ 对程序收敛性的影响,控制电压激励分别为Umax=9 V,调节松弛因子大小,记录程序是否收敛,以及程序收敛情况下一周期内的平均迭代次数和程序总计算时间。结果见表3,当λ=0.8时,有限元迭代过程不收敛,在确保收敛的情况下,采用较大的松弛因子可减少每时步的平均迭代次数和总计算时间,提升计算效率,最短总计算时间为228 s。通常,为同时满足稳定性和计算速度的要求,设置λ=0.5~0.7 即可。

表3 9 V 电压激励时不同松弛因子对迭代次数和总迭代时间影响
Tab.3 Effects of different relaxations on the iterations and computation time with voltage excitation amplitude of 9 V

松弛因子λ 平均迭代次数 一周期总计算时间/s 0.4 12 336 0.5 9 282 0.6 7 251 0.7 7 228 0.8 不收敛

4.2 磁通密度分布与瞬时损耗特性

为进一步检验算法的有效性,本文对爱泼斯坦方圈的磁通密度分布特性与瞬时功率损耗特性进行了分析与计算。在t=0.005 s 时的磁通密度分布计算结果如图12 所示。可以看出,主磁路上磁通密度分布较为均匀,侧壁上磁通密度幅值计算结果与实验测量值基本一致;而在叠片的内角处,磁通密度较大,在外角部分磁通密度较小。

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

根据4.1 节中磁通密度B 和动态磁场Htot 计算结果,本文进一步计算了Umax=9 V 和Umax=13 V 两种激励大小在参考单元M 处一周期内的瞬时损耗功率 tp 并与测量值进行对比,结果如图13 所示。在t=0.00 2 s 和t=0.012 s 附近,测量与计算的瞬时损耗达到极大值;在t=0.005 s 和t=0.015 s 附近,瞬时损耗达到极小值。损耗负值来源于可逆磁化的贡献。可以看出,本文所提方法在模拟瞬时功率方面具有较高的精度,瞬时损耗的最大误差不超过6%。

图13 动态磁滞瞬时损耗测量值与计算值对比
Fig.13 Comparison between measured and calculated instantaneous loss when considering dynamic hysteresis

完成一个周期有限元计算后,积分得到每一个单元磁滞回线的面积并求和即可得到铁心区域整体铁心总损耗PFe

式中,f 为激励的频率;ne 为铁心离散单元总数。图14 所示为不同磁通密度幅值时的铁心总损耗计算值和测量值大小。可以看出,通过提出的有限元算法计算得到的一个周期内的铁心损耗与实验测量结果基本一致。

图14 铁心损耗测量值与计算值对比
Fig.14 Comparison between measured and calculated power losses

5 结论

为解决时域有限元分析中考虑铁心材料的动态磁滞特性出现难以收敛的问题,本文提出了一种基于等效磁阻率的高效稳定改进固定点迭代策略,并编写了相应的有限元程序,实现了硅钢片的动态磁滞特性进行模拟计算,并通过实验测量结果对比,验证了该有限元算法的有效性与准确性,具体贡献如下:

1)采用解析方法生成一阶回转曲线,辨识逆Everett 函数,建立逆Preisach 磁滞模型,可准确模拟电工钢片静态磁滞特性。

2)采用固定点迭代技术耦合动态磁滞特性,通过分析涡流损耗和剩余损耗等效磁场对有限元迭代过程的影响,提出基于等效磁阻率的固定点策略,并在有限元迭代过程引入松弛因子,通过实验数据与计算结果的对比分析,验证了本文提出的方法可实现低频动态磁滞磁场的准确计算,并具有良好的计算速度和稳定性。

3)通过场-路耦合有限元算法与动态磁滞模型结合,可充分考虑铁心总损耗对电路和磁场的反馈效应,实现瞬时损耗的准确计算,最大计算误差不超过6%,可满足工程计算要求。

参考文献

[1] 刘刚, 荣世昌, 武卫革, 等.基于混合有限元法和降阶技术的油浸式变压器绕组2 维瞬态流-热耦合场分析[J].高电压技术, 2022, 48(5): 1695-1705.Liu Gang, Rong Shichang, Wu Weige, et al.Twodimensional transient flow-thermal coupling field analysis of oil-immersed transformer windings based on hybrid finite element method and reduced-order technology[J].High Voltage Engineering, 2022, 48(5):1695-1705.

[2] Cui Han, Ngo K D T.Transient core-loss simulation for ferrites with nonuniform field in SPICE[J].IEEE Transactions on Power Electronics, 2019, 34(1): 659-667.

[3] 李冰, 王泽忠, 刘恪, 等.特高压变压器直流偏磁对绕组电流的影响[J].电工技术学报, 2020, 35(7):1422-1431.Li Bing, Wang Zezhong, Liu Ke, et al.Research on winding current of UHV transformer under DC-bias[J].Transactions of China Electrotechnical Society, 2020,35(7): 1422- 1431.

[4] 陈龙龙, 魏晓光, 焦重庆, 等.混合式高压直流断路器分断过程电磁瞬态建模和测试[J].电工技术学报, 2021, 36(24): 5261-5271.Chen Longlong, Wei Xiaoguang, Jiao Chongqing, et al.Electromagnetic transient modeling and test of hybrid DC circuit breaker[J].Transactions of China Electrotechnical Society, 2021, 36(24): 5261-5271.

[5] 李永建, 闫鑫笑, 张长庚, 等.基于磁-热-流耦合模型的变压器损耗计算和热点预测[J].电工技术学报, 2020, 35(21): 4483-4491.Li Yongjian, Yan Xinxiao, Zhang Changgeng, et al.Numerical prediction of losses and local overheating in transformer windings based on magnetic-thermalfluid model[J].Transactions of China Electrotechnical Society, 2020, 35(21): 4483-4491.

[6] Yamazaki K, Fukushima N.Torque and loss calculation of rotating machines considering laminated cores using post 1-D analysis[J].IEEE Transactions on Magnetics, 2011, 47(5): 994-997.

[7] 迟青光, 张艳丽, 陈吉超, 等.非晶合金铁心损耗与磁致伸缩特性测量与模拟[J].电工技术学报,2021, 36(18): 3876-3883.Chi Qingguang, Zhang Yanli, Chen Jichao, et al.Measurement and modeling of lossand magnetostrictive properties for the amorphous alloy core[J].Transactions of China Electrotechnical Society, 2021, 36(18): 3876-3883.

[8] Lin D, Zhou P, Fu W N, et al.A dynamic core loss model for soft ferromagnetic and power ferrite materials in transient finite element analysis[J].IEEE Transactions on Magnetics, 2004, 40(2): 1318-1321.

[9] Zhou Ping, Lin Dingsheng, Lu Chuan, et al.A new algorithm to consider the effects of core losses on 3-D transient magnetic fields[J].IEEE Transactions on Magnetics, 2014, 50(2): 365-368.

[10] Lin D, Zhou P, Chen Q M, et al.The effects of steel lamination core losses on 3D transient magnetic fields[J].IEEE Transactions on Magnetics, 2010,46(8): 3539-3542.

[11] 赵志刚, 马习纹, 姬俊安.考虑频率效应的正弦及直流偏磁条件电工钢磁滞特性模拟及验证[J].中国电机工程学报, 2021, 41(23): 8178-8186.Zhao Zhigang, Ma Xiwen, Ji Junan.Simulation and verification on hysteresis characteristics of electrical steel under sinusoidal and DC bias conditions considering frequency effects[J].Proceedings of the CSEE, 2021, 41(23): 8178-8186.

[12] Xiao Meng, Zuo Yuxiang, Li Yang, et al.Core loss calculation of anode saturable reactor in damping oscillation state based on J-A theory[J].IEEE Transactions on Applied Superconductivity, 2021,31(8): 1-4.

[13] 刘任, 李琳, 王亚琦, 等.基于随机性与确定性混合优化算法的Jiles-Atherton 磁滞模型参数提取[J].电工技术学报, 2019, 34(11): 2260-2268.Liu Ren, Li Lin, Wang Yaqi, et al.Parameter extraction for Jiles-Atherton hysteresis model based on the hybrid technique of stochastic and deterministic optimization algorithm[J].Transactions of China Electrotechnical Society, 2019, 34(11): 2260-2268.

[14] 陈龙, 易琼洋, 贲彤, 等.全局优化算法在Preisach磁滞模型参数辨识问题中的应用与性能对比[J].电工技术学报, 2021, 36(12): 2585-2593, 2606.Chen Long, Yi Qiongyang, Ben Tong, et al.Application and performance comparison of global optimization algorithms in the parameter identification problems of the preisach hysteresis model[J].Transactions of China Electrotechnical Society, 2021, 36(12): 2585-2593, 2606.

[15] 赵小军, 刘小娜, 肖帆, 等.基于Preisach 模型的取向硅钢片直流偏磁磁滞及损耗特性模拟[J].电工技术学报, 2020, 35(9): 1849-1857.Zhao Xiaojun, Liu Xiaona, Xiao Fan, et al.Hysteretic and loss modeling of silicon steel sheet under the DC biased magnetization based on the preisach model[J].Transactions of China Electrotechnical Society, 2020,35(9): 1849-1857.

[16] Hussain S, Lowther D A.An efficient implementation of the classical preisach model[J].IEEE Transactions on Magnetics, 2018, 54(3): 1-4.

[17] 段娜娜, 徐伟杰, 李永建, 等.基于极限磁滞回线法的软磁复合材料磁特性模拟[J].电工技术学报,2018, 33(20): 4739-4745.Duan Nana, Xu Weijie, Li Yongjian, et al.Electromagnetic property modeling of the soft magnetic composite material based on the limiting loop method[J].Transactions of China Electrotechnical Society, 2018, 33(20): 4739-4745.

[18] 张长庚, 李永建, 杨庆新.一种考虑磁滞和微观涡流效应的电磁模拟方法[J].中国电机工程学报,2016, 36(21): 5966-5974, 6041.Zhang Changgeng, Li Yongjian, Yang Qingxin.An electromagnetic simulation method considering hysteresis and micro-eddy current effect[J].Proceedings of the CSEE, 2016, 36(21): 5966-5974, 6041.

[19] Yue Shuaichao, Anderson P I, Li Yongjian, et al.A modified inverse vector hysteresis model for nonoriented electrical steels considering anisotropy for FEA[J].IEEE Transactions on Energy Conversion,2021, 36(4): 3251-3260.

[20] Dlala E, Belahcen A, Arkkio A.Efficient magnetodynamic lamination model for twodimensional field simulation of rotating electrical machines[J].Journal of Magnetism and Magnetic Materials, 2008, 320(20): 1006-1010.

[21] Dlala E.Efficient algorithms for the inclusion of the preisach hysteresis model in nonlinear finite-element methods[J].IEEE Transactions on Magnetics, 2011,47(2): 395-408.

[22] Sadowski N, Batistela N J, Bastos J P A, et al.An inverse Jiles-Atherton model to take into account hysteresis in time-stepping finite-element calculations[J].IEEE Transactions on Magnetics, 2002, 38(2): 797-800.

[23] Fallah E, Moghani J S.A new approach for finiteelement modeling of hysteresis and dynamic effects[J].IEEE Transactions on Magnetics, 2006, 42(11): 3674-3681.

[24] Dlala E, Belahcen A, Arkkio A.A fast fixed-point method for solving magnetic field problems in media of hysteresis[J].IEEE Transactions on Magnetics,2008, 44(6): 1214-1217.

[25] Dlala E, Belahcen A, Arkkio A.Locally convergent fixed-point method for solving time-stepping nonlinear field problems[J].IEEE Transactions on Magnetics, 2007, 43(11): 3969-3975.

[26] Li Wei, Koh C S.Investigation of the vector Jiles-Atherton model and the fixed point method combined technique for time-periodic magnetic problems[J].IEEE Transactions on Magnetics, 2015, 51(4): 1-6.

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

An Efficient and Stable Time Domain Finite Element Method Considering Dynamic Hysteresis Effect

Wei Peng1 Chen Long1 Ben Tong1 Jing Libing2 Zhang Xian3

(1.College of Electrical Engineering and New Energy China Three Gorges University Yichang 443002 China 2.Hubei Provincial Engineering Technology Research Center for Microgrid China Three Gorges University Yichang 443002 China 3.State Key Laboratory of Reliability and Intelligence of Electrical Equipment Hebei University of Technology Tianjin 300130 China)

Abstract In the time-domain finite element analysis, considering the hysteresis properties of the iron core and the feedback effect of the loss equivalent magnetic field on the magnetic field distribution is of great significance for finely simulating the transient electromagnetic characteristics of electrical equipment.Aiming at the problems of long calculation time and divergence in the finite element analysis coupled with dynamic hysteresis characteristics, this paper proposes an improved fixed point iteration strategy based on equivalent reluctivity, and applies it to time domain finite element analysis implemented with Preisach model.With the modified reluctivity and the introduced relaxation factor, the finite element algorithm proposed in this paper has good convergence speed and stability with high accuracy.

Firstly, since numerical application of the inverse Preisach model is easy and has high accuracy, the hysteresis behavior of the core material is described by the inverse Preisach model.The Preisach model is established based on the analytical first order reversal curves, of which the parameters are identified by several symmetric hysteresis loops.Secondly, by analyzing the influence of the dynamic magnetic field component on the convergence,a fixedpoint iterative strategy based on equivalent reluctivity is proposed, and a relaxation factor is introduced in finite element iterations to enhance the convergence speed and stability of the algorithm.Finally, taking the Epstein frame as an example, experiments and finite element calculations are carried out on the magnetic field distribution and instantaneous loss characteristics of the iron core.There are 9 analytical polynomial parameters of B30P105 electrical steel sheet, identified by Genetic Algorithm.The static hysteresis loops simulated by the Preisach model is of high accuracy.In the 2-D dimensional computation model, the triangle located at the limb-yoke region is selected as an observation spot.The calculated dynamic hysteresis loops, exciting current and instantaneous loss are compared with measured ones, under voltage excitation amplitude of 9 V and 13 V, respectively.The maximum error of the instantaneous loss is less than 6%.After one period computation, the total iron loss is calculated by integrate the area of dynamic hysteresis loop in each element.The calculated iron losses under different flux density are in substantial agreement.In order to investigate the effect of relaxation factor on the convergence property of the program, various values of relaxation factor are chosen and the results show that the average number of iterations and computation time are inversely proportional to the magnitude of the relaxation factor.However, the iteration process doesn't converge with the relaxation factor set as 0.8.

The following conclusions can be drawn from the simulation analysis: (1) The Preisach hysteresis model based on analytical first order reversal curves can accurately simulate the static hysteresis characteristics of electrical steel sheets.(2) By incorporating the dynamic hysteresis effect into the improved fixed-point iteration process, the proposed finite element algorithm can realize precise calculation of low-frequency dynamic hysteresis magnetic field, with good computation speed and stability.(3) The error of the calculated instantaneous loss is small, which meets engineering application requirements.

Keywords:Time-domain finite element method, fixed-point technique,hysteresis modeling, iron losses

中图分类号:TM153

DOI: 10.19595/j.cnki.1000-6753.tces.221180

国家自然科学基金(52007102, 52207012, 51977147)和省部共建电工装备可靠性与智能化国家重点实验室(河北工业大学)开放课题基金(EERIKF2021015)资助项目。

收稿日期2022-06-21 改稿日期 2022-07-16

作者简介

魏 鹏 男,1996 年生,硕士研究生,研究方向为磁性材料磁特性、电磁场数值分析。

E-mail:m18805167298@163.com

陈 龙 男,1989 年生,博士,讲师,研究方向为磁性材料磁特性模拟、全局优化设计。

E-mail:chenlong@ctgu.edu.cn(通信作者)

(编辑 郭丽军)