（固体力学专业论文）边界元法反求解二维位势问题的边界条件.pdf

边界元法反求解二维位势问题的边界条件 摘要 在工程结构中经常会遇到这样的问题，即部分边界上的物理量不便直接测 量，而是通过其它边界上己知量来确定。从数学的角度来看，这属于反问题的 一种。对于复杂结构的反问题，通常采用数值方法求解，如有限元法、有限差 分法、边界元法等。因为边界元法仅利用边界量进行求解，故对于上述一类反 问题分析有较大优势。 过去边界元法求解薄体结构存在着几乎奇异积分的困难，目前几乎奇异积 分问题己得到成功解决，解析积分算法也己成功应用于分析二维位势问题常规 结构和薄体结构。本文在此基础上，首先介绍了几乎奇异积分及奇异值分解方 法的基本理论，引用并分析了二维位势边界元法中几乎奇异积分的解析积分算 法，再采用奇异值分解方法，对常规和薄体结构位势反问题实施了正则化。特 别对二维薄体位势边界元法，成功实现了几乎奇异积分和不适定反问题的双重 正则化。本文用边界元方法计算了二维位势各向同性常规结构和薄体结构以及 二维位势正交各向异性常规结构及薄体结构的c a u c h y 边值问题。文中给出了较 多数值算例，证明了该法的有效性。 关键词：反问题，边界元法，几乎奇异积分，解析积分，奇异值分解，位势问 题 i n v e r s ed e t e r m i n a t i o no ft h eb o u n d a l 了c o n d i t i o n s o f t w o d i m e n s i o n a lp o t e n t i a lp 1 o b l e m s b yu s i n g t h eb o u n d a r ye l e m e n tm e t h o d a b s t r a c t i ti so f t e ne n c o u n t e r e dw h hs o m ep r o b l e m si ne n g i n c e r i n g 血a t 也ep h y s i c a lp a r a m e t e r s o na p a n o f m eb o u n d a r ya r en o td e t e m i n e dd i r c c t l y 瓶di tc a nb cw o r k e do u tb ym e 肌so f t h ek n o w nc o n d i t i o no nt 1 1 eo t h e rb o u n d a r y f r o mt h ev i e w p o i n to fm a t h e m a t i c s ，t l l i si sa k i r l do fm v e r s ep r o b l e m ，w i t hr e g a r dt oc o m p i i c a t e di n v e r s ep r o b l e m s ，n u m e r i c a lm e t h o d s a r ea l w a y sa d o p t e dt od e a lw i mt h e m ，s u c ha sf i n i t ce l e m e n tm e t l l o d ，f i n i t ed i 毹r e n c e m e t h o da n db o u n d a r ye l e m e n tm e t h o d b o u n d a r ye k m e n tm e t h o dh a sag r e a ta d v a i l t a g e o v e ro t l l e r so nm ei n v e r s ep r o b l e m sa b o v ei nt h a tb o u n d a 町e l e m e n tm e t h o d 、v o r k so u tt h e r e s u l t so n l yb a s e do nt h eb o u n d a l yd a h m i nt 1 1 ep a s t ，t h e r ee x i s t sd i m c l l l t yf o rc a l c i l l a t i n gn e a r l ys i n g u l a ri m e g 豫1 sa r i s i n g 丘o m b o u n d a r ye l e m e n ta n a l y s i so f t h ep r o b l e m so f t l l i nb o d y h o w e v e r ，b yn ow ，也ed i 伍c u l t yo f n e a r l ys i n g u l a ri n t e g r a l sh a sb e e nr e s o l v e db ya n a l y t i c a 圭i n t e g r a la l g o r i t l l l n ，w h i c hh a sb e e n p u ti n t o u s e si n a i l a l y z i n gm ec o n v e n 廿o n a l s t n l 曲【l r e sa i l dt l l i nb o d ys 仃u 帆l r e si nt h e p o t e ：呲i a lp r o b l e m s b 硒e d o n 也i s t h e o r y ，t l l ep r e s e mp a p e rf i r s t l y i n t r o d u c e st h e 胁d 锄e n t a lt h e 嘶e so ft h en c a r l ys i n g u l a ri m e g r 出sa n dt h es i n g u l a rv a i u e d e c o m p o s i t i o n ， a 1 1 d q u o t e st h e 删”i c a la l g o r i t l l mo fm en e a r l ys i n g u l a ri n t e g r a l so ft w o d i m e n s i o n a l p o t e n t i a lb o u n 出【r y e l e m c n tm e m o d t h e nt h e p o t e n t i a l i n v e r s e p r o b l e m s o ft h e c o n v e n t i o n a ls m l c t u r e s 锄dm i nb o d ys t n l c t u r e sa r e r e g u l a r i z e db yt l l es i n g u l a rv a l u e d e c o m p o s i t i o n e s p e c i a l l yf o rt h ep r o b l e mo f t h i nb o d y d u p l er e g u l a r i z a t i o no fb o t l lt h e n e a r l ys 协g l l l a ri n t e g m l s a n d i l l - p o s e d i n 憎s e p r o b l e m s i sa l s os u c c e s s f u li nt h e t w o - d i m e n s i o n a lb o 啪d a r ye l e m e mm e t h o d b a s e do nt h cs t r a t e g y ，t h eb o u n d a r ye l e m e n t m 甜1 0 di s a p p l i e dt 0 瞰a l y z i n gt l l e i n w r s ep r o b l e m so ft h et 、v o d i l i l e n s i o n a l c a u c h y p 删a l 、i t l li s o t m p i ca n do n l l o 缸d p i cm a t e r i a l sf 研c o n v e n t i o n a l 曲m t i l r e sa 1 1 dt h i m b o d y s m l c t u r e s a g o o d n 眦l b e ro f n u m e r i c a le x 锄p l e s a r e 西v e n a n dd e m o n s 仃a t et 1 1 ev a l i d i t yo f 吐j ed n 悠曲tn l e 西d d k e y w o r d s ：i n v e r s ep r o b l e m ，b e m ，n e a r l ys i n g u l a ri n t e g r a l ，a n a l y t i c a li n t e g r a l s i n g u l a rv a l u ed e c o m p o s i t i o n ，p o t e n t i a lp r o b l e m 合肥工业大学 本论文经答辩委员会全体委员审查，确认符合合肥工业大学硕 士学位论文质量要求。 答辩委员会签名( 工作单位、职称) 主席：丝傀红 委员：易斓 构饷睁 ；膨罨 戡痧主、 导师 彩嶙一 j 星叛患砻如 沈王影走耋 插图清单 图3 1二维普通位势问题示意图1 2 图3 2线性元示意图1 4 图3 3普通位势问题1 6 图3 4c a u c h v 边值问题1 6 图3 5方形闭区域热流问题1 8 图3 6方形闭区域热流问题a d 边的位势计算结果19 图3 7方形闭区域热流问题a d 边的位势梯度计算结果19 图3 ，8圆环区域热传导2 0 图3 9环形域热传导问题的位势计算结果2 1 图3 1 0 环形域热传导问题的位势梯度计算结果2 1 图3 1 1 环形域热传导问题的计算结果误差2 2 图3 1 2 线性单元2 3 图3 13 接近度e 2 3 图3 1 4 薄体矩形区域热流2 5 图3 1 5 例3 3 位势计算结果( 6 口= 1 e _ 8 ) 2 5 图3 1 6 例3 3 位势梯度计算结果( 6 d = 1 e 一8 ) 2 6 图3 17 例3 36 ，口取不同值时位势梯度计算结果比较2 6 图3 18 例3 4 薄壁圆环2 7 图3 19 薄壁圆环热传导位势计算结果2 7 图3 ，2 0 薄壁圆环热传导位势梯度计算结果2 7 图4 1正交各向异性问题2 9 图4 2 边界线性单元3 0 图4 3正交各向异性方形区域温度场问题3 2 图4 4正交各向异性方形区域温度场a d 边位势计算结果3 2 图4 5正交各向异性方形区域温度场a d 边位势梯度计算结果3 2 图4 6正交各向异性方形区域温度场a d 边位势计算误差3 3 图4 7 正交各向异性方形区域温度场a d 边位势梯度计算误差3 3 图4 8圆形域温度场的边界条件3 4 图4 9圆形域热传导位势计算结果3 4 图4 1 0 圆形域热传导位势梯度计算结果3 4 图4 ，1 1 圆形域热传导位势梯度计算结果( 共3 6 0 单元) 35 图4 1 2 薄壁圆环区域热流3 8 图4 13 例4 3 不同厚度比位势计算结果比较3 9 图4 1 4 例4 3 不同厚度比位势梯度计算结果比较3 9 图4 15 例4 3 位势结果误差比较4 0 图4 1 6 例4 3 位势梯度结果误差比较4 0 表格清单 表3 1 方形闭区域热流问题a d 边计算误差2 0 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作及取得的研究成果。据我所 知，除了文中特别加以标注和致谢的地方外，论文中不包含其他人已经发表或撰写过的研究成果， 也不包含为获得 盒目王些盔堂或其他教育机构的学位或证书而使用过的材料。与我一同 工作的同志对本研究所做的任何贡献均已在论文中作了明确的说明并表示谢意。 学位论文作者签名：才为鸭努 签字日期：厂年翻日 学位论文版权使用授权书 本学位论文作者完全了解金鳇王些太堂有关保留、使用学位论文的规定，有权保留并向国 家有关部门或机构送交论文的复印件和磁盘，允许论文被查阅和借阅。本人授权金8 b 王些盍堂可 以将学位论文的全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或扫描等复制手 段保存、汇编学位论文。 ( 保密的学位论文在解密后适用本授权书) 学位论文作者签名：栖雾 签字日期：护f 年6 月f g 日 学位论文作者毕业后去向 工作单位： 通讯地址： 导师签名；饬 签字日期p 岁年月伽 电话： 邮编： 致谢 本论文是在导师牛忠荣教授及周焕林副教授悉心指导下完成的。 首先我要感谢牛老师三年来对我的辛勤教育和培养。在这三年里，不管是 在学习还是在生活上，牛老师都给了我极大的关心和帮助。在学习与论文工作 中，牛老师一方面严格要求，另一方面为我们创造了良好的学习与工作环境， 并对论文作了大量修改；在生活上，我能感受到牛老师就像父母关心孩子一样 关心自己的学生，千言万语难以表达我内心里的感激之情。 作者同时感谢周焕林副教授在论文阶段对我的工作以及论文的完成给予的 热情帮助及悉心指导，感谢周老师对我的辛勤付出，我在研究生阶段工作 的完成与周老师的帮助是分不开的，在此我向周老师表达深深的谢意! 在此我还要感谢程长征老师和胡宗军老师对我给予的无私帮助! 我还要对我的任课老师王建国教授、王左辉教授，王炯华教授、刘一华教 授、盛宏玉教授，蔡敏教授、白嘉楠副教授，赵速梅副教授等致以衷心的感谢 和诚挚的敬意! 同时，还要感谢合肥工业大学土建学院和工程力学系的所有老师和同学对 我生活上的关心和学业上的鼓励! 最后，我还要谢谢我的家人以及我所有的朋友! 本文研究工作得到国家自然科学基金( 1 0 2 7 2 0 3 9 ) 的资助。 作者：杨智勇 2 0 0 5 年6 月 1 1力学数值方法概述 第一章绪论 力学是一门基础学科，它所阐明的规律带有普遍的性质；力学又是一门技 术科学，它是许多工程技术的理论基础川。对于力学问题，其控制方程和边界 条件易于确定，但除了一些简单问题可以给出解析解外，绝大多数的问题不得 不求助于数值解。早在计算机出现之前，力学工作者就提出了若干种分析和设 计大型、复杂结构的计算方法，高速有效的数值计算是力学知识通向工程应用 的主要手段。 2 0 世纪中叶，随着计算机问世，计算理论和数值方法得到了迅速发展，有 限元法的诞生，标志着计算力学作为现代力学学科的重要分支之一，登上了国 际力学界的舞台。如今，有限差分法( f d m ) 、有限单元法( f e m ) 和边界单元法 ( b e m ) 是计算力学领域中的几种最重要的数值方法，它们在力学研究和工程应 用诸多领域得到了广泛应用。 有限差分法是一种历史悠久、理论简洁的数值解法【2 l 。它将所考虑的领域 划分适当的网格，在网格的结点上采用差分方程( 代数方程) 近似地代替控制 方程( 微分方程) ，在边界格子点上引用差分方程表示边界条件，把求解微分方 程的问题转换为求解代数方程的问题。当采用较多结点时，近似解的精度可以 得到改进。有限差分法可以求解相当复杂的问题，且能提供方便的误差估计方 法。对求解建立于空间坐标系的流体流动问题，优势明显，在流体力学领域， 它至今仍占主导地位。但是有限差分法很难处理复杂的边界或有应力、应变奇 异的问题，对三维问题也不容易处理。随着其他数值方法的深入发展，有限差 分法在许多问题上的应用已经没有优势。 有限元法把有限差分法的网格点离散改造成更为灵活的有限单元离散，利 用每一个单元内假设的近似函数来分片地表示全求解域上待求的未知场函数 【3 】。根据最小势能原理、虚功原理或者广义交分原理h l 建立有限元方程，充分 利用计算机能力求解未知变量以后，通过插值函数计算出各单元内场函数的近 似值，从而得到整个求解域上的近似解。理论证明，随着划分单元数目的增多 以及插值函数精度的提高，解的近似程度将不断提高。如果单元满足收敛要求， 近似解收敛于精确解。有限元法开拓了现代数值方法的新天地，现己推广渗透到 力学的各分支和许多其他学科领域，用于固体力学的有限元法已日趋完善，近 年来有限元法已发展到流体力学、温度场、磁场和声场等问题的求解，分析对 象也从弹性材料扩展到塑性、弹塑性、粘塑性和复合材料等。如今有限元法程 序有了方便的构模和前后处理功能，并将求解各种问题的程序归入一个程序系 统形成通用的程序包。现在功能齐全的大型通用程序包已经商品化，在科学研 究和工程应用中起到愈来愈大的作用，国际上流行很多功能强大的有限元分析 软件，主要有美国的n a s t r a n 、a b q u s ，德国的a s k a ，英国的p a f e c 等。 有限元法适用范围广，精度高，但也有它自身的缺陷。它和有限差分法一样， 需要在整个求解域上离散，问题的自由度和原始信息量大。对无限或半无限求 解域，有限元法只能采取将求解区域近似取为有限域的办法。 边界元法起源于上世纪六十年代，它继承了有限差分法和有限元法的离散 的思想，其特色在于待解量都限定在结构物的表面( 也叫边界) 上。边界元法 的优点在于：( 1 ) 其降维特点使待解量压缩到很小的规模，数据的前后处理特 别方便；( 2 ) 边界量全部确定后可以按照需要逐点逐向计算任意点的物理量， 而且各种计算值都可以达到很高的同阶计算精度；( 3 ) 由于奇性基本解常在无 限域中给出，当跟特殊构造的单元如无限元、奇异元等相结合，边界元法用来 分析无限域问题、半无限域问题和应力集中等奇性问题尤其有效”】。( 4 ) 最近 几年边界元法研究表明，对复杂几何区域问题，如含有薄体结构，只要有效处 理边界元法几乎奇异积分问题，就可从容处理边界单元的疏密过渡，从而有效 处理薄体等结构问题。 当然，边界元法的短处也是不容忽视的，其一，由于它的系数矩阵为非对 称满阵，因而3 0 以上的枫时都耗在系数矩阵的形成上，因此解题规模受到存 储量和计算时间两方面的制约1 6 】；其二，由于所用基本解和分析途径因问题而 异，这就不利于编写通用程序，而且它要用到较多较深的数学工具，因而大大 限制了它的工程应用。 目前，边界元法应用软件主要有b e a s y ，s u r f e s ，e z b e a ，h i b e m ， b e s m a p 等，但总体上，边界元法程序包的发展现状落后于边界元法的发展水 平，落后于软件技术的发展水平，边界元法及其软件的工程推广应用任重道远 【7 1 。 1 2 边界元法的发展 边界元法实际上是紧随着有限元法发展起来的，和有限元法相同，边界元 法的研究也始于2 0 世纪五、六十年代。最初f j r i z z o 和t ，a c r u s e 将该方 法命名为边界积分方程法( b o u n d a r yi n t e g r a le q u a t i o nm e t h o d ) ，虽然边界积分 方程之说法实际上是数学家提出的。历史上，b e t t i ( 1 8 7 2 ) 提出了贝蒂互等定 理，k e i v i n ( 1 8 4 8 ) 求出了集中力作用于无限大物体的弹性力学基本解。 s o m i g l i a n a ( 1 8 8 5 ) 结合k e i v i n 特解和b e t t i 互等定理得到了用边界位移和面力 表示的内点位移的边界积分公式，从而为弹性力学边界积分方程的建立作出了 理论准备。 边界积分方程法是通过许多俄国学者的工作逐步被人们了解的， m u s k h e i i s h v i l i ( 1 9 5 3 ) 将积分方程方法用于结构力学分析。k e l l o g g ( 1 9 5 3 ) 用 积分方程方法求解l a p l a c e 问题【8 1 。积分方程法在一些流体力学和位势问题中 被广泛采用并称为源法。它现在被称为“间接”( i n d i r e c t ) 边界元法，即未知 量不是问题的物理变量。其中，j a s w o n 对l a p l a c e 方程由势理论建立了边界积 分方程的数值方法，为间接边界元法的提出作出了重要的贡献。关于“直接” ( d i r e c t ) 边界元法，首先是由j a s w o n 和p o n t e r 【9 】( 1 9 6 3 ) 对弹性杆件扭转问 题进行的，他们利用翘曲函数数值求解了扭转问题的抗扭刚度和边界剪应力。 这是有效利用边值和法向导数之间积分关系的最早论文之一。 随着边界积分方程法的发展，1 9 7 7 年j a s w o n 和s y m m 出版了一本位势和 弹性力学积分方程法的专著【1 0 】。由于c a b r e b b i a 在有限元和边界元方面的工 作积累，使他在文【1 1 中第一次创造性地采用了“边界元法”这一名称，1 9 7 8 年b r e b b i a 编写了国际上第一本书名含有“边界元法”的边界元专著【1 2 】，该书 阐明了边界元法和其它数值方法特别是有限元法的关系，提出了如何用加权余 量法来建立边界积分方程，初步形成了边界元法的理论体系，确立了边界元法 作为一种数值方法的地位，标志着边界元法从此进入了系统性的研究时期。 从1 9 7 8 年以来，国际边界元法会议基本上每年召开一次，边界元法的研究 迅速在世界上普及，经过4 0 多年来的发展，边界元法的研究领域已经相当丰富。 固体和流体，线性和非线性，稳态和瞬态，普通材料和复合材料，均质和非均 质，正问题和反问题，常规结构和薄体结构等等都有涉猎。 1 3位势反问题及反热传导问题概述 1 3 1 位势问题 位势问题在数学分析上也属于椭圆型方程问题，它代表的是一种定常状态。 诸如稳定温度场，无旋稳恒流动，稳定浓度场，静电场等等都满足l a p l a c e 方程， 通称位势问题，可统一表示为 其中“= “( z ) 叫做位势，在具体问题上可以是温度、水头、浓度、压力、电位等 等，6 = 6 ( z ) 称为体源；q 是该问题的分析域，r 1 是位势给定边界，r 2 是位势法 向梯度给定边界，”是边界外法线。和弹性力学问题类似，位势问题也只有在 简单状况下才有解析解，对于工程结构中的位势问题一般都采用数值分析方法。 而边界元法对求解位势问题相当有效，这是因为边界元法兼有降维、高精度、 3 适解性强等优点，尤其是对薄体问题。早在1 9 6 3 年j a s w o n 【1 3 1 和s y m m 。4 1 就用间 接边界元法求解了位势问题；如今边界元法在位势问题的研究已相当成熟，从二 维问题到三维问题，从各向同性到各向异性。位势问题也分正问题和反问题，其 中反热传导问题是一种典型的位势反问题。 l ，3 ，2 反闻题及反热传导闻题概述 反问题总是伴随着正问题而产生，从数学角度来看，对于微分方程的定解 问题，亦即正问题，通常是寻求满足方程初始条件与边界条件的解。如果在微分 方程定解问题中，原来已知量中的某一个或几个条件变成未知量，而原来未知的 函数，现在或者是已知的，或者是与这未知函数有关的某些“信息”是已知的( 这 些“信息”称为“附加条件”) ，即由支配方程、初始条件以及附加条 牛去确定 原定解问题中的未知量，此即反问题。 反问题的研究涉及到未知边界条件识别，材料物理常数估算。边界几何形 状重建，裂纹孔洞和缺陷等的测评，还有结构残余应力分析等等。对几何区域 确定的物理问题，如果已知其控制方程，材料常数，但仅已知部分边界条件， 这类问题称为边界条件识别型反问题。 对反热传导问题，如果部分边界上的温度和热流己知，但剩余边界上的温 度和热流不明确，这类部分边界含有o v e r s p e c m e d 边界条件的问题通常称为 c a u c h y 型位势反问题。从上世纪6 0 年代初就有人对反热传导问题进行了研究， 经过4 0 多年的发展取得了很大进展。早在1 9 6 0 年g s t o l z 【1 5 】就首先给出了球体热 传导反问题的数值解，是研究反热传导问题的开创性工作；i f r a n k 1 6 1 首先提出 用最小二乘法解反热传导问题，将物体表面热流用时间的多项式表示，用最小 二乘法定出热流多项式中的系数。j v b e c k i 】研究了一维非线性反热传导问题， 他在该文中提出的数值方法被广泛引用，他采用c r a n k n i c o l s o n 有限差分格式离 散热传导方程，用非线性最小二乘法求解，他的优点是处理不精确的数据稳定 性好：h r b u s b y 等i l 驯用有限元法和动态规划方法研究了二维热传导反问题， 他们用动态规划方法使最小二乘法误差函数对未知的边界热流向量取极小值， 从而求出边界热流向量；后来n z a b a r a s 等【1 9 l 用边界元法求解了二维反热传导问 题。 以前的反热传导问题主要偏重于航空航天工程中，且研究的多是较简单的 问题1 20 1 ，而边界元法有降维、高精度、适解性强等优点，因此边界元法在反问 题中的应用及研究虽然时间不长，但发展很快，经过十几年的发展。取得了很 大的进展，文献也比较多。d l e s n i c 等【2 1 1 采用边界元法对边界进行离散，并用 最小二乘法及最小能量原理分别进行了求解了一维线性非稳态热传导反问题； n s m e r a 等【”1 用边界元法求解了各向异性常规结构的热传导率问题：s o m c h a r t c h a n t a s i r i w a n 【2 3 】研究了二维及三维反热传导问题；f d e l v a r e 等1 2 4 l 以及d l e s n j c 4 等【2 s 1 均采用了一种迭代的边界元方法求解了c a u c h y 型位势反问题；k h a y a s h i 等【2 6 】【2 7 】对c a u c h y 边值问题采用了一种非迭代的直接方法进行求解a 总体来说， 边界元法在位势反问题中的研究已取得了一定的成果，但是由于反问题的不适 定性，所研究的问题大部分都是各向同性问题或一维、二维常规结构，对于正 交各向异性问题、各向异性问题以及高维复杂的问题( 如薄体问题) 研究依然 不多，需要作进一步研究工作。 1 4 本文的工作背景及主要内容 1 4 1本文工作背景 1 、几乎奇异积分 边界元法的一个重要优势就是问题的维数降低一维，但由于边界积分方程 积分核函数含有奇性基本解，使得对常规边界元法，降维不得不付出处理奇异 积分的代价。然而常规边界积分方程的奇异特性是边界元法满足良好数值精度 的保证。当源点不在边界，但趋近于边界；或源点虽在边界上，但不在当前积 分单元上而趋近于当前积分单元时，这两种情况源点与场点的距离趋近于零， 但不等于零，这时即产生所谓的几乎奇异积分( n e a t l ys i n g u l 甜i n t e g t a l ， n e a r - s i n g u l a ri n t e g r a l ，q u a s i s i n g u l a ri n t e g r a l ) ，几乎奇异积分的第一种情况表现 在求近边界内点的物理参量的场合，反映了边界元法的“边界层效应”( b o u n d a r v l a y e r e f 艳c t ) 【2 8 ，2 9 ，30 1 ，第二种情况在边界单元疏密不均时易于发生，尤其表现 在薄体结构问题上。一般问题物理场不奇异，几乎奇异积分的积分核在数学上 也是有限的，即数学上也不奇异，几乎奇异积分计算的失真本质是由数值积分 引起的，由于积分函数在此积分区间上数值变化很大，所以常规高斯积分不能 有效计算肛”。文【2 9 】对几乎奇异积分研究作了简要综述。 2 、半解析及解析的正则化算法 牛忠荣等【2 8 r 3 1 】对二维边界元法中的几乎奇异积分提出了一种通用的半解 析正则化算法，该算法适于线性元，采用分部积分，直接处理常规边界积分方 程离散积分后的积分核函数，将引起几乎奇异的积分元素变换到积分号之外， 剩下的积分为几乎弱奇异积分，采用常规高斯积分计算。周焕林等【2 9 】沿此思想， 处理了二维常规和薄体结构位势和热弹性力学问题中的几乎奇异积分，文 3 2 ，3 3 】将此方法用于求解二维弹性力学薄壁结构。文 2 8 ，2 9 ，3 2 ，3 3 将此算法用于 二次元，只要把与源点邻近的二次元分作两段线性元，就可应用上述半解析积 分算法。 文【3 1 3 3 】中半解析积分算法没处理几乎弱奇异积分，当源点进一步趋近边 界，几乎弱奇异积分的影响逐步显示出来影响结果的精确度，牛忠荣等【3 4 1 对此 进一步推导，从而得到完全解析的积分算法，周焕林等【2 9 3 5 ”】导出了解析积分 算法的另一种形式，并将此解析积分算法用于位势和热弹性力学常规结构和薄 体结构问题1 2 9 ，叫，取得很好效果。 类似二维问题处理的恩想，牛忠荣等【3 7 】对三维问题几乎奇异积分提出了一 种通用的半解析正则化算法，采用分部积分，将线性三角形单元上的几乎奇异 面积分变换为沿单元边界围道上的一系列的线积分，对此线积分采用常规高斯 积分计算，并成功分析弹性力学三维薄体结构问题。周焕林1 2 9 】将此算法推广应 用于分析三维位势和热弹性力学常规结构和薄体结构问题。 此外牛忠荣等【3 8 l 还导出了一种新的位势梯度边界积分方程一一位势问题 自然边界积分方程。奇异性比常规边界积分方程降低一阶，结合几乎奇异积分 的解析正则化方法，可完全消去积分几乎奇异特性的影响。 1 4 2 本文主要内容 周焕林在文【2 9 】的总结展望中率先提出将几乎奇异积分正则化算法推广用 于反问题边界元法，此后作了大量的前期调研和预研工作，且撰写并提交了 2 0 0 5 年国家自然科学基金申请书。提出采用迭代法、最小二乘法和奇异值分解 技术对反问题进行正则化，研究加快收敛，减少迭代的松弛措施。本文在此思 想之下开展工作，结合奇异值分解算法( s i n g u l a rv a l u ed e c o m p o s i t i o n 以下简 称s v d ) 刈和几乎奇异积分正则化算法，对二维位势反问题进行了研究。对薄 体结构问题成功实现了几乎奇异积分和不适定反问题的双重正则化。 本文主要内容可概括为以下几个方面： 1 ) 引入s v d 方法，对s v d 法的基本原理作了阐述。 2 ) 引用二维位势问题及几乎奇异积分的基本理论，对二维各向同性常规 结构及薄体结构位势反问题进行了分析，并对反问题实施了正则化，实现了几 乎奇异积分和不适定反问题的双重正则化，同时给出了算例及结果。 3 ) 分析了二维正交各向异性常规结构及薄体结构位势反问题，并对几乎 奇异积分及反问题实施了正则化过程，给出了算例及结果。 最后。对全文进行了总结。 6 第二章奇异值分解方法( s i n g u l a r v a l u e d e c o m p o s i t i o n ) 本文在求解位势反问题时。边界元法是依据部分边界条件来求解其余未知 的边界物理量。此过程中遇到奇异的代数矩阵求解问题，本章从一般情况下叙 述矩阵的奇异值分解方法。 2 1矩阵的奇异值分解 在整个复数空间内，有如f 定理： 定理l ：设4 c ”，则存在酉矩阵p ，q ，使得 - p ”4 q = i 宇3 i 这里d = 硪昭( d i ，d 2 ，d ，) ，且吐d 2 d ， o ， d ，( f _ 1 ，2 ，r ) 称为4 的奇异值，而 4 = p l 宇3 l q ” 称为矩阵的奇异值分解式h 叭。其中矩阵p 盯，q8 为矩阵p ，q 的共轭转置矩阵。 证明略，可参见文【4 0 】。 虽然上面的定理是在复数空间内证明的，值它也可以推广到实矩阵，且当a 为实矩阵时，我们也可以证明对任一实矩阵4 且”“，则存在单位正交矩阵 p ，q ，使得 p 4 q = i 宇3 i 这里d = 硪裙( d i ，d 2 ，d ，) ，且d i d 2 矗， o ， 一( j _ l ，2 ，) 称为4 的奇异值，而 一= p 孑：i q 7 ( 2 1 ) 即为实矩阵4 的奇异值分解式。 其中矩阵p ，q7 为矩阵p ，q 的转置，且p ，q 均为实矩阵。 证明：设曰= 4 7 一，则占7 = 一7 爿) 7 = 爿7 7 ) = 4 7 一= 丑，即丑为h 阶实对称矩阵， 所以矩阵丑的特征值及特征向量均为实数及实向量，以下的证明与原定理类似， 参见文【4 0 】，故略。 t 7 2 2矩阵的奇异值分解的另一种形式 对于式( 2 1 ) ，一且，p 矗”，q 皿。当聊h 且4 的列向量满秩时， 则矩阵4 的奇异值分解式为 止p 阁q 7 ( 2 2 ) 其中d 为一n n 阶的对角阵，且对角线上元素均大于零，若爿的列向量不是满 秩的，则我们也可以将矩阵一的奇异值分解式表达成( 2 2 ) 式，只是在对角阵d 中的对角线上的元素含有零。 在式( 2 - 2 ) 中，p 仍为一埘m 阶的方阵，设p = ，p 2 ) ，其中e 为m 行阶， 罡为川伽一仃) 阶，则有 p ( 言) = ，只) ( 言) = d + 只o ) = 只d 令，= 置，= d = w w l ，矿= q 因此可将( 2 2 ) 式改写为如下形式： ( 2 3 ) 所以任意给一个聊n 阶的矩阵一，其中肌”，可以分解成( 2 3 ) 式，其中矩阵u 是m 行阶的，且列向量正交，矩阵矿是h h 阶的对角阵，且w ，o ，f - l ，2 ，h ， 矩阵y 是竹胛阶的单位正交矩阵，即 莩玩也 ；i ；主： ( 2 - 3 - a ) i = l 一。一。 亦即 7 = 6 。 ，= l l i 仃 1 ，珂 ” = - ( 2 3 - b ) ( 2 3 一c ) 由于矿是单位正交矩阵，所以也有y 矿7 = 1 。 s v d 方法对于m n 的情况也同样适用。在这种情况下，对于矩阵4 的奇 异值有h = 0 ，当_ ，= 卅+ l ， 时，并且矩阵u 相应的列向量也均为零，同时等 式( 2 3 - a ) 也仅适用于七，曼m 的情况。 2 3奇异值分解方法( s v d ) 【3 9 2 3 1对于方阵的s v d 法 如果矩阵4 是方阵，例如h 阶的，那么相应的矩阵【厂，矿和缈也是肿c n 阶的 方阵。它们的逆阵都是容易计算的：由于u ，y 是单位正交矩阵，所以它们的逆 阵就等于他们的转置；是对角阵，所以它的逆阵也是一个对角阵，只是它的 逆阵对角元素等于矿中对应元素的倒数。从( 2 3 ) 式则立即可以得到下式 一= y 【威口g ( 1 1 0 ) 】- u 7 ( 2 4 ) 使上式不能成立的情况是对角阵中的对角元素w ，中有一个或几个为零，或者 虽不为零但数值太小几乎接近于零，以至受舍入误差影响太大，因此不能计算， 或计算误差太大，也就是说矩阵爿奇异，逆阵不存在。奇异值分解方法( 以下简 称s v d ) 就能处理这种情况。 “ 通常矩阵的条件数是指矩阵特征值中最大的数与最小的数之间的比值，如 上面的矩阵4 ，它的条件数即w ，中最大的数与最小的数之间的比值。如果一个 矩阵4 的条件数太大或趋于无穷，也就是说其中某几个w 的倒数超出了计算机 的计算精度，我们就说矩阵一是奇异的( 例如对于单精度w 一般不应小于l o 一， 双精度w ，一般不应小于1 0 。2 ，当然，不同计算机的精度一般不同) 。 对于奇异的方阵，零空间( n u l l s p a c e ) 和值域的概念很重要。先来看我们所熟 悉的联立方程组 4 x = 6 ( 2 5 ) 其中4 为一方阵，西和x 为列向量。方程( 2 - 5 ) 定义了一个从向量x 到向量6 的线 性变换a 。如果a 是奇异的，则存在一x 的子空间，叫做零空间，满足映射： 4 - x = o 。零空间的维数被称为爿的退化阶数。 同时也存在一个有关6 的子空间，对于该空间中任非零向量6 ，存在某 个向量x ，使得4 x = 6 成立，这个子空间亦称做4 的值域。值域的维数也就是 矩阵爿的秩。如果4 不奇异，则一的秩为h ，如果4 奇异，则爿的秩小于h 。 事实上零空间的维数加上爿的秩等于n 。 特别应当指出的是在s v d 法中，矩阵u 相应于w o 的列向量是矩阵a 的 列向量的一组基，而矩阵y 相应于w ，= o 的列向量是零空间的一组基。 现在我们来看看方程组( 2 - 5 ) ，如果6 = o ，则任意一个相应于w = 0 的矿的 9 列向量均为原方程的一个解。 证明：设 y = “，v ：，) 则有 y 7 = 再设w = o ，则x = ”必为原方程爿工= 0 的解，因为y 为单位正交矩阵，所以 有 一x = 4 - q 7 = u = u = u 0 ： ， o ，” 1 ，“ r v l r v 2 ： r v n 0 = 0 证毕。 如果6 0 ，首先的问题是原方程有没有解，如果向量髓能由4 的列向量线 性组合，亦即原方程有解，设x 为原方程的解，显然有对相应于w = 0 的y 中 的列向量”经过任意的线性组合加上x 仍为原方程的解，即原方程有无穷多解。 这时一般想得到一个令h 2 最小的解，我们可以在w ，= o 时用零取代1 w ，然后 采用下式计算 x = y 【硪昭( 1 7 u ) 7 6 )( 2 6 ) 上式所求x 即为令h 2 最小的解。 证明：设4 x = o ，4 x = 6 ，x 由( 2 6 ) 式给出，则有 4 r o + x ) = 6 考虑l x + x i ， 卜+ x l = p - ，7 6 + x i 1 0 = l 矿妒u 7 西+ 矿7 x 】 = 阿- 7 6 + y 7 l 考虑到矿一- u - 6 仅当w ，o 时有非零元素， 素，因此有 ( 2 7 ) 而矿- x 仅当= 0 时有非零元 l 工+ x i l x i 证毕。 如果在方程组( 2 5 ) 中，6 不能由4 的列向量经线性组合得到，也就是说原 方程无解，这种情况在力学当中经常遇到，这通常是由于数值计算引起的误差 所致，这时我们仍然可以由( 2 6 ) 式构造一个原方程的解，虽然它并不是原方程 的精确解，但是在所有可能的向量工当中，它是使阻x 一6 l 取得最小值的解。 证明；假设我们对上面所求的x 用任选的某一x 来进行修正，即用x + 工取代x ， 则月，x 一6 将被加上一个量= 爿x ，显然6 能由4 的列向量线性组合，我们可 以得到 阻x 一6 + 6 i = i ( u y 7 ) ( 矿一u 7 6 ) 一西+ 西i = l ( u 矽u 1 一1 ) 6 + 6 l = j e ，- ( f y 孵一一1 ) u 6 + ，7 6 1( 2 8 ) = l ( 形一一1 ) - u 7 6 + u 7 圳 在上式中最后一行，缈- 形一一1 仅当w ，= o 时有非零元素，而u 7 6 仅当w 0 时 有非零元素，因此我们有( 2 - 6 ) 式所取x 使阻工一6 i 取得最小值，证毕。 2 ，3 2 对于非方阵的情况 对于非方阵的情形与方阵的证明类似，x 也可由( 2 6 ) 式给出，只是u 7 不再 是一个方阵，而是一个n m 阶的矩阵。 在我们对位势偏微分方程进行求解时，首先通过插值建立联立方程组，对 于位势反问题方程个数一般大于未知数，普通的高斯消除法不能解决这一类的 问题，而s v d 方法不受方程个数及未知量个数的限制。我们知道反问题一般是 不适定的，无论对原方程有无穷多解的不定解问题的情况还是对原方程无解的 超定解问题的情况，用s v d 方法都能给出一个方程的解，尤其对超定解问题， 用s v d 方法一般能够求出满足要求的近似解的。 使用s v d 法有个值得注意的地方就是如果w ，太小，我们在求未知量工 时就用零取代l w ，进行计算，但究竟多小，本文不作详细讨论，可以参阅文f 3 9 。 第三章二维各向同性位势反问题 在实际工程中，由于空间和环境条件等的限制，很难或不能测量物体某 部分边界上的物理量，而只能通过物体的其他部分边界上进行测量得到的已知 量来确定，这种情况在许多物体的热传导等研究中会遇到，比如在航空航天领 域、智能材料、涂层结构以及压力容器等等，因此对此类反问题的研究是很有 广泛意义的。 3 1 各向同性材誊斗常规结构位势反问题 3 1 1 位势问题的边界积分方程 对各向