三维有限元法计算过程.docx
- 文档编号:11019744
- 上传时间:2023-02-24
- 格式:DOCX
- 页数:30
- 大小:361.18KB
三维有限元法计算过程.docx
《三维有限元法计算过程.docx》由会员分享,可在线阅读,更多相关《三维有限元法计算过程.docx(30页珍藏版)》请在冰豆网上搜索。
三维有限元法计算过程
三维有限元法计算过程
三维有限元法的计算过程:
1)网格单元剖分;
2)线性插值;
3)单元分析;
4)总体刚度矩阵合成;
5)求解线性方程组等部分组成
、偏微分方程对应泛函的极值问题
矿井稳恒电流场分布示意图
主要任务是分析在给定边界条件下,求解稳定电流场的Laplace
方程或Poisson方程的数值解,即三维椭圆型微分方程的边值问题:
Lu三(▽)+(a)+(▽)=F
excycyczcz
du
二E岳=0cn
(空+◎严0
丄
F=I、(x-Xo、(y-yo)、(z-Zo)
上述微分方程边值问题等价于下面泛函的极小值问题:
--U2-U2-U212
J[U]=q[(龙)(片)(芯)]-fU}dxdydz匚』•门二UdS
、网格剖分
:
?
i:
:
hi
图2-4单巷道地质模型
1、网格单元的类型
2、网格单元剖分原则及其步长选择
因此,网格内的单元剖分应按以下剖分原则
1)、各单元节点(顶点)只能与相邻单元节点(顶点)重合,而
不能成为其它单元内点;
3)、
在场变化剧烈的区域网格剖分单元要密一些,在场变化平缓
的区域单元密度应小。
方向的两侧作对称变步长剖分,距o越远,步长应越大。
常用的变步长方法
有:
Xi勺一:
Xi=(iT)c
以上各式中c为常数,展「1、厶Xi为同一坐标轴上相邻步长值。
以X方向为例,可知,X正方向与负方向对称,只相差一负号。
若令,只要给出距
原点最近节点的坐标.^1,由上式即可求出其它相应的步长・叹。
同理可求得y、z方向上的变步长湖、迄。
3、网格剖分方法
kLtkLt+1
kLt+r
1LZ
kLt+(LX-1)L
kLt+m
iLz+1
kLt+(LX-1)L
z+1
kLt+Lz-1kLt+(m+1)Lz+1kLt+LxLz-1
=(k+1)Lt-1
图2-6平面内节点编号示意图
图2-7长方体单元编号示意图
在程序设计的过程中,将采用四面体单元对研究区域进行分析。
因此对上面的长方体网格单元还需作进一步的处理。
任取一编号为m的长方体单元,显然根据m的值可求出其八个节点的节点编号。
为讨论方便,将长方体单元的各顶点的节点号,将其重新编号为0、1、2、3、4、5、6、7,如图2-8所示。
图2-8长方体单元节点编号示意图
可知长方体单元可划分为六个四面体单元,其方法为用1、2、5、6顶点
组成的平面将长方体分成二个三棱柱体,每个三棱柱体又可分成三个四面体单元,其结果如下:
(0,1,2,4)、(1,2,4,5)、(2,4,5,6)、(1,2,3,5)、(2,3,5,6)、(3,5,6,7)。
其组合规律为:
(i+j,1+i+j(j+1)/2,
2ij(5-j)/2,4ij),i=0,1,j=0,1,2。
在具体的计算过程中,按长方体网格编号逐一计算,在每个长方体单元内再按以上的组合规律依次对六个四面体计算。
这样处理后四面体的总数为
6*L_n*(Ly_1)。
三、线性插值分析
任取一四面体单元如图2-9所示,设其顶点相应的节点序号为(i,j,l,m),并设各节点电位为Uk,k=i,j,I,m。
对应坐标为(Xi,yi,z),(xj,yj,zj),(xi,yi,zi),(xm,ym,zm),当单元足够小时,设四面体内部电导率为常数、其电位呈现线性变化。
即:
U=匕2x比3丫心4z(2-29)
对i,j,l,m四个节点有:
5"j〉2Xr3『i〉4N
』UjI2Xj+°3yji4Zj
U]〉2X「〉3y]r4zl(2-30)
.Um八1「2Xm〉3Ym:
4Zm
整理式(2-29)、(2-30)得:
U二NiUiNjUjN|U|NmUm(2-31)
式中u和Nk都是的(x,y,z)坐标函数,其中Nk的值为:
1
山=6何+5"^+0旧,ci」』(2-32)
可知:
akbsdk(k=i,j,l,m)、D只与节点坐标有关。
为了分析式(2-32)所对应Nk的几何意义,现引入面积坐标和体
积坐标的概念。
如图2-10所示为.:
ijl,(Xj,yj、(x「yj)、(xi,yi)为其顶点
上式值为图2-10中三角形.-:
Pj面积的两倍,所以:
具有以下性质:
Ni(x,y)汨
Nk(Xn”n)=丿
4)与节点i对边平行的直线上各点的N,x,y)值不变。
利用面积坐标的定义,可以方便地确定三角单元中任一点的面积
坐标,且可建立一个面积坐标网,如图2-11所示。
因为面积坐标与通常使用的直角坐标系的选择无关,可利用该坐标网构造任意性线、非线性的插值函数。
图2-11三角形面积坐标网
采用与上述相似的分析方法,可知Nk(x,y,z)(k",j,l,m)为任一点p(x,y,z)在四面体ijlm内的体积坐标,具有与三角单元面积坐标相似的性质。
pjlm
Ni(x,y,z)=D
Dijlm
(2-36)
Dpilm
Nj(x,y,z)=
Dijim
(2-37)
Dpijm
Ni(x,y,z):
Dijlm
(2-38)
Dpijl
Nm(x,y,z厂
(2-39)
Djim
上式中Djim为其对应节点所构成四面体的体积。
在式(2-31)假设四面体内电位为线性变化,P点的电位值是通过该点的体积坐标将各节点的电位值线性化计算得出。
如果四面体内的电位值为非线性变化,则只要改变相应点体积坐标的系数即可
形如:
U=SNj5ONjUjPN|U|QNmum
(2-40)
式中S、O、P、Q为P点体积坐标对应的系数,由于四面体内电位的非线性变化规律比较复杂,系数值的确定比较困难,因此在以下各章节的讨论中只考虑电位值的线性变化,对于非线性的变化规律将以后讨论。
四、网格单元分析
整个求解区域Q被剖分为一系列小四面体单元后,泛函J(u)也相应地被离散为各单元上的泛函Je(u),故有J(u)八Je(u)。
m为四面体m
单元总个数。
在任一四面体单元上泛函形式为:
对任意节点的微分得:
Je(u)e:
u:
zuu:
zuu:
“:
u「()()()]負CX「ukcxcy^ukcy已zcukcz
aa
-f—^}dxdydz亠门uds
ukeU
其中k二i,j,l,m
令门二"也•亠(卫)—•亠(d)d•亠
(2)]dXdydz
x:
uk:
x:
y:
uk:
y:
z:
uk:
z
则:
旦^.飞.IHf庇kQ
—dxdydziy'
-e
u^ds
:
'Uk
(2-43)
现对各项分别进行讨论:
1、对于:
•:
■■项
由前面的讨论可知:
u八NkUk,
1
Nk=D(akbkXCkydkZ)
所以
:
U°Nu
_:
xfx
:
:
Nj
丄5
-X
.NlNm
"l-Um
:
X:
X
又因为:
-:
Nkbk
同理可得:
.:
uk(Tx)
-:
Nk
=u
:
:
Uk
Nk
.:
u
■:
y
.:
u
.z
「x「Ukb
u(屯)=Dy(bkbiUi+bkbjUj+bkbUl+bkbmUm).x:
Uk:
xD
-:
Uk釣
1
D^CkCjUiCkCjUjCkClUlCkCmUm)
1
W訶d©UidkdjUjdkdlUldkdmUm)
(2-44)
(2-45)
(2-46)
(2-47)
(2-48)
UUU、
_()()
(二)
:
z:
Uk:
z
1
P(bibkUi
bjbkUjbibkUi
bmbkUm)(CiCkUi-CjCkUj
CiCkUi
CmCkUm)(didkUidjdkUjdidkUidmdkUm)]
1
^^[(bibkCiCkdidk)Ui(bjbk-CjCkdjdQUj
(2-49)
(b|bkclckdldk)Ul(bmbkcmckdmdk)Um]
记:
FrWmca96从gmCjSdjdQUj
(blbkClCkdldk)Ul(bmbk
CmCkdmdk)Um
所以:
匚Fkdv=
D2k
36VeFk
(2-50)
111fdxdydz
2、对于爲:
:
Uk项
由前面分析可知:
f=I、:
(x一x0)、;(y一y0)、•(z一z0),U=Nk所以有:
cUk
Eu
iiifdxdydz=I111Nk(x,y,z)、(x-Xo)、(y-y。
)、(z-zo)dxdydz(2-51)
Ve-ukVe
由:
函数的性质可知:
iiif—dxdydz=Nk(Xo,yo,z°)I
Ve:
uk
(2-52)
I供电点在点(xk,yk,zk)
=<
卫供电点不在点(Xk,yk,Zk)
由式(2-25)分析可知,该项只与研究区域i■■的第三类边界条件有关,所以
当四面体完全位于求解区域内部时(即非边界单元),该项值为零。
在边界处则须对该项进行计算。
1)关于值的计算
由电场理论知识可知,对于点源和均匀介质的情况,其电位函有如下的形
式:
U(x,y,z)=C(C为常数)(2-53)
r
贝S:
卫=(2-54)
£nrr
式中,为由点源到计算点的径向矢量r和外法线矢量n之间的夹角。
式(2-54)可改写成如下形式:
-UUcos二=0(2-55)
mr
对比混合边界条件(2-7)的公式可知:
设边界四面体单元嘔如图2-9所示,由节点i,j,l组成的面为边界
面,另一节点m位于研究区域内。
在三角形•冷内,值是坐标(x,y,z)及cost的函数,为简化计算,
当三角形的面积较小时,用其重心处的值来代替整三角形的值。
设三角形重心P的坐标为(xt,yt,zt),贝卩:
式中,nop为由供电点O到三角形重心的径向矢量,
nop=(Xt—x°,yt—y°,z-z。
)n=n”nij
其中nH、nj分别为节点l,i和节点l,j所在直线的径向矢量
nii十-X|,yi-yi,Zi-Z|)nij=(Xj-X|,yj-yi,Zj-zj
可得:
niinij
=[(yi-yi)(Zj-Zi)—(Zi—Zi)(yj—yi)]i+[(Z—zJ(Xj—xj
(2-59)
(2-60)
(2-61)
(2-62)
-(洛-xi)(Zj-z)]j[(人—X)(yj—yj-(%—yJ(Xj-xjk令A=(yi—yi)(Zj—Zi)-(乙-z)(yj-yi)
B=(乙—Zi)(Xj—xi)—(Xi—xJ(Zj-zi)
C=(Xi—Xi)(yj—yi)_(yi_y|)(Xj-X|)
则:
n=(A,B,C)
所以:
cos,A(Xt-Xo)B(yt-y°)C(Zt-z°)
J==一
r[(Xt-Xo)2(y^y。
)2(Zt-Zo)2]IA2B2C2
I!
r^uds
2)关于…'e:
uk值的计算
由式(2-31)计算可知丄匕-Nk,所以:
cUk
UdS=(NiNkUiNjNkUjNlNkUlNmNkUm)dS(2-64)
:
e;Uk-
为方便计算,设四面体单元-e的边界面r所在的平面为XOY平面,I为原点,li为X轴,并设节点i,j,l,m的坐标分别为(Xi,0,0),(Xj,yj,0),(0,0,0)和
(Xm,ym,Zm)。
再设三角形内的一动点P的坐标为(X,y,0)
因为N「Nj・N「Nm=1。
所以在△上N「NjNl=1。
故Ni,Nj,Ni中
只有两个自由量,设为Ni,Nj,则
11NiNjdxdy
A
yj,-
二odyXj
X_x
Xi+jiyj
yj
/XXjy、y
(-Xi
)dx
Wyj
XiYj
12
yjy
c[(Xi
2Xiyj
XiYj€
2412
_Xi2
-y)
yj
2
Xj2Xjy
-(—y)]-一(Xi
XiYj
yj
y)}dy
(2-65)
yjx
NjNjdxdy二0dy勺
△
:
_e
yj
yj
2
y2dXyj
(2-66)
因为M,Nj,N|中只有两个自由量,
而且三角形
△三个节点i,j,l
的选取排列是任意的。
因此:
12
kHt
心e
k=t
.6
JJNkNtdxdy=«A
t",j,l
(2-67)
所以:
uds
-■(NiNkuiNjNkuj
NiNkUiNmNkUm)ds
乂-.(NiNkUiNjNkUjNlNkUlNmNkUm)dxdy
%坐&竺&乞01
-Ui1
61212
&生&经&鱼0
Uj
12612
&担占色&$0
Ul
12126
0000
-Um一
k
k
k
k
i
j(2-68)
l
m
五、系数矩阵合成
由单元分析可知:
其中Tt项只对参于对边界单元的计算
由泛函极值的必要条件凹巴=0可知:
CUk
°Fk-INk+甲k=0
36Ve
(2-70)
即:
FkVIM
36Ve
供电点在节k点k,Zk)=(x°,y°,z。
)供电点不在(节,点,zQ=(x°,y°,z°)
对k=ijl,m,上式可写成如下的矩阵形式[4]:
kekekeke1uJ_INi1
iiijilimii
ieIeIe
IkjjkjlkjmiiW丨=iINj
(2-71)
对kek;muiini
I11lm丨丨1II1I
1称kmm」LUm」JNmJ
系数矩阵即为四面体单元Ve的单元刚度矩阵Ke,显然单元刚度矩阵为对称阵,其阶数等于单元节点数。
在具体求解计算过程中,应将单元系数矩阵合成总的系数矩阵K。
即
K八Ke,其中Ke=CTKeC为单元系数矩阵在总系数矩阵中的形式;Cm
为4n阶的转换矩阵,其作用将单元体系数矩阵变成总系数矩阵,其元素值为6=1,C2j-1,C31=1,C4m=1其余值全为零。
整个求解区域离散为四面体单元,每一个单元体都可求得其对应的系数矩阵,将各单元刚度矩阵叠加即可得总系数矩阵。
设整个求解区域有n个节
点,则K是一个nn阶的系数矩阵。
合成总矩阵后,可得由n个方程
组成的大型线性方程组,其矩阵表达形式为:
系数矩阵K的特点
总系数矩阵K的阶数为nn,其数据量较大,在计算时对计算机的内
性,得到一定程度的解决。
有限元总系数矩阵具有以下特点:
1)对称性总系数矩阵是一个对称矩阵,运行时系数矩阵可以只存贮矩阵的上三角部分或下三角部分,可节省一半的存贮量。
2)稀疏性由于有限元方法采用四面体单元结点的相互关联性,这决定
了总系数矩阵具有大型、稀疏的特征。
实际上,只有当结点i是结点j的相邻
点时kj=0即总系数矩中的某一行中只有几个非零元素,对于一个四面体单元的四个结点,与任一结点对应的整体刚度矩阵中最多只有15个非零元素。
当
空间结点数较多,网格剖分较密,总系数矩阵很大时其矩阵的稀疏性表现得越突出。
另一方面,总系数矩阵的结构(即非零元素的排列方式)依赖于结点的编号规律。
相邻结点的编号相差较大,其系数矩阵的稀疏性越明显。
3)带状分布特性对于有限元数值分析中的总系数矩阵,非零元素分布在
以主对角线为中心的斜带形区域中,由于对称性,存贮系数时,可只取上三角(或下三角)部分,因而每行元素只存贮半个带宽的元素就可以了。
4)正定性有限元数值分析中,方程离散后形成的系数矩阵具有正定性特点,即主对角线元素在所处行与列为最大且为正值。
六、求解线性方程组
有限元法是求解复杂电磁场问题的有效工具之一,大型工程电磁场有限元计算最后都归结为大型稀疏有限元方程的求解。
在求解区域上对泊松方程进行离散并设置各节点参数及其边界条件之后,便得到一个大型的线性方程组。
求解此方程组,就可以得到各节点的电位值,进而计算出不同装置形式下的视电阻率值。
目前求解线性方程组KU二F的方法很多:
1)、Gauss-Sidel迭代法;
2)、超松驰迭代法;
3))高斯消去法(GausS;
4)、乔列斯基分解法(LDLt)等。
这些方法在求解线性方程组方面已经取得了较为广泛的应用,并且具有计
算精度较高、稳定性强等优点。
本文在程序设计的过程中,主要采用了常规的Gauss消去法和适合大型稀疏矩阵方程的LDLt分解法,并对其计算精度进行对比。
1、Gauss消去法
Gauss消去法是计算机上常用的求解线性方程组的有效算法,此方法分为
消元过程和回代过程。
消元过程是将原方程组化为上三角形方程组的过程,而
回代过程是求解上三角形方程的过程
对于线性方程组:
(3-1)
其中:
K=(kj)nnU珂比…,
F=(f1,f2,…
计(
为方便起见,
3-1)
式为K
(1)U:
二F⑴
(1)
11
(1)■k12
(1)
21
(1)------k22
k2n
即
_k
k
KU=F
为对称正定型稀疏矩阵
••…,Un)T
…,fn)T
■uj
U2
k^n
第一次消元:
设k
(1)=0,至U第i(i=2,3,
记li^k(])/ki11),(i
■",n
=2,3;
)个方程上,得
n),k
(2)u=F
(2),
k⑴
k12
k⑵
k22
其中:
k(『-⑺谓
i⑵=fi
(1)-li1
第二次消元:
「f1⑴1
f2⑴
冗1)
(3-2)
以-1门乘式(3-2)中的第一个方程加
即:
VuJ
U2
9
f
(2)
12
9
a
lUn_
1】
a
f,2l
k⑴
k1n
k⑵
k2n
kn2)
(i,j=2,3,,n)
(i=2,3,,n)
设k;?
式0,记心=ki
(2)/k2?
,到第i(i=3,4,……,n)个方程上,
(i=3,4,得K(3)U
(3-3)
n),以-li2乘式(3-3)中的第二个方程加
_k
(1)[
(1)
11k12
k⑵
k22
I
(1)
k13
k⑵
k23
k⑶
k33
(3-4)
第k-1次消元完成后,得与式
(3-1)
等价的方程组,K(k)U二F(k),即:
■k
(1)
(1)
11k12
k⑵
k22
kkkk
knk)
(3-5)
其中:
ki(k)=k(j7
(kJ)
~-likjkkJj
(i,j
=k,k1,,n)
(k)_f(Z
i_fi
(i=k,k1,,n)
如此进行下去,形方程组,即:
共经过
n-1次消元计算后得到与式(3-1)
等价的上三角
I
(1)I
(1)
k11k12
k⑵
k22
I(i)
ki3
k⑵
k23
(3-6)
然后利用回代公式求解线性方程组的解为:
(n)(n)
n/knn
n
(i)v|(i)(i)
Ui=(fi…kisus)/kii
s兰十
Un
(i=n-1n-2,,2,1)(3-7)
2、LDLT分解法
对有限元方程组KU二F的系数矩阵K采用半带宽压缩存储技术,用二维数组KN(N,M)存放系数矩阵K,其中N为网格节点总数,M为半带宽。
将系数矩阵K矩阵分解为:
K二LDLt(3-8)
式中,L为带状的下三角矩阵,lj=O(i-j・M)
1
D为对角阵dH-
lii
jJ1l
*litljtij-
(i=1,2,3N—Ji)
(3-9)
t注m1tt
式中:
tg0i兰M+1
J—Mi>M+1
(3-10)
于是,解有限元方程KU=F等价于解方程LDLTU二F,Y二LtU,
则得方程组:
LDY=F
'T
LU=Y
J
(3-11)
yfi—(匹)(i—1,2,...,N)(3-12)
t4mltt
为比较以上两种方程组求解方法的计算精度及其稳定性,令式(2-73)中
的U向量值为1,求解右端项F。
在程序计算过程中,利用总系数矩阵K和所
得的右端项F,计算向量U的值,比较U与1的大小关系,即可判断比较方程组求解方法的计算精度及其稳定性。
具体计算结果如表3-2所示。
表3-2计算结果比较
节点数
1000(1QX1QX10)
2541(11X21X11)
求解方法
Gauss
LDLT
Gauss
ldlt
最大值
0.98272
0.9999855
0.97256
1
最小值
-1.368484
0.9999702
-10.25627
0.9999926
平均值
0.9168088
0.99997689
0.84958515
0.9999955
标准差
0.2080572
3.3699E-006
0.40298052
1.5569E-006
经过对上述两种方程组求解算法的计算结果比较可知,在相同条件下
LDLt法计算精度较高,稳定性好,但速度较慢;Gauss消去法计算速度快,
但其计算精度较低,稳定性差。
因此在计算过程中将采用ldlt分解算法。
3、LU分解法
高斯消元的过程相当于下述矩阵乘法运算
M⑵M⑴〔A|b-U|y1
因此由分块矩阵乘法可得
M⑵M⑴A=UM⑵M⑴b=y(3-2-15)
令L=(M
(2)M(f(M■
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 三维 有限元 计算 过程