基于传输线迭代的2D轴对称非线性静磁场模型的求解方法与流程
本发明涉及数值计算领域，具体而言，涉及一种基于传输线迭代法的非线性静磁场模型的有限元求解方法，该方法主要针对2D轴对称非线性静态电磁场进行求解。
背景技术：
有限元法是工业设计中最常用的数值计算方法，被诸多商用仿真软件采用，应用广泛。然而，随着求解模型的日益复杂化以及分网单元数目的不断增多，以传统的牛顿迭代法为核心的非线性有限元求解方法面临着求解耗时严重的问题，这直接关系到产品研发的速度和效率。
有限元问题的求解的核心在于求解线性方程组，而对于非线性问题来说，传统的牛顿迭代法每一步都要利用新的迭代结果重新生成有限元模型的全局矩阵，随着模型分网的不断增大，全局矩阵的维度不断变大，每一步矩阵的LU分解等消耗的时间会相应的增大，总体的求解时间可能随着分网的变密而成几何式增大。
因此，需要研究一种新的迭代方法，以解决牛顿迭代法求解有限元非线性问题时带来的求解时间长，效率低的问题。
技术实现要素：
本发明提供了一种基于传输线迭代法的2D轴对称非线性静磁场模型的有限元求解方法，用以解决牛顿迭代法求解有限元非线性问题时带来求解时间长，效率低的问题。
为了达到上述目的，本发明提供了一种基于传输线迭代法的2D轴对称非线性静磁场模型的有限元求解方法，其包括以下步骤：
S1：确定待求解的变量以及求解域，待求解的变量为一2D轴对称非线性静磁场的磁势，2D轴对称非线性静磁场由通电线圈中的电流产生，通电线圈周围的各元件均为铁磁材料，求解域为2D轴对称非线性静磁场所在的区域；
S2：建立一平面x-y坐标系，以2D轴对称非线性静磁场的对称轴为y轴，在y轴上选定其中一点为原点，并设定经过原点并与y轴垂直的直线为x轴，即x-y平面为2D轴对称非线性静磁场所在区域一过对称轴的截面所在的平面；
S3：列出2D轴对称非线性静磁场中的控制方程和边界条件式并组成一微分方程组，其控制方程为：
其中，J为电流密度变量，μ为三角单元的磁导率，A为磁势，
边界条件式为：
Γ1：A＝0，
Γ1表示磁势A在边界Γ1上的分布，Γ2表示磁势A沿边界的外法线方向的变化率，
S4：使用三角单元对求解域进行分网，得到包含多个三角单元的有限元网络，该有限元网络中的三角单元总个数为N，节点总个数为M，并分别对三角单元和节点进行1～N和1～M的编号，其中1000≤N≤3000；
S5：根据微分方程组的泛函形式，推导出每一个三角单元的单元矩阵[Ye]和激励源单元矩阵[Je]，其中，每一[Ye]均为3×3的矩阵，每一[Je]均为1×3的矩阵：
[Je]＝[Jl Jm Jn]，
l、m、n分别为推导每一三角单元的[Ye]和[Je]时，三角单元的三个顶点的编号，
r和s分别为三角单元的三个顶点编号1、m和n中的其中两个顶点编号，
x1、xm和xn分别为节点l、节点m和节点n在平面坐标系中的横坐标，y1、ym和yn分别为节点l、节点m和节点n在平面坐标系中的纵坐标，Δ为节点l、节点m和节点n组成的三角单元的面积；
S6：根据得到的每一个三角单元的单元矩阵[Ye]和激励源单元矩阵[Je]，对N个三角单元进行有限元装配，得到全局矩阵Y和J，其中Y为M×M矩阵，J为M×1矩阵；
S7：求解非线性方程组YA＝J，得到2D轴对称非线性静磁场中每个节点的磁势A，其中A为M×1的节点磁势矩阵，A＝[A1 A2 … AM]T；
S8：根据步骤S8中计算得到的节点磁势矩阵A，按照以下各式计算每一个三角单元的磁感应强度B，其中，
S9：根据铁磁材料的B-H曲线以及步骤S8中计算得到的每一个三角单元的磁感应强度B，并计算出每一个三角单元的磁导率μ；
S10：以步骤S4中的分网结果为基础，对求解域进行精细的三角分网，得到三角单元总个数为N'、节点总个数为M'的有限元网络，并分别对三角单元和节点进行1～N'和1～M'的编号；
S11：按照步骤S5中的方法，对步骤S10中得到的有限元网络再次计算每个三角单元的单元矩阵[Ye]和激励源单元矩阵[Je]；
S12：有限元网络转化为电路模型，将步骤S11中得到的单元矩阵[Ye]视为电路的导纳矩阵，激励源单元矩阵[Je]视为每个节点与地之间的电流源矩阵，对有限元网络中的每一个三角单元均建立一个等效电路网络，建立等效电路网络的方法如下：
将单元矩阵[Ye]对角线上的元素视为自导，非对角线上的元素视为互导，
对于非对角线上的元素，若Yrs＞0，则在三角单元对应的等效电路网络中的节点r和节点s之间设置一受控电流源，该受控电流源中的电流大小为UrsYrs,方向为从节点r流向节点s，其中Urs为节点r与节点s之间的磁势差，
对于非对角线上的元素，若Yrs＜0，则在三角单元对应的等效电路网络中的节点r和节点s之间设置一纯电阻，该纯电阻的导纳为|Yrs|，
若有限元单元矩阵[Ye]的第r行的所有元素之和不等于0，当第r行所有元素之和大于零时，在节点r与地之间设置一纯电阻，该纯电阻的导纳为Yrl+Yrm+Yrn，当第r行所有元素之和小于零时，则在节点r与地之间设置一受控电流源，该受控电流源中的电流大小为Ur0·|Yrl+Yrm+Yrn|,方向为从节点r流向地，其中Ur0为节点r与地之间的磁势差，
在每一节点与地之间均设置一电流源，节点l、节点m、节点n与地之间的电流源中的电流大小分别为Jl、Jm、Jn，电流方向为由地流向节点；
S13：组装电路，将步骤S12中建立的每一个三角单元对应的等效电路网络通过节点进行连接，组装成一个完整的非线性电路网络，该非线性电路网络等效为包含一线性网络与多个非线性待求元件的电路；
S14：对于步骤S13中得到的非线性电路网络，为了使用传输线迭代方法进行求解，需要在非线性元件与线性网络之间添加一段传输线，由于传输线对信号传输的延时作用，电路的非线性求解过程包括入射阶段和反射阶段，
入射阶段，非线性电路元件的电压信号向线性网络进行入射，等效为传输线导纳与虚拟电流源的并联，
反射阶段，电压信号由线性网络传向非线性元件，对非线性元件进行求解，如此不断迭代入射阶段和反射阶段，直到电路达到稳态，
(一)在线性部分与非线性元件之间添加一段传输线，传输线的导纳的计算方法如下：
(1)确定每一个三角单元的磁导率μ的估计值，检查经过步骤S10分网后得到的三角单元的重心对应的第一次分网的三角单元，并将对应的第一次分网的三角单元的磁导率设为三角单元的磁导率，
(2)非线性元件的导纳是一个关于磁导率μ的变量，将上一步得到的μ值代入到非线性元件表达式，得到的结果作为对应的传输线的导纳值，
(二)设迭代开始时每一节点的电压均为0，当第n个节点电压信号以Vin反射到线性网络时，每一非线性待求元件等效为一导纳和一电流源的并联电路，其中，导纳为对应的传输线导纳Yn，电流源中的电流值为2VinYn，对该等效电路进行求解，得到每一节点的电压值Vin，
(三)根据各个节点的电压值，利用非线性元件与电压之间的关系式，计算并更新非线性元件的导纳值，
(四)计算各个节点向非线性元件入射的电压值Vrn，节点n处的Vrn＝Vn-Vin，
(五)入射过程，每一非线性待求元件等效为一导纳与一电流源的并联电路，其中，导纳为对应的传输线导纳Yn，电流源中的电流值为2VrnYn，得到每一非线性元件两端的电压
(六)计算各个节点向线性网络入射的电压值Vin，节点n处的
(七)重复步骤(二)～(六)，直至相邻两次迭代中，步骤(二)所求的电压值Vn达到预设的收敛误差，此时计算得到的各节点的电压值Vn即为所求电压值，
S15：根据每一节点的电压值绘制2D轴对称非线性静磁场中的磁势云图。
在本发明的一实施例中，步骤S14中(二)对等效电路的求解为使用节点电压法进行求解，其步骤为：
(一)计算得到矩阵YV＝I，其中Y为电路导纳矩阵，由于迭代过程中，导纳矩阵Y保持不变，仅需计算一次，V为待求节点电压，I为节点电流；
(二)迭代第一步对矩阵Y进行LU分解，即Y＝LU，其中L为单位下三角矩阵，U为上三角矩阵，之后的每一次迭代，无需计算此步，直接计算步骤(三)；
(三)使用公式V＝U-1(L-1I)求解节点电压V。
在本发明的一实施例中，步骤S14中(五)，入射过程中，每一非线性元件两端电压的求解是独立的，在此使用多核并行技术同时对多个非线性元件两端电压进行求解。
本发明提供的基于传输线迭代法的2D轴对称非线性静磁场模型的有限元求解方法解决了牛顿迭代法求解有限元非线性问题时带来求解时间长，效率低的问题。本发明中的每一步骤均无需再次计算全局矩阵，只需要进行一次全局矩阵的LU分解，然后即可重复使用，从而节省计算时间；同时，该方法非常适合使用并行算法进行加速，能够进一步的加快有限元问题的求解。相对于传统的牛顿迭代法，本发明在求解时间上有着很大的优势，有着广阔的应用前景。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为产生一2D轴对称非线性静磁场的接触器机械结构示意图；
图2为与图1对应的静磁场的有限元模型区域示意图；
图3a为对求解域进行首次分网(粗糙分网)的示意图；
图3b为对求解域进行再次分网(精细分网)的示意图；
图4为求解域中的三角单元的示意图；
图5为等效电路网络的示意图；
图6为对三角单元进行组装的示意图；
图7为传输线迭代等效示意图；
图8为反射过程的等效示意图；
图9为入射过程的等效示意图；
图10为静磁场中的磁势云图；
图11为传统牛顿迭代法与传输线迭代法求解时间对比；
图12为牛顿迭代法与传输线迭代法不同分网大小计算时间对比；
图13为牛顿迭代法与传输线迭代法单步计算时间对比。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明提供的基于传输线迭代法的2D轴对称非线性静磁场模型的有限元求解方法包括以下步骤：
S1：确定待求解的变量以及求解域，待求解的变量为一2D轴对称非线性静磁场的磁势，2D轴对称非线性静磁场由通电线圈中的电流产生，通电线圈周围的各元件均为铁磁材料，求解域为2D轴对称非线性静磁场所在的区域；
图1为产生一2D轴对称非线性静磁场的接触器机械结构示意图，如图所示，图1中为一接触器，其中的线圈中通有电流，线圈周围的铁芯、推杆、衔铁等的材质均为铁磁材料，该接触器周围的静磁场以推杆为中心而轴对称，因此，可以先以静磁场的其中任意一个截面(如图中的方框部分)为研究对象。
图2为与图1对应的静磁场的有限元模型区域示意图，该有限元模型区域为图1中的俯视区域的右侧部分，该区域也就是本发明针对的求解域。
S2：建立一平面x-y坐标系，以2D轴对称非线性静磁场的对称轴为y轴，在y轴上选定其中一点为原点，并设定经过原点并与y轴垂直的直线为x轴，即x-y平面为2D轴对称非线性静磁场所在区域一过对称轴的截面所在的平面；
S3：列出2D轴对称非线性静磁场中的控制方程和边界条件式并组成一微分方程组，其控制方程为：
其中，J为电流密度变量，μ为三角单元的磁导率，A为磁势，
边界条件式为：
Γ1：A＝0，
Γ1表示磁势A在边界Γ1上的分布，Γ2表示磁势A沿边界的外法线方向的变化率，
S4：使用三角单元对求解域进行分网，得到包含多个三角单元的有限元网络，该有限元网络中的三角单元总个数为N，节点总个数为M，并分别对三角单元和节点进行1～N和1～M的编号，其中1000≤N≤3000；
图3为对求解域进行首次分网(粗糙分网)的示意图。
S5：根据微分方程组的泛函形式，推导出每一个三角单元的单元矩阵[Ye]和激励源单元矩阵[Je]，其中，每一[Ye]均为3×3的矩阵，每一[Je]均为1×3的矩阵：
[Je]＝[Jl Jm Jn]，
图4为求解域中的三角单元的示意图，l、m、n分别为推导每一三角单元的[Ye]和[Je]时，三角单元的三个顶点的编号，
r和s分别为三角单元的三个顶点编号1、m和n中的其中两个顶点编号，
x1、xm和xn分别为节点l、节点m和节点n在平面坐标系中的横坐标，y1、ym和yn分别为节点l、节点m和节点n在平面坐标系中的纵坐标，Δ为节点l、节点m和节点n组成的三角单元的面积；
S6：根据得到的每一个三角单元的单元矩阵[Ye]和激励源单元矩阵[Je]，对N个三角单元进行有限元装配，得到全局矩阵Y和J，其中Y为M×M矩阵，J为M×1矩阵；
S7：求解非线性方程组YA＝J，得到2D轴对称非线性静磁场中每个节点的磁势A，其中A为M×1的节点磁势矩阵，A＝[A1 A2…AM]T；
S8：根据步骤S8中计算得到的节点磁势矩阵A，按照以下各式计算每一个三角单元的磁感应强度B，其中，
S9：根据铁磁材料的B-H曲线以及步骤S8中计算得到的每一个三角单元的磁感应强度B，并计算出每一个三角单元的磁导率μ；
S10：以步骤S4中的分网结果为基础，对求解域进行精细的三角分网，图3b为对求解域进行再次分网(精细分网)的示意图，得到三角单元总个数为N'、节点总个数为M'的有限元网络，并分别对三角单元和节点进行1～N'和1～M'的编号；
S11：按照步骤S5中的方法，对步骤S10中得到的有限元网络再次计算每个三角单元的单元矩阵[Ye]和激励源单元矩阵[Je]；
S12：有限元网络转化为电路模型，将步骤S11中得到的单元矩阵[Ye]视为电路的导纳矩阵，激励源单元矩阵[Je]视为每个节点与地之间的电流源矩阵，对有限元网络中的每一个三角单元均建立一个等效电路网络，图5为等效电路网络的示意图，如图所示，建立等效电路网络的方法如下：
将单元矩阵[Ye]对角线上的元素视为自导，非对角线上的元素视为互导，
对于非对角线上的元素，若Yrs＞0，则在三角单元对应的等效电路网络中的节点r和节点s之间设置一受控电流源，该受控电流源中的电流大小为UrsYrs,方向为从节点r流向节点s，其中Urs为节点r与节点s之间的磁势差，
对于非对角线上的元素，若Yrs＜0，则在三角单元对应的等效电路网络中的节点r和节点s之间设置一纯电阻，该纯电阻的导纳为|Yrs|，
若有限元单元矩阵[Ye]的第r行的所有元素之和不等于0，当第r行所有元素之和大于零时，在节点r与地之间设置一纯电阻，该纯电阻的导纳为Yrl+Yrm+Yrn，当第r行所有元素之和小于零时，则在节点r与地之间设置一受控电流源，该受控电流源中的电流大小为Ur0·|Yrl+Yrm+Yrn|,方向为从节点r流向地，其中Ur0为节点r与地之间的磁势差，
在每一节点与地之间均设置一电流源，节点l、节点m、节点n与地之间的电流源中的电流大小分别为Jl、Jm、Jn，电流方向为由地流向节点；
S13：组装电路，图6为对三角单元进行组装的示意图，如图所示，将步骤S12中建立的每一个三角单元对应的等效电路网络通过节点进行连接，组装成一个完整的非线性电路网络，该非线性电路网络等效为包含一线性网络与多个非线性待求元件的电路；
S14：对于步骤S13中得到的非线性电路网络，为了使用传输线迭代方法进行求解，需要在非线性元件与线性网络之间添加一段传输线，由于传输线对信号传输的延时作用，电路的非线性求解过程包括入射阶段和反射阶段，图7为传输线迭代等效示意图，
入射阶段，非线性电路元件的电压信号向线性网络进行入射，等效为传输线导纳与虚拟电流源的并联，
反射阶段，电压信号由线性网络传向非线性元件，对非线性元件进行求解，如此不断迭代入射阶段和反射阶段，直到电路达到稳态，
(一)在线性部分与非线性元件之间添加一段传输线，传输线的导纳的计算方法如下：
(1)确定每一个三角单元的磁导率μ的估计值，检查经过步骤S10分网后得到的三角单元的重心对应的第一次分网的三角单元，并将对应的第一次分网的三角单元的磁导率设为三角单元的磁导率，
(2)非线性元件的导纳是一个关于磁导率μ的变量，将上一步得到的μ值代入到非线性元件表达式，得到的结果作为对应的传输线的导纳值，
(二)设迭代开始时每一节点的电压均为0，当第n个节点电压信号以Vin反射到线性网络时，每一非线性待求元件等效为一导纳和一电流源的并联电路，其中，导纳为对应的传输线导纳Yn，电流源中的电流值为2VinYn，对该等效电路进行求解，得到每一节点的电压值Vin，图8为反射过程的等效示意图，
(三)根据各个节点的电压值，利用非线性元件与电压之间的关系式，计算并更新非线性元件的导纳值，
(四)计算各个节点向非线性元件入射的电压值Vrn，节点n处的Vrn＝Vn-Vin，
(五)入射过程，每一非线性待求元件等效为一导纳与一电流源的并联电路，其中，导纳为对应的传输线导纳Yn，电流源中的电流值为2VrnYn，得到每一非线性元件两端的电压图9为入射过程的等效示意图，
(六)计算各个节点向线性网络入射的电压值Vin，节点n处的
(七)重复步骤(二)～(六)，直至相邻两次迭代中，步骤(二)所求的电压值Vn达到预设的收敛误差，此时计算得到的各节点的电压值Vn即为所求电压值，
S15：根据每一节点的电压值绘制2D轴对称非线性静磁场中的磁势云图，如图10所示为静磁场中的磁势云图。
在本发明中，步骤S14中(二)对等效电路的求解可以使用节点电压法进行求解，其步骤为：
(一)计算得到矩阵YV＝I，其中Y为电路导纳矩阵，由于迭代过程中，导纳矩阵Y保持不变，仅需计算一次，V为待求节点电压，I为节点电流；
(二)迭代第一步对矩阵Y进行LU分解，即Y＝LU，其中L为单位下三角矩阵，U为上三角矩阵，之后的每一次迭代，无需计算此步，直接计算步骤(三)；
(三)使用公式V＝U-1(L-1I)求解节点电压V。
在本发明中，步骤S14中(五)，入射过程中，每一非线性元件两端电压的求解是独立的，在此使用多核并行技术同时对多个非线性元件两端电压进行求解。
下面阐述本发明的有益技术效果：
相比传统牛顿迭代法求解非线性有限元问题，使用本发明提供的传输线迭代法，能够非常显著的减少计算所用的时间。图11～图13对比了传统牛顿迭代法与传输线迭代法的计算效率。图11中，在单核计算下，传统的牛顿迭代法计算时间几乎是传输线迭代法的6倍，使用本发明提供的方法，求解速度得到了很大的提升；图12中，随着求解模型的分网单元增加，牛顿迭代法与传输线迭代法的求解时间比值不断增大，说明出本发明能够有效的处理有限元模型变复杂的情况；图13显示的是这两种方法的单步求解时间对比，可见，本发明具有非常显著的时间优势，能大幅度提高计算效率。
本发明提供的基于传输线迭代法的2D轴对称非线性静磁场模型的有限元求解方法解决了牛顿迭代法求解有限元非线性问题时带来求解时间长，效率低的问题。本发明中的每一步骤均无需再次计算全局矩阵，只需要进行一次全局矩阵的LU分解，然后即可重复使用，从而节省计算时间；同时，该方法非常适合使用并行算法进行加速，能够进一步的加快有限元问题的求解。相对于传统的牛顿迭代法，本发明在求解时间上有着很大的优势，有着广阔的应用前景。
本领域普通技术人员可以理解：附图只是一个实施例的示意图，附图中的模块或流程并不一定是实施本发明所必须的。
本领域普通技术人员可以理解：实施例中的装置中的模块可以按照实施例描述分布于实施例的装置中，也可以进行相应变化位于不同于本实施例的一个或多个装置中。上述实施例的模块可以合并为一个模块，也可以进一步拆分成多个子模块。
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围。
- 还没有人留言评论。精彩留言会获得点赞！