基于传输线迭代的2D轴对称非线性静磁场模型的求解方法与流程

技术特征：

1.一种基于传输线迭代法的2D轴对称非线性静磁场模型的有限元求解方法，其特征在于，包括以下步骤：

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

S2：建立一平面x-y坐标系，以2D轴对称非线性静磁场的对称轴为y轴，在y轴上选定其中一点为原点，并设定经过原点并与y轴垂直的直线为x轴，即x-y平面为2D轴对称非线性静磁场所在区域一过对称轴的截面所在的平面；

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

其中，J为电流密度变量，μ为三角单元的磁导率，A为磁势，

边界条件式为：

Γ₁：A＝0，

Γ₂：

Γ₁表示磁势A在边界Γ₁上的分布，Γ₂表示磁势A沿边界的外法线方向的变化率，

S4：使用三角单元对求解域进行分网，得到包含多个三角单元的有限元网络，该有限元网络中的三角单元总个数为N，节点总个数为M，并分别对三角单元和节点进行1～N和1～M的编号，其中1000≤N≤3000；

S5：根据微分方程组的泛函形式，推导出每一个三角单元的单元矩阵[Y^e]和激励源单元矩阵[J^e]，其中，每一[Y^e]均为3×3的矩阵，每一[J^e]均为1×3的矩阵：

$\[Y^e\]=\left[\begin{array}{ccc}{Y}_{ll}& Y_{lm}& Y_{ln}\\ Y_{ml}& Y_{mm}& Y_{mn}\\ Y_{nl}& Y_{nm}& Y_{nn}\end{array}\right],$

[J^e]＝[J_l J_m J_n]，

l、m、n分别为推导每一三角单元的[Y^e]和[J^e]时，三角单元的三个顶点的编号，

r和s分别为三角单元的三个顶点编号1、m和n中的其中两个顶点编号，

$\left\{\begin{array}{ccc}{b}_l=y_m-y_n& b_m=y_n-y_l& b_n=y_l-y_m\\ c_l=x_n-x_m& c_m=x_l-x_n& c_n=x_m-x_l\end{array}\right.,$

$r_c=\frac{x_l+x_m+x_n}{3},$

$\mathrm{\Δ}=\frac{1}{2}(b_l{c}_m-b_m{c}_l),$

x₁、x_m和x_n分别为节点l、节点m和节点n在平面坐标系中的横坐标，y₁、y_m和y_n分别为节点l、节点m和节点n在平面坐标系中的纵坐标，Δ为节点l、节点m和节点n组成的三角单元的面积；

S6：根据得到的每一个三角单元的单元矩阵[Y^e]和激励源单元矩阵[J^e]，对N个三角单元进行有限元装配，得到全局矩阵Y和J，其中Y为M×M矩阵，J为M×1矩阵；

S7：求解非线性方程组YA＝J，得到2D轴对称非线性静磁场中每个节点的磁势A，其中A为M×1的节点磁势矩阵，A＝[A₁ A₂ … A_M]^T；

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

$B^2=B_x^2+B_y^2,$

$B_x=-\frac{1}{2\mathrm{\Δ}}(c_l{A}_l+c_m{A}_m+c_n{A}_n),$

$B_y=\frac{A_l+A_m+A_n}{x_l+x_m+x_n}+\frac{1}{2\mathrm{\Δ}}(b_l{A}_l+b_m{A}_m+b_n{A}_n);$

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

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

S11：按照步骤S5中的方法，对步骤S10中得到的有限元网络再次计算每个三角单元的单元矩阵[Y^e]和激励源单元矩阵[J^e]；

S12：有限元网络转化为电路模型，将步骤S11中得到的单元矩阵[Y^e]视为电路的导纳矩阵，激励源单元矩阵[J^e]视为每个节点与地之间的电流源矩阵，对有限元网络中的每一个三角单元均建立一个等效电路网络，建立等效电路网络的方法如下：

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

对于非对角线上的元素，若Y_rs＞0，则在三角单元对应的等效电路网络中的节点r和节点s之间设置一受控电流源，该受控电流源中的电流大小为U_rsY_rs,方向为从节点r流向节点s，其中U_rs为节点r与节点s之间的磁势差，

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

若有限元单元矩阵[Y^e]的第r行的所有元素之和不等于0，当第r行所有元素之和大于零时，在节点r与地之间设置一纯电阻，该纯电阻的导纳为Y_rl+Y_rm+Y_rn，当第r行所有元素之和小于零时，则在节点r与地之间设置一受控电流源，该受控电流源中的电流大小为U_r0·|Y_rl+Y_rm+Y_rn|,方向为从节点r流向地，其中U_r0为节点r与地之间的磁势差，

在每一节点与地之间均设置一电流源，节点l、节点m、节点n与地之间的电流源中的电流大小分别为J_l、J_m、J_n，电流方向为由地流向节点；

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

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

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

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

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

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

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

(二)设迭代开始时每一节点的电压均为0，当第n个节点电压信号以V_in反射到线性网络时，每一非线性待求元件等效为一导纳和一电流源的并联电路，其中，导纳为对应的传输线导纳Y_n，电流源中的电流值为2V_inY_n，对该等效电路进行求解，得到每一节点的电压值V_in，

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

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

(五)入射过程，每一非线性待求元件等效为一导纳与一电流源的并联电路，其中，导纳为对应的传输线导纳Y_n，电流源中的电流值为2V_rnY_n，得到每一非线性元件两端的电压

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

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

S15：根据每一节点的电压值绘制2D轴对称非线性静磁场中的磁势云图。

2.根据权利要求1所述的基于传输线迭代法的2D轴对称非线性静磁场模型的有限元求解方法，其特征在于，步骤S14中(二)对等效电路的求解为使用节点电压法进行求解，其步骤为：

(一)计算得到矩阵YV＝I，其中Y为电路导纳矩阵，由于迭代过程中，导纳矩阵Y保持不变，仅需计算一次，V为待求节点电压，I为节点电流；

(二)迭代第一步对矩阵Y进行LU分解，即Y＝LU，其中L为单位下三角矩阵，U为上三角矩阵，之后的每一次迭代，无需计算此步，直接计算步骤(三)；

(三)使用公式V＝U^-1(L^-1I)求解节点电压V。

3.根据权利要求书1所述的基于传输线迭代法的2D轴对称非线性静磁场模型的有限元求解方法，其特征在于，步骤S14中(五)，入射过程中，每一非线性元件两端电压的求解是独立的，在此使用多核并行技术同时对多个非线性元件两端电压进行求解。