一种基于区域分解的快速扫频仿真方法与流程

本发明涉及微波工程中的快速扫频仿真领域，尤其涉及一种基于区域分解的快速扫频仿真方法。
背景技术：
：在微波工程领域，一段宽频带内对电大尺寸以及几何结构复杂微波器件的精确电磁仿真是不可避免的任务。传统的有限元方法对器件模型进行扫频时，对每个频点都要进行重新建模并且重复求解，当模型尺寸过于复杂或者剖分精度要求过高时，未知量就会过多，采用直接解法计算机内存往往达不到要求，采用迭代解法也需要很长的时间才能收敛。前人将模型降阶法结合到有限元方法中，建立了有限元模型降阶法，将重复求解原模型矩阵方程转化成重复求解降阶后的小型的矩阵方程，大大提高了扫频的效率。但究其降阶的本质，是对与原始有限元矩阵同等规模的矩阵循环求解若干次，因此当遇到未知量较大结构复杂的模型时，直接解法不能胜任降阶的任务，只能采用迭代求解，因此降阶时间同样会非常慢。相比较于传统的有限元法，区域分解法是一种非常适合应用在大规模复杂结构尤其是周期型结构进行大规模数值计算的方法。该算法把所要求解的区域分解成若干子域，在各个子内进行单独计算，子域和子域之间通过边界连续性条件进行迭代传输运算，将原大规模矩阵的求解转化为众多子域独立求解，有缩小计算规模和预处理迭代效率高等优点，因此求解速度很快。为了能够使有限元模型降阶法适用更多结构，我们将区域分解法和有限元模型降阶法相结合，极大提高了对大规模复杂结构尤其是周期型结构模型降阶的效率，相比于单纯的有限元模型降阶法适应性更广，相比于传统的有限元法扫频速度更快，而且有时间精度都可控的特点。技术实现要素：本发明的目的在于提供一种基于区域分解的快速扫频仿真方法。这种方法利用共形分割面相邻子域之间拥有公共棱边信息，各相邻子域之间能够互相传输的特性对有限元模型进行降阶操作，在保证结果精度高的同时与传统的有限元法扫频相比速度更快，并且由于采用了区域分解的技术，提高模型降阶的效率，从而使有限元模型降阶适应性更广。实现本发明的技术解决方案是：一种基于区域分解的快速扫频仿真方法，其包括如下步骤：第一步：将要仿真的模型划分为若干个不重叠的子区域，从而将大规模复杂问题划分成了一系列小规模简单问题，便于提高求解效率。划分子域时，需尽可能地保证每个区域的尺寸相同以获得最佳求解效率，且子域尺寸不超过半个最小扫频波长，以获得性态良好的子域有限元矩阵。利用四面体单元对该模型进行离散；第二步：针对各个子域建立电场矢量波动方程，并在各子域交界面上引入电磁场传输条件，利用whitney基函数构建出各子域在频率展开点下的有限元系数矩阵方程，并求解出各子域内部关于棱边的电场未知向量；第三步：相邻子域之间传输关于棱边的向量信息，在该频率展开点下构建出全局有限元系统krylov子空间的每一列；第四步：利用构建出的krylov子空间，对全局有限元系统不含频率项的有限元系数矩阵进行降阶；第五步：将频率展开点对应频率区间内的各个频点加入降阶后的有限元矩阵方程，采用直接求解器求解，并利用krylov子空间依次恢复出频率展开点对应频带范围各个频点上的全局棱边电场值。进一步的，所述第二步中，根据精度和效率要求确定频率展开点的个数，若设置k个展开点，则将扫频范围等分为k个区间，每个展开点选取为该区间的中点，以确保在等宽的展开范围内精度一致。进一步的，所述第三步中，在有限元模型降阶的过程中，每形成krylov子空间的一列都涉及到对与原始有限元矩阵同等规模的矩阵求一次逆，其具体矩阵求逆的流程如下：首先在采用区域分解法将各个子域的未知向量求出；然后，相邻子域之间传输交换信息，将各子域的未知向量组合成全局向量并存储到krylov子空间对应列中；若要进行下一列的求解运算，则将krylov子空间当前列按照子区域棱边号分离并发送到各个子域，再经过各个子域向量运算，构建出新的各个子域的方程右边向量，在系数矩阵不变的情况下，同样采用区域分解法求解出各个子域的新的未知向量，经过各子域传输组合成下一列全局向量，依次进行即可构建出全局有限元系统的krylov子空间。本发明与现有技术相比，其显著效果是：(1)极大提高了对微波器件快速扫频仿真的效率，相比于传统方法具有更高的效率；(2)本发明方法极易在并行计算平台上实现，能够利用分布式并行计算机大大提高扫频仿真效率，尤其适用于大规模复杂的微波器件扫频仿真；(3)扫频仿真时间和精度具有良好的可控性，通过设置降阶模型的维度，从而可以在确保满足精度的前提下缩短扫频时间。附图说明图1是本发明方法的整体过程流程图；图2是本发明实施例1中仿真对象波导滤波器的具体结构示意图；图3是本发明实施例1中采用本发明方法与传统方法在10-15ghz频段计算出的s参数曲线对比。具体实施方式下面结合附图和实施例对本发明的技术方案作进一步的说明。图1给出了本发明所述方法的主要流程，下面结合图1对本发明所述方法的步骤作进一步详细描述：s1：将要仿真的模型划分为若干个不重叠的子区域，划分时要注意尽量保证各子域大小相同，来保证并行计算时不同进程之间的负载均衡；利用四面体网格单元对整个模型进行剖分，保证每个区域剖分密度合理均匀；s2：对全局棱边按区域号的顺序进行依次排序并将各个区域棱边信息发送到各个子域，使得每个子域拥有对应子区域所有棱边信息，而相邻子区域通过dirichlet边界条件和阻抗边界条件连接；s3：各个子域根据电场矢量波动方程作为控制方程，利用whitney基函数建立对应区域棱边的si矩阵和不含频率项的ti矩阵，对求解模型激励利用dirichlet边界条件构建每个区域不含频率项的右边向量bi；s3.1：通过whitney基函数建立对应区域棱边的si矩阵和不含频率项的ti矩阵的关键如下：其中，nprog表示划分的区域个数；i表示区域号；为第i个子域划分到的四面体总数；和分别为单元矩阵和中的矩阵元素，m,n为六面体棱边序号；s3.2：对求解模型激励利用dirichlet边界条件构建每个区域不含频率项的右边向量bi关键如下：其中，einc代表激励源；s4：通过si,ti,bi，选取合适的频率展开点，构建出各子域该频率展开点下的相关系数矩阵和右边向量，采用区域分解迭代求解器，求解出各个子域内棱边的未知矢量。s4.1：若令s＝jk0(s4.1.1)则子域i的系数矩阵可以表示为：ai＝si+s2ti(s4.1.2)子域i的右边向量可以赋值为：bi＝s·bi(s4.1.3)s4.2：根据区域分解法中的撕裂对接算法，利用dirichlet边界条件连接各子域内的线性系统，有如下形式：ad·xd＝bd(s4.2.1)其中，ad,xd,bd分别具有如下分块的形式：式中，各子域内部的棱边和交界面上与两个子域相连的棱边用上标r表示，与多个子域相连拐角上的棱边用上标c表示，因此，式(s4.2.2)中的表示的是子域i内部及边界上的棱边相关的矩阵，即式(s4.1.2)中的ai，是与拐角处的棱边相关的矩阵，表示两者之间的耦合，根据各子域棱边信息便可求出矩阵称作布尔矩阵，表示子域i交界面上的棱边在整个子域内部与交界面上的编号的映射关系，比如某一棱边在交界面上的编号为a，在子域i内的编号为b，则否则同理；表示子域i内与拐角有关的矩阵元素在整个拐角集合上的投影。式(s4.2.4)中的表示子域i的右边向量，即式(s4.1.3)中的bi。式(s4.2.3)中的是i子域内部的棱边和交界面上的棱边的电场值，λ是定义在子域交界面上的未知矢量，将参与矩阵方程联合求解。因此，在ai已知，右边向量已知以及布尔矩阵已知的情况下，只有xd是唯一未知矢量。利用区域分解迭代求解器即可求出交界面上的场值即λ继而求出每个子域的再将每个子域的剔除与第i-1个子域分界面棱边场值之外剩余已排序好的棱边场值，最后通过各子域之间的传输组合，即可形成真正有关全局的待求矢量xn，具体表现为如下形式：xn＝gather(xn+2)(s4.2.6)在之前步骤中已将所有棱边按照区域号排序，因此在分界面是共形的前提下，相邻的两个子域在分界面有着相同的棱边顺序，因而只要找到两个区域在分割面上边的对应关系然后就可以直接进行传值。式中表示除去与第i-1个子域分界面棱边场值之外剩余已排序好的棱边场值，是我们真正要求的全局待求矢量。同样，经过以下步骤：bn+2＝split(xn)(s4.2.8)其中bn+2是新的右边向量b，代入式(s4.2.1)即可形成新的矩阵方程。s5：基于共形网格分界面上的棱边是公共边的特性，因此相邻子域棱边信息可以利用相应子域传输收发组合成全局棱边信息，每列循环迭代实现对该频率展开点下全局有限元系统krylov子空间的构建；s5.1：将全局有限元系统线性化为只和频率相关的线性系统，可得如下形式：其中s0表示频率展开点，k和s0b分别表示加入频率展开点后形成的新的维度为n的系数矩阵和右边向量。a2n表示的是与k系数矩阵相关的一个维度为2n的矩阵。σ是一个和频率有关的极小量，未知向量y与待求矢量x相关，只是临时的一列无需求出；s5.2：构建全局有限元系统krylov子空间(维度[n×q])具体实施流程如下：输入量：频率展开点s0，各子域不含频率项的si,ti,bi，全局不含频率项的s,t,b；(1)利用频率展开点s0，形成每个子域的系数矩阵ki和右边向量bi，利用区域分解迭代求解器，求解出xn+2(2)令p0＝gather(bi)；q0＝xn＝gather(xn+2)；(3)β＝||q0||(4)(5)for(j＝1；j＜q；j++)(6)(7)for(i＝1；i＜j+1；i++)(8)(9)(10)endfor(11)(12)if(hj+1,j≡0)(13)breakdown(14)(15)endfor(16)vn×q＝{q1,q2,q3,…,qq}其中，将*式展开，有后式代入前式，有：qj+1＝-k-1(tqj-1+2s0tqj)(s5.2.2)为使式(s5.2.2)适用于区域分解求解器求解，具体步骤如下：(1)bd＝split(tqj-1+2s0tqj)(2)(3)qj+1＝gather(xd)s6：利用构建出的krylov子空间对全局不含频率项的相应系数矩阵和右边向量降阶。循环扫频，加入频率项形成降阶后的有限元矩阵方程并采用直接求解器求解，并利用krylov子空间依次恢复出该频点上的全局棱边电场值。实施例1采用本发明所述方法对一个波导滤波器进行仿真分析，波导内部不填充介质，除两个端口外，其余外表面均为理想金属面。模型尺寸分别为：a＝19.05mm，b＝9.525mm，c＝2.86mm，d＝2.01mm，h＝18.9mm，l＝17.23mm，如图2所示。将整个波导滤波器模型划分为三个子域，每个子域在仿真计算机上分配一个进程，在左端施加te10模的入射波，右端加1/4最大波长厚度的完美匹配层(pml)，以达到完全吸收无反射的目的；设置剖分精度为1/15最小波长；扫频范围为10ghz～15ghz，设置200个频点，间距25mhz；分别用本发明方法与传统有限元模型降阶方法计算该模型在10ghz～15ghz的s参数即散射参数，其中传统方法具体采用迭代解法求解；如图3所示给出了采用本发明方法与传统方法计算出的该波导滤波器的s参数，其中，s21为正向传输系数，也就是增益，s11为输入反射系数，也就是输入回波损耗，不难看出，两者吻合较好，证明了本发明方法的正确性；从下表1中可以看出，传统方法与本发明方法的求解效率不同，在求解相同模型和未知量数目情况下，本发明仅需1658秒，而传统方法则需要30398秒，通过比较发现本发明比传统方法的求解效率高出18.5倍，证明了本发明方法相比于传统方法的高效性。表1：方法未知量展开点数v矩阵列数降阶时间(秒)总求解时间(秒)传统方法429672203038430398本发明方法4296722016411658当前第1页12