基于棱边单元法处理电磁场中移动边界问题的耦合方法

技术领域

本发明涉及电磁场仿真技术领域，具体的说是涉及一种基于棱边单元法处 理电磁场中移动边界问题的耦合方法，该方法可用于解决电气装备中电磁场移 动边界问题，诸如处理电机中的转子旋转问题，处理IC（integrated circuit）装备 中的外加旋转磁场与腔室内磁场的耦合问题，处理电磁阀动态过程中的动铁和 阀体间的相对运动问题等问题。

背景技术

棱边单元法是一种在电磁场数值求解中广泛应用的数值方法，采用棱边单 元法处理电磁场中的移动边界一直是一个难点问题。通过对已有文献和技术资 料的详细查询，总结出目前采用棱边单元法处理此类问题时应用较为广泛的4 种技术手段：MB（moving band，移动带）方法，基于Ladder网格的MPC方法， 气隙单元法，以及基于CG法（或预处理CG法）的拉格朗日乘子法。

第一种方法是MB方法，这种方法在模型的运动区域和静止区域之间建立 一个带状单层网格区域，当运动区域的网格发生刚体移动，为避免带状区域网 格出现较大的剪切变形，需要对带状区域进行重新划分网格Remesh，用以重新 匹配运动和静止区域的网格，参见参考文献：Davat B,Ren Z,Lajoie-Mazenc M. The movement in field modeling[J].Magnetics,IEEE Transactions on,1985,21(6): 2296-2298.；第二种方法是基于Ladder（梯形）网格的MPC法（multipoint constraint method，多点约束法），这种方法首先在模型的运动区域和静止区域的边界处建 立滑移面，该滑移面采用Ladder网格，在运动过程中，为保证滑移面上磁矢势 自由度的切向连续性，采用MPC方法对其进行耦合，参见参考文献：ANSYS Manual13.0；第三种是气隙单元法，这种方法根据旋转电机中气隙的几何特点， 得到定子和转子上未知变量之间的解析关系，将气隙处理为一种特殊的宏单元 （macro-element），并与有限元方法相结合进行求解，获得整个电磁场求解域的 解，参见参考文献：Abdel-Razek A,Coulomb J,Feliachi M,et al.Conception of an  air-gap element for the dynamic analysis of the electromagnetic field in electric  machines[J].Magnetics,IEEE Transactions on,1982,18(2):655-659.；第四种方法 是基于CG法（conjugate gradient，共轭梯度法）（或预处理CG法）的拉格朗日 乘子法，这种方法首先在模型的运动区域和静止区域的边界处建立滑移面，然 后采用拉格朗日乘子法对滑移面上的单元建立约束矩阵，来保证滑移面上磁矢 势自由度的连续性，由于拉格朗日乘子约束矩阵使得生成的总体矩阵具有病态 特性，还需要结合CG法（或预处理CG法）来求解方程，以获得合理的解，参 见参考文献：Golovanov C,Coulomb J L,Marechal Y,et al.3D mesh connection  techniques applied to movement simulation[J].Magnetics,IEEE Transactions on, 1998,34(5):3359-3362.。

但是在上述方法在实际运用中都有一定的限制：采用MB方法需要在每个 时间步都进行网格重构、节点插值，以及重新生成树，因此求解精度和效率受 到较大影响；

气隙单元法使得有限元矩阵带宽变的较大，求解效率变低，此外，该方法 的应用仅限于旋转类的电磁问题；

基于Ladder网格的MPC方法：由于滑移面上磁矢势自由度具有方向性（沿 棱边的切向），为保证MPC方法能对其有效地约束，需要对滑移面处采用Ladder 单元；在运动过程中，当滑移面上的节点不再保持一一对应关系时，采用MPC 方法可以直接耦合Ladder网格节点的自由度，从而保证滑移面上自由度的连续 性；但是基于Ladder网格的MPC方法，对滑移面的网格有限制，只能处理梯 形的网格，因此难以应用到复杂模型中；

基于CG法（或预处理CG法）的拉格朗日乘子法：该方法在滑移面上通过 引入拉格朗日矢量乘子，建立磁矢势和拉格朗日矢量乘子的约束矩阵，从而保 证滑移面上磁矢势自由度的连续性；由于生成的总体矩阵具有病态特性，还需 要采用CG法（或预处理CG法）来迭代求解方程组，以在不唯一的解空间中获 得合理的解；基于CG法（或预处理CG法）的拉格朗日乘子法，在采用CG法 （或预处理CG法）求解方程组时，解的收敛性往往得不到保证，从而导致求解 失败。

发明内容

鉴于已有技术存在的缺陷，本发明的目的是要提供一种基于棱边单元法处 理电磁场中移动边界问题的耦合方法，该方法能够有效克服MB方法的求解效 率和精度问题，基于Ladder网格的MPC方法对网格的限制性问题，气隙单元 法的求解效率以及应用范围的问题，以及基于CG法（或预处理CG法）的格朗 日乘子法的方程组求解的收敛性问题等问题。

为了实现上述目的，本发明的技术方案：

基于棱边单元法处理电磁场中移动边界问题的耦合方法，其特征在于：所 述方法包括：

S1：建立待分析对象的三维模型，对三维模型的运动区域及静止区域剖分 网格，并在运动区域及静止区域的边界处形成滑移面；

S2：建立电磁场求解的基本方程，为保证磁矢势的唯一性，引入拉格朗日 乘子规范方案，即采用拉格朗日乘子法施加规范条件:计算各个区域 的棱边单元矩阵，并将所述棱边单元矩阵集成到集成到有限元总体方程，其中为散度（div）计算符号，A为磁矢势；

S3：当运动区域的网格沿着滑移面移动时，在所述三维模型的滑移面上， 采用拉格朗日乘子法对磁矢势A进行耦合，生成拉格朗日乘子约束矩阵，并将 所述约束矩阵集成到所述有限元总体方程中；

S4：在所述三维模型的滑移面上，对S2中施加规范条件时所引入的拉格朗 日标量乘子自由度，采用MPC方法进行耦合；

S5：对S2、S3及S4步骤所得到的有限元总体方程，施加边界条件，进行 求解，得到待分析对象的电磁场分析结果。

与现有技术相比，本发明的有益效果：根据本发明获得的有限元矩阵性态 良好，计算结果精度可靠。本发明方法对处理运动边界，以及不匹配网格耦合 界面的电磁场问题具有统适性，且实施简便，更合适科学实践和工程应用。

附图说明

图1为本发明所示耦合方法的流程示意图；

图2为本发明所示耦合方法的三维模型的滑移面示意图；

图3为本发明所示耦合方法的实施例模型示意图；

图4为本发明所示耦合方法的实施例的棱边单元示意图；

图5为本发明所示耦合方法的实施例的网络模型示意图；

图6为本发明所示耦合方法的实施例所述的磁感应强度矢量图一；

图7为本发明所示耦合方法的实施例所述的磁感应强度矢量图二。

图中：1、静止区域，2、运动区域，3、滑移面，M、运动方向，N、空隙 距离为0，A、滑移面3的A侧，B、滑移面3的B侧。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图，对 本发明进行进一步详细说明。

本发明基本设计原理：首先建立电磁场求解的基本方程，并在电磁场求解 域引入拉格朗日乘子规范方案以保证解的唯一性；然后在电磁场滑移面上，对 磁矢势自由度采用拉格朗日乘子法来保证切向的连续性；最后对拉格朗日乘规 范方案中所引入的拉格朗日标量乘子采用MPC（多点约束法）方法保证其连续 性。

如图1所示，本方法包括如下步骤：

S1：建立待分析对象的三维模型，分别对三维模型的运动区域及静止区域 剖分网格，并在运动区域及静止区域的边界处形成滑移面。

本步的目的是：在运动区域和静止区域的边界上建立一个滑移面，使得运 动区域网格可以在滑移面上移动。滑移面上的自由度通过以下步骤的耦合方法 进行耦合。

如图2所示，鉴于运动区2和静止区域1是两套网格，因此在滑移面3上， 是有两套网格发生重叠的，在滑移面3的A侧和B侧都有面网格，且这两层面 网格可以不匹配，这让我们在剖分网格的时候更加自由，受到的限制比较少。

S2：建立电磁场求解的基本方程，并在电磁场求解域引入拉格朗日乘子规 范方案，即采用拉格朗日乘子法施加规范条件:以保证解的唯一性，计 算各个区域的棱边单元矩阵，并将所述棱边单元矩阵集成到有限元总体方程中， 其中为散度（div）计算符号，A为磁矢势；

所述的S2包括：

S21、建立各个区域的电磁场求解的基本方程：

我们以静磁场为例进行说明，不导电空间中的静磁场方程：

$\left(\begin{array}{cc}▿\times \left[v\right]▿\times A=J_s& \mathrm{in\Omega }---\left(1\right)\end{array}\right)$

边界条件：

$\left(\begin{array}{cc}A\times n=0& \mathrm{on}\partial \Omega ---\left(2\right)\end{array}\right)$

其中[ν]:磁阻矩阵(磁导率矩阵的逆)、A:磁矢势、J_s:电流密度矢量源项、 旋度(curl)计算符号、Ω:不导电的三维空间领域、n:分析领域的外法向、 分析领域的边界；

S22、为保证电磁场有限元方程解的唯一性，需要引入拉格朗日标量乘子p来 对磁矢势A施加规范条件，所述的规范条件为：

$\left(\begin{array}{cc}▿·A=0& \mathrm{in\Omega }---\left(3\right)\end{array}\right)$

对于S21中的磁矢势A采用棱边单元法进行离散，即自由度是在棱边上， 方向沿着棱边方向，而对于所述拉格朗日标量乘子p采用节点元法进行离散，即 自由度在单元的角节点处：

A=[W]{A_e}      （4）

p={N}{p_e}       （5）

这里，A_e表示单元中的磁矢势自由度，p_e表示单元中的拉格朗日标量乘子自 由度，[W]是磁矢势A的形函数矩阵，{N}是拉格朗日标量乘子p的形函数矩阵。

根据（4）、（5），我们可得到对应于（1）、（2）、（3）的单元矩阵方程，并 集成到有限元总体方程：

$\left(\begin{array}{cc}\left(\begin{array}{cc}{K}^{\mathrm{AA}}& G^{\mathrm{AP}}\\ {\left(G^{\mathrm{AP}}\right)}^T& 0\end{array}\right)\left(\begin{array}{c}{A}_e\\ p_e\end{array}\right)=\left(\begin{array}{c}{J}_e^S\\ 0\end{array}\right)& \mathrm{in\Omega }\end{array}\right)---\left(6\right)$

其中，K^AA表示单元磁阻矩阵，G^AP表示拉格朗日标量乘子p对磁矢势A的单 元约束矩阵，(G^AP)^T表示G^AP的转置矩阵，表示施加在单元上的电流密度；

K^AA和G^AP，及的表达式如下：

单元磁阻矩阵$K^{\mathrm{AA}}={\int }_{\mathrm{Ve}}{(▿\times {\left[W\right]}^T)}^T\left[v\right](▿\times {\left[W\right]}^T)\mathrm{dV};$

单元拉格朗日乘子矩阵

单元电流密度源项矢量

S3：当运动区域的网格沿着滑移面移动时，在所述三维模型的滑移面上， 采用拉格朗日乘子法对磁矢势A进行耦合，生成拉格朗日乘子约束矩阵，并对 所述约束矩阵集成到有限元总体方程中。

为了确保滑移面上磁矢势A的切向连续性，在滑移面上，需要引入拉格朗 日矢量乘子λ_Γ耦合磁矢势A的切向量，包括：

S31：建立滑移面处磁矢势A约束方程的弱形式：为保证磁矢势A在滑移面处 的切向连续性，需要在滑移面处采用拉格朗日乘子法施加连续性条件，即在滑 移面引入拉格朗日矢量乘子λ_Γ来施加滑移面处约束条件，得到滑移面处磁矢势 A的约束方程的弱形式：

${\int }_{{\Gamma }_c}(A^A\Lambda n_{\Gamma }-A^B\Lambda n_{\Gamma })·{\lambda }_{\Gamma }\mathrm{dS}=0---\left(7\right)$

其中，A^A、A^B分别表示在滑移面A、B侧的磁矢势，如图2所示，n_Γ表示滑移 面的外法向，Λn_Γ表示与外法向n_Γ垂直的分量，Γ_c表示滑移面。

为了得到含拉格朗日乘子约束的有限元方程，还需要对式（7）进行变分， 变分表达式为：

${\int }_{{\Gamma }_c}(A^A\Lambda n_{\Gamma }-A^B\Lambda n_{\Gamma })·{\mathrm{\delta \lambda }}_{\Gamma }\mathrm{dS}=0---\left(8\right)$

${\int }_{{\Gamma }_c}(\delta A^A\Lambda n_{\Gamma }-\delta A^B\Lambda n_{\Gamma })·{\lambda }_{\Gamma }\mathrm{dS}=0---\left(9\right)$

S32：建立滑移面处有限元方程的矩阵形式：

对于滑移面处的拉格朗日矢量乘子λ_Γ采用棱边单元法进行离散，即自由度 是在棱边上，方向沿着棱边方向：

λ_Γ=[W_λ]{λ_e}      （10）

（10）式中，[W_λ]是拉格朗日矢量乘子的形函数矩阵，其表达形式和磁矢 势A的形函数矩阵相同；{λ_e}是单元中的拉格朗日矢量乘子自由度，方向沿单元 棱边的方向；

结合（8）、（9）、（10）和（6），获得滑移面处的单元矩阵，并进行集成， 得到有限元总体方程：

$\left(\begin{array}{ccc}{K}^{\mathrm{AA}}& G^{\mathrm{AP}}& K_{\Gamma }^{\mathrm{A\lambda }}\\ {\left(G^{\mathrm{AP}}\right)}^T& 0& 0\\ {\left(K_{\Gamma }^{\mathrm{A\lambda }}\right)}^T& 0& 0\end{array}\right)\left(\begin{array}{c}{A}_e\\ p_e\\ {\lambda }_e\end{array}\right)=\left(\begin{array}{c}{J}_e^S\\ 0\\ 0\end{array}\right)---\left(11\right)$

比较（10）和（4），我们在滑移面处引入了λ_e以及相应的拉格朗日乘子约 束矩阵和

$K_{\Gamma }^{\mathrm{A\lambda }}={\int }_{{\Gamma }_c}\left(\right[W^A]\Delta n_{\Gamma }-[W^B\left]\Delta n_{\Gamma }\right)\left[W_{\lambda }\right]\mathrm{dS}---\left(12\right)$

${\left(K_{\Gamma }^{\mathrm{A\lambda }}\right)}^T={\int }_{{\Gamma }_c}{\left[W_{\lambda }\right]}^T({\left[W^A\right]}^T\Lambda n_{\Gamma }-{\left[W^B\right]}^T\Delta n_{\Gamma })\mathrm{dS}---\left(13\right)$

S4：在三维模型的滑移面上，对S2中所引入的拉格朗日标量乘子p，采用 MPC（多点约束法）方法进行耦合。

对于拉格朗日标量乘子p，在滑移耦合界面处，采用MPC多点约束法耦合， 保证拉格朗日标量乘子的连续性。

S5：对S2、S3及S4步骤所得到的有限元总体方程，施加边界条件，并进 行求解，得到待分析对象的电磁场分析结果。

通过求解有限元总体方程，获得整个求解域的磁矢势A，并通过公式（14）、 （15），求得磁感应强度B及磁场强度H：

$B=▿\times A---\left(14\right)$

H=[ν]B       （15）

其中：旋度(curl)计算符号；[ν]:磁阻矩阵(磁导率矩阵的逆)。

下面以对旋转线圈的磁场分析的具体实施例为例，说明上述步骤：

本案例的模型由定子、转子，以及和转子共同旋转的线圈组成，我们需要 求解：当转子和线圈旋转时，整个区域中磁场的分布，这是一个比较典型的电 磁场移动边界问题。

其具体过程：

a)首先建立待分析对象的三维几何模型，包括旋转线圈、定子部件等（如 图3）；

b)对三维几何模型的运动区域和静止区域进行网格剖分，将所述三维几何 模型剖分为如图所示的六面体棱边单元（如图4），得到如图所示的三维 有限元模型（如图5）；

c)设定该三维有限元模型的计算总时长以及单个时间步大小，并对模型中 的内圆柱和外圆柱分别创建物理模型，赋予单元类型、材料属性等，并 对外圆柱设定角速度，对线圈施加恒定电流；

d)对内圆柱和外圆柱的有限元模型，设置规范方案为拉格朗日乘子规范方 案；

e)施加内圆柱和外圆柱的耦合条件：选择内圆柱的外环面和外圆柱内环面 作为滑移面，采用拉格朗日乘子法耦合磁矢势A，采用MPC方法耦合拉 格朗日标量乘子p自由度；

f)计算由d)规范的棱边单元矩阵并进行集成，计算滑移界面上的拉格朗日 矢量乘子约束矩阵并集成，最后获得总体有限元矩阵。

g)对整个有限元模型的上底面、下底面所有节点，和外圆柱的外环面上所 有节点的棱边电磁势A，以及相应的拉格朗日标量乘子p自由度进行约 束；

h)对施加边界约束的有限元方程进行求解,获得各个时间步的计算结果和 磁场分布。

最后求解结束，获得线圈旋转过程中电磁场的场量，图6和图7为线圈旋 转到不同位置时的磁感应强度的矢量分布图。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局 限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本 发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护 范围之内。