瞬态电磁-温度场耦合计算中自适应时间步长研究张宇娇;汪振亮;徐彬昭;黄雄峰;智李【摘要】瞬态温度场响应时间比电磁场长,在采用有限元法求解瞬态电磁-温度场顺序耦合时,若采用统一计算时间步长会导致温度场求解次数过多而增加计算时间.该文提出自适应时间步长耦合方法,采用指数平滑法预测电磁-温度场耦合时间间隔,并在两个相邻的耦合时间节点间,根据响应特征值和预测-校正法计算电磁场、温度场最佳离散步长.与传统等步长耦合方法对比,电磁、温度场均使用最佳的离散步长,同时避免了温度场的过频计算,减少了计算时间.最后,以通电铜导环为例,采用自适应时间步长耦合计算铜导体在交流电流作用下的瞬态温升,与传统等步长耦合方法相比,计算相对误差在 1%以内,总计算时间减少 26.6%,证明了该方法的有效性.【期刊名称】《电工技术学报》【年(卷),期】2018(033)019【总页数】8页(P4468-4475)【关键词】有限元法;电磁-温度场耦合;自适应时间步长;指数平滑预测;响应特征值【作者】张宇娇;汪振亮;徐彬昭;黄雄峰;智李【作者单位】三峡大学梯级水电站运行与控制湖北省重点实验室宜昌 443002;三峡大学梯级水电站运行与控制湖北省重点实验室宜昌 443002;三峡大学梯级水电站运行与控制湖北省重点实验室宜昌 443002;三峡大学梯级水电站运行与控制湖北省重点实验室宜昌 443002;三峡大学梯级水电站运行与控制湖北省重点实验室宜昌 443002【正文语种】中文【中图分类】TM155在线圈发射装置、继电器、电机等设备中，运行时产生过高的温升将影响金属导体的导电性、绝缘材料的介电性能，甚至会使材料出现破坏性的热膨胀，严重影响设备的安全运行[1-8]。
近年来，国内外学者针对电磁装备开展了大量关于电磁与热耦合问题求解方法的研究[9-11]。
采用有限元法分析上述问题时，瞬变的大电流使得导体温度在较短时间内迅速增加，因而在计算过程中不仅要考虑电磁场对温度场的影响，还需考虑温度变化对电磁场材料导电性能的影响，为此需要建立瞬态电磁-温度场强耦合模型。
目前采用顺序耦合的方法求解瞬态电磁-温度场强耦合模型时，电磁场和温度场往往采用相同的时间步长[3-5]。
然而在求解过程中，温度场响应较电磁场慢，若采用相同的时间步长计算，使得温度场求解次数过于频繁，增加了求解时间。
针对两物理场顺序耦合时响应时间不同的问题，文献[12]根据现有软件，提出代码耦合的概念，对各物理场采用不同的时间离散策略，并在时间的耦合点上进行载荷的传递。
文献[13]针对线性系统采用不同时间步长计算，并提出时间积分算法进行了稳定性分析。
文献[14]提出在流-固耦合时，对流体区域和固体区域采用不同时间步长的DFMT-SCSS（sub-cycled conventional staggered scheme）算法，选取固体求解区域时间步长为流体求解区域的倍数。
文献[15]提出电磁场计算自适应时间步长计算方法。
文献[16]针对微波加热问题，采用自适应时间步长计算温度场。
文献[17]针对两物理场顺序耦合问题，采用自适应时间步长算法使两物理场获取相同的时间步长。
然而上述方法在计算两物理场耦合时，即使提到了自适应时间步长，但最终各物理场仍然采用等步长计算，无法体现各物理场时间常数不同且实时变化的特点，造成两物理场未获得最佳离散策略的缺陷。
针对上述研究缺陷，本文提出电磁场与温度场耦合计算中自适应时间步长概念，即各物理场采用自适应时间步长算法获取最佳离散时策略，并在预测的耦合时间节点上进行自适应耦合。
最后以通电铜导环热分析为例，在电流初始值为零时，对铜导块施加最大值为10 000 A、0 Hz交流电，计算0.1 s内铜导体温度变化。
与传统的电磁-温度场等步长计算方法相比，计算时间减少20%，误差在1%以内，证明了该方法的有效性。
如图1所示，在tn时刻，首先通过指数平滑法预测下一耦合时间点tn+1。
然后在两耦合时间节点间，通过预测-校正法和响应特征值计算电磁场和温度场最佳离散步长，得到tn ~tn+1时间段内电磁场和温度场最佳时间离散步长为和。
电磁场、温度场分别采用、计算至耦合时间节点。
设有一物理场量u在区域W和边界G上满足一阶物理方程，即在区域W内（1）在边界G上（2）式中，A为求解区域W内满足微分算子；B为求解边界G上满足微分算子。
对区域离散，采用伽辽金或最小位能原理等，将式（2）和式（3）改写成有限元格式，即（3）式中，C为阻尼矩阵；K为系数矩阵；P为载荷矩阵。
对式（3）采用两点循环公式求解得式中，；；；q为积分参数，q =0、1/2、1分别为向前差分公式、中心差分公式和向后差分公式。
采用数值方法求解式（4）时，当载荷P或计算结果u变化较快时，为保证稳定性，需选取较小时间步长进行计算。
而当P或u变化较慢时，为节省计算时间，可选取较长步长计算。
求解中若采用恒定时间步长，则无法根据P或u的变化调整步长大小，通常为了保证稳定性，统一采用较小时间步长，导致计算时间增加[18]。
为此提出时间步长自适应控制算法，即根据载荷变化率及计算结果自动调整时间步长[16,19]。
在满足计算精度及稳定性的前提下减少计算时间。
如图2所示，首先由载荷离散误差确定时间步长为Dt1n+1，再由计算稳定性可确定离散步长为Dt2n+1，由此得出Tn时刻最大时间步长Dtn+1为在较小的离散时间段t∈(tn-1,tn)内，不考虑载荷正负号变化，式（4）中载荷矩阵P采用线性插值，等效载荷采用斜坡加载，如图3所示。
由载荷离散产生的误差为[16]式中，为第n个时间步内产生的离散误差。
对式（6）采用梯形公式积分，得到离散误差近似为[15]根据式（7），载荷离散误差近似正比于Dt2，可将下一步步长计算分为如下两步[16]：1）步长预测。
根据第n步计算误差，预测第n+1步步长。
式中，为安全系数，=0.9；为允许最大误差。
2）步长校正。
判断当第n+1步时间步长所产生误差是否满足<。
如果不满足，则采用式（9）进行修正迭代计算，直至满足<。
式中，为安全系数，=0.75。
J. R. Hughes[20]提出根据响应特征值确定计算稳定时间步长Dtn+1的方法，定义Dtn+1l为振荡限制条件。
当Dtn+1l>>1时系统处于振荡状态。
为保证计算稳定性，可取最大步长满足（10）（11）式中，f为稳定系数，f<1；Δun为tn-1到tn时间段场量u的变化。
电磁-温度场自适应时间步长耦合整体流程如图4所示，主要包括以下四个步骤，其中步骤1和步骤2属于初始化部分。
1）参数初始化。
建立电磁场和温度场计算有限元模型，并给定材料参数及载荷。
2）耦合时间点确定。
首先电磁-温度场采用较小步长耦合三次，用来获取温度场触发热量Qpre的预测数据。
然后采用指数平滑法预测温度场触发热量，通过触发热量判断耦合时间节点。
3）电磁场自适应步长计算。
获取电磁场自适应时间步长，启动电磁场计算，当电磁场累积热量达到温度场触发热量时，启动温度场计算。
4）温度场自适应时间步长计算。
获取温度场最适应步长，将电磁场计算平均热功率作为载荷加至温度场，当温度场计算完毕后，更新电磁场节点温度，并进行下一时间步长的温度场触发热量计算。
反复迭代计算至最终时间。
对于瞬态电磁-温度场顺序耦合，本文提出根据累积热量判断电磁-温度场耦合时间节点。
假设tn时刻输入电流密度为Jn。
为控制两场耦合的非平衡误差[17]，根据材料电阻率随温度变化特点，确定电磁场计算过程中允许的温度变化最大值，再通过指数平滑法预测累积发热量Qpre，最后启动电磁场计算，当tn+1时刻电磁场累积热量达到Qpre时，tn+1即为电磁-温度场耦合时间点。
2.1.1 最大温度变化计算线性材料热功率与电流密度的关系为（13）式中，V为单元体积，m3；Pn为tn时刻发热功率，W；Jn为tn时刻电流密度，A/m2；为电阻率，W·m；a为电阻随温度变化率，W·m/℃；为0℃时电阻率。
当输入热功率为Pn，由温度变化引起的功率计算误差应满足当式（14）取等号时，可得出tn时刻温度允许变化的最大值。
2.1.2 温度场触发热量计算根据tn时刻允许最大温度变化，由tn、tn-1、tn-2时刻输入热量及温度变化预测tn+1时刻温度场触发热量。
1）计算各时刻输入热量及温度变化，见表1。
2）计算各时刻温度随热量变化率，见表2。
3）温度场触发热量计算。
采用指数平滑方法预测tn+1时刻温度场变化率Kn+1为[21]温度场触发热量为（16）2.1.3 耦合时间判断当tn+1时刻电磁场累积热量满足以下条件之一时，tn+1为电磁-温度场耦合时间点，启动温度场计算。
1）单元最大热量变化达到单元预测热量阈值时自动启动温度场计算，即（17）2）总体热量变化达到预测阈值时自动启动温度场计算，即（18）式中，为自然系数，控制计算的精度。
根据似稳条件忽略位移电流，时变电磁场方程可写为式（19）形式的矢量磁位方程。
式中，W1为导体区域；W2为非导体区域；A为矢量磁位，Wb/m；V为电位，V；s为电导率，S/m；m为磁导率，H/m；Js为源电流密度，A/m2；为导体在磁场中运动速度，m。
采用伽辽金法将式（19）写成有限元格式[22]，即（20）式中，；R为电磁场阻尼矩阵；S为电磁场系数矩阵；J为磁场载荷矩阵。
采用上文中自适应时间步长算法，可得tn时刻电磁场最佳离散时间步长为。
在只考虑热传导和对流条件下，温度场导热微分方程可写成[19]式中，为密度，kg/m3；为比热容，J/(kg·K)；为发热功率，W。
初始及边界条件为式中，为初始温度分布；Ta为环境温度；为边界温度分布；为恒定温度边界条件；为散热量边界条件；q为边界发热功率；h为边界对流换热系数。
采用伽辽金法将式（21）和式（22）形成有限元格式，即（23）式中，C为温度场比热矩阵；K为温度场导热矩阵；P为能量载荷矩阵。
同理，采用第2节中自适应时间步长算法，可得tn时刻温度场最适合离散步长为DtnT。
然而实际计算温度场时，还需考虑DtnT与耦合时间间隔Dt（Dt=tn+1-tn）的大小。
可采用如下方式计算实际温度场离散步长DtT：1）当时，在tn~tn+1时间段内温度场离散步长DtT=Dt。
2）当时，在tn~tn+1时间段内温度场离散步长Dt为（24）式中，[ ]为取整符号，取大于或等于某数的最小整数。
本节以铜导环电磁-热分析为例说明自适应时间步长在电磁-温度场耦合中的应用。
如图5a所示，将一长度为10 mm、厚度为2 mm的铜导环置于空气中，铜导环的上下及内外表面对流换热系数均为5 W/(m2·K)。
对铜导环施加电流i=104sin(100πt)（单位：A），持续时间0.1 s，分析铜环的温度变化。
