基于面片拼接的等几何分析方法求解波导本征值问题

王　峰，林　皋，刘　俊，2，李建波

(1.大连理工大学水利工程学院，辽宁大连116024; 2.南京水利科学研究院水文水资源与水利工程科学国家重点实验室，江苏南京210098)



基于面片拼接的等几何分析方法求解波导本征值问题

王峰1，林皋1，刘俊1，2，李建波1

(1.大连理工大学水利工程学院，辽宁大连116024; 2.南京水利科学研究院水文水资源与水利工程科学国家重点实验室，江苏南京210098)

摘要：基于NURBS的等几何分析方法集成了计算机辅助设计(CAD)和有限元方法的优点，CAD模型、网格划分和数值仿真均采用同样的几何描述.然而，由于单个NURBS曲面片拓扑的局限性，单片等几何分析方法难以处理介质分布不均匀以及截面形式复杂的多连通区域问题.本文基于面片拼接，将等几何分析方法用来求解此类问题的波导本征值.细分前后，NURBS曲面片拼接处的控制点和网格必须匹配.通过Galerkin法来离散波导本征值问题的Helmholtz控制方程，计算结果表明该方法具有自由度消耗小、精度高、收敛速度快等优点.

关键词：波导本征值;等几何分析; NURBS;面片拼接;介质分布不均匀;多连通区域

1　引言

波导［1］是一种在微波或可见光波段中传输电磁波的装置，常用于无线电通讯、雷达、导航等无线电领域，其常用的截面形状有矩形波导、圆形波导、椭圆波导和环形波导［2］.按照电磁波在传输方向上有没有电场和磁场分量，可以分为横磁波(Transverse Magnetic，TM)、横电波(Transverse Electric，TE)、横电磁波(TEM).其中波导本征值问题是微波理论中最关键问题之一，因为它不仅与不同模式的电磁波传输特性有关，而且还是许多微波部件分析优化的基础.

求解波导本征值问题的方法主要有解析法和数值解法两类.解析法能求解的实例比较少，且截面形状较为简单正规，介质均匀规则.除了均匀波导外，不均匀波导元件目前已经大量的应用于微波系统［3］中，如可变移相器、抗匹配器等.

对于介质分布不均匀和截面复杂的波导，只能通过数值解法来求解.随着计算机科学技术的发展，出现了许多求解波导本征值问题的数值解法，如有限元方法［4～7］，边界元方法［8］，无网格法［9，10］，比例边界有限元方法［11，12］.

等几何分析方法(Isogeometric Analysis，简称IGA)是Hughes［13］于2005年提出来的，其将NURBS(Non-U-niform Rational B-spline)基函数引入到等参有限元中.该方法可以有效衔接计算机辅助设计(Computer Aided Design，CAD)和有限元方法(Finite Element Method，FEM)，将分析计算建立在精确的几何模型之上，消除了计算模型与几何模型之间的非一致性问题.它同时兼备了FEM和无网格方法的优点，其基本思想是，几何模型和网格计算模型采用相同的基函数空间.该方法具有以下独特的优势:(1)避免了FEM分段多项式逼近问题域所产生的误差;(2)由于NURBS的几何精确特性，即使是稀疏的参数网格也保留了模型的几何信息，同时细分过程中无需访问原CAD模型，这便于实现网格的自适应保形细分，节省了工作量.

目前该方法已经被用于电磁学问题［14～16］，张勇［17］将其与比例边界有限元(Scaled Boundary Finite Element Method，SBFEM)耦合来求解波导本征值问题.然而由于NURBS曲面片是定义在矩形参数域中的，受制于拓扑的局限性，单片等几何分析方法尚且难于处理介质分布不均匀以及截面形式复杂的波导本征值问题.由于NURBS曲面片的张量积特性，可以通过多面片拼接［18］的等几何分析方法来处理上述问题，即对NURBS曲面片进行“加法”操作，这充分借鉴了子结构法的基本思想［19］.由于等几何分析方法采用等参元的思想，电场或磁场纵向分量用NURBS基函数来构造，Helmholtz控制方程的弱形式通过Galerkin变分原理来实现，最后施加不同的边界条件，便可得横磁波(TM)和横电波(TE)的等几何分析方程，本文通过三个算例来验证其有效性.

2　等几何分析基本思想

2.1NURBS基函数

NURBS基函数是通过B样条基函数的有理化来实现的［20］，令Ξ= {ξ0，ξ1，…，ξn + p +1}是一个参数空间中的单调不减的坐标序列，即ξi≤ξi +1，i = 0，1，…，n + p.其中ξi是节点，Ξ为节点矢量，n +1和p分别是B样条基函数的个数和次数.对应于节点矢量Ξ的各阶基函数可由de Boor-Cox公式递归得到，即

当节点矢量Ξ端点的重复度等于p +1时，这样的节点矢量是开放型的，开放型节点矢量是CAD标准的，其定义的B样条基函数在端点处等于1，通常把节点矢量定义在标准区间［0，1］中.在等几何分析中，单元的细分有三种方法，即节点插入的h细分、基函数升阶的p细分和结合两者而成的k细分方法.图1给出了节点矢量为Ξ= {0，0，0，0.2，0.4，0.6，0.8，1，1，1}的基函数图形.

赋给每个B样条基函数Ni，p(ξ)权值wi(0＜wi≤1)，并进行加权平均，就得到NURBS基函数

若权值相等，则NURBS基函数退化为B样条基函数.二维NURBS基函数可通过一维NURBS基函数张量积得到，即

给出控制点Pi，便可构造NURBS曲线，即

式中: Pi是NURBS控制点.图2给出了一NURBS曲线，其采用的基函数如图1所示.

2.2波导本征值问题的变分方程

在等几何分析中，NURBS基函数将节点矢量所张成的参数域^Ω映射为控制点所在的物理域Ω，即

对于NURBS曲面也有类似的表达，且为

其中: Nk为NURBS基函数，Pk为控制点，NC为控制点个数，也是基函数的个数，x =(x，y)为笛卡尔坐标，ξ=(ξ，η)为参数坐标.等几何分析采用等参元的思想，即场变量uh(x)与几何形状x的近似采用相同的NURBS基函数，则任意场变量可近似表示为

式中: uk是与控制点Pk对应的控制点变量.

问题域Ω的网格可由节点矢量和参数域中的节点区间映射生成.设节点矢量为{η0，η1，…，ηm}，从两个方向节点矢量中去除重复节点可得子矢量，定义参数域中的节点区间Qi，j=通过式(6)可映射成物理域上的一个剖分.可以看出，网格的生成和细化均可以自动实现，且自动生成的网格能够准确描述计算物理域，不会引入离散误差，这实现了几何设计和数值分析的无缝连接.

波导内电磁场量满足的控制微分方程为

式中:▽2为Laplace算子;对于TM波，u表示为纵向电场分量，即u = Ez;对于TE波，u表示为纵向磁场分量，即u = Hz.kc为截止波数，为频率，μ为磁导率，ε为介电常数，β为相位常数.

不同波型满足的边界条件为

式中:∂Ω为问题域Ω的边界线，n为外法线方向.

由Galerkin变分原理，上述问题可推得如下代数方程

其中:

由于NURBS基函数具有局部支撑性和非负性，只有少量NURBS基函数在边界上非零.对于TM波，边界上为零的基函数为，不为零的基函数为，则有

由NURBS基函数的单位分解特性可知

2.3NURBS曲面片的拼接

在实际的波导装置中，波导横截面往往需要多个NURBS曲面片来共同描述，例如介质分布不均匀和截面形式复杂的波导.然而由于NURBS曲面片的分片连续性及有理形式，NURBS曲面片实现光滑拼接出现了瓶颈，这里我们只研究交接处C0连续下的NURBS曲面片的拼接问题.

多片等几何分析的基本思想为:求出每块NURBS曲面片的系数矩阵，然后再通过连通数组(connectivity arrays)组装成波导横截面的总体系数矩阵.这里我们以两个B样条曲面片拼接为例来说明，如图3所示，下标f代表交接处的控制点，下标n代表面片其他处的控制点.

其中

如图4所示，如果对面片2进行节点插入的h细分方法，我们可以得到一组新的控制点.

设片1和片2的控制点变量为

交接处解满足C0连续

3　数值算例

3.1非均匀介质的矩形波导

本节应用多片等几何分析方法求解波导本征值问题，为了面片拼接的实施，这里我们选取齐次控制点

同时为了验证方法的有效性，我们选取非均匀矩形波导、均匀弯曲型L形金属波导和均匀十字形金属波导进行求解，定义数值解的相对误差为

为了验证多片等几何分析方法的有效性，本文选取如图5所示的波导模型，相关材料参数为ε1=1.0，ε2=1.5，μ1=μ2=μ0.表1给出了通过多片等几何分析方法和有限元方法计算的前三阶模态的截止波数，通过与解析解对比可以看出，在采用更少自由度的情况下，等几何分析方法可以求得更准确的TE、TM波截止波数.图6给出了TE11模态的相对误差随自由度变化的收敛图，很明显，随着自由度的增加，等几何分析方法和有限元方法均能收敛于解析解，但是等几何分析方法误差小、收敛速度快.同时可以看出，等几何分析方法在1039个自由度的情况下，TE11模态的相对误差是0.14%，而有限元方法在2197个自由度的情况下，其相对误差为1.39%.

表1　截止波数(kca)对比

表2　计算自由度、计算内存占用和计算时间比较

为了与FEM在计算内存占用和计算时间作对比，本文选用同样的400个网格剖分，其对应的计算自由度数、计算内存占用和计算时间见表2.在此网格剖分的情况下，IGA计算得到的TE11模态的相对误差是0.16%，而FEM计算得到的TE11模态的相对误差是1.22%.这里很明显可以看出，IGA中的NURBS等参元之间可以共享更多的自由度，从而节省了自由度消耗，降低了计算所用的空间存储.

3.2均匀弯曲型L形金属波导

考虑弯曲型金属波导，其界面形式如图7所示，它通过两个NURBS曲面片来构造其图形，其齐次控制点见表3.表4给出了等几何分析方法和有限元方法计算的TE波前五阶模态的截止波数，通过与参考解对比可以看出，在自由度差不多的情况下，等几何分析方法求得的TE波截止波数误差远小于有限元方法得到的误差，甚至小三个数量级，这充分说明了等几何分析方法可以通过更少的自由度来求得TE波更准确的截止波数.

表3　两个曲面片的齐次控制点

表4　截止波数对比

为了与FEM在计算内存占用和计算时间作对比，本文选用同样的网格剖分，其对应的计算自由度数、计算内存占用和计算时间见表5.由于NURBS等参元之间可以共享更多的自由度，因此在同样数量单元的情况下，IGA所花费的计算自由度会少于FEM.图8给出了在第2个网格剖分的情况下前5阶截止波数相对误差变化趋势图，可以看出FEM相对误差要远远大于IGA的计算相对误差，甚至大两到三个量级.可见，IGA与传统的FEM相比，具有消耗自由度少、精度高、收敛速度快等优点.

表5　计算自由度、计算内存占用和计算时间比较

3.3均匀十字形金属波导

十字形金属波导［5］的截面形式如图9所示，这里可以通过五个NURBS曲面片拼接构造.为了与文献［5］取得的结果作比较，只给出了通过等几何分析方法求得的TE波第一、第三阶模态的截止波数，见表6.

通过与参考解对比可以看出，在自由度差不多的情况下，等几何分析方法求得的TE波截止波数与有限元方法计算得到的TE波截止波数相近，这验证多片等几何分析方法的有效性.

表6　TE波第一、第三阶模态的截止波数对比

4　结论

本文通过拼接将等几何分析方法由单片扩展到多片，并推导了波导本征值问题的等几何分析离散方程.通过对非均匀介质和截面形式复杂的波导求解，可以看出等几何分析方法具有计算自由度少、收敛速度快、精度高等优点，而且还具有非常好的细分保形特性.但是面片拼接需要定义多个曲面片，这不利于解决多孔结构，因此面片拼接和面片裁剪结合是未来等几何分析方法走向工程应用的关键.

参考文献

［1］俎栋林.电动力学［M］.北京:清华大学出版社，2006.162-246.Zu Dong-lin.Electrodynamics［M］.Beijing，China: Tsinghua University Press，2006.162 -246.(in Chinese)

［2］冯剑，张贵新，刘亮，等.一种环形波导微波等离子体装置［J］.高电压技术，2009，35(1):48 -53.Feng Jian，Zhang Gui-xin，Liu Liang，et al.A microwave plasma system with ring waveguide［J］.High Voltage Engineering，2009，35(1):48 -53.(in Chinese)

［3］张克潜，李德杰.微波与光电子学中的电磁理论［M］.北京:电子工业出版社，2001.179 -202.Zhang Ke-qian，Li De-jie.Electromagnetic Theory for Microwaves and Optoelectronics［M］.Beijing，China: Publishing House of Electronics Industry，2001.179 -202.(in Chinese)

［4］Schiff Bernard，Yosibash Zohar.Eigenvalues for waveguides containing re-entrant corners by a finite-element method with superelements［J］.IEEE Transactions on Microwave Theory and Techniques，2000，48(2):214 -220.

［5］徐善驾.波导本征值问题的有限元分析［J］.电子学通讯，1982，4(4):222 -234.Xu Shan-jia.The finite-element analysis of waveguide eigenvalue problem［J］.Journal of Electronics，1982，4(4): 222 -234.(in Chinese)

［6］李融林，倪光正，俞集辉.B样条有限元法解波导本征值问题［J］.电子学报，1997，25(3):5 -9.Li Rong-lin，Ni Guang-zheng，Yu Ji-hui.B-spline finite element solution of waveguide eigenvalue problems［J］.Acta Electronica Sinica，1997，25(3):5 -9.(in Chinese)

［7］李建兵.基于有限元法的填充非均匀介质脊波导传输特性研究［D］.兰州:兰州交通大学，2012.1 -61.Li Jian-bing.Research on transmission characteristics of inhomogeneous dielectric loaded ridge waveguide by the finite element method［D］.Lanzhou，China: Lanzhou Jiaotong University，2012.1 -61.(in Chinese)

［8］占腊民.两类电磁场本征值问题的研究［D］.武汉:华中科技大学，2004.79 -92.Zhan La-min.Study on two kinds of electromagnetic field eigenvalue problems［D］.Wuhan，China: Huazhong University of Science and Technology，2004.79 -92.(in Chinese)

［9］Ooi B L，Zhao G.Element-free method for the analysis of partially-filled dielectric waveguides［J］.Journal of Electromagnetic Waves and Applications，2007，21(2):189 -198.

［10］张淮清，俞集辉.波导本征问题分析的径向基函数方法［J］.电子学报，2008，36(12):2433 -2438.Zhang Huai-qing，Yu Ji-hui.Radial basis function method for the eigen analysis of waveguide［J］.Acta Electronica Sinica，2008，36(12):2433 -2438.(in Chinese)

［11］Lin G，Liu J，Li J B，et al.Scaled boundary finite element approach for waveguide eigenvalue problem［J］.IET Microwaves，Antennas＆Propagation，2011，5(12):1508 -1515.

［12］林皋，刘俊.波导本征问题分析的比例边界有限元方法［J］.计算力学学报，2013，30(1):1 -9.Lin Gao，Liu Jun.Scaled boundary finite element method for the eigen analysis of waveguide［J］.Chinese Journal of Computational Mechanics，2013，30(1):1 -9.(in Chinese)

［13］Hughes T J R，Cottrell J A，Bazilevs Y.Isogeometric analysis: CAD，finite elements，NURBS，exact geometry and mesh refinement［J］.Computer Methods in Applied Mechanics and Engineering，2005，194(39):4135 -4195.

［14］Buffa A，Sangalli G，Vázquez R.Isogeometric analysis in electromagnetics: B-splines approximation［J］.Computer Methods in Applied Mechanics and Engineering，2010，199(17):1143 -1152.

［15］Zhang Y，Lin G，Hu Z Q，et al.Isogeometric analysis for elliptical waveguide eigenvalue problems［J］.Journal of Central South University，2013，20(1):105 -113.

［16］张勇，林皋，刘俊，等.波导本征问题的等几何分析方法［J］.应用力学学报，2012，29(2):113 -119.Zhang Yong，Lin Gao，Liu Jun，et al.The isogeometric analysis for eigen problem of waveguide［J］.Chinese Journal of Applied Mechanics，2012，29(2):113 -119.(in Chinese)

［17］张勇，林皋，胡志强.比例边界等几何分析方法Ⅰ:波导本征问题［J］.力学学报，2012，44(2):382 -392.Zhang Yong，Lin Gao，Hu Zhi-qiang.Scaled boundary isogeometric analysis and its application I: eigenvalue problem of waveguide［J］.Chinese Journal of Theoretical and Applied Mechanics，2012，44(2):382 -392.(in Chinese)

［18］Cottrell J A，Hughes T J R，Bazilevs Y.Isogeometric Analysis: Toward Integration of CAD and FEA［M］.New York: Wiley，2009.87 -92.

［19］张亚辉，林家浩.结构动力学基础［M］.大连:大连理工大学出版社，2007.217 -227.Zhang Ya-hui，Lin Jia-hao.Fundamentals of Structural Dynamics［M］.Dalian，China: Dalian University of Technology Press，2007.217 -227.(in Chinese)

［20］Piegl L，Tiller W.The NURBS Book［M］.Berlin: Springer Verlag，1997.47 -139.

王峰男，1987年5月出生，山东莱阳人，大连理工大学水利工程学院博士研究生.主要从事无网格法及等几何分析的研究及其应用于力学和电磁学分析.

E-mail: wangfengdut@gmail.com

林皋男，1929年1月出生，江西丰城人，大连理工大学水利工程学院教授，博士研究生导师，中国科学院院士.主要从事大坝、核电结构、地下结构地震响应分析与安全评价，以及等几何分析在力学边值问题中的应用等方面的研究.

E-mail: gaolin@dlut.edu.cn

Isogeometric Analysis for the Waveguide Eigenvalue Problem Based on the NURBS Patch Splicing Technique

WANG Feng1，LIN Gao1，LIU Jun1，2，LI Jian-bo1
(1.School of Hydraulic Engineering，Dalian University of Technology，Dalian，Liaoning 116024，China; 2.State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering，Nanjing Hydraulic Research Institute，Nanjing，Jiangsu 210098，China)

Abstract:Isogeometric analysis(IGA)based on the non-uniform rational B-splines(NURBS)has integrated the advantages of Computer Aided Design(CAD)and the finite element method(FEM).The main feature of IGA is the usage of one common geometry representation for creating CAD models，for meshing，and for numerical simulation.However，IGA based on one single NURBS patch is difficult to deal with the inhomogeneous mediums and complex multiply connected domains because of the limitation of the NURBS patch topology.In this paper，IGA based on the patch splicing is used to solve the waveguide eigenvalue problem of this kind.The control points and meshes of different patches must coincide on the interface，even after refinement.The Helmholtz governing equation can be discreted using the Galerkin procedure.Numerical examples are presented to show that IGA possesses the advantages of better convergence on a per-degree-of-freedom and high accuracy.

Key words:waveguide eigenvalue; isogeometric analysis; NURBS; patch splicing; inhomogeneous mediums; multiply connected domains

作者简介

基金项目:国家自然科学基金重点项目(No.51138001);国家自然科学基金委创新研究群体基金(No.51121005);国家重大科技专项(No.2011ZX06002 -10，No.2013ZX06002001-09);水文水资源与水利工程科学国家重点实验室开放研究基金(No.2012491611);上海交通大学海洋工程国家重点实验室开放基金(No.1202);中国博士后科学基金(No.2013M530919，No.2014T70251)

收稿日期:2014-05-12;修回日期: 2014-07-30;责任编辑:李勇锋

DOI:电子学报URL：http: / /www.ejournal.org.cn10.3969/j.issn.0372-2112.2016.01.029

中图分类号：TN814

文献标识码：A

文章编号：0372-2112(2016)01-0200-06