有限单元法的数学基础
有限单元法的数学基础
1、引言
有限元方法归根结底是一种数值计算方法,它有严格的数学证明作为其近似的客观性和合理性的保证。力学问题最终归结为一组微分方程的边值问题或者初值问题抑或是混合问题。比如弹性静力学最终归结为L-N 方程的微分提法。在很难或者根本不可能得到所得方程的理论解的情况下,究竟用什么样的方法才能得到方程的近似解(这种近似解已经能够满足实际工程的需要),在这种情况下,二十世纪五六十年代由结构力学家进而由数学家提出和证明了这种思想方法的合理性。有限元方法产生于力学计算,但是,它本质上并不是力学的专利。世间万物的变化过程很多都可以通过微分方程特别是偏微分方程来描述,也就是说,微分方程是很多现象和过程的数学结构,而大多数的微分方程是不能得到理论解的,这时候就可以使用有限元方法来求其近似解,因为有限元方法是求解微分方程(组)的数值计算方法。它适用于力学的微分方程,也同样适用于其它领域的相应的微分方程的数值求解。 2、有限元方法数学根源
对于一个给定的微分方程定解问题,为了求其近似解,我们可以使用Ritz 方法和Galerkin 方法。下面分别阐述这两种方法,然后讨论有限元方法和他们的关系。
(1) Ritz 法
Ritz 法源于最小势能原理,设H 是可分的Hilbert 空间,在H 中取有限维空间Sn ,它是由N 个线性无关向量φ1, φ2, , φN 张成,即:
S N
⎧
≡⎨ωn , ωn =⎩
N
∑C φ, ∀(C
i i
i =1
⎫
, C , C ) ∈R 12N N ⎬
⎭
用S N 代替H ,在S N 上求泛函J(w)的极值, 即求U N ∈S N ,使得 J (U N ) = m in ω
∈S N
N
J (ωN )
实际上寻求U N 只需通过解一个线性方程组 J (ω) =
12
(D ω, ω) -F (ω) ≥0
D--------双线性形式 F--------线性泛函
N
ωN =
∑C φ
i
i =1
i
N
N
N
J (ωN ) = =
121
D (∑C i φi , ∑C i φi ) -F (∑C i φi )
i =1N
i =1
N
i =1
i
i
∑2
D (φi , φj ) C i C j -
∑F (φ) C
i =1
i , j =1
因此,J (ωN ) 是一个以C 1, C 2, , C N 为未知数(自变量)的二次多项式
j (C 1, C 2, , C N ) ,如果二次项的系数矩阵
φi φ, [D (
j
) , i ]=j 1,
2
,
是正定的,那么j =j (C 1, C 2, , C N ) 在N+1维空间是一个N ,
开口向上的椭球抛物面,它有且只有一个极(最)小值点,所谓在S N 上求J (ωN ) 的极值,就是确定C 01, C 02, , C 0N ,使得:
j (C 1, C
02
, , C
0N ) =min C
1, , C N
∈R
j (C 1, C
00
, , C 2
0N
)
极值条件:
∂j ∂C i
|C 0
1, C 2
, , C
0N
=0 (i =1, , N )
得:
n
∑D (φφ
i
i =1
) C j
0i =F (φi ) (i =1, , N )
000T
即:C =[C 1, C 2, , C N ] 适合方程组:
KC=F
T
F =[F (φ1), , F (φ1)]
,D (φN , φ) 1⎫⎛D (φ1, φ1) ,D (φ, 2φ) ,1
⎪D (φ1, φ2) ,D (φ, 2φ) , ,D (φN , φ) 2
2⎪ K =
⎪ ⎪D (φ, φ), D (φ, φ), , D (φ, φ) 1N 2N N N ⎝⎭
。。。。。。。
总结:我们可以发现,Ritz 法主要体现在两点:第一,用有限维试解函数空间的基
函数生成真解的一个近似函数,这个试验函数的系数是基本未知量;第二,在求解过程中进行了积分弱化,也就是说,由原来方程的高阶可微要求,通过一次积分和分步积分,得到了弱解形式;第三,使用势能驻值原理。 (2) Galerkin 法
Galerkin 法源于虚位移原理,虚位移原理的本质只是一个数学恒等式。Galerkin 方法也是用S N 代替H, 求U N ∈S N ,使得对任意W N ∈S N 有
D (U N , ω
) =
(f ω,
N N
)
且满足
n
U N =
∑C φ
i
i =1
i
仍然只需要解一个代数方程组,事实上,设:
n
ωN =
∑C φ
i
i =1
n n
0i
i
U N =
∑C
i =1
φi , ωN =
∑C φ
i
i =1
i
则:
N
N
0 =
∑
i , j =1
D (φi , φj ) C j C i -∑C i F (φi )
i =1
N N
i
j
N
=
∑[∑D (φ, φ
i =1
j =1
) C
0j
-∑F (φi )]C i
i =1
由于C i =(i =1, 2, , N ) 是任意的,必须有:
N
∑D (φ, φ
i
j =1
j
) C
0j =F (φi ), (i =1, 2, N )
因此得到与Ritz 方法相同的代数方程组,此时U 就是在S N 取的Galerkin 意义下的广义解。
总结:我们可以发现,Galerkin 法主要体现在两点:第一,用有限维试解函数空间的基函数生成真解的一个近似函数,这个试验函数的系数是基本未知量;第二,由高阶可微变换成弱解积分形式;第三,使用虚位移原理。
(3) 有限元法
有限元方法和Ritz 法、Galerkin 法有什么样的关系呢?
通过上面的分析可以看出,除去由于方法不同可能带来的数值运算差别之外,仅从原理上来说,Ritz 法和Galerkin 法位移的区别就在于前者使用了势能驻值原理,而后者使用了虚位移原理。而有限元法属于Ritz —Galerkin 法的范畴,但又是改进,关键在于有限元法选取分块多项式函数类作为子空间S N 。也就是说有限元方法可以从Ritz 法导出,也可以从Galerkin 法导出(即既可以基于势能驻值原理导出也可以基于虚位移原理导出),但是唯一的区别是它的试解函数空间选取不同。它采用分片插值函数空间,这正好体现了有限元剖分,而上面两种方法虽然也是近似,但是没有分片的概念。而这正是有限元的核心。
具体说,就是先把Ω区间剖分为一系列单元,然后在每个单元内做多项式插值,而使他们在单元的公共点、边、面上满足一定的连续性条件,以保证这种函数组成的有限维空间S N 是H 的子空间。
所以,有限元的基础是变分原理和分片多项式插值。
下面,以二维Poisson 方程的第一边值问题为例来说明有限元方法的这种思想。 Poisson 方程的边值问题:
- u =f (x , y ) i n Ω满足边值条件:
u =α(x , y ) in ∂Ω 第一边值条件
即
⎧- u =f (x , y ) in Ω
⎨
⎩u =0 in ∂Ω
求广义解,就是求 U ∈H
01
(Ω) ,使得:
D (u , ω) =(f , ω), ∀ω∈H
01
(Ω)
其中:D (u , ω) =
⎰⎰(u
Ω
x
w x +u y w y ) dxdy
(f , ω) =
⎰⎰
Ω
f ωd x d y
对Ω进行三角剖分,记三角形单元为e 1, e 2, , e N E ,Ω的内部节点为
P 1, P 2, , P N P ,边界节点为P NP +1, P NP +2, , P NP +M ,所有单元e i 的最大边长是h ,对每一个
单元e=△P i P j P m 做线性插值,使它在三个顶点各取已知值ωi , ωj , ωm , 当节点满足∂Ω上时,给定值为0,这样在Ω上就构成了一个多面体函数ωh (x , y ) ,当(ω1, ω2, , ωN P ) 遍历空间
R N P 所有向量时,函数ωh (x , y ) 就构成了一个N 维函数空间S N 。显然,S N 是一个线性空间,
且是H 01(Ω) (C 10(Ω) 在∙意义下完备化)。
S N 作为H
01
(Ω) 的NP 维线性子空间,它在NP 个基函数
φ1(x , y ), φ2(x , y ), , φN P (x , y ) ,使得ωn (x , y ) 可以表示成:
N P
ωn (x , y ) =
∑ωφ(x , y )
i i
i =1
有限元方法就是求U 1, U 2, , U N P ,使得:
N P
U n (x , y ) =
∑U φ(x , y )
i i
i =
满足:D (u n , w n ) -(f , w n ) =0, ∀w n ∈S N
N P
设:ωn =
N P
∑ωφ(x , y ) ,以U n , ωn 的表达式代入得:
i i
i =1
N P
i
j
∑ωu
i , j =1
D (φi , φj ) -∑ωi (f , φi ) =0
i =1
把D (φi , φj ) ,(f , φi ) 分别写成:
NE
D (φi , φj ) =
∑⎰⎰(∂x
n =1e
n
∂φi ∂φj
∂x
+
∂φi ∂φj ∂y ∂y
) dxdy
N E
(f , φi ) =
∑⎰⎰
n =1e
n
f φi dxdy
简记:D n (φi , φj ) =
⎰⎰(∂x
e n
∂φi ∂φj
∂x
+
∂φi ∂φj ∂y ∂y
) dxdy
(f , φi ) n =
⎰⎰
e n
f φi d x d y
设单元e n = P k P l P m ,因此在单元e n 上只有三个基函数φk , φl , φm 不为零,其它所以基函数在e n 上的值都为零。从而有:
N P N P
i
j
N E i
j
n
∑ωu
i , j =1
D (φi , φj ) =
∑ωu ∑D
i , j =1
n =1
N P
i
j
(φi , φj )
N E
=
∑∑ωu
n =1i , j =1N E
D n (φi , φj )
=
∑∑
N E
ωi u j D n (φi , φj )
n =1i , j =k l , m ,
= ∑
n =1
⎡ωk ⎤⎛D (φk , φk ) ,D (φk , φl ) ,D (φk , φm ) ⎫⎡u i ⎤
⎪⎢⎥⎢⎥
ωl D (φl , φk ) ,D (φl , φl ) ,D (φl , φm ) ⎪⎢u j ⎥ ⎢⎥
D (φ, φ), D (φ, φ) , D (φ, φ) ⎪⎢⎥⎢⎥⎝m k m l m m ⎣ωm ⎦⎭⎣u m ⎦
由于在单元e n 上:
φk =N k , φl =N l , φm =N m 得:
N P
N E
D (u n , ωn ) =
∑ω
i , j =1
i
u j D (φi , φj ) =∑{δ*}e
n =1
T
n
[K ]e {δ}e
n
n
其中:
⎡ωk ⎤⎡u k ⎤
⎢⎥⎢⎥=ωl {δ}e =u l ⎢⎥⎢⎥n ⎢⎢⎣ωm ⎥⎦⎣u m ⎥⎦
{δ*}e
n
虚伪移向量 位移向量 单元上的刚度矩阵
⎛D n (N k , N k ) ,D n (N k , N l ) ,D n (N k , N m ) ⎫
⎪=D n (N l , N k ) ,D n (N l , N l ) ,D n (N l , N m ) ⎪ D (N , N ), D (N , N ), D (N , N ) ⎪
m k n m l n m m ⎝n ⎭
[K ]e
n
N P N P
i
N E i
i n
N E
同理: (f , ωn ) =
∑ω
i =1
(f , φi ) =
∑ω∑(φ, φ)
i =1
n =1
=∑∑
ωi (f , φi ) n
n =1i =k , l , m
N E
=
∑
n =1
⎡ωk ⎤⎢⎥ωl ⎢⎥⎢⎣ωm ⎥⎦
T
⎡(f , φk ) n ⎤
⎢⎥(f , φl ) n ⎢⎥⎢⎣(f , φm ) n ⎥⎦
NE
T
e n
=
∑{δ*}{F }
n =1
e n
{F }e :单元荷载向量
n
NE
T
NE
T
最终:∑{δ*}
n =1
e n
[K ]e {δ}e
n
n
-∑{δ*}
n =1
e n
{F }e
=0
n
意义:
∂φi ∂φj ∂x ∂x
∂φi ∂φj ∂y ∂y
1.⎰⎰(
Ω
+) dxdy
⎰⎰φαfdxdy
Ω
Ω上的计算积分,转换到离散单元上计算,再叠加。 2.稀疏性:不相邻的单元,(φi , φj ) =0 (4) 有限元法的其他解释
进一步上升和总结,可以将Ritz 法、Galerkin 理解为加权残差法的特例。具体见《土
工计算机分析》(龚晓楠 等1999)。 3、补充
实际上,有限元方法的产生应该是一件很自然的事情。因为,对于一个一般的定义于连续区域上(物质的或空间的)微分方程,将其积分是一件很自然的事情,对其利用分部积分也是一件很自然的事情,然后对积分区域剖分也是很自然的,因为根据高等数学的知识我们就已经知道积分对于连续区域的可加性了。这些在探索其理论解的过程中都应该是常见的,唯一的“不自然”的地方就是用插值函数来代替精确函数。因此,有限元方法最先由力学家而不是数学家发现却是“不自然”的。