地球物理中的有限单元法分解Word格式文档下载.docx
《地球物理中的有限单元法分解Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《地球物理中的有限单元法分解Word格式文档下载.docx（21页珍藏版）》请在冰点文库上搜索。
是单元的u值列向量；
是单元系数矩阵，
，s，t=i，j，m，
是对称矩阵。
4总体合成:
将单元的场值列向量ue扩展成全体节点的场值向量u，u=（u1u2…ui…uj…um…und），按照节点的总体序号，将单元系数矩阵中的各元放在
的相应行与列的交叉位置上，其余位置的元为零，这样单元积分可写成
各单元积分相加时，只要将ke相加即可；
5求变分：
通过以上四个步骤，已经将连续函数g的泛函，离散成各节点g值的多元函数：
泛函的极值等于多元函数的极值，用多元函数求极值的方法，对上式求微分，
因为K是对称矩阵，有
，所以
，由于
，所以由上式得
。
得到含有ND个元的ND个方程联立的线性代数方程组；
6解线性代数方程组：
首先将第一类边界条件代入，通过定带宽储存的对称带型线性方程组，解方程组，得到各节点的u值，至此有限单元法的求解过程结束。
三、实现过程
为了验证有限单元法的有效性，我们设计一个规则形状的地下矿体，给出模型：
1、模型
密度均匀的水平半无限空间，一个均匀球体S，球体半径R=10m，剩余密度σ=1g/cm3，球心坐标（a，b）=（200，-100）。
对于均匀球体来说，它与将其全部剩余质量集中在球心处的点的质量所引起的异常完全一样。
若球心的埋藏深度为D，球的半径为R，剩余密度为σ，则它的剩余质量为M=（4πR3σ）/3,通过原点的任意水平剖面上则重力异常的解析解表达式为：
设测线长400m，高程变化（-200,300），地形设为一曲线：
y=20*sin（0.02*x）+30，其中x为测线离原点的水平距离，y为高程，则S引起的重力异常为：
2、剖分
通过Matlab建模，得到地形曲线图，再用三角单元对划出的区域进行剖分（程序附后），剖分后对单元和节点进行编号，并将节点的xy坐标和单元节点号列表（表1和表3），分别放在XY（2，ND）和I3（3，NE）两个二维数组中。
剖分图如下：
图1：
剖分图
根据重力异常的解析解表达式算出区域边界上及区域内部的节点值。
编制计算程序，根据边界上节点的场值，算出区域内部节点的场值，将计算出来的内部节点场值与解析解比较，在有效数字内，若计算值与理论值一致或相差不大，则验证了有限单元法的可效性。
3、剖分后的节点分布及单元编号见下表：
（1）节点总数ND=33；
（2）单元总数NE=40；
（3）单元节点编号
表1
单元
序号
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
i
j
m
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
（4）节点坐标
表2
节点
x
80
120
y
165
300
44.3
4712
172.1
7355
49.9
9147
174.9
9573
43.5
0927
171.7
5464
160
200
240
280
28.8
3252
164.4
1626
14.8
6395
157.4
3198
10.0
7671
155.03836
17.3
7467
320
360
400
158.6
8733
32.3
3098
166.1
6550
45.8
7336
172.9
3668
49.7
8717
174.8
9359
（5）第一类边界节点数ND1=24；
（6）第一类边界节点号和场值：
表3
边界节点序号
边界节点场值
（*103）
2.676
8204
2.023
8112
1.249
8540
4.031
5208
1.398
0970
5.914
4786
1.534
9158
9.042
7687
1.646
9267
14.66
69839
1.720
8469
21.182479
1.7467240
19.149463
1.7208469
11.445633
1.6469267
6.4876173
1.5349158
4.0165414
1.3980970
2.6832691
1.9555145
1.2498540
根据书中给出的第一类边界条件、三角元剖分、线性插值的位场延拓得有限单元程序框图：
编制计算程序，并将已知的边界节点数据输入，来计算区域内部的节点值（5、8、11、14、17、20、23、26、29点）。
运行程序后，得计算结果如下：
内部节点号
数值解
解析解
误差
0.0028377837
0.0024170650
0.0004207187
0.0036490047
0.0028453949
0.0008036098
0.0045046713
0.0033407882
0.0011638831
0.0054020132
0.0038639188
0.0015380944
0.0061408007
0.0042171525
0.0019236482
0.0060967710
0.0041428832
0.0019538878
0.0049695643
0.0036416100
0.0013279542
0.0037122145
0.0029888200
0.0007233946
0.0027385908
0.0024087478
0.0003298430
小结
通过计算出的解与用解析公式得到的解进行比较，误差相差不大，证明了有限元方法是一种有效的正演方法。
附录：
计算机程序
地形及图形剖分程序：
clc,clear
x=0:
400;
y=-200:
0;
[XY]=meshgrid（x,y）;
y1=20*sin（0.02*x）+30;
%地形%
Model=50*exp（-（200-X）.^2/120-（-100-Y）.^2/120）;
%模型%
x1=40*x（1:
11）;
y2=20*sin（0.02*x1）+30;
%观测点%
%绘图
figure
plot（x1,y2,'
Rv'
'
MarkerFaceColor'
r'
MarkerSize'
8）
legend（'
剖分图'
）
holdon
area（x,y1,'
linestyle'
none'
contourf（X,Y,Model,400,'
set（gca,'
ticklength'
[0.00.0],'
FontName'
Verdana'
FontWeight'
Bold'
fontsize'
12）%
colormapjet（32）
forj=1:
length（x1）
fori=1:
X1（i,j）=x1（j）;
end
Y1（1,j）=300;
Y1（2,j）=300/2+y2（j）/3;
Y1（3,j）=y2（j）;
end
plot（40*x（1:
11）,y2,'
fori=1:
holdon
plot（x1,Y1（i,:
）,'
ko'
x1,Y1（i,:
k'
linewidth'
1.2）;
plot（zeros（3,1）+x1（j）,Y1（:
j）,'
zeros（3,1）+x1（j）,Y1（:
length（x1）-2
plot（x1（i:
i+1）,[Y1（1,i）,Y1（2,i+1）],'
i+1）,[Y1（3,i）,Y1（2,i+1）],'
plot（x1（length（x1）-1:
length（x1））,[Y1（1,length（x1）-1）,Y1（2,length（x1））],'
length（x1））,[Y1（3,length（x1）-1）,Y1（2,length（x1））],'
k=0;
fori=3:
-1:
k=k+1;
text（X1（i,j）+2,Y1（i,j）+10,num2str（k）,'
color'
10,'
fontname'
华文中宋'
）;
axis（[0400-200400]）
节点坐标值计算及比较的Fortran程序：
PROGRAMsecond
DIMENSIONX（1:
3）,Y（1:
3）,B（1:
3）,C（1:
3）
DIMENSIONXX（1:
3,1:
11）,YY（1:
11）,XY（1:
2,1:
33）
DIMENSIONDG（1:
11）
DIMENSIONUO（40）,NDB（9）
INTEGER:
:
ND,NE,ND1,IE
REAL,ALLOCATABLE:
KE（:
:
）,SK（:
）,A（:
）,I3（:
P（:
）,NB1（:
）,U1（:
）,U（:
ND=33
NE=40
ND1=24
ALLOCATE（U（1:
ND）,P（1:
ND）,NB1（1:
ND1）,U1（1:
ND1）,I3（1:
NE））
!
地质体的模型参数
PI=3.14159!
PI常量
G=6.672E-11!
万有引力常量
X1=200.;
Y1=-100.!
球心S坐标（x1,y2）
R=10.!
半径
P1=1.!
密度
DOI=1,11
XX（1,I）=40.*（I-1）!
观测点坐标xx
YY（1,I）=20.*sin（0.02*XX（1,I））+30.!
观测点坐标yy
XX（2,I）=XX（1,I）!
网格点坐标xx1
YY（2,I）=（300.+YY（1,I））/2.!
网格点坐标yy1
XX（3,I）=XX（1,I）!
网格点坐标xx2
YY（3,I）=300.!
网格点坐标yy2
ENDDO
计算异常场值
DOI=1,3
DOJ=1,11
异常体S异常值
DG（I,J）=4.*PI*R**3.*P1*G*（YY（I,J）-Y1）
&
/（3.*（（XX（I,J）-X1）**2.+（YY（I,J）-Y1）**2.）**（3.0/2.0））*1.0E9
WRITE（*,*）DG（I,j）
ENDDO
INPUTDATA
K=3*（I-1）+1
XY（1,K）=XX（1,I）
XY（1,K+1）=XX（2,I）
XY（1,K+2）=XX（3,I）
XY（2,K）=YY（1,I）
XY（2,K+1）=YY（2,I）
XY（2,K+2）=YY（3,I）
U（K）=DG（1,I）
UO（K）=DG（1,I）
U（K+1）=DG（2,I）
UO（K+1）=DG（2,I）
U（K+2）=DG（3,I）
UO（K+2）=DG（3,I）
输入参数I3，存放单元节点编号
OPEN（2,FILE="
I3.txt"
READ（2,*）
DOI=1,NE
READ（2,*）I3（1,I）,I3（2,I）,I3（3,I）
CLOSE
（2）
WRITE（*,*）"
SUCCESS:
InputI3.txt"
J=0
NDB=（/5,8,11,14,17,20,23,26,29/）
DO111I=1,ND
DOK=1,9
IF（I.EQ.NDB（K））GOTO111
J=J+1
NB1（J）=I
U1（J）=U（I）
111CONTINUE
CALLFUNCTION函数
CallMBW
CALLMBW（NE,I3,IW）
CallMBW"
ALLOCATE（KE（1:
3）,SK（1:
ND,1:
IW）,A（1:
IW））
CallUK1
CALLUK1（ND,NE,IW,I3,XY,SK）
WRITE（*,*）'
Success:
CallUK1!
'
CallUB1
CALLUB1（ND1,NB1,U1,ND,IW,SK,U）
CallUB1!
CallLDLT
CALLLDLT（SK,ND,IW,U,IE）
CallLDLT!
NB1U1
OPEN（3,FILE='
NB1U1.txt'
WRITE（3,"
（2A15）"
）"
NB1"
"
U1"
DOI=1,ND1
WRITE（3,"
（I15,f15.5）"
）NB1（I）,U1（I）
CLOSE（3）
WRITE（*,*）'
OutputNB1U1.txt!
XY
OPEN（3,FILE="
XY.txt"
X"
Y"
DOI=1,ND
（2F15.5）"
）XY（1,I）,XY（2,I）
OutputXY.txt"
Result
OPEN（4,FILE='
Result.txt'
WRITE（4,"
（4A15）"
节点号"
数值解"
解析解"
误差"
DOI=1,ND
（I15,f15.10,f15.10,f15.10）"
）I,U（I）,
UO（I）,U（I）-UO（I）
CLOSE（4）
OutputResult.txt!
DEALLOCATE
DEALLOCATE（I3,KE,A,U,P,NB1,U1）
ENDPROGRAMsecond
SUBROUTINE
MBW计算总体系数矩阵的半带宽
SUBROUTINEMBW（NE,I3,IW）
INTEGERM
REALI3（1:
NE）
IW=0
M=MAX（ABS（I3（1,I）-I3（2,I））,ABS（I3（2,I）-I3（3,I））,
ABS（I3（3,I）-I3（1,I）））
IF（（M+1）.GT.IW）IW=M+1
END
UK1集成总体矩阵
SUBROUTINEUK1（ND,NE,IW,I3,XY,SK）
NE）,XY（1:
ND）,SK（1:
IW）
REALX（1:
REALKE（1:
DOJ=1,IW
SK（I,J）=0
DOL=1,NE
DOJ=1,3
I=I3（J,L）
X（J）=XY（1,I）
Y（J）=XY（2,I）
CALLUKE1（X,Y,KE）
DO1J=1,3
NJ=I3（J,L）
DO1K=1,J
NK=I3（K,L）
IF（NJ.LT.NK）GOTO11
NK=（NK-NJ+IW）
SK（NJ,NK）=SK（NJ,NK）+KE（J,K）
GOTO1
11NJ=（NJ-NK+IW）
SK（NK,NJ）=SK（NK,NJ）+KE（J,K）
NJ=NJ+NK-IW
1CONTINUE
RETURN
END
UKE1计算单元系数矩阵KE
SUBROUTINEUKE1（X,Y,KE）
C
（1）=Y
（2）-Y（3）
C
（2）=Y（3）-Y
（1）
C（3）=Y
（1）-Y
（2）
B
（1）=X（3）-X
（2）
B
（2）=X
（1）-X（3）
B（3）=X
（2）-X
（1）
S=2.*（C
（1）*B
（2）-C
（2）*B
（1））
DOJ=1,I
KE（I,J）=（C（I）*C（J）+B（I）*B（J））/S
UB1加上第一类边界条件
SUBROUTINEUB1（ND1,NB1,U1,ND,IW,SK,U）
REALNB1（1:
ND1）,SK（1:
IW）,U（1:
ND）
U（I）=0.
J=NB1（I）
SK（J,IW）=SK（J,IW
升级会员 