CN106649939B - 基于传输线迭代的2d轴对称非线性静磁场模型的求解方法 - Google Patents

基于传输线迭代的2d轴对称非线性静磁场模型的求解方法 Download PDF

Info

Publication number
CN106649939B
CN106649939B CN201610859242.5A CN201610859242A CN106649939B CN 106649939 B CN106649939 B CN 106649939B CN 201610859242 A CN201610859242 A CN 201610859242A CN 106649939 B CN106649939 B CN 106649939B
Authority
CN
China
Prior art keywords
node
matrix
nonlinear
triangular unit
admittance
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201610859242.5A
Other languages
English (en)
Other versions
CN106649939A (zh
Inventor
杨文英
彭飞
郭久威
贾楠
翟国富
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Institute of Technology
Original Assignee
Harbin Institute of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Institute of Technology filed Critical Harbin Institute of Technology
Priority to CN201610859242.5A priority Critical patent/CN106649939B/zh
Publication of CN106649939A publication Critical patent/CN106649939A/zh
Application granted granted Critical
Publication of CN106649939B publication Critical patent/CN106649939B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]

Abstract

本发明提供了一种基于传输线迭代法的2D轴对称非线性静磁场模型的有限元求解方法，该方法对求解域进行粗糙分网和精细分网两次分网，并在非线性元件与线性网络之间添加相应的传输线，以通过多次迭代入射阶段和反射阶段对电路进行求解，进而得到2D轴对称非线性静磁场中的磁势云图。与现有的牛顿迭代法相比，本发明在求解时间上有着很大的优势，有着广阔的应用前景。

Description

基于传输线迭代的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分解，然后即可重复使用，从而节省计算时间；同时，该方法非常适合使用并行算法进行加速，能够进一步的加快有限元问题的求解。相对于传统的牛顿迭代法，本发明在求解时间上有着很大的优势，有着广阔的应用前景。
本领域普通技术人员可以理解：附图只是一个实施例的示意图，附图中的模块或流程并不一定是实施本发明所必须的。
本领域普通技术人员可以理解：实施例中的装置中的模块可以按照实施例描述分布于实施例的装置中，也可以进行相应变化位于不同于本实施例的一个或多个装置中。上述实施例的模块可以合并为一个模块，也可以进一步拆分成多个子模块。
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围。

Claims (3)

1.一种基于传输线迭代法的2D轴对称非线性静磁场模型的有限元求解方法，其特征在于，包括以下步骤：
S1：确定待求解的变量以及求解域，待求解的变量为一2D轴对称非线性静磁场的磁势，2D轴对称非线性静磁场由通电线圈中的电流产生，通电线圈周围的各元件均为铁磁材料，求解域为2D轴对称非线性静磁场所在的区域；
S2：建立一平面x-y坐标系，以2D轴对称非线性静磁场的对称轴为y轴，在y轴上选定其中一点为原点，并设定经过原点并与y轴垂直的直线为x轴，即x-y平面为2D轴对称非线性静磁场所在区域一过对称轴的截面所在的平面；
S3：列出2D轴对称非线性静磁场中的控制方程和边界条件式并组成一微分方程组，其控制方程为：
其中，J为电流密度变量，μ为三角单元的磁导率，A为磁势，
边界条件式为：
Γ1：A＝0，
Γ2：
Γ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分网后得到的三角单元的重心对应的第一次分网的三角单元，并将对应的第一次分网的三角单元的磁导率设为经过步骤S10分网后得到的三角单元的磁导率，
(2)非线性元件的导纳是一个关于磁导率μ的变量，将上一步得到的μ值代入到非线性元件表达式，得到的结果作为对应的传输线的导纳值，
(二)设迭代开始时每一节点的电压均为0，当第n个节点电压信号以Vin反射到线性网络时，每一非线性待求元件等效为一导纳和一电流源的并联电路，其中，导纳为对应的传输线导纳Yn，电流源中的电流值为2VinYn，对该等效电路进行求解，得到每一节点的电压值Vin，
(三)根据各个节点的电压值，利用非线性元件与电压之间的关系式，计算并更新非线性元件的导纳值，
(四)计算各个节点向非线性元件入射的电压值Vrn，节点n处的Vrn＝Vn-Vin，
(五)入射过程，每一非线性待求元件等效为一导纳与一电流源的并联电路，其中，导纳为对应的传输线导纳Yn，电流源中的电流值为2VrnYn，得到每一非线性元件两端的电压
(六)计算各个节点向线性网络入射的电压值Vin，节点n处的
(七)重复步骤(二)～(六)，直至相邻两次迭代中，步骤(二)所求的电压值Vn达到预设的收敛误差，此时计算得到的各节点的电压值Vn即为所求电压值，
S15：根据每一节点的电压值绘制2D轴对称非线性静磁场中的磁势云图。
2.根据权利要求1所述的基于传输线迭代法的2D轴对称非线性静磁场模型的有限元求解方法，其特征在于，步骤S14中(二)对等效电路的求解为使用节点电压法进行求解，其步骤为：
步骤(1)：计算得到矩阵YV＝I，其中Y为电路导纳矩阵，由于迭代过程中，导纳矩阵Y保持不变，仅需计算一次，V为待求节点电压，I为节点电流；
步骤(2)：迭代第一步对矩阵Y进行LU分解，即Y＝LU，其中L为单位下三角矩阵，U为上三角矩阵，之后的每一次迭代，无需计算此步，直接计算步骤(3)；
步骤(3)：使用公式V＝U-1(L-1I)求解节点电压V。
3.根据权利要求书1所述的基于传输线迭代法的2D轴对称非线性静磁场模型的有限元求解方法，其特征在于，步骤S14中(五)，入射过程中，每一非线性元件两端电压的求解是独立的，在此使用多核并行技术同时对多个非线性元件两端电压进行求解。
CN201610859242.5A 2016-09-28 2016-09-28 基于传输线迭代的2d轴对称非线性静磁场模型的求解方法 Active CN106649939B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610859242.5A CN106649939B (zh) 2016-09-28 2016-09-28 基于传输线迭代的2d轴对称非线性静磁场模型的求解方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610859242.5A CN106649939B (zh) 2016-09-28 2016-09-28 基于传输线迭代的2d轴对称非线性静磁场模型的求解方法

Publications (2)

Publication Number Publication Date
CN106649939A CN106649939A (zh) 2017-05-10
CN106649939B true CN106649939B (zh) 2019-11-01

Family

ID=58853474

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610859242.5A Active CN106649939B (zh) 2016-09-28 2016-09-28 基于传输线迭代的2d轴对称非线性静磁场模型的求解方法

Country Status (1)

Country Link
CN (1) CN106649939B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107609274B (zh) * 2017-09-14 2020-03-13 哈尔滨工业大学 基于传输线与级别调度法的二维静磁场并行有限元方法
CN109408927B (zh) * 2018-10-13 2022-01-18 哈尔滨工业大学 一种基于黑盒传输线模型的二维静磁场并行有限元加速方法
CN113326642B (zh) * 2021-05-08 2024-01-26 英特工程仿真技术（大连）有限公司 一种含薄气隙结构的轴对称电磁场气隙力计算方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102054070A (zh) * 2009-10-30 2011-05-11 新思科技(上海)有限公司 非线性电路直流工作点的支路电流计算方法与装置
CN104298809A (zh) * 2014-08-27 2015-01-21 天津大学 一种基于矩阵指数电磁暂态仿真的非线性建模求解方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102054070A (zh) * 2009-10-30 2011-05-11 新思科技(上海)有限公司 非线性电路直流工作点的支路电流计算方法与装置
CN104298809A (zh) * 2014-08-27 2015-01-21 天津大学 一种基于矩阵指数电磁暂态仿真的非线性建模求解方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
新型等效源法求解轴对称非线性静磁场;闫照文 等;《华中理工大学学报》;20001130;第28卷(第11期);57-59 *
轴对称静磁场的矢量边界元法;郭航 等;《电工技术学报》;19950531;第1995年卷(第2期);33-37 *
轴对称非线性磁场分布的有限元分析;谢干权 等;《计算数学》;19810228;第1981年卷(第1期);44-56 *

Also Published As

Publication number Publication date
CN106649939A (zh) 2017-05-10

Similar Documents

Publication Publication Date Title
Zheng et al. New strategies for some issues of numerical manifold method in simulation of crack propagation
CN106649939B (zh) 基于传输线迭代的2d轴对称非线性静磁场模型的求解方法
Li et al. A comparison of 3D particle, fluid and hybrid simulations for negative streamers
Faganello et al. Magnetic reconnection and Kelvin–Helmholtz instabilities at the Earth's magnetopause
Wiengarten et al. Implementing turbulence transport in the CRONOS framework and application to the propagation of CMEs
Bu et al. A note on block representations of the group inverse of Laplacian matrices
Palmroth et al. Preliminary testing of global hybrid-Vlasov simulation: Magnetosheath and cusps under northward interplanetary magnetic field
Ryang Three-spin giant magnons in AdS5× S5
Huang et al. Kinetic simulations of secondary reconnection in the reconnection jet
Matyjasek Extremal limit of the regular charged black holes in nonlinear electrodynamics
Bufferand et al. Implementation of multi-component Zhdanov closure in SOLEDGE3X
Guo et al. Re‐reconnection processes of magnetopause flux ropes: Three‐dimensional global hybrid simulations
Goulart et al. Disformal invariance of Maxwell’s field equations
Pasko et al. Mesosphere‐troposphere coupling due to sprites
Ngobeni et al. Cosmic ray anisotropies in the outer heliosphere
Fu et al. A three‐dimensional MHD simulation of the multiple X line reconnection process
Hesse et al. On the topology of flux transfer events
He et al. A magnetic null geometry reconstructed from Cluster spacecraft observations
Elder et al. Magnetic nulls in interacting dipolar fields
Cherkis et al. The't Hooft–Polyakov monopole in the presence of a't Hooft operator
Chernicoff et al. Thin-shell wormholes in AdS5 and string dioptrics
Zalesak et al. Finite temperature effects on the evolution of ionospheric barium clouds in the presence of a conducting background ionosphere: 1. The simplest case: Incompressible background ionosphere, equipotential magnetic field lines, and an altitude‐invariant neutral wind
Rapoport et al. Excitation of planetary electromagnetic waves in the inhomogeneous ionosphere
Aveiro et al. Equatorial spread F studies using SAMI3 with two-dimensional and three-dimensional electrostatics
Xie et al. Darwin model in plasma physics revisited

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant