（电磁场与微波技术专业论文）有限元方法在二维散射问题中的应用.pdf

摘要 有限元是计算电磁学的主流方法之一，对复杂结构和非均匀介质问题有很 强的描述能力。在开域电磁散射问题中，有限元需要用虚拟边界截断无限大空 间，本文用完全匹配层作为主要的开域散射边界条件，对二维散射问题进行研 究。对比了吸收边界条件与共性匹配层的吸收性能，并对后者的使用条件及参 数设置进行了深入地探讨。 文中几个典型散射体的计算结果与其他文献的对照表明了这种方法的有 效性和准确性。 关键词有限元、共形匹配层、电磁散射 a b s t r a c t t h ef i n i t ee l e m e n tm e t h o d ( f e m ) i so n eo ft h em a j o rn u m e r i c a lm e t h o d si n c o m p u t a t i o n a le l e c t r o m a g n e t i c s d u et oi t sv e r s a t i l i t ya n df l e x i b i l i t y , t h ef e mi s c a p a b l eo fm o d e l i n gc o m p l i c a t e ds t r u c t u r e sa sw e l la si n h o m o g e n e o u sm a t e r i a l s w h e ns o l v i n ge l e c t r o m a g n e t i c ( e m ) s c a t t e r i n gp r o b l e mw i t ht h ef e m ，t h ei n f i n i t e r e g i o ne x t e r i o rt ot h es c a t t e r e rm u s tb et r u n c a t e db yaf i c t i t i o u sb o u n d a r y f o rm o s t a p p l i c a t i o n si nt h i sp a p e r , t h i st r u n c a t i o nb o t m d a r y - - t h ep e r f e c t l ym a t c h e d l a y e r ( p m l ) i sa d o p t e df o rt h ea n a l y s i so ft w od i m e n s i o n a ls c a t t e r i n gp r o b l e m w e c o m p a r et h ea b s o r b i n ga b i l i t yb e t w e e na b s o r b i n gb o u n d a r yc o n d i t i o n ( a b c ) a n d p e r f e c t l ym a t c h e dl a y e r ( p m l ) ，a n dw ea l s od i s c u s sd e e p l yo nt h ea p p l i c a t i o n c o n d i t i o na n dp a r a m e t e rs e r i n go fp m l s e v e r a le x a m p l e sa r ea n a l y z e d ，t h er e s u l t sc o m p a r e dt ot h o s eo b t a i n e db y o t h e ra u t h o r ss h o wt h ev a l i d i t ya n d a c c u r a c yo ft h i sm e t h o d k e yw o r d sf e m ，p m l ，e l e c t r o m a g n e t i c ss c a t t e r i n g 西北工业大学硕士学位论文第一章绪论 第一章绪论 1 1 电磁场有限元理论的发展 有限元建立了一条连续系统离散逼近的自然途径，也许这就是有限元方 法最重要的特性1 1 】。到目前为止，作为一种十分有效的数值算法，有限元法 已广泛的应用于各类工程技术领域，例如流体力学、空气动力学、结构应力、 应变分析、以及各种场变量温度、压力、电磁场的计算等等。本节主要 介绍有限元理论的一个分支，电磁学有限元理论的发展。 1 9 4 3 年，c o u r a n t 的论文中首次明确地出现了有限元的思想最小位能 原理及分片插值的离散形式。然而许多学者认为，有限元初期最重要的贡献 应来自于结构工程师l “，他们根据工作需要提出许多启发性的构想并付诸实 践。1 9 5 3 年，工程师们首次【3 j 用数字计算机实现了有限元刚度矩阵的求解。 1 9 6 0 年，有限元的名称 4 1 正式出现。 六十年代中期，数值分析学家认识到有限元思想的重要性，将偏微分方 程理论、泛函分析、逼近论等引入有限元的理论体系，建立了有限元方法的 数学基础。自此，有限元开始在工程计算中得到广泛应用。直至成为某些学 科数值算法的支柱。 1 9 6 9 年，p p s i l v e s t e r 有关波导模求解的文章是一个里程碑，标 志着电磁学有限元方法的出现。不久，有限元被用于求解电机磁场、静电场、 波导本征值、涡流场、散射与辐射场等等，涉及电磁学的各个领域。有关电 磁学有限元方法的论文数量逐年递增【“，六十年代末仅十几篇，七十年代末 达到1 0 0 篇左右，八十年代末已接近1 0 0 0 篇；到九十年代，有限元已经同矩 量法，时域有限差分法共同构成计算电磁学的三大主流算法。许多商业有限 元软件也包含了电磁场训算功能，例盘n m a r c 、a n s o f t 、a n s y s l 6 】等，专为天线， 微波器件设计服务的有限元软件也开始出现。 大量利技文献与商业化软件的出现表明，在电磁学相关领域，有限元理 沧已经趋向成熟，然而有限元方法t 分适用于闭域问题的求解，对于散射等 西北工业大学硕士学位论文第一章壹苦论 开域问题则较为困难，三十年来这个领域新的方法层出不穷，但是一直未发 展到接近工程应用的水平。显然，这是一个广阔而复杂，值得深入研究的领 域。 1 2 有限元理论在开域散射问题中的应用 有限元求解电磁散射问题的困难之处在于，有限元是一种区域性方法， 受计算机存储空间和运算速度的限制；而电磁散射是一个开域问题，有限元 方法不可能将离散域扩展到无限空间，这就需要一个特殊的边界条件开 域边界条件为有限元确定一个有限大的计算区域，将开域问题转化为闭域问 题。 有限元开域边界条件可以分为【7 j 局域边界条件( l o c a lb o u n d a r y c o n d i t i o n ) 与全域边界条件( g l o b a lb o u n d a r yc o n d i t i o n ) 。全域边界条件 准确地描述了散射体外部空间的电磁特性，是精确的边界条件，然而由于考 虑了边界上所有节点间的相互作用，全域边界条件离散化之后会形成一个稠 密矩阵，使有限元与全域边界条件的耦合矩阵失去稀疏性，不利于存储和求 解。相反，局域边界条件只考虑相邻节点间的相互作用，因此保持了系统矩 阵的稀疏性；不足之处是局域边界条件一般为精确辐射边界条件的某种近似， 精度和应用范围有限。 电磁散射问题的有限元理论与其它学科相比有许多共同之处，其个性部 分都直接或问接地与开域边界条件联系在一起，这一点是理解其发展趋势的 关键。所以，本节对这一领域研究进展的论述以开域边界条件为核心，以有 限元方法的实施步骤为框架，分为边界条件的处理、系统方程组的求解、物 理量与几何量的插值形式三部分。 1 边界条件的处理 有限元应用于散射问题之初，研究人员就 - 分清楚局域与全域边界条件 的缺陷，并且一商努力从各种途径去弥补它们，新的开域边界条件往往给研 究工作带来巨大的进展。 ( 1 ) 全域边界条件 电磁散射问题早期的开域边界条张二几乎都是全域边界条件，最早引入的 西北工业大学硕士学位论文第一章绪论 是边界积分方程( b o u n d a r yi n t e g r a le q u a t i o n ) ，先驱者是s i l v e s t e r 、 m c d o n a l d 等，他们在七十年代初的工作奠定了全域边界条件的基础，使有限 元具备了求解散射问题的能力。 单矩法( u n i m o m e n tm e t h o d ) i “9 j 是稍后出现的全域边界条件，在单矩法的 边界上，散射场展开为本征函数和的形式，待求变量是各个本征函数的展开 系数。由于这种方法的边界形状仅限于圆形( 2 d ) 、球形( 3 d ) ，对于狭长的散 射体效率不高，故未能得到推广。单矩法最新的进展是将边界形状扩展为椭 圆形( 2 d ) 。 再后来出现的有限元矩量法，有限元边界元1 1 3 j 等方法，在本质上都可 以归入有限元边界积分方程法【5 】或统称为有限元混合方法【1 0 】( h y b r i d m e t h o d ) 。七十年代由于计算机尚未普及，混合方法没能得到重视；至八十年 代中期，该方法逐渐得到关注，被广泛地应用于二维1 n , 1 2 , 1 3j 、三维电磁散射 计算以及天线辐射特性的计算，是公认的精确数值方法。 如前所述，全域边界条件离散化后产生稠密矩阵，破坏了系统矩阵的稀 疏性，因此限制了有限元求解电大尺寸散射问题的能力，在三维计算中这种 现象更为严重。现阶段主要的解决途径有两个：一是建立专门的耦合矩阵算 法，具体的研究进展见系统方程组求解部分：二是采用局域边界条件。 ( 2 ) 局域边界条件 吸收边界条件a b c ( a b s o r b i n gb o u n d a r yc o n d i t i o n ) 是局域边界条件的典 型代表，七十年代末出现，在有限元【“j 、时域有限差分6 1 等方法中得到广 泛应用。吸收边界条件模拟了散射体外部自由空间的电磁特性，使散射波能 够自由穿过边界，或者说被边界完全吸收。 吸收边界条件的建立方法是将局域算子应用于波动方程，或者使散射波 展开式的前几项在边界处被吸收。常用的是e n g q u i s t - m a j d a 和 b a y l i s s - t u r k e 型【”j 吸收边界条件，前者适用于平面边界，后者适用于圆 形、球形边界。二维的吸收边界条件已经有了完整的理论体系【1 8 ，l ，研究人 员对其较为熟悉，并在实践中广为采用。三维矢量吸收边界条件目前仍在发 展中，有许多问题尚待解决，是研究的前沿方向之一。由于提高吸收边界条 件的阶数或者采用“共形”技术能够减小边界与物体之间的距离，提高计算效 率，所以高阶和共形p “吸收边界条件也是研究热点。现有的高阶和共形吸收 西北工业大学硕士学位论文 第一章绪论 边界条件公式繁琐，计算迭代次数多，未能达到节约计算时间的目的，稳定 性与吸收效果亦不理想【2 l 】，总体上说是不成功的。高阶a b c 必须与高阶有限 元配合，为克服这一限制，s t u p f e l 等建立了数值吸收边界条件 瞄2 3 i n a b c ( n u m e r i c a la b s o r b i n gb o u n d a r yc o n d i t i o n ) 。n a b c 将边界点场强 的法向导数表示为其相邻节点场强的线性组合。 a 。”。) = c ，”，) j = 0 因为计算散射场需要知道的a 。”值，所以n a b c 带来的另一好处是便于 r c s 计算。相对于a b c ，n a b c 有一些独特的优点，但还不足以使其得到广 泛应用。 m e i ( m e a s u r e de q u a t i o no fi n v a r i a n c e ) 方法 2 q 的思路与n a b c 相似， 由m e i 等在1 9 9 2 年提出，最初用于截断有限差分网格，1 9 9 3 年开始被用 于有限元开域问题。在m e i 方法中，m e i 系数将边界某一点的场与其相邻节 点处的场联系起来： “。) ：羔叩缸，) j = 1 m e i 系数通过测度子( m e t r o n ，实际上是某种虚拟的源分布) 来确定。m e l 方法不会改变有限元矩阵的稀疏性，且具有能够贴近散射体放置的优点。但 是用m e i 方法不能直接得到a 。”的值，用差商得到的a 。“在精度上比”低一 阶。 迄今为止，m e i 方法的主要应用限于二维，公开文献中还没有见到有关其 稳定性的结论，m e i 方法的许多优点也需要在理论上进一步证实。与数值吸收 边界条件不同的是，尽管有关m e i 方法的疑问也许更多，但它的应用以及所引 发的研究热情却要广泛得多。 完全匹配层p m l ( p e r f e c t i ym a t c h e dl a y e r ) 是近年出现的，最出色的局 域边界条件。在理想状态下，p m l 对入射波的吸收与频率及入射角无关，不 会改变有限元系统矩阵的稀疏性，而且更好的吸收效果只需增) j h p m t 胸厚度就 可以实现。大量的研究成果 2 s - 3 1 】表明，p m t 。能够设置在距散射体很近的距离 i ：( 零点几个波氏) ，且吸收效果优于a b c 。 完全匹配层的概念由b e r e n g e r 于1 9 9 4 年首先以分离场变量的形式提出 4 西北工业大学硕士学位论文第一章绪论 【3 2 】，这种形式的p m l 对麦克斯韦方程做了较大的改动，仅适用于有限差分法。 有关这种形式p m l 的吸收特性及应用参见文献阿。“。 1 9 9 5 年，s a c k s 等提出了张量形式的p m l m j ，由于后者对有限元、有限 差分法来说只是一种各向异性的媒质，因此获得了广泛的应用，其中文献【鼍”】 在理论一k 对比了b e r e n g e r 的p m l 与张量形式的p m l ，表明二者在理论基础和 吸收效果上是等价的。张量形式的柱面、球面p m l 不久也出现了，其中柱面 p m l 在旋转体的散射计算中表现了极佳的吸收效果【“矧。最近的研究包括将 平面、柱面、球面p m l 组合使用mj ，使匹配层的外形更加灵活。 1 9 9 8 年，t e i x e i r a 等在理论上提出了共形匹配层，并将平面、柱 面、球面匹配层归结为共形匹配层的特殊形式，完善了张量匹配层的理论体 系。二维共形匹配层由h w a n g 等用f d t d 实现，对于椭圆柱，机翼等散射 体，反射误差( r e f l e c t i o ne r r o r ) 低于一4 0 d b 。l i u 等【3 9j 指出共形匹配层 要求其内部的等距曲面具有g 2 连续性，现有单元类型不能满足这一要求，文 中将主曲率用法曲率近似，用曲边六面体矢量单元实现了三维共形匹配层。 ( 3 ) 其它形式的边界条件 在静态、准静态电磁场问题中，最早的局域边界条件是把离主场区较远 的开域场强设为定值或零，例如在静电场问题中，标量电位以设为零；在静 磁场或准静态问题中，矢量磁位或标量磁位可以设为一个定值。显然，这种 方法的效率与精度都不高。气球法( b a l l o o n i n gt e c h n i q u e ) 、无限元 ( i n f i n i t ee l e m e n t ) 等也是首先应用在静态问题中的方法。文献对静态、 准静态的开域边界条件做了较全面的总结。场反馈法( f i e l df e e d b a c k f o r m u l a t i o n ) 、自适应吸收边界条件( a d a p t i v ea b s o r b i n gb o u n d a r y c o n d i t i o n ) 、透明吸收边界( t r a n s p a r e n ta b s o r b i n gb o u n d a r y ) 等则是近年 来发展的、针对散射问题的开域边界条件。出于种种原因，这些方法现阶段 未获推广；但它们各具特色且富有启发性，未来的发展潜力，以及是否会成 为主流方法尚不能妄下结论。 ( 4 ) 小结 局域边界条件和全域边界条件的绝大部分特点是互补的：局域边界条件 保持了有限元系统矩阵的稀疏性，但或多或少地引入了误差；全域边界条件 是精确的边界条件，可是破坏_ 有限元系统矩阵的稀疏性；局域边界条件必 西北工业大学硕士学位沧文 第一章绪论 须距散射体一定距离，保持一定的形状，不兔扩大了网格区域；全域边界条 件则可以紧贴散射体设置，形状也比较灵活。局域边界条件的发展方向是缩 小网格区域，降低反射误差；全域边界条件的研究则侧重于有限元系统矩阵 的处理和求解。但二者最根本的目标是相同的，即提高有限元求解电大尺寸 散射问题的能力。 2 系统方程组的求解 经过半个世纪的发展，有限元系统方程组的求解理论已经相当成熟，有 专门的系数矩阵存储方案及配套的解法。系统方程组的解法可以分为直接法 和迭代法两类 4 1 4 24 3 】：迭代法用有限步的迭代运算得到具有指定精度的近似 解，例如高斯一赛德尔迭代法；直接法则可以在预定的运算次数内求出方程 组的精确解，例如乔列斯基分解法。中小规模的方程组一般用直接法，大规 模、超大规模方程组宜于使用迭代法。对于不同的开域边界条件，散射问题 产生的系统方程组不同，相应地求解方法亦不同。 ( 1 ) 全域边界条件 全域边界条件离散化之后产生稠密矩阵。如果边界上有n 个未知量，则 生成n 阶方阵，直接求逆需要d ( 3 ) 次运算和o ( 2 ) 的存储量；对于电大尺 寸三维问题的全域边界条件，今天的超级计算机也难以处理。因此，这一领 域主要的研究方向是减少有限元一全域边界条件耦合矩阵的存储量和计算量 。 早期的耦合矩阵处理方法主要有内视法、外视法，二者本质上都是耦 合矩阵的分块解法；对于特殊的边界形状( 如圆形) ，由于边界积分是卷积， 可采用f f t 计算：显然，这几种方法并未从根本上解决问题。 快速多极子方法f m m ( f a s tm u t i p o l em e t h o d ) 陋删是很有前途的方法， 用这种方法可以将矩阵的存储量及矩阵与向量乘积的运算量降至o ( 15 ) 。如 果结合快速多级方法( m u l t i l e v e l ) ，上述两项指标可以降至o ( n l o g n l 。利 用小波的特性也可以将积分方程产生的稠密矩阵稀疏化，加快积分方程的求 解速度。 ( 2 ) 局域边界条件 局域边界条件保持了有限元系统矩阵的稀疏性，系统方程组的求解实际 j 二是稀疏矩阵的求解。对j i 大型方程组迭代法几乎是唯一可行的解法f 4 7 j 。 西北工业大学硕士学位论文 第一章绪论 现阶段流行的迭代法有c g ，b i c g 。q m r ，g m r e s 等，文献【艚】对这些方法做了 较全面的介绍。 迭代法的收敛速度主要取决于方程组系数矩阵的条件数，因此局域边界 条件系统方程组的求解主要体现在迭代法的预处理技术上。有限元、有限差 分等网格方法产生的系统方程组本身就是病态的，局域边界条件，特别是p m l 的引入进一步破坏了系统矩阵条件数，对预处理技术的要求也更高。根据作 者的经验，未加预处理的b i c g 在求解含p m l 的系统方程组时基本上不收敛。 因此，含p m l 的有限元系统方程组解法是一个前沿课题。最近，文献【4 9 】 研究了采用g m r e s 时的预处理技术，文献 5 0 】研究t b i c g 的预处理技术及并 行计算技术。 ( 3 ) 小结 局域边界条件和全域边界条件产生不同类型的系统矩阵，因此研究工作 的侧重点也不同：前者需要降低系统矩阵的条件数，提高收敛速度，主要研 究方向是预处理技术：后者则不仅要加快求解速度，而且要降低存储量，主 要研究方向是稠密矩阵的稀疏化方法。 3 物理量与几何量的插值形式 有限元用网格单元来描述对象的物理性质与几何参数，多样化的单元形 式是有限元的特色与优势之。电磁学有限元方法的单元类型可以分为节点 单元和矢量单元两大类。 节点单元是一种古老的单元形式，单元类型丰富，单元特性均有手册可 查，多数有限元软件使用的也是节点单元。节点单元描述标量场是很理想的， 因此适用于二维散射问题：用节点单元描述三维问题必须先将矢量场标量化， 再对标量方程进行离散。 一般认为 1 0 1 ，由于节点单元描述的物理量在单元间总是连续的，在三 维矢量场的计算中存在三个难以处理的问题： ( 1 ) 介质交界面上场的不连续性： ( 2 ) 伪解( s p u r i o u ss o lu t i o n ) 川的存在； ( 3 ) 导体、介质边角处场的奇异性。 这些难点几乎都可以通过引入矢量单元来克服。矢量单元f 5 2 】( v e c t o r e l e m e n t s ) 亦称棱边单元( e d g ee l e m e n t s ) ，它将物理量的自由度贼予单元的 7 西北工业大学硕士学位论文第一章绪论 棱边而不是节点。在相邻矢量单元的交界面上，切向场连续，法向场可以不 连续，因此难点( 1 ) 很自然地得到了解决。关于矢量单元能否消除伪解存在争 议。由于应用矢量单元后绝大部分算例确实没有伪解出现，而电磁散射问题 一般不涉及伪解( 伪解与静态场有关，静态场对远场没有贡献) ，所以难点( 2 ) 在本文可以搁嚣起来。 电磁场在导体、介质的边角处存在两类奇异性，一类是场矢量的方向发 生突变，另一类是场矢量的幅值变得无限大。如果不采用含有奇异性的基函 数，无论节点单元还是矢量单元都对第二类奇异性无能为力 8 2 ；但是矢量 单元能够描述第一类奇异性，即难点( 3 ) 可以部分解决。另外，矢量单元处理 普通的d i r i c h l e t 边界条件也比节点单元方便。 基于以上几个原因，大部分学者认为l ”j ，矢量单元在三维问题中的表现 优于节点单元。在电磁学领域，矢量单元的应用始于八十年代初期，九十年 代开始被研究人员广泛接受，成为三维计算电磁学的主流单元。矢量单元在 散射问题中的应用可参见文献【7 , 1 0 】。 矢量单元也存在某些缺点。首先，矢量单元对场变量的描述精度低于同 阶的节点单元：其次，矢量单元基函数问的正交性不如节点单元，系统矩阵 的条件数因此增大，影响迭代法的收敛速度；另外，矢量单元是一种新型的 单元，与绝大部分软件包的网格生成器、前后处理程序等不兼容，给实际应 用带来困难。 发展节点单元与矢量单元的混合单元则可以取长补短。混合单元的形式 有三种：一种是在矢量方程描述的区域采用棱边单元，在标量方程描述的区 域采用节点单元，两种区域的交界面建立相应的边界条件；另一种形式是在 结构或介质不连续的区域采用棱边单元，在媒质特性均匀的区域采用节点单 元；第三种形式是构造棱边一节点复合单元，法向使用节点基函数，切向使 用棱边基函数。混合单元的共同缺点是构建过程复杂，通用性差。 矢量单元和节点单元都在发展，最近的研究表明，通过建立适当的阻抗 边界条件，传统的节点单元电能够解决导体边角的奇异性问题。除了普通的 四面体、六面体矢量单元外，曲边高阶矢量单元也已经构造出来，这有助 _ i 二提高矢量单元的插值精度和汁算效率。目前高阶矢量单元的应用比较少见。 这说明其理论还不够成熟，网格生成器的研究r = 作血相对滞后。 西北工业大学硕士学位论文第一章绪论 有限元法弱形式的变分表达式仅要求单元问c o 连续。因此大部分研究人 员对单元拼接连续性没有足够的重视。我们在研究中发现，共形匹配层等较 复杂的局域边界条件要求单元间有高阶的拼接连续性，这一连续性一旦得到 保证。散射问题的计算精度会得到大幅度的提升。 2 小结 三维矢量场的插值问题是电磁学有限元方法所独有的，矢量单元的出现 解决了多年遗留的难题，促进了电磁学有限元方法的发展。从文献数量上看， 现阶段应用最广的是三角形( 2 d ) 、四面体单元( 3 d ) 等非结构单元，这类单 元结构简单，有限元刚度矩阵的元素能够解析计算，对复杂的几何结构有很 强的适应性。这种现象一方面说明工程应用需要简单、成熟的单元，另一方 面说明其它单元的特性尚不理想，网格生成的研究也相对滞后。 要实现算法的高效就必须进行分工，特殊的单元形式为特殊的目标服务。 因此，与系统方程组解法的研究一样，我们也需要适应不同开域边界条件的 单元形式，例如适于共形匹配层的、具有高阶连续性的单元。 4 总结 有限元的广泛应用，首先来自于工程分析与设计的需要。随着工程项目 的大型化、复杂化以及对性能指标要求的日益提高，越来越多的问题需要用 有限元等数值方法解决。有限元的数学基础已经比较稳固，研究方向集中在 具体问题的应用，以及同其它理论、方法的结合。数值方法受到计算机存储量 与运算速度的限制，对于散射问题的有限元方法，矛盾集中于如何处理开域 边界条件；其它方面的研究工作，如系统方程组的求解等都是围绕这个核心 展开的。 对于电磁散射的有限元方法，前沿的研究领域是： 1 、开域边界条件方面：全域边界条件的局域化方法：局域边界条件的误 差控制；共形局域边界条件的研究；新开域边界条件的研究。 2 、系统方程组的求解方面：稠密矩阵的稀疏化、快速解法；局域边界条 件下稀疏、病态有限元方程组的解法及其预处理技术的研究。 3 、单元插值形式方面：矢量单元插值精度和单元效率的提高：矢量一节 点混合单元的研究：满足高阶拼接连续性的单元；层次单元、样条单元的研 究。 西北工业大学硕士学位论文第一章绪论 4 、有限元与其它方法的结合，例如有限元与高频方法；有限元在具体问 题中的应用，如飞机进气道的散射特性研究等：时域有限元方法。 0 西北工业大学硕士学位论文第二章有限元的前期处理 第二章有限元的前期处理 2 1 引言 有限元法不仅适用于求解不规则的几何形状和边界条件，而且能计算解决 复杂的多介质问题；这是其它数值方法难于处理之处，而对于有限元法来说却 很简单。 有限元法求解的基本思想是把边值问题的求解转化为条件变分问题的求 解。具体过程是先将待求域剖分成有限个单元，然后通过分片插值把能量泛涵 求极值的问题转化为多元函数求极值的问题，进一步根据极值条件，得到多元 线性代数方程组，解此方程组就得到问题的近似解。 对于复杂的几何形f b 计算之前有限元法需要做大量的原始数据准备工 作。为避免手工输入出现差错，数据准备的工作应尽量减少。完成这项工作的 辅助程序主要有自动网格剖分、数据校对诊断、节点编号优化、绘图与显示等。 本文完成了自动网格剖分和绘图两部分。 2 2 平面有限元网格剖分概述 剖分作为有限元离散过程的第一步，其结果直接影响最终的计算精度及运 算时间。总体上说，剖分越细致，计算结果就越精确，但同时也需要更多的存 储容量和运算时问，处理好这两者之间的矛盾，减少人为干预，全面准确反映 物体的相关特性，这几点是网格自动剖分的主要任务。 剖分的基本原则是使得削分后的网格模型尽可能的逼近原始物体的形状。 从形成网格的区域来划分，有限元网格的剖分方法可分为整体法和分区法两大 类。 整体法是把整个结构当作一+ 个区域进行剖分，一次性形成整体网格，这种 方法适合规则、均匀材料组成的单联体结构。 分区法则足把整个结构分作若j f 个予区域，在各子区域进行剖分，最后将 西北工业大学硬士学位论文 第二章有限元的前期处理 各子区域拼接形成整体网格。这种方法适合不规则、材料非均匀的复杂结构。 从形成网格的坐标空间来说，剖分方法又可分为直接法和变换法两种。 根据边界的已知节点坐标，直接法可直接得出区域内网格节点的整体坐 标，因此能保涯网格在实际空间的形状、位置。但剖分难度大。 变换法则是先求规则区域网格节点的局部坐标，然后通过某种变换求出节 点的整体坐标。变换法的节点坐标容易计算，但难以照顾内部网格的实际位置。 本文采用的平面超单元变换剖分法分别属于上述的分区法和变换法。 2 3 平面超单元变换剖分法 以平面八节点等参数单元为例，本节简要介绍一下平面单元变换剖分法的 具体步骤： ( 1 ) 根据己知的结构信息，将有限元的分析对象划分成一定数目的子区 域，每个子区域称为个超单元。划分的原贝j 是每个超单元内部材料性质单一， 各边的边界能用二次曲线较好的描述。对应于与节点等参数单元，超单元的形 状也是八节点四边形，各边为二次曲线，其形状完全由八个节点确定，如图 2 1 。 ( 2 ) 确定超单元节点坐标，横竖两边( 非严格意义上) 需要生成的网格 节点数目等。如果生成不均匀网格，还要考虑剖分因子奎、以。对直线边， 只要确定改变所在的弧个端节点坐标就可以了，中间节点的坐标可不必输入， 通过自动插值得到。因为实际情况比较复杂，很难全面考虑，1 ，2 两步的部 分工作需手工完成。 ( 3 ) 将超单元变换成规则区域正方形，求出规则区域下各网格节点 的局部坐标。 ( 4 ) 利用参数变换求出超单元内对应节点的整体坐标，同时定义单元编 号，节点编号和单元信息。 西北工业大学硕士学位论文第二章有限元的前期处理 ( b ) 局部坐标 图2 - 1 超单元网格变换 若对局部坐标采用等间距剖分，则图2 - l ( b ) 中0 的局部坐标为 2 白2 1 + ( i - 1 ) 二n - 。唧1 ( 2 3 1 a ) ， 2 1 + ( j - i ) 嘉 ( 2 3 1 6 ) n 为只只、曼只边，即正方形f 方向要生成的节点数；m 为鼻只、置只边， 也即正方形玎方向要生成的节点数。剖分网格为八节点等参数单元时n 和m 只 能为奇数，等值线将正方形在古方向划分为n 一1 个单元，在，7 方向划分为m 一1 个单元。 西北工业大学硕士学位论文 第二章有艰元的前期处理 对应超单元中图2 1 ( a ) 中r 的整体坐标为 8 x = 。皤，叩) x 【。 ( 2 3 2 a ) 口= 1 3 y = n 。( 喜，r ) y 。 ( 2 3 2 b ) 口= l 善、印为正方形中0 对应网格节点的局部坐标，r 。、y 。为八个已知的超 单元节点坐标。为八节点超单元的形状函数，定义如下： 。= 三( 1 十甄x l + 呷吼沁幺+ 叩玩一1 ) 口= l ，2 ，3 ，4 虬：! ( i + 毒乞) ( j 玎z ) 口：5 ，7 口= 委( 1 + ，7 巩x 1 一亏2 ) a = 6 ，8 ( 2 - 3 3 ) 幺、为图2 - 1 ( b ) ，中八个节点的局部坐标，例如 善，= 1 ( 2 - 3 4 a ) 叩l = 1 ( 2 3 4 b ) 若对局部坐标采用不等间距剖分，则图2 - l ( b ) 中的r 的局部坐标为 2 掣 岛= 一l + 弓l 掣 p = l 2 7 7 。= 一1 + 鲁l 五； 西北工业大学硕士学位论文第五章稀疏矩阵的存储及求解 整刚矩阵主兀在一维矩阵中存储的序号为( 1 ，3 ，6 ，1 0 ，1 4 ，1 8 ) ，第i 行第j 列矩阵元素k 。( f j ) 在一维数组中的位置是：在第i 行主对角线元 素序号之前( f _ ，) 个。如第4 行第2 列元素k 。，在第4 行主对角线元素 足。之前( 4 2 ) = 2 个，即一维数组的笫1 0 2 = 8 个数组元素。 可见这种存储方式是以主对角线元素在一维数组中的位置为参考基准 的。实际编程中主元在一维矩阵中的位置存在数组岍】中，岫】称为主元指 示矩阵。建立主元指示矩阵i m 】后。k 】中每个元素都可通过【m 】来确定。 如图5 - 3 ：第i 行主元k i i 在一维存储方式下为黝私硒m 【丑：第i 行第j 列 元素k 。应是主元x 。向前移( 一j ) 个，即为黝铂嬲m 【f 】一( 一s ) l ；第i 行第 一个非零元素在该行第，一( n a l l 一m d 一1 d + l 列。 o 0 n a ( i 一1 ) 0zo n a ( i ) l 卜 n a ( i ) - r a ( i 一1 ) 图5 3 【k l 中元素的位置 形成主元指示矩阵m 】的核心问题是如何确定每行主元左边带内元素 的个数。对于未作任何处理的原始单元刚度矩阵，第i 行主元左边带内元素 的个数等于i 点与其相关节点号差的最大值；m a ) ( ( 一) ，j c ，j 为i 点 的相关节点。 5 3 整体刚度矩阵的集合 有限元法中所有单元刚度矩阵按“一一定规则合并，组成表示整个有限元 求解区域特征的整刚矩阵。在几个单- = z ，土、- + - - d i 处，场量的值对于所有包 含陔节点的单元都是相同的这是集合过程的基础。 西北工业大学硕士学位论文第五章稀疏矩阵的存储及求解 整刚集合的简要步骤如下： ( 1 ) 初始状态整刚一维数组置零： ( 2 ) 对全部单元逐个进行考察，计算出单元刚度矩阵，形成主元指示矩 阵脚】 ( 3 ) 将每个单刚矩阵的下三角部分集合到整刚矩阵一维数组中，实现 k 。】+ k 】j 区】 ( 5 31 ) 如上节所述，蟛应叠加到r m h r s m ( o o j ) 】中，这里下标i 、j 的 定义如图5 4 3 7 ( 1 1 5 1 4 8 2 圈5 - 4 单元节点位置的定义 在上图所示的八节点等参数中，( 1 ，2 8 ) 是局部坐标系下八个节点的 编号。假如节点6 在整体网格中的编号为1 0 0 ，节点7 在整体网格中的编号为 1 1 5 ；则单刚矩阵计算出后，他的第7 行第6 列元素叠加进一维正刚矩阵时i 为 1 1 5 ，j 为1 0 0 。 5 4l d f f 分解法 为便于说明有限元方程组在其系数矩阵一维变带宽存储下的求解过程， 本节简要介绍一t l d l r 分解法。 用正定对称矩阵d 】表示整刚矩阵k 】，将口】分解为下三角阵 l l 及上z 角阵眇】( c m u t 分解) l 】= 嗍旺，】 ( 5 4 1 ) 其中 西北工业大学硕士学位论文第五章稀疏矩阵的存储及求解 g l = m = 因为阻】是对称的，因此有 e r 可写成 0 o l l ，2 2 ，3 。如 ： | nl 。2 1 ”1 2 1 0 0 ( 5 - 4 2 ) ( 5 - 4 3 ) 1 4 - - k r = u l g l ( 5 - 4 4 ) 1 ，l 1 1 1 o l 。氓1 l 。! f l n k ，3 3 l ( 5 44 ) 或 阱= p 】i v , 】( 5 - 4 6 ) 则 阻】= 眇r 【d 1 瓯】 ( 5 4 7 ) 由于妙1 、 u 1 】都是上三角形矩阵，陋】、 u r 都是下三角形矩阵，则由 分解的唯一性有 叫= 妙r 【d 】 ( 5 - 4 8 a ) 陋，】= n 】 ( 5 4 8 b ) 比较眇】及眇】可见 u o = t ， ( f j ) ( f = ，) ( 5 - 4 1 2 ) 对于一般为变带宽带状、对称、正定的有限元整刚矩阵区1 ，l d l r 分解 具有以下特点： ( 1 ) 因为屯、k 在带外为零，同时带外订。为零，所以屯= 0 ，故带外 元素不必进行分解，也无须存储；但带内零元素要参与分解，必须存储，因 此正好适应了变带宽存储的特点； ( 2 ) 对其他带内元素分解时，右端项的化约和回代皆不参与运算。这一 特点对右端为多维向量时的求解有利，分解完成后重复进行化约和回代过程 即可。 5 5 整刚一维存储时的三上膨分解法 将l d l r 分解法用于整刚按一维变带宽存储的有限元方程求解时，还需要 对基本的l d l r 公式做一些改动，以提高求解和存储的效率。 ( 1 ) 医】的下三角带内元素按一维变带宽存储时，( 5 4 1 2 ) 式可修改为 不。 式中： k q 一t 。k l * 姬 i n t ” ，= 】2 ( 5 - 5 1 ) j m ，m 十1 ， o ( i j ) m = m i n ( m ，月，) ( 5 - 52 ) 所，、m ，是k 】中第i 行和第j 行最左边非零元素所占列号，如图5 5 所 西北工业大学硕士学位论文第五章稀疏矩阵的存储及求解 图5 5 变带宽存储中的整刚元素 图s 一5 中，o 为j 行最左边元素，。，为i 行最左边元素。，k 是待 修改为屯的元素，b 是尚未修改的主元。k 】矩阵的下三角部分毛之前的元 素己全部变换为陋】的相应元素( | cj ) ，位于0 之后的元素尚未修正。故计 算( 5 - 5 1 ) 式时，可自上而下从左到右逐行把k 】的下三角部分变换为下三 角阵团。这样，仍可用原来k 】的内存存储变换好的m 。 矩阵化约： 对如下形式的有限元整体方程 医】口 = ) 令 k 】= 陋p r l 陋r 则有 l i d 】“【三】7 口 = p 令 矿) = 【d 】_ l 时口j 则 弛弦】= 即 ( 5 53 ) ( 5 - 5 4 ) ( 5 - 55 ) ( 5 - 56 ) ( 5 57 ) 西北工业大学硕士学位论文 第五章稀疏矩阵的存储及求解 矩阵回代： 令 。 ，1 。，： o | 1 l | n ： | n l | n 2 此时眇】为上三角阵，则 展开后为 ， 一 l 。3 ，删 置 趁 b ： 玩 b 。 ( 5 5 8 ) y 。= 骂i 。 y ：= o 毛一，。- y 1 ) ，： ri - 1、 咒= f 骂一o y i ， ( 5 5 9 ) _ 呐 厂。一l、 y 。= l 或一o y 。l ，m ， 窿，】= d i - 1 e r 眇】伍) = 口) y l y 2 ： m ： y 。 由于化约后z 仍在原丸的位置，故可自下而上求得 x ( 5 51 0 ) ( 5 51 1 ) ( 5 51 2 ) r ，= y ，一虬h ( 5 51 3 ) 同样由于化约后屯仍在原j i 。的位置，故 所以有 t 1 1 j = l 。 ( 5 51 4 ) 如：；h viiiiiniiii儿 ：；， l ” 甜 o i 西北工业大学硕士学位论文第五章稀疏矩阵的存储及求解 一= y ，一乇坼t , ， i m t ，i = ”，拧一1 2 , 1 ( 5 51 5 ) i r = i + l 对上式的条件说明如下： ( 1 ) 规定，m 。是由于满足这一条件时，0 限制在k 行的带宽内取值； 否则k 处于带外，是零元素，在k 】中没有存储： ( 2 ) 上限b 与行号及半带宽d 有关 b = 7 + d ( 7 + d ) 这是因为陋1 按下三角存储，故上三角眇】的k 由f + l 至即，变为下 三角的i ，由i + 1 ( i 卅 ) 至b ，如图5 6 ： ( 3 ) m k 是第k 行最左边的列号，转置后变为k 的行号。 图5 - 6 循环上限b 的值 归纳以上分析，将l d l r 法用于求解有限元方程时具有以下特点： ( 1 ) 不改变k 】的带状结构，带外的零元素分解后仍为零元素； ( 2 ) 一维数组 删a r s 】中的任意一个元素仅在计算岛时使用一次， 以后不再使用，故( 或) 可存储在原的单元内，这样节省了大量的工 作单元，这是该方法一个很突出的优点； 圈5 - 7 计算时的有关元素 1_1 游 旧刘i p 州 一 h 一卜 一 一h 西北工业大学硕士学位论文第五章稀疏矩阵的存储及求解 ( 3 ) 在计算任意一个0 时仅与i 行半带宽谚范围内上部三角区域内 的元素有关，如图5 7 ： ( 4 ) 另外还有普通l d i ? 法本身所具有的一些优点，如速度快，精度高， 避免了开平方根的运算，程序简单等。 4 7 西北工业大学硕士学位论文第六章计算结果 第六章计算结果 为证明方法的有效性，下面给出几个典型二维散射体的计算结果。实际应 用p m l 方法过程中，可尽量缩小虚拟边界，以便减少节点数目。 1 金属圆柱 6 1 金属圆柱及椭圆柱 ( a ) t m 波 ( b ) t e 波 图6 1 金属圆柱双站r c s 西北工业大学硕士学位论文第六章计算结果 说明：，= 五，旯= 0 0 l m ，入射角庐= 0 0 ，计算结果与矩量法进行比较，两者 符合得较好。 2 金属椭圆柱 ( a ) 计算结果 ( b ) 参照结果 图6 2 金属椭圆柱的双站雷达散射截面( t m 波) 说明：金属椭圆柱面与虚拟边界之间剖分两层单元，参照结果摘自文献 厂 瓣吨乏 西北工业大学硕士学位论文 第六章计算结果 6 2 金属方柱 ( a ) t m 波 ( b ) t e 波 图6 3 金属方柱的双站r c s 说明：金属方柱，边k a = 2 2 ，五= o 0 1 m 入射角西= 0 t 1 西北工业大学硕士学位沦文 第六章计算结果 6 3 涂层柱体 1 ，涂层金属圆柱b = 2 5 9 ，r 。= 0 4 2 ，r 2 = 0 5 2 ) 雪 了 釜 d e i g e e ( a ) 计算结果 ( b ) 参照结果 图6 4 涂层金属圆柱的双站雷达散射截面 说明：涂层内分两层单元。参照结果摘自文献。 西北工业大学硕士学位论文第六章计算结果 2 涂层金属方柱 ( a ) 计算结果 ，如掣 ( b ) 参照结果 图6 5 涂层金属方柱双站雷达散射截面( t m 波) 说明：曲线a 即为无涂层导体方柱的t m 波散射截弼，参照结果摘自文献 西北工业大学硕士学位论文第六章计算结果 6 4 介质壳 无限长圆介质壳“= 4 ，r 。= 0 2 5 2 ，r 二= 0 3 2 ) d 6 9 r e e ( a ) 计算结果 簪d e g r e e s ) ( b ) 参照结粜 图6 6 无限长吲介质壳层的双站雷达散射截面 西北工业大学硕士学位论文 第六章计算结果 说明：参照结果摘自文