﻿功能:无限区域中动态弹性波水平振动位移解

格式:
DynamicWaveU(f,t,x,y,Vp,Vs,Density,n)
DynamicWaveU(f,t,x,y,Vp,Vs,Density)
DynamicWaveU(f,t,x,y,Vp,Vs)
DynamicWaveU(f,t,x,y,Vp)
DynamicWaveU(f,t,x,y)

f:为符号变量存储或自定义函数,表示产生垂直力源的函数
t:求解的时间,为数值或矩阵变量
x:观测点水平坐标
y:观测点竖直坐标
Vp:纵波波速,默认4200
Vs:横波波速,默认2700
Density:介质密度,默认2490
n:积分复化数,默认500

原理:

如下图,在无限大弹性介质中,有一垂直力源F作用于原点,求在此作用下点(x0,y0)的水平位移
{<http://img1.ph.126.net/VMjgKZwr1wdRdGX1Ie7C9A==/6598263836286190876.png>}
参考:赵海波,王秀明.一种优化的交错变网格有限差分法及其在井间声波中的应用[J].科学通报,2007,6,52(12):1387-1395

例子:
这里我们假设F为三角形加载,如下图
{<http://img1.ph.126.net/IDgC4-1nxcSw0T3WxT6cfA==/6608261695515360805.png>}
这里算例取
/*
T1=0.001s
T2=0.015s
K1=1000
Vp=4200m/s
Vs=2700m/s
x=32m
y=25m
*/
先将如下力源函数存储在自定义文件里,文件名称为"Ftest.vb"
/*
Public Function Ftest(ByRef t as Double) as Double
if t<0.001
Return 1000*t
ElseIf t<0.015
Return (0.015-t)/0.014
else
Return 0
End If
End Function
*/
//依次执行如下几条命令可以看到
t=0:0.001:0.06;
t=t';
x=DynamicWaveU(Ftest,t,32,25,4200,2700);
x=StandardOne(x)//将数据无量纲化,注意,之前t求转置就是为这一步做准备的,如果不无量化就不需要转置;
plot(t,x)
{<http://img2.ph.126.net/h21A1lWZc6VBKISbnGtkgg==/1394708509702172114.png>}