地球物理中的有限单元法解剖Word下载.docx
《地球物理中的有限单元法解剖Word下载.docx》由会员分享,可在线阅读,更多相关《地球物理中的有限单元法解剖Word下载.docx(18页珍藏版)》请在冰豆网上搜索。
![地球物理中的有限单元法解剖Word下载.docx](https://file1.bdocx.com/fileroot1/2022-10/12/d58a3c2c-9341-4e1c-a94c-5f3fc3c8c6ef/d58a3c2c-9341-4e1c-a94c-5f3fc3c8c6ef1.gif)
3、用优先单元法解决泛函极值问题其步骤大致如下:
把研究区域剖分成有限个小单元,在每个单元上,把函数简化成线性函数、二次函数或高次函数,这称为单元上函数的插值。
用简化后的函数计算每个单元上的泛函。
各单元之间,通过单元间节点上的函数值相互联系起来。
对各单元的泛函求和,获得整个区域上的泛函。
这样,有限单元法将连续函数的泛函,离散成各单元节点上函数值得泛函。
根据泛函取极值的条件,得到各节点的函数值应满足的线性代数方程组。
解代数方程组,得到各节点的函数值。
有限单元法的主要优点是,适用于物性复杂分布的地球物理问题,而且,其解题过程也比较规范化。
这些优点是有限单元法在地球物理中获得广泛的应用。
但是,有限单元法是区域性方法,必须在全区域中进行剖分。
剖分后的单元和节点数目多,最后得到的线性代数方程组很大。
特别是三维问题和地球物理中常遇到的无解区域问题,需要中、大型计算机,才能完成有限单元法的计算。
这是有限单元法的主要缺点。
二、基本步骤
1用三角单元对区域进行剖分:
用三角单元对整个区域进行剖分,因为三角单元的边容易拟合地形线的形状。
三角形的顶点称为节点,用节点上的离散场值来近似场值的连续分布。
剖分后对单元和节点进行编号,次序号(i,j,m)按逆时针方向排列。
第一类边界条件的节点上的场值是已知的,其余节点上的场值是待求的;
2线性插值:
在三角单元内假定场值u是线性变化的,则u=ax+by+c,u还可表示为
,其中
是形函数,
3单元分析:
将全区域的积分分解为单元e的积分之和
,
继续推得单元e上的积分
其中
是单元的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)+3