用径向基函数求微幅波的势函数
一、 径向基函数
1.径向基函数介绍
径向基函数是无网格方法的一种，我们知道有限元计算中网格畸变会带来困难，相对与有限元，无网格方法不受网格的影响，取点具有随意性，能更准确地获得获得更加复杂系统的近似解，对于流体，以拉格朗日建立的光滑粒子动力学方法（SPH ）现阶段解决复杂的问题得到广泛的应用。
径向基函数是一类以点x 到xi 的距离di=||x-xi||为自变量的函数。
它具有形式简单、空间维数无关、各向同性的优点，由于径向基函数有许多表达式，在这里采用吴宗敏（1995）提出的正定紧支径向基函数，其表达式为
())5307282366(1)(54326
r r r r r r x I +++++-=+φ （1）
()⎩⎨⎧≤≤-=-+
其他时当 010116
r r r mI I d d r =，mI d 是定义在节点I x 处的径向基函数的支撑域半径，I I x x d -=是以点x 到节点I x 的距离；
2.对于径向基函数的插值
设一个函数为
()a x x u T )(φ= （2）
[]T N a a a a ,,,21 = （3）
()[]T N x x x x )(,,),()(21φφφφ = （4） 式（4）有N 个未知数，令近似函数u(x)在节点I x 处的值等于函数u(x)在该节点处的值I u ，即I I u x u =)(，可得到N 个线性方程组：
u Aa = （5）
式中 ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=)()()()()()()
()()()(....)()(21222121211121N N N N N N N T T T x x x x x x x x x x x x A φφφφφφφφφφφφ （6） []T N u u u u ,,,21 = （7）
这里u 一般是已知的；’
由式（5）解出系数矩阵（3）式a ，代入式（2）中，得
u x N u A x x u T )()()(1==-φ
式中形函数矩阵N 为
1)()(-=A x x N T φ （8）
二、 微幅波的势函数
微幅波势函数的表达式：
)sin(cosh )(cosh t kx kh
h z k Ag σσφ-+= （9） 其中kh gk tanh 2=σ （10）
kh gT L tanh 2π
=
（11） L
k π2= （12） tanh(x)和cosh （x ）分别为双曲正切和双曲余弦函数。
三、 计算步骤及过程
1.已知条件
设周期T=6s ，波高A=1.06m ，水深h=4m ，带入（4）可以得到波长L=34.75m ，如图
2.利用势函数
在t=0的条件下势函数为kx kh
h z k Ag sin cosh )(cosh +=
σφ （13） 在势函数的范围内随意取五个点内，分别为（4.344，-3.5），（13.031 ，-3.3），（17.375 ，-3.2），（21.719 ，-3.1），（30.406 ，-3）。
3.利用径向基函数
由（2）式是理论解，这里用径向基函数得到的一个数值解来模拟理论解，同时证明径向基函数的精确性。
利用I I x x d -=得到五个点及自身的距离，mI d 是定义在节点I x 处的径向基
函数的支撑域半径，在这里取mI d =L=34.5m ，得到mI
I d d r =，代入吴宗敏公式（1）
得)(j i
x φ再将)(j i x φ代入（6）式得到55⨯的矩阵，其值为 ⎥⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=63.80 3.80 6 2.115.350.863.800.030.872.115.3565.352.110.873.805.3563.800.0290.872.113.806)()()()()()()()()(552515522212512
111x x x x x x x x x A φφφφφφφφφ 将A 求逆得1-A
⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=- 0.63 1.72- 1.72- 7.70 2.22 11.57- 1.07- 6.37 0.14 1.07- 2.22 11.57-
19.26 11.57- 2.22 1.07- 6.37 11.57- 7.70 1.72- 0.14
1.07-
2.22 1.72- 0.631A 另外将五个点代入（13）式中可以得到五个势函数的值i u ，代入式（7）得到
[]T u 5.60- 5.58- 0.00 5.55 5.53=
用（5）式u Aa =求逆u A a 1-=，得到
[] 0.93 3.87- 0.19 3.65 0.90-=a 。
四、 求微幅波势函数的数值解 根据mI
I d d r =可以 设L z x r /))5.3(() 4.344(221--+-=
L z x r /))3.3(() 13.031(222--+-=
L z x r /))2.3(() 17.375(223--+-=
L z x r /))1.3(() 21.719(224--+-=
L z x r /))3(() 30.406(225--+-=
得到数值解函数的表达式
()∑===n
i i i T
r a a y x y x u 1)(),(,φφ ()=y x u ,)1)(5307282366()1(5
43265
1i i i i i i i i r r r r r r r a -+++++-∑= 五、 精度验算
在波长L=34.75m 范围内取点与真实值比较有
值解是具有可行的。
六、 计算程序
在在这里使用MATLAB 来进行计算
程序如下：
l=34.75;
x=[ l/8
3*l/8
l/2
l*5/8
l*7/8];
y=[ -3.5
-3.3
-3.2
-3.1
-3];
A=1.06,h=4,g=9.8,t=0;
k=2*pi/l;
zi=(g*k*tanh(k*h))^0.5;
% zi 为σ
for i=1:5
for j=1:5
fai(i,j)=((x(i)-x(j))^2+(y(i)-y(j))^2)^0.5/l;
faix(i,j)=(1-fai(i,j))^6*(6+36*fai(i,j)+82*fai(i,j)^2+72*fai(i,j)^3 ...
+30*fai(i,j)^4+5*fai(i,j)^5);
end
end
faix=inv(faix);
%faix 为求逆矩阵
for i=1:5;
ceita(i)=A*g/zi*cosh(k*(y(i)+h))/cosh(k*h)*sin(k*x(i)-zi*t);
end
ceita=ceita';
%a 为待定系数
a=faix*ceita;
for i=1:5;
%faii 为所求的径向基函数
r(i)=((x(i)-10*l/16)^2+(y(i)+3.4)^2)^0.5/l
end
faii=a(1)*(1-r(1))^6*(6+36*r(1)+82*r(1)^2+72*r(1)^3 ...
+30*r(1)^4+5*r(1)^5)+a(2)*(1-r(2))^6*(6+36*r(2)+82*r(2)^2 ...
+72*r(2)^3 +30*r(2)^4+5*r(2)^5)+a(3)*(1-r(3))^6*(6+36*r(3) ...
+82*r(3)^2+72*r(3)^3 +30*r(3)^4+5*r(3)^5)+a(4)*(1-r(4))^6*(6+36*r(4) ...
+82*r(4)^2+72*r(4)^3+30*r(4)^4+5*r(4)^5)+a(5)*(1-r(5))^6*(6+36*r(5) ...
+82*r(5)^2+72*r(5)^3+30*r(5)^4+5*r(5)^5);
七、 对于求时间的势函数
设)(),(),,(t y x t y x ςφφ=
此方程为强解的形式，优点是计算方便，不复杂，但精确度差，特别是计算高阶的函数。
我们上面所做的是关于t=0时得到的一个函数),(y x φ，可以设
N b b b t ,,,21 =分别得到),(,),,(),,(),(21y x y x y x y x N φφφφ =
利用径向基函数，可以得到函数的表达式
),()(),()(),()(),,(2211y x t N y x t N y x t N t y x N N
φφφφ+++= )(t N i 为形函数。
参考文献
张雄，刘岩，无网格法<M>，清华大学出版社，2005.。