（无线电物理专业论文）固体介质中声波散射的有限元数值模拟.pdf

固体介质中声波散射的有限元数值模拟 李树榜( 无线电物理) 指导教师：李书光教授 摘要 超声波入射到固体中裂纹、孔洞等缺陷上时将发生散射，散射波 中包含了缺陷的特征信息，了解缺陷和超声波的相互作用对缺陷的检测 和定征起着关键性的作用。现有的解析方法只适用于几何形状简单的散 射体，由于实际散射体形状各异，除十分简单的情形外，传统的解析方 法对它们难以奏效，只能借助于数值方法。本文建立了各向同性固体的 具有吸收边界的二维平面应变有限元模型，在研究了脉冲超声波在无缺 陷固体中传播的基础上，首先模拟计算了脉冲超声点源加载时，超声波 在固体内部平行于表面的裂纹、垂直于表面裂纹、与表面成4 5 0 角斜裂 纹、四分之一圆周的弧形裂纹和背表面裂纹等5 种形式裂纹上的散射。 然后分别研究了平面脉冲超声纵波和横波入射时，在固体中横向圆柱形 孔上的散射。通过对纵、横波的分离和在全波场中对散射波的分离，清 晰地显示了缺陷对脉冲超声波的散射和波型问复杂的模式转换。模拟结 果以位移场、散度场( 纵波) 、旋度场( 横波) 波场快照和a 扫描曲线 形式给出，部分结果与现有动态光弹实验照片进行了比较，并做出了评 价。 关键词：固体介质，声波，散射，有限元方法，数值模拟 f i n i t ee l e m e n ts i m u l a t i o no fu l t r a s o n i c 矸白y e s c a t t e r i n gi ns o l i d l is h u b a n g ( r a d i o p h y s i c s ) d i r e c t e db yp r o f e s s o rl is h u - g a n g a b s t r a c t d e f e c t ss u c ha sc r a c k sa n dh o l e sa c ta ss o u r c e so fw a v es c a t t e r i n g w h e ni l l u m i n a t e db ya r ti n c i d e n tu l t r a s o n i cp u l s e t h e s c a t t e r i n gw a v e s c o n t a i na l lt h en e c e s s a r yi n f o r m a t i o no ft h ef l a w s u n d e r s t a n d i n gt h e u l t r a s o n i cw a v e f l a wi n t e r a c t i o ni sa ne s s e n t i a lf e a t u r ei nd e t e c t i n ga n d c h a r a c t e r i z i n gt h ed e f e c t s m o s to ft h ea v a i l a b l ea n a l y t i c a lm e t h o d sf o r s c a t t e r i n gp r o b l e me o l l f j m e t ot h ef l a w sw i t hs a m p l eg e o m e t r i e su n d e r i d e a l i z e dc o n d i t i o n s h e n c e ，n u m e r i c a lt e c h n i q u e sa r ec a l l e di n t op l a yf o r t h ec o m p u t a t i o no fs c a t t e r e dw a v ef i e l d s i nt h et h e t i s ，t w o - d i m e n s i o n a l p l a i ns t r a i nf i n i t ee l e m e n tm o d e l sw i t ha b s o r b i n gb o u n d a r yc o n d i t i o n sh a v e b e e n d e v e l o p e dt oi n v e s t i g a t et h e u l t r a s o n i cw a v ep r o p a g a t i o na n d s c a t t e r i n gi i li s o t r o p i cs o l i dm e d i a f i r s t l y , s c a t t e r i n gb yf i v ek i n d so fc r a c k s h a v eb e e ns t u d i e da sf o l l o w s ：p a r a l l e lt ot h es u r f a c e p e r p e n d i c u l a rt ot h e s u r f a c e , 4 5 0a n g l et ot h es u r f a c e ，q u a r t e r - c i r c l e ，b a c k - f a c ec r a c kw i t ha p u l s e d - u l 仃a s o n i c - p o i n t - s o u r c el o a d e d s e c o n d l y , t h ei n t e r a c t i o no fp u l s e d - u l t r a s o n i c - - p l a i n - c o m p r e s s i o n s h e a r - w a v ew i t hac y l i n d r i c a lh o l eh a sb e e n i n v e s t i g a t e d i no r d e rt os h o wt h es c a t t e r i n ga n dc o m p l e xm o d e c o n v e r s i o n c l e a r l y , t h et o t a lf l i e di sd e c o m p o s e di n t ol o n g i t u d i n a la n d 缸 a n s v e r s ep a r t s ， a n dt h es c a t t e r e df i e l di so b t a i n e db ys u b t r a c t i n gi n c i d e n tw a v ef r o mt h e t o t a lf i e l d s n 抡r e s u l t sa r cp r e s e n t e di n t h ef o r m o fa - s c a np l o t so r d i s p l a c e m e n t , d i v e r g e n c ea n dc u r ls n a p s h o t s ，a n d s o m eo ft h e ma r e c o m p a r e dw i t ht h ep r e v i o u se x p e r i m e n t a lo b s e r v a t i o n so b t a i n e db yt h e d y n a m i cp h o t o e l a s t i c i t ym e t h o d ，a n de v a l u a t i o n s a l ed o n e k e yw o r d s ：s o l i dm e d i a , a c o u s t i cw a v e ，s c a r e r i n g ，f i n i t ee l e m e n tm e t h o d ， n u m e r i c a ls i m u l a t i o n 独创性声明 本人声明所呈交的论文是我个人在导师指导下进行的研究工作及 取得的研究成果。尽我所知，除了文中特别加以标注和致谢的地方外， 论文中不包含其他人已经发表或撰写过的研究成果，也不包含为获得中 国石油大学或其它教育机构的学位或证书而使用过的材料。与我一同工 作的同志对本研究所做的任何贡献均已在论文中作了明确的说明并表 示了谢意。 签名；垄丝盗纱。年华月f 。日 关于论文使用授权的说明 本人完全了解中国石油大学有关保留、使用学位论文的规定，即： 学校有权保留送交论文的复印件及电子版，允许论文被查阅和借阅；学 校可以公布论文的全部或部分内容，可以采用影印、缩印或其他复制手 段保存论文。 ( 保密论文在解密后应遵守此规定) 学生签名： 导师签名： 中国石油大学( 华东) 硕士论文第1 章前言 第l 章前言 无损检测是现代工业许多领域中保证产品质量与性能、稳定生产 工艺的重要手段。而超声检测是一种重要的无损检测技术，由于它穿透 能力强、对人体无害，已较广泛应用于工业及高技术产业中。但到目前 为止，工业用超声无损检测大多还停留在了解材料与构件内是否有缺 陷，或凭经验大致判断缺陷的大小与位置的水平。近年来人们对超声 无损检测提出了更高的要求，不仅要知道缺陷的有无及其位置，还希望 知道缺陷的大小、形状、取向和内含物等，这样就有可能确定缺陷的危 害性大小和缺陷所在部件的可使用寿命。这些额外信息，只有试图从接 收到的散射波中汲取。这时，就需要首先了解不同形状、大小等障碍物 超声散射波各具有什么特征。超声散射的研究便由定性上升到定量的范 畴。而在无损检测中，声波的载体绝大多数是固体，缺陷则是固体中的 声波传播的障碍物，会引起声波的散射。于是，对固体中超声散射的研 究变得突出出来，受到了人们的极大关注，也取得了不少进展。目前， 对一些形状比较简单、有旋转对称性的散射体进行了实际计算，求算了 大量的解。然而，流体中的散射问题本已很复杂，和流体相比，固体中 声波的传播和散射规律有其独特的复杂性，且实际散射体形状各异，绝 大多数无法求得散射问题的解析解嘲。因此，一些数值分析方法( 如有 有限元、有限差分、边界元等) 在处理散射问题中得到广泛应用，出现 了大量的科技文献。特别是有限元法，经过了近半个世纪的发展，其理 论基础己趋向成熟，已在一些实际工程问题中得到运用。有限元波动模 拟理论及应用的研究有着广阔前景，在理论和实际应用中都有着重要意 义，值得进一步探讨。 声学中讨论散射现象已有很长的历史。约从上世纪4 0 年代开始， 中国石油大学( 华东) 硕士论文第1 章前言 超声在固体中的散射及应用研究得到了发展。1 9 5 6 年，y m g 和t r u e l l “ 对固体中异质球形散射体( 颗粒) 的声散射问题作了开创性的研究，使用 波函数展开法从理论上详细分析了平面入射纵波在不同边界条件下的 散射截面，并在r a y l e i g h 近似条件下，得到一系列解析表达式，为固体 中的声散射问题的研究奠定了基础。随着波动问题研究的不断深入和研 究领域的不断扩展，经历了几十年的共同努力，国内外学者己提出多种 方法用于声波散射的远近场分析。其中能求出散射场准确解的主要有波 函数展开法和积分方程法“1 。 波函数展开法的基本思想是对标量或矢量波动方程作变量分离， 利用一些曲线坐标和与之对应的特殊函数，如b e s s e l ，w e b e r , m a t h i e u 函数来构造波动方程的解，并进一步满足边界条件，求解边值问题。对 波函数展开法，只有目标表面与坐标面相一致时才能使用这种方法。而 曲线坐标的种类有限，到目前为止只有1 1 种( 直角坐标、圆柱坐标、椭 圆柱坐标、球坐标、圆锥坐标、抛物线坐标、扁长球面坐标、扁球面坐 标、椭圆球面坐标、抛物面坐标) ，且对于矢量波动方程，分离坐标的 数量仅限于前6 种。因此，该方法在使用上受到很大限制嘲。所涉及文 献主要集中在对不同性质、不同边界状态的球状。”和柱状“8 1 等几何形 状规则的散射体的研究。但作为一种经典的方法，它所提供的原则和解 答，对求解声波散射问题起到了非常重要的作用。 积分方程法一般是从弹性动力互等原理出发导出积分方程。根据 积分表达式定理，用一个面积分和一个体积分来描述由任意形状的物体 产生的散射波。对于材料性质不连续的夹杂物，可以通过等价的体积力 与周围介质g r e e n 函数的体积分来描述由不连续区域引起的散射波。由 于体积分中的未知函数含有总波场，所以我们得到的是关于散射波的积 2 中国石油大学( 华东) 硕士论文第1 章前言 分方程，其类型主要为f r e d h o l m 型积分方程和具有c a u c h y 型积分核的 奇异积分方程，可用渐近方法或迭代方法进行求解。如果第一次试函数 用入射波场代替总波场，所得的解被称为b o r n 嘲近似。如果用静态场代 替未知的总波场，所得的解被称为静态近似“。积分方程对于二维定常 波和瞬态波问题很有效。 解析方法能够求出准确解，且在理论上已相当成熟，但只适用于 几何形状简单的散射体，由于实际散射体形状各异，除十分简单的情形 外，传统的解析方法难以奏效“。 一般来说声波是看不到的。人们在研究声波时，通常用传声器、 超声换能器等进行测量。在研究固体中超声波的传播及散射时，则只能 在固体表面接收声波，得到局部信息。声波在固体内部的行为及变化过 程则要推论，而不能直接测量。为了显示固体中的声波，从上世纪7 0 年代开始使用的光弹法是目前最有效的方法。它利用超声波在固体中传 播时所引起的应力，用来显示透明固体中应力波的传播。实验观测了超 声波在裂纹和空腔等缺陷上的散射“”，取得了一批可喜的成果。但现 有的光弹法亦有其很大的局限性，一是光弹法只能显示二维声场；二是 显示声波的灵敏度不高，不易观察某些散射细节。另外，由于在固体中 制作多种多样的散射体不但在人力物力上难以承受，而且相应的理论难 度也很大，所以目前只限于几何形状比较规则的散射体。 近年来，一种返朴归真的研究方法愈来愈受到国内外研究者的重 视。这就是采用有限差分、有限元方法将连续问题转变为离散问题，并 运用计算机直接进行波动数值模拟“”。这是因为，有限区域内的介质无 论在几何上和力学性质上多么复杂，其中的波动在原则上可以用时空离 散技术加以模拟。波动计算机模拟可视为研究近场波动问题的一种具有 中国石油大学( 华东) 硕士论文第1 章前言 普适性的试验方法。这种数值试验类似于原型和物理模型试验，但是， 在模拟复杂近场波动的能力上，在模型“制作”，参数选取和变动的灵 活性上，较之实物模型试验皆具有无比的优越性。同时，取得准确而完 整的数据无需任何“测量”设备，并可灵活而直观地显示试验结果。此 方法己应用于很多学科领域。这些应用研究表明，在许多情况下研究者 能够判断数值模拟结果的合理性，大大推进了研究工作。研究声波散射 常用的数值模拟方法主要有有限差分、边界元和有限元等，这些数值方 法皆在工程中得到广泛应用，已经取得了大量有使用价值的研究成果。 有限差分方法是将偏微分方程化为有限差分方程的方法。它首先 把微分方程的定解区域剖分成由有限个格点组成的网格，然后用各离散 点上待求函数的差商来近似代替该点的微商，从而把需要求解的微分方 程转化为一组差分方程的线性代数方程组，并用该差分方程所得的各离 散点上的待求函数的数值解来近似原来微分方程的连续解。它的优点是 概念和计算格式比较简单“，对于求解具有规则几何形状和均匀材料特 性的问题特别有效，且程序设计简单，收敛性也较好。不过，由于它局 限于规则的差分网格，如正方形、矩形或正三角形网格等，几何上复杂 的区域难于准确表达，不适于非均匀的网格，且沿非线性边界难以引入 边界条件，使它的适应性受到极大影响。 对于边界元方法，1 9 6 3 年，j a s w o n 和s y m m “”首次发表了二维 位势问题边界积分方程解法。他们把边界离散成许多小段( 单元) ，把每 个小段内的边界变量视为常值，采用配点法将边界方程化为一组代数方 程。方程的系数，对于非奇异积分采用s i m p s o n 积分法，对于奇异积分 采用解析积分，这便是现代边界元中的常元法。 声学中最早采用边界元方法研究声辐射的是c h e n 和s c h w e i k e r t 4 中国石油大学( 华东) 硕十论文第1 章前言 ( 1 9 6 3 ) “”。由于边界元方法在解决无限域问题上显示出巨大的优势，它 受到了美国海军的大力支持。到上世纪六十年代末七十年代初，边界元 方法已广泛应用于声辐射和散射的研究中，并取得了极大的发展。在以 后的几十年中，声学中的边界元方法得到了进一步发展，直至今日它仍 是计算声学中的研究主题之一。运用边界元方法，k i t a h a r a ( 1 9 8 9 ) “8 1 在 频域中讨论了致密非均匀体对弹性波的散射；h s a t o 和y s h i n d o ( 2 0 0 2 ) “”研究了单个多边形夹杂对平面纵波和横波的散射，并将结果推广到多 个多边形夹杂的情况。 尽管边界元方法是一种强大的数值计算工具，但还存在很多不足。 它事先无法估计计算精度，难以考虑如何通过自动加密来修正误差：网 格重新划分极为费时，重新划分后网格后的边界元系数都必须重新计 算；计算结果的可靠性往往取决于人的理论素质和工程经验；现有的边 界元数值方法的求解效率、精度和求解工程问题的规模都有待提高。另 外，边界元的系数矩阵是非对称满阵，对大型工程问题，边界元将丧失 其优点，所需存储量和求解所需时间可能比有限元方法还多1 。 有限元方法的思想起源于上世纪4 0 年代一些力学专家对飞机结构 的设计。到上世纪6 0 年代，w c l o u g h 进行了系统总结而成为一种比较 完整的方法，并将其命名为“有限元法”。而直到1 9 7 2 年，才由b a b u s k a 和a z i z 全面系统地揭示了偏微分方程的现代数学理论和有限元法的内部 联系，从而为该方法奠定了严格的数学基础。 由于有限元法广泛适用于各类偏微分方程的边值问题，并特别适合 于计算机编程实现，因此自有限元方法问世以来，人们很快就认识到它 的潜力，特别是在工程技术领域。随着计算机技术的迅速发展，有限元 方法的应用也得到了迅猛地发展，并日趋完善。如今，有限元方法的应 中国石油大学( 华东) 硕士论文第1 章前言 用已从固体力学领域扩展到流体力学、热传导、电磁场、地球物理和声 学等领域。 有限元方法求解基本方程( 常微或偏微分方程) 和相应定解条件数 值解的基本思想是将连续的求解域离散为有限个、并按一定方式相互连 接在一起的单元组合体。在每个单元内，假设近似函数，来分片地表示 整个求解域上的未知场函数。由于在联接相邻单元的结点上，场函数应 具有相同的数值，因而将它们用作数值求解的基本未知量。这样一来， 求解原来待求场函数的无穷多自由度问题转换成求解场函数结点值的 有限自由度问题o ”。 有限元法吸取了有限差分法中离散化处理的内核，又继承了变分 计算中选择插值函数关于区域积分的合理方法，使它具有很大的灵活性 和通用性。t e m p l e ( 1 9 8 8 ) 。4 和b o n d ( 1 9 9 0 ) 详细对比了现有的各种数 值方法，得出的结论是：因有限元方法可以方便灵活地模拟不同类型基 材中任意形状的缺陷而优于其它方法。对于固体中波动数值模拟，有限 元方法的主要优点有：( 1 ) 它的结点配置比较任意，对于具有复杂形状 的散射体，可以使边界结点完全落在边界上，在边界上给出较好的逼近； 对于由几种材料组合而成的系统，可以通过把结点取在分界面上而得到 较好的处理。( 2 ) 在有限元分析中，几何边界条件可以通过在装配一级 上对离散方程简单的修改而很容易给定，自然边界条件和不同介质分界 面则无需考虑。( 3 ) 有限元法用统一的观点对区域内的结点及边界结点 列出计算格式，使各结点的精度在总体上比较协调。最后的有限元离散 矩阵具有对称正定的性质，给问题的求解带来方便。( 4 ) 有限元方法可 以得到连续的近似，任一点场变量及其导数的计算简单方便。 一 有限元近场波动数值模拟的基本问题可以概括为：近场区域内波 6 中国石油大学( 华东) 硕士论文第l 章前言 动的离散模拟和近场区域边界的处理。对前者，有限元方法中通常采用 三种离散模型，即集中质量、一致质量和混合质量模型。廖振鹏 ( 1 9 9 2 ) 叫详细对比三种模型，得出的结论是，从实用观点看，集中质量 模型是波动模拟的最佳选择。对于后者，从上世纪6 0 年代末以来，国内 外进行了大量关于人工边界条件的研究工作。如s m i t h ( 1 9 7 4 ) u s j 消除界 面反射的叠加技术，取得了很好的效果，但要大大加重计算量。l y s m e r 和k u h l e m e y e r ( 1 9 6 9 ) ”1 提出了粘性阻尼边界，可以衰减掉大部分的反 射。我国学者廖振鹏“”对人工边界条件的建立、稳定性等也做了大量 研究。直到今天，对于人工边界条件的研究依然没有停止，仍有新的边 界条件的提出啪。总之，各种人工边界条件都有一定的适用范围，依据 不同情况而选取。 运用有限元方法，国外许多学者对固体中声波散射进行了模拟。 k m a j a l e e l ( 1 9 9 3 ) 嘲计算了脉冲超声波在固体中裂纹及双层复合材料 界面裂纹上的散射，结果以a 扫描曲线和波场快照的形式给出。d d a t t a ( 1 9 9 6 ) ”1 研究了超声波在各向同性和正交各向异性固体介质中的传播 过程，模拟了实验的超声回波技术，获取了a 扫描数据并对其在频域 做了特征分析。n n k i s h o r e ( 2 0 0 0 ) 啪1 计算模拟了各向同性固体中裂纹和 柱形孔等孤立缺陷对超声波的散射，得到散射波能量的角分布。t i t o ( 2 0 0 0 ) 。”研究了脉冲超声波在固体中柱形夹杂上的散射。 国内早在上世纪八十年代初，有限元方法开始在工程声学中得到 应用嘲。例如，分析不规则房间的声场，分析电磁换能器的振动问题和 评价抗性消声瓦的特征等。但此后一直发展缓慢，虽然人们越来越认识 到声波散射数值模拟对于超声检测的重要性溉”，但目前我国计算机技 术在超声检测中的应用主要集中在超声信号的采集、量化和处理及超声 7 中国石油大学( 华东) 硕士论文第1 章前言 成像系统上，在超声波散射数值模拟方面的研究进行的还很少，是一 个极其薄弱的环节，急待加强。 总之，运用有限元方法，许多国外学者对固体中声波的散射数值 模拟研究已经取得了若干可喜成果，有些已经在生产实际中得到应用， 但成果主要集中在较为规则的散射体上，还有很多问题没有解决。我国 声波散射有限元数值模拟的起步较晚，研究成果还很少。本文的主要工 作是建立固体介质中任意形状散射体的二维有限元模型，编写通用有限 元计算程序，研究声波在固体中裂纹和空腔等缺陷上的散射。 全文内容安排如下： 第二章从弹性力学的基本方程入手，导出各向同性弹性固体中波 动模拟的有限元运动方程；对计算中使用的线性三角形单元和双线性四 边形单元做单元特性分析：讨论并选定有限元人工边界条件；给出中心 差分法和n e w m a r k 方法两种积分格式。 第三章介绍所建立的二维有限元模型和计算程序，对比相同计算 条件下中心差分法和n e w m a r k 方法的计算结果，选定使用n e w m a r k 方 法；模拟计算超声脉冲点源激发声波在固体板中的传播，验证所建立的 有限元模型和编写的程序的有效性。 第四章模拟计算板内部平行于表面的不同长度直裂纹、垂直于表 面的裂纹、与表面成4 5 。角斜裂纹、四分之一圆周弧形裂纹和背表面裂 纹对超声脉冲点源激发声波的散射；讨论裂纹长度对超声脉冲回波、透 射波幅度的影响；分析各裂纹散射声场的特征。 第五章研究脉冲超声平面纵波和横波在固体中圆柱形孔上的散 射，计算结果与现有动态光弹实验照片进行对比。 第六章总结全文，展望以后可以开展的研究方向和研究内容。 8 中国石油大学( 华东) 硕士论文第2 章弹性固体中声波波动问题的有限元法 第2 章弹性固体中声波波动问题的有限元法 2 1 弹- 性力学基础和固体中声波方程 2 1 1 弹性力学基本方程幢1 o x ，q ，o r z ，勺，f ，吃来表示。其中吒，q ，吒为正应力；，吃为剪切 盯= 吒q 巳吃r ( 2 1 1 ) 甜= 习= 缸v w r g - 2 ， 弹性体内任意一点的应变，可以由6 个分量t ，o ，乞， 来表示。其中毛，勺，巳为正应变；，如为切应变。表示为矩阵形式 。= ko 乞岛比l r ( 2 1 3 ) 9 中国石油大学( 华东) 硕士论文第2 章弹性固体中声波波动问题的有限元法 抛 巳2 面 却 占，2 万 却却 2 瓦+ 万。岛 矩阵形式为 占= 三馨 其中上为矩阵算子，即 l = 跏 占= 。 瑟 抛a v 如2 万+ 瓦2 孙0 w 。西+ 瓦2 k ( 2 1 4 ) ( 2 1 5 ) ( 2 1 6 ) ( 2 ) 物理方程应力一应变关系 对于各向同性的线弹性材料，应力与应变关系的矩阵表达式为 盯= 觑 ( 2 1 。7 ) 其中d 称为弹性矩阵 d ：墨! ! 二尘 ( 1 + u ) ( 1 - 2 0 ) l 旦旦 l 一，l 一， 1 j l 1 一u l 对 称 oo 0 0 l 一2 0 2 ( 1 一们 1 0 o 0 0 1 2 0 2 0 一d 1 o o o o 0 1 2 0 2 ( 1 一们 ( 2 1 8 ) a一如o a 一良 o a一七a一钞 a锣a一缸o o o a 一瑟 o a一钞o a良o o 中国石油大学( 华东) 硕士论文第2 章弹性固体中声波波动问题的有限元法 弹性体v 的一部分边界上外力q ，、q ，、q ：已知，称为力的边界 f z = n 一。+ n ，口+ n r ： q 1 1 0 0 其中珂，、n ，、也为边界外法线的方向余弦。式( 2 1 1 0 ) 的矩阵形式为 以= 以曼台n y 0 耋蔓0 q11100 ， 以= i b取 l ( 2 ) 1 o b 仇1 只= q ，0 = 9 ，只= q ： ( 2 1 1 3 ) “= 石，v ；，w = 石 ( 2 1 1 4 ) 主里互迪盔堂l 兰垄! 堡主丝苎 箜! 童堂丝旦堡主主婆婆垫塑里堕互堡! 旦墨 弹性介质在应力应变状态下运动时，其运动应力平衡微分方程为 等+ 拿o y + 誓+ q = p 窘m 睨“ 等等+ - 甜- g ”- + 岛= p 矿0 2 v 鲁+ 鲁+ 誓+ q = p 窘 耶， 呶咖 o z 0 i 其中，q ，g ，q 为单位体积的体积力在善，) ，z 方向的分量。将( 2 1 1 5 ) 式中未知函数统一用位移函数表示，整理简化后可以得到各向同性弹性 固体介质运动平衡的位移方程 训瓦o o + 2 “+ q = p 争删a ( 砷) 石o o + 肺+ g = p 争洲 ( 撕) 西o o + 2 w 也= p 窘 ( 2 1 1 6 ) 眩 “ 其中旯和为l 砌e 常数= 导+ 导+ 等，口= v 一罢+ 多+ 尝 为体应变。 式( 2 1 1 6 ) 可写为如下矢量形式 ( 旯荆v v 科胛2 + q = p 窘 ( 2 1 1 7 ) 根据场论中有关拉普拉斯运算的恒等式 可2 一v 目n v x v 越 ( 2 1 1 鼬 中国石油大学( 华东) 硕士论文第2 章弹性固体中声波波动问题的有限元法 式( 2 1 1 7 ) 可以表不为 ( 2 + 2 t ) v v u + t v x v x u + q = p 窘( 2 1 1 9 ) 式( 2 1 1 9 ) q h 左边第一项对应于纵波，第二项对应于横波，第三项为体 力项，它表示震动源。 进一步对上式中的位移矢量邸和体力q 分别做亥姆霍兹 ( h e l m h o l t z ) 分解，令 以= v 矿+ v 妒，q - - v o , + v x q 。( 2 1 2 0 ) 其中伊和y 称为位移位，缈为标量位，y 为矢量位，分别对应于纵波和 横波；q 和q ，为体力q 的标量位和矢量位，分别对应于纵波震源和横 波震源，且满足 v x ( v 力= 0 ，v ( v 妒) = 0 v ( v g ) = 0 ，v ( v 岛) = o ( 2 1 2 1 ) 可以得到各向同性弹性固体介质中的纵、横波方程为 俨缈一岿一专g v 2 缈专鲁= 一虿1g (2122) 辄，= j 孕， ，l = 居删表示黼黼传播航 对k ( 2 1 2 0 ) 中左式两侧分别取散度和旋度并利用( 2 1 2 1 ) 式，可得 v = v ( v 妒) + v ( v 妒) = v 2 妒 中国石油大学( 华东) 硕士论文第2 章弹性同体中声波波动问题的有限元法 v 却= v x ( v 伊) + v x ( v x y ) = v ( v 妒)( 2 1 2 3 ) 由上式可见，位移脚矢量场的散度场和旋度场分别为只与纵波和横波有 关的量。如果求得位移矢量场，分别求其散度场和旋度场，即可将总 波场分离为纵、横波分量。 2 2 弹- 陛固体中声波波动问题的有限元法 2 2 1 动力学有限元方程。蚓 固体中声波的传播为动力学问题，固体中质点位移、速度、应力、 应变和载荷等物理量均与时间有关。运用有限元方法求解此类问题时， 处理的是四维 ，y ，：，f ) 问题。这里采用部分离散的方法，即只对空间域 进行离散，将区域划分为有限个单元。 由于只对空间域进行离散，在单元口内选择位移插值函数u ( x ，y ，z ，t ) u ( x ，y ，z ，f ) = l v ( x ，y ，z ) a 扣( t )( 2 2 1 ) i “( x ，y ，z ，f ) l 其中髓 y , z ，t ) = iv ( x , y ，z ，f ) l ，n ( x , y ，z ) 为单元插值函数矩阵，4 扣( f ) 【以薯y ，z ，r ) j 为与时间有关的单元结点位移向量。 设西8 ( f ) ，矗8 ( f ) 分别为单元结点速度和加速度向量， 示为 i k x , y ，：，f ) = n ( x , y , z ) a 9 ( f ) j i ( y ，z ，f ) = n ( x , y ，z ) 矗8 ( f ) 将( 2 2 1 ) 式代入( 2 1 5 ) 式，可得 1 4 分别可以表 ( 2 2 2 ) ( 2 2 1 3 ) 中国石油大学( 华东) 硕士论文第2 章弹性固体中声波波动问题的有限元法 s = l n ( x ，y , z ) a ( 。( f ) 令b = l n ( x , y ，z ) ，称为应变矩阵，则可将应变应力表示为 暑= b a 。( f ) 盯= d e = d b a 。( f ) ( 2 2 5 ) ( 2 2 6 ) 有限元运动方程可应用拉格朗e l ( l a g r a n g e ) 方程或哈密顿( h a m i l t o n ) 原理得到。这里应用拉格朗日方程，拉格朗日方程为 其中 旦丝一丝：o 席甜加 ( 2 2 7 ) 三= t 一 ( 2 2 8 ) 上式中上为拉格朗日函数，r 为动能，万，为势能。在单元内可分别表 示为 r = 百1 肛7 谢矿 一 万争吉肛7 耐矿一肛7 缈一肛7 酗 - r ”r t ”舒” ( 2 2 9 ) ( 2 2 1 0 ) 其中矿。为单元体积，p 为密度。将式( 2 2 1 ) - - ( 2 2 6 ) 分别代入( 2 2 9 ) 和 ( 2 2 1 0 ) 5 戈，可得r 和为 r = ep = 寻饥壹肛r n d n d v i zr = p = 寺五7 【肛7 ( 2 2 1 1 ) 中国石油大学( 华东) 硕士论文第2 章弹性固体中声波波动问题的有限元法 = 萋e 石，= i 14 7 【善e 胪7 d b d 矿k 一口r 【窆j r 烈矿 。1 。“” 8 1 ”“ ( 2 2 1 2 ) + 窆肛7 酗】- a r p 。( t ) 其中只( ，) 为结点集中力向量。 定义 m 扣) _ 肛7 n d v ，k 扣= 胪7 d b d v 掣= 7 酗，群虹7 q d v ( 2 2 1 3 ) 其中m 扣为单元质量矩阵，k 。为单元刚度矩阵，碟。为由单元表面力 产生的单元结点载荷向量，群曲由单元体积力产生的单元结点载荷列向 量。 将式( 2 2 1 1 ) 和( 2 2 1 2 ) 代入c 2 7 ) 式，得到计算区域整体有限元运动 m i i ( t ) + k a ( t ) = p ( f ) r 2 2 1 4 ) f 式( 2 2 1 4 ) 为关于时间f 的二阶常微分方程组。其中m = m ，为整 p = i f 体质量矩阵； 置= 置， 为整体刚度矩阵； 归l f j ( f ) = ( 。+ 碍砖) + 只o ) ，为整体结点载荷向量。 e = l 2 2 2 弹性平面问题的单元分析嘲 本文中建立的有限元模型均为二维平面模型，使用了三结点线性三 角形单元和四结点双线性四边形单元。下面对这两种单元分别予以介 ! 堂燮塑兰奎塑圭堡窒笙! 童登丝旦堡主重丝垫垫塑望盟宣堡歪鲨 绍。 ( 1 ) 线性三角形单元 在x y 平面内的三角形单元如图2 1 所示， 有i 、，七三个结点，结点顺序按逆时针排列， 每个结点有两个位移分量，共6 个结点位移， 即 口。= 陋jm u j0 唯】7 位移插值多项式取为 u ( x ，力= q + x + 嘞y 图2 i 线性三角形单元 v t x , y ) 2 a 4 + a 5 x + 1 2 6 y q 2 1 5 ) 其中嘶为广义坐标，q ，口：，鸭可以通过结点位移甜，甜j ，表示， o t 4 ，瓯用结点位移q ，_ ，h 表示，可得 嘶= 击( q + a f l t j + 吼) = 西1 ( 6 f 珥+ 屯哟+ 以蜥) 1 口3 = 互j 约+ q q + c k u k ) 嘞= 砑1 ( q b + q 。+ 咏唯) = 击( 6 i v j + 屯叶+ 玩吒) = 去 巧+ 。_ + q 雌) ( 2 2 1 6 ) 1 7 中国石油大学( 华东) 硕七论文第2 章弹性固体中声波波动问题的有限元法 ，1 1 一只 肚三1 1 _ 乃 1 1 磁以 = 1 + v k + x 幽_ x l y k x l y l x k y ) a 1 2 xj y k x k y j oj 2 x k y t x t o k 2 x i y i x j y b i = y i y k b i = y t y b k = y i y | c 2 x k xj ，c i 2 x l x k c i 。x j x tq 2 1 7 ) 将( 2 2 1 6 ) 式代入( 2 2 1 5 ) 式，可得 秘q ，协= n 1 珏i + n j + n 弘k v ( x , y 1 = n y l + n j + n i ，k 其中j ，m ，帆为插值函数，其表达式为 m = 击( q + t 工+ q y ) n j = 紊j + b x + c j 订 以= 击( 吼+ 以x + q j ，) 将式( 2 2 1 8 ) 写成矩阵形式为 口化力= n a 力 ( 2 2 1 8 ) ( 2 2 1 9 ) ( 2 2 2 0 ) 其中 删= 暖期 = 阳力m 品m _ 也也：一 1 8 ( 2 2 2 1 ) 中国石油大学( 华东) 硕士论文第2 章弹性固体中卢波波动问题的有限元法 对平面问题，应变向量8 变为 。= k 勺】r = 窿考考+ 刳7 将( 2 2 2 0 ) 式代入( 2 2 2 2 ) 式得 8 = 口口。) 其中 曰：上 2 4 b s 0 b j 0 c ，0 c lb j c i 对平面应力问题，弹性矩阵d 为 0 钆0 c j 0 c k b j c ib k 。= f v o 旦 f 2 2 2 2 ) ( 2 2 2 3 ) ( 2 2 2 4 ) ( 2 2 2 5 ) 其中e 为介质杨氏模量，d 为泊松比。对平面应燹问题，弹性矩阵d 只 需将( 2 2 2 5 ) 式中e 换成_ 冬，u 换成二l 。 l d 。1 一d 设单元厚度为t ，单元刚度矩阵和质量矩阵按( 2 2 1 3 ) 计算，可得 k 和) - 肛7 d b d v = b 7 d b tp = t a b 7 d b m 扣) - 肛7 n d v ( 2 2 2 6 ) o班。班。啦 ”o班。班o o班。啦。班 班。班。班o o啦。班。班 啦。班。班。 嘞 中国石油大学( 华东) 硕士论文第2 章弹性固体中声波波动问题的有限元法 ( 2 ) 双线性四边形单元嘲 如图2 2 所示，双线性四边形单元有4 个结点，顺序按逆时针排列， 所示坐标为整体坐标，每个结点有两个自由度。用图2 3 所示的自然坐 标f 和叮将单元映射成矩形，其坐标为局部坐标。在局部坐标系下插值 函数为 图2 2 双线性四边形单元 l = 去( 1 一烈l 一玎) 2 = 丢( 1 + 绷刊 3 = 三( 1 均( 1 4 = 三( 1 一手) ( 1 + ，7 ) 应变矩阵b 由下式确定 拈南慨口：曰，b a 式( 2 2 2 8 ) 中每个置( i = 1 , 2 ，3 ，4 ) 可表示为 图2 3 自然坐标系下的叔 线性四边形单元 中国石油大学( 华东) 硕士论文第2 章弹性同体中声波波动问题的有限元法 置= 口盟一6 盟。 西叻 oc 盟一b o n , a 玎a f 。盟一6 盟a 盟一6 盟 a qa88q ( 2 2 2 9 ) 爹毅a , b , c , d 给定如f 口2 言队( 善一1 ) + y 2 - 1 - 掌) + y 3 ( 1 + d + y 4 ( 1 一剜 6 = 丢叭玎柚奶( 1 一，7 ) + 乃( 1 + 卅y 。( 一1 刊】 。2 茅而( 叩一1 ) + x 2 ( 1 一，7 ) + 乃( 1 + 刁) + y 4 - 1 - 叩) 】 d 2 i _ ( f 一1 ) + x 2 - 1 - 善) + x 3 ( 1 + 0 4 - x 4 ( 1 一善) 】 ( 2 2 3 0 ) 坐标变换的雅克比( j a b i ) 行列式i ，| 由下式确定 m ；k 而屯毛 善要铫r y r 1 y0 o 善+ l 一善一 0 2 一孝一1 o + 0 3 善+ ，7一，7 1 0y 4 ( 2 2 3 1 ) 弹性矩阵d 同( 2 2 2 5 ) 式。 设单元厚度为t ，则双线性四边形单元的单元刚度矩阵和质量矩阵 可以表示为 k = f 1f l 口7 伽f i ，眵玎 ( 2 2 3 2 ) m ( f ) = f lf l p t n 7 l v l i i d 当a , 1 ( 2 2 3 3 ) 2 2 3 人工边界条件 由于计算机存储量的限制，选取的有限元求解域往往只是实际区 2 l 中国石油大学( 华东) 硕士论文第2 章弹性同体中声波波动问题的有限元法 域非常有限的一部分，计算区域和实际区域的其它部分间存在人工截断 边界。为保证波动模拟的可靠性与真实性，必须在人工截断边界上加上 适当的人工边界条件，以吸收到达边界处的应力波。迄今为止，人们已 构造出了许多种人工边界条件，综合考虑现有的各种边界条件对各种波 型的吸收效果，本文中选用l y s m e r ( 1 9 6 9 ) t ”) 粘性阻尼边界和 s a m m ( 1 9 9 8 ) 吸收边界的组合。 l y s m e r 粘性阻尼边界为在边界上设置一系列的粘性阻尼器，来衰 减掉大部分的反射。这一人工边界的作用效果与波源的频率和入射的角 度无关。对于二维问题，人工边界上单元l y s m e r 秸性阻尼边界可以表 示为 仃= a 0 3 - 。 ( 2 2 3 4 ) 西 、_ 其中 4 = 活：， 叫s ， 式( 2 2 3 5 ) e e 口和6 为无量纲常数，计算中均取为1 ；v p 和v s 分别为介质 中纵波和横波波速。边界单元的阻尼矩阵可以表示为 c l = p 1 a l v d g q 2 ，3 6 ) s a r m a 吸收边界就是在边界向计算区域内若干层单元上添加瑞利 ( r a y l e i g h ) 阻尼，形成边界阻尼区。瑞利阻尼矩阵可以表示为 g = n 射+ 户暂 ( 2 2 3 7 ) 其中口和为瑞利阻尼系数，满足 a + p 砰z 2 c o t 。 中国石油大学( 华东) 硕士论文第2 章弹性固体中卢波波动问题的有限元法 哆为结构自振动固有频率；毒为相对于一振型的临界阻尼比，对于 不同材料，取值约在0 0 5 0 3 之间。取两种固有振型联立方程，计算 可以确定口和口的值。 加入边界阻尼项的整体有限元运动方程( 2 2 1 4 ) 变为 m i i ( t ) + c 谤( f ) + k a ( t ) = p ( f ) ( 2 2 3 9 ) 其中 c = c l + c 2( 2 2 4 0 ) 2 2 4 整体矩阵存储和质量矩阵模型 同一般线性代数方程组相比，有限元方程( 2 2 3 9 ) 式的系数矩阵有着 自身显著的特点，皆具有对称性、正定性、高稀疏性。对于具有上万阶 系数矩阵的大型、超大型问题，这种零元素在系数矩阵中所占比例将很 大，常常超过9 0 。如不能对其有效存储，将会造成计算机存