基于传输线与级别调度法的二维静磁场并行有限元方法

摘要

基于传输线与级别调度法的二维静磁场并行有限元方法，属于电器数值计算领域，该方法主要针对二维非线性静态电磁场进行求解，包括二维平面和二维轴对称情形。本发明的优点是：采用传输线迭代法和级别调度法进行有限元的迭代求解，在迭代求解过程中，全局矩阵Y能够保持不变，在矩阵求解过程当中，采用LU分解法，只需要在计算的第一步进行LU分解，由于LU分解一般占用矩阵求解的95％左右的时间，使用这种方法，在每一个迭代步当中不需要再次执行全局矩阵的LU分解过程，能够节约95％的时间。同时，我们将级别调度法运用到了LU分解之后的矩阵三角求解过程当中，该算法能有效的加速三角求解过程。

说明书

技术领域

本发明属于电器数值计算领域，具体涉及一种传输线迭代法与级别调度法结合的二维静磁场有限元并行加速求解方法，该方法主要针对二维非线性静态电磁场进行求解，包括二维平面和二维轴对称情形。

背景技术

有限元法是工业设计中最常用的数值计算方法，被诸多商用仿真软件采用，应用广泛。然而，随着求解模型的日益复杂化以及分网单元数目的不断增多，以传统的牛顿迭代法为核心的非线性有限元求解方法面临着求解耗时严重的问题，这直接关系到产品研发的速度和效率。

有限元问题的求解的核心在于求解线性方程组，而对于非线性问题来说，传统的牛顿迭代法每一步都要利用新的迭代结果重新生成有限元模型的全局矩阵，随着模型分网的不断增大，全局矩阵的维度不断变大，每一步矩阵的LU分解等消耗的时间会相应的增大，总体的求解时间可能随着分网的变密而成几何式增大。

因此，需要研究一种新的迭代方法，以解决牛顿迭代法求解有限元非线性问题时带来的求解时间长，效率低的问题。

发明内容

本发明的目的是为了解决现有牛顿迭代法求解有限元非线性问题时带来求解时间长、效率低的问题，提供一种将传输线迭代法和级别调度法结合来实现对二维非线性静磁场模型进行有限元并行加速求解的方法，该方法能够对所求问题并行地进行求解计算。

为实现上述目的，本发明采取的技术方案如下：

基于传输线与级别调度法的二维静磁场并行有限元加速方法，所述方法具体步骤如下：

步骤一：建立一个二维平面坐标系，建立所求静磁场问题的几何模型；

步骤二：对于二维非线性静磁场当中的控制方程和边界条件，得到一组微分方程组，其控制方程为：

其中，A为待求变量磁势，μ₀为空气磁导率，M为磁化强度矢量，α_m为M与x轴正方向之间的夹角，J为电流密度；边界条件为：

Γ₁：A＝0，Γ₁表示磁势A在边界上的分布；

Γ₂表示磁势A沿边界的外法线方向的变化率；

对于二维轴对称结构的静磁场,得到一组微分方程组，其控制方程为：

令μ'＝xμ，A'＝xA，得到

其中，A'为待求变量磁势，μ'为空气磁导率，M为磁化强度矢量，α_m为M与x轴正方向之间的夹角，J为电流密度；边界条件为：

Γ₁：A'＝0，Γ₁表示磁势A在边界上的分布；

Γ₂表示磁势A沿边界的外法线方向的变化率；

二维和二维轴对称问题，均可采用以下步骤进行求解，

步骤三：进行求解区域的粗糙分网，使用三角单元对求解域进行离散分网，得到包含多个三角形单元的有限元网络，该有限元网络中的三角单元总个数为N1，节点总个数为M1，并分别对三角单元和节点进行1～N1和1～M1的编号；

步骤四：在每一个三角单元当中，计算出单元系数矩阵[Y^e]，激励电流源单元矩阵[J^e]，磁化强度矩阵[M^e]，其中，[Y^e]为3×3的矩阵，[J^e]和[M^e]均为3×1的矩阵，设K、L、N为三角形单元三个节点按照逆时针的编号；其中，二维情况下，矩阵[Y^e]中的每一个元素为(i，j＝K，L，N)，i，j为变量的下标，b，c为变量表达式，具体形式为

Δ为三角形面积，μ为材料的磁导率，x，y为点的横、纵坐标值；

二维轴对称情况下，矩阵[Y^e]中的每一个元素为(i，j＝K，L，N)，i，j为变量的下标，b，c为变量表达式，具体形式为

Δ为三角形面积，μ为材料的磁导率，x，y为点的横、纵坐标值；

步骤五：根据得到的每一个三角单元的单元系数矩阵[Y^e]和激励电流源单元矩阵[J^e]以及磁化强度矩阵[M^e]，对N1个三角单元进行有限元装配，得到全局矩阵Y、J和M，其中Y为M1×M1矩阵，J、M为M1×1矩阵；

步骤六：根据步骤五得到的全局矩阵Y、J、M，求解非线性方程组YA＝J+M，得到二维轴对称非线性静磁场中每个节点的磁势A，其中A为M1×1的节点磁势矩阵，A＝[A₁>2>M1]^T；

步骤七：根据步骤六中计算得到的节点磁势矩阵A，按照以下各式计算每一个三角单元的磁感应强度B，其中，

其中，B_x为磁感应强度B在x方向上的分量，B_y为磁感应强度B在y方向上的分量，A_K、A_L、A_N分别为ΔKLN三点的磁势值；

步骤八：根据铁磁材料的B-H曲线以及步骤七中计算得到的每一个三角单元的磁感应强度B，并计算出每一个三角单元的磁导率μ；

步骤九：进行精细分网，以步骤三中的分网结果为基础，对求解域进行精细的三角分网，得到三角单元总个数为N2、节点总个数为M2的有限元网络，并分别对三角单元和节点进行1～N2和1～M2的编号；

步骤十：按照步骤四中的方法，对步骤九中得到的有限元网络再次计算每个三角单元的单元系数矩阵[Y^e]₁和激励电流源单元矩阵[J^e]₁、磁化强度矩阵[M^e]₁；

步骤十一：有限元网络转化为电路模型，将步骤十中得到的单元系数矩阵[Y^e]₁视为电路的导纳矩阵，激励电流源单元矩阵[J^e]₁和磁化强度矩阵[M^e]₁视为每个节点与地之间的电流源矩阵，对有限元网络中的每一个三角单元均建立一个等效电路网络；

步骤十二：组装电路，将步骤十一中建立的每一个三角单元对应的等效电路网络通过节点进行连接，组装成一个完整的非线性电路网络，该非线性电路网络等效为包含一线性网络与多个非线性待求元件的电路；

步骤十三：对于步骤十二中得到的非线性电路网络，为了使用传输线迭代方法进行求解，需要在非线性元件与线性网络之间添加一段传输线，由于传输线对信号传输的延时作用，电路的非线性求解过程包括入射阶段和反射阶段，

入射阶段，非线性电路元件的电压信号向线性网络进行入射，等效为传输线导纳与虚拟电流源的并联，

反射阶段，电压信号由线性网络传向非线性元件，对非线性元件进行求解，如此不断迭代入射阶段和反射阶段，直到电路达到稳态，

(一)在线性部分与非线性元件之间添加一段传输线，传输线的导纳的计算方法如下：

(1)确定每一个三角单元的磁导率μ的估计值，检查经过步骤九分网后得到的三角单元的重心对应的第一次分网的三角单元，并将对应的第一次分网的三角单元的磁导率设为三角单元的磁导率，

(2)非线性元件的导纳是一个关于磁导率μ的变量，将上一步得到的μ值代入到非线性元件表达式，得到的结果作为对应的传输线的导纳值，

(二)设迭代开始时每一节点的电压均为0，当第n个节点电压信号以入射电压V_in反射到线性网络时，每一非线性待求元件等效为一导纳和一电流源的并联电路，其中，导纳为对应的传输线导纳Y_n，电流源中的电流值为2V_inY_n，对该等效电路进行求解，求解表达式为YV＝J+M+I_c+2V_inY_n，其中I_c为电路当中的压控电流源，Y为电路的导纳矩阵，迭代过程当中Y保持不变，当得到每一节点的电压值V_in，YV＝J+M+I_c+2V_inY_n的求解使用LU分解法，由于矩阵Y在迭代过程当中保持不变，只需要在迭代的第一步将Y分解为LU，即Y＝LU，LU分别为下三角矩阵和上三角矩阵，之后每次求解都直接进行三角求解LUV＝J+M+I_c+2V_inY_n，求解采用并行的级别调度法进行；

(三)根据各个节点的电压值，利用非线性元件与电压之间的关系式，计算并更新非线性元件的导纳值，

(四)计算各个节点向非线性元件入射的电压值V_rn，节点n处的V_rn＝V_n-V_in，

(五)入射过程，每一非线性待求元件等效为一导纳与一电流源的并联电路，其中，导纳为对应的传输线导纳Y_n，电流源中的电流值为2V_rnY_n，得到每一非线性元件两端的电压V_NEn，由于各个元件的计算能够单独计算，这一步采用并行算法进行计算，

(六)计算各个节点向线性网络入射的电压值V_in，节点n处的V_in＝V_NEn-V_rn，

(七)重复步骤(二)～(六)，直至相邻两次迭代中，步骤(二)所求的电压值V_n达到预设的收敛误差，此时计算得到的各节点的电压值V_n即为所求电压值；

步骤十四：根据每一节点的电压值绘制二维非线性静磁场中的磁势云图。

本发明相对于现有技术的有益效果是：采用传输线迭代法和级别调度法进行有限元的迭代求解，在迭代求解过程中，全局矩阵Y能够保持不变，在矩阵求解过程当中，采用LU分解法，只需要在计算的第一步进行LU分解，由于LU分解一般占用矩阵求解的95％左右的时间，使用这种方法，在每一个迭代步当中不需要再次执行全局矩阵的LU分解过程，能够节约95％的时间。同时，我们将级别调度法运用到了LU分解之后的矩阵三角求解过程当中，该算法能有效的加速三角求解过程。

附图说明

图1为下三角矩阵以及级别图；

图2为下三角矩阵按照级别排序图；

图3为一个大小为4096×4096的下三角矩阵按照级别大小进行排序之后的矩阵；

图4为静磁场几何模型图；

图5为利用牛顿迭代法计算CPU核心数与分网单元数目下每一个迭代步的计算时间曲线图；

图6为利用本发明方法计算CPU核心数与分网单元数目下每一个迭代步的计算时间曲线图；

图7为利用牛顿迭代法计算CPU核心数与分网单元数目下总计算时间曲线图；

图8为利用本发明方法计算CPU核心数与分网单元数目下总计算时间曲线图。

具体实施方式

下面结合附图对本发明的技术方案作进一步的说明，但并不局限于此，凡是对本发明技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围，均应涵盖在本发明的保护范围中。

具体实施方式一：本实施方式记载的是基于传输线与级别调度法的二维静磁场并行有限元加速方法，所述方法具体步骤如下：

步骤一：建立一个二维平面坐标系，建立所求静磁场问题的几何模型，如图4所示；

步骤二：对于二维非线性静磁场当中的控制方程和边界条件，得到一组微分方程组，其控制方程为：

其中，A为待求变量磁势，μ₀为空气磁导率，M为磁化强度矢量，α_m为M与x轴正方向之间的夹角，J为电流密度；边界条件为：

Γ₁：A＝0，Γ₁表示磁势A在边界上的分布；

Γ₂表示磁势A沿边界的外法线方向的变化率；

对于二维轴对称结构的静磁场,得到一组微分方程组，其控制方程为：

令μ'＝xμ，A'＝xA，得到

其中，A'为待求变量磁势，μ'为空气磁导率，M为磁化强度矢量，α_m为M与x轴正方向之间的夹角，J为电流密度；边界条件为：

Γ₁：A'＝0，Γ₁表示磁势A在边界上的分布；

Γ₂表示磁势A沿边界的外法线方向的变化率；

二维和二维轴对称问题，均可采用以下步骤进行求解，

步骤三：进行求解区域的粗糙分网，使用三角单元对求解域进行离散分网，得到包含多个三角形单元的有限元网络，该有限元网络中的三角单元总个数为N1，节点总个数为M1，并分别对三角单元和节点进行1～N1和1～M1的编号；

步骤四：在每一个三角单元当中，计算出单元系数矩阵[Y^e]，激励电流源单元矩阵[J^e]，磁化强度矩阵[M^e]，其中，[Y^e]为3×3的矩阵，[J^e]和[M^e]均为3×1的矩阵，设K、L、N为三角形单元三个节点按照逆时针的编号；其中，二维情况下，矩阵[Y^e]中的每一个元素为(i，j＝K，L，N)，i，j为变量的下标，b，c为变量表达式，具体形式为

Δ为三角形面积，μ为材料的磁导率，x，y为点的横、纵坐标值；

二维轴对称情况下，矩阵[Y^e]中的每一个元素为(i，j＝K，L，N)，i，j为变量的下标，b，c为变量表达式，具体形式为

Δ为三角形面积，μ为材料的磁导率，x，y为点的横、纵坐标值；

步骤五：根据得到的每一个三角单元的单元系数矩阵[Y^e]和激励电流源单元矩阵[J^e]以及磁化强度矩阵[M^e]，对N1个三角单元进行有限元装配，得到全局矩阵Y、J和M，其中Y为M1×M1矩阵，J、M为M1×1矩阵；

步骤六：根据步骤五得到的全局矩阵Y、J、M，求解非线性方程组YA＝J+M，得到二维轴对称非线性静磁场中每个节点的磁势A，其中A为M1×1的节点磁势矩阵，A＝[A₁>2>M1]^T；

步骤七：根据步骤六中计算得到的节点磁势矩阵A，按照以下各式计算每一个三角单元的磁感应强度B，其中，

其中，B_x为磁感应强度B在x方向上的分量，B_y为磁感应强度B在y方向上的分量，A_K、A_L、A_N分别为ΔKLN三点的磁势值；

步骤八：根据铁磁材料的B-H曲线以及步骤七中计算得到的每一个三角单元的磁感应强度B，并计算出每一个三角单元的磁导率μ；

步骤九：进行精细分网，以步骤三中的分网结果为基础，对求解域进行精细的三角分网，得到三角单元总个数为N2、节点总个数为M2的有限元网络，并分别对三角单元和节点进行1～N2和1～M2的编号；

步骤十：按照步骤四中的方法，对步骤九中得到的有限元网络再次计算每个三角单元的单元系数矩阵[Y^e]₁和激励电流源单元矩阵[J^e]₁、磁化强度矩阵[M^e]₁；

步骤十一：有限元网络转化为电路模型，将步骤十中得到的单元系数矩阵[Y^e]₁视为电路的导纳矩阵，激励电流源单元矩阵[J^e]₁和磁化强度矩阵[M^e]₁视为每个节点与地之间的电流源矩阵，对有限元网络中的每一个三角单元均建立一个等效电路网络；

步骤十二：组装电路，将步骤十一中建立的每一个三角单元对应的等效电路网络通过节点进行连接，组装成一个完整的非线性电路网络，该非线性电路网络等效为包含一线性网络与多个非线性待求元件的电路；

步骤十三：对于步骤十二中得到的非线性电路网络，为了使用传输线迭代方法进行求解，需要在非线性元件与线性网络之间添加一段传输线，由于传输线对信号传输的延时作用，电路的非线性求解过程包括入射阶段和反射阶段，

入射阶段，非线性电路元件的电压信号向线性网络进行入射，等效为传输线导纳与虚拟电流源的并联，

反射阶段，电压信号由线性网络传向非线性元件，对非线性元件进行求解，如此不断迭代入射阶段和反射阶段，直到电路达到稳态，

(一)在线性部分与非线性元件之间添加一段传输线，传输线的导纳的计算方法如下：

(1)确定每一个三角单元的磁导率μ的估计值，检查经过步骤九分网后得到的三角单元的重心对应的第一次分网的三角单元，并将对应的第一次分网的三角单元的磁导率设为三角单元的磁导率，

(2)非线性元件的导纳是一个关于磁导率μ的变量，将上一步得到的μ值代入到非线性元件表达式，得到的结果作为对应的传输线的导纳值，

(二)设迭代开始时每一节点的电压均为0，当第n个节点电压信号以入射电压V_in反射到线性网络时，每一非线性待求元件等效为一导纳和一电流源的并联电路，其中，导纳为对应的传输线导纳Y_n，电流源中的电流值为2V_inY_n，对该等效电路进行求解，求解表达式为YV＝J+M+I_c+2V_in>n，其中I_c为电路当中的压控电流源，Y为电路的导纳矩阵，迭代过程当中Y保持不变，当得到每一节点的电压值V_in，YV＝J+M+I_c+2V_inY_n的求解使用LU分解法，由于矩阵Y在迭代过程当中保持不变，只需要在迭代的第一步将Y分解为LU，即Y＝LU，LU分别为下三角矩阵和上三角矩阵，之后每次求解都直接进行三角求解LUV＝J+M+I_c+2V_inY_n，求解采用并行的级别调度法进行；

(三)根据各个节点的电压值，利用非线性元件与电压之间的关系式，计算并更新非线性元件的导纳值，

(四)计算各个节点向非线性元件入射的电压值V_rn，节点n处的V_rn＝V_n-V_in，

(五)入射过程，每一非线性待求元件等效为一导纳与一电流源的并联电路，其中，导纳为对应的传输线导纳Y_n，电流源中的电流值为2V_rnY_n，得到每一非线性元件两端的电压V_NEn，由于各个元件的计算能够单独计算，这一步采用并行算法进行计算，

(六)计算各个节点向线性网络入射的电压值V_in，节点n处的V_in＝V_NEn-V_rn，

(七)重复步骤(二)～(六)，直至相邻两次迭代中，步骤(二)所求的电压值V_n达到预设的收敛误差，此时计算得到的各节点的电压值V_n即为所求电压值；

步骤十四：根据每一节点的电压值绘制二维非线性静磁场中的磁势云图。

具体实施方式二：具体实施方式一所述的基于传输线与级别调度法的二维静磁场并行有限元方法，步骤十一中，等效电路网络的建立方法如下：

将单元矩阵[Y^e]对角线上的元素视为自导，非对角线上的元素视为互导，

对于非对角线上的元素，若表示矩阵[Y^e]第r行，第s列的元素，则在三角单元对应的等效电路网络中的节点r和节点s之间设置一压控电流源，该受控电流源中的电流大小为U_rsY_rs，方向为从节点r流向节点s，其中U_rs为节点r与节点s之间的磁势差，

对于非对角线上的元素，若则在三角单元对应的等效电路网络中的节点r和节点s之间设置一纯电阻，该纯电阻的导纳为

在每一节点与地之间均设置一电流源，节点K、节点L、节点N与地之间的电流源中的电流大小分别为J_K、J_L、J_N，电流方向为由地流向节点。

具体实施方式三：具体实施方式一所述的基于传输线和级别调度法的二维静磁场并行有限元方法，步骤十三的(二)中，级别调度法的建立方法如下：

以下三角矩阵Lx＝b的求解为例，L为一个大小为n×n的下三角矩阵，x，b均为一个n×1的矩阵，其计算步骤为：

A、计算矩阵L当中第i变量x(i)的级别level(i)，对于第i(i＝1,2,…,n)行中的所有元素,如果矩阵L的第i行，第j(j＝1,2,…,n)列的元素L_ij不为零，那么更新level(i)＝1+maxlevel(j)，直到遍历所有元素，

B、计算得到每一个变量的级别之后，将矩阵L按照变量x的级别大小，对其进行行列变换，使得变量的级别为从小到大的顺序，新排序之后的矩阵变为L'，x和b也跟着发生变换，得到L'x'＝b'，如图1，2分别为某矩阵按照级别大小排序前后，小方框内的数字为该列变量x的级别；

C、如图3，为一个大小为4096×4096的下三角矩阵按照级别大小进行排序之后的矩阵，其中黑色小点代表该位置元素非零，级别不同的变量使用虚线分隔开，在同一级别的变量能够同时进行求解，在级别level 1中，通过x'(j)＝b'(j)/L'(j)计算变量的值，在同一级别内，将各个变量的求解放到不同的计算单元中，实现并行计算，

在level 1中对角线上正方形区域下的长方形区域内的变量值已经在level 1当中进行了求解，在等式L'x'＝b'中，将已求出的x'代入，然后将常数移到右侧，更新b'，

与上面的计算方法一样，接下来计算level 2当中的正方形区域内的变量值，使用公式x'(j)＝b'(j)/L'(j)，求解变量值，在等式L'x'＝b'中，将已求出的x'代入，然后将常数移到右侧，再次更新b'，

之后的级别计算与上述一样，直到计算结束，完成未知变量x的求解；

上三角矩阵Ux＝b的求解方法与L的求解过程类似，只不过是从相反的方向做一遍。

从图5、图6、图7、图8通过时间对比发现，与传统有限元软件当中所使用的牛顿迭代法相比，本发明所用计算时间有很大的优势。