（地质资源与工程专业论文）地震波有限元集中质量矩阵及并行算法模拟研究.pdf

中文摘要 地震波数值模拟是地震勘探中的一个重要课题，模拟结果的好坏直接影响着地震 资料处理和解释的质量。地震波数值模拟算法较多，其中有限元法是重要算法之一， 其优势在于：有限单元形状具有任意性，对于不规则界面能很好的逼近等。在有限元 地震波正演模拟中，采用s 釉a 衰减边界条件时，在介质与衰减层的交界上会产生人 为反射。衰减层内引入的阻尼矩阵使得衰减层与介质之间存在明显的物性差异，这种 差异是人为反射产生的根本原因，计算结果证实了人为反射的存在。为了消除人为反 射对地震波正演模拟的影响，在衰减层内采用变比例系数法，从而在不增加衰减层厚 度的前提下，既减缓介质与衰减层之间物性参数的突变，压制交界面上的人为反射， 又保证波在衰减层内被充分地吸收。 由变分原理得到的质量矩阵计算公式，与刚度矩阵计算方法相同，采用相同的插 值函数，在有限单元内积分得到，质量矩阵又称为一致质量矩阵。一致质量矩阵能真 实的反映有限单元内部质量分布，但一致质量矩阵是带状稀疏矩阵，求解有限元方程 时，数据量和运算量巨大，串行算法难以满足需求。当采用四边形有限单元时，在同 一方向上节点数相同的条件下，将一致质量矩阵和刚度矩阵看作是块状三对角矩阵， 采用并行求解块状三对角线性方程组的方法，实现了一致质量矩阵有限元方程的并行 求解。通过对地震波有限元模拟结果分析，证明并行算法的正确性和可行性。 集中质量矩阵是对一致质量矩阵的近似，忽略了相连单元内各个节点之间的加速 度差异，这些节点被看作具有相同的加速度。集中质量矩阵为对角阵，则可以对稀疏 的刚度矩阵采用压缩存储的方式，在求解有限元常微分方程过程中，无需将压缩矩阵 解压，可以直接做运算操作，如此，大大降低了数据量和运算量。 三维地震波数值模拟能更真实的描述实际地震波的传播，提供三维空间的波场特 征。采用集中质量矩阵，实现了三维波动方程的有限元数值模拟，较好的模拟了不同 地质条件下的波场特征。在求解过程中，对刚度矩阵采用压缩存储的方式。采用并行 算法实现了多节点区域模拟。 关键字：地震波动方程有限元数值解并行算法集中质量矩阵衰减边界条件 a b s t r a c t s e i s m i cn u m e r i c a ls i m u l a t i o ni so i o ft h em o s ts i g n i f i c a n ts u b j e c t si nt h es e i s m i c e x p l o r a t i o n , a n dt h er e s u l t sa f f e c tt h eq u a l i t i e so f t h es e i s m i cp r o c e s s i n ga n di n t e r p r e t a t i o m t h e r ea r es o m em e t h o d s a n df i n i t e - e l e m e n ti st h em o r ei m p o r t a n to n e f i n i t e - e l e m e n th a s s o m ea d v a n t a g e s ，t h ea r b i t r a r yf o r m so f t h ef i n i t ee l e m e n t , n u m e r i c a la p p r o x i m a t i o no nt h e i r r e g u l a ri n t e r f a c e sa n ds oo n ，h e ns e i s m i cs i m u l a t i o ni sc a r r i e do u tv i af i n i t e - e l e m e n t m e t h o d ，u s i n gs a r m ad a m p i n gb o u n d a r yc o n d i t i o nw i l lc a l k t h es p u r i o u sr e f l e c t i o n so n t h ei n t e r f a c e sb c 帆懒t h em e d i u m r e g i o na n dt h ed a m p i n gr e g i o n 1 1 l es p u r i o u sr e f l e c t i o n s w e r ei n d u c e db yt h ea b r u p tc h a n g eb e t w e e nt h et w or e g i o n s i no r d e rt oa v o i d i n gt h e s p u r i o u sr e f l e c t i o n s ，ap r o p o r t i o nc o e f f i c i e n tw a si n t r o d u c e da n di t sv a l u ei si n c r e a s i n g f r o mz e r ot oac o n s t a n ti nt h et r a n s i t i o nr e g i o nw h i c hi sl o c a t e di nt h ed a m p i n gr e g i o n b y t h i sm e a nw ew e a k e nt h ea b r u p t n e s so nt h ei n t e r f a c e s , m e a n w h i l et h ew a v ei nd a m p i n g r e g i o ni sa b s o r b e dp e r f e c t l yw i t h o u te x p a n d i n gt h ed a m p i n gr e g i o n t h ec o m p u t a t i o n a lf o r m u l ao fc o n s i s t e n tm a s sm a t r i xd e r i v 甜f r o mt h ev 捌o n a l p r i n c i p l e si st h e 毹哪a st h es t i f f n e s sm a t r i x s - t h e r e f o r et h ec o n s i s t e n tm 螂m a t r i x 眦 d e s c r i b et h em a s sd i s t r i b u t i n go ft h ef i n i t ee l e m e n t b u tt h ec o n s i s t e n tm a s sm a t r i xi sa s p a r s em a t r i x , s ot h e r ei st o om u c hd a t aa n dv o l u m ec a l c u l a t i o n w h e nt h em a s sm a t r i xa n d t h e5 倘e n e s sm a t r i xa mr e g a r d 鹪t h et h r e eb l o c k - t r i d i a g o n a lm a t r i x e s , f i n i t e - e l e m e n t e q u a t i o np a nb es o l v e db yt h ep a r a l l e l 姆f i t h mf o rb l o c k - t r i d i a g o n a ll i n e a rs y s t e m s 1 1 1 e n u m e r i c a ls i m u l a t i o nr e s u l t sp r o v et h a tt h ep a r a l l e la l g o r i t h mo nm u l t i - c p ui sc o r r e c ta n d f e a s i b l e h e nt h ef i n i t ee l e m e n ti ss os m a l lt h a tt h en o d e so f t h et a m ef i n i t ee l e m e n to a n $ e e l n a st oh a v et h es a m ea c c e l e r a t i o n s ，t h ec e n t r e dm a s $ m a t r i xw a su s e dt os u b s t i t u t ef o rt h e c o n s i s t e n tm a s sm a t r i x , a n dt h ep r e c i s i o ni s n tb e l o wt h ec o n s i s t e n tm a s sm a t r i x s 1 1 抢 c e n t r e x ln 埔船m a t r i xi sa d i a g o n a l - m a u i x , a n dt h es t i f f n e s sm a t r i xc a nb ec o m p r e s s e d t h u s t h ed a t aq u a n t i t ya n dt h ec a l c u l a t i o na m o u n ta r es a v e d 3 ds e i s m i cn u m e r i c a ls i m u l a t i o n 咖t r u l ys i m u l a t et h ea c t u a ls e i s i m i cw a v ea n dp r o v i d e3 dc h a r a c t e r i s t i c so f t h ew a v ef i e l d f i u l t e - e l e m e n te q u a t i o nu s i n gt h ec e n t r e dm a s sm a t r i xi ss o l v e di nt h e3 d n 圮n u m e r i c a l s i m u l a t i o nb yt h ep a r a l l e la l g o r i t h mo nm u l t i c p ue x p a i l d st h es i m u l a t i o ns p a c e k e y w o r d s ：s e i s m i cn u m e r i c a ls i m u l a t i o n , f i n i t e - e l e m e n tn u m e r i c a ls o l u t i o n s ，p a r a l l e l a l g o r i t h m ，c e n t r e dm a s sm a t r i x , d a m p i n gb o u n d a r yc o n d i t i o n , 独创性声明 我呈交的学位论文是在导师指导下个人进行的研究工作及取得的研 究成果。尽我所知，除了文中特别加以标注和致谢的地方外，论文中不 包含其他人已经发表或撰写过的研究成果，也不包含为获得其它学位或 证书而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均 已在论文中作了明确的说明并表示了谢意。特此声明。 声明人( 签名) ：王目墓塑! z 年岁月3 0 日 关于论文使用授权的说明 本人完全了解中国石油大学有关保留、使用学位论文的规定，即： 学校有权保留送交学位论文的复印件，允许学位论文被查阅和借阅；学 校可以公布学位论文的全部或部分内容，可以采用影印、缩印或其他复 制手段保存学位论文。特此说明。 说明人( 签名) ：王目墓指导教师( 签名)岁a 孓i ? t 地震波有限元集中质量矩阵及并行算法模拟研究 创新点摘要 1 、在分析衰减吸收边界条件的基础上，对s a n n a 衰减边界条件做了改进。通过 在衰减层内设置过渡带，并在过渡带内采用变比例系数的方法，来减缓衰减层内阻尼 参数的突变，克服了介质与衰减层之间界面上产生的反射波干扰，压制了人为截断边 界上的边界反射( 见第3 章) 。 2 、将并行求解块状三对角线性方程组算法应用到地震波有限元方程求解中，实 现了地震波一致质量矩阵有限元数值解的并行求解。通过对算法的分析，根据求解过 程中数据的变化，将系数矩阵的三条对角线上的子矩阵分开存储，在求解过程中，仅 保存有效的数据。一致质量矩阵对应的有限元方程数据量和运算量巨大，通过并行运 算，实现了串行算法难以实现的多节点区域模拟。通过对不同地质模拟结果的分析， 证实了算法的正确性和有效性( 见第4 章、第6 章) 。 3 、当有限单元的尺寸远小于地震波波长时，一致质量矩阵可以近似为集中质量 矩阵。采用质量矩阵近似的方法和并行算法，实现了三维空间地震波方程的有限元数 值解，并模拟了不同地质模型的波场特征。( 见第5 章、第6 章) 。 第l 章弓i 言 第1 章引言 地震数值模型是在已知地下结构和地质介质性质的条件下，研究地震波的传播规 律，并模拟地表观测信号，在地震勘探中具有重要意义。地震反演则是在己知地表观 测信号，根据地震波传播的规律，通过分析地表信号，来反推地下未知构造和未知介 质的物性。因此，能否取得一个好的正演模拟结果，直接影响着地震勘探方法的发展、 地震资料处理的效果和资料解释的质量。随着石油勘探技术的不断发展，地震勘探开 始由构造勘探转向岩性勘探和隐蔽性油气藏勘探，对地震波正演模拟技术也提出了更 高的要求，高精度的数值模拟技术和三维数值模拟技术研究有着重要意义。 1 1 论文选题的目的及其意义 地震波模拟是研究地震波传播问题的重要方法。自上世纪6 0 年代以来，地震模 拟技术得到了很好的发展。目前，比较常用的方法主要有：有限元法【n 、射线法、有 限差分法、反射率法、傅立叶方法、虚谱法和边界元等。射线法是建立在高频近似的 假设之上，适用于光滑的非均匀结构和弱各向异性问题。反射率法适用于横向均匀和 层状介质，对于复杂的介质模拟有一定的局限性。傅立叶方法具有较高的精度，且产 生的频散较小，但仅适用于计算光滑变化的模型。虚谱法的基本思想是通过付氏变换， 在时间域上作差分，避免对空间坐标的求导运算。虚谱法在频率域分辨率高，但在时 间域分辨率较低。有限差分是通过有限差分算子，将波动方程离散化以后，以差分来 代替微分。有限差分只考虑有限个离散点上的函数值，而不考虑点的邻域函数的变化； 有限差分法在空间域内具有较高的分辨率，可以很好地适应剧烈变化的地下介质情 况但是在频率域内分辨率却很低，此外，有限差分法的最大不足是只适用于规则的 网格剖分，不能处理不规则形状的边界条件；对于自由边界条件，使用有限差分法和 傅立叶法来实现都比较困难。 基于变分原理的有限元法是采用剖分插值的思想，对求解空问分段近似【3 】。首先 将求解空间离散为一系列有限单元，由单元的节点波场值插值得到求解空间内所有节 点的波场值。有限单元剖分具有任意性，对于不规则的地质体和弯曲界面，可以通过 不同形状的有限单元来逼近。因此，比较适合几何和物理条件比较复杂的问题。在推 导有限元方程时，使边界条件同时满足。因此，有限元处理边界问题相对比较简单。 所以，有限元法具有不受边界几何条件的限制，单元网格剖分灵活，且处理边界条件 第1 章弓l 言 比较方便的独特优势。 有限元法也存在不足，其需要巨大的数据量和运算量。因此，有限元很难实现多 节点区域的模拟。在研究过程中，也会因数据量和运算量巨大而影响研究的进展。面 对这些难题，人们大多将一致质量矩阵近似为集中质量矩阵，变稀疏的对称质量矩阵 为对角阵，来降低求解运算的难度。一致质量矩阵是非对角阵，假设空间有一个节点， 整体一致质量矩阵就是为n n 阶对称矩阵。采用一致质量矩阵，求解空间内的节点 之间相互耦联，当一个节点存在惯性力时，其余节点也存在惯性力；而采用集中质量 矩阵，求解空间内的节点之间则是解耦的【4 1 。一致质量矩阵与刚度矩阵采用相同的位 移插值函数，反映了单元质量的实际分布，因此，一致质量矩阵能提供精确的结果， 但运算量巨大。当前，计算机技术有了突破性的发展，在新的运算环境下，采用一致 质量矩阵有限元来实现地震波数值模拟是完全有可能的。 随着石油勘探技术的发展，为了适应实际生产研究的需要，三维地震波数值模拟 已经成为一个新的研究热点。在实际地震勘探中，地震波是在一个半无限空间中传播， 在三维空间内模拟地震波的传播问题，相对二维空间更真实一些。在三维空间内，传 播的是一个球面波，波的能量是按照球面l r 2 来衰减的；而在二维空间内，将三维空 间中的球面波描述为二维空间的圆柱，波的能量是按1 ，来衰减的 5 1 。三维空间的数 值解提供一个立体的波场，不局限在单一平面上；而二维空间的数值解是一个平面上 的波场，不能全面地描述来自空间上的波场。因此，三维地震波数值模拟要优于二维 的数值模拟。目前，在三维空间中。比较常用的地震波数值解方法有；有限差分、虚 谱法、边界元和有限元等方法，其中有限元法数值模拟研究并不多。采用有限元法求 解三维空间的地震波动方程，有望体现有限元的独特优势，为地震勘探提供更丰富的 波场信息，更好地推动地震勘探技术的发展。 以往的有限元数值解主要是串行算法，难以解决多节点区域正演模拟中的巨大数 据量和运算量问题。随着计算机硬软件计算的发展，微机群和基于消息传递的h 伊i 并行算法开发平台的出现，在一定程度上解决了巨大数据量和运算量的难题。有限元 的基本思想就是：先化整为零，再积零为整，并行处理的基本原则是：化整为零，分 而治之，有限元的基本思想与并行处理的原则一致嘲。因此，有限元法适合并行处理。 第l 章引言 1 2 地震有限元数值解发展现状 有限元法【l 】为一种有效的数值计算方法，在工程求解运算中占有重要地位，应用 领域广泛。有限元法最初在上世纪五十年代主要用于飞机结构分析，到上世纪六十年 代开始推广到弹性力学平面问题，后来开始应用到热传导、电磁场、流体力学等连续 性问题中。地球物理数据模拟是有限元法和边界元法的重要使用领域1 7 一。 在上世纪七十年代，哥伦比亚大学将有限元法开始应用到地震波勘探领域。此后， 有限元法以其独特的优势，吸引了众多科研工作者畔圳的注意。m a r f u r t 1 5 】( 1 9 8 4 ) 研 究对比了传统的有限差分法和有限元法的精度，研究表明，在各向均匀介质中，对于 标量波动方程和泊松比小于o 3 的弹性波动方程，有限元和有限差分效果相当；对于 弹性波动方程，当泊松比大于0 3 时，有限元法的计算精度要高于有限差分法 p a d o v a n i 【1 6 】( 1 9 9 4 ) 研究了在地震波有限元模拟中的低阶有限元和高阶有限元法的精 度。s e r o n ( 1 9 9 6 ) 等研究了弹性介质中的有限元法数值解。8 a r m a t l 7 1 研究了弹性波传 播有限元模拟中衰减吸收边界问题，给出了应用于有限元的无反射吸收边界条件。有 限元法在射线追踪、逆时偏移和多分量模拟等方面，也得到了一定的研究和应用 0 8 - 2 2 1 。 在国内，上世纪八十年代是有限元法在地震勘探领域发展的一个重要时期。崔力 科 i m 2 和牟永光( 1 9 8 7 ) 等田。瑚人对有限元数值解的稳定性、精确度等问题做了详 细的研究，求解了无阻尼地震模型、三维阻尼地震模型的有限元数值解以及有限元偏 移成像技术。杜世通b l o l ( 1 9 8 2 ) 从变分原理出发，推导了弹性波动方程的有限元方 程式，并分别给出了二维变速不均匀介质、非理想弹性介质中波动方程有限元数值解 法和有限元偏移方法，并将该技术应用到井间观测系统中( 1 9 8 6 ) 。但在上世纪八十 年代，由于计算机的处理运算性能有限，制约了有限元的进一步发展和应用。 到了上世纪九十年代，计算机技术有了一定的改善，人们继续着对有限元法的深 入研究 2 6 - 3 6 。王尚旭t 3 2 1 ( 1 9 9 0 ) 利用有限元方法研究了双相各向同性介质中弹性波 传播规律。邵秀民等【3 3 i ( 1 9 9 3 ) 研制了地震波传播数值模拟有限元程序系统，将不规 则四边形有限单元应用到了地震波正演模拟中。王润秋瞰】( 1 9 9 5 ) 等将求解空间有二 维扩展到了三维空间，通过三维空间一维化的思路，将三维空间映射到一维空间中， 并给出了对应的有限元解三维弹性波方程的并行算法。周辉和徐世浙【3 6 】( 1 9 9 7 ) 对各 第1 章引言 向异性介质中波动方程的有限元数值解的稳定性做了深入探讨。此时，巨型机和并行 机价格昂贵，无法普及到各个研究和生产领域。 到了本世纪，高性能的计算机开始出现，计算环境有新的改善 3 3 - - 4 2 1 。周删”1 等 人利用六面体单元和三线性插值，实现了各向异性介质中三维三分量波动方程的有限 元数值解。柯本喜1 3 8 】( 2 0 0 1 ) 联合使用三角形和四边形有限单元，求解二维起伏地表 条件下标量波动方程响应，并给出了较理想的人工边界条件，达到了较好的效果。张 剑锋1 3 9 】采用集中质量矩阵实现了各向异性介质的弹性波有限元数值解。张美根和王妙 月【柏，a 2 ( 2 0 0 0 ) 实现各向异性介质的有限元正演模拟、有限元偏移技术和有限元参 数反演技术。杨顶辉【4 3 】( 2 0 0 2 ) 对孔隙各向异性介质中基于b i s q 模型的弹性波传播 进行了有限元方法模拟。马德堂和朱光明m l ( 2 0 0 4 ) 把有限元法和伪谱法相结合，求 解起伏地表条件下的解弹性波动方程。黄自萍 4 7 1 等( 2 0 0 4 ) 则采用区域分解的方法， 将有限元和有限差分相结合求解弹性波动方程，在起伏地表采用有限元法来模拟，在 规则区域使用有限差分法。在地震勘探中，有限元法还广泛地应用到油藏描述和应力 分析中。张帆【船】( 2 0 0 0 ) 采用有限元的方法构造应力场，来预测地下构造裂隙的发育 情况。在研究过程中，人们主要还是采用集中质量矩阵来实现地震波数值解。 有限元法在工程上的贡献巨大，自有限元法诞生以来，有限元法仍然在不断地发 展，先后出现了广义协调元、理性有限元和无单元等新型有限元。贾晓峰和胡天跃【4 1 】 将移动最小二乘和无单元g a l e r k i n 法开始应用到地震波正演模拟和成像中，并取得了 很好的效果。无单元g a l e r k i n 法无需剖分单元网格，但在积分求解中需要划分背景网 格辅助。该方法是采用拟合的思想，节点处的位移不满足6 函数性质，对畸变有一定 的抵抗性，但在边界处理上存在一定的难度，且数据量和运算量巨大。 当前有限元结构分析程序库发展很快，已经出现了如s a p 系列、a d i n a 、 a n y s y s 等比较成熟的大型程序库，但是这些软件都是基于传统的串行计算机体系 结构来设计，并在串行计算机上运行。由于受到计算机硬件水平的限制，在应用规模 和运算速度上也存在一定的局限性，约束了大型和超大型有限元的应用。自上世纪八 十年代以后，有限元的并行算法在力学分析上发展较快，开始出现各种有限元并行分 析的方法【5 2 1 。采用多c p u 的并行处理算法，可以突破单一c p u 的物理限制5 2 】，提高 求解问题的效率。目前，地震波有限元数值解还是传统的串行算法，每一台机器的物 4 第1 章引言 理空间和处理能力是一定的，单靠单一c p u 来是实现有限元数值解，求解空间有限。 随着计算机硬软件技术的进一步发展，高性价比的微机群开始进入各研究所和生产部 门在地震勘探领域，对于海量数据的求解运算不再是难题。地震资料处理技术开始 由串行求解转向多c p u 的并行化处理，数值运算环境有了重大变革，有限元并行处 理技术可以在新的环境下，获得进一步的发展。 1 3 研究思路和主要内容 地球介质是一个半无限空间，采用数值模拟只能实现一个有限求解区域【辨5 9 1 。 人为边界产生的响应( 如反射等) 是实际介质中不存在的，也是模拟中所不希望的干 扰。消除或者减弱人为边界反射，是地震波正演模拟中的主要内容。在有限元正演模 拟中，比较常用的边界条件是衰减边界条件，其具有稳定性好，吸收角度全的优点。 衰减边界条件是根据粘滞性介质的阻尼形成机理来构建，阻尼效果的好坏与阻尼因子 有关。通常，整体阻尼矩阵由整体质量矩阵和整体刚度矩阵通过比例系数来构建。因 此，在衰减层内直接添加衰减边界条件，则使得介质与衰减层的物性参数存在突变， 从而产生干扰反射波。通过采用变比例系数法，来缓解衰减边界区域内的物性突变， 避免干扰反射，同时也保证了边界条件的吸收效果。 一致质量矩阵的计算方法与刚度矩阵相同，采用相同位移插值函数，在有限单元 空间内积分计算。采用一致质量矩阵能准确地描述地震波传播特征，因此，一致质量 矩阵较集中质量矩阵更准确一些。但是一致质量矩阵与刚度矩阵一样，也是稀疏矩阵， 当空间中有玎个节点时，一致质量矩阵是n x n 阶矩阵，所以有限元方程求解过程中 数据量和运算量巨大，限制了其应用和发展。有限元法的基本思想是通过将求解空间 剖分成一系列有限单元，然后再集合有限单元参数矩阵为整体参数矩阵，得到整体有 限元方程 6 2 1 。并行处理技术基本原则是将一个整体任务划分成几个子任务，并分别处 理解决嘲因此，有限元法的基本思想和并行处理技术的基本原则相一致，所以有限 元法可以采用并行处理的技术来解决。有限元并行处理主要包括两个部分：参数矩阵 的并行计算和有限元常微分方程的并行求解。其中有限元常微分方程是通过递推求解 线性方程组来实现，并行求解系数矩阵为稀疏矩阵的线性方程组，算法实现难度较大， 可以将系数矩阵看作是块状三对角矩阵。通过并行求解块状三对角线性方程组的算 法，来实现一致质量矩阵对应的有限元方程并行求解。 第1 章引言 三维空间的数值解描述的是球面波场，能得到立体空间上的波场信息和传播规 律，所以三维空间的数值模拟能接近实际地震波传播的性质。采用有限元法实现三维 空间波动方程正演模拟，可以发挥有限元法的优势，提供丰富的地震波场信息。采用 一致质量矩阵对应的有限元，数据量和运算量巨大。当有限单元的尺寸与地震子波波 长之比足够小时，可以将质量矩阵近似为集中质量矩阵，变稀疏对称矩阵为对角阵， 此时，刚度矩阵可以采用压缩存储的方式。如此，大大的降低运算量和较少运算量， 在此基础上，采用并行求解的方法，扩展了求解空间。 课题的主要内容及成果有： l 、研究有限元正演模拟中的衰减边界条件。在分析衰减吸收边界条件机理的基 础上，证实了交界面上容易产生人为反射波，这种人为反射是由于衰减层与介质之间 介质的物性突变所引起的。通过在衰减层内设置过渡带，对衰减边界条件做了改进。 在过渡带内对阻尼因子采用变比例系数的方法，来减缓衰减层内阻尼参数的突变，克 服了人为反射波的产生，消除了截断边界反射。 2 、采用并行算法，实现了二维空间地震波一致质量矩阵有限元数值模拟。一致 质量矩阵能更准确的描述有限单元的质量分布，当采用四边形有限单元，且同一方向 的节点个数相同时，一致质量矩阵和刚度矩阵可以看作块状三对角矩阵，将并行求解 块状三对角线性方程组的算法应用到地震波有限元方程求解中，可以实现地震波有限 元并行求解。根据对求解过程中数据量、运算量和通信量的分析，在均衡数据量和运 算量的基础上，可增加首末c p u 上的求解空间，减少c p u 个数与模型空间大小之间 的限制；在算法实现过程中，将系数矩阵在三条对角线上的矩阵分开存储，并仅保留 有意义数据，对无意义数据不记录，如此节约了数据空间。通过对地质模拟结果韵分 析，证实了算法的可行性和正确性。 3 、采用并行求解的方法和质量矩阵近似的方法，实现三维空间的地震波有限元 数值模拟。在求解过程中，对刚度矩阵采用压缩存储的方法，仅记录相关节点的信息， 减少了数据量；在数值运算时，直接对压缩矩阵做运算操作，不必解压到原先的数据 形式，减少了运算量。 6 第2 章波动方程有限元数值解法 第2 章波动方程有限元数值解法 2 1 有限元数值方法概述 有限元法是一种重要的数值算法。将求解空间划分为一系列有限单元，对每个单 元采用插值函数来近似场函数在单元内的局部特征。在理论上讲，当插值函数越接近 场函数，数值解越接近精确解，通常场函数比较复杂，一般采用有限项多项式来代替 场函数。当增加有限单元的节点数，或减小有限单元的尺寸时，数值解将逐渐逼近精 确解。 l 、有限元法的主要步骤【3 l 采用有限元法来求解物理问题，需要遵循以下步骤： 求解区域离散化。将求解区域划分成一系列的有限单元，单元的类型要根据求 解问题确定。选择适当的插值函数。根据求解精度要求和单元类型，确定插值函数； 计算单元质量矩阵和刚度矩阵。根据平衡条件，推导有限元p 的单元质量矩阵肘函 和单元刚度矩阵墨赫，计算单元质量矩阵和刚度矩阵，得到单元有限元方程。集合 所有单元质量矩阵和刚度矩阵，得到整体质量矩阵和刚度矩阵，来获得整体有限元方 程。求解节点的未知位移。联合边界条件和载荷向量，根据求解问题要求，选定求 解有限元方程的方法，求解有限元方程。 2 、微分方程近似解法 偏微分方程的求解要求满足一定的初始条件和边界条件嘲，称为初边界条件解。 边界条件的类型有两个：一个是基本边界条件( 又称狄里赫利边界条件) ；另一个是 自然边界条件( 又称纽曼边界条件) 。在研究应力应变问题时，这两类边界条件相对 应于位移边界条件和应力边界条件。变分法将导出自然满足狄里赫利边界条件，即未 知函数“在边界上具有给定值；g a l e k i n 剩余量法将要求自然满足纽曼边界条件，即 未知函数的导数在边界上具有给定值。首先根据问题的边界条件性质，选择不同的 近似方法。此外，r a y l e i g h r i t z 法要求所求物理问题的变分泛函必需存在嘲，且有明 确的泛函表达式。有些物理问题的泛函可能不存在，或者是难以得到难明确的泛函关 系，对于这类问题，采用c r a l e k i n 剩余量法比较方便。g a l e k i n 剩余量法是一种误差最 小化的近似方法，直接对控制方程求解，而不要求问题的变分泛函一定存在，其加权 7 第2 章波动方程有限元数值解法 函数取插值函数的基函数。 3 、有限单元类型及插值函数 在二维空间中，有限单元的形状主要有三角形有限单元和四边形有限单元。三角 形有限单元大多是采用线性插值，计算简单方便，可以逼近任意的弯曲界面，这是有 限单元法优于其他求解方法的优势。虽然三角形网格有较好的适应性，但当内角很小 时，计算的精度和稳定性很差州。四边形有限单元采用的是曲线插值函数，采用不规 则四边形有限单元也可以近似不规则地质体，四边形有限单元的稳定性好，即使采用 较扁的四边形有限单元，也能得到较好的计算效果哗】。 4 、有限元收敛的要求 在有限元法中，场函数的总体泛函是由单元泛函集合而成。如果采用完全多项式 作为单元的插值函数，则有限元解在一个有限尺寸的单元内可以等于精确解。但是实 际上有限元的插值函数只能取有限项多项式，所以只能是接近精确解的近似解。有限 元趋于精确解，需要满足的收敛准则网： 完备性要求 如果在泛函中场函数的最高阶导数是坍阶，则有限元解收敛的条件之一是单元 场函数的插值函数至少是所次完备多项式。或者说插值函数中必须包括本身和直至埘 阶导数为常数的项。当插值函数是多项式时，满足单元内部函数的连续性要求。主要 是单元交界面上的连续性。 协调性要求 如果泛函中的最高阶导数是m 阶，则插值函数在单元界面上必须具有c 叫连续 性，即在相邻单元的交界面上函数应有直至脚- l 阶的连续性导数。选取的单元具备完 备性和协调性时，有限元才是收敛的，即当单元尺寸趋于零时，有限元解趋于精确解。 2 2 波动方程有限元数值解方法原理 传统的位移元主要是在单元内假定位移场，通过变分原理建立有限元方程。变分 问题就是求解泛函的极值，即能量的极值。根据哈密尔顿原理，一个保守力场系统中， 系统的动能与位能之差的时间平均值为最小嘲。有限元法是把求解空间离散为一系列 有限单元，在每一个单元内，根据一定条件用简单的函数来近似化未知函数的局部特 征。只要有限单元的大小尺度合适，选择一定的插值函数，近似函数的精度就满足要 3 第2 章波动方程有限元数值解法 求。只要在单元内位移满足于完备性、协调性和稳定性要求，有限元方程就是收敛的 【蚓。以三维的有限元推导为例，说明有限元方程的推导。 在三维空间上，直角坐标系下的标量波动方程为 窘- i - 雾4 - 窘= 吉鲁o t c 2 一t ， 缸。 却2昆2 v 2 2 却为波场函数l f “o 膨f ) 。在二维空间中，标量方程式没有( 2 - - 1 ) 式中y 项。采用 r a y l e i g h - r i t z 法，取三维空间波动方程的变分方程： 胍p 2 窘+ 雾+ 害一窘 = 0( 2 2 ) 使用分部积分，化简上式： v 2m 罡警+ 万o u 万o & + 西3 u i a # u + 窘面 蚴 = j l i 罡& 卜+ 压陪面卜+ j l ，罡擞k c 2 吲 考虑到边界岛、岛和岛上的给定值，在边界变分为零。得到三维问题的变分积分， 为： ，= 埘惝饼+ 盼2 到 对于二维问题，则为： ( 2 - - 4 a ) ，= 三l 2 障) 2 + 2 + 2 窘“ ) 掀( 2 - - 4 b ) 对有限单元e 推导有限元方程。用单元节点波场值表示单元内任意点的波场值。取试 验函数： ”。= 眵版 ( 2 5 ) 二维空间上，主要的有限单元类型有三角形有限单元和四边形有限单元。三节点 的三角形有限单元插值函数移) 为线性函数，可以表示为： 9 第2 章波动方程有限元数值解法 彩 = 1 x il 毛 z 船心x 2 z ：2 ， 【ix 4 毛 而毛 x 2 2 2 屯毛 x 4 2 4 ( 2 6 a ) ( 2 6 b ) 其中，毛) 为有限单元的第i 个节点坐标。在三维空间中，以八节点的六面体有 限单元为例，选插值函数为彩) 可以表示为： 铆= ( 1z ) ，zx yy z 荔x y z 其中k ，以，毛) 为有限单元的第f 个节点坐标。可以得到有限单元方程： y i 毛 x 2 j ，2 9 2 黾y s z j x 4 y 4 2 4 x 5 儿毛 x b y 6 2 6 x 7 y 7 2 7 h 地2 l ( 2 - - 6 e ) j 【谚办施+ u m n y 2 胍( 警誓+ 等等+ 警警卜= 。c 2 一， 其中v 为速度。通常上式表示为： 肘豇赫+ 眉赫“赫= 0 ( 2 - - 8 ) 其中 = j 【谚力施 ( 2 9 a ) 碥f ( 芸誓+ 等等+ 警警 拉( 2 - - 9 b ) 1 0 为示表以1，j可 毛龟乃 而屯b 粉 i目屯船 x 插 a 例 = 旬 辨 融单限有形边四的点节 个四以 碱矾矾撕张碱辅矾砌勿乃“办办新靠w妇乃儿儿趵虮m耽均儿儿儿乃儿 “砌巧鼽职靠新鲰 而如乃“办办新力朋儿儿凡儿儿乃以 而砌彤所斯靠新鲰 第2 章波动方程有限元数值解法 称朋赫为单元质量矩阵，置知为单元刚度矩阵，屹，为单元位移向量，西知为“知对 时间t 的二阶导数。 为了计算研究的方便，在二维空间采用长方形有限单元，单元节点编号如图2 一 l a 所示。将有限单元节点坐标代入( 2 6 a ) ，得到插值函数： b z 识= ( 1 - 故一言) 丸= ( t 一言) 詈 。2 州， 九= 舡云) b y c a ) 二维空问有限单元( b ) 三维空间有限单元 图2 1 空间中的有限单元 将插值函数带入单元的质量矩阵矩阵和单元刚度矩阵二维表达式，得到： m k ：氅 j o 4 2 24 12 21 12 21 4 2 24 ( 2 一l l a ) 第2 章波动方程有限元数值解法 x ：妥 j b 口 l l 2 l 2 1 l 2 l 一1 l 2 1 2 一l l 1 2 一l 1 一 2 1 2 l 口 _ i b 1 1 1 2 1 2 一l l l 2 1 2 l 2 1 2 l 一1 l 2 1 2 一l l ( 2 一l l b ) 在三维空间，采用长方体有限单元( 图2 一l b ) ，通过手算得到单元参数矩阵。 设有限单元在x 方向的长度为口，y 方向的长度为b ，z 方向的长度为c ，长方体单元 节点的坐标为：( o ，0 ，o ) 、( 0 ，b ，o ) 、( 口，b ，o ) 、( 口，0 ，o ) 、( 0 ，0 ，c ) 、( o ，b ， c ) 、( 口，b ，c ) 和( 口，0 ，c ) ，将节点坐标代入( 2 - - 6 b ) ，得到插值函数： 魂= ( 一言 ( ，一詈 ( t 一詈) 如= ( ，一激一句 九= 言詈( 1 一詈) 九2 争a k 歌b a 一习c ) 。2 吨， 九= ( 一言) ( ，一詈) 詈 。1 。 九= ( 一珧詈 2 一 c 将插值函数( 2 一1 2 ) 代入( 2 9 ) ，得到单元一致质量矩阵( 2 - - 1 3 a ) 式和单元刚度 矩阵( 2 - - 1 3 b ) 式： 1 2 嚣爿口* 办 = 第2 章波动方程有限元数值解法 = 寨 确= 轰 们 + 一 3 6 b 4 4 2 2 2 2 一l l a b j 一 3 6 c 84 2 4 4 2l2 484 2 2 4 2l 2 484l24 2 4 2 482l2 4 4 2l28424 2 4 2l48 4 2 l2422 484 2l2 4 4 2 48 422 4 2 l一1 2 44 2l22 2 442 一l 一2 2 4 224 2 11 2ll一242 2 l221244 一l一22l一2 44 2一ll24 22 4 4 2 2 2 2 l l 一2 2 4 4 一l l 2 2 42l 2 4 2 124 212 - 4 之一l _ 2 _ 4 - 2 一l - 2 - 4 - 2 1 - 2 222一l 一222l 一4一ll2 4ll一2 14 4 2 一l一442 2224 222 4 2 4 _ 2 一l _ 2 1 - 2 _ 4 - 2 1 2 一l _ 2 _ 4 - 2 4 _ 2 一l - 2 - 4 - 242l2 - 1242l - 2l242 - 42 l24 ( 2 1 3 a ) ( 2 一1 3 b ) 乞o 。2 o 也2 4 。o也：2也o 4 第2 章波动方程有限元数值解法 集合空间内所有有限单元方程，得到整体有限元方程： j l f 疗+ k u ：0( 2 1 4 ) 称m 为整体质量矩阵，置为整体刚度矩阵，u 为整体位移向量， 为u 关于时间的 二阶导数。因此，有限元方程( 2 1 4 ) 是关于时间t 的二阶常微分方程。 在质量矩阵有两种计算方法1 4 , 6 5 - - 6 9 3 ，一致质量矩阵( 协调质量矩阵) 和集中质量 矩阵。根据变分原理推导得到的质量矩阵表达式( 2 - - 9 a ) ，称为一致质量矩阵。一致 质量矩阵对应的有限元方程是空间勰联的，即当内部任一节点承受一个力时，内部其 余节点也会瞬时产生一个加速度，节点间距离远近的不同，加速度大小不i 司t 4 1 。单元 一致质量矩阵是对称矩阵，集合得到的整体一致质量矩阵也是对称矩阵，若空间上的 节点数为, ，则一致质量矩阵膨与剐度矩阵量一样，都是n x n 阶对称矩阵。集中质 量矩阵是一致质量矩阵的一种近似形式，只有当有限单元的尺寸相对于地震波波长度 之比足够小，单元节点之间位移对时间的二阶导数之间的差异可以忽略的时候，这种 近似才可行。当单元质量分布均匀时，单元的各节点上分担的质量可以按照线段、面 体或者体积来确定【6 5 阐。两种质量矩阵类型不同，求解对应有限元方程在数据量和运 算量上的差异也较大。 2 3 有限元方程差分求解 偏微分方程求解初值问题，要求给定t = o 时刻全域上的未知函数值，由此推定f = o 时的求解未知函数。要求在时间域上的求积分。在地震波有限元数值解中，通常 在空间上采用有限元，得到关于时间的二次常微分方程式( 2 一1 4 ) ，采用有限差分法 来求解。因此，在数值求解运算中，也需要保证有限差分法求解有限元方程的算法稳 定性。 有限差分法主要分为无条件稳定和条件稳定两种 6 6 1 。在任意初始条件下，对于任 意的网格间距血和时间步长4 f ，解都不会无限的增大或者振荡，则这种差分格式为 无条件稳定；当对网格间距彳x 和时间步长加的取值在一定范围内时，才不会出现解 的无限增大或者振荡，则这种差分格式属于条件稳定。对于条件稳定差分格式，确定 稳定条件的方法把当系数矩阵的谱半径要小于1 ，即庐1 以此来确定空间网格间距 出和时间步长小的取值范围。在满足条件限制之下，算法才能稳定，避免发散，称 1 4 第2 章波动方程有限元数值解法 这种条件限制为稳定条件。 有限差分按照差分格式可以分为隐式和显式两类。隐式法是根据当前时刻以及前 几时刻的波场值，来建立下一时刻波场值的方程组，根据方程组来确定下一时刻波场 值；而显式法是根据当前时刻以及前几时刻的波场值，直接递推下一时刻的波场值。 比较常用的无条件稳定隐式差分方法有n e w m a r k 方法。大多数隐式有限差分都属于 无条件稳定，但需要构建耦联的方程组，对于节点数多的大型有限元问题，采用隐格 式求解有限元方程运算量增大。显式有限差分格式直接递推，并不需构建耦联方程， 但是显式有限差分格式大多是条件稳定，若想稳定求解有限元方程，需要对空间网格 大小和时间采样间隔做限制。 下面分别选取隐式有限差分格式和显式差分格式中比较有代表性的求解方法： n e w m a r k 差分和中心差分为例，推导地震波有限元方程求解的表达式： l 、n e w m a r k 法 n e w m a r k 法是比较有代表性的隐式差分格式，是