地震波动方程(教学知识)

第三章 地震波动方程
现在，我们用前一章提出的应力和应变理论来建立和解在均匀全空间里弹性波传播的地震波动方程。

这章涉及矢量运算和复数，附录2对一些数学问题进行了复习。

3.1 运动方程（Equation of Motion ）
前一章考虑了在静力平衡和不随时间变化情况下的应力、应变和位移场。

然而，因为地震波动是速度和加速度随时间变化的现象，因此，我们必须考虑动力学效应，为此，我们把牛顿定律（ma F =）用于连续介质。

3.1.1一维空间之振动方程式
质点面上由于应力差的存在而使质点产生振动。

如图1-3所示，考虑一薄棒向x 轴延伸，其位移量为u ：
Fig3-1
则其作用力为“应力”X “其所在的质点面积”，所以其两边的作用力差为
()()()dxds x
x dx x ds ∂∂=
-+σ
σσ 惯量﹙inertia ﹚为
22t
u dxds ∂∂ρ
所以得出
x
t u ∂∂=∂∂σρ22 ……………………………………………………... （3-1）
其中ρ为密度﹙density ﹚，σ为应力﹙stress ﹚=x
u
E ∂∂。

3-1式表示，物体因介质中的应力梯度﹙stress gradient ﹚而得到加速度。

如果ρ与E
为常数，则3-1式可写为
2
22221t u
c x u ∂∂=∂∂ …………………………………………………… （3-2）
其中ρ
E
c =
运用分离变量法求解（3-2）式，设u=F(x)T(t),(3-2)式可以变为
T X c T X ''=
''2
1
设22ω-='
'=''T
T X X c 则可得：c
x i
t
i e
X e
T ωω±±∝∝,
考虑欧拉公式：)sin()cos(),sin()cos(t i t e t i t e
t i t
i ωωωωωω-=+=-
()
()
()
()
ct x c
i
ct x c
i
ct x c
i
ct x c
i
De
Ce
Be
Ae
u ---+-++++=ωωωω （3-3）
其中A,B,C,D 为根据初始条件和边界条件确定的常数。

考虑到ω可正可负，方程式的解具有()()ct x g ct x f u ++-=的形式，其中f 及g 为波的函数，以c 的波行速度向+x 与-x 方向传递。

我们可以采用如下程序模拟地震波的传播。

平面波在均匀介质里沿x 方向传播，
剪切波的齐次微分方程可表达为：
2
2
222x u t u ∂∂=∂∂β 这里u 是位移。

对100公里的波长和假定秒公里/4=β的情况，我们写出用有限差分法解这方程的计算机程序。

用长度间距公里1=dx ，时间间距1.0=dt 秒。

假
定在u （50公里）震源时间函数的形式为：
()()5sin 250t t u π= 0＜t ＜5秒
用u （0公里）的应力自由边界条件和u （100公里）的固定边界条件。

用有限差分图解来近似二次导数：
2
1
1222dx
u u u x u i i i -++-=∂∂ 以4秒的间隔画出1-33秒的图。

M = moviein(101);
dx=1;dt=0.1;tlen=3;beta=4; %初始化变量,tlen为震源持续时间，beta为波传播的速度
u1=zeros(101,1);u2=u1;u3=u1;
%u1为前一个时刻的各点的位移，u2为当前时刻的位移，u3为下一个时刻的位移值，开始均假定为零
t=0;
jj=0;
while (t<=33) %模拟的最长时间为33秒
for ii=2:100
rhs=beta^2*(u2(ii+1)-2*u2(ii)+u2(ii-1))/dx^2; %方程的解
u3(ii)=dt^2*rhs+2*u2(ii)-u1(ii); %对时间求导数
end
%左边为自由边界条件，右边为固定边界条件
u3(1)=u3(2); %左边为自由边界条件
u3(101)=0.0; %右边为固定边界条件
%左右两边为自由边界条件
% u3(1)=u3(2); %左边为自由边界条件
% u3(101)=u3(100); %右边为自由边界条件
%左右两边为固定边界条件
% u3(1)=0.0; %左边为固定边界条件
% u3(101)=0.0; %右边为固定边界条件
if(t<=tlen)
u3(51)=(sin(pi*t/tlen)).^2; %地震震源时间函数
end
for ii=1:101
u1(ii)=u2(ii);u2(ii)=u3(ii); %时刻的更新
end
plot(u2); %绘制目前的波形图
ylim([-1.2 1.2]);
M(:,jj+1) = getframe; %获得当前的图像
t=t+dt; %时间延长
end
movie(M) %演示波形传播
3.1.2三维空间之振动方程式
推导三维空间之振动方程式的过程，与上节中所采用的一维空间讨论方式类似，如图3-2所表示，先探讨在x 方向之位移量u ：
Fig3-2
在y-z 面上的作用力差为：()()[]dydz x dx x xx xx σσ-+ 在x-z 面上的作用力差为：()()[]
dxdz y dy y yx yx σσ-+ 在x-y 面上的作用力差为：()()[]dxdy z dz z zx zx σσ-+ 惯量为：
22t
u
dxdydz ∂∂ρ
得出
x zx yx xx f z
y x t u +∂∂+∂∂+∂∂=∂∂σσσρ22 ………………………………….. ﹙3-4﹚
其中σxx 、σyx 及σzx 分別为stress tensor 在xx ﹙x 面方向、x 力方向﹚，yx ﹙y 面方向、x 力方向﹚及zx ﹙z 面方向、x 力方向﹚方向的分量。

注意，在本讲义中有关stress tensor 的两个下标﹙indexes ﹚之定义，依序为面的方向与力的方向。

将σxx 、σyx 及σzx 与其对应的应变之关系代入3-4式可推导得出三维空间之振动方程式如下：
()x f u x
t u +∇+∂∂+=∂∂222μθ
μλρ …………………………………. ﹙3-5a ﹚
其中λ及μ为常数，而2
∇为Laplacian operator ，代表22
2
222z
y x ∂∂+∂∂+∂∂。

以相同的方法，可以得出在y 及z 方向的振动方程式，若其位移量分別为v 与w ，。