波动方程的变步长有限差分数值模拟

合集下载

5第五讲 典型模型方程-波动方程的差分格式

5第五讲 典型模型方程-波动方程的差分格式

修正方程:
u u c x 2u c x2 3u 1 2 2 1( 1) 3 c x x t x 2 6
a t x
稳定条件 :

1
耗散及频散特性. G 1 c c cos ic sin
Euler差分方法(如下)不稳定
为了使其为稳定差分方法,修正如下

因此,此差分方程不一定总相容于PDE,当 则相容,为条件相容。
t 0, x 0 x
t 0, x 0
2
t
0
?
时,v保持为常数,
Lax差分方法的耗散性和色散性分析
当v不为1时,耗散性很强。
理解下图???
Euler差分方法(如下)不稳定
初值问题的精确解
由数学物理方程得:
(3) (4)
由(4)分别求t、x导数,代入(2)式可以验证解的正确性, 即: cF cF 0
欧拉显格式:Euler Explicit Method
都为一步显格式
1、 2、
: :
1 un j
时间、空间前差
时间前差、空间中心差
两种差分方程的时间项都为一阶精度。 每个方程中未知量只有 ,都为一步显格式 显格式: 递推;而不同求解方程组
蛙跳格式的耗散性和色散性分析
蛙跳格式的缺点
1、初始值需要给两个时间层上的值。解决:可以 在初始计算时,仅给n=0时刻的初始值,采用两层 差分格式来得出n=1层上的值,然后在采用蛙跳格 式循环计算。
n u 不依赖于 j ,导致了 2 n 1 3 1 2n 4 2 两个不相关的解。 u j ...u j u j ,u j ...u j u j
一阶迎风格式的放大因子:

地震波波动方程数值模拟方法(严选优质)

地震波波动方程数值模拟方法(严选优质)

地震波波动方程数值模拟方法地震波波动方程数值模拟方法主要包括克希霍夫积分法、傅里叶变换法、有限元法和有限差分法等。

克希霍夫积分法引入射线追踪过程,本质上是波动方程积分解的一个数值计算,在某种程度上相当于绕射叠加。

该方法计算速度较快,但由于射线追踪中存在着诸如焦散、多重路径等问题,故其一般只能适合于较简单的模型,难以模拟复杂地层的波场信息。

傅里叶变换法是利用空间的全部信息对波场函数进行三角函数插值,能更加精确地模拟地震波的传播规律,同时,利用快速傅里叶变换(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技术将系数矩阵转化为三个一维数组来存储,分别存储系数非零元素、非零元素对应所在的行以及每列起始非零元素所在位置.经CSC技术压缩存储后显著减少了内存空间及计算量,在计算时只有少许的非零元素参加计算,且根据三个一维数组可以简便地找到对应的非零元素,进而采用LU分解快速而精确地求解.本文基于Jo等提出的最优化9点差分方法,首次应用CSC技术在频率空间域进行二维声波方程数值模拟.通过对Corner-edge模型和二维Marmousi模型进行试算,可以显著节省内存需求,明显提高计算速度,进而得到精度较高的正演结果.【期刊名称】《石油地球物理勘探》【年(卷),期】2014(049)002【总页数】7页(P288-294)【关键词】频率—空间域;CSC;系数矩阵;一维数组;LU分解;数值模拟【作者】吕晓春;顾汉明;成景旺;周丽【作者单位】中国地质大学(武汉)地球物理与空间信息学院,湖北武汉430074;中国地质大学(武汉)构造与油气资源教育部重点实验室,湖北武汉430074;中国地质大学(武汉)地球物理与空间信息学院,湖北武汉430074;中国地质大学(武汉)构造与油气资源教育部重点实验室,湖北武汉430074;中国地质大学(武汉)地球物理与空间信息学院,湖北武汉430074;中国地质大学(武汉)构造与油气资源教育部重点实验室,湖北武汉430074;中国地质大学(武汉)地球物理与空间信息学院,湖北武汉430074【正文语种】中文【中图分类】P6311 引言为了反演地下介质参数或研究地震波在各种复杂介质中的传播机制,需要进行波场数值模拟。

波动方程的数值模拟方法包括有限差分法、有限元法、伪谱法等。

在时间—空间域的数值模拟技术已较成熟,并广泛应用于复杂介质的正演模拟中。

一维波动方程的有限差分法

一维波动方程的有限差分法

学生实验报告实验课程名称偏微分方程数值解开课实验室数统学院学院数统年级2013 专业班信计02班学生姓名学号开课时间2015 至2016 学年第 2 学期数学与统计学院制开课学院、实验室:数统学院实验时间:2016年6月20日五.实验结果及实例分析1、u(x,t)在t=0.5,1.0,1.5,2.0时刻的数值解、精确解以及绝对误差表1 u(x,t)在t=0.5,1.0,1.5,2.0时刻的数值解时刻tt=0.5,1.0,1.5,2.0时刻的数值解t=0.5 0 -0.0059 -0.0113 -0.0155 -0.0182 -0.0192 -0.0182 -0.0155 -0.0113 -0.0059 0 t=1.0 0 -0.3090 -0.5877 -0.8090 -0.9510 -0.9999 -0.9510 -0.8090 -0.5877 -0.3090 0 t=1.5 0 0.0020 0.0038 0.0052 0.0061 0.0064 0.0061 0.0052 0.0038 0.0020 0 t=2.0 0 0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090 0.5878 0.3090 0表2 u(x,t)在t=0.5,1.0,1.5,2.0时刻的精确解时刻tt=0.5,1.0,1.5,2.0时刻的精确解t=0.5 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0 t=1.0 0 -0.3090 -0.5878 -0.8090 -0.9511 -1.0000 -0.9511 -0.8090 -0.5878 -0.3090 0 t=1.5 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0 t=2.0 0 0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090 0.5878 0.3090 0表3 u(x,t)在t=0.5,1.0,1.5,2.0时刻的绝对误差时刻tt=0.5,1.0,1.5,2.0时刻的绝对误差t=0.5 0 0.0059 0.0113 0.0155 0.0182 0.0192 0.0182 0.0155 0.0113 0.0059 0 t=1.0 0 0.0000 0.0000 0.0001 0.0001 0.0001 0.0001 0.0001 0.0000 0.0000 0 t=1.5 0 0.0020 0.0038 0.0052 0.0061 0.0064 0.0061 0.0052 0.0038 0.0020 0 t=2.0 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0说明:在t=0.5时刻的绝对误差最大,t=1.5时刻次之,t=1与t=2时刻的绝对误差均较小,由于0.11r hτ==<,该格式稳定,由数值计算得到的矩阵不难看出,数值解符合理论解。

有限差分法地震波传播数值模拟

有限差分法地震波传播数值模拟

=
kΔx
=
2π λ
Δx ≤ 1
,即只要一个波长包含几个空间步
长,随着差分精度2M的提高,上述高阶差分解法产生
的数值频散会逐渐减小。
不同差分精度空间频散曲线
不同差分精度时间频散曲线
五、边界问题
自由边界条件
内部边界条件 吸收边界条件
设计吸收边界条件的目标:
z 方程+边界条件数学上是非病态的
连续
问题 z 方程+边界条件可以近似描述无限介质中的物理过程 z 边界条件和内部点的计算方式是相容、不冲突的
-----------J.M. Carcione
地震波传播数值模拟应用领域
地震波传播理论
数据采集
理论指导 物性参数
研究传播规律
正演模拟
指导设计 观测系统
验证
地震解释
提供理论数据 试验处理流程
数据处理
提供正演方法
岩石物理
参数反演
断层下覆界面反射能量强
炮点
T=2000ms
炮点 T=2300ms
炮点位于11km处的单炮记录
?21?4?1?6?1?m?2m?1222m246333m24lllol622m32m?m4?m6??m?m?2m?m?2?c1m??m??c2?m??c3???m??cm??m??1??0????0???m????0??35?1?133335?555?135?mm?m2n?12n?12n?1?35?1llloln?2n?1?c1??1??n???3?2n?1??c2??0?n?5??0?2n?1?c??3???m??m??m?n?2n?1????2n?10???cn????ox2ox数值频散试验dxdz10mdt1ms10高阶差分为何会消除数值频散

基于波动方程有限差分数值模拟及FCT消频散分析

基于波动方程有限差分数值模拟及FCT消频散分析
第 36卷第 2期
中 国 锰 业
2018年 4月
CHINA′SMANGANESEINDUSTRY
Vol.36No.2 Apr.2018
基于波动方程有限差分数值模拟及 FCT消频散分析
赵 威,李伟波,桂志先,周 游,于晓东
(长江大学 油气资源与勘探技术教育部重点实验室,湖北 武汉 430100)
Δx=Δz=h
(4)
1.2 一阶速度—应力弹性波方程
在二维各向同性介质中,不考虑体力的影响下,
弹性波方程可以化为一阶速度—应力的形式:



!"
$ ""
% "&

#
"
&




!& '
$ "& "
% && &


""

'
$()*
!" "
高,占用内存小等特点,但由于其本身时间和空间都 采用差分 离 散 形 式,往 往 会 出 现 数 值 频 散 问 题[2]。 怎样减少频散,提高模拟精度始终是有限差分数值 模拟的关键问题。为此,学者们做了大量工作,当信 号在一个波长范围内,具有多个离散节点是,频散现 象基本消失。一般随着网格间距增大,频散问题越 严重;反之若减小网格间距相应会增加网格点数,将 会大大增加计算量,造成计算时间过长。既能控制 模拟精度,又能节省计算时间的方法显得尤为重要。 Boris和 Book等人 首 先 提 出 了 通 量 校 正 传 输 的 方 法,它最早起源于求解流体力学连续方程中,随后才 开始被运用于声波方程数值模拟的求解;国内杨顶 辉等[3]将其 应 用 于 弹 性 波 动 方 程 正 演 数 值 模 拟 频 散处理,取得一定效果。为了更好的对比分析未加 FCT数值模拟和加了 FCT方法的区别,本文基于声 波和弹性波两方面,借助于交错网格的模拟优势,对 各向同性介质波场进行了有限差分数值模拟,并对 其相应的波场快照和计算效率进行对比分析。

基于GID有限元前处理的波动方程数值模拟

基于GID有限元前处理的波动方程数值模拟

基于GID有限元前处理的波动方程数值模拟刘静;文山师;黄晶晶【摘要】在地震波数值模拟计算过程中,缺乏简单易行的有限元前处理方法,使得复杂构造模型较难建立和分析.本文以二维声波方程为例结合GID软件,网格剖分部分采用三角形单元模拟速度界面,把单元内的场和波速均看作单元上的线性函数;GID 软件可以方便地进行网格剖分和设置网格控制节点,通过编写用户自定义”问题类型”,建立并输出已有的有限元计算程序的初始模型.将GID软件前处理与有限元计算程序整合,提高了方法的效率,简单易行.【期刊名称】《工程地球物理学报》【年(卷),期】2014(011)002【总页数】7页(P243-249)【关键词】数值模拟;有限元;GID;声波方程;三角形单元【作者】刘静;文山师;黄晶晶【作者单位】山西省煤炭地质115勘查院,山西大同037003;中石化西北油田分公司勘探开发研究院,新疆乌鲁木齐830011;中石化石油工程地球物理有限公司河南分公司,河南南阳473000【正文语种】中文【中图分类】P631.41 引言地震波场的数值模拟技术是在已知地下介质结构和参数的情况下,利用理论计算的方法研究地震波在地下介质中的传播规律,合成地震记录的一种技术。

地震勘探中的数值模拟方法主要以射线理论和波动方程理论为基础,有射线追踪法,柯西霍夫积分法,有限元法,有限差分法和伪谱法[1~6]。

有限差分法直接用差分代替微分,因其方法简单、精度高,在地震模拟中而得到了广泛的研究和应用。

但其固有缺陷是不能准确模拟具有复杂几何形态的物性界面,有限元法则是求解原问题等价泛函的变分或原问题的等效积分方程的弱解(当等价泛函不存在时),因而能够适应较有限差分更为剧烈的物性变化,加之种类繁多的插值形函数,使其能够模拟很复杂的几何界面。

有限元法的主要缺点是计算和存储量都很大,效率相对较低。

建立有限元分析模型比较复杂且存在困难,因此可以用一些成形的软件作为有限元网格剖分的工具,建立并输出可用于已有有限元计算程序的初始模型,将大大提高方法的效率[7]。

声波方程有限差分数值模拟的变网格步长算法

声波方程有限差分数值模拟的变网格步长算法

声波方程有限差分数值模拟的变网格步长算法声波方程有限差分数值模拟是一种常用的声波传播模拟方法,可以在计算机上通过数值计算求解声波传播的过程。

在进行这种数值模拟时,常常需要选择合适的网格步长,以保证计算结果的准确性和计算效率。

本文将介绍一种变网格步长算法,用于优化声波方程有限差分数值模拟的计算。

声波方程可以用下面的形式表示:∂^2p/∂t^2=c^2∇^2p其中p是声场变量,t是时间,c是声速,∇^2是Laplace算子。

为了将声波方程用有限差分方法进行离散化计算,我们需要将空间和时间分别离散化。

首先,将空间离散化为网格,在每个网格点上计算声场的值。

其次,将时间离散化为离散的时间步长,通过迭代计算不同时间步长上的声场分布。

为了保证计算结果的准确性,网格步长应当满足Nyquist采样定理的要求。

即网格步长应小于声波的最小波长的一半。

根据声波方程的性质,我们可以通过声速和最高频率来估计声波的最小波长。

然后,我们可以根据最小波长来选择合适的网格步长。

然而,在实际的声波传播计算中,声场的变化往往不是均匀的。

有些区域的声场变化较大,而其他区域的声场变化较小。

如果我们在整个计算区域都采用较小的网格步长,将会造成计算资源的浪费。

因此,需要一种方法能够根据声场的变化情况来自适应地调整网格步长。

变网格步长算法就是一种能够根据声场变化情况自动调整网格步长的算法。

其基本思想是根据声场在不同网格上的变化率来决定每个网格上的网格步长。

具体的算法步骤如下:1.初始化:选择一个合适的初始网格步长。

通常可以选择根据声波的最小波长来确定。

2.计算网格步长:在每个时间步长上,对于每个网格点,计算其周围网格点上的声场变化率。

常用的方法是计算声场在三个相邻时间步长上的差分值,然后取绝对值并求平均。

根据声场变化率,调整当前网格点上的网格步长。

变化率大的网格点应该有更小的网格步长,而变化率小的网格点则可以有更大的网格步长。

3.更新声场:根据调整后的网格步长,更新所有网格点上的声场值。

三维非均匀介质波动方程有限差分python开源代码

三维非均匀介质波动方程有限差分python开源代码

文章标题:探索三维非均匀介质波动方程有限差分python开源代码1. 简介在地质勘探、医学成像和地震监测等领域,对三维非均匀介质波动方程的研究与应用日益重要。

而有限差分方法在数值求解波动方程中具有广泛的应用。

在本文中,我们将探讨如何利用Python编程语言实现三维非均匀介质波动方程的有限差分方法,并开源共享相应的代码,以便更多人能够深入理解和应用这一重要领域。

2. 三维非均匀介质波动方程简介三维非均匀介质波动方程描述了波在非均匀介质中的传播规律,是地震勘探、医学成像等领域中常见的数学模型之一。

该方程的数值求解通常采用有限差分方法,通过离散网格化空间和时间来逼近连续的微分方程,从而得到数值解。

3. 有限差分方法有限差分方法是数值求解微分方程的一种常见方法,其基本思想是将微分方程中的导数用差分近似代替,从而将连续的问题转化为离散的问题。

在三维非均匀介质波动方程中,有限差分方法可以有效地模拟波的传播过程,并得到波场的数值解。

4. Python编程实现利用Python编程语言实现三维非均匀介质波动方程的有限差分方法具有许多优势,如简洁易读的代码、丰富的科学计算库等。

在实现过程中,我们可以利用NumPy库进行数组操作,使用Matplotlib库进行波场可视化,并通过SciPy库进行数值求解等。

5. 开源代码共享在本文中,我们将共享我们编写的三维非均匀介质波动方程有限差分Python开源代码,包括空间离散化、时间离散化、边界条件处理、波场更新等关键部分。

我们也会附上详细的注释和使用说明,以便感兴趣的读者能够下载并运行我们的代码,深入理解和学习有限差分方法在波动方程中的应用。

6. 个人观点和理解通过编写三维非均匀介质波动方程的有限差分Python开源代码,我深刻体会到数值模拟在地质勘探、医学成像等领域中的重要作用。

Python作为一种强大的科学计算语言,为我们提供了丰富的工具和库,使得数值模拟变得更加高效和灵活。

波传播和地震模型的部分求和有限差分法

波传播和地震模型的部分求和有限差分法

波传播与地震模型分析: 部分求和有限差分法的应用地球物理场景中,波的传播和地震模型的研究是基于对地质结构及其物理性质的精确描述。

部分求和有限差分法是一种强有力的数值分析方法,用于解决这些模型在变化复杂的介质中传播的问题。

本文旨在从该方法的理论基础、实现方案到最终的应用效果三个方面进行全面阐述。

一、波动方程与有限差分法的原理波的传播现象可以通过波动方程来描述,而有限差分法则提供了一种强大的数值分析工具,对这一方程进行离散化处理。

(1)波动方程概览波动方程作为一种偏微分方程,能够刻画波在介质中传播的速度、方向和形态。

地震波动方程则特别体现了地震波在地球介质中的传播规律。

(2)有限差分法导论有限差分法通过将连续域离散化,使用格点近似连续函数及其导数,将微分方程转化为代数方程组。

从而可以利用计算机进行数值求解。

二、部分求和有限差分法的理论与实施在传统有限差分法的基础上,部分求和技术的引入,有效提升了数值计算的稳定性和精确度。

(1)部分求和技术的理念部分求和技术(PS-FDM)在于对邻近的离散点进行部分求和处理,减少高频数值震荡,这对于模拟高频率的地震波尤其有效。

(2)实施方案实施部分求和有限差分法时,构建精确的离散网格是首要步骤。

之后,确定合适的时间步长和空间步长,依靠合适的算法框架,将求和处理的策略嵌入其中。

三、地震模型下的波传播分析运用部分求和有限差分法对实际地震模型进行分析,可以模拟出地震波在不同地质结构中的传播效果。

(1)模型设定根据地质数据建立地震模型,包括地层的分布、物理属性以及断层等结构特征。

(2)波场模拟通过部分求和有限差分法模拟不同类型地震波的传播过程,提取波场数据。

(3)结果分析与应用最终得到的波场数据可以用于评估地震对建筑结构的影响,或用于地震预警系统。

通过对波场模拟结果的详细分析,可以优化地震模型,提高地震预测的准确性。

延续前文的论述,对部分求和有限差分法(PS-FDM)的理论框架进行了详细展开,本节将着重讲述该方法的具体实施步骤,以及如何针对地震模型进行高效的波传播模拟。

声波方程有限差分数值模拟的变网格步长算法

声波方程有限差分数值模拟的变网格步长算法
第4卷 第3期 2007 年 6 月
工程地球物理学报
CHIN ESE JO U RN A L O F EN GI NEERIN G G EOP HY SICS
文章编号: 1672 7940( 2007) 03 0207 06
V ol 4, N o 3 Jun , 2007
声波方程有限差分数值模拟的 变网格步长算法
缺乏稳定性。本文对解决上述问题有较高优越性的变网格算法进行了介绍, 对传统的有限差分法与变网格差
分算法在内存需求、计算速率等方面的差别进行 了比较, 并对变网格算法中的边界条件、时间积分的快速展开
算法进行了阐述, 总结了变网格算法的优点。
关键词: 正演模拟; 变网格; 边界条件; 网格步长
中图分类号: P631
算[ 5] 。为简便起见, 仅对基本方程( 1) 的快速展开
法做一些介绍, 但对变密度的情况可用类似的方
法进行推导。
对震源项 s = ( x , z , t ) = g( x , z ) h( t ) 边界条
件为吸收边界时, 式( 1) 的常规解为
P(x , z, t) =
[
sin( L) h( t 0L
( 7)
这里 2N 是算 子长度[ 7] 。如 果在深 度 z 0 =
l z 的步长是双倍的, 式( 6) 在长度为 2N 的有限
差分法可被使用到深度为( l - N ) z , 式( 7) 可在
深度 l z 以下使用。在剩余的( N - 1) 个网格点,
可结合式( 6) 和式( 7) 来计算其导数。在深度( l -
文献标识码: A
收稿日期: 2007 03 19
Acoustic Equation Numerical Modeling on a Grid of Varying Spacing

变步长激发极化法有限元数值模拟

变步长激发极化法有限元数值模拟

其 中 u = ( 1 8 u …u )

警 等 c
(i) d ̄ c Ny
警d叼 d

( 3
d警( ・ N) i警
3 3 2 边 界 积 分 ..
如果 单元 的 一个 面 口 落 在 无 穷 远 边 界 上 ,

则边 界积 分 为
图 1 模 型剖 分 示 意 图
维普资讯
第3 卷 第3 0 期
物探 化探 计 算 技 术
28 月 0 年5 0
文章 编号 :10 一 l4 2 0 ) 3 - 2 2 4 0 1 7 9(0 8 O - 0 l —o
变 步 长 激发 极 化 法有 限元数 值 模 拟
杨晓弘, 何继善, 童孝忠
率 ; 为 电 流强 度 ; 为 模 型 区域 ;凡为边 界 的外 , 力 法 向方 向 ;
测 点到 电源点 A、 B的距 离 ; 为无 穷远边 界 。 F
在剖分单元上进行插值 , 合理存储总刚度矩阵 , 用 常规 L L D T方法计算大型稀疏矩阵 , 等效视电阻 用 率 法来计 算模 型视极 化率 , 在缩 短计算 时 间的情况
下 , 证 了计算 结果 的精度 。 保
2 变 分 问题
与 边 值 I 题 式 ( ) 价 的 变 分 I 为 口 l ] 1等 司题
1 三维模 型的边值问题
中间梯 度 法是 激 发极 化 法 中最 常用 到 的物 探 方法 之一 , 特别适 合 寻找直 立 的高阻脉 和水平 的 良
导层 , 在三 维构 造 中 , 点 电源 电场 模 型 的 电位 双
中 图分 类号 :P6 13 2 3 . 4 文献 标识码 :A 边值 问题 为

有限差分法及其应用

有限差分法及其应用

有限差分法及其应用1有限差分法简介有限差分法(FDM)是计算机数值模拟最早采用的方法,至今仍被广泛运用。

该方程将解域划分为差分网格,用有限个网络节点代替连续的求解域。

有限差分法通过泰勒级数展开等方法,把控制方程中的导数用网格节点上的函数值得差商代替进行离散,从而建立以网格节点上的值为未知数的代数方程组。

该方法是一种直接将微分问题变为代数问题的近似值解法,数学概念直观,表达简单,是发展较早且比较成熟的数值方法。

2有限差分法的数学基础有限差分法的数学基础是用差分代替微分,用差商代替微商而用差商代替微商的意义是用函数在某区域内的平均变化率来代替函数的真是变化率。

而根据泰勒级数展开可以看出,用差商代替微商必然会带来阶段误差,相应的用差分方程代替微分方程也会带来误差,因此,在应用有限差分法进行计算的时候,必须注意差分方程的形式,建立方法及由此产生的误差。

3有限差分解题基本步骤有限差分法的主要解题步骤如下:1)建立微分方程根据问题的性质选择计算区域,建立微分方程式,写出初始条件和边界条件。

2)构建差分格式首先对求解域进行离散化,确定计算节点,选择网格布局,差分形式和步长;然后以有限差分代替无线微分,以差商代替微商,以差分方程代替微分方程及边界条件。

3)求解差分方程差分方程通常是一组数量较多的线性代数方程,其求解方法主要包括两种:精确法和近似法。

其中精确法又称直接发,主要包括矩阵法,高斯消元法及主元素消元法等;近似法又称间接法,以迭代法为主,主要包括直接迭代法,间接迭代法以及超松弛迭代法。

4)精度分析和检验对所得到的数值进行精度与收敛性分析和检验。

4商用有限差分软件简介商用有限差分软件主要包括FLAC、UDEC/3DEC和PFC程序,其中,FLAC是一个基于显式有限差分法的连续介质程序,主要用来进行土质、岩石和其他材料的三维结构受力特性模拟和塑性流动分析;UDEC/3DEC是针对岩体不连续问题开发,用于模拟非连续介质在静,动态载荷作用下的反应;PFC是利用显式差分算法和离散元理论开发的微、细观力学程序,它是从介质的基本粒子结构的角度考虑介质的基本力学特性,并认为给定介质在不同应力条件下的基本特征主要取决于粒子之间接粗状态的变化,适用于研究粒状集合体的破裂和破裂发展问题,以及颗粒的流动(大位移)问题。

有限差分法不同边界条件下的数值模拟

有限差分法不同边界条件下的数值模拟

有限差分法不同边界条件下的数值模拟文章介绍了地震数据处理中所使用的数值模拟法,对采用有限差分法所使用不同边界条件处理方式进行了数值模拟,通过波场快照直观的得出了不同的边界吸收条件的吸收效果,对结果进行了对比,分析总结了各种方法的优缺点。

标签:数值模拟;有限差分;边界条件随着近年来国家宏观经济调控,经济增长的速度逐步减缓,能源行业受此影响最为严重,许多煤矿是在亏损的情况下生产,直接导致了地质行业投入的减少。

物探行业压力也越来越大,物探行业应该抓紧发展先进技术,提高能源勘探的效率。

在物探行业中,地震勘探作为一个重要的手段,发挥着巨大的作用。

数据处理作为地震勘探的一大重要环节,所采用的各种方法和技术手段也一直在更新和进步。

在地震勘探处理方法研究中,地震数值模拟技术可以在室内完成地震数据模型的建立,并对其地震数据进行各种方法的处理,查看处理方法的效果和数据的好坏,另一方面,地震数值模拟进行正演获得的数据也可以作为反演的基础进行比对。

在地震数据处理的过程中,如何模拟地震波的传播便是需要解决的问题。

在二十世纪70年代开始采用显示差分格式来模拟地震波的传播。

由于有限差分法适用条件广,计算速度比较快,占用计算机内存少,编程比较容易实现,模拟精度相对较高而得到广泛应用。

但是有限差分法模拟地震波场时,由于计算机运算核心的限制,有限差分方法只能得到有限的数据点,地震波动方程只能是在有限差分方程中求得近似解,这时就考虑到人工边界问题,如果不对边界进行处理,波在通过边界时会产生反射,因此我们希望对添加的边界进行处理来消除这些反射。

在20世纪70年代,地球物理学界陆续采用了不同的边界条件来实现削弱地震波在通过边界时的反射,比如reynold边界、clayton边界、cerjan边界,以及后来提出的PML层边界条件,每种边界条件都在不同程度上实现了地震波通过边界时的衰减。

为了验证以上边界条件在数值模型的效果,在文章中,我们设计了一些简单的数值模型,给出了不同的边界条件,通过波形在通过不同边界条件时反射进行比较,观察每种方法衰减反射的效果。

波动方程求解方法

波动方程求解方法

常用的波动方程求解方法主要有以下几种:有限差分法、有限元法和伪谱法、积分方程法等。

1、有限差分方法由于适应性强,计算快速,因此是最先发展起来而且使用范围最广的数值方法,有限差分方法最大的弱点之一就是会产生数值频散。

有限差分法采用差分算式近似逼近偏导数运算,从而使波动方程的偏导数运算问题转化成差分代数问题,最后通过求解差分代数方程组得到近似解结果。

有限差分法的差分算式本身就是一种局部点运算,不需要考虑原函数中所求点值在邻域范围上的函数的变化情况,而只需要用到所求点值附近点上的值,所以能够很好的适用于复杂情况, 但是难保模拟精度。

有限差分方法有较高的空间域分辨率,而在频率域上分辨率反而会极低,稳定性同时还受到网格间距和时间步长的影响。

同时,虽然有限差分法还伴随有数值频散的问题,但是计算速度较快。

有限差分法目前主要有以下三大类:规则网格方程、弹性方程和交错网格方程。

有限差分法的具体操作可以分为两个部分:(1)用差分代替微分方程中的微分,将连续变化的变量离散化,从而得到差分方程组的数学形式:(2)求解差分方程组。

在第一步中,通过网格剖分法,将函数定义域分成大量相邻而不重合的子区域。

通常采用的是规则的剖分方式,最常用的是正方形网格。

这样可以便于计算机自动实现和减少计算的复杂性。

网格线划分的交点称为节点。

若与某个节点P 相邻的节点都是定义在场域内的节点,则P 点称为正则节点;反之,若节点P 有处在定义域外的相邻节点,则P 点称为非正则节点。

在第二步中,数值求解的关键就是要应用适当的计算方法,求得特定问题在 所有这些节点上的离散近似值。

目前最常用的两种有限差分方法包括:基于位移 波动方程的二阶中心差分法和基于一阶速度-应力波动方程的高阶交错网格法, 前者算法简单,易于实现,但差分精度具有局限性,最后得到的是节点上z x ,分量的位移离散近似值,后者算法稍复杂,但可以提高差分精度,最终得到的是节点上的位移速度离散近似值。

有限差分法模拟地震波场

有限差分法模拟地震波场

有限差分法模拟地震波场
1.原理
本次模拟中,根据《裂隙介质地震波场特征理论与应用》中的一种高精度有限差分算法。

通过建立一个200*200的数字网格,作为地下地质模型,并将其看作是一个均匀各向同性介质,其速度为2000m/s,将震源放在模型中心(100,100),以高阶有限差分为数值模拟方法实现地震波场模拟。

根据应力分析可得波动方程微分形式:
∂2P ∂x2+ ∂2P
∂z2
= 1
V2

2P
∂t2
其中,时间导数项采用二阶中心差分来逼近:
∂2P ∂t2=
P t+Δt−2P t+P t−Δt
Δt2
+o(Δt2)
空间导数项应用泰勒展开来构造,对于函数P以x为变量,将P x+ζΔx与P x−ζΔx在x点用Taylor级数展开,并将两式相减,消除高于二阶导数项,可得到函数的空间一阶,二阶导数,其中二阶导数的2n阶差分逼近式为:
∂2P ∂x2=1
Δx2
c i
d i
n
i=1
+ o(Δx2n)
式中c i为常数,d i=P x+ix−2P x+P x−ix,在本次模拟中,为了提高数值模拟的精度,在数值计算中,使用的是高阶有限差分,本次精度达到8阶,其中c i取c i=(1.6 -0.2 0.25 -0.01786),带入计算公式即可。

此次
实验精度为8阶,其有限差分格式满足稳定条件:V maxΔt
max⁡(Δx,Δz)=2000∗0.002
20

0.555。

2.结果
图1 地震波场传播过程
图2.3.1 震源子波图2.3.2 单炮地震记录。

地震波场数值模拟方法

地震波场数值模拟方法

第42卷第2期2003年6月石 油 物 探GE OPHY SIC A L PROSPECTI NG FOR PETRO LE UMV ol.42,N o.2Jun.,2003文章编号:100021441(2003)022*******地震波场数值模拟方法张永刚(中国石油化工股份有限公司科技发展部,北京100029)摘要:简要总结了地震波场数值模拟的各种方法的基本原理及其主要特点,对最近在该领域出现的一些方法和研究结果做了简要的阐述,并对比了各种方法的优缺点。

在此基础上提出了运用波动方程数值模拟作为基础,结合射线方法辅助识别波场类型,用于分析异常波的产生机理和出现特点的基本思想,这对复杂条件下的地震勘探具有指导和借鉴意义。

关键词:地震波场;数值模拟;射线追踪;有限元;伪谱法;正演模拟中图分类号:P63114+1 文献标识码:AOn numerical simulations of seismic w avefieldZhang Y onggang(Department of Science and T echnology Development,SI NOPEC,Beijing100029,China)Abstract:This paper reviews the principles and characteristics of various numerical simulations of seismic wavefield,and com2 pares the merits and defects of the simulations.S ome newly emerged methods and results are briefly discussed.The author pro2 poses to study the generation mechanism and characteristics of abnormal waves based on wave equation numerical simulation supplemented by ray tracing.K ey w ords:seismic wavefield;numerical simulation;ray tracing;finite element;pseudo2spectrum;forward m odeling 地震波场数值模拟是研究复杂地区地震资料采集、处理和解释的有效辅助手段,地震波场数值模拟的主要方法包括2大类,即波动方程法和几何射线法。

油气地球物理2007年总目次

油气地球物理2007年总目次
维普资讯






2 07年 0
总 目次

综 述・
地震 勘探 技术 概 述 … … … … … … … … … … … … …… …… … … … … … …… … 王有新 王延光 ( ) 1
浅谈 L u 集群在国内石油企业的应用现状及发展趋势 ………………………………… 戴 猛( ) ix n 1
快 照复 制技术 在 分布式 环 境 中应 用 … … … … … … … … …… … …… … … … … … … … … 侯 飞( ) 1
火成岩侵入形成的构造物理模拟实验 ………………… ………………………………… 孙志信( ) 2 横向各向同性介质纵波非双曲线 时差速度分析 ……… ……………… 杜启振 孙晶波 刘莲莲( ) 2
‘不口研 究 . 皇 综 合 九 ‘ 丌

准噶尔盆地车排子地区“ 三多” 成藏特征分析 …………………………………………… 宋传春( ) 1 济阳坳陷东部走滑构造形成机制 ……………………………………… 胡贤根 谭明友 张明振( ) 1 垦 7 断块东二段油藏建模研究 …………………… …………………………… 徐 强 宋建国( ) l 1 东营凹陷深层构造剖析及其地质意义 ………………………………… 贾红义 于建国 王金铎( ) 1 济阳坳陷深层构造层序划分 ……………………………………………………………… 贾红义( ) 3 济阳坳陷桩海地区下古生界深埋藏溶蚀特征及储集意义
广义希尔伯特变换在含噪信号边缘检测中的应用效果分析 …… … … 党志敏 贺振华 黄德济( ) 2
三维转换波资料处理方法研究及其应用 … … 毕丽飞 王延光 王 慧 乔玉雷 石建新 王秀玲( ) 2 叠前地震反演技术在气层预测中的应用 … …………………………… 王玉梅 孟宪军 慎 国强( ) 2

起伏地表弹性波传播有限差分法数值模拟

起伏地表弹性波传播有限差分法数值模拟

起伏地表弹性波传播有限差分法数值模拟董良国;郭晓玲;吴晓丰;马在田【期刊名称】《天然气工业》【年(卷),期】2007(027)010【摘要】有限差分是进行地震波传播数值模拟的最常用方法,但该方法处理起伏的自由边界比较困难.为此,通过对不同地形起伏情况下自由边界的具体分析,将整个二维空间离散点划分为24类,对每一类自由边界处的网格点选择了合理的表现方式,实现了起伏地表自由边界条件的数值化.该方法可以模拟出地表起伏情况下弹性波复杂的传播现象,为进行起伏地表地震波传播规律研究、山地地震勘探野外观测系统设计、山地地震勘探干扰波分析和识别、以及静校正研究提供了正演模拟工具.模拟实例表明,地形起伏引起面波、体波等地震波型之间的相互转化,产生了大量的散射P波、散射S波和散射面波,尤其是由沿地表传播的强能量面波,在地表起伏及近地表物性突变处产生了大量的强能量散射面波,同时也产生了相对较弱的散射P 波和散射S波,这是造成山地地震资料信噪比低的主要原因.【总页数】4页(P38-41)【作者】董良国;郭晓玲;吴晓丰;马在田【作者单位】同济大学海洋地质国家重点实验室;同济大学海洋地质国家重点实验室;同济大学海洋地质国家重点实验室;同济大学海洋地质国家重点实验室【正文语种】中文【中图分类】P61【相关文献】1.双相各向异性介质弹性波传播交错网格高阶有限差分法模拟 [J], 裴正林2.任意起伏地表弹性波方程交错网格高阶有限差分法数值模拟 [J], 裴正林3.起伏地表条件下2.5维声波方程有限差分法数值模拟 [J], 齐鹏;孙建国4.起伏地表弹性波传播的间断Galerkin有限元数值模拟方法 [J], 薛昭;董良国;李晓波;刘玉柱5.有限积分法与有限差分法在弹性波数值模拟中的对比分析 [J], 李明智;熊章强;张大洲因版权原因,仅展示原文概要,查看原文内容请购买。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

·8·
油气地球物理
2007年 7 月
为压制频散, 提高模拟精度, 只有减小采样间隔, 如 用常网格 4m 步长算法进行模拟, 内存需求将会扩 大 1 倍, 以上述( 1600m ×1600m ) 的模拟区域为例, 计算时间将会增大 1.2 倍左右。用文中所述变步长 算法进行模拟, 低速地表用 4m 小步长, 地表以下用 8m 大步长, 变网格步长算法的内存需求仅是常网格
" $ 2
!


=-


!t
P !




L P+PL

+s
!
!
( 3)
% $ 2


L=
!c



!
2+
!

!x !z
( 4)
收稿日期: 2007-03-23; 修订日期: 2007-04-27 作者简介: 李胜军, 男, 在读硕士研究生, 研究方向为地震 波 传 播 理 论 。 联 系 电 话 :( 0546) 8392055, E-mail: hdpulis@126.com, 通 讯 地 址 : ( 257061) 中国石油大学( 华东) 地球资信与信息学院。 * 中国石油大学( 华东) 研究生创新基金资助, 编号: S2006—06。
0 0.0
800
1600 ( m)

0.0
800
1600 ( m)
0.5
0.5
1.0
1.0
1.5
( s)
( a) 均匀网格
1.5
( s)
( b) 变网格
图 7 常网格与变网格得到的炮集记录
4 结论
上述分析与理论模型试算结果表明, 变网格步 长差分算法是一种高效的正演模拟算法。它能较好 地实现对低速地表、高速层夹层、复杂界面的数值模 拟, 成功地解决小网格步长差分算法内存需求量大、 计算效率低等问题。在特殊地区用增大( 减小) 步长 的模拟方法减小了内存需求, 提高了计算效率; 在特 定位置减小步长, 解决了对超薄层地质体、倾斜界面 的数值模拟效果不理想的问题: 解决了用常规高阶 有限差分或傅立叶变换法模拟存在的问题。
1 时间积分
均 匀 介 质 中 的 二 维 声 波 方 程 可 用 下 式 表 示[2]

!



=- L P+s
( 1)
!t
" # 2

22
- L =c
!2+
!
2( x, z, t) , 代表压 力项; c=c( x, z) , 代表速 度; s=s( s, z) , 代表震源函数; L2 为差分算子。 在 密 度 !=!( x, z) 变 化 的 情 况 下 , 常 用 的 是 V idale 给 出 的 公 式[5]
$# % P !x, z, t "= t sin !!L "h !t- ! "d ! g !x, z "( 5) 0L 该常规解通过契比雪夫展开近似为
$ ! "% &N/2
P !x,
z,

"


b2k+1

R iL!
Q2k+1
iL! R
g !x, z " ( 6)
式中: ! 为变量; L! 2 是空间差分算子 L2 的数值 逼近;
2007 年 7 月
油气地球物理 PETROLEUM GEOPHYSICS
第5卷 第3期
波动方程的变步长有限差分数值模拟 *
李胜军 1, 2) 孙成禹 1) 张玉华 1) 倪长宽 1) 1) 中国石油大学地球资源与信息学院; 2) 中石油勘探开发研究院西北分院
摘要: 有限差分算法是常用的正演模拟方法之一, 其包含的地震信息丰富, 且实现简单。传统的有限差分方法通常都 采用均匀网格步长, 在对含低速 / 高速介质、薄层 / 厚层介质的模型进行波场模拟时往往缺乏稳定性。文章介绍了一 种可以有效解决上述问题的变网格算法, 对常规有限差分法与变网格差分算法在内存需求、计算速率等方面的差别 进行了比较, 对变网格差分算法中的边界条件、时间积分的快速展开算法作了阐述, 进而总结了变网格算法的优点。 关键词: 变步长; 边界条件; 计算时间; 快速展开法; 数值模拟
·6·
油气地球物理
2007年 7 月
波动方程 的 时 间 积 分 可 通 过 K osloff等 提 出 的 快速展开法( R E M ) 来计算[6]。为 简便起见 , 仅对基 本方程( 1) 的快速展开法做一些介绍。对变密度的 情况可用类似的方法进行推导。
对震源项 s( x, z, t) =g( x, z) h( t) , 边界条件为 吸收边界时, 式( 1) 的常规解为
4m 步长算法的 65% 左右, 计算时间也缩短 45% 左 右。与常网格 8m 步长算法相比, 计算时间是常网格 8m 步长的 1.2 倍, 当然, 内存需求也相应有所扩大。 图 7( b) 是变网格步长算法的模拟结果( 为了对比 明显, 两图使用了相同的增益) , 可以看出, 频散得 到了较好的压制, 分辨率也明显提高。
在图 2 显示的网格中, 可使用这种变网格算法 来改变空间分辨率。在细网格区域不需要用傅立叶 变换法或高阶有限差分算子, 四阶或六阶精度的算 法就可以满足精度的要求。

如果 Z 方向上的步长通过其他参数来调整, 类 似的方程可被推导。一般通过函数 m( k) 来产生, 有

% 2
Dz
!P
" i, j
1 #2!z
$2 k

!

# $ P - 2P +P k i, j+2# k+n $ i, j i, j- 2# k+n $ ( 11)
为保证这种差分算法的稳定性, 压力必须内插 在水平方向两倍步长边界上的一些点上。对参数长 度 2N= 10 的那些点在图 5 中已列出。内插计算时 可使用特别精确且没有错误发生的三角内插算法。
$x 在 X 方向上是连续变化的, 用常规的有限差分 算法就可计算。但是在 Z 方向上, %z 是不连续的, 基 于连续变化的傅立叶变换法不再适用, 所以 Z 方向上的导数要用高阶有限差分法近似。


z0
图 1 低速地表模型差分网格分布示意图 x

图 2 低速夹层模型差分网格分布示意图 x 1500m /s
图 3 给出了简单的模型: 用固定小网格步长 !x 和 "z 采样, 可以很好地满足低速层的采样需求, 但 是这样造成了对高速层的过采样。为避免这种过采 样, 可以从某一深度 z0( 图 1) 开始采用双倍的步长 来采样。在修改的网格上, X 方向的导数用傅立叶法 或有限差分法计算, #x 通过采样定理来确定。由于
列表( 黑圆点位置属于精细网格循环; 方框点属于 双倍的步长循环点; 五角星位置点的导数被计算) 。 在每一列用不同的差分模式来计算在不同深度的二
阶导数。针对十阶精度有限差分相应的 m( k) 的值列 于表 1, 其他阶精度的 m( k) 值可用相应的方法计算。
z0 z
图 5 必须内插压力的网格点
3 算例分析
Qk 指修改了的契比雪夫多项式; 系数 bk 为
# bk=

J2k+1
!!R
" h
!t-
!
"d!
0R
( 7)
式中: Jk 为 k 阶精度的贝赛尔函数。 如果空间参数 - L! 2 的特征值很靠近实轴、R2 比
最大特征值大, 当 N>R·t 时, ( 6) 、( 7) 式按指数规 律收敛。参数 - L! 2的最大特征值由 X 方向上的差分

1 #!z $2 k

!

#P - k i, j+m(k)
$ 2Pi, j +Pi, j- m(k)
( 12)
差分因子
!
必须根据函数

m(
k)
计算, 这可通过
Fornberg( 1988a) 提出的一种算法来实现[9]。图 4 表
示了方程( 12) 在双倍步长区域差分算子 2N= 10的
情况下有限差分参数的变化情况。显示了 6 个网格
参数的最大特征值 Ax 和 Z 方向上的相应 特征值 Az 的和来决定。这里, R2 可由( 8) 式计算
! " 2 2
R =c
Ax !!x
"2 +
Az !!z "2
max
( 8)
2 空间导数计算
通常情况下, 勘探区域的地表较为复杂且介质 中波的传播速度较低。为较好地压制数值频散并提 高模拟结果的分辨率, 就要在整个计算区域采用小 网格步长来做差分数值模拟。这样就造成了对高速 区域的过采样, 浪费了计算机资源。为减小内存需 求和节省计算时间, 可以用图 1 所示方法进行模 拟。对于地层中有低速层的情形可以采用图 2 所示 方法来达到提高分辨率并且计算量增加不大的差 分算法。
为验证上述算法的正确性、有效性, 设计了图 6 所示的低速地表地质模型。数值模拟时, 震源设置在 ( 800m , 24m ) 处。
ab
c def

800
1600 ( m)

1000m /s
z0
图 4 过渡带的差分格式
表 1 图 4 显示的有限差分参数 m( k) 值
m( k) a





m( 0) 0

"
( 9)
对应的双步长网格定义如下


& 2
Dz
$P
% i, j

"
相关文档
最新文档