一种用于电磁场求解的自适应区域分解有限元方法与流程

1.本发明涉及电气设备及仿真计算的技术领域，尤其涉及一种用于电磁场求解的自适应区域分解有限元方法。


背景技术：

2.随着我国工业生产自动化水平的不断提高，电气设备在工业应用中的地位越来越重要，准确高效的数值分析方法对于电气设备的优化设计和安全运行起着关键的作用。自适应有限元方法，即在有限元法基础上引入误差分析和离散调节，以其准确性、灵活性和高效性，已经成为当前电磁分析的主导数值方法，被广泛应用在求解实际工程问题中。基于特定问题求解过程中误差分析的结果，该方法能够产生满足精度要求的最优网格，其中所含单元数目尽可能少且误差近乎均匀分布。根据所使用的离散调节方法，自适应有限元法被分成三类，即r
‑
法，p
‑
法和h
‑
法。在h
‑
法中，按照误差分析的结果，网格被有序的加密或稀疏，最终生成一个误差接近均匀分布的最优网格；在求解过程中，针对不同的区域，加密和稀疏操作都会有涉及，以实现用尽可能少的计算资源得出足够精确的数值解。
3.由于在电磁分析中多涉及非线性材料，且部分问题结构复杂，时间跨度长，导致计算规模很大，需要大量的计算资源来处理网格剖分以及后续的代数系统求解，最终影响到电磁分析及设计的效率。随着计算机硬件和并行技术的进步，大规模工程问题得以通过把区域分解理念引入到有限元法来求解；区域分解方法的理念是把求解域分解成若干个重叠的或不重叠的子区域，并在对各子问题的边界条件进行合理处置后，使用并行处理单元分别求解各子问题；在区域分解有限元方法形成的代数系统求解中，多使用krylov子空间法结合合适的预处理算子以提高计算效率；区域分解法已经被应用在求解低频电磁场问题中，有效地提高了问题的求解效率，但是传统的区域分解有限元方法的仍存在计算精度相对不足、求解效率相对有限等问题。


技术实现要素：

4.本部分的目的在于概述本发明的实施例的一些方面以及简要介绍一些较佳实施例。在本部分以及本技术的说明书摘要和发明名称中可能会做些简化或省略以避免使本部分、说明书摘要和发明名称的目的模糊，而这种简化或省略不能用于限制本发明的范围。
5.鉴于上述现有存在的问题，提出了本发明。
6.因此，本发明解决的技术问题是：传统区域分解有限元方法的计算精度相对不足、求解效率相对有限等问题。
7.为解决上述技术问题，本发明提供如下技术方案：基于对所求解电磁场问题无先验知识，对所述电磁场整个求解域进行一致剖分，生成初始离散网格；根据目标子区域的数目和所述初始离散网格，对求解任务进行均匀分配并调用区域分解求解器对离散后得到的电磁问题进行求解；基于求解结果，计算用来描述误差信息和三角形网格单元形状尺寸特征的hessian矩阵和度量张量矩阵；对计算数值结果进行误差分析，若误差不满足预设值，
则对网格做自适应处理并再次计算，直到误差满足预设值并输出数值结果。
8.作为本发明所述的用于电磁场求解的自适应区域分解有限元方法的一种优选方案，其中：对所述电磁场整个求解域进行区域划分包括，基于求解问题离散网格按计算量将求解域合理的划分为多个子求解域。
9.作为本发明所述的用于电磁场求解的自适应区域分解有限元方法的一种优选方案，其中：所述电磁场其基于磁势函数的暂态控制方程包括，
[0010][0011]
其中，ν为磁阻率，u为磁势，σ为电导率，t为时间，j为电流密度。
[0012]
作为本发明所述的用于电磁场求解的自适应区域分解有限元方法的一种优选方案，其中：基于伽辽金有限元离散后得到的所述电磁问题完全展开的代数方程包括，
[0013][0014]
其中，对角线块矩阵(a
ij
,i＝j)表征子区域i内自由度对于系数矩阵的贡献，非对角线块矩阵(a
ij
,i≠j)表征子区域j内自由度对于子区域i系数矩阵的贡献，u
i
为子区域i的未知量矩阵，f
i
为子区域i的载荷矩阵。
[0015]
作为本发明所述的用于电磁场求解的自适应区域分解有限元方法的一种优选方案，其中：基于有限元策略生成的离散方程按照子区域进行分解，各子区域和整个定义域内未知量的迭代方程包括，基于雅克比法，所述未知量的迭代方程表示为：
[0016][0017]
该方程等价于：
[0018][0019]
其中：
[0020][0021]
基于：
[0022][0023][0024]
其中，上标n 1/n为迭代计算步数，r
i
为从整个定义域到子区域i的限制算子，为m
i
×
m的矩阵，m
i
和m分别为各子区域和整个定义域内未知量的个数；
[0025]
所述未知量的迭代方程最终写为：
[0026][0027]
其中，
[0028][0029]
r
n
＝f
‑
au
n
[0030]
其中，a为电磁场刚度矩阵，u为磁势未知量矩阵，f为电磁场载荷矩阵，上标t为矩阵转置。
[0031]
作为本发明所述的用于电磁场求解的自适应区域分解有限元方法的一种优选方案，其中：所述hessian矩阵包括，hessian矩阵h
k
表示为：
[0032]
h
k
＝{h
k,ij
}；i/j＝1,2
[0033][0034]
其中，u
k
为磁势的第k次迭代解，x1、x2分别代表二维空间x和y两个方向的自变量。
[0035]
作为本发明所述的用于电磁场求解的自适应区域分解有限元方法的一种优选方案，其中：所述h
k
对应的度量张量矩阵m
k
包括，
[0036]
m
k
＝v
t
diag(λ)v
[0037]
其中，λ和v分别是h
k
的特征值与特征向量。
[0038]
作为本发明所述的用于电磁场求解的自适应区域分解有限元方法的一种优选方案，其中：对所述计算数值结果进行误差分析包括，在由误差信息决定的所述度量张量矩阵m
k
的几何空间内，网格质量由下式判断：
[0039][0040]
其中，|δ|是单元的面积，p(δ)是单元的周长，n*是网格单元期望数目，h*(n*)是网格单元的平均尺寸，f(t)(0≤f≤1)是极值点在t＝1处的任意凸函数。
[0041]
本发明的有益效果：本发明提出了一种将自适应方法和区域分解方法优点结合起来的自适应区域分解有限元方法，能够在保证求解精度的基础上大幅提高计算效率；使用灵活的加性schwarz区域分解法作为krylov子空间法的预处理算子，并引入基于度量张量
的各向异性网格自适应方法，对网格质量进行分析后可以根据磁场分布生成最优网格；与传统区域分解法相比，本发明所提出的方法在保证计算精度的同时，降低了所需自由度数目，求解效率有显著的提高。
附图说明
[0042]
为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其它的附图。其中：
[0043]
图1为本发明一个实施例提供的一种用于电磁场求解的自适应区域分解有限元方法的基本流程示意图；
[0044]
图2为本发明一个实施例提供的一种用于电磁场求解的自适应区域分解有限元方法的c型执行器的结构和尺寸示意图；
[0045]
图3为本发明一个实施例提供的一种用于电磁场求解的自适应区域分解有限元方法的c型执行器的磁场能量收敛曲线示意图；
[0046]
图4为本发明一个实施例提供的一种用于电磁场求解的自适应区域分解有限元方法的c型执行器的网格自适应过程示意图；
[0047]
图5为本发明一个实施例提供的一种用于电磁场求解的自适应区域分解有限元方法的c型执行器的磁势填充等势线示意图；
[0048]
图6为本发明一个实施例提供的一种用于电磁场求解的自适应区域分解有限元方法的永磁同步电机模型示意图；
[0049]
图7为本发明一个实施例提供的一种用于电磁场求解的自适应区域分解有限元方法的永磁同步电机模型的离散网格及区域划分示意图；
[0050]
图8为本发明一个实施例提供的一种用于电磁场求解的自适应区域分解有限元方法的永磁同步电机模型的磁势填充等势线示意图。
具体实施方式
[0051]
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合说明书附图对本发明的具体实施方式做详细的说明，显然所描述的实施例是本发明的一部分实施例，而不是全部实施例。基于本发明中的实施例，本领域普通人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明的保护的范围。
[0052]
在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是本发明还可以采用其他不同于在此描述的其它方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似推广，因此本发明不受下面公开的具体实施例的限制。
[0053]
其次，此处所称的“一个实施例”或“实施例”是指可包含于本发明至少一个实现方式中的特定特征、结构或特性。在本说明书中不同地方出现的“在一个实施例中”并非均指同一个实施例，也不是单独的或选择性的与其他实施例互相排斥的实施例。
[0054]
本发明结合示意图进行详细描述，在详述本发明实施例时，为便于说明，表示器件结构的剖面图会不依一般比例作局部放大，而且所述示意图只是示例，其在此不应限制本
发明保护的范围。此外，在实际制作中应包含长度、宽度及深度的三维空间尺寸。
[0055]
同时在本发明的描述中，需要说明的是，术语中的“上、下、内和外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一、第二或第三”仅用于描述目的，而不能理解为指示或暗示相对重要性。
[0056]
本发明中除非另有明确的规定和限定，术语“安装、相连、连接”应做广义理解，例如：可以是固定连接、可拆卸连接或一体式连接；同样可以是机械连接、电连接或直接连接，也可以通过中间媒介间接相连，也可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。
[0057]
实施例1
[0058]
参照图1，为本发明的一个实施例，提供了一种用于电磁场求解的自适应区域分解有限元方法，包括：
[0059]
s1：基于对所求解电磁场问题无先验知识，对电磁场整个求解域进行一致剖分，生成初始离散网格；
[0060]
s2：根据目标子区域的数目和初始离散网格，对求解任务进行均匀分配并调用区域分解求解器对离散后得到的电磁问题进行求解；
[0061]
s3：基于求解结果，计算用来描述误差信息和三角形网格单元形状尺寸特征的hessian矩阵和度量张量矩阵；
[0062]
s4：对计算数值结果进行误差分析，若误差不满足预设值，则对网格做自适应处理并再次计算，直到误差满足预设值并输出数值结果。
[0063]
s1～s4步骤具体包括：
[0064]
对电磁场整个求解域进行区域划分包括：
[0065]
在加性schwarz区域分解法中引入基于度量张量的各向异性网格自适应方法，将h
‑
型自适应有限元方法与区域分解方法的优点结合起来，在求解过程中生成对场分布响应更好的网格；
[0066]
具体的，采用基于求解问题离散网格的方法对计算量进行合理的划分，如metis工具；由于不受几何尺寸的限制，该方法能更容易适应实际问题离散中存在的网格尺寸差异；此外，metis工具能够更加充分的利用计算资源，具有更高的并行率。在生成的非重叠子域基础上，可以通过对每个子域向外增加一层或者几层单元来生成重叠的子域；因此，针对非重叠区域分解算法所做的区域划分也可以轻松的适用于重叠区域分解算法。
[0067]
进一步的，加性schwarz区域分解法包括：
[0068]
在使用有限元方法求解电磁场问题时，多通过引入势函数并结合材料的本构方程对原有的麦克斯韦方程进行变形以简化计算，而后基于定义域的离散网格及其形状函数对势函数进行有限元离散，并结合边界条件形成最后的代数系统；在区域分解有限元方法中，各子区域构建的子问题经过边界未知量及相应的边界条件耦合，并最终借助不同的求解方法实现并行求解。
[0069]
基于加性schwarz的区域分解法是一种通用、高效且稳定的区域分解方法，可以针对不同的问题对子区域的重叠面积进行调整，并能实现稳定收敛，接下来结合电磁场控制
方程对其基本原理进行介绍；对于二维问题，其基于磁势函数的暂态磁场控制方程可以表示为：
[0070][0071]
其中，ν为磁阻率，u为磁势，σ为电导率，t为时间，j为电流密度。
[0072]
通过伽辽金有限元离散后得到电磁问题完全展开的代数方程可以表示为：
[0073]
au＝f
[0074]
其中，a为电磁场刚度矩阵，u为磁势未知量矩阵，f为电磁场载荷矩阵。
[0075]
在原求解域被分解为n个子求解域后，势函数向量被相应的分成n个子向量，同时基于有限元方法生成的离散方程可以按照子区域进行分解，各子区域和整个定义域内未知量的个数分别用m
i
和m来表示，因而通过伽辽金有限元离散后得到电磁问题完全展开的代数方程可以重写为：
[0076][0077]
其中，对角线块矩阵(a
ij
,i＝j)表征子区域i内自由度对于系数矩阵的贡献，非对角线块矩阵(a
ij
,i≠j)表征子区域j内自由度对于子区域i系数矩阵的贡献，u
i
为子区域i的未知量矩阵，f
i
为子区域i的载荷矩阵。
[0078]
基于雅克比法，未知量的迭代方程表示为：
[0079][0080]
该方程等价于：
[0081][0082]
其中：
[0083][0084]
基于：
[0085][0086][0087]
其中，上标n 1/n为迭代计算步数，r
i
为从整个定义域到子区域i的限制算子，为m
i
×
m的矩阵，m
i
和m分别为各子区域和整个定义域内未知量的个数；
[0088]
未知量的迭代方程最终写为：
[0089][0090]
其中，
[0091][0092]
r
n
＝f
‑
au
n
[0093]
其中，上标t为矩阵转置。
[0094]
在计算过程中，每个子区域对应的对角线块矩阵的求逆矩阵计算可以利用并行处理器实现，为提高未知量的迭代方程的求解效率以及避免潜在的收敛问题，经区域分解后形成的代数系统可使用krylov子空间法结合预处理算子进行求解，在本发明中采用共轭梯度法。
[0095]
由于电磁分析中众多自由度的准确解是不可知的，对于数值结果进行有效的误差评估是必不可少的，同时也需要针对评估结果对离散进行调整以获得满足精度的结果及避免计算资源浪费；因此，本发明在加性schwarz区域分解算法中引入网格自适应算法来保证求解的精度及提高计算效率；为避免子区域网格独立调整可能导致的悬点，在每步计算完成后基于计算结果统一对整个求解域的网格进行调整；为得到对磁场分布响应更好的网格，本发明采用基于度量张量的各向异性网格自适应方法来对求解过程中的网格进行调整，其他技术与经典自适应有限元方法保持一致。
[0096]
其中，基于求解结果，计算用来描述误差信息和三角形网格单元形状尺寸特征的hessian矩阵和度量张量矩阵包括：
[0097]
hessian矩阵h
k
表示为：
[0098]
h
k
＝{h
k,ij
}；i/j＝1,2
[0099][0100]
其中，u
k
为磁势的第k次迭代解，x1、x2分别代表二维空间x和y两个方向的自变量；
[0101]
进一步的，h
k
对应的度量张量矩阵m
k
包括：
[0102]
m
k
＝v
t
diag(λ)v
[0103]
其中，λ和v分别是h
k
的特征值与特征向量。
[0104]
如果系统误差满足期望值，求解结束输出数值结果，否则对网格做进一步自适应
处理并再次计算；在由误差信息决定的度量张量m
k
的几何空间内，网格质量由下式判断：
[0105][0106]
其中，|δ|是单元的面积，p(δ)是单元的周长，n*是网格单元期望数目，h*(n*)是网格单元的平均尺寸，f(t)(0≤f≤1)是极值点在t＝1处的任意凸函数。
[0107]
更进一步的，基于度量张量，通过移动节点、交换边、增删节点对网格进行自适应处理，以生成在度量张量定义的几何空间内的最优网格剖分；经过网格自适应后，数值解精度和自由度数量的比率得到显著的提高，而后返回对求解域进行区域分解并求解步骤进行再次计算，直到误差满足预设值并输出数值结果；对于暂态问题，在上述循环的基础上增加步进计算步骤，并可与自适应时间步长调整方法结合，在满足精度要求的前提下用尽可能少的时间步数完成计算。
[0108]
本发明在加性schwarz区域分解法中引入基于度量张量的各向异性网格自适应方法，将h
‑
型自适应有限元方法与区域分解方法的优点结合起来，在求解过程中生成对场分布响应更好的网格；采用基于度量张量的各向异性网格自适应策略，可以基于当前网格的数值解对网格质量进行判断，并进行后续的各向异性调整以获得最优离散系统；计算过程中的多次子区域划分通过metis实现，按照子区域数量和当前网格对求解任务进行近似平均的分配；为进一步提高求解效率，在代数方程求解中将区域分解法用作krylov子空间法的预处理算子；使用有限元求解缺乏先验知识的电磁问题，可以先对整个求解域进行粗略的均匀剖分，后续再依据误差分析的评估结果对离散网格进行自适应调整；本发明提出的方法遵循同样的逻辑，在自适应离散后通过区域分解有限元方法进行求解；为简化计算步骤及避免分散网格操作可能导致的悬点(子区域间不匹配网格)，区域分解求解后的误差评估以及离散网格自适应调整均在全局层面进行；因此，本发明提出一种将自适应方法和区域分解方法优点结合起来的自适应区域分解有限元方法，能够在保证求解精度的基础上大幅提高计算效率；使用灵活的加性schwarz区域分解法作为krylov子空间法的预处理算子，并引入基于度量张量的各向异性网格自适应方法，对网格质量进行分析后可以根据磁场分布生成最优网格；与传统区域分解法相比，本发明方法在保证计算精度的同时，降低了所需自由度数目，求解效率有显著的提高。
[0109]
实施例2
[0110]
参照图2～8为本发明另一个实施例，为对本方法中采用的技术效果加以验证说明，本实施例采用传统技术方案与本发明方法进行对比测试，以科学论证的手段对比试验结果，以验证本方法所具有的真实效果。
[0111]
如图2所示为一个c型执行器的结构和尺寸图，浅色区域为流通1000安匝电流的绕组，在数值问题中采用一个包含执行器的求解域，并在该求解域内生成一套均匀的初始网格，其中包含3076个三角形单元和1617个自由度，在后续的计算中将分别采用传统区域分解法和自适应区域分解法来求解该问题，并就效率和精度对两种方法进行比较。
[0112]
在使用传统区域分解法进行求解的过程中，通过整体限制单元的尺寸来增加自由度的个数，进而提高解的精度，在经过八次整体加密后，所求得的磁场能量基本达到稳定，如图3所示，最终生成的网格含有33760个三角形单元和17141个自由度，由于采用一致的剖
分原则，加密过程同时会对计算结果误差小的区域进行加密，不可避免会造成一定的计算资源浪费。
[0113]
使用自适应区域分解法求解该问题，在初始网格的基础上经过四次自适应操作后所求解的磁场能量与区域分解法的结果基本相当，如图3所示，在求解过程中，自适应方法基于误差分析结果对当前离散网格进行调整，在误差大的区域增加节点并在误差小的区域移除节点，以保证在避免计算资源浪费的前提下得出满足精度要求的结果。具体的网格自适应过程及其对应的子区域划分在图4中给出，最终生成的网格有16787个三角形单元和8446个自由度，其中自由度数目约减少50％，观察生成的网格可以发现，磁导率突变边界及尖角周围的网格得到加密，而磁场梯度较小区域的网格尺寸相对较大；对于该电磁问题来说，自适应区域分解法只需传统区域分解法所需计算资源的二分之一左右就能得出接近相同精度的结果；以区域分解法所得磁势结果作为参考值，自适应区域分解法所得磁势的l2误差为2.4
×
10
‑3；为更好的展示求得的结果，磁势的填充等势线在图5中给出。
[0114]
图6为一个四对极二十四槽的表贴式永磁同步电机的二维模型，分别采用传统区域分解有限元方法和本发明所提的自适应区域分解有限元方法对其空载运行状态进行分析；在初始均匀网格基础上经过七次均匀网格加密收敛后的离散网格及其相应的区域划分如图7(a)所示，其中包含11788个三角形单元及5970个自由度，求解域基于离散网格被划分成十五个子求解域，基于该网格求解的磁势分布如图8(a)所示；采用自适应区域分解有限元方法求解时经过五次自适应调整生成的网格及所求结果分别如图7(b)和8(b)所示，通过比较可以发现两种方法所得结果具有良好的一致性，且以区域分解法所得结果作为基准值，自适应区域分解法计算所得磁势的l2误差为3.1
×
10
‑3，验证了本发明所提方法的准确性；最终自适应调整后的网格单元数和自由度数分别为7512和3801，远小于均匀加密网格各相应参数，其中自由度数目约减少36％，计算规模得到显著降低。
[0115]
应说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。