CN109408927B - 一种基于黑盒传输线模型的二维静磁场并行有限元加速方法 - Google Patents

Abstract

本发明公开了一种基于黑盒传输线模型的二维非线性静磁场模型的有限元求解方法，所述方法包括如下步骤：一、确定待求解的变量以及求解域；二、建立rz直角坐标系；三、列出二维轴对称非线性静磁场中的控制方程和边界条件式并组成一微分方程组；四、对求解域进行分网；五、计算有限元的系数矩阵；六、建立等效的黑盒电路模型；七、在黑盒电路模型与线性电路之间***一段传输线线段；八、进行传输线法迭代；九、重复进行迭代，直到迭代结果收敛到固定误差，结束求解。本发明能够对所求问题并行地进行求解计算，并且能够用于复杂的有限元分网模型求解当中，从而解决牛顿迭代法求解有限元非线性问题时带来的求解时间长、效率低的问题。

Description

一种基于黑盒传输线模型的二维静磁场并行有限元加速方法

技术领域

本发明属于电磁场数值计算技术领域，涉及一种二维非线性静磁场模型的有限元求解方法，具体涉及一种将二维非线性静磁场的有限元问题转化为等价的黑盒电路问题，并借助传输线迭代法实现对二维非线性静磁场的有限元并行加速求解的方法。

背景技术

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

有限元问题的求解的核心在于求解线性方程组，而对于非线性问题来说，传统的牛顿迭代法每一步都要利用新的迭代结果重新生成有限元模型的全局矩阵，随着模型分网的不断增大，全局矩阵的维度不断变大，每一步矩阵的LU分解等消耗的时间会相应的增大，总体的求解时间可能随着分网的变密而成几何式增大。使用传输线迭代法能够避免重复的LU分解过程，求解时间。然而，传统的传输线迭代法只能用于一阶三角形单元的有限元模型当中，无法用于复杂的分网单元模型。因此，需要研究一种新的迭代方法，来提高传输线法的适用范围和应对复杂模型的能力，从而解决牛顿迭代法求解有限元非线性问题时带来的求解时间长、效率低的问题。

发明内容

本发明的目的是提供一种基于黑盒传输线模型的二维非线性静磁场模型的有限元求解方法，该方法能够对所求问题并行地进行求解计算，并且能够用于复杂的有限元分网模型求解当中，包括三角形单元和四边形单元等，从而解决牛顿迭代法求解有限元非线性问题时带来的求解时间长、效率低的问题。

本发明的目的是通过以下技术方案实现的：

一种基于黑盒传输线模型的二维非线性静磁场模型的有限元求解方法，包括如下步骤：

一、确定待求解的变量以及求解域，其中：所述待求解的变量为一二维轴对称非线性静磁场的磁势A，所述求解域为二维轴对称非线性静磁场所在的区域；

二、建立一个rz直角坐标系；

三、列出二维轴对称非线性静磁场中的控制方程和边界条件式并组成一微分方程组，其中：

所述控制方程为：

式中，J为电流密度变量，μ′＝rμ，μ为三角单元的磁导率，A′＝rA，A为磁势，r为横坐标，z为纵坐标；

所述边界条件式为：

Γ₁:A′＝0；

式中，m为求解边界上的法向量；

四、采用分网程序对求解域进行分网，将求解域离散为三角形单元或者四边形单元，在每一个有限元单元中，变量

式中，n为有限元单元的节点数目，N_j为单元当中的形函数，

为对应节点的变量值的大小；

五、按照下式计算有限元的系数矩阵：

其中，每一项的计算表达式如下：

六、将步骤五当中的矩阵K^e看作电路的导纳矩阵，b^e看作电流源向量，建立等效的黑盒电路模型；

七、在黑盒电路模型与线性电路之间的连接导线任意位置***一段传输线线段；

八、进行传输线法迭代：

将电压信号入射到线性网络中，求解该电路，即方程组：(Y_linear+Y_TL)A′＝b+2V_iY_TL，在整个迭代过程当中，矩阵(Y_linear+Y_TL)保持不变，在迭代的第1步执行一次矩阵的LU分解操作，后续步骤无需再次执行，求解完成后，计算反射电压V_r＝A′-V_i；

将电压信号入射到非线性有限元单元中，求解该电路，即非线性方程组：

采用牛顿迭代法进行求解，第k步的迭代公式如下：

其中：

每个非线性单元内的方程组求解都单独放到独立的计算核心中求解，实现并行计算，求解完成后，再次计算反射电压

式中，Y_linear为入射到线性网络内的线性电路的导纳矩阵；Y_TL为入射到线性网络内的传输线电路的导纳矩阵；A′为待求节点电压；V_r为反射到非线性电路的反射电压；V_i为入射到线性电路的入射电压；K^e为非线性单元的系数矩阵；

为反射回非线性单元内的传输线的导纳矩阵；

为反射回非线性单元内的反射电压；J^e为雅克比矩阵；

为第k次迭代时的雅克比矩阵；f为计算当中的中间变量；A^e′为单元内每个节点的磁势值。

九、重复进行步骤八的迭代，直到迭代结果收敛到固定误差，结束求解。

相比于现有技术，本发明具有如下优点：

1、采用有限元方法对静磁场进行离散建模，将有限元中的每个单元的单元系数矩阵等效为一个密封的黑盒电路模型，从而实现将有限元方程组等效为一个非线性的电路网络问题。

2、在非线性单元的黑盒模型与线性单元的网络之间***传输线线段，利用传输线的迭代规律，实现传输线迭代，从而实现静磁场问题的求解。

3、在传输线迭代过程中，每个非线性黑盒模型都是相互独立的，它的求解计算可以采用单独的计算资源(CPU或者GPU核心)进行计算，实现并行加速。

4、采用传输线法和黑盒模型能够进行有限元的迭代求解，能够处理多种复杂的有限元模型，例如三角形单元和四边形单元。

5、在迭代求解过程中，全局矩阵Y能够保持不变，在矩阵求解过程当中，采用LU分解法，只需要在计算的第一步进行LU分解，由于LU分解一般占用矩阵求解的95％左右的时间，使用这种方法将大幅减小每一步的迭代时间。

6、使用传输线法能够将非线性单元与线性求解区域进行隔离，实现并行加速效果。

附图说明

图1为黑盒电路模型；

图2为***传输线后的黑盒电路网络；

图3为电压入射到线性网络后的等效电路图；

图4为电压反射回非线性单元后的等效电路图；

图5为迭代求解流程图；

图6为某型号接触器的结构图；

图7为某型号接触器的剖面图；

图8为接触器简化模型；

图9为接触器计算模型；

图10为分网结果，(a)三角形单元，(b)四边形单元；

图11为采用传输线法求解得到不同状态下的磁场分布情况；

图12为单步计算时间加速效果。

具体实施方式

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

具体实施方式一：本实施方式提供了一种基于黑盒传输线模型的二维非线性静磁场模型的有限元求解方法，所述方法包括如下步骤：

一、确定待求解的变量以及求解域，待求解的变量为一二维轴对称非线性静磁场的磁势A，二维轴对称非线性静磁场由通电线圈中的电流产生，通电线圈周围的各元件均为铁磁材料，求解域为二维轴对称非线性静磁场所在的区域。

二、建立一个rz直角坐标系，r为横轴，z为纵轴。

三、列出二维轴对称非线性静磁场中的控制方程和边界条件式并组成一微分方程组，其控制方程为：

其中，J为电流密度变量，μ′＝rμ，μ为三角单元的磁导率，A′＝rA，A为磁势，r为横坐标，z为纵坐标，边界条件为：Γ₁:A′＝0；

m为求解边界上的法向量。

四、采用分网程序对求解区域进行分网，将求解域离散为三角形单元或者四边形单元。在每一个有限元单元中，变量

其中，n为有限元单元的节点数目，N_j为单元当中的形函数，

为对应节点的变量值的大小。

五、计算有限元的系数矩阵。根据伽辽金法，对控制方程进行加权积分，得到方程：

写成矩阵形式K^eA^e′＝b^e，即为：

其中，每一项的计算表达式如下：

六、将步骤五当中的矩阵K^e看作电路的导纳矩阵，b^e看作电流源向量，建立等效的黑盒电路模型，如图1所示。

七、黑盒模型与线性电路之间是通过导线相连接的，在黑盒模型与线性电路之间的这些导线的任意位置***一段传输线线段，如图2所示。

八、进行传输线法迭代。

电压信号入射到线性网络中，如图3所示，求解该电路，即方程组：(Y_linear+Y_TL)A′＝b+2V_iY_TL，在整个迭代过程当中，矩阵(Y_linear+Y_TL)保持不变，在迭代的第1步执行一次矩阵的LU分解操作，后续步骤无需再次执行，求解完成后，计算反射电压V_r＝A′-V_i。

电压信号入射到非线性有限元单元中，如图4所示，求解该电路，即非线性方程组：

该矩阵采用牛顿迭代法进行求解，第k步的迭代公式如下：

其中：

由于每个非线性单元相互之间都是隔离的，每个单元内的方程组求解都单独放到独立的计算核心中求解，实现并行计算。求解完成后，再次计算反射电压

式中，Y_linear为图3当中线性电路的导纳矩阵；Y_TL为图3虚线框内传输线电路的导纳矩阵；A′为待求节点电压；V_r为反射到非线性电路的反射电压；V_i为入射到线性电路的入射电压；K^e为非线性单元的系数矩阵；

为图4中传输线的导纳矩阵；

为反射回单元内的反射电压；J^e为雅克比矩阵；

为第k次迭代时的雅克比矩阵。

不断地进行迭代，直到迭代结果收敛到固定误差，结束求解。求解流程图如图5。

具体实施方式二：本实施方式提供了一种基于黑盒传输线模型的二维非线性静磁场模型的有限元求解方法，所述方法具体实施步骤如下：

一、确定待求解的变量以及求解域，待求解的变量为一二维轴对称非线性静磁场的磁势A，二维轴对称非线性静磁场由通电线圈中的电流产生，通电线圈周围的各元件均为铁磁材料，求解域为二维轴对称非线性静磁场所在的区域。本实施方式针对某型号轴对称结构接触器进行计算，该模型的外观图如图6所示，剖面图如图7所示。对其进行简化，并建模，简化模型如图8所示，计算模型如图9所示。

二、建立一个rz直角坐标系，r为横轴，z为纵轴。

三、列出二维轴对称非线性静磁场中的控制方程和边界条件式并组成一微分方程组，其控制方程为：

其中，J为电流密度变量，μ′＝rμ，μ为三角单元的磁导率，A′＝rA，A为磁势，r为横坐标，z为纵坐标，边界条件为：Γ₁:A′＝0；

m为求解边界上的法向量。

四、采用分网程序对求解区域进行分网，将求解域离散为三角形单元或者四边形单元，分网结果如图10所示。在每一个有限元单元中，变量

其中，n为有限元单元的节点数目，N_j为单元当中的形函数，

为对应节点的变量值的大小。

对于三角形单元，

Δ为三角形单元的面积，

对于四边形单元，

五、计算有限元的系数矩阵。根据伽辽金法，对控制方程进行加权积分，得到方程：

写成矩阵形式K^eA^e′＝b^e，即为：

其中，每一项的计算表达式如下：

六、将步骤五当中的矩阵K^e看作电路的导纳矩阵，b^e看作电流源向量，建立等效的黑盒电路模型，如图1所示。

七、在黑盒模型与线性电路之间***一段传输线线段，如图2所示。

八、进行传输线法迭代。

电压信号入射到线性网络中，如图3所示，求解该电路，即方程组：(Y_linear+Y_TL)A′＝b+2V_iY_TL，然后，计算反射电压V_r＝A′-V_i。

电压信号入射到非线性有限元单元中，如图4所示，求解该电路，即非线性方程组：

该矩阵采用牛顿迭代法进行求解，第k步的迭代公式如下：

其中：

求解完成后，再次计算反射电压

不断地进行迭代，直到迭代结果收敛到固定误差，结束求解。求解流程图如图5。

采用传输线法求解得到不同状态下的磁场分布情况，结果如图11所示。

表1不同的求解模型

case 1-5为三角形单元分网；case 6-10为四边形单元分网。

表2 N-R迭代法和本发明BB-TLM迭代法的误差对比

No.1-5为表1case 5模型当中任取五个点的误差对比；No.6-10为表1case 10模型任取五个点的误差对比。

对表1当中case10模型分别使用1、4、8、16、20个CPU核心进行加速，对比牛顿迭代法与传输线法的单步计算时间，由图12可见，传输线法的单步计算时间比牛顿迭代法快几倍。

对表1当中case10模型分别使用1、4、8、16、20个CPU核心进行加速，对比牛顿迭代法与传输线法的总计算时间，当CPU核心增大时，传输线法的加速效果比牛顿迭代法更加明显。

表3针对表1中的模型分别采用牛顿迭代法和传输线法求解(CPU核心为20)

由表3可知，当分网模型增大时，传输线法的计算时间比牛顿迭代法快。

Claims (1)

1.一种基于黑盒传输线模型的二维非线性静磁场模型的有限元求解方法，其特征在于所述方法包括如下步骤：

一、确定待求解的变量以及求解域，其中：所述待求解的变量为一二维轴对称非线性静磁场的磁势A，所述求解域为二维轴对称非线性静磁场所在的区域；

二、建立一个rz直角坐标系；

三、列出二维轴对称非线性静磁场中的控制方程和边界条件式并组成一微分方程组，其中：

所述控制方程为：

式中，J为电流密度变量，μ′＝rμ，μ为三角单元的磁导率，A′＝rA，A为磁势，r为横坐标，z为纵坐标；

所述边界条件式为：

Γ₁:A′＝0；

Γ₂:

式中，m为求解边界上的法向量；

四、采用分网程序对求解域进行分网，将求解域离散为三角形单元或者四边形单元，在每一个有限元单元中，变量

式中，n为有限元单元的节点数目，N_j为单元当中的形函数，

为对应节点的变量值的大小；

五、按照下式计算有限元的系数矩阵：

其中，每一项的计算表达式如下：

六、将步骤五当中的矩阵K^e看作电路的导纳矩阵，b^e看作电流源向量，建立等效的黑盒电路模型；

七、在黑盒电路模型与线性电路之间的连接导线任意位置***一段传输线线段；

八、进行传输线法迭代：

将电压信号入射到线性网络中，求解该电路，即方程组：(Y_linear+Y_TL)A′＝b+2V_iY_TL，在整个迭代过程当中，矩阵(Y_linear+Y_TL)保持不变，在迭代的第1步执行一次矩阵的LU分解操作，后续步骤无需再次执行，求解完成后，计算反射电压V_r＝A′-V_i；

将电压信号入射到非线性有限元单元中，求解该电路，即非线性方程组：

采用牛顿迭代法进行求解，第k步的迭代公式如下：

每个非线性单元内的方程组求解都单独放到独立的计算核心中求解，实现并行计算，求解完成后，再次计算入射到线性电路的入射电压

式中，Y_linear为入射到线性网络内的线性电路的导纳矩阵；Y_TL为入射到线性网络内的传输线电路的导纳矩阵；A′为待求节点电压；V_r为反射到非线性电路的反射电压；V_i为入射到线性电路的入射电压；K^e为非线性单元的系数矩阵；

为反射回非线性单元内的传输线的导纳矩阵；

为反射回非线性单元内的反射电压；

为第k次迭代时的雅克比矩阵；

九、重复进行步骤八的迭代，直到迭代结果收敛到固定误差，结束求解。

Images