地球物理中的有限单元法解剖.docx

上传人:b****1 文档编号:1326810 上传时间:2022-10-20 格式:DOCX 页数:18 大小:231.33KB
下载 相关 举报
地球物理中的有限单元法解剖.docx_第1页
第1页 / 共18页
地球物理中的有限单元法解剖.docx_第2页
第2页 / 共18页
地球物理中的有限单元法解剖.docx_第3页
第3页 / 共18页
地球物理中的有限单元法解剖.docx_第4页
第4页 / 共18页
地球物理中的有限单元法解剖.docx_第5页
第5页 / 共18页
点击查看更多>>
下载资源
资源描述

地球物理中的有限单元法解剖.docx

《地球物理中的有限单元法解剖.docx》由会员分享,可在线阅读,更多相关《地球物理中的有限单元法解剖.docx(18页珍藏版)》请在冰豆网上搜索。

地球物理中的有限单元法解剖.docx

地球物理中的有限单元法解剖

地球物理算法技术(论文)

 

地球物理中的有限单元法

 

院系:

地球物理与信息技术院

姓名:

刘雅宁

学号:

2010120053

任课老师:

张贵宾

地球物理中的有限单元法

一、有限单元法的介绍

在地球物理理论计算中,存在着两类基本问题:

正问题和反问题。

给定场源的分布,求解场值的大小,这是正问题,或者称为正演问题。

地球物理正演的数值计算方法,种类很多,最常用的有:

有限差分法和有限单元法。

有限单元法是50年代首先在弹性力学中发展起来的方法。

主要优点是,适用于物性参数复杂分布的区域,但计算量大。

随着计算机技术的发展,有限单元法在解决各个工程领域的许多数学物理问题中,得到了广泛的应用,称为一种高效、通用的计算方法。

地球物理中的一些边值问题,也采用了有限单元法,解决了许多从前无法计算的地球物理问题。

有限单元法解决数学物理边值问题的基本思路和过程如下:

1、给出地球物理边值问题中的偏微分方程和边界条件(及初始条件)。

这一点看起来似乎容易,但做起来并不容易,特别是边界条件的给定。

只有对地球物理方法的原理和问题有深入的理解,才能给边值问题中的偏微分方程和边界条件以正确的描述。

2、将地球物理边值问题转变为有限元方程。

实现这种转变的主要数学工具是变分法,用变分法得到的有限元法方程称为泛函极值问题。

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

5

5

5

5

8

8

8

8

11

11

11

11

14

14

14

14

17

17

17

j

1

2

3

6

4

5

6

9

7

8

9

12

10

11

12

15

13

14

15

m

4

1

2

3

7

4

5

6

10

7

8

9

13

10

11

12

16

13

14

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

17

20

20

20

20

23

23

23

23

26

26

26

26

29

29

29

29

32

32

32

32

18

16

17

18

21

19

20

21

24

22

23

24

27

25

26

27

30

28

29

30

33

15

19

16

17

18

22

19

20

21

25

22

23

24

28

25

26

27

31

28

29

30

(4)节点坐标

表2

节点

1

2

3

4

5

6

7

8

9

10

11

x

0

0

0

40

40

40

80

80

80

120

120

y

30

165

300

44.3

4712

172.1

7355

300

49.9

9147

174.9

9573

300

43.5

0927

171.7

5464

12

13

14

15

16

17

18

19

20

21

22

120

160

160

160

200

200

200

240

240

240

280

300

28.8

3252

164.4

1626

300

14.8

6395

157.4

3198

300

10.0

7671

155.03836

300

17.3

7467

23

24

25

26

27

28

29

30

31

32

33

280

280

320

320

320

360

360

360

400

400

400

158.6

8733

300

32.3

3098

166.1

6550

300

45.8

7336

172.9

3668

300

49.7

8717

174.8

9359

300

(5)第一类边界节点数ND1=24;

(6)第一类边界节点号和场值:

表3

边界节点序号

1

2

3

4

6

7

9

10

12

13

15

边界节点场值

(*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

16

18

19

21

22

24

25

27

28

30

31

32

33

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点)。

运行程序后,得计算结果如下:

内部节点号

数值解

解析解

误差

5

0.0028377837

0.0024170650

0.0004207187

8

0.0036490047

0.0028453949

0.0008036098

11

0.0045046713

0.0033407882

0.0011638831

14

0.0054020132

0.0038639188

0.0015380944

17

0.0061408007

0.0042171525

0.0019236482

20

0.0060967710

0.0041428832

0.0019538878

23

0.0049695643

0.0036416100

0.0013279542

26

0.0037122145

0.0029888200

0.0007233946

29

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

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置:首页 > 幼儿教育 > 幼儿读物

copyright@ 2008-2022 冰豆网网站版权所有

经营许可证编号:鄂ICP备2022015515号-1