声波方程有限差分法数值模拟
声波方程数值模拟实验报告
声波方程数值模拟实验报告.基础理论知识需要的已知条件包括：1）震源函数地层速度（波速） 边界条件.2.2.2苇=v 2(諾二2)S(t)一 t :x :zv （x, Z ）是介质在点（x , z ）处的纵波速度，u 为描述速度位或者压力的波场，s （t ）为震 源函数。
为求式（4-1）的数值解，必须将此式 离散化，即用有限差分来逼近导数，用差商代替 微商。
为此，先把空间模型网格化（如图4-1所示）。
设x 、z 方向的网格间隔长度为h ，厶t 为时间采样步长，则有:X 二i^h（i 为正整数） z = j.）h（j 为正整数）^n ：t（n 为正整数）u ：j 表示在（i,j ）点，k 时刻的波场值。
k1. 2)3)>2u 2f cu〒=V p ^~2 t :x22w 2 w 厂二 V （—T :t ；x 声波方程的有限差分法数值模拟 对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述:2. 弹性波方程:二b s (t) 一 z22_w )；z 2)(4-1)将u i,j在（i,j）点k时刻用Taylor展式展开：k 1 kjUu i,j 7,j ■—ct将U i k j 」在(i,j)点k 时刻用Taylor 展式展开:k k4 k k5 k[U i2jui 2,j] 3[ui 4,j ui ・1,j] —?U i,j }S(t)*、(i -i °)**「(j - j °)(4-7)式中v(i, j)为介质速度的空间离散值，：h 是空间离散步长，=t 为时间离散步长，s(k)为震源函数，关于 s(k) 一般使用一个理论的雷克型子波代替，即:上式中，t 为时间，f 为中心频率，一般取为20-40HZ ， 为控制频带宽度的参数,(4-2)k 1k ；'UUi,jUi,2*「7 2 ；:t 2(4-3)将上两式相加，略去高阶小量, 整理得(i,j)点k 时刻的二阶时间微商为:2k 1kk J；：u u i,j -2u i,j u i,j.:t 242(4-4)对于空间微分，采用四阶精度差分格式,(以X 方向为例)即将U i*；j 、*Tj 分别在(i,j)点k 时刻展开到四阶小量，消除四阶小量并解出二阶微分得:-2~ u 11 k k 4 k k 5 k{—T7[u i 二 j +u id2,j ]+：[u i4j +U i*,j ] —=u i,j } 12 3 2:x L X 2(4-5)同理可得:2；=u 1.1 kk4 kk5 k一 2 人 2{一 石[Ui,j ，+Ui,j~2]+：[Ui,j 」+Ui,j^]—；U i,j}:z-z12 32(4-6)这就实现了用网个点波场值的差商代替了偏微分方程的微商，将上三个式子代入 (4-1)式中得:k 1 Ui,j kk 4= 2u i,j —Ui,jV :-1 1 k k 4 k k—{p [u=j U i/ yij —j：h 25 k H-U i,j }h 2 {12s (t )二 e(-2 f / )2t2cos2 二 ft(4-8)般取3-5。
有限差分法地震波传播数值模拟
2k
2 x
ω2
⎟⎟⎠⎞2
⎪⎫ ⎬ ⎪⎭
∂3P ∂z∂t 2
−
v2 4
∂3P ∂x2∂z
+
3v 4
∂3P ∂x2∂t
−
1 v
∂3P ∂t 3
=
0
Elastic Wave: (Bottom Boundary)
Utt = α 2U xx + β 2U zz + (α 2 − β 2 )Wxz Wtt = β 2Wxx + α 2Wzz + (α 2 − β 2 )U xz
2
MM 4
⎤ ⎥ ⎥
⎢⎢⎡CC12((MM
) )
⎤ ⎥ ⎥
M = M
6
⎥ ⎥ ⎥
⎢⎢C3(M
)
⎥ ⎥
⎢M⎥
2M ⎥
M⎦
⎢⎣CM(M ) ⎥⎦
⎡1⎤ ⎢⎢0⎥⎥ ⎢0⎥ ⎢⎢ M ⎥⎥ ⎢⎣0⎥⎦
⎡1
⎢ ⎢
13
⎢ 15
⎢ ⎢
M
⎢⎣12N −1
3 33 35
M 32 N −1
5 53 55
M 52 N −1
1992，1994年，Tessmer et al.在模拟二维以及三维不规则 地表面波时，同样也是使用的如上吸收边界。
采取的措施
z 高阶差分解法--提高计算精度，减小数值频散 z 采用基于特征分析方法得到的吸收边界条件
数值频散问题------高阶差分解法
声波：
∑ ∂2 f
∂x 2
=1 Δx 2
M
Cm(M )[ f (x + mΔx) − 2 f (x) +
m=1
f (x − mΔx)] + o(Δx2M )
声波方程有限差分正演
题目：使用Ricker 子波，刚性边界条件，并且初值为零，在均匀各向同性介质条件下，利用交错网格法求解一阶二维声波方程数值解。
解：一阶二维声波方程：22222221zPx P t P c ∂∂+∂∂=∂∂ (1)将其分解为：21P c t Px P z x z x z V V x z V tV t ∂∂∂⎧=+⎪∂∂∂⎪∂∂⎪=⎨∂∂⎪∂∂⎪=⎪∂∂⎩(2)对分解后的声波方程进行离散，可得到：112211,-1,,,122[]N n n n n m i m j i m j xi j xi j m t VVc P P h +-+---=∆=+-∑ 112211,1,,,122[]Nn n n n m i j m i j m zi j zi j m t VV c P P h +-++---=∆=+-∑ 1111212222,,m 1,,,,11[]Nn n n n n n i ji jmxi j xi m j zi j m zi j m m tc PP cVVVVh+++++++-+--=∆=+-+-∑h z x =∆=∆针对公式(1)，使用二阶中心差商公式：2P(,,1)2(,,)(,,1)i j n P i j n P i j n t +-+-∆222(1,,)2(,,)(1,,)(,1,)2(,,)(,1,)P i j n P i j n P i j n xc P i j n P i j n P i j n z +-+-⎧⎫+⎪⎪⎪⎪∆=⎨⎬+-+-⎪⎪⎪⎪⎩∆⎭(3)变形：P(,,1)=2(,,)(,,1)i j n P i j n P i j n +--2222(1,,)2(,,)(1,,)t (,1,)2(,,)(,1,)P i j n P i j n P i j n xc P i j n P i j n P i j n z +-+-⎧⎫+⎪⎪⎪⎪∆+∆⎨⎬+-+-⎪⎪⎪⎪⎩∆⎭(4)对离散格式作时间和空间三重Fourier 变换：0P(,,)(,,)x z i j n P k k w ↔ ，0P(,,1)(,,)*exp()x z i j n P k k w iw t +↔∆0P(1,,)(,,)*exp(k )x z x i j n P k k w i x +↔-∆，0z P(,1,)(,,)*exp(k )x z i j n P k k w i z +↔-∆对公式(4)进行Fourier 变换：2222exp()2exp()h exp()2()exp()2exp()h x x z z ik x ik x iw t iw t t c ik z ik z -∆-+∆⎡⎤+⎢⎥∆=--∆+∆⎢⎥-∆-+∆⎢⎥⎢⎥⎣⎦2222exp()2exp()h exp()2()=exp()2exp()h x x z z ik x ik x iw t iw t t c ik z ik z -∆-+∆⎡⎤+⎢⎥∆-+-∆∆⎢⎥-∆-+∆⎢⎥⎢⎥⎣⎦222222sin sin 22sin (2x z k x k zw tt c h∆∆+∆=∆） (5) 公式(5)右端必须满足下列条件：22222sin sin 220(x z k x k zt c h∆∆+≤∆≤）1 取x k 和z k 最大值，即=x x k π∆，z =k z π∆，则有：22220t c h≤∆≤1因此tc ∆≤即为所求得的稳定性条件。
二维声波方程交错网格有限差分数值模拟研究
在对介质模型进行离散化处理的过程中，网格
是一种常用手段。对波动方程进行网格离散，可以
利用交错网格的差分形式。交错网格就是把速度和
应力分配到两套不同的网格中，这样可以使速度、应
力得到很好的耦合［２］。利用交错网格有限差分法对
一阶速度—应力波动方程进行求解时，应力、速度等
分量在模型交错网格节点中的位置分布如图 １所示。
－Ｕ ］ ｋ＋１／２ ｉ，ｊ－ｎ＋１
其中 ｘ＝ｉΔｘ，ｚ＝ｊΔｚ，ｔ＝ｋΔｔ，ｉ、ｊ、ｋ分别表示空 间和时间网格点。Ｕｋｉ，＋ｊ１／２，Ｗｋｉ＋＋１１／／２２，ｊ＋１／２，Ｐｋｉ，＋ｊ＋１１／２／２，Ｑｋｉ，ｊ＋１／２ 和 Ｓｋｉ＋１／２，ｊ分别是速度 Ｖｘ、Ｖｚ与应力 σｘｘ、σｚｚ、σｘｚ的离 散值。
数的问题转化为求解网格节点上的差分方程组的问
题，得到数值解。在波动方程网格离散化的过程中，
可以利用交错网格的差分形式。
１．１ 均匀各向同性介质二维声波方程
均匀各向同性介质二维声波方程可表示为：
２ｕ（ｘ，ｘｙ２，ｚ，ｔ）＋２ｕ（ｘ，ｚｙ２，ｚ，ｔ）＝
ｖ２（１ｘ，ｚ）２ｕ（ｘ，ｔｙ２，ｚ，ｔ）
（１）
应力 Ｐｘｘ： Ｐｋｉ＋＋１１／２，ｊ＝Ｐｋｉ＋１／２，ｊ＋Ｃ１１ΔΔｘｔｎΣＮ＝１ＣＮｎ［Ｕｋｉ＋＋ｎ１，／ｊ２－Ｕｉｋ－＋ｎ１＋／２１，ｊ］＋
Ｃ１３ΔΔｚｔｎΣＮ＝１ＣＮｎ［Ｗｉｋ＋＋１１／／２２，ｊ＋ｎ－１／２
－Ｗ ］ ｋ＋１／２ ｉ＋１／２，ｊ－ｎ１／２
应力 Ｑｚｚ： Ｑｉｋ＋＋１１／２，ｊ＝Ｑｋｉ＋１／２，ｊ＋Ｃ１３ΔΔｘｔｎΣＮ＝１Ｃ（ｎＮ）［Ｕｉｋ＋＋ｎ１，／ｊ２－Ｕｉｋ－＋ｎ１＋／２１，ｊ］＋
Ｃ３３ΔΔｚｔｎΣＮ＝１Ｃ（ｎＮ）［Ｗｋｉ＋＋１１／／２２，ｊ＋ｎ－１／２－Ｗｉｋ＋＋１１／／２２，ｊ－ｎ＋１／２］
声波方程有限差分数值模拟的变网格步长算法
声波方程有限差分数值模拟的变网格步长算法声波方程是描述声波在介质中传播的方程，在进行数值模拟时需要使用有限差分方法来近似求解。
有限差分方法将连续的空间和时间离散化，通过在离散的网格节点上计算声压场的数值来模拟声波的传播过程。
在有限差分数值模拟中，选择合适的网格步长对数值结果的精度和计算效率具有重要影响。
变网格步长算法是一种能够根据需要在不同区域自适应地调整网格步长的技术。
在声波方程有限差分数值模拟中，声波会在不同介质中以不同的速度传播，因此在网格步长选择上需要考虑介质的变化。
当声波在介质变化剧烈的区域传播时，使用较小的网格步长可以提高模拟的精度。
而在介质变化缓和的区域，使用较大的网格步长能够减少计算量。
变网格步长算法的基本思想是通过对声波传播区域进行划分，然后根据介质的变化情况调整不同区域的网格步长。
具体步骤如下：1.对传播区域进行划分：根据介质的变化情况，将声波传播区域划分为多个区域，每个区域具有不同的网格步长。
2.确定初始网格步长：在初始时，可以根据经验或者初步模拟结果选择一个合适的网格步长。
3.计算声波传播：在每个区域内分别进行声波方程的有限差分数值模拟，使用当前区域的网格步长进行计算。
4.判断误差与精度：通过计算得到的声压场数值和预期精度进行比较，如果达到要求，则结束模拟；否则，进行下一步。
5.调整网格步长：根据当前区域的模拟误差情况，调整该区域的网格步长，使其更适应声波在当前区域的传播特性。
可以根据误差大小和梯度等因素进行调整。
6.重复模拟：根据调整后的网格步长，重新进行声波方程的有限差分数值模拟，并返回步骤4通过以上步骤，可以实现声波方程有限差分数值模拟的变网格步长算法。
这种算法能够根据介质的变化情况自适应地调整网格步长，提高模拟的精度和计算效率。
然而，变网格步长算法的实现并不简单，需要合理设置划分区域和调整步长的策略，并在计算过程中不断优化和调整，才能达到较好的数值模拟效果。
声波有限差分数值模拟
声波方程的离散化
时间离散化
将时间轴划分为一系列离散的时 间点，以便在每个时间点上求解 声波方程。
空间离散化
将连续的空间区域划分为一系列 离散的网格点，以便在每个网格 点上求解声波方程。
声波方程的有限差分近似
声波速度的近似
使用有限差分法近似声波的速度，以 便在离散化的空间和时间网格上求解 声波方程。
调试技巧
使用断点、单步执行和变量监视等技巧，定位和 修复代码中的错误。
04 声波有限差分法的应用
在地震波模拟中的应用
总结词
地震波模拟是声波有限差分法的重要应用领域之一，通过模拟地震波在地壳中的传播过程，可以用于地震勘探、 地震工程和地震安全性评估等领域。
详细描述
声波有限差分法能够模拟复杂的地震波场，包括反射、折射、散射等现象，为地震勘探提供精确的数据。在地震 工程中，该方法可用于评估结构的抗震性能和设计地震防护措施。此外，在地震安全性评估中，声波有限差分法 可以模拟地震对地下设施的影响，为设施的安全设计和运营提供依据。
谱方法
将谱方法与有限差分模拟 相结合，提高数值计算的 稳定性和精度。
无网格方法
探索与无网格方法（如移 动最小二乘法、径向基函 数等）的结合，实现更加 灵活的离散方式。
THANKS FOR WATC限差分法需要将连续的空间离散化为网格，这可能导致 结果的精度和可靠性受到网格大小和形状的影响。
边界条件处理
在处理具有复杂边界条件的声波问题时，有限差分法可能 需要采用特殊的边界条件处理方法，增加了模拟的复杂性 和难度。
数值色散误差
有限差分法在模拟声波传播过程中会产生数值色散误差， 即模拟结果与真实解之间的误差随模拟步数的增加而增大， 这可能影响模拟结果的精度和可靠性。
超声波在各向同性固体中传播的数值模拟_肖开丰
收稿日期:2007-08-09 修回日期:2007-08-17 第25卷 第9期计 算 机 仿 真2008年9月 文章编号:1006-9348(2008)09-0292-03超声波在各向同性固体中传播的数值模拟肖开丰,宋文爱(中北大学信息与通信工程学院,山西太原030051)摘要:介绍以应用动量守恒和应力—应变关系为各向同性弹性介质进行微变形建立超声波传播方程,并利用一种有限元方法和差分法,分析超声波(p 波)在各向同性介质中传播,及在裂缝上的散射结果。
并且通过计算机程序模拟二维超声波的传播,散射的数值模拟与几何理论结果是相符合的。
因此超声检测仿真软件是可以预测超声检测过程中的波形。
超声探头的建模发展趋势是利用F D M 和F E M 方法研究超声检测的仿真软件,分析结果将广泛应用到检测过程中。
关键词:有限元;超声波;波传播中图分类号:T P 391.9 文献标识码:BN u m e r i c a l S i m u l a t i o no f t h e P r o p a g a t i o n o f U l t r a s o n i cWa v e i na nI s o t r o p i c S o l i dX I A OK a i -f e n g ,S O N GWe n -a i(S c h o o l o f I n f o r m a t i o n a n d C o m m u n i c a t i o n E n g i n e e r i n g ,N o r t h U n i v e r s i t y o f C h i n a ,T a i y u a n S h a n x i 030051,C h i n a )A B S T R A C T :T h i s p a p e r p r e s e n t s t h e a p p l i c a t i o n o f t h e e q u a t i o n s o f c o n s e r v a t i o n o f m o m e n t u ma n d t h e s t r e s s -s t r a i n r e l a t i o n s f o r a n i s o t r o p i c e l a s t i c m e d i u mu n d e r g o i n g i n f i n i t e s i m a l d e f o r m a t i o nt o e s t a b l i s ha n u l t r a s o n i c w a v e p r o p a g a -t i o ne q u a t i o n ,a n du s e st h ef i n i t ee l e m e n t m e t h o d(F E M)a n dt h ed i f f e r e n c em e t h o df o r a n a l y s i i n gu l t r a s o n i c (P-w a v e )w a v e p r o p a g a t i o ni ni s o t r o p i c m e d i a ,a n dt h ec r a c k s s c a t t e r i n gr e s u l t s .A n dt h et w o -d i m e n s i o n a l u l t r a s o n i c w a v e p r o p a g a t i o n h a s b e e ns i m u l a t e db y c o m p u t e r p r o g r a m a n d t h e s c a t t e r e d w a v e i s w e l l e x p l a i n e db ya g e o m e t r i c a l t h e o r y .T h e r e f o r e u l t r a s o n i c t e s t i n g s i m u l a t i o n s o f t w a r e c a n p r e d i c t t h e c o u r s e o f u l t r a s o n i c t e s t i n g .T h e m o d e l s o f u l -t r a s o n i c p r o b e s w e r e d e v e l o p e d ,a n dt h e y w e r e c o m b i n e dw i t h t h e F E M a n d t h e F D M t o d e v e l o p a c o m p l e t e s i m u l a t o r o f u l t r a s o n i c t e s t i n g (U T ).T h e s e f e a t u r e s e n a b l e t h e m o d e l i n g o f a w i d e r a n g e o f u l t r a s o n i c t e s t i n g .K E Y WO R D S :F i n i t e e l e m e n t ;U l t r a s o n i c w a v e ;P r o p a g a t i o n1 引言固体中超声的散射是声学的基本问题之一,同时是超生无损检测的核心问题。
声波方程有限差分数值模拟的变网格步长算法
声波方程有限差分数值模拟的变网格步长算法声波方程有限差分数值模拟是一种常用的声波传播模拟方法，可以在计算机上通过数值计算求解声波传播的过程。
在进行这种数值模拟时，常常需要选择合适的网格步长，以保证计算结果的准确性和计算效率。
本文将介绍一种变网格步长算法，用于优化声波方程有限差分数值模拟的计算。
声波方程可以用下面的形式表示：∂^2p/∂t^2=c^2∇^2p其中p是声场变量，t是时间，c是声速，∇^2是Laplace算子。
为了将声波方程用有限差分方法进行离散化计算，我们需要将空间和时间分别离散化。
首先，将空间离散化为网格，在每个网格点上计算声场的值。
其次，将时间离散化为离散的时间步长，通过迭代计算不同时间步长上的声场分布。
为了保证计算结果的准确性，网格步长应当满足Nyquist采样定理的要求。
即网格步长应小于声波的最小波长的一半。
根据声波方程的性质，我们可以通过声速和最高频率来估计声波的最小波长。
然后，我们可以根据最小波长来选择合适的网格步长。
然而，在实际的声波传播计算中，声场的变化往往不是均匀的。
有些区域的声场变化较大，而其他区域的声场变化较小。
如果我们在整个计算区域都采用较小的网格步长，将会造成计算资源的浪费。
因此，需要一种方法能够根据声场的变化情况来自适应地调整网格步长。
变网格步长算法就是一种能够根据声场变化情况自动调整网格步长的算法。
其基本思想是根据声场在不同网格上的变化率来决定每个网格上的网格步长。
具体的算法步骤如下：1.初始化：选择一个合适的初始网格步长。
通常可以选择根据声波的最小波长来确定。
2.计算网格步长：在每个时间步长上，对于每个网格点，计算其周围网格点上的声场变化率。
常用的方法是计算声场在三个相邻时间步长上的差分值，然后取绝对值并求平均。
根据声场变化率，调整当前网格点上的网格步长。
变化率大的网格点应该有更小的网格步长，而变化率小的网格点则可以有更大的网格步长。
3.更新声场：根据调整后的网格步长，更新所有网格点上的声场值。
声波方程有限差分数值模拟的变网格步长算法
工程地球物理学报
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
声波在非线性介质中传播行为的数值模拟与分析
声波在非线性介质中传播行为的数值模拟与分析声波是一种机械波，它是通过介质中分子的振动传播的。
当声波传播过程中介质的特性发生变化，例如介质的非线性性质，声波的传播行为将会发生一系列的变化。
为了更好地理解声波在非线性介质中的传播行为，数值模拟和分析成为一种重要的研究方法。
一、声波在非线性介质中的基本特性声波在非线性介质中的传播行为与线性介质有所不同。
在线性介质中，声波的传播速度与频率无关，而在非线性介质中，声波的传播速度将与频率有关。
这种频率依赖性导致了声波在非线性介质中的传播现象，例如声速的失真和声波的各种非线性效应。
二、数值模拟方法为了模拟声波在非线性介质中的传播行为，数值模拟方法成为一种有效的研究手段。
数值模拟方法可以通过计算机模拟声波在非线性介质中的传播过程，从而得到声波的传播速度、频率特性以及非线性效应等信息。
1. 声波方程的数值解声波方程是描述声波在介质中传播的基本方程。
在非线性介质中，声波方程可以表示为非线性的偏微分方程。
为了求解这个方程，可以使用数值方法，例如有限差分方法或有限元方法，将偏微分方程离散化为代数方程组，然后通过迭代求解得到声波的数值解。
2. 非线性介质参数的数值计算非线性介质的参数对声波的传播行为有重要影响。
为了进行数值模拟，需要准确地计算非线性介质的参数。
这可以通过实验测量或者理论计算来获得。
例如，可以利用声压级对声波的非线性参数进行测量，或者使用分子动力学模拟等方法计算非线性介质的分子结构和相互作用。
三、声波在非线性介质中的数值模拟与分析通过数值模拟和分析，我们可以研究声波在非线性介质中的传播行为。
例如，可以模拟声波在非线性介质中的传播速度随频率变化的情况。
结果显示，随着频率的增加，声波的传播速度将减小，这是由于非线性介质中的频率依赖性导致的。
此外，数值模拟还可以研究声波在非线性介质中的非线性效应，例如声波的失真、声波的非线性色散和声波的非线性吸收等。
通过数值模拟，我们可以定量地分析这些非线性效应对声波传播行为的影响，从而更好地理解声波在非线性介质中的传播机制。
声波有限差分数值模拟
2P
x2
1 x2
[
5 2
Pn i, j,k
4 3
P P n i 1, j,k
n i 1, j,k
1 12
P P n i2, j,k
n i2, j,k
]
2P y2
1 y2
[
5 2
Pn i, j,k
4 3
P P n i, j 1,k
n i, j 1,k
1 12
P P n i, j 2,k
u(xΒιβλιοθήκη a)a2 2!2u x2
a4 4!
4u x4
(n为偶数)
an n!
nu xn
o(an )
（2-4）
二、声波方程的差分近似
（2-4）中分别取 a x a 2 x a 3 x
并做 o(a6 ) 误差截断得：
u(x x) 2u(x) u(x x) x2 2u x4 4u x6 6u
用类似的方法可以求得
2u 的四阶误差截断表达式：
x2
2u
x2
1 x2
{
5 2
u(x)
4 3
u(x
x) u(x
x) 1 u(x 2
12
x) u(x 2
x)}
2u 的二阶误差截断表达式：
x2
（2-8） （2-9）
2u
x2
1 x
2
u
(
x
x) 2u(x) u(x x)
（2-10）
二、声波方程的差分近似
2u x2
b1 x2
1 x2
3r1
3 10
r2
1 45
r3
二、声波方程的差分近似
2u
声学波传播过程的数值模拟分析
声学波传播过程的数值模拟分析声学波传播是研究声波在不同介质中传播规律的一门学科。
通过数值模拟分析声学波的传播过程，我们可以更好地理解和预测声波在不同介质中的行为，为声学相关领域的研究和应用提供有力支持。
声学波传播的数值模拟分析首先需要确定所研究的问题，如声源的特性、介质的物理参数以及边界条件等。
然后，通过建立合适的数学模型和方程组，利用计算机进行数值计算和解析。
最后，根据模拟结果对声波传播过程进行分析和评估。
在声学波传播的数值模拟分析中，常用的方法包括有限差分法（FDM）、有限元法（FEM）和边界元法（BEM）等。
这些方法各有特点，可以根据具体问题和需求选择合适的方法进行模拟分析。
以有限差分法为例，它是一种离散化计算的方法。
首先，将声波传播问题的连续域转化为离散的有限差分网格，将时间和空间分割成小块。
然后，根据声学波动方程将声场的变化量用差分的形式表示。
最后，通过数值计算和迭代求解差分方程组，得到声场在各个时间和位置的数值解。
有限差分法的数值模拟分析具有一定的深度。
通过改变差分网格的分辨率，我们可以探究声波传播过程中的细节和特征。
例如，在分析声波在不同介质中的传播速度和衰减率时，可以通过调节网格大小和时间步长的方法来探讨它们对声波传播的影响。
此外，还可以研究声波在复杂介质结构中的传播规律，如声波在不同形状和密度的障碍物中的散射和衍射现象。
声学波传播的数值模拟分析还可以应用于声波在医学成像和工程设计中的研究。
例如，在医学领域中，数值模拟分析可以用于研究超声波在人体组织中的传播规律，以帮助医生进行准确的诊断和治疗。
在工程设计中，数值模拟分析可以用于研究声波在复杂环境中的传播特性，如建筑物中的声学设计和噪音控制。
当然，声学波传播的数值模拟分析也存在一些挑战和限制。
首先，模拟的精确度和计算效率之间存在着一定的平衡。
增加模拟的精度会导致计算量的增加，而过于追求计算效率可能会牺牲模拟的准确性。
其次，模拟结果往往需要与实验数据进行对比验证，以确保模拟的可靠性。
声波方程数值模拟矩形网格有限差分系数确定法
声波方程数值模拟矩形网格有限差分系数确定法梁文全;王彦飞;杨长春【摘要】压制数值频散是有限差分方法的关键问题之一.目前压制数值频散的方法大多假设不同方向空间偏导数的空间步长相同,导致算法精度低,计算效率低.为此,提出使用线性方法压制声波方程矩形网格有限差分算子的数值频散,并进行了稳定性分析、频散分析和数值模拟.通过频散分析和数值模拟,验证了本文方法能够有效压制矩形网格有限差分数值频散,相较于泰勒展开方法和最小二乘方法,线性方法计算有限差分系数的效率更高,可以替代传统的正方形有限差分网格和相应的系数用于声波方程数值延拓.【期刊名称】《石油地球物理勘探》【年(卷),期】2017(052)001【总页数】7页(P56-62)【关键词】声波模拟;时间—空间域;有限差分格式;矩形网格;数值频散【作者】梁文全;王彦飞;杨长春【作者单位】龙岩学院资源工程学院,福建龙岩364000;中国科学院地质与地球物理研究所,中国科学院油气资源研究重点实验室,北京100029;中国科学院地质与地球物理研究所,中国科学院油气资源研究重点实验室,北京100029【正文语种】中文【中图分类】P631地震波数值模拟延拓是地震学和勘探地震学的重要基础，在理论和应用中都发挥了重要的作用［1］。
地震波数值模拟延拓的方法包括有限差分法、伪谱法、传统有限元法、谱元法、有限体积法、间断Galerkin法等［2，3］。
有限差分方法因为计算效率高、所需内存小、实现简单而广泛应用于地震正演研究，同时是叠前逆时深度偏移成像、全波形反演、速度分析技术迅速发展的基础［4，5］。
数值频散是对波动方程的时间和空间偏导数的离散化造成的，使相速度变成了网格间距的函数，降低了模拟的精度。
在地震波数值模拟中，一般使用正方形（正方体）差分网格进行数值模拟［6，7］。
在实际地震资料处理中，需要使用矩形网格有限差分格式。
Chen［8］基于平均导数法提出一种新的9点有限差分算子，可以解决频率—空间域数值模拟中不同方向空间采样间距不同的问题，并提高了频率域数值模拟的精度。
一维声波方程有限差分模拟
一维声波方程有限差分模拟一维声波方程是$$\frac{\partial^2 p}{\partial t^2} = c^2 \frac{\partial^2p}{\partial x^2}$$其中 $p(x,t)$ 是在位置 $x$ 和时间 $t$ 上的压力，$c$ 是声速。
为了进行数值模拟，我们需要使用有限差分方法来近似空间和时间上的偏导数。
一个常用的方法是使用二阶中心差分，即$$\frac{\partial^2 p}{\partial t^2} \approx \frac{p_j^{n+1} -2p_j^n + p_j^{n-1}}{\Delta t^2}$$$$\frac{\partial^2 p}{\partial x^2} \approx \frac{p_{j+1}^n -2p_j^n + p_{j-1}^n}{\Delta x^2}$$其中 $p_j^n$ 表示在位置 $j\Delta x$ 和时间 $n\Delta t$ 上的压力。
将这两个式子代入原方程得到$$\frac{p_j^{n+1} - 2p_j^n + p_j^{n-1}}{\Delta t^2} = c^2\frac{p_{j+1}^n - 2p_j^n + p_{j-1}^n}{\Delta x^2}$$整理后可得$$p_j^{n+1} = 2p_j^n - p_j^{n-1} + c^2 \frac{\Delta t^2}{\Deltax^2} (p_{j+1}^n - 2p_j^n + p_{j-1}^n)$$这是通过有限差分方法得到的一维声波方程的数值模拟公式。
我们可以从初始状态 $p_j^0$ 和 $p_j^1$ 开始迭代，按照上述公式计算 $p_j^{n+1}$ 直到达到所需的时间步数。
需要注意的是，为了保证数值稳定性，需要满足 $\Delta t \leq \frac{\Delta x}{c}$ 的条件。
有限差分法不同边界条件下的数值模拟
有限差分法不同边界条件下的数值模拟文章介绍了地震数据处理中所使用的数值模拟法，对采用有限差分法所使用不同边界条件处理方式进行了数值模拟，通过波场快照直观的得出了不同的边界吸收条件的吸收效果，对结果进行了对比，分析总结了各种方法的优缺点。
标签：数值模拟；有限差分；边界条件随着近年来国家宏观经济调控，经济增长的速度逐步减缓，能源行业受此影响最为严重，许多煤矿是在亏损的情况下生产，直接导致了地质行业投入的减少。
物探行业压力也越来越大，物探行业应该抓紧发展先进技术，提高能源勘探的效率。
在物探行业中，地震勘探作为一个重要的手段，发挥着巨大的作用。
数据处理作为地震勘探的一大重要环节，所采用的各种方法和技术手段也一直在更新和进步。
在地震勘探处理方法研究中，地震数值模拟技术可以在室内完成地震数据模型的建立，并对其地震数据进行各种方法的处理，查看处理方法的效果和数据的好坏，另一方面，地震数值模拟进行正演获得的数据也可以作为反演的基础进行比对。
在地震数据处理的过程中，如何模拟地震波的传播便是需要解决的问题。
在二十世纪70年代开始采用显示差分格式来模拟地震波的传播。
由于有限差分法适用条件广，计算速度比较快，占用计算机内存少，编程比较容易实现，模拟精度相对较高而得到广泛应用。
但是有限差分法模拟地震波场时，由于计算机运算核心的限制，有限差分方法只能得到有限的数据点，地震波动方程只能是在有限差分方程中求得近似解，这时就考虑到人工边界问题，如果不对边界进行处理，波在通过边界时会产生反射，因此我们希望对添加的边界进行处理来消除这些反射。
在20世纪70年代，地球物理学界陆续采用了不同的边界条件来实现削弱地震波在通过边界时的反射，比如reynold边界、clayton边界、cerjan边界，以及后来提出的PML层边界条件，每种边界条件都在不同程度上实现了地震波通过边界时的衰减。
为了验证以上边界条件在数值模型的效果，在文章中，我们设计了一些简单的数值模型，给出了不同的边界条件，通过波形在通过不同边界条件时反射进行比较，观察每种方法衰减反射的效果。
有限差分法地震波场数值模拟中的吸收边界条件
（7）
⎧P(0, j, k +1) = (2 − 2A − A2 )P(0, j, k) + 2A(1+ A)P(1, j, k)
⎪ ⎪
− A2P(2, j, k) + (2A −1)P(0, j, k −1) − 2AP(1, j, k −1)
⎪⎪P(M , j, k +1) = (2 − 2A − A2 )P(M , j, k) + 2 A(1+ A)P(M −1, j, k)
Chen Jingguo
China University of Geosciences (Beijing) (100083) E-mail: chenjg_cugb@ Abstract
The Finite Difference Method (FDM) is one of the most important methods of seismic wave field numerical modeling. But when we model seismic wave field with FDM, absorption boundary problem should be considered because the machine can not model the geological body of infinite boundary. This paper presents 2-D acoustic wave equation sample, through contrasting effectiveness of non-absorbing with the 1st order and 2nd order absorbing boundary conditions, illuminate the artificial boundary reflection (pseudo-reflection) problem and show the importance of the absorbing boundary conditions. Key words: finite difference method, seismic wave field, numerical modeling, absorbing boundary conditions, pseudo-reflection
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
}
main() {
FILE *fp; int i,j,k; float v[I][J],Dt[I][J];
//定义文件
double **u1,**u2,**u3,**u4,**u5,w[KMAX];
通过以上的计算，可以将 k=130 和 k=240 时波场的数据进行导出，用 Matlab 绘出图像。其中，在 Matlab 里用到的成像程序如下（这里只列举第一张波前快照的程序，第二张的与之相同）： clear all; load wavefront1.dat;
m=wavefront1(1); n=wavefront1(2); len=length(wavefront1)-2; nn=len/m; if n~=nn; disp('The data file length error!'); return; end
u k 1 i, j
2uik, j
uik,
1 j
v2t 2 h 2
{
1 12
[
uk i2,
j
uk i2,
j
]
4 3
[uik1,
j
uk i 1,
j
]
5 2
uik,
j
}
v 2t 2 h 2
{
1 12
[
uik,
j2
uik, j2 ]
4 3
[uik,
j
1
主要集中在前部，大致上 t>35 的时候 S(t)都为零值， 故在绘制 S(t)图像时 k 取 100，来比较清晰得呈现雷 克子波。本实验 不涉及吸 收边界条件 的使用。这样
声波方程数值模 拟所需的（ 1）震源函 数（2）地层
速度（波速）（3）边界条件均定义完全。
为求（1）式的数值解，必须将此式离散化（包 括时间离散和空间离散，这里 dt 取 0.002，dh 取 4， 即用有限差分来逼近导数，用差商代替微商。为此，
图4
图5
通过调节增益（x(j)=(i-1)*dx+wavefront((i-1)*n+j+2)*12*dx，dx 前面的 12），对波场
值放大一定倍数，可以使能量小的增大，使成像清晰。从下图中可清楚看到波前的图像，数层
的地震波随时间的推移由震源开始向外扩散。k=130 时，大约是时间为 260ms 时的图像，此时
关键词：声波方程 数值模拟 地震波 有限差分 引言：地震波场模拟即地震正演，是指已知模型结构，通过物理或数值计算的
方法模拟该地质结构下的地震波的传播，最终合成地震记录，也可以认为其是野 外数据采集过程的室内再现。物理模拟花费昂贵，人们一般采用比较经济的数值 模拟技术。地震波场数值模拟是在给定数学模型（如弹性波方程，声波方程等）、 震源和地下几何界面、物性参数（岩层密度、速度等）情况下，研究弹性波或声 波的传播规律。地震波波动方程数值模拟方法主要包括克希霍夫积分法、傅里叶 变换法、有限元法和有限差分法等。相对于上述几种方法，有限差分法是一种更 为快速有效的方法。虽然其精度比不上有限元法，但因其具有计算速度快，占用 内存较小的优点，在地震学界受到广泛的重视与应用。本实验给定大小两个模型， 对于小模型，整个区域的速度值可设为常数，即只有一种介质，将震源点放在模 型中间，分别记录两个时刻的波前快照（即该时刻区域内所有网格点的波场值）。
由上面的程序可得震源离散后的图像，如下图：
图2
小模型波前快照
选用一个小模型，即 200*200 的范围，视为完全弹性模型，为一层介质，在介质中取速 度为定值，本实验中取为 2500，将震源放在模型中央，分别记录两个时刻的波前快照，（即 区域内所有网格点的波场值）。第一时刻为地震波还未传播到边界上的某时刻，第二时刻为 地震波已经传播到边界上的某时刻，体会其人工边界反射。
别对应哪些波。
实验内容：
对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述：
2u t 2
v
2
(2u x 2Fra bibliotek2u z 2
)
S (t )
（1）
其中 v(x,z)是介质在点(x , z)处的纵波速度，本实验中 v 为常数 2500，u 为描述速度位或者压力的波场，
s(t) 为震源函数，本实验选用雷克子波（零相位子波）
e s(t) (2f / )2t2 cos 2ft
上式中，t 为时间，fm 为中心频率，一般取为 20-40HZ，本实验取 20，γ为控制频带宽度的参数， 一般取 3-5,这里取 4，在实际计算过程中，需把此震 源函数离散，参 与波场计 算。由图像 可知子波能量
图3
根据网格设计可知，网格间距为 10m，则 200*200 的矩阵实际上大致是个 2000*2000 平方米的区域， 这样震源在（100，100）即中心点，即在地质剖面上，宽度为 1000m，深度为 1000m，假设介质均匀各向 同性，且边界为自由界面，所以在介质中速度 v 以一常数 2500m/s 的传播速度传播，应该在 400ms 左右时 到达边界。
uik, j1 ]
5 2
uik,
j
}
s(t) * (i i0 ) * ( j j0 )
（2）
一般来讲，差分时差分 格式阶数越高， 得到的波场值精 度越高，但稳定性 会受影响，差分 格式的稳定 性不仅与差分格式本身有 关，也随着差分 精度的不同而不同 ，而且与网格步 长之比的大小有 关，值得指出
u1[i][j]=0;u2[i][j]=0;u3[i][j]=0; //对速度，u1，u2，u3 赋初值
}
for(k=0;k<KMAX;k++)
{
w[k]=exp(-pow(2*PI*FM/R,2)*pow(k*DT,2))*cos(2*PI*FM*k*DT);
}
//波场值计算
for(i=0;i<I;i++)
������������ 为 Nyquist 频率，一般取为主频的两倍，G 为每个波长所占的网格点数，时间、空间为两阶差分的 情况 G 取 8，而时间、空间为四阶差分的情况 G 取 4。
则由上面两式得 v∈[1000，3061.8]
由上面定义的数值，用 Matlab 将震源函数 s(t) 成图，下面是成图的源程序：
的地震波的振幅愈来愈小。
以下为小模型波前快照的源程序： /********************************** *声波方程正演模拟——小模型波前快照* ***********************************/ #include <stdio.h> #include <math.h> #include <stdlib.h> #define PI 3.14159265 #define FM 20 #define R 4 #define KMAX 1000 #define DT 0.002 #define I 200 #define J 200 #define SI 100 #define SJ 100 #define DH 10
第一时刻为地震波还未传播到边界上的某时刻，第二时刻为地震波已经传播到边
界上的某时刻，体会其人工边界反射；对于大模型，定义为水平层状速度模型（至
少两层）；做两个实验，一是将震源点放在区域表层任一点，记录下某些时刻的
波前快照，体会地震波在两种介质的分界面上传播规律；二是合成一个地震记录，
即记录下与震源同一深度点的各点所有时刻的波场值，并指出记录上的同向轴分
m[i]=(double*)malloc(y*sizeof(double));
return m;
}
void exchange(double **u1,double **u2,double **u3) //循环交换赋值 {
int i,j; for(i=0;i<I;i++)
for(j=0;j<J;j++) u1[i][j]=u2[i][j];
3
vmax * t / h
的是差分格式的稳定性与微分方程无关，本实验所采用的（2）式为条件稳定，稳定条件：
8，
（2）式由差分方程近似替代微分方程所得，所以如果空间和时间采样间隔不当，就会产生频散现象，导致
波形畸变，甚至派生出多个同相轴，为减少频散，Δh 需满足 Dablain 的经验公式： h vmin /(GfN )
声波方程有限差分法数值模拟
学 院：海洋地球科学学院 专业：地球信息科学与技术 年级：2012 级 姓 名：王昊 指 导教师：宋鹏
声波方程有限差分法数值模拟
2012 级 地球信息科学与技术 12040032037 王昊
摘要：本实验应用声波方程作为正演模拟的波动方程，将所提供震源函数离散
后绘图，并将给定两个二维速度-深度模型（一个小模型；一个大模型），绘出图 形。通过模拟地震波在介质中的传播，理解实际勘探中地震波在地层中的传播规 律，在模拟水平层状速度模型中，体会地震波在两种介质分界面的传播规律，并 能够从地震记录中识别出反射波，透射波，多次波，折射波和绕射波。并通过模 拟人工合成的地震记录，体会地震勘探基本原理和方法，验证地震波传播能量波 形变化趋势。
在小模型中，可将介质看做是单一介质，即在整个介质中，波的传播速度是相同的，在 这里设地震波速度为 2500m/s，以下是 Matlab 成图的源程序：
clf for x=1:200
f(x)=2500; end plot(f);axis([1,200,1,3000]);title('小模型的速度-深度模型'); figure; x=[0 0 200 200]; y=[0 200 200 0]; patch(x,y,'g'); title('小模型'); axis ij; 用 matlab 编译后得图像，见下图：