一种基于自适应时间步长的瞬变电磁‑温度场耦合计算方法与流程

本发明一种基于自适应时间步长的瞬变电磁-温度场耦合计算方法，涉及电磁温度计算领域。

背景技术：

在线圈发射装置、继电器、电机等设备中，运行时产生过高的温升将影响线圈和接头导电性、材料的绝缘性能，甚至会产生破坏性的热膨胀，对于设备运行性能及安全性产生一定的影响。在上述问题中，瞬变的大电流使得导体的温度在较短的时间内迅速增加，因而在计算过程不仅考虑电磁场对温度场的影响，还需考虑温度变化对电磁场材料导电性能的影响，为此需要建立瞬变电磁-温度场顺序强耦合模型。目前采用有限元法求解瞬变电磁-温度场间接耦合时，电磁场和温度场往往采用相同的时间步长。然而在求解过程中，温度场响应较电磁场慢，若采用相同的时间步长计算，使得温度场求解过频，造成了求解时间的增加。

对于两物理场瞬态间接耦合时，物理场响应时间不同的问题。有方法根据现有软件，提出代码耦合概念，对各物理场采用不同的时间离散策略，并在时间的耦合点上进行载荷的传递。也有方法针对在流-固耦合时，对流体区域和固体区域采用不同的时间步长的DFMT-SCSS算法，选取固体求解区域时间步长为流体求解区域的倍数计算。然而上述方法在计算两物理场耦合时，两物理场时间步长只是存在简单的倍数关系，而物理场都采用恒定步长计算，存在两物理场未获得最佳离散策略的情况。为此提出自适应步长耦合概念，即各物理场采用自适应步长算法获取最佳离散时间策略，并在预测的耦合时间节点上自动耦合。可以使各物理场获得最佳的时间离散策略，也避免了响应较慢的物理场过频的计算。

技术实现要素：

为解决上述技术问题，针对瞬态电磁-温度场间接耦合问题，本发明提出一种基于自适应时间步长的瞬变电磁-温度场耦合计算方法。在T_n时刻，如图1所示，首先通过指数平滑法预测下一耦合时间点T_n+1。然后在两耦合时间节点间，通过预测-校正法和响应特征值计算电磁场和温度场最佳离散步长。得到T_n至T_n+1时间段内电磁场和温度场最佳时间离散步长为△t_n^E和△t_n^T。电磁场、温度场分别采用△t_n^E、△t_n^T计算至耦合时间节点。使各物理场在两耦合时间节点上，以最佳的时间步长计算，避免了响应较慢的物理场过频的计算。

本发明所采用的技术方案是：

一种基于自适应时间步长的瞬变电磁-温度场耦合计算方法，包括以下步骤：

1)、初始化。首先建立分析对象的电磁场和温度计算有限元模型，并给定电磁场计算中的磁导率、电阻率和温度场计算的比热容、导热率等材料参数。最后施加初始条件、边界及载荷；

2)、耦合时间点确定。采用温度场触发热量来判断电磁温度场耦合时间节点。首先电磁-温度场采用较小步长耦合三次，用来获取温度场触发热量Q_pre的预测数据。然后采用指数平滑法预测温度场触发热量，通过触发热量判断耦合时间节点；

3)、电磁场自适应步长计算。通过载荷离散误差和响应特征值确定电磁场载荷最佳离散的步长△t_n^E。启动电磁场，采用最佳离散步长计算。当电磁场累积热量达到温度场触发热量时，暂停电磁场计算；

4)、温度场自适应步长计算。同3)中，通过载荷离散误差和响应特征值确定电磁场载荷最佳离散的步长△t_n^T。将电磁场计算平均热功率作为载荷加至温度场，当温度场计算时间与电磁场同步时，暂停温度场计算。更新电磁场节点温度，并进行下一时间步的温度场触发热量计算。反复迭代计算至最终时间。

本发明一种基于自适应时间步长的瞬变电磁-温度场耦合计算方法，优点在于：

1)、由于顺序耦合时，在计算电磁场的单个时间步内，并未考虑温度场的变化，因而存在累积的非平衡误差。采用温度场触发热量来判断耦合时间节点，可以有效控制非平衡误差。

2)、与传统等步长耦合方法相比，避免了由于温度场响应时间比电磁场长，而采用统一计算时间步长导致温度场求解次数过多而增加计算时间的问题，减小了计算时间。

附图说明

图1是电磁-温度场自适应耦合示意图。

图2是电磁-温度场自适应步长耦合流程图。

图3是载荷离散误差图。

图4(a)是铜导环结构示意图；

图4(b)是有限元模型图。

图5(a)等步长温度计算温度分布图；

图5(b)自适应步长温度分布图。

图6(a)是1号单元计算结果图；

图6(b)是140号单元计算结果图。

具体实施方式

一种基于自适应时间步长的瞬变电磁-温度场耦合计算方法，包括以下步骤：

1)、初始化。首先建立分析对象的电磁场和温度计算有限元模型，并给定电磁场计算中的磁导率、电阻率和温度场计算的比热容、导热率等材料参数。最后施加相应的初始条件、边界及载荷；

2)、耦合时间点确定。采用温度场触发热量来判断电磁温度场耦合时间节点。首先电磁-温度场采用较小步长耦合三次，用来获取温度场触发热量Q_pre的预测数据。然后采用指数平滑法预测温度场触发热量，通过触发热量判断耦合时间节点；

3)、电磁场自适应步长计算。通过载荷离散误差和响应特征值确定电磁场载荷最佳离散的步长△t_n^E。启动电磁场，采用最佳离散步长计算。当电磁场累积热量达到温度场触发热量时，暂停电磁场计算；

4)、温度场自适应步长计算。同3)中，通过载荷离散误差和响应特征值确定电磁场载荷最佳离散的步长△t_n^T。将电磁场计算平均热功率作为载荷加至温度场，当温度场计算时间与电磁场同步时，暂停温度场计算。更新电磁场节点温度，并进行下一时间步的温度场触发热量计算。反复迭代计算至最终时间。

具体步骤如下：

步骤1)：建立电磁场计算有限元模型，并给定磁导率、电阻率材料参数，施加载荷。由似稳条件忽略位移电流，时变电磁场方程可写为如下形式的矢量磁位方程：

式中：Ω₁为导体区域，Ω₂为非导体区域。A为矢量磁位(Wb/m)；V为电位(V)；σ为电导率(S/m)；μ为磁导率(H/m)；J_s为源电流密度(A/m²)。

采用伽辽金法将上式写成有限元格式：

式中：R为电磁场阻尼矩阵；S为电磁场系数矩阵；J为磁场载荷矩阵。

步骤2)：建立温度场场计算有限元模型，并给定比热容、导热率材料参数。在只考虑热传导和对流条件下，温度场导热微分方程可写成如下形式：

式中：ρ为密度(Kg/m)，C_p为比热容(J/(Kg.K))，Q为发热功率W。

初始条件及边界条件为：

式中：G(x,y,z)表示初始温度分布；F(x,y,z)表示恒定温度边界条件；Γ_q表示散热量边界条件，q为边界发热功率，h为边界对流换热系数。

采用伽辽金法将式(3)，(4)形成有限元格式如下：

式中：C为温度场比热矩阵，K为温度场导热矩阵，P为载荷矩阵。

步骤3)：采用温度场触发热量来判断电磁温度场耦合时间节点。电磁-温度场采用较小步长耦合三次，用来获取温度场触发热量Q_pre的预测数据。

步骤4)：获取电磁场计算过程中允许材料变化最大温度。线性材料热功率P与电流密度为J关系：

P_n＝∫_VJ_n²ε(T)dV (6)

ε＝aT+ε₀ (7)

式中：V为单元体积，P_n为t_n时刻发热功率。J_n为T_n时刻电流密度，ε为电阻率，a为电阻随温度变化率，ε₀为0℃时电阻率。

当输入热功率为P_n，由温度变化引起功率计算误差应满足：

|P_n(T_n)-P_n(T_n+△T)|≤γP_n(T_n) (8)

式(8)取等号时，可得出T_n时刻温度允许变化最大值△T_max。

步骤5)：温度场触发热量计算。根据T_n时刻最大温度变化△T_max。由T_n，T_n-1，T_n-2时刻输入功率及温度变化，预测T_n+1温度场触发热量

首先计算各时刻输入热量及温度变化，如表1。

表1各时刻热量、温度

再计算各时刻温度随热量变化率，如表2。

表2时刻温度变化率

最后计算温度场触发热量计算。采用指数平滑方法预测t_n+1时刻温度场变化率K_n+1：

K_n+1＝αK_n+(1-α)K_n-1+(1-α)²K_n-2 (9)

温度场触发热量为：

Q^pre＝△T_maxK_n+1 (10)

步骤6)：启动电磁场，根据载荷离散误差确定电磁场时间步长t∈(t_n-1,t_n)时，当载荷矩阵P采用线性插值时，等效载荷采用斜坡加载，如图3所示。由载荷离散产生的误差为如式(13)。

对式(13)采用梯形公式积分，得离散误差近似如式(14)。

根据式(14)，载荷离散误差近似正比于(△t)²，可将下一步长计算可分为如下两步：

(a)步长预测。根据第n步计算误差，预测第n+1步长△t¹_n+1。

式中：为安全系数，e_tolerance为允许最大误差；为第n时间步内产生的离散误差。

(b)步长校正。判断当第n+1时间步长△t¹_n+1所产生误差是否满足e^k+1<e_tolerance。如不满足采用(8)进行修正迭代计算，直至满足e^k<e_tolerance。

式中：为安全系数，

步骤7)：根据响应特征值确定电磁场时间步长Hughes提出根据响应特征λ_r值确定计算稳定时间步长△t_n+1的方法^[16]，定义△t_n+1λ为震荡限制条件。当△t_n+1λ>>1时系统处于震荡状态。为保证计算稳定性可取最大步长需满足：

式中：f<1，f为稳定系数；λ_r为响应特征值；△u_n为t_n-1到t_n时间段场量u的变化。

步骤8)：电磁场在t_n时刻离散时间步长为：

步骤9)：耦合时间判断。当T_n+1时刻电磁场累积热量满足以下条件之一时，T_n+1为电磁-温度场耦合时间点，启动温度场计算。

(a)单元最大热量变化达到步骤5)中单元预测热量阈值时自动启动温度场计算：

(b)总体热量变化达到步骤5)中单元预测热量阈值时自动启动温度场计算：

步骤10)：启动温度场计算，计算温度场自适应时间步长温度场自适应时间步长计算同步骤6)、步骤7)、步骤8)；

步骤11)：当温度场计算时间到达步骤9)中耦合时间T_n+1时，停止温度场计算，更新电磁场节点温度。

步骤12)：反复迭代步骤4)～步骤10)，直至计算总时间T_total。

实施例：

以通电铜导环瞬态温升为例:

本节以电铜导环热分析为例说明自适应时间步长在电磁-热耦合中的应用，该模型广泛存在于继电器、线圈发射等装置中。将一长10mm，厚2mm铜导环置于空气中，如图4(a)所示)，铜导环的上下及内外表面对流换热系数均为5W/m²。对铜导环施加电流i＝10⁴sin(50πt)，持续时间0.1s，分析铜环的温度变化。

建立轴对称模型，如图5(b)所示，电磁场计算区域包含铜环、空气区域，采用三角划分网格，共988个节点，2061个单元。温度场计算区域为铜环区域，也采用三角形划分网格，共121个节点，200个单元。

当式(15)中电磁场载荷离散误差取为e_tolerance＝0.5％，式(8)、式(21)中取γ＝1％，β＝1。电磁-温度场前三个时间步内采用1ms耦合，随后采用自适应时间步长耦合。最后将自适应时间步长计算结果与定步长计算结果对比，如图5(a)、图5(b)所示，分别为铜导环施加交流载荷0.1s时，采用等步长和自适应时间步长计算所得的温度分布云图。

由图6(a)、图6(b)可见，采用自适应时间步长所得的铜环温度分布与定步长温度分布基本一致。其中1号单元温度最高而相对误差最小为0.35％，140号单元温度最低而相对误差最大为0.53％。为此选取1号单元和140号单元，分析其在整个加热过程中温度及误差的变化规律。

1号和140号在0.1s内自适应时间步长与定步长计算所得单元温度对比如图6(a)、图6(b)所示，铜导块温度近似以0.02s为周期上升。以0.06s～0.08s为例，在区域a和c中，电流幅值较小，温度场升较慢，对应的温度场计算步长较大。而区域b，电流幅值较大，温度上升较快，对应温度场计算时间步长较小，说明电磁-温度场自动步长耦合算法的有效性。在整个计算过程中，自动步长计算与定步长计算误差在在0.7％以内，明了该方法的准确性。

电磁-温度场自适应时间步长与定步长耦合计算性能对比如表3所示。在0.1s内，采用定步长耦合，选取时间步长为1ms时，电磁场和温度场均需计算100次，整体计算时间为235s。而采用自适应步长求解时，电磁场计算93次，温度场计算29次，计算时间为184s，计算时间减小21％。

表3计算性能对比