一种基于非结构网格求解避雷针电场的方法及系统

本发明涉及计算电磁学领域，具体涉及一种基于非结构网格求解避雷针电场的方法及系统。

背景技术：

避雷针广泛应用于建筑、电力设备的防雷保护当中。其基本原理是在避雷针的顶端形成畸变的空间电场，从而引导雷电向避雷针放电，然后通过接地引线将雷电引入大地，保护周围物体免遭雷击。避雷针周围电场分布的计算，对于防雷保护、避雷针结构设计等具有重要意义。因此如何能够高效、准确地计算避雷针电场，是需要解决的现实需求。

自21世纪初法国学者tonti提出元胞法以来，元胞法凭借其对于各类网格的普适性和多物理场建模的统一性，已被众多学者广泛应用到电场、磁场、温度场等多种物理场问题的求解当中。但是，当元胞法应用于避雷针电场的求解时，其现有技术仍存在一些明显缺陷尚需解决。首先，元胞法的主要瓶颈在于其只有二阶精度，这导致了该方法难以适应于一些有高精度要求的问题。如果为了追求计算精度，采用较细的网格，则会大量耗费计算内存和时间；而如果为了追求计算效率，采用较粗的网格计算，则在避雷针电场分布非均匀的情况下，本构方程离散误差将显著增加，从而严重影响避雷针电场的计算精度，因此，目前尚缺乏一种能有效兼顾计算效率和精度的数值方法。

技术实现要素：

为了克服现有技术在计算避雷针电场时的缺陷，本发明提供一种基于非结构网格求解避雷针电场的方法及系统，本发明能够显著提升原有方法的计算精度，且在修正过程中网格剖分保持不变，既不增加系数矩阵规模又不增加矩阵带宽，仅需对元胞法本构方程进行重构与修正，从而可以在不加密网格的情况下实现避雷针电场的高精度求解，有效地兼顾了计算效率和精度。

为达到上述目的，本发明采用如下技术方案：

一种基于非结构网格求解避雷针电场的方法，包括如下过程：

利用元胞法计算避雷针电场初始数值解；

利用避雷针电场初始数值解获得避雷针电场求解域内电场的分布函数；

利用电场分布函数对避雷针电场的本构方程中的本构矩阵进行修正，得到全局本构矩阵；

利用全局本构矩阵和避雷针电场的拓扑矩阵更新避雷针电场的系数矩阵，然后重新计算避雷针电场，得到避雷针电场的数值解。

优选的，利用元胞法计算避雷针电场初始数值解的过程包括：

根据避雷针电场中的变量关系，基于元胞法对避雷针进行离散建模，得到对应的tonti图(即托特图)；

确定避雷针电场的求解域，对求解域进行网格剖分，获得主网格和对偶网格信息；

根据主网格和对偶网格信息，计算避雷针电场的拓扑矩阵和本构矩阵；

根据tonti图，将拓扑矩阵和本构矩阵组装为系数矩阵；

根据避雷针电场的边界条件，建立元胞方程组；

求解元胞方程组，获得避雷针电场初始数值解。

优选的，利用避雷针电场初始数值解获得避雷针电场求解域内电场的分布函数的过程包括：

将避雷针电场初始数值解作为离散样本，插值获得求解域内的电位分布函数；

对电位分布函数进行求导，获得求解域内每个网格单元内超收敛点处电位的导数样本；

使用超收敛导数恢复方法，利用导数样本恢复各网格单元内的电场分布函数。

优选的，利用电场分布函数对避雷针电场的本构方程中的本构矩阵进行修正，得到全局本构矩阵的过程包括：

将分布函数代入避雷针电场的本构方程对本构方程中的本构矩阵进行修正，得到修正后的局部本构矩阵；其中，修正后的局部本构矩阵中元素的表达式如下：

式中，表示修正后的局部本构矩阵中第j行第i列的元素，表示恢复得到的电通量，ui表示恢复得到的电位差；ε表示空气介电常数，e^*表示恢复得到的电场分布函数，表示第j个对偶网格面，li表示第i条主网格棱边，ds表示对面积分，dl表示对棱边积分；

将修正后的局部本构矩阵组装为全局本构矩阵。

优选的，利用全局本构矩阵和避雷针电场的拓扑矩阵更新避雷针电场的系数矩阵，然后重新计算避雷针电场，得到避雷针电场的数值解的过程包括：

将避雷针电场的拓扑矩阵与全局本构矩阵组装为新的系数矩阵；

根据避雷针电场已知的边界条件和新的系数矩阵，重新建立元胞方程组；

求解重新建立的元胞方程组，得到避雷针电场的数值解。

本发明还提供了一种基于非结构网格求解避雷针电场的系统，包括：

初始数值解计算模块：用于利用元胞法计算避雷针电场初始数值解；

分布函数计算模块：用于利用避雷针电场初始数值解计算避雷针电场求解域内电场的分布函数；

修正模块：用于利用电场分布函数对避雷针电场的本构方程中的本构矩阵进行修正，得到全局本构矩阵；

电场计算模块：用于利用全局本构矩阵和避雷针电场的拓扑矩阵更新避雷针电场的系数矩阵，然后重新计算避雷针电场，得到避雷针电场的数值解。

优选的，初始数值解计算模块利用元胞法计算避雷针电场初始数值解的过程包括：

根据避雷针电场中的变量关系，基于元胞法对避雷针进行离散建模，得到对应的tonti图；

确定避雷针电场的求解域，对求解域进行网格剖分，获得主网格和对偶网格信息；

根据主网格和对偶网格信息，计算避雷针电场的拓扑矩阵和本构矩阵；

根据tonti图，将拓扑矩阵和本构矩阵组装为系数矩阵；

根据避雷针电场的边界条件，建立元胞方程组；

求解元胞方程组，获得避雷针电场初始数值解。

优选的，分布函数计算模块利用避雷针电场初始数值解计算避雷针电场求解域内电场的分布函数的过程包括：

将避雷针电场初始数值解作为离散样本，插值获得求解域内的电位分布函数；

对电位分布函数进行求导，获得求解域内每个网格单元内超收敛点处电位的导数样本；

使用超收敛导数恢复方法，利用导数样本恢复各网格单元内的电场分布函数。

优选的，修正模块利用电场分布函数对避雷针电场的本构方程中的本构矩阵进行修正，得到全局本构矩阵的过程包括：

将分布函数代入避雷针电场的本构方程对本构方程中的本构矩阵进行修正，得到修正后的局部本构矩阵；其中，修正后的局部本构矩阵中元素的表达式如下：

式中，表示修正后的局部本构矩阵中第j行第i列的元素，表示恢复得到的电通量，ui表示恢复得到的电位差；ε表示空气介电常数，e^*表示恢复得到的电场分布函数，表示第j个对偶网格面；li表示第i条主网格棱边，ds表示对面积分，dl表示对棱边积分；

将修正后的局部本构矩阵组装为全局本构矩阵。

优选的，电场计算模块利用全局本构矩阵和避雷针电场的拓扑矩阵更新避雷针电场的系数矩阵，然后重新计算避雷针电场，得到避雷针电场的数值解的过程包括：

将避雷针电场的拓扑矩阵与全局本构矩阵组装为新的系数矩阵；

根据避雷针电场已知的边界条件和新的系数矩阵，重新建立元胞方程组；

求解重新建立的元胞方程组，得到避雷针电场的数值解。

与现有技术相比，本发明具有如下有益效果：

本发明基于非结构网格求解避雷针电场的方法通过将避雷针电场初始数值解作为先验信息，获得求解域内的电场分布函数，然后利用电场分布函数对避雷针电场的本构方程中的本构矩阵进行修正，最后利用全局本构矩阵和避雷针电场的拓扑矩阵更新避雷针电场的系数矩阵，然后重新计算避雷针电场，获得高精度数值解。从上述可以看出，相比于传统的提高精度的自适应方法，本发明既不增加系数矩阵规模又不增加矩阵带宽，仅需对本构矩阵进行重构更新，修正了方程组中的权重系数，从而可以在不加密网格的情况下实现高精度求解，有效地兼顾了计算精度和效率。并且，该方法可同时适用于结构网格和非结构网格。

附图说明

图1本发明基于非结构网格求解避雷针电场的方法的计算流程示意图；

图2本发明实施例的避雷针模型示意图；

图3本发明实施例对应的元胞法tonti图；

图4本发明实施例计算出的电位分布图；

图5本发明实施例中采用软件仿真的电位分布图。

具体实施方式

下面结合附图和实施例将详细描述本发明的实施例，需要指出的是，所述实施例是示例性的，旨在用于解释本发明，而并非是对本发明的实施方式的限定。

本发明基于非结构网格求解避雷针电场的方法是一种基于非结构网格下元胞法本构矩阵自适应修正的避雷针电场计算方法，参照图1，包括如下步骤：

s1，基于元胞法对避雷针电场边值问题进行离散建模，求解元胞方程组获得模型初始数值解；具体包括以下步骤：

s1-1，根据避雷针电场中的变量关系，基于元胞法进行离散建模，得到对应问题的tonti图，参照图3；

s1-2，确定避雷针电场的求解域，对求解域进行网格剖分，获得主网格和对偶网格信息；

s1-3，依据网格剖分信息(即主网格和对偶网格信息)，计算避雷针电场的拓扑矩阵g和本构矩阵m。此时的本构矩阵m中的元素计算公式如下：

其中mj,i表示本构矩阵中第j行第i列的元素，φj表示穿过第j个对偶网格面的电通量，δui表示沿第i条主网格棱边的电位差；ε表示空气介电常数；表示第j个对偶网格面的面积；li表示第i条主网格棱边的长度。

s1-4，依据tonti图，将拓扑矩阵和本构矩阵组装为系数矩阵；

s1-5，根据已知的避雷针电场边界条件，强加边界条件，建立元胞方程组；

s1-6，求解元胞方程组，获得避雷针电场初始数值解。

s2，以元胞法求得的避雷针电场初始数值解作为先验信息，使用超收敛导数恢复技术获得求解域内电场的分布函数；具体包括以下步骤：

s2-1，以元胞法求得的避雷针电场初始数值解作为离散样本，插值获得求解域内的电位分布函数；

s2-2，对获得的电位分布函数进行求导，获得求解域内每个网格单元内超收敛点处电位的导数样本；

s2-3，使用超收敛导数恢复方法，利用上述导数样本恢复各网格单元内的电场分布函数。

s3，将恢复得到的电场分布函数代入到元胞法本构方程(即雷针电场本构方程)中进行重构与修正；具体包括以下步骤：

s3-1，将恢复得的电场分布函数代入到元胞法本构方程中的本构矩阵进行修正，该修正是基于本构关系是元胞法离散误差主要来源，而修整过程既不增加系数矩阵规模又不增加矩阵带宽；本构矩阵中的局部本构矩阵修正后的元素计算公式如下：

其中表示修正后的局部本构矩阵中第j行第i列的元素，表示恢复得到的电通量，δui表示恢复得到的电位差；ε表示空气介电常数，e^*表示恢复得到的电场分布函数，ds表示对面积分，dl表示对棱边积分。

s3-2，将修正后的局部本构矩阵组装为全局本构矩阵。

s4，更新元胞方程组的系数矩阵，重新计算得到避雷针电场的高精度数值解，具体包括以下步骤：

s4-1，将拓扑矩阵与全局本构矩阵组装为新的系数矩阵；

s4-2，根据已知的避雷针电场边界条件，强加边界条件，并利用系数矩阵重新建立元胞方程组；

s4-3，求解重新建立的元胞方程组，得到避雷针电场的高精度数值解。

从上述可以看出，本发明通过将元胞法的初始数值解作为先验信息，使用恢复技术获得求解域内的电场分布函数，然后代入到本构矩阵当中进行重构修正，最后重新计算更新后的方程组，获得高精度数值解。相比于传统的提高精度的自适应方法，本发明既不增加系数矩阵规模又不增加矩阵带宽，仅需对本构矩阵进行重构更新，修正了方程组中的权重系数，从而可以在不加密网格的情况下实现高精度求解，有效地兼顾了计算精度和效率。并且，该方法可同时适用于结构网格和非结构网格。

实施例

本实施例以图2所示的避雷针为例，计算该避雷针周围的电场分布。其中边界条件为避雷针表面和底部都接地，满足顶部处电位设为恒定值右侧边界使用渐进边界条件，求解域内的控制方程为

本实施例基于非结构网格求解避雷针电场的方法的具体过程包括以下步骤：

(1)基于元胞法对避雷针电场边值问题进行离散建模，求解元胞方程组获得模型初始数值解；具体包括以下步骤：

1-1)根据该实施例避雷针电场中所涉及的变量关系，基于元胞法进行离散建模，得到对应的tonti图，如图3所示。图中的表示电位，u表示电位差，φe表示电通量，g表示梯度拓扑矩阵，d％表示散度拓扑矩阵，m表示本构矩阵；

1-2)确定避雷针电场的求解域。对求解域进行网格剖分，获得主网格和对偶网格信息；

1-3)依据网格剖分信息，计算该避雷针电场的拓扑矩阵g和本构矩阵m。此时本构矩阵m中的元素计算公式如下：

其中，mj,i表示本构矩阵m中第j行第i列的元素，表示穿过第j个对偶网格面的电通量，ui表示沿第i条主网格棱边的电位差；ε表示介电常数；表示第j个对偶网格面的面积；li表示第i条主网格棱边的长度。

1-4)依据tonti图，将拓扑矩阵和本构矩阵组装为系数矩阵。

1-5)根据已知的避雷针电场边界条件，强加边界条件，建立元胞方程组如下：

1-6)求解元胞方程组，获得避雷针电场初始数值解。

(2)以元胞法求得的避雷针电场初始数值解作为先验信息，使用超收敛导数恢复技术获得求解域内电场的分布函数；具体包括以下步骤：

2-1)以元胞法求得的避雷针电场初始数值解作为离散样本，插值获得求解域内的电位分布函数此时插值获得的电位分布函数表达式为：

其中表示多项式的系数，r1为x多项式的幂次参数，r2为y多项式的幂次参数，x为插值函数的横坐标自变量，x0为离散样本的横坐标，y为插值函数的纵坐标自变量，y0为离散样本的纵坐标。

2-2)，对获得的电位分布函数进行求导，获得求解域内每个网格单元内超收敛点处电位的导数样本

2-3)，使用超收敛导数恢复技术，利用上述导数样本恢复获得求解域内电场强度的分布函数如下：

其中，e^*为恢复得到电场强度分布函数，r为导数恢复算子。

(3)将恢复得到的电场分布函数代入到元胞法本构方程中进行重构与修正；具体包括以下步骤：

3-1)将恢复得的电场分布函数代入到元胞法本构方程的本构矩阵中进行修正，此时本构矩阵中的局部本构矩阵修正后的元素计算公式如下：

其中，表示修正后的局部本构矩阵中第j行第i列的元素，表示恢复得到的电通量，ui表示恢复得到的电位差；ε表示空气介电常数，e^*表示恢复得到的电场分布函数，ds表示对面积分，dl表示对棱边积分。

3-2)，将修正后的局部本构矩阵组装为全局本构矩阵。

(4)更新元胞方程组的系数矩阵，重新计算得到避雷针电场的高精度数值解，具体包括以下步骤：

4-1)将拓扑矩阵与更新后的本构矩阵(即全局本构矩阵)组装为新的系数矩阵；

4-2)根据已知的避雷针电场边界条件，强加边界条件，并利用系数矩阵重新建立元胞方程组；

其中，m^*表示修正后的本构矩阵，g^t为拓扑矩阵g的转置矩阵。

4-3)求解重新建立的元胞方程组，输出结果，得到避雷针电场的高精度数值解。

本实施例的计算结果如图4所示，商业软件仿真结果如图5所示。由图4和图5可以看出，本实施例的避雷针电场计算结果与仿真结果高度一致。这说明本发明在没有大幅加密网格的情况下，既不增加系数矩阵规模又不增加矩阵带宽，仅需对元胞法本构方程进行重构与修正，便实现避雷针电场的高精度求解，有效地兼顾了计算效率和精度，从而证明了本发明的先进性。