（车辆工程专业论文）基于参数曲面的边界面法研究.pdfVIP

摘要 本文介绍了边界元法的一种新的实现方法，称为边界面法。在传统的边界元 法中，网格单元不仅用来进行边界积分和变量插值，而且用来近似几何体。当离 散网格较稀疏时，会引起较大几何误差，因而影响计算精度。在边界面法中，边 界积分和场变量插值都是在以边界表征的实体边界曲面的参数空间里进行。在边 界积分过程中，积分点的几何数据，如物理坐标、雅可比、外法向量都是直接由 曲面算得，而不是通过单元插值近似，从而避免了几何误差。 文中首先基于格林公式简洁地推导了满足拉普拉斯控制方程的位势问题的正 则边界积分方程。然后讨论了在曲面参数空间利用曲面单元实现边界积分和变量 插值方案。曲面单元定义在曲面的二维参数空间，这与边界元法中定义在三维物 理空间的曲面单元不同。为了有效地分析薄型和细长结构，在参数空间里基于曲 面单元细分方法开发了自适应近奇异和奇异积分技术。单元细分直接在曲面的参 数空间里进行，简单而有效。根据推进波前法实现的一般流程，对该方法进行了 改进并开发了自适应曲面网格生成技术。在网格生成过程中，利用自适应四叉树 结构控制网格尺寸相应地生成曲面单元。 边界面法的实现是直接基于边界表征的c a d 模型，因此该方法具有与c a d 软 件无缝连接的天然优势，能使几何设计与工程分析集成于统一的框架。我们利用 c + + 语言开发了边界面法与u g 的接口，具有复杂几何的数值实例说明了这两者的 集成是可行的，这将为实现自动化分析迈出了重要的一步。 具有不同几何、边界条件类型的大量三维位势问题的数值实例表明，边界面 法具有良好的收敛性和可以获得较精确的数值结果。通过与传统边界元法比较表 明，该方法不仅比边界元法具有更高的精度，而且数值结果对网格密度的敏感性 较小。另外，边界面法可以方便有效地分析具有细小特征和薄型结构的三维位势 问题。 关键词：边界元法；参数曲面；曲面单元；边界面法；c a d 模型 a b s t r a c t t h i sw o r kp r e s e n t san e wi m p l e m e n t a t i o no ft h eb o u n d a r ye l e m e n tm e t h o d ( b e m ) ， h e r ec a l l e dt h eb o u n d a r yf a c em e t h o d ( b f m ) t h ec o n v e n t i o n a lb e m u s e st h es t a n d a r d e l e m e n t sf o rb o u n d a r yi n t e g r a t i o na n da p p r o x i m a t i o no ft h eg e o m e t r y t h eg e o m e t r i c e r r o r si n t r o d u c e db yc o a r s em e s hc a nl e a dt oa c c u r a c yp r o b l e m s h o w e v e r , i nt h eb f m b o t hb o u n d a r yi n t e g r a t i o na n dv a r i a b l ea p p r o x i m a t i o n a r ep e r f o r m e di nt h e2 - d p a r a m e t r i cs p a c eo ft h eb o u n d a r ys u r f a c e o fs o l i dr e p r e s e n t e db yt h eb o u n d a r y r e p r e s e n t a t i o nd a t a ( b r e p ) f o rb o u n d a r yi n t e g r a t i o n ，t h eg e o m e t r i cd a t aa t g a u s s i a n i n t e g r a t i o np o i n t s ，s u c ha st h ep h y s i c a lc o o r d i n a t e s ，t h ej a c o b i a n sa n dt h eo u tn o r m a l s ， a r ec a l c u l a t e dd i r e c t l yf r o mt h ef a c e sr a t h e rt h a nf r o me l e m e n t s ，a n dt h u sn og e o m e t r i c e r r o rw i l lb ei n t r o d u c e d i nt h eb e g i n n i n go ft h i sp a p e r ，w ed e d u c e ds e l f - r e g u l a rb o u n d a r yi n t e g r a le q u a t i o n o ft h ep o t e n t i a lp r o b l e m sg o v e r n e db yl a p l a c e se q u a t i o ns i m p l yb a s e do nt h eg r e e n f o m u l a t h e nt h es c h e m e sf o rt h eb o u n d a r yi n t e g r a t i o na n dv a r i a b l ea p p r o x i m a t i o n b a s e do nt h es u r f a c ee l e m e n ti np a r a m e t r i cs p a c ea r ed i s c u s s e d e a c hs u r f a c ee l e m e n t i sd e f i n e db yas e to fn o d e sg i v e ni n2 - dp a r a m e t r i cs p a c eo ft h es u r f a c ew h i c hi s d i f f e r e n t 行o mt h es u r f a c ee l e m e n td e f i n e di n3 dp h y s i c a ls p a c ea su s e di nt h eb e m i no r d e rt od e a lw i t ht h i na n ds l e n d e rs t r u c t u r e se f f e c t i v e l y , a na d a p t i v ei n t e g r a t i o n s c h e m ef o rn e a r l ys i n g u l a ri n t e g r a l sa n ds i n g u l a ri n t e g r a l sh a sb e e nd e v e l o p e db a s e d o n2 - ds u r f a c ee l e m e n ts u b d i v i s i o n t h ep r o c e s so fs u b d i v i s i o ni sp e r f o r m e dd i r e c t l y i n p a r a m e t r i es p a c em a k i n g i t si m p l e m e n t a t i o ns i m p l ea n de f f e c t i v e w e h a v e i m p r o v e dt h ei m p l e m e n t a t i o no ft h ea d v a n c i n gf r o n tm e t h o d ( a f m ) a c c o r d i n gi t s g e n e r a lp r i n c i p l e a n da na d a p t i v em e t h o df o rg e n e r a t i n gs u r f a c ee l e m e n t sh a sa l s o b e e nd e v e l o p e d s u r f a c ee l e m e n t sa r eg e n e r a t e da l o n gw i t ha na d a p t i v eq u a d t r e e p r o c e d u r eu s e dt oc o n t r o lt h ee l e m e n ts i z e t h eb f mh a sr e a lp o t e n t i a lt os e a m l e s s l yi n t e r a c tw i t hc a ds o f t w a r e ，i n t e g r a t i n g e a s i l yg e o m e t r i cd e s i g na n de n g i n e e r i n ga n a l y s i si n t oac o m p l e t e l yu n i f i e df r a m e w o r k ， b e c a u s ei t si m p l e m e n t a t i o nc a nb ed i r e c t l yb a s e do nac a d m o d e lt h r o u g hi t sb r e p d a t a w eh a v ed e v e l o p e da ni n t e r f a c eb e t w e e nb f ma n du g - n xw i t hc + + c o d e n u m e r i c a le x a m p l e si n v o l v i n gc o m p l i c a t e dg e o m e t r yh a v ed e m o n s t r a t e dt h a t t h e i n t e g r a t i o no fb f ma n du g - n x i ss u c c e s s f u l ，w h i c hm a yp r o v i d ea ni m p o r t a n ts t e p t o w a r da u t o m a t i cs i m u l a t i o n t h eb f mh a sb e e nv e r i f i e dt h r o u g han u m b e ro fn u m e r i c a le x a m p l e sw i t h d i f f e r e n tg e o m e t r i e sa n db o u n d a r yc o n d i t i o nt y p e i ti sd e m o n s t r a t e dt h a to u rm e t h o d h a sh i g hc o n v e r g e n c er a t ea n dy i e l d sv e r ya c c u r a t en u m e r i c a lr e s u l t s c o m p a r e dw i t h t h ec o n v e n t i o n a lb e m ，t h eb f mc a nn o to n l yp r o v i d em o r ea c c u r a t er e s u l t st h a nt h e b e m ，b u ta l s oi sl e s ss e n s i t i v et ot h ec o a r s e n e s so f t h em e s h w h a t sm o r e ，o u rm e t h o d c a nd e a lw i t ht h es t r u c t u r e sw i t hf e a t u r e si ns m a l ls i z eo rt h i ns t r u c t u r e st os i m u l a t e 3 一dp o t e n t i a lp r o b l e m se f f i c i e n t l y k e yw o r d s ：b e m ；p a r a m e t r i cs u r f a c e ；s u r f a c ee l e m e n t ；b o u n d a r yf a c em e t h o d ；c a d m o d e l i n g l v 1 1研究背景 第1 章绪论 绝大多数工程实际问题是不可能得到解析结果，只能依靠数值计算方法得到 所要求精度的各种结果。因此，以数值计算方法为基础的计算机辅助工程( c a e ) 软件被广泛用来分析实际工程问题。利用c a e 软件成功地推动工程的发展和产品 研发，面临两大难点【l ，2 】：l 、如何有效地将任意复杂形状的实体离散成有限个单 元，为后续相应的数值计算方法建立分析模型，比如有限元法( f e m ) 【3 】、边界元法 ( b e m ) 阶】。2 、在实际问题中不可避免存在大规模的数值计算。 当今有限元法发展比较成熟，相应地产生了以该方法为基础的大量商业软件， 如a n s y s 、n a s t r a n 和a b a q u s 等。这些软件在分析实际工程问题和研发产品 中发挥了重要的作用。但是，由于有限元法需要将整个求解区域离散成有限个单 元，因此将会导致一个很大的代数方程组，特别是求解三维( 3 d ) 问题。对三维问 题，需要将用来分析的c a d 几何模型离散成若干个体单元网格。网格单元用来能 量积分和变量插值，还用来插值逼近原始的c a d 几何。对于复杂的实体并且含有 细小特征，离散成为具有形状良好的实体单元往往是比较困难的。此外，有限元 法【2 】的实现是基于所求问题控制方程和边界条件的等效积分“弱 形式，其试函数 ( 一般为插值函数) 要求具有一阶连续性。在求解力学问题时，应力解总是比位移解 低一阶，导致应力值精度不高，但在实际问题中更关注于应力值，比如应力集中 部位及最大值。 相比于有限元法，以边界积分方程为基础的边界元法( b e m ) 【4 5 】是一种很有效 的数值方法，它能弥补有限元法存在的上述缺点。例如，边界元法只需要对物体 的边界进行离散，使求解问题域降低一级，这不仅便于模拟复杂的几何形状，而 且在很大程度上简化了分析和计算程序。在边界元法中采用了所分析物理问题的 解析基本解，通常具有更高的计算精度。对力学问题，数值计算结果的位移与应 力具有同等的精度。另外，边界元法可以方便地分析无限域和奇异性等问题。 边界元法具有上述优良特点，可现在为什么没有得到广泛的应用呢? 因为， 一直以来有下面的三大问题限制着边界元法的发展。l 、需要所求物理问题的解析 基本解，而对非线性问题很难找到相应的解析基本解。2 、要处理由基本解带来的 奇异积分和近奇异积分。3 、方程系数矩阵是满阵，限制了问题的求解规模。但近 几年，双向互异法、快速算法的提出和应用，精确求解奇异积分和近奇异积分技 术的发展和完善，为解决以上三大瓶颈问题提供了适合的途径。这也为边界元的 基十参数曲面的边界面泫研究 发展提供了良机。 在传统的边界元法中，把三维c a d 几何模型离散成边界元分析模型后，c a d 模型的原始几何信息基本丢掉。边界元分析的几何是基于网格单元通过l a g r a n g e 或h e r m i t e 插值方法近似的。基于网格单元的几何插值引起的几何误差，从根本上 会导致计算精度问题【6 1 ，甚至对有些计算起着决定性的影响。边界元分析模型和 c a d 几何模型的分离，使设计和分析成为两个互相独立的过程。在边界元自适应 网格细分过程中，需要反复地与c a d 系统进行交互，而每个阶段的交互是很复杂 的。因此，尽管自适应网格细分技术理论已经成熟，但主要局限于学术研究，很 少深入到工程分析应用领域【6 】。 上述问题引起了工程师和数值研究者的极大兴趣，国内外许多学者对其进行 了研究。国内，张见明等人【7 。1 5 】提出了杂交边界点法( h b n m ) ，边界积分和变量插 值直接在基于边界表征的实体边界参数曲面上进行，该法成功地实现了与c a d 几 何的无缝连接。国外，l w a n g 1 6 】通过几何细分方法使边界元法与c a d 几何融合一 体，但几何细分方法是以单元为基础的一种插值逼近的几何建模方法，从根本上 会带来几何误差。类似地在有限元法中，t j r h u g h e s 等人【6 】提出了基于非均匀有 理b 样条( n u r b s ) 的等几何分析方法，已经运用在线弹性固体和结构力学及流体力 学中。这种方法现只运用在完整的b 样曲面上，但还未见相关文献报道运用在复杂 的裁剪曲面上。b h a v y a a g g a r w a l r 7 】利用b 样单元分析板弹性问题。y b a z i l e v s 等人 【”】提出了基于t s p l i n e 等几何分析方法，并成功地运用于二维及三维流体和结构分 析。j u r i j sb a z i l e v s 1 9 l 利用基于n u r b s 等几何分析方法分析流固耦合问题。f c i r a k 等人【2 0 】利用曲面单元分析薄板变形问题等。 等几何分析方法已成功地运用在有限元法的诸多领域中，但在有关边界元法 中的报道较少。本文基于参数曲面，提出了边界元法的一种新的实现方法，在这 称为边界面法( b o u n d a r yf a c em e t h o d ，b f m ) 2 1 - 2 3 】。该方法继承了传统边界元法的 特点。重要的是，这种方法具有等几何分析的许多优良特点，而且可以做到与c a d 软件无缝集成。 1 2 边界元法简介 边界积分方程一边界元法( b o u n d a r yi n t e g r a le q u a t i o n b o u n d a r ye l e m e n t m e t h o d ) 简称边界元法( b e m ) 。该方法的产生可追溯到1 0 0 多年前。1 9 世纪就有人提 出了一些积分等式和位势理论，可把线性偏微分方程的边值问题转化为等价的边 界积分方程求解，而积分解的存在性和唯一性等问题早已得到解决。随后边界元 发展较缓慢。在2 0 世纪7 0 年代后，随着有限元法的同趋成熟，以及其缺陷同益明 显，人们将目光转向了边界元，并将有限元法的离散技术运用在边界元法中，从 而成为一种比较成熟的工程分析工具。应该指出的是，英国科学家c a b r e b b i a 4 】 所编著的边界元方法著作对该方法的推广和应用取得了很大的促进作用。 该方法以边界积分方程为数学基础，同时采用了有限元法的单元离散技术。 它通过将边界离散为边界单元，将边界积分方程离散为线性代数方程组，再由数 值方法求解线性代数方程组，从而得到原问题的边界积分方程的解。边界元法具 有边界积分方程的深厚数学根基，它又是在计算机飞速发展的前提下，在有限元 法之后发展起来的一种有效的工程数值分析方法。它不仅可以用于分析弹性力学 等固体力学问题，同时也应用于流体力学、热传导、以及电磁场等其它物理领域。 边界元法的最大特点是降低了求解问题的维数，将二维问题化为其边界线上的 一维问题，将三维问题化为边界面上的二维问题。例如，分析二维带孔矩形板问 题时，图1 1 和图1 2 比较了边界元法和有限元的离散模型，边界元法只需要求解域 的边界离散因而简单的多。该方法只以边界变量为基本变量，域内未知量可以在 需要时根据边界变量求出。这种方法通常具有较高的精度，而且在一些情况下比 有限元法更为有效。因此这种方法受到国内外许多专家学者的重视，并且随着相 应软件的发展与推广，逐步在工程上得到了应用。 图1 1带孔板边界元离散模型 图1 2 带子l 板有限元离散模型 我国边界元法的研究大约开始于2 0 世纪8 0 年代，清华大学杜庆华院士首先在 国内丌创了工程边界元法的研究。3 0 年以来，随着国际边界元法的发展，我国学 者也取得了很大的成绩。自19 8 5 年由杜院士发起并主持全国第一届边界元会议以 来，到目前为止，我国已举行了多次这样的会议及相关的国际会议。这些会议的 召开促进了我国边界元法的研究工作。 1 3研究意义 边界元法作为以计算机为基础的数值计算方法运用于工程分析起源于6 0 年 代。随后，边界元法广泛地运用在各种工程数值分析领域，并在8 0 年代基本成熟。 我们知道，c a d 起源还晚点，在7 0 、8 0 年代。尽管几何是分析的支撑基础，但为 什么用于边界元法分析的模型跟c a d 模型如此不同。许多边界元法的计算程序在 技术上成熟远早于c a d 模型的广泛应用。这很容易让人提出这样一个问题，是否 可以将分析模型与c a d 几何模型融为一体呢? 答案当然是肯定的。 本文基于参数曲面，提出了边界元法的一种新的实现方法一边界面法( b f m l 2 1 - 2 3 】。其最主要的目的是使分析模型与c a d 几何模型融为一体。不管以多么粗糙 的网格离散，分析模型在几何上是精确的，以避免由几何误差的影响降低计算精 度。另外一个重要的目的是，在自适应过程中消除与c a d j l 何模型的反复交互。 一旦初始网格形成，自适应网格细分变得简单。而另一个目的是在c a d 几何模型 上使网格生成更加轻松简单。 边界面法是建立在以数学解析表达的参数曲面基础上。不论是对边界的积分 和还是对场变量的插值都是在曲面的参数空间里进行。直接在曲面的参数空间内 进行边界积分，直接利用c a d 造型系统中参数曲面的几何信息是本文方法与当今 主流c a e 软件的重要不同之处。有限元法和传统边界元法中，函数插值和数值积 分( 能量积分或边界积分) 都是在单元内进行，且必须依赖于单元。本文方法中需要 分析计算的几何变量直接来自c a d 几何模型，因而自然地与实体造型系统融为一 体。 1 4 研究的主要内容 本文围绕边界面法的成功实现，研究的主要内容分以下四个方面： 1 改进传统边界元法的积分和插值方法 ( 1 ) 实现插值区域与积分区域的相互独立 在传统的边界元中，首先将边界离散成单元，数值积分和几何插值都在单元 里进行。在单元内利用l a g r a n g e 或h e r m i t e 多项式对物理场变量进行插值逼近。本 文尝试让数值积分区域与变量插值区域分离，使积分和插值域相互独立。这样， 变量插值不一定完全依懒于单元，使插值方法更加灵活。既可以做到基于单元用 分段多项式插值，也可以基于整个曲面利用移动最小二乘法( m l s ) 或非均匀有理b 样条( n u r b s ) 等方法进行插值逼近；也可以方便地融合连续插值和非连续性插值。 ( 2 ) 实现在曲面参数空间的积分和插值 本文尝试在实体边界的曲面参数空间旱进行边界积分和变量插值。试图避免 因基于单元用多项式插值近似几何引起的几何误差。本文方法的几何数据直接通 过参数曲面的参数变量计算获得。在参数曲面上，可以利用f h j 面单元、移动最小 二乘法或非均匀有理b 样条等方法对变量进行插值逼近。本文着重研究基于曲面单 元在曲面的参数空间里对变量插值的方法。 2 精确计算奇异积分和近奇异积分 本文拟采用局部坐标变换的方法尝试消除奇异性，提高积分精度。为处理近 奇异积分，拟采用基于参数空间的单元细分技术开发自适应积分方案。该方法是 不依赖于所求解的物理问题，是一种具有一般性的通用方法。 3 实现任意曲面网格生成方法 本文的方法必须要对实体的边界曲面在参数空间里离散为曲面网格。能否生 成质量好的网格是实现边界面法的关键。本文拟选择当今流行的曲面网格生成技 术之一一推进波前法( a d v a n c i n gf r o n tm e t h o d ，a f m ) ，生成较好的曲面网格。在 掌握该技术的一般原理后，研究确定相应的数据结构和网格质量参数，使网格能 满足边界面法所需边界离散的要求。 4 实现与实体造型c a d 系统的无缝连接 边界面法只需要边界表征的几何实体的参数曲面，而大部分c a d 商业软件( 如 u g 、p r o e ) 中的几何实体采用边界表征。因此，边界面法可以方便地做到与c a d 软件的无缝连接，实现c a d 与c a e 的集成。本文尝试基于u g ：次开发，将边界面 法的算法程序接入u g 实体建模环境中，调用u g o p e na p i 的相关函数获得实体边 界曲面的几何信息，为后继边界面法分析所利用。这样，可以方便地分析具有复 杂形状几何的实际工程问题。 皋于参数曲面的边界面法研究 2 1引言 第2 章边界积分方程及变量插值 描述工程物理问题的控制方程一般都是偏微分方程，称为数学物理方程，如 拉普拉斯( l a p l a c e ) 方程、泊松( p o i s s o n ) 方程等。这些微分方程是泛定的，只是同一 类物理现象的共性，要获得这些方程具体问题的解，还要加上定解条件一一初始 条件和边界条件。 本文的边界面法跟边界元法一样，第一步是根据微分方程的定解问题，推出 相应的边界积分方程。c a b r e b b i a 4 】为了与有限差分法、有限元法相联系，并统 一在相同的数学基础上，提出了以加权余量法( w e i g h t e dr e s i d u a lm e t h o d ) 作为边界 元法的出发点推导边界积分方程。在数学上，把一个区域上的积分转化为区域边 界上的积分根本方法是基于格林公式的应用。本章基于格林公式中的基本原理， 简洁推导出与拉普拉斯控制方程等价的自身正则化的边界积分方程。 在有限元法和传统的边界元法中 3 - 5 】，对场变量是基于单元用分段多项式插值 近似的，例如l a g r a n g e 或h e r m i t e 插值方法。该插值方法实现必须依赖于二维或三 维物理空间的单元。有限元法中的能量积分和边界元法的边界积分也是在物理空 间的单元里实施的。也就是说，积分和变量插值区域相同，都在同一单元内进行。 由张见明等人【7 - 1 5 】提出的杂交边界点法( h d b n m ) ，t j r h u g h e s 等人f 6 1 8 】提出 的基于非均匀有理b 样条的等几何分析方法和a m i ts h a w 等人【2 4 】提出的基于非均 匀有理b 样条的无网格方法等内容的启发，本文中对变量的插值，直接在边界曲面 的参数空间里进行，并且使积分区域和插值区域相互独立。两者区域的独立，可 以自由得选择插值函数。可以引用无网格法中的变量近似方法构造形函数，使边 界面法具有无网格法的特点；也可以利用传统边界元法和有限元法的分段多项式 插值方法，使边界面法具有基于网格插值方法的特点。从而，边界面法提供了一 个自由而统一的变量插值框架，使该方法的实现灵活自由。边界面法已成功在曲 面参数空间里利用移动最小二乘变量近似方法实现了【2 1 1 ，使该方法具有像边界点 法( b n m ) 等无网格法【7 - 1 5 , 2 5 - 2 9 】的特点。 另外，由本章2 2 2 节中离散的边界积分方程( 2 2 4 ) 可知，形函数m ( s ) 不要求连 续性，因而不必保证在边界曲面上处处连续。但是，在有限元法中必须保证形函 数在整个求解域是一阶连续性的。这就在理论上允许我们可以很自由地用非连续 的非结构网格离散曲面边界。如果用参数空间的曲面单元插值变量，在相邻曲面 上共享棱边的单元不要求插值节点共用，可以在各自的曲面边界上利用非连续单 元( 图2 1 ) 。因此，可以对不同的曲面用不同密度的网格来离散实体边界，如图2 2 所示的立方体边界离散。这有利于求解薄型或含有细小特征结构的问题，也利于 避免边界奇异性处理多域问题。 图2 1平面上的节点分布图2 2 不同的单元密度网格 在参数曲面上，对物理场变量的插值逼近有多种方法。张见明等人【7 - 1 5 1 提出了 基于参数曲面上的移动最小二乘插值方案。a m i ts h a w 等人利用非均匀有理b 样 条插值逼近未知量。在边界元法中，节点只分布在边界曲面上，以上两种方案只 需分别在每个独立的边界曲面上实施，不依赖单元。用以上两种方法开发出的边 界类型的无网格法，虽然插值不需单元，但仍然需要积分单元。不管是利用移动 最小二乘法还是非均匀有理b 样条插值，都能获得较高的精度，但这两种方法很 难用于复杂的裁剪曲面和具有狭长曲面的薄型结构。另外，需要选择经验参数， 计算效率较低，不适宜用于求解大规模复杂的工程问题。为此，本文利用分段多 项式在二维参数曲面单元内对未知量进行插值，类似于有限元法中的二维平面单 元【3 1 。 2 2 边界积分方程及积分方案 2 2 1正则边界积分方程的推导 对任意d ( d = 1 ，2 ，3 ) 维域q ，满足拉普拉斯控制方程的位势边界值问题可以表 示为： 地2 = 喜筹- o ，沌q u = f i - ，v x f 。 ( 2 1 ) u ，fn i 兰q = - 4 ，v x f 。 式中q 的整个边界为f = f 。+ f 。，厅和虿分别是位势已知边界r 。和法向流已知边界 f 。的边界值，n 是边界外法向矢量，n i 是法矢量分量，x i 是d 维物理空间坐标分量， 基于参数曲面的边界面法研究 i = 1 ，d ，如图2 3 所示。 u r = c n l = 9 把上述问题转化为等效积分形式： l 刃2 谢q = l v 蔷a 犁0 2 u q = 。( 2 2 ) 其中y 为任意试函数，对以上方程的左端项做一次分部积分变换，可得： 嘻争一喜毒挚+ 瑶孙珈 泣3 ， 一喜善挚叫r v 鼢如智挑魄 玎锄 对( 2 3 ) 方程的右端的第一项继续用分部积分， 喜瓦钆o ut q = 唾喜缸n 喀椭 泣4 ， 将( 2 4 ) 式代入方程( 2 3 ) 得： 谢q 一4 ( 未妻户f = 。 位5 ， 由于v 是任意的，取v 为满足拉普拉斯方程基本解矿，满足 俨! i 。 包6 ， ，= 0 s y l l 式中，万( 厂) 为狄拉克函数，y ，s 分别为边界上的源点和场点，r 为原点和场点间的 欧氏空间距离，。矿是基本解，它表示在无限域中源点y 处有单位点源时在任意场 点s 处产生的势，如图2 4 所示。当y s 满足拉普拉斯方程，当y = s 基本解奇异， 因此基本解又称奇异解。对二维问题普拉斯方程的基本解为： “s ：上l n j l ( 2 7 ) 2 万 r ( s ，y ) 对三维问题普拉斯方程的基本解为： 矿( s ，y ) _ 去志 q 沼 图2 4 位势问题的源点y 与场点s 考虑8 ( r ) 函数的性质 肌y 胂= 心谠宇 方程( 2 5 ) 可以表达成以下两个方程 ( 2 9 ) 吣) = 唾( 掣坤一( y ，s ) 警卜) ，v y e f 2 ( 2 1 0 ) 。= 4 ( 掣心一( y ，s ) 警一v y 磊 ( 2 1 1 ) 以上两方程给出了点y 在q 和域外壶的积分方程。如果使点y 位于边界上，把方程 ( 2 5 ) 中所有的未知变量转化到边界上，就可以得到边界积分方程。下面给出得到 正则化边界积分方程的推导。 幽2 5 推导边界积分方程示意图 如图2 5 所示，对边界上任取一点p ，方程( 2 1 0 ) 等价于下面的方程 珩，2 “掣【u ( s ) - u ( p ) 渤掣卜q + “( p ) 掣r ( s ) 如果y 趋近与p ，并与p 重合时，可以把所有的未知变量转化到边界上， 积分方程。现在考虑方程( 2 1 2 ) 右端的第二项边界积分，对二维问题有 唾挚f = 唾去挚r = 去4 象r 如图2 6 ，露为边界外法向量，两向量尹和万的夹角余弦值为 ( 2 1 2 ) 得到边界 ( 2 1 3 ) 皋于参数曲面的边界面法研究 又因为 所以有 n 一r 一元尹 咖归丽2 可 红a n 知考2 骨强+ 骨嘭2 首 一2 瓦+ 瓦2 矿强+ 1 矿嘭2 两 堡订：c o s 口d f o n 则 唾言r = 唾吾c o s 口d r = 唾扣o = 唾d o = 2 万 将式( 2 1 7 ) 代入( 2 1 3 ) ，可得 唾塑玺堂订= 1 1 1 1 2 6 计算积分4 要订示意图 同理，对三维位势问题也可以推出该结论也成立。 所以，当p 与y 重合时，可以得到 o = d f l 型o n 叫y 帅y ，s ) 掣卜) ( 2 1 9 ) 令v ( y ，s ) = ”5 ( y ，s ) ，q = a u a n ，可得以下正则化边界积分方程 0 = f ( 甜( s ) 一“( y ) ) 9 5 ( s ，y ) d f f g ( s 5 ( s ，y ) d r ( 2 2 0 ) 式中q s ( s ，y ) 为基本解的单层势梯度，对三维位势问题表达为 g坳)=掣删=谤-而1丽ar(s,y)on ( 2 2 1 ) q 万厂一i s v i 砌i s i 方程( 2 2 0 ) 之所以称为自身正则化边界积分方程，是因为方程右端的第一项，当场 点s 趋近源点y 时，q s ( s ，y ) 项的分母趋近于o ，而“( s ) 一u ( y ) 的值也趋近于0 ，从而消 除了该项被积函数的奇异性。 到此，满足拉普拉斯控制方程位势边界值问题( 2 1 ) 转化为边界积分方( 2 2 0 ) 。 本文基于边界表征实体造型方法提出的边界面法，就是通过数值方法求解该边界 积分方程。本文讨论该方法只考虑三维情况，事实上在真实的物理世界中，都是 旧 q 眨 亿 三维的，二维情况只是针对三维空间的一些特殊问题作降低维数的合理简化。 在边界值问题中，边界条件有以下三种：给边界施加位势( 比如，温度、位移 等) ，称为本质边界条件( d i r i e h l e t ) ；给边界施加通量，即法向流量( 比如，温度法 向流量，力等) ，称为自然边界条件( n e u m a n n ) ；在一部分边界施加本质边界条件， 而另一部分施加自然边界条件( c a u c h y ) ，或在给定位势和通量的线性组合 ( r o b i n ) ，都为混合边界条件( m i x e d ) 。 2 2 2 边界积分方程离散 用边界面法求解边界积分方程( 2 2 0 ) 式，首先需萼将区域q 的边界f 离散成若 干个( 以) 边界单元f j ，和相应的个插值节点。然后，在边界节点或给定点上施 加相应的边界条件。设节点七上的位势和法向流值分别为，g 。对任意节点k ， 要么u 。已知，要么g 。已知。将边界上任意一点s 的位势“( s ) 和法向流q ( s ) 分别表示 为个节点的值插函数，即 甜( s ) = “( 叩) = m 魄 仨1 ( 2 2 2 ) g ( s ) = “( “，1 ，) = m 吼 式中，( “，1 ，) 为边界参数曲面的二维参数坐标，三维物理空间坐标( x ，y ，z ) 满足 x = x ( “，1 ，) ，y = y ( “，1 ，) 和z = z ( u ，d ，m ( s ) 为形函数。将边界离散后，将( 2 2 2 ) 式代入 方程( 2 2 0 ) ，可以把边界积分方程离散为关于个节点未知量的方程，即 n 。nn 。w o = 一芝f 。q s ( s ，y ) ( m ( s ) 一m ( y ) ) 蚝d r + 芝f “5 ( s ，y ) m ( s ) 吼d r ( 2 2 3 ) 方程( 2 2 3 ) 可以组装成矩阵形式为： q 。 h n l q ： 马： h n 2仨 g i g 2 。 ： g ， g i ： g 2 ： ： g n 2 g i g 2 g 懈 或hu=gq 式中， 纸：笠f ，9 5 ( s ，y 从以( s ) 一m ( y d r 瓯= n u mf ，“5 ( s ，y ，) m ( s ) d r 上式中，f ，为形函数m ( s ) 值不为。的边界单元，数量记为n u m 。 将方程( 2 2 5 ) 进行变换，使未知量移到左边，已知量移到右边， a x = b ( 2 2 4 ) ( 2 2 5 ) ( 2 2 6 ) ( 2 2 7 ) 形成线性组， ( 2 2 8 ) w； 一 一 一 基于参数曲面的边界面法研究 式中，x 是“。或g 。的未知数向量。求解方程( 2 2 8 ) ，就可以得到所有节点 k ( k = 1 ，2 n ) 上未知量u 。或g 女的值。根据节点上的值，相应地可以求得边界上和域 内任意一点的位势和流量值【4 ，5 1 。 2 2 3 边界积分方案 用数值方法求解( 2 2 0 ) 式边界积分方程，首先需要将边界1 1 离散成单元，然后 考虑源点y 对每个单元进行相应的数值积分。方程( 2 2 0 ) 右端的第一项是正则的， 因此，对任何边界单元可以用正则高斯积分计算这一项。考虑积分方程中位势基 本解，当场点s 趋近源点y 时，r ( s ，y ) 一0 ，从而“5 ( s ，y ) - - - o o ，这就导致在局部小范 围内，值的变化相当大。因此不能用正则高斯积分计算方程右端的第二项，需要 特殊的积分技术。 当源点y 位于积分单元的顶点或单元的内部，方程( 2 2 0 ) 右端的第二项产生弱 奇异性。在不同的边界单元内，当场点s 接近源点y 时，方程( 2 2 0 ) 右端的第二项 就产生近奇异性。当求解薄形结构或边界单元分布很不规则时，就存在这种情况。 当场点s 与源点y 相隔一定的距离，就可以用正则高斯积分计算，该距离由场点s 所 在单元的尺寸和场点与源点间距离等因素决定。处理弱奇异性和近奇异性的积分 技术将在第三章给出详细介绍。因此，考虑场点s 与源点y 的相对位置对基本解带 来的奇异性，边界积分可分为正则化积分、奇异积分和近奇异积分。 本文方法的实现是直接基于边界表征的实体造型数据结构，实体的边界曲面 都以参数形式表达的。边界积分直接在曲面二维参数空间里完成，没有依赖传统 的边界单元。在数值积分过程中，被积函数的几何变量，比如高斯积分点的坐标、 雅可比、外法向量是直接通过参数曲面单元中的曲面参数变量计算获得，而不是 通过单元插值近似的。直接在曲面的参数空间内进行边界积分，直接利用c a d 造 型系统中参数曲面的几何信息是本文方法与当今主流c a e 软件的重要不同之处。 有限元法和传统边界元法中，函数的插值和数值积分( 能量积分或边界积分) 都是在 单元内进行，且必须依赖于单元。本文方法中需要的分析计算的几何变量直接来 自c a d 几何模型，因而自然地与实体造型系统融为一体。 2 2 4 域内和边界上点位势和流量 求方程( 2 2 8 ) 后，边界上插值节点未知量u 。或q 。值已知，则边界上的位势和法 向流量可以由( 2 2 2 ) 式近似计算得到。为计算边界上位势的梯度，可以由以下方程 1 1 0 1 计算： f 牙( 一c ) 1 刀。 ，l ： 以， f a a x , 1 a c # a u = i 钆锄锄魏锄i a “ ( 2 2 9 ) 【a 历o vjl 钆加西a x , l o v 儿a u l o x , j 其中，c 为传导率( 或渗透率) ，在2 2 1 节边界积分方程推导过程中，取c = 1 0 ： 吃，i = 1 ，2 ，3 是边界外法向量分量；g f = 加挑为位势的梯度；魄a , 和阮a , 由参数 曲面方程直接计算获得；笳l a s , 可以由下式计算得到 驯锄= m ，。，驯加= m ，u k ( 2 3 0 ) 域内任意一点y 的势u 和流量q 通过以下的积分方程计算 “( j ，) = l ”o 5 ( s ，y ) d r l g o 5 ( j ，y ) d f = 蔓h w “石r 一扣仅y ) d r q 3 d e l e m e n t s e l e m e n t s 吣卜c ，a u - ( y - 2 = - c 工哿州州r - f 等吲州r 2亿3 2 ， 叫磊盥f 。哿似蜩一磊盥工。帮吲归 其中，e l e m e n t s 为边界离散后的边界单元，数值积分在每个单元内进行。 如果当点y 接近边界，会产生近奇异积分，不能用正则高斯积分直接对边界单 元进行积分。对相应的单元应采用自适应近奇异积分技术计算其积分，该技术的 实现在第三章详细给出。 2 3基于曲面单元变量插值 2 3 1线性三角形曲面单元 首先以三角形曲面单元( s t 3 ) 为例，说明在曲面二维参数空间上基于单元怎样 实现变量插值。为方便处理边界奇异性及多域问题，曲面边界上单元插值节点向 曲面内做较小的偏移，生成非连续单元( 图2 7 b ) ，而在曲面的内部都采用连续单元 ( 图2 7 a ) 。图中空心圆点表示单元顶点，而实心圆点表示插值节点，该表示在文中 以下的情况相同。 元顶点 值节点 图2 7 三角形曲面单元：a 连续单元，b 非连续单元 假设三角形曲面单元三个插值点( “，v ) 上场变量值为谚( i = l ，2 ，3 ) ； 点变量值g ( u ，v ) 由以下具有完备基线性多项式函数近似 矽 ，1 ，) = a n + b v + c 三个插值节点上的变量满足方程( 2 3 3 ) ，所以有 单元内任意 ( 2 3