三维有限元法计算过程
三维有限元法计算过程
三维有限元法的计算过程：
1）网格单元剖分；
2）线性插值；
3）单元分析；
4）总体刚度矩阵合成；
5）求解线性方程组等部分组成。
一、偏微分方程对应泛函的极值问题
矿井稳恒电流场分布示意图
主要任务是分析在给定边界条件下，求解稳定电流场的Laplace 方程或Poisson方程的数值解，即三维椭圆型微分方程的边值问题：
)
()((0)(0
)()()(000z z y y x x I F u n u
n
u F z u z y u y x u x Lu w D ---=????
?????=+??=??=????+????+????≡ΓΓ+Γδδδγσσσ 上述微分方程边值问题等价于下面泛函的极小值问题：
dS U dxdydz fU z U y U x U U J w D ?????Γ+Γ+ΓΩ
+-??+??+??=222221
}])()()[(2{][γσσ
二、网格剖分
∞1
ρi i h ρ.........
...
1、网格单元的类型
图2-5 网格单元类型
2、网格单元剖分原则及其步长选择 因此，网格内的单元剖分应按以下剖分原则
1）、各单元节点（顶点）只能与相邻单元节点（顶点）重合，而
不能成为其它单元内点；
2）、如果求解区域对称，那么单元剖分也应该对称；
3）、在场变化剧烈的区域网格剖分单元要密一些，在场变化平缓
的区域单元密度应小。
4）、网格单元体的大小变化应逐步过渡。
根据上述剖分原则，以x 、y 、z 坐标轴原点o 为中心，分别向x 、y 、z 方向的两侧作对称变步长剖分，距o 越远，步长应越大。常用的变步长方法有：
c i x x i i )1(1+=?-?+ c x x i i =??+/1(i ≠0)
c x x i i =?-?+1
1
1(i ≠0) 以上各式中c 为常数，1+?i x 、i x ?为同一坐标轴上相邻步长值。以x 方向为例，可知，x 正方向与负方向对称，只相差一负号。若令00=?x ，只要给出距原点最近节点的坐标1x ?，由上式即可求出其它相应的步长i x ?。同理可求得y 、
z 方向上的变步长i y ?、i z ?。
3、网格剖分方法
图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）。其组合规律为：（j i +，2/)1(1+++j j i ，2/)5(2j j i -++，j i ++4）
，1,0=i ，2,1,0=j 。 在具体的计算过程中，按长方体网格编号逐一计算，在每个长方体单元内再按以上的组合规律依次对六个四面体计算。这样处理后四面体的总数为
1)(L *L_n *6Y -。
三、线性插值分析
任取一四面体单元如图2-9所示，设其顶点相应的节点序号为),,,(m l j i ，并设各节点电位为k u ，m l j i k ,,,=。对应坐标为),,(i i i z y x ，),,(j j j z y x ，),,(l l l z y x ，
),,(m m m z y x ，当单元足够小时，设四面体内部电导率为常数、其电位呈
现线性变化。即：
z y x U 4321αααα+++= （2-29）
图2-9 四面体单元示意图
对m l j i ,,,四个节点有：
??????
?+++=+++=+++=+++=m
m m
m l l l l j j j j i i i
i z y x u z y x u z y x u z y x u 4321432143214321αααααααααααααααα （2-30）
整理式(2-29)、（2-30）得：
m m l l j j i i u N u N u N u N U +++= （2-31）
式中u 和k N 都是的),,(z y x 坐标函数，其中k N 的值为：
)(1
z d y c x b a D
N k k k k k +++=
，m l j i k ,,,= （2-32） 可知：k k k k d c b a ,,,（m l j i k ,,,=）、D 只与节点坐标有关。
为了分析式（2-32）所对应k N 的几何意义，现引入面积坐标和体
积坐标的概念。如图2-10所示为ijl ?，),(i i y x 、),(j j y x 、),(l l y x 为其顶点坐标，则三角形ijl ?的面积为：
图2-10 三角形平面图
l
l j j i
i
ijl y x y x y x 1112
1
=
? 设),(y x p 点为三角形中任一点，则 l
j l
j l
l
j j i i i x x z
y y x
y x y x z c x b a 111
1+-=
++l
l
j j y x y x y
x 111= （2-33）
上式值为图2-10中三角形p j l ?面积的两倍，所以：
ijl
pjl
ijl i i i i y c x b a y x N ??=?++=
2)(),( （2-34）
同理可得：ijl
pil j y x N ??=
),(，ijl
pij l y x N ??=
),(
定义),(y x N i ，),(y x N j ，),(y x N l 为),(y x p 在三角形i j l ?内的面积坐标，具有以下性质：
1） 1),(),(),(=++y x N y x N y x N l j i （2-35） 2） 1),(≤y x N i 3） ??
?≠==n
k n k y x N n n k 0
1
),(
4） 与节点i 对边平行的直线上各点的),(y x N i 值不变。
利用面积坐标的定义，可以方便地确定三角单元中任一点的面积坐标，且可建立一个面积坐标网，如图2-11所示。因为面积坐标与通常使用的直角坐标系的选择无关，可利用该坐标网构造任意性线、非线性的插值函数。
图2-11 三角形面积坐标网
采用与上述相似的分析方法，可知),,(z y x N k (m l j i k ,,,=)为任一点
),,(z y x p 在四面体ijlm 内的体积坐标，具有与三角单元面积坐标相似的性质。
ijlm pjlm i D D z y x N =
),,( （2-36）
ijlm pilm j D D z y x N =
),,( （2-37）
ijlm pijm l D D z y x N =
),,( （2-38）
ijlm
pijl m D D z y x N =
),,( （2-39）
上式中ijlm D 为其对应节点所构成四面体的体积。
在式（2-31）假设四面体内电位为线性变化，P 点的电位值是通过该点的体积坐标将各节点的电位值线性化计算得出。如果四面体内的电位值为非线性变化，则只要改变相应点体积坐标的系数即可。
形如：
m m l l j j i i u QN u PN u ON u SN U +++= （2-40）
式中S 、O 、P 、Q 为P 点体积坐标对应的系数，由于四面体内电位的非线性变化规律比较复杂，系数值的确定比较困难，因此在以下各章节的讨论中只考虑电位值的线性变化，对于非线性的变化规律将以后讨论。
四、网格单元分析
整个求解区域Ω被剖分为一系列小四面体单元后，泛函J （u ）也
相应地被离散为各单元上的泛函)(u J e ，故有∑=m
e u J u J )()(。m 为四面体
单元总个数。
在任一四面体单元上泛函形式为：
?????ΓΩ+-??+??+??=e
e
dS u dxdydz fu z u y u x u
u J e 2222
21
}])()()[(2{)(σγσ
对任意节点的微分得：
?????ΓΩ??+??-???????+???????+???????=??e e
ds
u u
u dxdydz u u f z u u z u y u u y u x u u x u u u J k
k k k k k e σγσ})]()()([{)( （2-42）
其中m l j i k ,,,= 令 d x d y d z z
u
u z u y u u y u x u u x u k k k e
)]()()([
???????+???????+???????=Φ???Ωσ 则：
?????ΓΩ??+??-Φ=??e
e
ds u u
u dxdydz u u f u u J k k k e σγ)( （2-43）
现对各项分别进行讨论： 1、对于Φ项
由前面的讨论可知：∑=k
k k U N U ，)(1
z d y c x b a D
N k k k k k +++=
所以
m m l l j j i i
U x N U x N U x N U x N x U ??+??+??+??=?? （2-44） x
N x U
U k k ??=????)( （2-45） 又因为：
D
b x N k k =
??，k k N U U =?? 所以：
)(1)(2m m k l l k j j k i i k k U b b U b b U b b U b b D
x U U x U +++=??????? （2-46） 同理可得：
)(1)(2m m k l l k j j k i i k k U c c U c c U c c U c c D y U U y U +++=??????? （2-47） )(1)(2m m k l l k j j k i i k k U d d U d d U d d U d d D
z U U z U +++=??????? （2-48） 则：
)()()(z U
U z U y U U y U x U U x U k k k ???????+???????+??????? )]
()()[(1
2
m k m l k l j k j i k i m k m l k l j
k j i k i m k m l k l j k j i k i U d d U d d U d d U d d U c c U c c U c c U c c U b b U b b U b b U b b D +++++++++++=
]
)()()()[(1
2
m k m k m k m l k l k l k l j k j k j k j i k i k i k i U d d c c b b U d d c c b b U d d c c b b U d d c c b b D
+++++++++++=
（2-49） 记：
m
k m k m k m l k l k l k l j
k j k j k j i k i k i k i k U d d c c b b U d d c c b b U d d c c b b U d d c c b b F )()()()(+++++++++++=
所以：
k
e
k F V dv F D
e
362
σ
σ
=
=Φ???
Ω （2-50）
2、对于
dxdydz u u
f
e
k
???
Ω??项
由前面分析可知：)()()(000z z y y x x I f ---=δδδ，
k k
N U U
=??所以有： dxdydz z z y y x x z y x N I dxdydz u u
f
Ve
k Ve
k
)()()(),,(000---=????????
δδδ （2-51） 由δ函数的性质可知：
I z y x N dxdydz u u
f
k Ve
k
?=?????
),,(000 ???=)
,,(0
),,(k k k k k k z y x z y x I 供电点不在点供电点在点 （2-52）
3、对于
??
Γ??e
ds u u u
k σγ项
由式（2-25）分析可知，该项只与研究区域Ω的第三类边界条件有关，所以当四面体完全位于求解区域内部时（即非边界单元），该项值为零。在边界处则须对该项进行计算。
1） 关于γ值的计算
由电场理论知识可知，对于点源和均匀介质的情况，其电位函有如下的形式：
r
C
z y x U =
),,( （C 为常数） （2-53） 则：θc o s 3r
U n r r C n U -=?=?? （2-54）
式中θ为由点源到计算点的径向矢量r 和外法线矢量之间的夹角。式（2-54）可改写成如下形式：
0cos =+??θr
U
n U （2-55） 对比混合边界条件（2-7）的公式可知：
θγcos 1r
= （2-56）
设边界四面体单元e Ω如图2-9所示，由节点l j i ,,组成的面为边界面，另一节点m 位于研究区域内。
在三角形ijl ?内，γ值是坐标),,(z y x 及θcos 的函数，为简化计算，当三角形的面积较小时，用其重心处的γ值来代替整三角形的γ值。设三角形重心P 的坐标为),,(t t t z y x ，则：
???
?
??
???++=++=++=)(31)(31
)(31t j i t t j i t t j i t z z z z y y y y x x x x （2-57）
可知：
202020)()()(z z y y x x r t t t -+-+-= （2-58）
由空间解析几何与向量代数理论可知，op
op
n n n n ??=θcos
式中，op n
为由供电点O 到三角形重心的径向矢量，
),,(000z z y y x x n t t t op ---=
（2-59）
lj li n n n ?=
其中li n 、lj n
分别为节点i l ,和节点j l ,所在直线的径向矢量。
),,(l i l i l i li z z y y x x n ---=
（2-60） ),,(l j l j l j lj z z y y x x n ---=
（2-61）
可得：lj li n n
?
k
x x y y y y x x j z z x x x x z z i y y z z z z y y l j l i l j l i l j l i l j l i l j l i l j l i
)])(())([()])(()
)([()])(())([(-----+-----+-----= （2-62） 令 ))(())((l j l i l j l i y y z z z z y y A -----=
))(())((l j l i l j l i z z x x x x z z B -----=
))(())((l j l i l j l i x x y y y y x x C -----=
则： ),,(C B A n =
所以：
222202020000])()()[()()()(cos C
B A z z y y x x z z
C y y B x x A r t t t t t t ++?-+-+--+-+-=
=θ
γ 2） 关于
??
Γ??e
ds u u u
k
值的计算
由式（2-31）计算可知
k k
N U U
=??，所以： dS U N N U N N U N N U N N dS U U
U
e
e
m k m l k l j k j i k i k
????ΓΓ+++=??)( （2-64） 为方便计算，设四面体单元e Ω的边界面ijl ?所在的平面为XOY 平面，l 为原点，li 为X 轴，并设节点m l j i ,,,的坐标分别为)0,0,(i x ，)0,,(j j y x ，)0,0,0(和
),,(m m m z y x 。再设三角形内的一动点P 的坐标为)0,,(y x 。
可知：m j i m
m j m j i z y x z y y x x x D ==
0001111，2
j i y x e =
?
由式（2-36）~（2-39）所示的体积坐标可知：
j i j i m
m j m j i y x y x x x z y y y x x x
D N -
=
=
000001111
1,j
m
m m i
j y y z y y x x x D N =
=0
00
00
01111
1 i
j
i i j m
m j m j i l x x
y y x x x z y y y x x x x D N -
+-==
1000
01111
1,00
00
001111
1==
y y x x x D N j j i m
因为1=+++m l j i N N N N 。所以在Δ上1=++l j i N N N 。故i N ，j N ，l N 中
只有两个自由量，设为i N ，j N ，则
dxdy N
N j
i ???
dx y y y x y x x x dy j
j i j y y
y x x x y x i j
j
i j i j
j
)(
-=??-+
dy y x x x y x y x y y x y y x x x y x y
i
i i j i j j j j i j i y j i j
)}(])()[(2{22220
----+=?
12
24
e
y x j i ?=
=
（2-65） dx y y dy dxdy N N
j
j
i j i j
j
y y
y x x x y x j
j j
????-+
?
=0
22
6
12
e
y x j i ?=
=
（2-66） 因为i N ，j N ，l N 中只有两个自由量，而且三角形Δ三个节点i,j,l 的选取排列是任意的。因此：
l j i t t
k e
t k e dxdy N N t k ,,6
12
=???????=?≠?=???
（2-67）
所以：
ds u N N u N N u N N u N N ds u u
u
m k m l k l j k j i k i k
)(+++=????ΓΓ
σγσγ dxdy u N N u N N u N N u N N m k m l k l j k j i k i )(+++=??
σγ
m
k l k j k i k u u u u e e e e e e e e e m
l j i ====???????
?
???????
??????????
??????????????????=0000061212012612012126σγσγσγσγσγσγσγσγσγ （2-68）
记上式为k ψ
五、系数矩阵合成
由单元分析可知：
k k k k e N I F Ve
u u J ψ+?-=??36)(σ
（2-69） 其中k ψ项只对参于对边界单元的计算。 由泛函极值的必要条件
0)
(=??k
e u u J 可知： 036=ψ+?-k k k N I F Ve
σ
即：
k k k N I F Ve
?=ψ+36σ
（2-70）
式中 ??
?≠==)
,,(),,(0
),,(),,(1
000000z y x z y x z y x z y x N k k k k k k k 供电点不在节点供电点在节点
对m l j i k ,,,=，上式可写成如下的矩阵形式[4]：
?????
???????=?????????????????????
????
?m l j i m l j i e mm e lm e ll e jm e
jl e jj e
im
e il e ij e ii IN IN IN IN u u u u k k k 对k k k k k k k 称 （2-71） 系数矩阵即为四面体单元Ve 的单元刚度矩阵e K ，显然单元刚度矩阵为对称阵，其阶数等于单元节点数。
在具体求解计算过程中，应将单元系数矩阵合成总的系数矩阵K 。即
∑=m
e K K ,其中 C K C K e T e =为单元系数矩阵在总系数矩阵中的形式；C
为n ?4阶的转换矩阵，其作用将单元体系数矩阵变成总系数矩阵，其元素值为11=i C ，12=j C ，13=l C ，14=m C 其余值全为零。整个求解区域离散为四面体单元，每一个单元体都可求得其对应的系数矩阵，将各单元刚度矩阵叠加即可得总系数矩阵。设整个求解区域有n 个节
点，则K 是一个n n ?阶的系数矩阵。合成总矩阵后，可得由n 个方程组成的大型线性方程组，其矩阵表达形式为：
???
??
?
????
??????????=???????????
?????????????????????????????----------1n i 101n i 101n 1,n j 1,n 1,1n 1,0
n 1n i,ij i1i01n 1,1j 11101n 0,0j 0100IN IN IN IN u u u u k k k k k k k k k k k k k k k k
（2-72） 可记作：
F
KU = （2-73）
系数矩阵K 的特点
总系数矩阵K 的阶数为n n ?，其数据量较大，在计算时对计算机的内存要求较高、计算时间较长。为解决这一问题，可以利用总系数矩阵的某些特性，得到一定程度的解决。有限元总系数矩阵具有以下特点：
1)对称性 总系数矩阵是一个对称矩阵,运行时系数矩阵可以只存贮矩阵的上三角部分或下三角部分,可节省一半的存贮量。
2）稀疏性 由于有限元方法采用四面体单元结点的相互关联性，这决定了总系数矩阵具有大型、稀疏的特征。实际上，只有当结点i 是结点j 的相邻点时0k ≠ij 即总系数矩中的某一行中只有几个非零元素，对于一个四面体单元的四个结点，与任一结点对应的整体刚度矩阵中最多只有15个非零元素。当空间结点数较多，网格剖分较密，总系数矩阵很大时其矩阵的稀疏性表现得越突出。另一方面，总系数矩阵的结构（即非零元素的排列方式）依赖于结点的编号规律。相邻结点的编号相差较大，其系数矩阵的稀疏性越明显。
3)带状分布特性 对于有限元数值分析中的总系数矩阵,非零元素分布在
以主对角线为中心的斜带形区域中,由于对称性,存贮系数时,可只取上三角(或下三角)部分,因而每行元素只存贮半个带宽的元素就可以了。
4)正定性有限元数值分析中,方程离散后形成的系数矩阵具有正定性特点,即主对角线元素在所处行与列为最大且为正值。
六、求解线性方程组
有限元法是求解复杂电磁场问题的有效工具之一，大型工程电磁场有限元计算最后都归结为大型稀疏有限元方程的求解。在求解区域上对泊松方程进行离散并设置各节点参数及其边界条件之后，便得到一个大型的线性方程组。求解此方程组，就可以得到各节点的电位值，进而计算出不同装置形式下的视电阻率值。
目前求解线性方程组F
KU 的方法很多：
1）、Gauss-Sidel迭代法；
2）、超松驰迭代法；
3）、高斯消去法（Gauss）；
4）、乔列斯基分解法(LDL T)等。
这些方法在求解线性方程组方面已经取得了较为广泛的应用，并且具有计算精度较高、稳定性强等优点。本文在程序设计的过程中，主要采用了常规的Gauss消去法和适合大型稀疏矩阵方程的LDL T分解法，并对其计算精度进行对比。
1、Gauss消去法
Gauss消去法是计算机上常用的求解线性方程组的有效算法，此方法分为消元过程和回代过程。消元过程是将原方程组化为上三角形方程组的过程，而
回代过程是求解上三角形方程的过程
对于线性方程组：
F KU = （3-1）
其中：n n )k (K ?=ij 为对称正定型稀疏矩阵
T n 21)u ,......,u ,(u U = T n 21)f ,......,f ,(f F =
为方便起见，计（3-1）式为(1)(1)F U K =，即：
??
?
??
??
?????????=????????????????????????????????(1)n (1)2(1)1n 21(1)nn (1)n2
(1)n1
(1)2n (1)22(1)21(1)
1n
(1)12
(1)11f f f u u u k k k k k k k k k （3-2）
第一次消元：
设0k (1)
11≠，记(1)11(1)11/k k l i i =，
（n ,,3,2 =i ），以1l i -乘式（3-2）中的第一个方程加到第i （n ,,3,2 =i ）个方程上，得(2)(2)F U K =，即：
??
?
??
??
?????????=???????????????????????????????
?(2)n (2)2(1)1n 21(2)nn (2)n2
(2)2n (2)22(1)
1n
(1)12
(1)11f f f u u u k k k k k k k （3-3）
其中：)
1(11)1()2(k l k k j i ij ij -= （n ,,3,2, =j i ）
)1(11)1()2(f l f f i i i -= （n ,,3,2 =i ） 第二次消元：
设0k )
2(22≠，记)2(22)2(22k /k l i i =，（n ,,4,3 =i ），以2l i -乘式（3-3）中的第二个方程加
到第i （n i ,,4,3 =）个方程上，得(3)(3)F U K =，即：
??
?
????
?????????=???????????????????????????????
?(3)n (3)3(2)2(1)1n 21(3)nn (3)n3
(3)3n
(3)
33
(2)2n (2)23
(2)22(1)
1n
(1)13
(1)12(1)11f f f f u u u k k k k k k k k k k k （3-4） 第1-k 次消元完成后，得与式（3-1）等价的方程组，)()(F U K k k =，即：
???
?
??
?????????????
?=???????????????????????????????????????
?)(n )((2)2(1)1n 21)
(nn )(n )(n (2)2n (2)22
(1)
1n
(1)
12
(1)11f f f f u u u u k k k k k k k k k k k k k k k k
k k (k)kk
（3-5） 其中：)
1(11)1()(k l k k -----=k j k ik k ij k ij （n ,,1,, +=k k j i ） )
1(11)1()(f l f f -----=k k ik k i k i （n ,,1
, +=k k i ） 如此进行下去，共经过1-n 次消元计算后得到与式（3-1）等价的上三角形方程组，即：
??
?
??
??
?????????=???????????????????????????????
?(n)n (2)2(1)1n 21(n)nn (2)2n (2)23
(2)22(1)
1n
(1)13
(1)12
(1)11f f f u u u k k k k k k k k （3-6） 然后利用回代公式求解线性方程组的解为：
??
???-==∑+=)(n 1)()()
()(k /)u k f (u k /f u i ii
i s s i is i i i n nn n n n （1,2,,2n ,1n --=i ） （3-7） 2、LDL T 分解法
对有限元方程组F KU =的系数矩阵K 采用半带宽压缩存储技术，用二维数
组)M KN(N,存放系数矩阵K ，其中N 为网格节点总数，M 为半带宽。
将系数矩阵K 矩阵分解为：
T LDL K = （3-8） 式中，L 为带状的下三角矩阵，0l =ij )M (>-j i ；D 为对角阵ii
ii l 1d =
则： ∑
-=-=1
l l l k l j t t tt
jt it ij ij m
)t ;N 3,2,1(i j i m == （3-9）
式中：??
?+>-+≤=1
M M
1M 0
t i i i m （3-10）
于是，解有限元方程F KU =等价于解方程F U LDL T =，U L Y T =，
则得方程组： ???==Y
U L F
LD Y T （3-11）
)l f l (
f y 1
∑-=-=i t t tt
t
it i i m
)N ,...,2,1(-i （3-12） 式中m t 值与式（3-10）相同，最后，矩阵U 各元素的计算公式为：
ii i t t ti
i i m
l /)u l
y (u t 1
∑+=-
= 1,,1N N, -=i （3-13）
??
?+-≤+++->-=1
M N 1
M 1M N 1N t i i i m （3-14）
为比较以上两种方程组求解方法的计算精度及其稳定性，令式（2-73）中的U 向量值为1，求解右端项F 。在程序计算过程中，利用总系数矩阵K 和所得的右端项F ，计算向量U 的值，比较U 与1的大小关系，即可判断比较方程组求解方法的计算精度及其稳定性。具体计算结果如表3-2所示。
表3-2 计算结果比较
经过对上述两种方程组求解算法的计算结果比较可知，在相同条件下LDL T 法计算精度较高，稳定性好，但速度较慢；Gauss 消去法计算速度快，但其计算精度较低，稳定性差。因此在计算过程中将采用LDL T 分解算法。
3、LU 分解法
高斯消元的过程相当于下述矩阵乘法运算
[][](2)(1)||M M A b U
y =
因此由分块矩阵乘法可得
(2)(1)M M A U =(2)(1)M M b y = （3-2-15）
令 (2)(1)1(1)
1()()()
L M M M M ---== 可见只要消元过程能进行到底，就有以下等价关系
?
??==??→←=??→←===y Ux b
Ly b LUx b Ax Ux
y LU
A （3-2-16）
消元过程相当于分解A 为单位下三角阵L 与上三角阵U 的乘积，解方程组Ly=b ，回代过程就是解方程组Ux=y 。
以上分析与结论完全可以推广到n 阶线性方程组，相应的式子（3-2-16）中的L 为n 阶单位下三角阵,U 为上三角阵
2112111n n l L l l ??????=?????? 11
12
1(1)
(1)22
2(1)n n n nn a a a a a U a -??????=??
???
?
解线性方程组Ax=b 可利用矩阵A 的LU 分解转化为解两个三角形方程，这种解方程组的方法称为LU 分解法或三角状分解法。关于矩阵A 的LU 分解的存在唯一条件有如下两个定理。
由前面的分析可知，可以用顺序消去法求矩阵A 的LU 分解，为了简化分解，下面通过比较元素法导出一个直接求L 与U 的元素的公式。令
A=LU （3-2-17） 其中
matlab有限元分析实例
MATLAB: MATLAB是美国MathWorks公司出品的商业数学软件，用于数据分析、无线通信、深度学习、图像处理与计算机视觉、信号处理、量化金融与风险管理、机器人，控制系统等领域。 MATLAB是matrix&laboratory两个词的组合，意为矩阵工厂（矩阵实验室），软件主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式。 MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等。MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C++，JAVA的支持。 MATLAB有限元分析与应用:
《MATLAB有限元分析与应用》是2004年4月清华大学出版社出版的图书，作者是卡坦，译者是韩来彬。 内容简介： 《MATLAB有限元分析与应用》特别强调对MATLAB的交互应用，书中的每个示例都以交互的方式求解，使读者很容易就能把MATLAB用于有限分析和应用。另外，《MATLAB有限元分析与应用》还提供了大量免费资源。 《MATLAB有限元分析与应用》采用当今在工程和工程教育方面非常流行的数学软件MATLAB来进行有限元的分析和应用。《MATLAB有限元分析与应用》由简单到复杂，循序渐进地介绍了各种有限元及其分析与应用方法。书中提供了大量取自机械工程、土木工程、航空航天工程和材料科学的示例和习题，具有很高的工程应用价值。
有限元理论方法
关于有限元分析法及其应用举例 摘要：本文主要介绍有限元分析法，作为现代设计理论与方法的一种，已经在 众多领域普遍使用。介绍了它的起源和国内外发展现状。阐述了有限元法的基 本思想和设计方法。并从实际出发，例举了有限元法的一个简单应用———啤 酒瓶的应力分析和优化,表明了利用有限元分析法的众多优点。随着计算机的 发展，基于有限元分析方法的软件开发越来越多。本文也在其软件开发方面进 行阐述，并简单介绍了一下主流软件的发展情况和使用范围。并就这一领域的 未来发展趋势进行阐述。 关键词：有限元分析法软件啤酒瓶 Abstract：This thesis mainly introduces the finite element analysis, as a modern design theory and methods used widely in in most respects. And this paper introduces its origins and development in world. It also expounds the basic thinking and approach of FEM..Proceed from the actual situation,this text holds the a simple application of finite-element method———the analysis and optimized of an beer bottle and indicate the the numerous benefits of finite element analysis .As computers mature and based on the finite element analysis of the software development is growing. This article introduces its application in the software development aspects as well, and briefly states the development and scope of the mainstream software. And it’s also prospect future development tendency in this area . Key: Finite Element Analysis Software Beer bottle 0 绪论 有限元法(Finite Element Method，FEM)，是计算力学中的一种重要的方法，它是20世纪50年代末60年代初兴起的应用数学、现代力学及计算机科学相互渗透、综合利用的边缘科学。有限元法最初应用在工程科学技术中，用于模拟并且解决工程力学、热学、电磁学等物理问题。对于过去用解析方法无法求解的问题和边界条件及结构形状都不规则的复杂问题，有限元法则是一种有效的分析方法。有限元法的基本思想是先将研究对象的连续求解区域离散为一组有限个且按一定方式相互联结在一起的单元组合体。由于单元能按不同的联结方式进行组合，且单元本身又可以有不同形状，因此可以模拟成不同几何形状的求解小区域；
三维有限元法计算过程
三维有限元法计算过程 三维有限元法的计算过程： 1）网格单元剖分； 2）线性插值； 3）单元分析； 4）总体刚度矩阵合成； 5）求解线性方程组等部分组成。 一、偏微分方程对应泛函的极值问题 矿井稳恒电流场分布示意图 主要任务是分析在给定边界条件下，求解稳定电流场的Laplace 方程或Poisson方程的数值解，即三维椭圆型微分方程的边值问题：
) ()((0)(0 )()()(000z z y y x x I F u n u n u F z u z y u y x u x Lu w D ---=???? ?????=+??=??=????+????+????≡ΓΓ+Γδδδγσσσ 上述微分方程边值问题等价于下面泛函的极小值问题： dS U dxdydz fU z U y U x U U J w D ?????Γ+Γ+ΓΩ +-??+??+??=222221 }])()()[(2{][γσσ 二、网格剖分 ∞1 ρi i h ρ......... ... 1、网格单元的类型 图2-5 网格单元类型 2、网格单元剖分原则及其步长选择 因此，网格内的单元剖分应按以下剖分原则 1）、各单元节点（顶点）只能与相邻单元节点（顶点）重合，而
不能成为其它单元内点； 2）、如果求解区域对称，那么单元剖分也应该对称； 3）、在场变化剧烈的区域网格剖分单元要密一些，在场变化平缓 的区域单元密度应小。 4）、网格单元体的大小变化应逐步过渡。 根据上述剖分原则，以x 、y 、z 坐标轴原点o 为中心，分别向x 、y 、z 方向的两侧作对称变步长剖分，距o 越远，步长应越大。常用的变步长方法有： c i x x i i )1(1+=?-?+ c x x i i =??+/1(i ≠0) c x x i i =?-?+1 1 1(i ≠0) 以上各式中c 为常数，1+?i x 、i x ?为同一坐标轴上相邻步长值。以x 方向为例，可知，x 正方向与负方向对称，只相差一负号。若令00=?x ，只要给出距原点最近节点的坐标1x ?，由上式即可求出其它相应的步长i x ?。同理可求得y 、 z 方向上的变步长i y ?、i z ?。 3、网格剖分方法 图2-6 平面内节点编号示意图
有限元法的基本思想及计算 步骤
有限元法的基本思想及计算步骤 有限元法是把要分析的连续体假想地分割成有限个单元所组成的组合体，简称离散化。这些单元仅在顶角处相互联接，称这些联接点为结点。离散化的组合体与真实弹性体的区别在于：组合体中单元与单元之间的联接除了结点之外再无任何关联。但是这种联接要满足变形协调条件，即不能出现裂缝，也不允许发生重叠。显然，单元之间只能通过结点来传递内力。通过结点来传递的内力称为结点力，作用在结点上的荷载称为结点荷载。当连续体受到外力作用发生变形时，组成它的各个单元也将发生变形，因而各个结点要产生不同程度的位移，这种位移称为结点位移。在有限元中，常以结点位移作为基本未知量。并对每个单元根据分块近似的思想，假设一个简单的函数近似地表示单元内位移的分布规律，再利用力学理论中的变分原理或其他方法，建立结点力与位移之间的力学特性关系，得到一组以结点位移为未知量的代数方程，从而求解结点的位移分量。然后利用插值函数确定单元集合体上的场函数。显然，如果单元满足问题的收敛性要求，那么随着缩小单元的尺寸，增加求解区域内单元的数目，解的近似程度将不断改进，近似解最终将收敛于精确解。 用有限元法求解问题的计算步骤比较繁多，其中最主要的计算步骤为： 1）连续体离散化。首先，应根据连续体的形状选择最能完满地描述连续体形状的单元。常见的单元有：杆单元，梁单元，三角形单元，矩形单元，四边形单元，曲边四边形单元，四面体单元，六面体单元以及曲面六面体单元等等。其次，进行单元划分，单元划分完毕后，要将全部单元和结点按一定顺序编号，每个单元所受的荷载均按静力等效原理移植到结点上，并在位移受约束的结点上根据实际情况设置约束条件。 2）单元分析。所谓单元分析，就是建立各个单元的结点位移和结点力之间的关系式。现以三角形单元为例说明单元分析的过程。如图1所示，三角形有三个结点i，j，m。在平面问题中每个结点有两个位移分量u，v和两个结点力分量F x，F y。三个结点共六个结点位移分量可用列
有限元实例分析
作业一：有限元分析实例 实例：请对一个盘轴配合机构进行接触分析。轴为一等直径空心轴，盘为等厚度圆盘，其结构及尺寸如图所示。盘和轴为一种材料，材料参数为：弹性模量Ex=2.5E5，泊松比NUXY=0.3，摩擦系数MU=0.25，试采用有限元计算方法分析轴和盘在过盈配合时的应力应变分布以及将轴从盘心拔出时轴和盘的接触情况。 问题分析说明 （1）本题主要分析装配过程中结构的静态响应，所以分析步选择通用静态分析步。由于为过盈配合，属于大变形，故应考虑几何 非线性的影响。 （2）模型具有轴对称性，所以可以采取轴对称模型来进行分析，先建立二维模型计算，再转换为三维模型计算，这样可以节省计
算时间。分析过程由两个载荷步组成, 第一个载荷步为过盈分 析, 求解过盈安装时的情况。第二个载荷步为将轴从盘心拔出 时的接触分析, 分析在这个过程中盘心面和轴的外表面之间的 接触应力。它们都属于大变形问题, 属于非线性问题。在分析 时需要定义一些非线性选项来帮助问题的收敛。 （3）接触面之间有很大的相对滑动，所以模型要使用有限滑移。 模型建立的分析说明 （1）进定义单元类型此项实例分析的问题中涉及到大变形, 故选用So li d185 单元类型来建立本实例入部件模块，的模型。盘 轴接触问题属于面面接触, 目标面和接触面都是柔性的, 将使用接触单元T ARGET 170 和CO NTAT17 4来模拟接 触面。分别创建名为为part1、part2的部件。 （2）定义材料属性，在线性各向同性材料属性对话框中的EX (弹性模量) 文本框中输入 2 . 5E5,PRX Y (泊松比) 文本框中输入 0 . 3,并将定义的材料属性赋予给part1和part2。如下图所示。 （3）进入装配模块，创建两者间的装配关系。
电磁仿真算中的有限元法
1电磁仿真算法中的有限元法 1.1常规的电磁计算方法简介 从上世纪50年代以来,伴随着计算机技术的进步,电磁仿真算法也蓬勃发展起来,这其中主要包括:单矩法、矩量法和有限元法等属于频域技术的算法; 传输线矩阵法、时域积分方程法以及时域有限差分法等属于时域技术的算法。除了这些以外, 还有属于高频技术的集合衍射理论等。本文根据国内外计算电磁学的发展状况,对日常生活中比较常用的电磁计算方法做了介绍,并对有限元法做了重点说明。 ⑴矩量法 矩量法属于电磁场的数值计算方法中频域技术的一种, 它的基本原理是利用把待解的微积分方程转化成的算子方程, 然后将由一组线性组合表示的待求函数代入第一步中的算子方程, 然后将算子方程转化成矩阵方程, 最后再通过计算机进行大量的数值计算从而得到数值结果。该方法在求解非均勻和不规则形状对象时,面很广,但会生成病态矩阵,所以会在一定程度上受到限制。矩量法的特点就是适用于求解微积分方程, 并且求解方法统一简单。但缺点就是会占用大量计算机内存,影响计算速度。 (2)单矩法 单矩法是一种解析方法和数值方法相结合的混合数值算法法,该方法的关键在于,如何合理的选择一个球面最小的半径,使得能够将分析对象的结构全部包含在内,以便将内外场进行隔离。外边的散射场单独使用其他函数表示,而包围的内部区域使用有限元法亥姆赫兹(Helmholtz)方程。此方法对于计算复杂形体乃至复杂埋入体内的电磁散射是种极为有效的手段。 (3)时域有限差分法 时域有限差分法(FDTD）近几年来越来越受到各方的重视, 因为一方面它处理庞大的电磁福射系统方面和复杂结构的散射体时很突出,另外一方面则在于它不是传统的频域算法, 它是种时域算法, 直接依靠时间变量求解麦克斯韦方程组,可以在有限的时间和体积内对场进行数据抽样, 这样同时也能够保证介质边界
有限单元法与有限元分析
有限单元法与有限元分析 1.有限单元法 在数学中，有限元法（FEM，Finite Element Method）是一种为求解偏微分方程边值问题近似解的数值技术。求解时对整个问题区域进行分解，每个子区域都成为简单的部分，这种简单部分就称作有限元。它通过变分方法，使得误差函数达到最小值并产生稳定解。类比于连接多段微小直线逼近圆的思想，有限元法包含了一切可能的方法，这些方法将许多被称为有限元的小区域上的简单方程联系起来，并用其去估计更大区域上的复杂方程。它将求解域看成是由许多称为有限元的小的互连子域组成，对每一单元假定一个合适的(较简单的）近似解，然后推导求解这个域总的满足条件(如结构的平衡条件），从而得到问题的解。这个解不是准确解，而是近似解，因为实际问题被较简单的问题所代替。由于大多数实际问题难以得到准确解，而有限元不仅计算精度高，而且能适应各种复杂形状，因而成为行之有效的工程分析手段。 随着电子计算机的发展，有限单元法是迅速发展成一种现代计算方法。它是50年代首先在连续体力学领域--飞机结构静、动态特性分析中应用的一种有效的数值分析方法，随后很快广泛的应用于求解热传导、电磁场、流体力学等连续性问题。 1.1.有限元法分析本质 有限元法分析计算的本质是将物体离散化。即将某个工程结构离散为由各种单元组成的计算模型，这一步称作单元剖分。离散后单元与单元之间利用单元的节点相互连接起来；单元节点的设置、性质、数目等应视问题的性质，描述变形形态的需要和计算精度而定（一般情况单元划分越细则描述变形情况越精确，即越接近实际变形，但计算量越大）。所以有限元中分析的结构已不是原有的物体或结构物，而是同新材料的由众多单元以一定方式连接成的离散物体。这样，用有限元分析计算所获得的结果只是近似的。如果划分单元数目非常多而又合理，则所获得的结果就与实际情况相符合。 1.2.特性分析 1）选择位移模式： 在有限单元法中，选择节点位移作为基本未知量时称为位移法；选择节点力作为基本未知量时称为力法；取一部分节点力和一部分节点位移作为基本未知量时称为混合法。位移法易于实现计算自动化，所以，在有限单元法中位移法应用范围最广。 当采用位移法时，物体或结构物离散化之后，就可把单元总的一些物理量如
对有限元方法的认识
我对有限元方法的认识 1有限元法概念 有限元方法（The Finite Element Method, FEM）是计算机问世以后迅速发展起来的一种分析方法。每一种自然现象的背后都有相应的物理规律，对物理规律的描述可以借助相关的定理或定律表现为各种形式的方程（代数、微分、或积分）。这些方程通常称为控制方程（Governing equation）。 针对实际的工程问题推导这些方程并不十分困难，然而，要获得问题的解析的数学解却很困难。人们多采用数值方法给出近似的满足工程精度要求的解答。 有限元方法就是一种应用十分广泛的数值分析方法。 有限元方法是处理连续介质问题的一种普遍方法，离散化是有限元方法的基础。 这种思想自古有之：古代人们在计算圆的周长或面积时就采用了离散化的逼近方法：即采用内接多边形和外切多边形从两个不同的方向近似描述圆的周长或面积，当多边形的边数逐步增加时近似值将从这两个方向逼近真解。 近年来随着计算机技术的普及和计算速度的不断提高，有限元分析在工程设计和分析中得到了越来越广泛的重视，已经成为解决复杂的工程分析计算问题的有效途径，现在从汽车到航天飞机几乎所有的设计制造都已离不开有限元分析计算，其在机械制造、材料加工、航空航天、汽车、土木建筑、电子电器、国防军工、船舶、铁道、石化、能源、科学研究等各个领域的广泛使用已使设计水平发生了质的飞跃。 国际上早在 60 年代初就开始投入大量的人力和物力开发有限元分析程序。“有限单元”是由Clough R W于1960年首次提出的。但真正的有限元分析软件是诞生于 70 年代初期，随着计算机运算速度的提高，内、外存容量的扩大和图形设备的发展，以及软件技术的进步，发展成为有限元分析与设计软件，但初期其前后处理的能力还是比较弱的，特别是后处理能力更弱。
有限元例题
【1】图示弹性力学平面问题，采用三角形常应变元，网格划分及单元、节点编号如图1所示。试求： (1) 计算系统刚度矩阵的最大带宽； (2) 根据图中结构的边界约束状态，给出约束节点位移值。 【解】 (1) 相邻节点号的最大差为d = 4； 所以，半带宽为B = 2 ? (4 + 1) = 10。 (2) u1 = 0，v1 = 0，u4 = 0，v4 = 0。 【2】弹性力学平面问题4节点等参元，其单元自由度是多少？单元刚度矩阵是多少阶的？单元刚度矩阵有多少个元素？ 【解】平面问题4节点等参元，其单元自由度是4 ?2 = 8个；单元刚度矩阵是8 ? 8 阶的，单元刚度矩阵有64个元素。
【3】平面刚架结构梁单元（考虑轴向和横向变形）的自由度是多少？单元刚度矩阵是多少阶的？单元刚度矩阵有多少个元素？ 【解】平面刚架结构梁单元（考虑轴向和横向变形）的自由度是2 ? 3 = 6个；单元刚度矩阵是6 ? 6阶的；单元刚度矩阵有36个元素。 【4】已知一等截面直杆中某一微段的起始点坐标为0.5m，终点坐标为0.6m，起始点的位移为0.2mm，终点的位移为0.3mm。假定直杆内的位移是线性分布的。求该微段等截面直杆的位移表达式f(x)。 【解】已知：x i = 0.5m, x j= 0.6m, u i = 0.2mm = 0.2?10-3m, u j= 0.3mm = 0.3?10-3m。 即
【5】已知4节点一维问题中单元①(1, 2)的应力矩阵为 结构总体位移列阵为 求单元①的应力（用矩阵计算）。 【解】由总体结构位移列阵知，单元①的位移列阵为 由{σ} = [C] {?}e可求得单元①的应力
三维有限元方法-为一种新型的研究方法,是利用数学的形式概括事件的各种条件和性能
三维有限元方法-为一种新型的研究方法，是利用数学的形式概括事件的各种条件和性能 三维有限元方法-为一种新型的研究方法，是利用数学的形式概括事件的各种条件和性能，并进行重复分析计算的研究方法。有限元数值模型分析技术将现代数学、力学的基础理论与有限元分析技术、计算机图形学和优化技术相结合，具有丰富、完善的单元库、材料模型库和求解器，可利用数值模拟技术高效求解各类结构动力、静力和线性、非线性问题。将其应用于骨科领域，可以更好的进行各种骨科生物力学分析，对各种生物力学强度进行数值模拟分析，较精确地掌握各点的受力情况，了解内部应力应变的分布规律，获得应力应变分布图等，从而更好的指导临床治疗。 学术术语来源—— 锁骨中段骨折修复：重建钢板前置与上置的生物力学差异 文章亮点： 1 文章结果显示，不论怎样的载荷条件，骨折断端均会存在一定的应力。而且，前置位和上置位不同内固定方式对骨折端愈合的影响方面不存在明显差别，但在骨折断端应力和内固定应力方面，前置位均显著大于上置位。即提示，较之上置位，前置位固定具有更明显的应力集中效应。 2 临床对锁骨中段骨折进行修复的过程中，利用不同重建钢板位置进行内固定修复会产生不同的生物力学情况。其中，较之重建钢板上置内固定，前置内固定修复效果更佳，是一种较为可靠的治疗方法。 3 文章仅对不同重建钢板位置的内固定效果进行了分析研究，并未考虑到不同骨折类型力学特性以及不同钢板类型等因素的影响。并在研究过程中假设螺钉为圆形杆，因此最终的研究结果可能存在导致内固定装置最大等效应力下降的情况，计算精确度存在一定的误差。另外，文章中对所使用的各种生物材料的力学特性均进行了假设，与客观情况存在较大的差异。因此，文章还存在一定的不足之处，还需要在今后的研究中不断予以完善，以提高研究结果的准确性和可信性，更好的为临床治疗提供可参考的依据。 关键词：
有限元分析方法
百度文库- 让每个人平等地提升自我 第1章有限元分析方法及NX Nastran的由来 有限元分析方法介绍 计算机软硬件技术的迅猛发展，给工程分析、科学研究以至人类社会带来急剧的革命性变化，数值模拟即为这一技术革命在工程分析、设计和科学研究中的具体表现。数值模拟技术通过汲取当今计算数学、力学、计算机图形学和计算机硬件发展的最新成果，根据不同行业的需求，不断扩充、更新和完善。 有限单元法的形成 近三十年来，计算机计算能力的飞速提高和数值计算技术的长足进步，诞生了商业化的有限元数值分析软件，并发展成为一门专门的学科——计算机辅助工程CAE（Computer Aided Engineering）。这些商品化的CAE软件具有越来越人性化的操作界面和易用性，使得这一工具的使用者由学校或研究所的专业人员逐步扩展到企业的产品设计人员或分析人员，CAE在各个工业领域的应用也得到不断普及并逐步向纵深发展，CAE工程仿真在工业设计中的作用变得日益重要。许多行业中已经将CAE分析方法和计算要求设置在产品研发流程中，作为产品上市前必不可少的环节。CAE仿真在产品开发、研制与设计及科学研究中已显示出明显的优越性： ?CAE仿真可有效缩短新产品的开发研究周期。 ?虚拟样机的引入减少了实物样机的试验次数。 ?大幅度地降低产品研发成本。 ?在精确的分析结果指导下制造出高质量的产品。 ?能够快速对设计变更作出反应。 ?能充分和CAD模型相结合并对不同类型的问题进行分析。 ?能够精确预测出产品的性能。 ?增加产品和工程的可靠性。 ?采用优化设计，降低材料的消耗或成本。 ?在产品制造或工程施工前预先发现潜在的问题。 ?模拟各种试验方案，减少试验时间和经费。 ?进行机械事故分析，查找事故原因。 当前流行的商业化CAE软件有很多种，国际上早在20世纪50年代末、60年代初就投入大量的人力和物力开发具有强大功能的有限元分析程序。其中最为著名的是由美国国1
有限元法基础试题
有限元法基础试题（A ） 一、填空题（5×2分） 1.1单元刚度矩阵e T k B DBd Ω = Ω? 中，矩阵B 为__________，矩阵D 为___________。 1.2边界条件通常有两类。通常发生在位置完全固定不能转动的情况为_______边界，具体指定有限的非零值位移的情况，如支撑的下沉，称为_______边界。 1.3内部微元体上外力总虚功： ()(),,,,e x x xy y bx xy x y y by d W F u F v dxdy δστδτσδ??=+++++??+(),,,,x x y y xy y x u v u u dxdy σδσδτδδ??+++??的表达式中，第一项为____________________的虚功，第二项为____________________的虚功。 1.4弹簧单元的位移函数1N +2N =_________。 1.5 ij k 数学表达式：令j d =_____，k d =_____，k j ≠，则力i ij F k =。 二、判断题（5×2分） 2.1位移函数的假设合理与否将直接影响到有限元分析的计算精度、效率和可靠性。（ ） 2.2变形体虚功原理适用于一切结构（一维杆系、二维板、三位块体）、适用于任何力学行为的材料（线性和非线性），是变形体力学的普遍原理。 （ ） 2.3变形体虚功原理要求力系平衡，要求虚位移协调，是在“平衡、协调”前提下功的恒等关系。 （ ） 2.4常应变三角单元中变形矩阵是x 或y 的函数。 （ ） 2.5 对称单元中变形矩阵是x 或y 的函数。 （ ） 三、简答题（26分） 3.1列举有限元法的优点。（8分） 3.2写出有限单元法的分析过程。（8分） 3.3列出3种普通的有限元单元类型。（6分） 3.4简要阐述变形体虚位移原理。（4分） 四、计算题（54分） 4.1对于下图所示的弹簧组合，单元①的弹簧常数为10000N/m ，单元②的弹簧常数为20000N/m ，单元③的弹簧常数为10000N/m ，确定各节点位移、反力以及单元②的单元力。（10分） 4.2对于如图所示的杆组装，弹性模量E 为10GPa ，杆单元长L 均为2m ，横截面面积A 均为2×10-4m 2，弹簧常数为2000kN/m ，所受荷载如图。采用直接刚度法确定节点位移、作用力和单元②的应力。（10分）
5.3 三维静磁场的有限元分析
5.3 三维静磁场的有限元分析 5.3.1 边值问题 以标量磁位m ?表示的无源区磁场的边值问题与电位的拉普拉斯边值问题的数学表达形式完全一样，可以如前节所述的有限元分析。在此，考虑有电流存在以矢量磁位A 作为待求变量的有限元分析。 设在线性媒质中，磁场满足的边界条件：边界1S 面上有0A A =，在边界S 2面上取某种形式的对称面作为第二类齐次边界，在该面上磁场强度H 的切向分量为零： ()0=???=?n m n e A e H γ，有如下边值问题： ()()??? ??∈=???∈=∈=????2 10 s s V n m m e A A A J A γγ 5.3.2 场域剖分与插值 对于求解场域V ，根据其形状和场的定性分布，选择合适的单元(例如四面体单元)，进行场域剖分，得到0Z 个单元、0N 个节点。在单元e 内，对位函数A 进行插值。若采用四面体单元： ∑ == 4 1 j j e j N A A ~ 式中e j N 是单元形状函数，分量形式 z j zj e j y j yj e j x j xj e j z z y y x x A N A N A N A A A e e e e ~ e ~e ~A ~??? ? ? ?+ ???? ? ?+ ???? ? ?=++=∑ ∑ ∑ ===4 1 4 1 4 1 以矩阵表示磁矢量位A 在单元节点上的各分量 [][] ()z y x l A A A A A T l l l l e l ,,==43 21 ),,(][][~4 1 z y x l A N A N A e l T e j lj e j l ===∑=
ansys有限元建模与分析实例-详细步骤
《有限元法及其应用》课程作业ANSYS应用分析 学号： 姓名： 专业：建筑与土木工程
角托架的有限元建模与分析 一 、模型介绍 本模型是关于一个角托架的简单加载，线性静态结构分析问题，托架的具体形状和尺寸如图所示。托架左上方的销孔被焊接完全固定，其右下角的销孔受到锥形压力载荷，角托架材料为Q235A 优质钢。角托架材料参数为：弹性模量366E e psi =；泊松比0.27ν= 托架图（厚度：0.5） 二、问题分析 因为角托架在Z 方向尺寸相对于其在X,Y 方向的尺寸来说很小，并且压力荷载仅作用在X,Y 平面上，因此可以认为这个分析为平面应力状态。 三、模型建立 3.1 指定工作文件名和分析标题 （1）选择菜单栏Utility Menu → 命令.系统将弹出Jobname(修改文件名)对话框,输入bracket （2）定义分析标题 GUI ：Utility Menu>Preprocess>Element Type>Add/Edit/Delete 执行命令后，弹出对话框，输入stress in a bracket 作为ANSYS 图形显示时的标题。 3.2设置计算类型 Main Menu: Preferences … →select Structural → OK 3.3定义单元类型 PLANE82 GUI ：Main Menu →Preprocessor →Element Type →Add/Edit/Delete 命令，系统将弹出Element Types 对话框。单击Add 按钮，在对话框左边的下拉列表中单击Structural Solid →Quad 8node 82，选择8节点平面单元PLANE82。单击ok ，Element Types 对话框，单击Option ，在Element behavior 后面窗口中选取Plane strs w/thk 后单击ok 完成定义单元类型。 3.4定义单元实常数 GUI ：Main Menu: Preprocessor →Real Constants →Add/Edit/Delete ，弹出定义实常数对话框，单击Add ，弹出要定义实常数单元对话框，选中PLANE82单元后，单击OK →定义单元厚度对话框，在THK 中输入0.5.
有限元法与有限差分法的主要区别
有限元法与有限差分法的主要区别 有限差分方法(FDM)是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以Taylor级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。对于有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式。从差分的空间形式来考虑，可分为中心格式和逆风格式。考虑时间因子的影响，差分格式还可以分为显格式、隐格式、显隐交替格式等。目前常见的差分格式，主要是上述几种形式的组合，不同的组合构成不同的差分格式。差分方法主要适用于有结构网格，网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。构造差分的方法有多种形式，目前主要采用的是泰勒级数展开方法。其基本的差分表达式主要有三种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等，其中前两种格式为一阶计算精度，后两种格式为二阶计算精度。通过对时间和空间这几种不同差分格式的组合，可以组合成不同的差分计算格式。有限元方法的基础是变分原理和加权余量法，其基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解。采用不同的权函数和插值函数形式，便构成不同的有限元方法。有限元方法最早应用于结构力学，后来随着计算机的发展慢慢用于流体力学的数值模拟。在有限元方法中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线形组合来逼近单元中的真解，整个计算域上总体的基函数可以看为由每个单元基函数组成的，则整个计算域内的解可以看作是由所有单元上的近似解构成。在河道数值模拟中，常见的有限元计算方法是由变分法和加权余量法发展而来的里兹法和伽辽金法、最小二乘法等。根据所采用的权函数和插值函数的不同，有限元方法也分为多种计算格式。从权函数的选择来说，有配置法、矩量法、最小二乘法和伽辽金法，从计算单元网格的形状来划分，有三角形网格、四边形网格和多边形网格，从插值函数的精度来划分，又分为线性插值函数和高次插值函数等。不同的组合同样构成不同的有限元计算格式。对于权函数，伽辽金(Galerkin)法是将权函数取为逼近函数中的基函数；最小二乘法是令权函数等于余量本身，而内积的极小值则为对代求系数的平方误差最小；在配置法中，先在计算域内选取N个配置点。令近似解在选定的N个配置点上严格满足微分方程，即在配置点上令方程余量为0。插值函数一般由不同次幂的多项式组成，但也有采用三角函数或指数函数组成的乘积表示，但最常用的多项式插值函数。有限元插值函数分为两大类，一类只要求插值多项式本身在插值点取已知值，称为拉格朗日(Lagrange)多项式插值；另一种不仅要求插值多项式本身，还要求它的导数值在插值点取已知值，称为哈密特(Hermite)多项式插值。单元坐标有笛卡尔直角坐标系和无因次自然坐标，有对称和不对称等。常采用的无因次坐标是一种局部坐标系，它的定义取决于单元的几何形状，一维看作长度比，二维看作面积比，三维看作体积比。在二维有限元中，三角形单元应用的最早，近来四边形等参元的应用也越来越广。对于二维三角形和四边形电源单元，常采用的插值函数为有La grange插值直角坐标系中的线性插值函数及二阶或更高阶插值函数、面积坐标系中的线性插值函数、二阶或更高阶插值函数等。对于有限元方法，其基本思路和解题步骤可归纳为(1)建立积分方程，根据变分原理或方程余量与权函数正交化原理，建立与微分方程初边值问题等价的积分表达式，这是有限元法的出发点。(2)区域单元剖分，根据求解区域的形状及实际问题的物理特点，将区域剖分为若干相互连接、不重叠的单元。区域单元划分是采用有限元方法的前期准备工作，这部分工作量比较大，除了给计算单元和节点进行编号和确定相互之间的关系之外，还要表示节点的位置坐标，同时还需要列出自然边界和本质边界的节点序号和相应的边界值。(3)确定单元基函数，根据单元中节点数目及对近似解精度的要求，选择满足一定插
有限元法中的几个基本概念
诚信·公平·开放·共赢 Loyalty Fair Opening Win-win 有限元法中的几个基本概念 有限元法是把要分析的连续体假想地分割成有限个单元所组成的组合体，简称离散化。 这些单元仅在顶角处相互联接，称这些联接点为结点。 离散化的组合体与真实弹性体的区别在于：组合体中单元与单元之间的联接除了结点之外再无任何关联。但是这种联接要满足变形协调条件，即不能出现裂缝，也不允许发生重叠。显然，单元之间只能通过结点来传递内力。 通过结点来传递的内力称为结点力，作用在结点上的荷载称为结点荷载。当连续体受到外力作用发生变形时，组成它的各个单元也将发生变形，因而各个结点要产生不同程度的位移，这种位移称为结点位移。 在有限元中，常以结点位移作为基本未知量。并对每个单元根据分块近似的思想，假设一个简单的函数近似地表示单元内位移的分布规律，再利用力学理论中的变分原理或其他方法，建立结点力与位移之间的力学特性关系，得到一组以结点位移为未知量的代数方程，从而求解结点的位移分量。然后利用插值函数确定单元集合体上的场函数。显然，如果单元满足问题的收敛性要求，那么随着缩小单元的尺寸，增加求解区域内单元的数目，解的近似程度将不断改进，近似解最终将收敛于精确解。 附：FELAC 2.0软件简介 FELAC 2.0采用自定义的有限元语言作为脚本代码语言，它可以使用户以一种类似于数学公式书写和推导的方式，非常自然和简单的表达待解问题的微分方程表达式和算法表达式，并由生成器解释产生完整的并行有限元计算C程序。 FELAC 2.0的目标是通过输入微分方程表达式和算法之后，就可以得到所有有限元计算的程序代码，包含串行程序和并行程序。该系统采用一种语言(有限元语言)和四种技术(对象技术、组件技术、公式库技术生成器技术)开发而成。并且基于FELAC 1.0的用户界面，新版本扩充了工作目录中右键编译功能、命令终端输入功能，并且丰富了文本编辑功能，改善了用户的视觉体验，方便用户快速便捷的对脚本或程序进行编辑、编译与调试。其中并行版在前后处理上进行了相应的改进。
有限元计算原理与方法..
1.有限元计算原理与方法 有限元是将一个连续体结构离散成有限个单元体，这些单元体在节点处相互铰结，把荷载简化到节点上，计算在外荷载作用下各节点的位移，进而计算各单元的应力和应变。用离散体的解答近似代替原连续体解答，当单元划分得足够密时，它与真实解是接近的。 1.1. 有限元分析的基本理论 有限元单元法的基本过程如下： 1.1.1.连续体的离散化 首先从几何上将分析的工程结构对象离散化为一系列有限个单元组成，相邻单元之间利用单元的节点相互连接 而成为一个整体。单元可采用各种类 型，对于三维有限元分析，可采用四 面 体单元、五西体单元和六面体 单元等。在Plaxis 3D Foundation 程序中，土体和桩体主要采用包 含6个高斯点的15节点二次楔 形体单元，该单元由水平面为6 节点的三角形单元和竖直面为四 边形8节点组成的，其局部坐标 下的节点和应力点分布见图3.1，图3.1 15节点楔形体单元节点和应力点分布界面单元采用包含9个高斯点的 8个成对节点四边形单元。 在可能出现应力集中或应力梯度较大的地方，应适当将单元划分得密集些；
若连续体只在有限个点上被约束，则应把约束点也取为节点：若有面约束，则应 把面约束简化到节点上去，以便对单元组合体施加位移边界条件，进行约束处理； 若连续介质体受有集中力和分布荷载，除把集中力作用点取为节点外，应把分布 荷载等效地移置到有关节点上去。 最后，还应建立一个适合所有单元的总体坐标系。 由此看来，有限单元法中的结构已不是原有的物体或结构物，而是同样材料 的由众多单元以一定方式连接成的离散物体。因此，用有限元法计算获得的结果 只是近似的，单元划分越细且又合理，计算结果精度就越高。与位移不同，应力 和应变是在Gauss 积分点(或应力点)而不是在节点上计算的，而桩的内力则可通 过对桩截面进行积分褥到。 1.1. 2. 单元位移插值函数的选取 在有限元法中，将连续体划分成许多单元，取每个单元的若干节点的位移 作为未知量，即{}[u ,v ,w ,...]e T i i i δ=，单元体内任一点的位移为{}[,,]T f u v w =。 引入位移函数N (x,y,z )表示场变量在单元内的分布形态和变化规律，以便用 场变量在节点上的值来描述单元内任一点的场变量。因此在单元内建立的位移模 式为： {}[]{}e f N δ= （3-1） 其中：12315[][,,......]N IN IN IN IN =，I 为单位矩阵。 按等参元的特性，局部坐标(,,)ξηζ到整体坐标,,x y z （）的坐标转换也采用 与位移插值类似的表达式。经过坐标变化后子单元与母单元(局部坐标下的规则 单元)之间建立一种映射关系。不管内部单元或边界附近的单元均可选择相同的 位移函数，则为它们建立单元特性矩阵的方法是相同的。因此，对于15节点楔 形体单元体内各点位移在整体坐标系,,x y z （）下一般取：
有限元法基础重点归纳(精)
1、有限元这种数值计算方法起源于20世纪50年代中期航空工程中飞机结构的矩阵分析。 2、有限单元法的基本思想:在力学模型上将一个原来连续的物体离散成为有限个具有一定 大小的单元,这些单元仅在有限个节点上相连接,并在节点上引进等效力以代替实际作用于单元上的外力。 3、节点:网格间相互连接的点。 4、边界:网格与网格的交界线。 5、有限元的优点:①理论基础简明,物理概念清晰,且可在不同的水平上建立起对该法的 理解②具有灵活性和适用性,应用范围极为广泛③该法在具体推导运算中,广泛采用了矩阵方法。 6、有限单元法分类(从选择基本未知量的角度:位移法(以节点位移为基本未知量,通用 性广、力法(以节点力、混合法(一部分以节点位移,另一部分以节点力 7、有限元法分析计算的基本步骤:①结构的离散化②单元分析(选择位移模式,建立单元 刚度方程,计算等效节点力③整体分析④求解方程,得出节点位移⑤由节点位移计算单元的应变与应力。 8、单元划分:将某个机械结构划分为由各种单元组成的计算模型。 9、有限元法基本近似性------几何近似。
10、弹性力学的任务:分析弹性体在受外力作用并处于平衡状态下产生的应力、应变和位移状态及其相互关系等。 11、弹性力学假设所研究的物体是连续的、完全弹性的、均匀的、各向同性的、微小变形的和无初应力的 12、外力:体力(分布在物体体积内的力---重力、惯性力、电磁力面力(分布在物体表面上的力---流体压力、接触力、风力 13、应力:物体受外力作用,或由于温度有所改变,其内部发生的内力。σ={ σx σy σz τx τy τz } = [σx σy σz τx τy τz ]T 14、应变:物体受到外力作用时,其形状发生改变时的形变。---长度和角度。 ε={ εx εy εz γx γy γz } = [εx εy εz γx γy γz ]T 15、位移:弹性体在载荷作用下,不仅会发生形变,还将产生位移,即弹性体位置 的移动。 δ={u v w }=[u v w ]T 16:、变形协调条件:设想在变形前,把弹性体分为许多微小立方单元体。变形后,每个单元体都产生任意变形而变成一些六面体。可能发生这样的情况,这些六面体
CATIA有限元分析计算实例
CATIA有限元分析计算实例 11.1例题1 受扭矩作用的圆筒 11.1－1划分四面体网格的计算 （1）进入【零部件设计】工作台 启动CATIA软件。单击【开始】→【机械设计】→【零部件设计】选项，如图11－1所示，进入【零部件设计】工作台。 图11－1单击【开始】→【机械设计】→【零部件设计】选项 单击后弹出【新建零部件】对话框，如图11-2所示。在对话框内输入新的零件名称，在本例题中，使用默认的零件名称【Part1】。点击对话框内的【确定】按钮，关闭对话框，进入【零部件设计】工作台。 （2）进入【草图绘制器】工作台 在左边的模型树中单击选中【xy平面】, 如图11-3所示。单击【草图编辑器】工具栏内的【草图】按钮，如图11-4所示。这时进入【草图绘制器】工作台。 图11－2【新建零部件】对话框图11－3单击选中【xy平面】 （3）绘制两个同心圆草图 点击【轮廓】工具栏内的【圆】按钮，如图11-5所示。在原点点击一点，作为圆草图的圆心位置，然后移动鼠标，绘制一个圆。用同样分方法再绘制一个同心圆，如图11-6所示。 图11－4【草图编辑器】工具栏图11－5【轮廓】工具栏 下面标注圆的尺寸。点击【约束】工具栏内的【约束】按钮，如图11-7所示。点击选择圆，就标注出圆的直径尺寸。用同样分方法标注另外一个圆的直径，如图11-8所示。
图11－6两个同心圆草图图11－7【约束】工具栏 双击一个尺寸线，弹出【约束定义】对话框，如图11－9所示。在【直径】数值栏内输入100mm，点击对话框内的【确定】按钮，关闭对话框，同时圆的直径尺寸被修改为100mm。用同样的方法修改第二个圆的直径尺寸为50mm。修改尺寸后的圆如图11-10所示。 图11－8标注直径尺寸的圆草图图11－9【约束定义】对话框 （4）离开【草图绘制器】工作台 点击【工作台】工具栏内的【退出工作台】按钮，如图11-11所示。退出【草图绘制器】工作台，进入【零部件设计】工作台。 图11－10修改直径尺寸后的圆图11－11【工作台】工具栏 （5）拉伸创建圆筒 点击【基于草图的特征】工具栏内的【凸台】按钮，如图11-12所示。弹出【凸台定义】对话框，如图11-13所示。在【第一限制】选项组内的【长度】数值栏内输入50mm，点击对话框内的【确定】按钮，生成一个圆筒体，如图11-14所示。在左边的模型树上出现【填充器.1】元素。