[发明专利]基于传输线与级别调度法的二维静磁场并行有限元方法

摘要

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

主权项

基于传输线与级别调度法的二维静磁场并行有限元加速方法，其特征在于：所述方法具体步骤如下：步骤一：建立一个二维平面坐标系，建立所求静磁场问题的几何模型；步骤二：对于二维非线性静磁场当中的控制方程和边界条件，得到一组微分方程组，其控制方程为：∂∂y[1μ0∂A∂y+Mcosαm]+∂∂x[1μ0∂A∂x-Msinαm]=-J]]>其中，A为待求变量磁势，μ0为空气磁导率，M为磁化强度矢量，αm为M与x轴正方向之间的夹角，J为电流密度；边界条件为：Γ1：A＝0，Γ1表示磁势A在边界上的分布；Γ2：Γ2表示磁势A沿边界的外法线方向的变化率；对于二维轴对称结构的静磁场,得到一组微分方程组，其控制方程为：∂∂y[1xμ∂(xA)∂y+Mcosαm]+∂∂x[1xμ∂(xA)∂x-Mcosαm]=-J,]]>令μ'＝xμ，A'＝xA，得到∂∂y[1μ′∂A′∂y+Mcosαm]+∂∂x[1μ′∂A′∂x-Mcosαm]=-J,]]>其中，A'为待求变量磁势，μ'为空气磁导率，M为磁化强度矢量，αm为M与x轴正方向之间的夹角，J为电流密度；边界条件为：Γ1：A'＝0，Γ1表示磁势A在边界上的分布；Γ2：Γ2表示磁势A沿边界的外法线方向的变化率；二维和二维轴对称问题，均可采用以下步骤进行求解，步骤三：进行求解区域的粗糙分网，使用三角单元对求解域进行离散分网，得到包含多个三角形单元的有限元网络，该有限元网络中的三角单元总个数为N1，节点总个数为M1，并分别对三角单元和节点进行1～N1和1～M1的编号；步骤四：在每一个三角单元当中，计算出单元系数矩阵[Ye]，激励电流源单元矩阵[Je]，磁化强度矩阵[Me]，其中，[Ye]为3×3的矩阵，[Je]和[Me]均为3×1的矩阵，设K、L、N为三角形单元三个节点按照逆时针的编号；其中，二维情况下，矩阵[Ye]中的每一个元素为(i，j＝K，L，N)，i，j为变量的下标，b，c为变量表达式，具体形式为bK=yM-yNbM=yN-yKbN=yK-yMcK=xN-xMcM=xK-xNcN=xM-xK,]]>Δ为三角形面积，μ为材料的磁导率，x，y为点的横、纵坐标值；二维轴对称情况下，矩阵[Ye]中的每一个元素为(i，j＝K，L，N)，i，j为变量的下标，b，c为变量表达式，具体形式为Δ为三角形面积，μ为材料的磁导率，x，y为点的横、纵坐标值；步骤五：根据得到的每一个三角单元的单元系数矩阵[Ye]和激励电流源单元矩阵[Je]以及磁化强度矩阵[Me]，对N1个三角单元进行有限元装配，得到全局矩阵Y、J和M，其中Y为M1×M1矩阵，J、M为M1×1矩阵；步骤六：根据步骤五得到的全局矩阵Y、J、M，求解非线性方程组YA＝J+M，得到二维轴对称非线性静磁场中每个节点的磁势A，其中A为M1×1的节点磁势矩阵，A＝[A1 A2 L AM1]T；步骤七：根据步骤六中计算得到的节点磁势矩阵A，按照以下各式计算每一个三角单元的磁感应强度B，其中，B2=Bx2+By2,]]>Bx=-12Δ(clAl+cmAm+cnAn),]]>By=Al+Am+Anxl+xm+xn+12Δ(blAl+bmAm+bnAn);]]>其中，Bx为磁感应强度B在x方向上的分量，By为磁感应强度B在y方向上的分量，AK、AL、AN分别为ΔKLN三点的磁势值；步骤八：根据铁磁材料的B‑H曲线以及步骤七中计算得到的每一个三角单元的磁感应强度B，并计算出每一个三角单元的磁导率μ；步骤九：进行精细分网，以步骤三中的分网结果为基础，对求解域进行精细的三角分网，得到三角单元总个数为N2、节点总个数为M2的有限元网络，并分别对三角单元和节点进行1～N2和1～M2的编号；步骤十：按照步骤四中的方法，对步骤九中得到的有限元网络再次计算每个三角单元的单元系数矩阵[Ye]1和激励电流源单元矩阵[Je]1、磁化强度矩阵[Me]1；步骤十一：有限元网络转化为电路模型，将步骤十中得到的单元系数矩阵[Ye]1视为电路的导纳矩阵，激励电流源单元矩阵[Je]1和磁化强度矩阵[Me]1视为每个节点与地之间的电流源矩阵，对有限元网络中的每一个三角单元均建立一个等效电路网络；步骤十二：组装电路，将步骤十一中建立的每一个三角单元对应的等效电路网络通过节点进行连接，组装成一个完整的非线性电路网络，该非线性电路网络等效为包含一线性网络与多个非线性待求元件的电路；步骤十三：对于步骤十二中得到的非线性电路网络，为了使用传输线迭代方法进行求解，需要在非线性元件与线性网络之间添加一段传输线，由于传输线对信号传输的延时作用，电路的非线性求解过程包括入射阶段和反射阶段，入射阶段，非线性电路元件的电压信号向线性网络进行入射，等效为传输线导纳与虚拟电流源的并联，反射阶段，电压信号由线性网络传向非线性元件，对非线性元件进行求解，如此不断迭代入射阶段和反射阶段，直到电路达到稳态，(一)在线性部分与非线性元件之间添加一段传输线，传输线的导纳的计算方法如下：(1)确定每一个三角单元的磁导率μ的估计值，检查经过步骤九分网后得到的三角单元的重心对应的第一次分网的三角单元，并将对应的第一次分网的三角单元的磁导率设为三角单元的磁导率，(2)非线性元件的导纳是一个关于磁导率μ的变量，将上一步得到的μ值代入到非线性元件表达式，得到的结果作为对应的传输线的导纳值，(二)设迭代开始时每一节点的电压均为0，当第n个节点电压信号以入射电压Vin反射到线性网络时，每一非线性待求元件等效为一导纳和一电流源的并联电路，其中，导纳为对应的传输线导纳Yn，电流源中的电流值为2VinYn，对该等效电路进行求解，求解表达式为YV＝J+M+Ic+2VinYn，其中Ic为电路当中的压控电流源，Y为电路的导纳矩阵，迭代过程当中Y保持不变，当得到每一节点的电压值Vin，YV＝J+M+Ic+2VinYn的求解使用LU分解法，由于矩阵Y在迭代过程当中保持不变，只需要在迭代的第一步将Y分解为LU，即Y＝LU，LU分别为下三角矩阵和上三角矩阵，之后每次求解都直接进行三角求解LUV＝J+M+Ic+2VinYn，求解采用并行的级别调度法进行；(三)根据各个节点的电压值，利用非线性元件与电压之间的关系式，计算并更新非线性元件的导纳值，(四)计算各个节点向非线性元件入射的电压值Vrn，节点n处的Vrn＝Vn‑Vin，(五)入射过程，每一非线性待求元件等效为一导纳与一电流源的并联电路，其中，导纳为对应的传输线导纳Yn，电流源中的电流值为2VrnYn，得到每一非线性元件两端的电压VNEn，由于各个元件的计算能够单独计算，这一步采用并行算法进行计算，(六)计算各个节点向线性网络入射的电压值Vin，节点n处的Vin＝VNEn‑Vrn，(七)重复步骤(二)～(六)，直至相邻两次迭代中，步骤(二)所求的电压值Vn达到预设的收敛误差，此时计算得到的各节点的电压值Vn即为所求电压值；步骤十四：根据每一节点的电压值绘制二维非线性静磁场中的磁势云图。