波动方程数值模拟的三种方法及对比
地震波波动方程数值模拟方法(严选优质)
地震波波动方程数值模拟方法地震波波动方程数值模拟方法主要包括克希霍夫积分法、傅里叶变换法、有限元法和有限差分法等。
克希霍夫积分法引入射线追踪过程，本质上是波动方程积分解的一个数值计算，在某种程度上相当于绕射叠加。
该方法计算速度较快，但由于射线追踪中存在着诸如焦散、多重路径等问题，故其一般只能适合于较简单的模型，难以模拟复杂地层的波场信息。
傅里叶变换法是利用空间的全部信息对波场函数进行三角函数插值，能更加精确地模拟地震波的传播规律，同时，利用快速傅里叶变换(FFT)进行计算，还可以提高运算效率，其主要优点是精度高，占用内存小，但缺点是计算速度较慢，对模型的适用性差，尤其是不适应于速度横向变化剧烈的模型.波动方程有限元法的做法是：将变分法用于单元分析，得到单元矩阵，然后将单元矩阵总体求和得到总体矩阵，最后求解总体矩阵得到波动方程的数值解；其主要优点是理论上可适宜于任意地质体形态的模型，保证复杂地层形态模拟的逼真性，达到很高的计算精度，但有限元法的主要问题是占用内存和运算量均较大，不适用于大规模模拟，因此该方法在地震波勘探中尚未得到广泛地应用。
相对于上述几种方法，有限差分法是一种更为快速有效的方法。
虽然其精度比不上有限元法，但因其具有计算速度快，占用内存较小的优点，在地震学界受到广泛的重视与应用。
声波方程的有限差分法数值模拟对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述：)()(2222222t S zu x u v t u +∂∂+∂∂=∂∂ （4-1） (,)v x z 是介质在点（x , z ）处的纵波速度，u 为描述速度位或者压力的波场，)(t s 为震源函数。
为求式（4-1）的数值解，必须将此式离散化，即用有限差分来逼近导数，用差商代替微商。
为此，先把空间模型网格化（如图4-1所示）。
设x 、z 方向的网格间隔长度为h ∆，t ∆为时间采样步长，则有：z∆,i j1,i j +2,i j+1,i j-h i x ∆= （i 为正整数）h j z ∆= (j 为正整数)t n t =∆ (n 为正整数)k j i u , 表示在(i,j)点，k 时刻的波场值。
CSC频率—空间域波动方程数值模拟
CSC频率—空间域波动方程数值模拟吕晓春;顾汉明;成景旺;周丽【摘要】针对频率空间域波动方程数值模拟需要巨大内存空间的现状,提出了利用列索引压缩存储(CSC)技术存储大型稀疏非对称复数型的矩阵系数.CSC技术将系数矩阵转化为三个一维数组来存储,分别存储系数非零元素、非零元素对应所在的行以及每列起始非零元素所在位置.经CSC技术压缩存储后显著减少了内存空间及计算量,在计算时只有少许的非零元素参加计算,且根据三个一维数组可以简便地找到对应的非零元素,进而采用LU分解快速而精确地求解.本文基于Jo等提出的最优化9点差分方法,首次应用CSC技术在频率空间域进行二维声波方程数值模拟.通过对Corner-edge模型和二维Marmousi模型进行试算,可以显著节省内存需求,明显提高计算速度,进而得到精度较高的正演结果.【期刊名称】《石油地球物理勘探》【年(卷),期】2014(049)002【总页数】7页(P288-294)【关键词】频率—空间域;CSC;系数矩阵;一维数组;LU分解;数值模拟【作者】吕晓春;顾汉明;成景旺;周丽【作者单位】中国地质大学(武汉)地球物理与空间信息学院,湖北武汉430074;中国地质大学(武汉)构造与油气资源教育部重点实验室,湖北武汉430074;中国地质大学(武汉)地球物理与空间信息学院,湖北武汉430074;中国地质大学(武汉)构造与油气资源教育部重点实验室,湖北武汉430074;中国地质大学(武汉)地球物理与空间信息学院,湖北武汉430074;中国地质大学(武汉)构造与油气资源教育部重点实验室,湖北武汉430074;中国地质大学(武汉)地球物理与空间信息学院,湖北武汉430074【正文语种】中文【中图分类】P6311 引言为了反演地下介质参数或研究地震波在各种复杂介质中的传播机制，需要进行波场数值模拟。
波动方程的数值模拟方法包括有限差分法、有限元法、伪谱法等。
在时间—空间域的数值模拟技术已较成熟，并广泛应用于复杂介质的正演模拟中。
基于GID有限元前处理的波动方程数值模拟
基于GID有限元前处理的波动方程数值模拟刘静;文山师;黄晶晶【摘要】在地震波数值模拟计算过程中,缺乏简单易行的有限元前处理方法,使得复杂构造模型较难建立和分析.本文以二维声波方程为例结合GID软件,网格剖分部分采用三角形单元模拟速度界面,把单元内的场和波速均看作单元上的线性函数;GID 软件可以方便地进行网格剖分和设置网格控制节点,通过编写用户自定义”问题类型”,建立并输出已有的有限元计算程序的初始模型.将GID软件前处理与有限元计算程序整合,提高了方法的效率,简单易行.【期刊名称】《工程地球物理学报》【年(卷),期】2014(011)002【总页数】7页(P243-249)【关键词】数值模拟;有限元;GID;声波方程;三角形单元【作者】刘静;文山师;黄晶晶【作者单位】山西省煤炭地质115勘查院,山西大同037003;中石化西北油田分公司勘探开发研究院,新疆乌鲁木齐830011;中石化石油工程地球物理有限公司河南分公司,河南南阳473000【正文语种】中文【中图分类】P631.41 引言地震波场的数值模拟技术是在已知地下介质结构和参数的情况下，利用理论计算的方法研究地震波在地下介质中的传播规律，合成地震记录的一种技术。
地震勘探中的数值模拟方法主要以射线理论和波动方程理论为基础，有射线追踪法，柯西霍夫积分法，有限元法，有限差分法和伪谱法[1～6]。
有限差分法直接用差分代替微分，因其方法简单、精度高，在地震模拟中而得到了广泛的研究和应用。
但其固有缺陷是不能准确模拟具有复杂几何形态的物性界面，有限元法则是求解原问题等价泛函的变分或原问题的等效积分方程的弱解（当等价泛函不存在时），因而能够适应较有限差分更为剧烈的物性变化，加之种类繁多的插值形函数，使其能够模拟很复杂的几何界面。
有限元法的主要缺点是计算和存储量都很大，效率相对较低。
建立有限元分析模型比较复杂且存在困难，因此可以用一些成形的软件作为有限元网格剖分的工具，建立并输出可用于已有有限元计算程序的初始模型，将大大提高方法的效率[7]。
波动方程正演模拟边界条件的比较分析
波动方程正演模拟边界条件的比较分析付小波;韩超;原健龙;余嘉顺【摘要】通过数值模拟研究了透明边界、Clayton-Engquist边界和完全匹配层边界的吸收效果,得出如下结论:在反射角和频率相同的情况下,完全匹配层边界条件效果最好,Clayton-Engquist边界效果次之,而透明边界条件的效果最差.以边界条件对100 Hz模型边界垂直反射的吸收效果来衡量,完全匹配层边界条件与Clayton-Engquist边界条件的效果分别是透明边界条件的16.5倍和3.5倍.在＜40 Hz的低频范围内,或者在反射角＞65°的情况下,Clayton-Engquist边界相对透明边界的吸收效果相对优势显著变弱.而完全匹配层边界的吸收效果则在150 Hz频率范围内和75°反射角范围内始终保持稳定的相对优势.【期刊名称】《成都理工大学学报（自然科学版）》【年(卷),期】2015(042)004【总页数】8页(P492-499)【关键词】波动方程;正演模拟;有限差分;边界条件【作者】付小波;韩超;原健龙;余嘉顺【作者单位】成都理工大学地球物理学院,成都610059;成都理工大学地球物理学院,成都610059;成都理工大学地球物理学院,成都610059;成都理工大学地球物理学院,成都610059;新西兰皇家地质与核科学研究所,惠灵顿【正文语种】中文【中图分类】P631.4边界条件是地震波数值模拟方法技术中的一项重要内容。
许多专家学者从不同角度提出多种构造边界条件的方法。
1977年，Clayton与Engquist［1］根据旁轴近似理论（Claerbout［2］；Claerbout与Johnson［3］），提出利用一系列不同近似精度的单程波动方程来吸收模型边界的反射能量，称作Clayton-Engquist （下文采用简略记号CE来表达）边界条件。
1978年，Reynolds通过对波动方程的分解得到了透明边界条件［4］（下文采用简略记号TBC来表达）。
波动方程及其解法
波动方程及其解法波动方程是常见的偏微分方程之一，它描述的是波的传播和变化。
而在实际问题中，如声波、光波、电磁波等的研究中，波动方程的解法是被广泛使用的。
本文将介绍波动方程的基本概念及其解法。
一、波动方程的基本概念波动方程最基本的形式是一维波动方程，其数学表达式如下：$\frac{\partial^2 u}{\partial t^2}=c^2\frac{\partial^2 u}{\partial x^2}$其中，$u(x,t)$表示波的位移，$c$是波的速度。
可以看出，波动方程是一个描述时间和空间之间关系的方程。
在这个方程中，偏微分算子表达了波动的传播和变化的规律。
二、波动方程的解法1. 分离变量法分离变量法是解波动方程的最常见方法之一。
其主要思想是，将变量$x$和$t$分离出来，分别让它们满足不同的微分方程。
如一维波动方程可以假设其解为$u(x,t)=X(x)T(t)$，将其代入波动方程可得：$XT''=c^2X''T$进一步变形，可得：$\frac{T''}{c^2T}=\frac{X''}{X}$由此得到两个方程：$\frac{T''}{c^2T}=-\omega^2$$X''=-\omega^2X$其中，$\omega$为角频率，$-\omega^2$为分离出来的常数倍。
对于这两个微分方程，可以分别求解。
2. 叠加原理在叠加原理中，可以将波看做是多个波的叠加。
这种方法可以用于特定场合下的波动方程求解。
例如，在弹性绳的研究中，可以将弹性绳的振动看作是多个波的叠加。
在这种情况下，可以对不同的波求解，并把它们的解加起来成为最终的解。
3. 直接积分法直接积分法是一种基本的解微分方程的方法，同样也适用于波动方程的求解。
在直接积分法中，可以通过对波动方程进行积分，逐步求解出波的变化规律。
这种方法的实现需要考虑初值条件的限制，而条件的不同可能导致问题的复杂性。
2007射线追踪与波动方程正演模拟方法对比研究
由于地震波在整条路径上满足同一个射 线参数，因此射线路径上任意连续三点也将满 足同一个参数， 而三点间的射线表现形式为 Ｓｎｅｌｌ 定律。按照 Ｓｎｅｌｌ 定律，可导出一个求
２ 设计依据：
１ ） 根据《混凝土结构加固技术规范》
ＣＥＣＳ１４６：２００３，设计图纸和该根据工程检测 报告编号 ＢＯｂＬＯＪＧ０３３，本工程采用加大截面 加固法、外粘钢加固法等， 其工作程序如下：
可靠性稳定→加固方案→加固设计→施 工→验收
２）材料：外包钢采用 Ｑ２３５ 材料 Ｌ８０ × ８０ × ５ 的 Ｂ 型角钢， Ａ ｓ ＝ ７ ９ １ ｍ ｍ ２、ｆ ｙ ＝ ２ １ ５ Ｎ ／ ｍｍ２ 加固箍筋用扁铁 ４０ × ４，外包混凝土用 Ｃ ２ ５ （ ｆ Ｃ ＝ １ １ ． ９ Ｎ ／ ｍ ｍ ２） ， 外包钢加固后的尺寸 ｂ × ｈ＝５００ × ７００，角钢与扁钢的连接采用焊 接， ｈ ｆ ＝ ５ ｍ ｍ 焊缝饱满， 焊条 Ｅ ４ ３ 、Ｅ α ＝ ２ ． １ × １０５Ｎ／ｍｍ２，粘钢采用改性环氧树脂胶粘 剂。
１ 基于射线追踪的合成地震响应
射线追踪法的主要理论基础是，在高频近
似条件下， 地震波的主能量沿射线轨迹传播。 基于这种认识，运用惠更斯原理和费马原理来 重建射线路径，并利用程函方程来计算射线的 旅行时。在旅行时计算中应用有限差分等方 法， 以获得快速的解。射线法的主要优点是 概念明确，显示直观，运算方便，适应性强；其 缺陷是应用有一定限制条件，计算结果在一定 程度上是近似的，对于复杂构造进行两点三维 射线追踪往往比较麻烦。为了计算波沿射线 的旅行时和波的传播路径， 叙述如下。
用于波动方程模拟的Chebshev谱元法
１引言
在瞬态分析，工程地震学，计算声学等领域 （ 如无损检测，油气勘探等） ，如何利用数值计算的方法更
精确地得到弹性波动方程的解一直是国内 外研究者的工作重点。 着计算机技术的发展， 随 一些原来影响数值
计算方法应用的瓶颈一一被克服， 但对于大型的复杂二维问题或三维问题的研究， 仍然对原有的数值模拟方 法提出了挑战。 当前普遍使用的数值方法，如有限体积法 （ｉｔＶｌ ｅ ｔ ｄ，有限元 （ｉｔＥｅ ｅｔ ｈｄ， Ｆ ｉ ｏ ｍ Ｍｅ ｏ） ｎｅ ｕ ｈ Ｆ ｉ ｌ ｎＭｅ ｏ） ｎｅ ｍ ｔ
Ｃｅ ｈ 正 多 式 ｅ ｎｒ多 展开。３ 伽 金 法 解正 题的 分 式， 全 近 ｈ ｓ ｖ 交 项 或Ｌｇ ｄ 项式 ｂｅ ｅ ｅ （ 用 辽 方 求 交问 变 格 得到 局的 似 ）
解。 有关谱元法的详细数学表述请参看文献 ７ ０ 我们这里采用 Ｃｅｓｅ 正交多项式，它是如下奇异性 Ｓ ｒ－ｉ ｖｌ方程的特征函数 ｈｂｈｖ ｔｍＬｏ ｉ ｕ ｕｌ ｅ
似函数能最佳地逼近偏微分方程的精确解，测试函数 （ｅ Ｆｎｔｎ Ｔｓ ｕｃｏ）被引进用于验证近似解带来的余量是 ｔ ｉ 否达到最小。对基函数和测试函数的不同选择导致了上述这几种数值方法。
２ ｅｓｅ 谱元法 Ｃ ｂｈｖ ｈ
谱元法（ Ｅ ） 最早 Ｍ ， 在由Ｐｔａ ‘ 并 应用于流体动力学。 把有限 （ Ｓ Ｐ ａｒ提出２ 主要 ｅ ］ ， 它 元法和谱方法相结 合，
Ｃｅ ｓ ｖ ａ ｓｏ ｔ 配置点 权重 ｈ ｙ ｅＧｕ－ ｂ ｏ ｂ ｈ － ｓＬ ａ 及其 定义如 下，
（） ２
ｘｏ， ｉｓ 二务 Ｃ
ｗ二 ：—
ｒ 一 ， ， ７ ： 、 ‘ ． ， 一二 。
ｊ０Ｎ ｗ 二 ＝， ； ； －
波动方程与热传导方程的解法
波动方程与热传导方程的解法波动方程与热传导方程是物理学中常见的偏微分方程，它们描述了波动和热传导的过程。
在实际问题中，解这两个方程可以帮助我们了解和预测物理现象，例如声波传播、电磁波传播和热量传导等。
本文将介绍波动方程和热传导方程的解法及其应用。
一、波动方程的解法波动方程描述了波的传播和干涉。
通常表示为：∂²u/∂t² = v²∇²u其中，u代表波的振幅，t代表时间，v代表波速，∇²u是u的拉普拉斯算子。
1. 分离变量法分离变量法是求解偏微分方程的常用方法。
对于波动方程，我们可以假设u(x, t)的解为u(x, t) = X(x)T(t)，其中X(x)和T(t)是仅与x和t相关的函数。
将u(x, t)的表达式带入波动方程，我们可以得到两个关于X(x)和T(t)的普通微分方程。
通过求解这两个方程，我们可以得到波动方程的解。
2. 傅里叶变换法傅里叶变换法也是求解偏微分方程的重要方法。
通过将波动方程进行傅里叶变换，我们可以将其变换为关于频率和空间变量的代数方程，进而求解得到波动方程的解。
二、热传导方程的解法热传导方程描述了热量在物质中的传导过程。
通常表示为：∂u/∂t = α∇²u其中，u代表温度分布，t代表时间，α代表热扩散系数，∇²u是u 的拉普拉斯算子。
1. 分离变量法与波动方程类似，热传导方程也可以通过分离变量法求解。
我们可以假设u(x, t)的解为u(x, t) = X(x)T(t)，其中X(x)和T(t)是只与x和t有关的函数。
将u(x, t)的表达式带入热传导方程，我们可以得到两个关于X(x)和T(t)的普通微分方程。
通过求解这两个方程，我们可以得到热传导方程的解。
2. 球坐标系或柱坐标系下的解法对于具有球对称性或柱对称性的问题，我们可以将热传导方程转换为径向方程和角向方程，并通过求解这些方程得到热传导方程的解。
三、波动方程和热传导方程的应用波动方程和热传导方程广泛应用于物理学、工程学和其他领域中。
波动方程数值模拟技术及其应用
波动方程数值模拟技术及其应用作者姓名: 陈睿专业班级: 2008050603指导教师: 熊晓军摘要波动方程数值模拟技术在地震勘探中的应用非常广泛，特别是对于碳酸盐岩这一类重要的油气储集层。
本文主要介绍了声学波动方程的基本理论，相位移波动方程数值模拟方法，相位移加插值波动方程数值模拟方法的原理，并且采用相位移加插值的方法进行实际碳酸盐岩模型的数值模拟，根据实际区域的地质剖面猜测初始的地震模型，通过波动方程对该猜测的初始模型进行正演与偏移，再把通过偏移的地震剖面与实际的地震记录剖面对比，反复调整其中的相关参数，更新地质剖面，从而获得更加正确的地质解释模型。
对比地质模型与原始的地震资料，从而确定了猜测的正确性，为该地区以后的储层预测、地震资料解释提供了一定的参考价值。
综上的论述，本次研究为相同地震、地质条件下礁滩储层的波场特征认识积累了一些经验，为准确地进行礁滩储层预测奠定了一定的基础。
关键词：相位移波动方程数值模拟偏移Numerical Simulation Technology Of Wave Equation And Its ApplicationAbstract:The numerical simulation of wave equation is widely used in seismic exploration.Especially important to carbonate oil and gas reservoir.This paper introduces the basic theory of the acoustic wave equation, the phase shift of the wave equation numerical simulation method, the phase shift plus interpolation wave equation numerical simulation of the principle, and the phase shift plus interpolation, numerical simulation model of the actual carbonate, according to the geological profile of the actual region to guess the initial seismic model, forward modeling and migration by the wave equation of the initial model, and then offset seismic profiles with the actual seismic record section contrast, which repeatedly adjust the relevant parameters update the geological section, to obtain a more accurate geological interpretation model. Comparing the geological model and the original seismic data, in order to determine the correctness of the speculation, after the regional reservoir prediction, seismic data interpretation to provide a certain reference value.Comprehensive discourse on this study for the same earthquake, geological conditions, the wave field characteristics of the reef reservoir understanding gained some experience, and laid a foundation for the reef reservoir prediction. Keywords：Phase shift Wave equation Numerical Simulation Offset目录摘要 (I)第1章前言 (1)1.1 研究背景及意义 (1)1.2 研究内容 (1)1.3 研究方法 (1)第2章波动方程数值模拟的基本理论 (2)2.1 声学波动方程的基本理论 (2)2.1.1运动方程和应力位移方程 (2)2.1.2声学近似方程 (2)2.2 波动方程数值模拟的方法原理 (4)2.3 相位移波场延拓方法 (5)2.4 相位移加插值波场延拓方法 (7)2.5 点脉冲的实验 (12)2.5.1原理 (12)2.5.2实际中的问题及其解决办法 (13)2.5.3自激自收点脉冲的实验 (15)第3章碳酸盐岩模型的数值模拟 (16)3.1 二维模型的建模方法 (16)3.1.1二维地质模型描述 (16)3.1.2建立地质模型所面临的问题 (18)3.2 数值模拟的计算流程 (19)3.3 实例计算 (21)结论及建议 (24)致谢 (25)参考文献 (26)第1章前言1.1 研究背景及意义波动方程数值模拟技术在地震勘探中的起着重要作用。
python 波动方程
python 波动方程波动方程是描述波动现象的数学模型，它可以用来描述各种波动现象，如声波、光波、电磁波等。
在物理学、工程学、数学等领域中，波动方程都有着广泛的应用。
在本文中，我们将以Python为工具，介绍波动方程的基本概念和求解方法。
波动方程的基本概念波动方程是一个偏微分方程，它描述了波动现象的传播和变化。
波动方程的一般形式为：∂²u/∂t² = c²∇²u其中，u是波动的位移，t是时间，c是波速，∇²是拉普拉斯算子。
这个方程的意义是，波动的加速度与波动的曲率成正比。
也就是说，波动的加速度越大，波动的曲率也越大。
求解波动方程的方法求解波动方程的方法有很多种，其中比较常用的方法是有限差分法和有限元法。
有限差分法是一种数值求解方法，它将波动方程离散化为一组差分方程，然后通过迭代求解差分方程来得到波动的解。
有限差分法的优点是简单易懂，计算速度快，但是精度相对较低。
有限元法是一种更加精确的数值求解方法，它将波动方程离散化为一组有限元方程，然后通过求解有限元方程来得到波动的解。
有限元法的优点是精度高，但是计算速度相对较慢。
Python求解波动方程的实现在Python中，我们可以使用NumPy和Matplotlib等库来实现波动方程的求解和可视化。
下面是一个简单的Python程序，用于求解一维波动方程：import numpy as npimport matplotlib.pyplot as plt# 定义常量c = 1.0dx = 0.1dt = 0.01L = 10.0T = 10.0# 定义初始条件x = np.arange(0, L+dx, dx)u0 = np.exp(-0.5*(x-5)**2)# 迭代求解差分方程u = u0.copy()for n in range(int(T/dt)):u[1:-1] = u[1:-1] - c**2*dt/dx**2*(u[2:]-2*u[1:-1]+u[:-2])# 可视化结果plt.plot(x, u0, 'k--', label='Initial')plt.plot(x, u, 'r-', label='Final')plt.legend()plt.show()在这个程序中，我们首先定义了一些常量，如波速c、空间步长dx、时间步长dt、空间范围L和时间范围T。
波动方程ppt课件
＝
2π
2π
d ＝Cd
C (本题结束)
判断各点运动 方向的技巧
上坡下行
例题：有一列横波向右
下坡上行
传播， 画出波形曲线上 A、B 、C 、D 、E 、F 各 点的运动方向和四分之
y C· B· ·D
u
一周期后的波形曲线。
· A 0
T 4
E·
·F
x
特别要注意：波的传播方向，这是关键。
例题：图（a）中所表示的x ＝0 处质点振动的初相位
y(m) 0.04
0
－0.04
u＝0.08 m/s
.a
b.
0.2
0.4
x (m)
例题：一列沿x 正向传播的简谐波，已知t1=0和
t2=0.25s时的波形如图。
试求： （1）振动方程 （2）波动方程 （3）作出波源振动曲线。
（练习册P32计算题3·版书）
y(m)
u
0.02m
t1 t2
..
.
0
P
x (m)
置的位移（坐标为 y）随时间t 的变化关系。
波波
函 数
动 方 程
y
=A
cos ω
(
t－
x
u
)
+j
波向x 轴正方向 传播也称右行波
波向x 轴负方向 传播也称左行波
y
=A
cos
ω
(
t
＋
x
u
) +j
物理意义：波线上任一点（距原点为 x）处 的质点任一瞬间相对其平衡位置的位移。
当波向x 轴正方向传播而且已知 距离0点为xo的Q点振动方程为：
与图（b）所表示的振动的初相位分别为：
波动数值模拟的常加速度显式算法
波动数值模拟的常加速度显式算法孙明社;曲淑英;侯兴民【摘要】给出了一种将隐式时域逐步积分算法转换为显式时域逐步积分算法的方法，避免了求解耦联方程组，提高了计算效率.对于暂态波源作用下的弹性半空间，利用有限单元法划分网格、建立结构动力方程，并应用Fortran语言对中心差分算法和平均常加速度显式算法编程，求解脉冲荷载作用下的出平面运动.2种算法计算结果对比表明，常加速度显式算法可以较好地应用于工程波动数值模拟中.%Explicit time integral method of numerical simulation for engineering wave problem is an important topic in both national and international research work. An explicit time integral method transformed by corresponding im-plicit time integral method, avoiding solving coupled equations and improving computational efficiency, is presented in this paper. Structural dynamic equation of elastic half space under transient wave is integrated using the finite el-ement method, and the out-plane response under impulse load is computed by central difference method and con-stant acceleration explicit method respectively. Comparison between the calculated results of the two methods shows that the constant acceleration explicit method can be used in engineering wave numerical simulation effectively.【期刊名称】《烟台大学学报（自然科学与工程版）》【年(卷),期】2015(000)001【总页数】5页(P61-65)【关键词】常加速度显式算法;中心差分算法;出平面运动;工程波动【作者】孙明社;曲淑英;侯兴民【作者单位】烟台大学土木工程学院，山东烟台264005;烟台大学土木工程学院，山东烟台264005;烟台大学土木工程学院，山东烟台264005【正文语种】中文【中图分类】TU470在土木工程和地震工程领域有很多问题都可以归结为波动问题，如强震地面运动、土－结构相互作用、结构的无损检测与探伤等问题.因而，研究工程波的传播与振动具有非常重要的意义.通过建立力学模型，工程波动问题一般就可转化为偏微分方程的求解.由于实际工程中许多工程结构建立的动力方程为非线性，很难获得精确的解析结果，因而，数值方法便得到广泛应用.时域逐步积分法是结构动力方程求解中的一种有效方法.根据是否需要求解耦联方程组，时域逐步积分法又可分为显式算法和隐式算法.隐式时域逐步积分算法的研究成果较多，如常平均加速度方法、Houblt方法、Newmark方法、Gurtin方法、Wilson-θ方法、Park方法以及 a方法等[1－5].但是随着结构自由度数目的增大，求解耦联方程组的工作量巨大，由于显式时域逐步积分方法无需求解耦联方程组，因而具有明显的优势.廖振鹏、李小军、周正华、侯兴民等[6－20]对显式时域逐步积分算法都进行了研究.本文在国内外学者研究成果的基础上探讨了平均常加速度隐式算法，即通过矩阵级数展开转变为显式算法，最后将该显式算法在工程波动数值模拟中进行应用，并与文献[21]的中心差分算法比较，说明常加速度显式算法可以较好的应用于工程波动数值模拟中.1 隐式算法显式化推导1.1 基本方法对于时域逐步积分方法的广义线性加速度算法，每向前计算一步都需求解耦联线性方程组式中:M为质量矩阵;C为阻尼矩阵;K为刚度矩阵;{u}i+1为ti+1时刻的位移;{}i+1为ti+1时刻的广义荷载;a，b为常数，不同的加速度取不同的值.将式(1)中左侧{u}i+1的系数矩阵提出矩阵，可得对式(2)矩阵求逆矩阵，有将式(1)进行矩阵求逆、级数展开、级数截断最终得到式(9)，完成了隐式算法转化为显式算法.式(9)可不必求解耦联线性方程组，提高了计算速度，节省了计算时间.用上述方法，以平均常加速度方法为例，推导其隐式显式化过程，并进行出平面运动的波动数值模拟.1.2 平均常加速显式化当Newmark－β逐步积分方法中的参数γ、β取不同值时，该方法退化为平均常加速度法、线性加速度法、中心差分法等.当γ=0.5、β=0.25时，为平均常加速法，即式(1)中的参数取值:a=4，b=2.公式(1)可记作将式(10)与式(1)对比，应用前述推导过程，可以得到平均常加速方法对应的显式化公式如下应用式(13)、(14)可以分别得到离散时间点上的速度、加速度式(12)、(13)、(14)为平均常加速算法的显式化公式，若给定初始时刻的初始位移{u}0、初始速度{˙u}0、初始加速度{¨u}0，即循环使用公式(12)、(13)、(14)，得到所有离散时间点上的位移、速度、加速度.平均常加速方法显式化的收敛条件，由式(7)可知取向量的“1”范数，求解上式，可得对于单自由度无阻尼体系ξ=0，带入上式，得2 出平面波动数值模拟2.1 出平面问题如图1，在直角坐标系oxyz下，弹性半空间内，剪切模量μ1=1，介质密度ρ1=1，不考虑阻尼.波源为沿y轴方向深度hs处作用的暂态荷载:图1 出平面波源问题Fig.1 Problem of out-plane wave暂态波源荷载时间分布函数为三角形脉冲T(t):2.2 平面波动数值模拟波动数值模拟的精度需要将有限单元的尺寸划分的与波长相比足够小，因而在建立动力方程时才可以将每一有限单元内惯性体积力视作不随空间变化的恒定量.为了实现动力方程的空间解耦，采用了集中质量有限元模型.采用有限单元法将连续弹性半空间介质离散化，用有限个单元组合体代替原本的连续介质.用平面正方形单元将弹性半空间介质离散化，假定各个单元只在公共节点上相互铰接.离散的正方形单元边长即空间离散步长为:Δx=Δy=0.05;单元节点编号如图2;单元刚度[15]如矩阵(23).图2 正方形单元节点编号Fig.2 Node number of square unit将弹性半空间沿y轴向下划分n个网格;y轴左右两侧，沿x轴正负方向各划分m个网格，见图3.图3 弹性半空间网格离散Fig.3 Discrete mesh of elastic half space观测点位移峰值如表1.有限元网格划分m=40，n=40;时步Δt=0.0075 s;暂态波源作用位置hs=0即表面荷载;考虑二阶透射边界.采用Fortran程序语言进行中心差分算法和平均常加速隐式转换显式算法编程计算，取(0，0)(0.5，0)(1.0，0)(0，0.5)(0，1.0)5个观测点，得到暂态波源作用下的弹性半空间计算结果.相同时步下中心差分算法和常加速度显式算法不同观测点的位移曲线如图4.表1 观测点位移峰值Tab.1 Displacement peak of observation point观测点中心差分算法常加速度显式算法(0，0) 0.346 0.346(0.5，0) 0.194 0.194(1.0，0) 0.137 0.137(0，0.5) 0.223 0.223(0，1.0) 0.175 0.175图4 观测点位移曲线Fig.4 Displacement curve of observation point将式(12)中系数逆矩阵取不同展开阶数下的各观测点位移曲线如图5.图5 不同阶数下的观测点位移曲线Fig.5 Displacement curve of observation point under different order number3 结论通过对比中心差分和常加速显式方法的计算结果分析，得到如下结论.(1)相同的网格数目、时步、边界条件下，2种算法的出平面波动的主要波形完全吻合，峰值相同.曲线后段虽然存在一定的差别，但是，工程中我们主要关注主要波形，因而可以满足工程要求.(2)考虑阻尼的情况下，中心差分算法的显式优势不明显，而常加速度显式算法仍为高效率的显式算法，具有一定优势.(3)逆矩阵在单位矩阵处级数展开，阶数不同，但观测点的位移曲线基本重合，说明常加速显式算法只需取前几阶就可满足精度要求，且计算效率较高.参考文献:[1]Newmark N M.A method of computation for structural dynamics[J].Proc ASCE，1959，85(3):69－94.[2]Bathe K J，Wilson E L.Stability and accuracy analysis of direct integration method[J].Earthq Eng and Struct Dynam，1973(1):283－291. [3]Hiber H M，Hughes T J R，Tayler R L.Improned numerical dissipationfor time integration algorithins in structural dynamics[J].Eathq Eng and Struct Dynam，1977，5(3):283－292.[4]Hiber H M，Hughes T J R.Collocation dissipation and overshoot for time integration schemes in structural dynamics[J].Eathq Eng and Struct Dynam，1978，6(1):99－177.[5]Subbara J K，Dokainish M A.A survey of direct time integration methods in computational structural dynamics(II)[J].Computers and Structures，1989，32(6):1387－1401.[6]廖振鹏.近场波动问题的有限元解法[J].地震工程与工程振动，1984，4(2):1－14.[7]李小军，廖振鹏，杜修力.有阻尼体系动力问题的一种显式差分解法[J].地震工程与工程振动，1992，12(4):74－79.[8]钟万勰.结构动力方程的精细时程积分法[J].大连理工大学学报，1994，34(2):131－136.[9]李小军.地震工程中动力方程求解的逐步积分方法[J].工程力学，1996，13(2):110－118.[10]周正华，李山有，侯兴民.阻尼振动方程的一种显式直接积分方法[J].世界地震工程，1999，15(1):41－44.[11]杜修力，王进廷.阻尼弹性结构动力计算的显式差分法[J].工程力学，2000，17(5):37－43.[12]王进廷，杜修力.有阻尼体系动力分析的一种显式差分法[J].工程力学，2002，19(3):109－112.[13]陈学良，金星，陶夏新.求解加速度反应的显式积分格式研究[J].地震工程与工程振动，2006，26(5):60－67.[14]高小科，邓子辰，黄永安.基于三次样条插值的精细积分法[J].振动与冲击，2007，26(9):75－82.[15]廖振鹏，刘恒，谢志南.波动数值模拟的一种显式方法－一维波动[J].力学学报，2009，41(3):350－359.[16]刘恒，廖振鹏.波动数值模拟的一种显式方法－二维波动[J].力学学报，2010，42(6):1104－1116.[17]谢志南，廖振鹏.波动方程数值模拟的一种显式方法[J].力学学报，2011，43(1):154－161.[18]李长青，楼梦麟，余志武，等.近似平衡多项式加速度动力显式算法[J].应用力学学报，2011，28(5):475－480.[19]李长青，楼梦麟.中心－偏心差分法[J].同济大学学报:自然科学版，2011，39(2):179－186.[20]李长青，楼梦麟，蒋丽忠.结构动力方程求解中隐式格式向显式格式的转换[J].振动与冲击，2012，31(13):91－94.[21]廖振鹏.工程波动理论导论[M].北京:科学出版社，2002.。
三维波动方程有限差分正演方法
三维波动方程有限差分正演方法摘要：本文主要介绍了三维波动方程有限差分正演方法的基本原理和实现过程，并对其优缺点进行了分析。
通过数值实验验证了该方法的可行性和准确性，为地震勘探、地下水文学等领域的研究提供了参考。
关键词：三维波动方程、有限差分、正演方法、数值实验一、引言三维波动方程是地震勘探、地下水文学等领域的基础理论，其求解方法对于相关领域的研究具有重要意义。
有限差分正演方法是求解三维波动方程常用的数值方法之一，其具有计算速度快、精度高等优点，在实际应用中得到了广泛的应用。
本文将介绍三维波动方程有限差分正演方法的基本原理和实现过程，并对其优缺点进行分析，同时通过数值实验验证该方法的可行性和准确性。
二、三维波动方程有限差分正演方法的基本原理三维波动方程可以表示为：其中，u(x,y,z,t)为波动场，c(x,y,z)为介质速度，ρ(x,y,z)为介质密度，t为时间。
有限差分正演方法通过将空间和时间离散化，将三维波动方程转化为差分方程，进而求解波动场在不同时刻的数值解。
具体而言，有限差分正演方法将空间和时间分别离散化，将空间网格点和时间网格点相结合，构成一个四维网格空间，其中每个网格点对应一个波动场的数值解。
有限差分正演方法的求解过程可以分为以下几个步骤：1. 空间离散化对于三维空间，可以将其分为x、y、z三个方向，分别进行离散化。
假设在x方向上，空间步长为Δx，在y方向上，空间步长为Δy，在z方向上，空间步长为Δz。
则可以将空间网格点表示为(xi,yj,zk)，其中i=1,2,...,nx，j=1,2,...,ny，k=1,2,...,nz，nx、ny、nz分别为x、y、z方向的网格数。
2. 时间离散化假设时间步长为Δt，则在t时刻波动场的数值解可以表示为ui,j,k^n，其中n表示时间步数，i、j、k表示空间网格点的索引。
3. 有限差分近似将三维波动方程中的导数项用有限差分近似表示，例如：其中，ui,j,k^n表示在t时刻，在(xi,yj,zk)处的波动场数值解，c(xi,yj,zk)表示在(xi,yj,zk)处的介质速度，ρ(xi,yj,zk)表示在(xi,yj,zk)处的介质密度，Δx、Δy、Δz、Δt分别为空间和时间步长。
基于OpenMP+SIMD的波动方程数值模拟并行方法
基于OpenMP+SIMD的波动方程数值模拟并行方法曹丹平5 10 15 20 25 30（中国石油大学（华东）地球科学与技术学院,青岛266580）摘要：地震勘探中的大规模并行计算在科研与生产中发挥了重要的作用，但是大规模并行中基于MPI 和GPU 的并行方法对编程的要求较高。
本文以多核处理器中的小规模并行计算为例，在多核微机处理器上将支持多线程并行的OpenMP 方法与支持数据并行的SIMD 指令相结合，针对二维波动方程在时间和空间方向上的循环嵌套特点同时实现了基于OpenMP+SIMD 并行的波动方程数组模拟方法。
模型测试表明在普通的双核笔记本电脑上采用OpenMP+SIMD 的并行方法即可达到6 倍的加速比，从而在常规多核处理器上实现了计算效率的有效提高，同时为进一步提高大规模并行计算的效率奠定了基础。
关键词：波动方程；数值模拟；多核并行；单指令流多数据流中图分类号：P631.4Parallel computing of wave equation numerical modelling based on the OpenMP + SIMD techniqueCAO Danping(School of Geosciences, China University of Petroleum (Huadong), Qingdao 266580) Abstract: The large scale parallel computing is very important in the seismic exploration industry and scientific research, which is a very professional skill for the MPI and GPU programming. The small scale parallel computing problem is considered in this paper, two parallel computingmethods are combined to improve the compute efficiency. The OpenMP method that s upport the multithread parallel computing is used to manage the multicore resources of the CPU, the SIMD (Single Instruction, Multiple Data) method that support the data parallel computing is used to compute more data with a single CPU instruction. OpenMP and SIMD are combined together to compute the wave equation parallel, the cycle characteristics of wave equation in time step andspace step are considered. The model test shows that the parallel computing method ofOpenMP+SIMD improved the compute efficiency significantly, the speedup of dual coreprocessor is up to 6, which is very important to enhance the large scale parallel computing.Key words: Wave equation; numerical modelling; single instruction multiple data;parallel computationmulticore35 40 0引言开展波动方程数值模拟有助于准确认识和掌握油气勘探中的地震资料反射特征，同时正演模拟也是全波形反演和逆时偏移方法的重要组成部分，但是波动方程数值模拟存在耗时较大的问题，如何提高波动方程的正演模拟效率对于深入开展地震勘探的相关方法研究具有重要意义。
FDTD使用说明文档
FDTD使用说明文档FDTD（Finite-Difference Time-Domain）是一种计算电磁波动方程的数值模拟方法。
它通过将空间和时间离散化，将整个问题转化为了差分方程的求解。
FDTD方法适用于计算二维和三维空间中的电磁波的传播和辐射问题，广泛应用于大气物理、电磁学、光学和电磁兼容等领域。
下面是FDTD的使用说明文档，包括基本原理、步骤和参数设置等。
一、基本原理：FDTD方法基于麦克斯韦方程组，将空间和时间划分为网格进行离散化，通过差分形式的麦克斯韦方程进行求解。
具体步骤如下：1.空间离散化：将计算区域划分为网格，每个网格点上都有电场和磁场分量。
2.时间离散化：使用时间步长Δt，将时间进行离散化。
3.更新电场：根据麦克斯韦方程组的电场更新公式，根据磁场的值更新电场的值。
4.更新磁场：根据麦克斯韦方程组的磁场更新公式，根据电场的值更新磁场的值。
5.边界条件：设置适当的边界条件，如吸收边界条件、周期性边界条件等。
6.重复步骤3-5，直到模拟结束。
二、步骤：使用FDTD方法进行模拟一般可分为以下步骤：1.设定计算区域的大小和网格划分，根据模拟需求确定网格节点数和间距。
2.初始化电场和磁场，设置初始场分布。
3.根据模拟需求设置时间步长Δt，以及计算的总时间或模拟步数。
4.迭代更新电场和磁场，按照FDTD的原理进行计算。
5.设置边界条件和吸收边界条件，确保计算区域的边界不会对计算结果产生影响。
6.输出结果，根据需求选择输出电场、磁场以及网格中其他物理量的数值。
7.模拟结束。
三、参数设置：在使用FDTD方法进行模拟时，一些重要的参数需要进行合理的设置，以保证模拟结果的准确性和稳定性：1.网格分辨率：根据模拟的需求和计算资源，设置合适的网格划分和节点数，以充分捕捉到目标问题的细节。
2.时间步长：时间步长Δt决定了模拟的时间分辨率，需要根据模拟的频率范围和计算精度要求设置。
3.边界条件：选择适当的边界条件，可以是吸收边界条件、周期性边界条件等，以避免计算区域的边界对计算结果的影响。
波动数值模拟的稳定性
波动数值模拟的稳定性廖振鹏;谢志南【摘要】依据数值解收敛方式的不同首先将Lax稳定性区分为强稳定性和弱稳定性,进而将稳定性分析方法归纳为强和弱2种方法,并指出两者的关系.对前者回顾了谱分析和正则模态分析研究工作的主要进展,评论了在这一领域中Godunov和Ryabenkii的原创性工作以揭示其对局部稳定性研究的价值.初步论证了有界面弹性杆中波动数值模拟的经验稳定准则.对后者分别阐明了将其应用于有限差分和有限元模拟的前提,特别是用于论证有限元模拟的稳定性的价值.最后讨论了强、弱稳定性分析对进一步发展波动数值模拟技术的含意.%Lax-stability is first differentiated into the strong and weak categories according to different manners of convergence, techniques of the stability analysis are then classified into the strong and weak stability analysis and relationship between the two techniques is pointed out. For the former, important advances in the spectral analysis and the normal mode analysis are reviewed first, and comments on the original work by Godunov and Ryabenkii in this field are then made to reveal its value for studying the local stability. An elementary argument is presented for verifying the empirical stability criterion for numerical simulation of wave motion in an elastic bar with an interface. For the latter, its applicable premises are respectively clarified for simulation of finite differences and finite elements. In particular, its value in demonstrating the stability of the finite element simulation is stressed. Implication of the strong stability analysis and theweak one is finally discussed for further developing techniques of the numerical simulation of wave motion.【期刊名称】《哈尔滨工程大学学报》【年(卷),期】2011(032)009【总页数】8页(P1254-1261)【关键词】Lax稳定性;von Neumann分析;谱分析;正则模态分析;GKS定理;局部稳定性【作者】廖振鹏;谢志南【作者单位】中国地震局工程力学研究所,黑龙江哈尔滨150080;哈尔滨工程大学船舶学院,黑龙江哈尔滨150001;中国地震局工程力学研究所,黑龙江哈尔滨150080【正文语种】中文【中图分类】O241.82波动数值模拟是力学、电磁学、声学、地球物理和多个工程学科共同关注的领域，研究者在基础和应用研究方面皆取得了丰硕成果.但是，边界引入的失稳问题尚未彻底解决，特别是由边界引入的局部失稳的机理和消除方法，即使在线性范围内亦待进一步研究.这里“边界”泛指人工边界、物理边界或不同介质的分界面，“局部失稳”则指边界及其邻近空间区域在有限时间内由数值解格式不当所引发的误差放大.这一放大可能导致数值解灾难性发散而终止计算，亦可能仅在数值解中引入附加误差.例如，模拟辐射条件的人工边界格式与相邻内域节点的格式匹配不当就可能引发局部失稳［1-2］.迄今为止实用的稳定性分析方法为von Neumann分析［3-4］和谱分析［5］，前者从无限模型(含周期性边界模型)导出，未考虑非周期边界影响;后者从有限模型导出，在一定程度上考虑了边界影响，但未考虑边界可能引发的局部失稳.科技工作者通常采用经验和半经验的方法处理稳定问题［1-2］，数学家则对偏微分方程初边值问题数值解的稳定性作了深入的研究［6-11］.但是，数学家所采用的简化模型与各专业领域涉及的真实模型之间存在差距，同时，其分析方法与实际应用所要求的可操作性之间亦有差距.因此，应用后者提供的有益思想研究真实模型中发生的局部失稳问题，并由此发展实用的稳定性分析方法值得注意.本文试图在线性范围内考察和梳理稳定性的基本概念和主要研究方法，以利于进一步发展波动数值模拟技术，特别是解决人工边界引发的局部失稳问题.1 强稳定性和弱稳定性从收敛性和Lax等价定理出发，依据收敛方式的不同将Lax稳定性区分为强稳定性和弱稳定性.Lax等价定理是分析线性波动数值模拟问题的合理切入点［12］.该定理可表述为:如果连续模型的初边值问题适定和数值解格式相容，则收敛性等价于Lax稳定性.Lax稳定性涵义如下:就线性波动数值模拟的理论分析而言，数值解可写成式中:矩阵A为数值解算子，取决于内域场方程和边界条件的时空离散格式，因此，数值解方案的所有性质由A确定.A的阶数J与空间步距h相关，J阶向量vn和v0分别表示t=nk和t=0时刻的数值解，n为时步数，k为时间步距.Lax稳定性由下式定义:式中:T为计算时间长度，C为与k，h和T无关的常数，k0为某一正数，lnt表示取整.式(2)可表述为A的幂的范数一致有界.因此，稳定性分析是一个代数问题.稳定性分析的目的为确定条件式(2)对k和h取值所施加的限制，即k和h的取值应在k-h平面第一象限的某一区域之上.此区域称为稳定区，考虑到收敛性，稳定区以k-h平面坐标原点为聚点［11］.式(2)作为一个代数条件与时、空步距趋于零的方式无关，但是，从数值解收敛的角度分析，k和h趋于零的不同方式将导致条件式(2)本质上的差异.由此出发，这里依据k和h趋于零的2种不同方式将Lax稳定性区分为强稳定性和弱稳定性.1.1 强稳定性设在网格节点上连续模型的解为u(t)，强稳定性可以通过直接考察数值解vn对连续模型解un=u(nk)的收敛性引入.收敛性由下式定义:在波动数值解稳定性论著中通常引入参数Δτ=k/h，则稳定区可用参数Δτ和k描述.在数值解格式仅涉及参数Δτ的情况下，由稳定条件式(2)导出的稳定区简化为允许Δτ取值的实数区间，例如式中:Δτ*为稳定区的上界.对于稳定区内任一Δτ，h取决于k，故式(2)应在时、空步距同步趋于零时成立.数值解格式的强稳定性定义为在k和h同步趋于零的条件下式(2)成立.因此，数值解格式是否具有强稳定性，只需建立如式(4)所示类型的稳定准则，并检查Δτ=k/h是否满足此准则.由于k和h同步趋于零，当k→0时J→∞，即A的阶数随k→0而趋于无穷.由此可知，强稳定性涉及对一簇阶数变化的矩阵A的分析.对于模拟真实世界的模型，如何依据强稳定性确定稳定区，就作者所知，在一些重要情况下实用的分析方法尚待完善.1.2 弱稳定性弱稳定性可以通过分别考察数值解vn的空间和时间收敛性引入.首先对连续模型作空间离散得离散解w(t)所满足的常微分方程组.就稳定性分析而言，此方程组可以写成式中:J阶矩阵B为常微分算子.将式(5)作时间离散可导得式(1).于是，式(1)对于连续模型解的收敛问题简化为相互独立的2个问题:1)式(5)的初值问题的收敛性:2)时空离散模型的解vn对wn=w(nk)的收敛性.在保证前一收敛性的前提下，数值解格式的弱稳定性定义为:在给定h的条件下当k趋于零时式(2)成立.因为h已给定，A的阶数J与k无关，故弱稳定性分析仅涉及阶数不变的一簇算子A.依据弱稳定性确定稳定区已有成熟的谱分析方法.就保证数值解收敛而言，强、弱稳定性要求并无差别.两者的区别在于:1)前者直接分析时空离散格式，故可用于分析一般数值解格式可能出现的所有失稳现象，特别是局部失稳现象;后者则限用于分析离散模型中时、空变量相互独立的情形，特别是适用于波动的有限元模拟.2)就分析方法的可操作性而言，后者易于前者，且可利用有限模型的经典谱判据建立稳定准则，而前者在一些情况下建立稳定准则的实用方法尚待研究.3)就数值试验而言，强稳定性分析提供了稳定区间的界限，便于选取时间步距k;而弱稳定性保证Δτ→0时式(2)成立，为在数值模拟中通过逐步减小时间步距获得稳定的数值结果提供了依据.虽然弱稳定性难以为k的取值提供依据，这一缺点并不会成为数值试验的障碍，可以参考von Neumann分析结果或CEL条件试选k的实验数据［11］.下面进一步阐明这两类稳定性分析的基本概念及其演进.2 强稳定性分析回顾强稳定性谱分析和正则模态分析的主要进展，着重评论 Godunov和Ryabenkii的开创性工作［6］，特别是揭示在这项工作中隐含的对解决局部失稳问题有价值的思想.2.1 谱分析的进展及其局限性有限模型的谱分析建立在A的特征值和特征向量概念之上.A的全部特征值λ1，λ2…，λJ的集合{λi}称为A的谱，与谱{λi}对应的特征向量的集合记为{vi}.直接用A 的谱给出的稳定性判据称为谱判据.为了获得保证式(2)成立的谱判据可利用任一方阵的约当变换:式中:Q为J阶非奇异阵，Λ为A的约当标准型.假定式中:常数C与k和h无关，则得如下谱判据:1)对任意A，稳定的必要条件为稳定的充分条件为式中:谱半径ρ(A)=max|λi|.2)当A非缺损(non-defective)时，即当A有J个线性无关的特征向量时，稳定的充要条件为式(9).3)当A缺损(defective)时，稳定的充要条件为式(9)加上如下条件:{λi}中不含特征值λi，其模为1且代数重数(algebraic multiplicity)大于几何重数(geometric multiplicity)，即所有约当块皆为一阶.以上判据可称为经典谱判据.不仅当A非缺损时，而且当A缺损时，只要约当块的阶数有限，经典谱判据皆可用于判别强稳定性.经典谱判据的实用价值在于，对许多有意义的模型A为正规(normal)，此时使用这一判据的前提式(8)成立［5］.不过，对某些有意义的模型A为非正规，式(8)可能不成立并难以验证［6］.这是经典谱判据的局限性.为消除这一局限性，Kreiss和Buchanon使用酉变换替换约当变换将任一方阵A展开为［13］式中:U为J阶酉阵，Γ为J阶上三角矩阵:式中:Γ的对角元素按模的大小作升序或降序排列.由此可得保证式(2)成立的充要谱判据式中:K为与算子簇A无关的常数.式(13)提供了无前提条件的谱判据.上述谱判据可用于鉴别给定数值模拟方案的强稳定性.但是，有限模型的谱判据建立在排除所有内域节点控制方程与边界条件联合支持的失稳模态之上，因而不能揭示局部失稳的机理，即揭示某一节点与其邻近节点相互作用产生的失稳模态.因此，对于出现了局部失稳的数值模拟方案，谱分析无助于为改进方案提供线索和论证局部稳定性.这是有限模型谱分析的局限性.2.2 Godunov-Ryabenkii的工作和正则模态分析Godunov和Ryabenkii首先通过简单算例讨论了有限模型谱分析引入的矛盾，然后提出消除矛盾的如下思想:将有限模型分解为若干基本的无限模型，并用正则模态分析研究用谱分析难以处理的无限模型的稳定问题［6］.这一工作开启了稳定性研究的新阶段.为了进一步揭示这项工作的意义，特别是对研究局部失稳问题的重要性，本节首先对文献［6］所研究的算例作略为不同的分析，并介绍他们的和后续的主要研究结果;然后阐明此项成果为研究局部失稳问题所提供的启示.2.2.1 算例分析和相关研究结果考察单向波动边值问题:假定在x=0处连续模型的精确解已知.设，采用时、空截断误差皆为二阶的Friedrichs格式得到式(14)的离散方程［13］为就稳定性分析而言，式(15)的解可写成式(1)的形式，其中，A 的特征值为［14］由于A非缺损，利用式(9)得稳定准则式(4)，其中就此算例而言，不难证明Δτ*的正确值为1.这一不一致的结果表明，对于如式(16)所示的非正规矩阵A，依据经典谱判据导出的稳定准则不能保证式(2)成立.下面考虑式(14)的另一离散方案:利用单边差分公式作时空离散:若将式(18)的解写成式(1)的形式，其中，依据2.1节提供的缺损矩阵A的谱判据可得Δτ*＜2，但依据强稳定性要求得到正确结果是Δτ*=1，再次出现矛盾.为消除有限模型谱分析引入的矛盾，Godunov和Ryabenkii将上述有限模型分解为3个无限模型，即全无限和分别具有左、右端边界条件的2个半无限模型，并提出保证有限模型稳定的如下G-R条件:不允许这3个无限模型中的任一个具有特征值的模大于1的模态［6］.由于在数值模拟方案的常规设计中通常已保证von Neumann稳定条件，问题归结为保证半无限模型的稳定性.沿着这一方向后续研究工作表明G-R条件仍为稳定的必要条件，因为它未考虑当特征值的模为1时可能出现的失稳模态.为排除此模态需要补充一个条件，并针对一维一阶双曲型方程组证明了在一定附加条件下G-R条件加上补充的条件为稳定的充要条件，这个结果称为GKS定理［7-9］.GKS定理所补充的条件的涵义如下:边界条件和相邻内域节点运动方程的离散格式皆不支持能量从边界向内域流动的数值解，即离散格式不允许群速度从边界指向内域［10］.2.2.2 Godunov-Ryabenkii工作的评论1)研究人员通过数值实验早已观察到局部失稳源于边界与相邻内域节点的相互作用，意识到需要将这一相互作用分离出来加以研究，并试图实现这一分离以阐明失稳机理和制定消除局部失稳的措施［1-2，15］，但是，这些研究未能将这一“分离”的思想贯彻到底，未建立反映这一分离的基本数学模型，从而导致其研究结果的局限性.另一方面，在上述两例中算子A的非正规性皆来自内域离散格式的非对称性，即算例仅涉及内域稳定问题.实际上，若离散方案同时满足von Neumann稳定条件，则有限模型的谱判据并不引入矛盾.后继研究工作还表明，同时满足有限模型谱判据和von Neumann稳定条件在许多情况下可以保证数值模拟方案的稳定性［15-16］.但是，这补充了附加条件的谱分析亦不能揭示边界引入的局部失稳的机理.因此，Godunov等所研究的算例并未涉及本文提出的局部失稳问题.但是，他们随后提出的用基本无限模型替换有限模型以及用正则模态分析阐明其失稳模态的思想对指导真实模型局部失稳机理的研究及寻求消除局部失稳的数值模拟方案具有普遍价值.这是因为，1)Godunov和Ryabenkii提出的基本无限模型完全实现了将控制局部失稳的因素从众多影响因素中分离出来，使研究者的注意力集中到相互耦合的不多的几类节点运动方程上.消除局部失稳的尝试则归结为改变这些方程或其中之一的结构或系数，以避免形成失稳模态的影响因素之间的匹配.这为探寻稳定方案提供了中肯的线索.2)他们首先在稳定性研究中使用的正则模态分析则为阐明真实模型的局部失稳机理提供了技术思路.2)Godunov和Ryabenkii关于稳定性分析的原创性观点源于谱判据引入的矛盾.应当指出，他们所指出的矛盾并非谱判据本身的矛盾，而是未正确使用经典谱判据所引入的矛盾.如前所述，经典谱判据只能在式(8)成立的条件下使用.就上述2个例子而言，式(8)皆不成立.首先考察第1个例子.对此例难以直接验证式(8)，但可证明其不成立.若式(8)在Δτ＞0时成立，则谱判据式(9)为‖An‖一致有界的充要条件.因此，若谱判据式(9)成立，而‖An‖非一致有界，则式(8)不成立.下面证明‖An‖非一致有界.由于0≤n≤T/k和J=1/h=Δτ/k，当n≥J-2(此时T≥Δτ-2k)时，矩阵An右上角的元素为引入矢量v*=(0，0，…0，1)T，则Δτ在区间之内取值时，有当k→0和h→0时，由上式知‖An‖趋于无穷.因此，在引入矛盾的这一Δτ区间内，经典谱判据不成立.其次，看第2个例子.将式(18)式给出的A展开成式(7)的形式，则注意到 J随k→0而趋于无穷，对于在区间0＜Δτ＜1和1＜Δτ＜2内任一给定的Δτ，式(8)皆不成立.更有意思的是，应用无前提条件的谱判据式(13)可得稳定准则0＜Δτ＜1.因此，谱判据并未引入矛盾.不过，瑕不掩瑜，这一评论仅为科学的非逻辑发展补充一例.3 一维界面模型的稳定准则为阐明正则模态分析的基本思想，本小节应用这一分析方法初步论证文献［17］提出的一维分界面格式的经验稳定准则，并作数值检验.考虑无限弹性杆中剪切波的传播.设分界面位于x=0处，分界面两侧波动方程为式中:波速为剪切刚度以及密度，i=1，x＞0;i=2，x＞0.分界面连续条件和初值条件为设分界面右侧和左侧步距分别为h和(c2/c1)h，离散网格节点坐标为 xj=jh，j=0，1，2，…;xj=j(c2/c1)h，j=-1，-2，…(图 1).分界面及其邻近节点的递推公式(文献［14］式(15))，其二阶格式为式中Δτ=c1Δt/h，v(x，t)= ∂u/∂t.当 j=0 时式(23)给出界面节点的递推公式，当j≠0时给出界面两侧内域节点的递推公式.式(25)中参数p可以独立地选取而不影响公式的精度阶，但其取值与离散格式的稳定性相关.对于内节点递推公式［18］，依据 von Neumann稳定条件得到p=3/2.对于式(24)，文献［17］建议的经验准则为:对界面节点和两侧内域节点皆取p=3/2.图1 一维非均匀网格示意图Fig.1 One-dimensional non-uniform grid diagram下面应用正则模态分析对格式(24)的稳定性进行初步论证，即该格式不支持特征值的模超过1的如下模态:式中:φj(z)为非零模态向量，特征值z为复数.将式(26)代入式(24)并注意到式(25)给出的系数ξ1、η1、γ1、ζ1和 p=3/2 可得式中:差分方程式(27)的通解为式中:κ为复数，φ为二维非零常数向量.将式(28)代入式(27)整理后得到其中，式(27)成立的充要条件为式(30)的展开式一般为κ的四阶方程，但因四阶项的系数和常数项为零而退化为二阶方程.求解此二阶方程得到κ 的 2 个根κ1、κ2，再将κ =κ1，κ2代入式(29)得到对应的向量φ=φ1，φ2，其中，由于内域稳定准则为0＜Δτ≤1，假定在此区间内 |z|＞ 1，不难证明|κ1|＜ 1，|κ2|＜ 1.注意到模态式(26)在界面节点两侧的表达式为式中:b1、b2为常数.根据界面位移连续条件得到b1φ1(z)=b2φ2(z)，上式可写成因此，若格式(24)支持失稳模态式(26)，则将发生如式(32)所示失稳(图2).图2 一维界面模型的失稳形式Fig.2 The form of G-R instability of the one-D model with an interface关于式(24)支持上述失稳模态的论证是在|z|＞1的假定下做出的.但是，若采用前述经验稳定准则，这一假定不成立.将式(32)代入递推式(24)，若b1不为零得到当p=3/2时从式(33)给出的2个相容的方程得到解z=-1.这与假定|z|＞1矛盾，即式(24)不支持|z|＞1的模态式(32).设有限弹性杆中界面(x=0)右侧和左侧的波速分别为c1=100 m/s和c2=50 m/s，左端(x=-300 m)固定，右端(x=600 m)自由，界面阻抗比α=2.给定初始速度为零，初始位移f(x)为3次B样条函数，其非零部分分布于区间［-25 m，0］.界面节点和内域节点的递推公式由式(24)给出.两端节点的递推公式可作为界面节点递推公式的极限情形得到，例如，若j=0的右侧为自由面，令j=0和α=0，则式(24)退化为自由端递推公式.若j=0的右侧为固定端，令j=0和α→+∞，式(24)则退化为固定端递推公式.为检验上述稳定准则，在数值实验中仅变动界面节点和两端点的p值，而对内域节点取用固定值p=3/2.实验结果表明，当p≠3/2发生失稳.图3给出界面节点的一组数值解(Δτ=1，界面左侧和右侧空间步距分别为3 m和6 m). 图3 稳定方案与失稳方案的界面节点位移时程Fig.3 Displacement time history of interface node in stability plan and instability plan4 弱稳定性分析弱稳定性分析可用于线法(method of lines)建立的时空离散格式，亦可用于其他离散方法建立的数值解格式，只要当k→0时该格式退化为一组常微分方程，且在h→0时后者收敛于原偏微分方程问题的解.应用弱稳定性分析的主要问题是检查常微分方程组式(5)初值问题的解是否收敛于连续模型的解.下面针对有限差分和有限元方法分别讨论这一问题.4.1 有限差分解法就弱稳定性分析在有限差分解法中的应用而言，Strikwerda指出:若空间离散格式至少具有一阶精度，则上述收敛性等价于如下适定条件［11］.式中:C为与t无关的常数.假定式(5)中B为常系数方阵且维数不变，则适定条件式(34)可用B的谱判据检验［19］.但是，式(5)的收敛性涉及维数变化的一簇矩阵B，文献［19］所提供的判据不再有效.除对十分简单的模型可以证明式(34)成立，例如，Strikwerda通过类似Godunov和Ryabenkii采用的用无限模型代替有限模型的方法，对式(14)的一个简单的空间离散格式证明了式(34)成立［11］.在一般情形下，式(34)是否成立则难以论证.因此，弱稳定性分析一般不能用于有限差分格式的稳定性分析.它的另一缺陷来自回避了实际的时空离散模型，从而排除了分析在实际数值计算过程中出现的某些现象的可能，例如，阐明局部失稳机理的可能.不过，在弱稳定性分析的前提成立的条件下，弱稳定性能消除所有可能的失稳形态，从而达到保证数值解格式稳定的目的.这是弱稳定性分析的价值所在，下面通过分析有限元解法阐明这一价值.4.2 有限元解法相对于具有局域近似特征的有限差分方法，基于泛函变分原理的有限元方法具有全局近似特征.这一特征使得用有限元方法建立的常微分方程组式(5)的收敛性易于得到保证.以线性粘弹性波的数值模拟为例，设边界条件为给定位移、速度或应力，或给定应力与位移和速度的线性关系，则用有限元法得到的式(5)可写成二阶形式: Mw..(t)+Cw·(t)+Kw(t)=0. (35)式中:M、C、K分别为质量、阻尼、刚度矩阵，其中，M对称正定，K对称半正定，C正定.Dupout证明了式(35)的解收敛于原连续模型的解［20］.除粘弹性波动外，浅水波等问题有限元离散所得常微分方程组的收敛性亦已获得证明(参看文献［20］引用的文献).对于若干具有实际意义的情形，式(35)时间离散格式的稳定性已得到证明.首先，若C为Rayleigh阻尼，则式(35)等价于解耦的一组单自由度方程，可采用常规谱判据判断其稳定性［5］.其次，当C为一般正定阵时，用P层线性多步法离散式(35)，得式中:ακ、βκ、γκ为实数，αp=1，Z 为时间平移算子，Zκvn=vn+κ.Gekeler将由 M、C、K 和ρ(z)、σ(z)、τ(z)规定的离散格式的稳定问题归结为［21］:是否存在实参数η∈-s，[]0，s＞0，使方程ρ(z)-ησ(z)=0的所有根z满足|z|≤1且不存在|z|=1的复根?肯定回答是格式稳定的一个必要条件.在此条件下，稳定准则为式中:λi为对称阵M-1/2KM-1/2的特征值.5 结束语本文将Lax稳定性分析区分为强和弱2种稳定性分析方法.它们对进一步发展波动数值模拟技术各有其用途.强稳定性概念可用于任何时空离散模型的稳定性分析，特别适用于研究局部稳定性.它不仅可用于阐明实际模型中发生的典型局部失稳模态，而且对局部失稳模态取决于少数因素之匹配的认识，为探寻稳定的模拟方案提供了切实的技术路线.合理地使用弱稳定性分析可排除包含局部失稳在内的所有失稳模态，这一分析方法特别适用于完善波动的有限元模拟技术，例如，建立和论证稳定的、与内域数值解格式精度阶匹配的实用人工边界.参考文献:【相关文献】［1］LIAO Z P ，LIU J B.Numerical instabilities of a local transmitting boundary［J］.Earthquake Engineering ＆Structural Dynamics，1992，21(1):65-77.［2］LIAO Z P，ZHOU Z H ，ZHANG Y H.Stable implementation of transmitting boundary in numerical simulation of wave motion［J］.Chinese Journal of Geophysics，2002，45(4):554-568.［3］CRANK J，NICOLSON P.A practical method for numerical evaluation of solutions of partial differential equations of heat conduction type［J］.Proc Cambridge Philos Soc，1947，43:50-67.［4］HOFFMAN J D.Numerical methods for engineers and Scientists［M］.NewYork:McGraw-Hill Inc，1992:657-691.［5］HUGHES J R.Analysis of transient algorithms with particular reference to stability behavior［C］//Computational Methods for Transient Analysis.North-Holland:Elsevier Science，1983:67-155.［6］GODUNOV S，RYABENKII V.Spectral stability criteria for boundary value problems for non-self-adjoint difference equations［J］.Russian Mathematical Surveys，1963，18(3):1-12.［7］KREISS H O.Stability theory for difference approximations of mixed initial boundary value problems［J］.Mathematics of Computation，1968，22(104):703-714.［8］GUSTAFSSON B，KREISS H O，SUNDSTRÖM A.Stability theory of difference approximations for mixed initial boundary value problems［J］.Mathematics of Computation，1972，26(119):649-686.［9］GUSTAFSSON B，KREISS H O，OLIGER J.Time dependent problems and difference methods［M］.New York:John Wiley＆ Sons，1995:496-620.［10］TREFETHEN L.Group velocity interpretation of the stability theory of Gustafsson，Kreiss and Sundstr［J］.Journal of Computational Physics，1983，49(2):199-217. ［11］STRIKWERDA J C.Finite difference schemes and partial differential equations ［M］.2nd ed.［S.1］:SIAM，2004:1-34.［12］RICHTMYER R D，MORTON K W.Difference methods for initial-value problems ［M］.2nd ed.New York:Interscience，1967:39-58.［13］STRIKWERDA J C，WADE B A.A survey of the Kreiss matrix theorem for power bounded families of matrices and its extensions，in linear operators［C］//Linear Operators.Warszaw:Inst Math Pol Acad Sciences，1997:329-360.［14］PARTER S V.Stability，convergence，and pseudo-stability of finite-difference equations for an over-determined problem［J］.Numerische Mathematik，1962，4(1):277-292.［15］关慧敏，廖振鹏.局部人工边界稳定性的一种分析方法［J］.力学学报，1996，28(3):376-380.GUAN H M，LIAO Z P.A method for the stablitity analysis of local artificial boundaries ［J］.Acta Mechanica Sinica，1996，28(3):376-380.［16］SOUSA E.On the edge of stability analysis［J］.Applied Numerical Mathematics，2009，59(6)，1322-1336.［17］谢志南，廖振鹏.波动方程数值模拟的一种显式方法-界面节点递推公式的构建［J］.力学学报，2011，43(1):154-161.XIE Z N，LIAO Z P.An explicit method for numerical simulation of wave motion-constructing recursion formula for interface point［J］.Acta Mechanica Sinica，2011，43(1):154-161.［18］廖振鹏，刘恒，谢志南.波动数值模拟的一种显式方法——维波动［J］.力学学报，2009，41(3):350-360.LIAO Z P，LIU H，XIE Z N.An explicit method for numerical simulation of wave motion:1-D wave motion［J］.Acta Mechanica Sinica，2009，41(3):350-360.［19］HAIRER E，NØRSETT S P，WANNER G.Solving ordinary differential equationsI:nonstiff problems［M］.2nd ed.Berlin:Springer，1993:80-89.［20］DUPONT T.L2-estimates for Galerkin methods for second order hyperbolic equations［J］.SIAM Journal on Numerical Analysis，1972，10(5):880-889.［21］GEKELER E.Linear multistep methods for stable differential equations［J］.Mathematics of Computation，1982，39(60):481-490.。
波动方程求解方法
常用的波动方程求解方法主要有以下几种:有限差分法、有限元法和伪谱法、积分方程法等。
1、有限差分方法由于适应性强，计算快速，因此是最先发展起来而且使用范围最广的数值方法，有限差分方法最大的弱点之一就是会产生数值频散。
有限差分法采用差分算式近似逼近偏导数运算，从而使波动方程的偏导数运算问题转化成差分代数问题，最后通过求解差分代数方程组得到近似解结果。
有限差分法的差分算式本身就是一种局部点运算，不需要考虑原函数中所求点值在邻域范围上的函数的变化情况，而只需要用到所求点值附近点上的值，所以能够很好的适用于复杂情况, 但是难保模拟精度。
有限差分方法有较高的空间域分辨率，而在频率域上分辨率反而会极低，稳定性同时还受到网格间距和时间步长的影响。
同时，虽然有限差分法还伴随有数值频散的问题，但是计算速度较快。
有限差分法目前主要有以下三大类:规则网格方程、弹性方程和交错网格方程。
有限差分法的具体操作可以分为两个部分:(1)用差分代替微分方程中的微分，将连续变化的变量离散化，从而得到差分方程组的数学形式:(2)求解差分方程组。
在第一步中，通过网格剖分法，将函数定义域分成大量相邻而不重合的子区域。
通常采用的是规则的剖分方式，最常用的是正方形网格。
这样可以便于计算机自动实现和减少计算的复杂性。
网格线划分的交点称为节点。
若与某个节点P 相邻的节点都是定义在场域内的节点，则P 点称为正则节点;反之，若节点P 有处在定义域外的相邻节点，则P 点称为非正则节点。
在第二步中，数值求解的关键就是要应用适当的计算方法，求得特定问题在 所有这些节点上的离散近似值。
目前最常用的两种有限差分方法包括:基于位移 波动方程的二阶中心差分法和基于一阶速度-应力波动方程的高阶交错网格法， 前者算法简单，易于实现，但差分精度具有局限性，最后得到的是节点上z x ,分量的位移离散近似值，后者算法稍复杂，但可以提高差分精度，最终得到的是节点上的位移速度离散近似值。
内域波动数值模拟的显式方法的开题报告
内域波动数值模拟的显式方法的开题报告1. 研究背景和意义内域波动是指在流体内部产生的小尺度涡旋和涡动，其在工程中有重要的应用，如湍流控制、燃烧、混合、传热等领域。
因此，精确模拟内域波动对工程应用有着重要的意义。
液体或气体的流动都是由基本物理方程所描述，其中包括质量守恒、动量守恒、能量守恒等方程，在大篇幅条件下方程的数值求解是十分困难的，因此需要运用一些适当的计算方法来解决这些方程。
2. 研究本体和内容本文将研究内域波动的数值模拟方法，并重点分析显式方法的数值模拟。
首先，将讨论数值方法的基本原理和数值算法，然后介绍数值模拟的实现过程和算法，包括在不同边界条件下的求解方式。
在此基础上，将对内域波动的数值模拟与实验结果进行对比和分析，以验证模拟方法的有效性。
最后，对模拟结果进行总结，并对模拟方法的优化和改进提出一定的建议。
3. 研究方法研究方法主要包括理论分析、数值模拟和实验对比。
首先通过理论分析，确定数值模拟的基本原理和算法。
接着采用数值模拟的方法对内域波动进行模拟，并对不同的边界条件进行讨论和对比。
最后将模拟结果进行与实验对比，以验证模拟方法的有效性。
4. 研究进度计划第1-2周：研究文献，了解内域波动数值模拟的基本方法和原理第3-5周：确定数值模拟的基本方法和算法，编写程序进行验证第6-7周：进行数值模拟，对不同的边界条件进行讨论第8-9周：将模拟结果进行与实验对比，验证模拟方法的有效性第10-11周：对模拟结果进行总结，提出模拟方法优化和改进的建议第12-14周：完成毕业论文的写作和修改，准备答辩5. 研究预期成果通过本次研究，预计可以得到以下成果：1）研究内域波动的数值模拟方法和显式方法的数值模拟；2）对不同边界条件下内域波动的数值模拟进行讨论和对比；3）验证模拟方法的有效性，对模拟结果进行总结并提出建议；4）编写毕业论文并顺利完成答辩。
6. 研究难点和解决思路研究内域波动数值模拟的关键是数值方法的选择和算法实现。