2015-04-23 有限元方法(Fluids)
流固耦合计算力学
有限元方法（有限元方法（计算流体）计算流体）
刘谋斌
北京大学工学院
第一部分：第一部分：数学基础回顾
第二部分：第二部分：
第三部分：第三部分：
第四部分：第四部分：
：
第六部分：第六部分：一维问题FEM 分析二维问题FEM 分析非定常问题FEM 分析FEM 分析讨论
第五部分：第五部分非线性问题
函数逼近(function approximation)—采用较多数量的简单函数的组合来“近似”代替非常复杂的原函数, 函数逼近—一个复杂的函数，一个复杂的函数，可以通过一系列的基底函数(base function)的组合来“近似”。
（1）基于全域的展开（基于全域的展开（如傅立叶级数展开）如傅立叶级数展开）
Fourier, 1807，“任何周期信号都可用正弦函数的级数表示”；1822，“任何非周期信号都可用正弦函数的积分表示”
（2）基于子域(sub-domain)的分段函数(piecewise functions) 组合
基于子域的函数逼近：基于子域的函数逼近域，形式非常简单；：采用的基函数可以定义在任意相邻子形式非常简单；基函数为线性函数，基函数为线性函数，需要较多分段才能得到较好的逼近精度，得到较好的逼近精度，计算量较大。计算量较大。
求解微分方程的试函数(trial function)法
设方程为：设方程为：L (u ) =f L 为微分算子
∧M
试解为：试解为：u =∑c i φi
i =1
ϕi 为线性无关的基函数；为线性无关的基函数；c i 为待定系数。为待定系数。
基函数选定后，基函数选定后，关键问题是要确定系数
两类方法：两类方法：
（1）变分法
（2）加权余量法
1 数学基础回顾1 数学基础回顾—数学基础回顾—虚功原理与势能最小原理 对描述同一个物理问题的偏微分方程，对描述同一个物理问题的偏微分方程，可以通过对偏微分方程的各种求解方式求解，微分方程的各种求解方式求解，也可以从变分原理出发，化为虚功的形式，或者能量泛函求极值的问题。的问题。 例如，例如，理论力学中质点平衡方程
n ∑F =0i
i =1作用在质点上的合力为0
质点平衡的充要条件：质点平衡的充要条件：当质点
做任意的虚位移时，做任意的虚位移时，所有作用
在质点上的力所作的虚功为0. 虚功原理：虚功原理：∑F ⋅δr =0i i =1n
势能原理：势能原理：若作用在质点上的力是有势的，若作用在质点上的力是有势的，则在稳定平衡时，势能应处于最小，势能应处于最小，即势能u(x, y, z)应取极值（应取极值（极小值）。极小值）。 很多物理定律是通过很多物理定律是通过“对某种积分取极值”的简明的数学原则推求出来的。推求出来的。当被积函数中含有未知函数时，当被积函数中含有未知函数时，有关确定积分极值的问题属于变分计算。极值的问题属于变分计算。
1 数学基础回顾1 数学基础回顾—数学基础回顾—Rayleigh-Rayleigh -Ritz法Ritz 法
Rayleigh-Ritz 方法提供了一个求解边值问题的重要方法—用变分原理求边值问题的近似解析解，由于构造了近似试函数，由于构造了近似试函数，使一些不易求解的微分方程得到了近似解。不易求解的微分方程得到了近似解。
试函数的选择非常重要，试函数的选择非常重要，直接影响精度。直接影响精度。
用变分法求解微分方程，用变分法求解微分方程，相当于把一个含有较高阶导数的微分方程化为一个包含较低阶导数的积分方程，程化为一个包含较低阶导数的积分方程，数学处理上相对方便。数学处理上相对方便。 试函数的选择需要适合整个区域，试函数的选择需要适合整个区域，很难找到适合全区域的试函数，尤其是几何及边界条件复杂。尤其是几何及边界条件复杂。
有限差分启示：有限差分启示元上选择适当的试函数去近似，元上选择适当的试函数去近似：可以把整个求解域分成很多小单元，可以把整个求解域分成很多小单元，即分段逼近。即分段逼近。 有限元，在每一个单 用变分法求解微分方程，用变分法求解微分方程积分形式的变分提法—，但是很多微分方程的变分形式难以找到，但是很多微分方程的变分形式难以找到必须找出相应的泛函，必须找出相应的泛函，把微分方程改为，或者不存在！或者不存在！ 加权余量
1 数学基础回顾1 数学基础回顾—数学基础回顾—加权余量法
引入一组权函数Wi , (i =1, 2,…,n ) ；构造余量ε构造余量ε与全函数的内积，数的内积，并令内积为0并令内积为0（及余量与全函数正交）：及余量与全函数正交）：待解微分方程为：：=0即待解微分方程为∫ε⋅W d Ω=0i
加权余量法就是把域内每一点的余量ε加权余量法就是把域内每一点的余量ε乘以一系列权函数W i , 然后求和，然后求和，要求和为0要求和为0；
相当于强制使微分方程的误差在平均的含义上为0相当于强制使微分方程的误差在平均的含义上为0； n 个全函数，个全函数，n 个方程，个方程，从而可以确定式（从而可以确定式（2.3）2.3）中的n 个待定系数。个待定系数。
权函数有不同的选择方式
GM 选用选择试函数ϕi 本身作为权函数
物理意义：物理意义：力在可能的位移平面上的投影为0，即做功为0 内积为：内积为：=0
在四种加权余量方法(配点法、辽金法) 中，伽辽金法比较好配点法、最小二乘法、。最小二乘法、矩量法及伽
凡是存在变分的问题，凡是存在变分的问题，由伽辽金法求出的近似解与Ritz 变分法近似解完全相同。分法近似解完全相同。
对于不存在变分的问题，对于不存在变分的问题，伽辽金法是一个比较好的选择。伽辽金法是一个比较好的选择。 伽辽金法1915年就已经提出，年就已经提出，由于所选择的试函数必须在全部区域内满足边界条件，在全部区域内满足边界条件，长期难以得到实用，长期难以得到实用，尤其是对于几何形状及边界条件比较复杂的问题。是对于几何形状及边界条件比较复杂的问题。
伽辽金法与子域的离散概念，产生了强大的生命力 伽辽金有限元方法辽金有限元方法
变分原理—将微分方程的求解转换成一个变分形式的积分方程近似解的求解，分方程近似解的求解，要求微分算子对称正定、微分算子对称正定、线性，在流体力学问题中往往难以满足.
伽辽金方法可以视为变分原理的推广—广义变分原理 固体力学—泛函极值—虚功原理
n ∑F =0i
i =1n ∑F ⋅δr =
0i i =1
流体力学—伽辽金方法—虚功率原理
动量方程
虚功率原理：虚功率原理：
2 一维问题2 一维问题FEM一维问题FEM分析FEM 分析
固体力学FEM 分析基本流程
(1)节点编号及单元划分：节点编号及单元划分：就是将原整体结构按几何形状的变化性质划分节点并进行编号，性质划分节点并进行编号，然后将其分解为一个个小的构件（即：单元）单元）--单元特征分析；
(2)建立单元刚度方程：建立单元刚度方程：基于节点位移，基于节点位移，建立每一个单元的节点平衡关系（平衡关系（叫做单元刚度方程）叫做单元刚度方程）--直接法或者虚功原理等；
(3)建立系统刚度方程：建立系统刚度方程：然后将各个单元进行组合和集成，然后将各个单元进行组合和集成，以得到该结构的整体平衡方程（到该结构的整体平衡方程（也叫做整体刚度方程），也叫做整体刚度方程），
(4)边界条件处理及方程求解：边界条件处理及方程求解：按实际情况对方程中一些节点位移和节点力给定相应的值，移和节点力给定相应的值，就可以求解出所有的节点位移和支反力，支反力，
(5)相关变量求解：相关变量求解：最后在得到所有的节点位移后，最后在得到所有的节点位移后，就可以计算每一个单元的其它力学参量（每一个单元的其它力学参量（如应变、如应变、应力）应力）
2 一维问题2 一维问题FEM一维问题FEM分析FEM 分析
流体力学FEM 分析基本流程
(1)节点编号及单元划分：节点编号及单元划分：就是将原复杂区域划分为不重叠、就是将原复杂区域划分为不重叠、无
间隙的有限个特征子域（间隙的有限个特征子域（即：单元）单元）--单元特征分析；(2)建立单元刚度方程：基于节点，基于节点，建立单元有限元方程（建立单元有限元方程（也统称单元刚度方程）称单元刚度方程）--变分原理或伽辽金方法；
(3)建立系统刚度方程：然后将各个单元有限元方程进行组合和集成，集成，以得到整体区域的整体有限元方程（以得到整体区域的整体有限元方程（也统称整体刚度方程）；方程）；
(4)边界条件处理及方程求解：边界条件处理及方程求解：按实际情况处理边界条件，按实际情况处理边界条件，并求解整体方程；解整体方程；(5)相关变量求解：相关变量求解：
11
2 一维问题2 一维问题FEM一维问题FEM分析FEM 分析
例1：求解下述一维问题的伽辽金有限元解
基本方程边界条件
u |x =0=0
du
精确解
dx
|x =1=0 求解流程
（1）单元划分—连续体离散为有限个特征子域（连续体离散为有限个特征子域（单元），单元），对），对
节点及单元编号
12
2 一维问题2 一维问题FEM一维问题FEM分析FEM 分析
节点在单元中的局部编号与整体编号的关系
布尔矩阵
单元分析—插值函数（插值函数（与维数、与维数、单元节点个
数及分布、数及分布、节点类型相关）节点类型相关）试函数
u
(e )
=a 1+a 2x
(e )
1
a 1, a 2为待定系数
局部坐标中对于单元中每个节点有
x =0:u (e ) x =h :u 2
a u =1=a 1
1(e ) (e )
=a 1+a 2h
a 2=(u 2−u 1)
h 13
(e ) 1
2 一维问题2 一维问题FEM一维问题FEM分析FEM 分析
回代可得即
u
(e )
=u
(e )
1
x (e ) x (1−) +u 2h h
下标重复表示项数叠加
u
(e )
=u
(e ) N
⋅ϕ
(e ) N
ϕ
(e ) N
（局部）局部）基函数、基函数、插值函数、插值函数、形状函数
与固体力学有限元分析中杆单元类似！
14
2 一维问题2 一维问题FEM一维问题FEM分析FEM 分析
u
(e )
N
=∆⋅u i
(e )
N 2×1
(e ) Ni
u =∑u
e =1
E
(e )
=∑[ϕ]
(e )
N 1×2
⋅[u ]
=∑[ϕ]
(e )
N 1×2
⋅[∆]
(e )
Ni 2×4
⋅[u i ]4×1
=[ϕi ]1×4⋅[u i ]4×1=[u i ]1×4⋅[ϕi ]4×1
整体形状函数
15
2 一维问题2 一维问题FEM一维问题FEM分析FEM 分析
整体形状函数的性质
近似函数分布图
16
2 一维问题2 一维问题FEM一维问题FEM分析FEM 分析
（2）伽辽金有限元方程建立
代入待求微分方程
余量为
构造与整体权（构造与整体权（形状）形状）函数的内积，函数的内积，令内积为0，可得
分部积分，分部积分，可得
2 一维问题2 一维问题FEM一维问题FEM分析FEM 分析
整理、整理、改写为
一维情况，一维情况，边界为两个端点
或表示为
整体有限元方程
18
也可对单元写出伽辽金有限元方程
代入待求微分方程，代入待求微分方程，得单元余量
构造内积，构造内积，令内积为
分部积分，分部积分，整理可得
或写出单元有限元方程
19
单元有限元方程
整体有限元方程
通过布尔矩阵、通过布尔矩阵、局部插值函数及整体插值函数展示了整体与局部有限元
方程中各系数的关系—可组装、可组装、合成
20
（3）单元（单元（系数）系数）计算
21
（3）单元（单元（系数）系数）计算
取4个单元，个单元，
h=1/4
同理系数B 、C 、F 、G 可以求出，可以求出，A 、B 、C 各个单元相同
F= ，各个单元不同, x m 为总体坐标G 表示边界条件，表示边界条件，x =1.0出存在自然边界条件各单元线性代数方程组
22
（4）总体合成
把各单元的系数矩阵进行整体合成，整体合成，各单元系数矩阵在整
体系统矩阵中的位置由布尔矩阵决定布尔矩阵的作用：布尔矩阵的作用：使相邻单元的贡献在共同的节点上叠加。
[A ij ]5×5=∆⋅A
(1)
iN (1)NM
⋅∆+∆
(1)Mj (2)iN
⋅A
(2)NM
⋅∆
(2)Mj
+∆⋅A
(3)iN (3)NM
⋅∆
(3)Mj
+∆
(4)iN
⋅A
(4)NM
⋅∆
(4)Mj
23
（5）边界处理与方程求解
边界处理
对称、对称、三对角矩阵
24
（5）边界处理与方程求解
二次插值收敛性远强于线性插值，强于线性插值，收敛速度按h 2k 增长, k 为插值函数阶次
25
例2：求解下述一维问题的伽辽金有限元解
基本方程边界条件
精确解
求解流程
（1）单元划分—简单分为2个单元，个单元，单元长度h
=0.5（2）插值函数—选择两节点线性插值函数，选择两节点线性插值函数，同上
26
（3）建立伽辽金有限元方程组
近似函数代人微分方程，微分方程，产生残差（余量）余量） 用插值函数对残差加权，差加权，并令内积为0
分部积分，分部积分，整理、整理、改写，改写，可得伽辽金有限元方程 该伽辽金有限元方程建立过程对单元和整体都适用
（4）单元分析
及同样方法可以求出
27
（5）组装合成
（6）边界处理与方程求解（边界处理与方程求解（略）
28
例3：求解下述二维问题的伽辽金有限元解
基本方程边界条件
求解流程
（1）单元划分—把二维计算域划分为有限个特征单元（2）插值函数—选择合适的插值函数
（3）伽辽金有限元表达式—代人微分方程、微分方程、产生残差，以插
值函数为权函数，值函数为权函数，构造内积，令内积为0
29
（3）伽辽金有限元表达式—代人微分方程、微分方程、产生残差，以插
值函数为权函数，值函数为权函数，构造内积，令内积为
分部积分、分部积分、改写整理
边界项
二维问题边界为一曲线，沿该边界曲线的插值函数可借
鉴一维插值函数
l : 边界线长度s : 沿边界项距离
伽辽金有限元方程
30
3 二维问题3 二维问题FEM二维问题FEM分析FEM 分析
伽辽金有限元方程
（4）单元分析
（5）整体合成
（6）边界处理及方程求解
31
3 二维问题3 二维问题FEM二维问题FEM分析FEM 分析
二维有限元的插值函数—类似固体力学有限元中二维形状函数 3节点三角形单元—形式
3个节点
确定系数
与固体力学中3节点三角形单元形函数一致
系数a i , b i , c i 仅取决于单元的几何形状32
二维有限元的插值函数—类似固体力学有限元中二维形状函数 3节点三角形单元—性质
33
多节点三角形单元—高阶插值高阶插值函数需要更多节点
往往采用面积比形成的自然坐标
二次插值节点分布二维完全多项式排列
34
4节点矩形单元
质心坐标
面积坐标
双线性插值矩形 正方形
35
4节点矩形单元
双二次插值
9节点矩形单元8节点矩形单元4端点、端点、4中点加质心4端点、端点、4中点完全二次多项式非完全二次多项式36
例4：求解下述问题的伽辽金有限元解基本方程
边界条件
方法一：方法一：插值函数与时
空同时关联（空同时关联（全离散）全离散）
方法二：方法二：插值函数仅与
空同时关联（空同时关联（半离散）半离散）
37
求解策略（半离散法）：离散法）：在空间上的离散完全依照定常问题伽辽金方法，伽辽金方法，从而把偏微分方程转化为函数节点值对时间t 的常微分方程，常微分方程，时间上可以用差分法求解（时间上可以用差分法求解（隐式、显式）。 求解流程
（1）单元划分—把计算域划分为有限个特征单元
（2）插值函数—选择合适的插值函数
（3）伽辽金有限元表达式—代人微分方程、微分方程、产生残差，以插
值函数为权函数，值函数为权函数，构造内积，令内积为
0分部积分、积分、改写整理
38
分部积分、积分、改写整理解耦时间
39
可写成
其中
（4）单元分析
（5）整体合成：整体合成：整体系统有限元方程
40
4 非定常问题4 非定常问题FEM非定常问题FEM分析FEM 分析
（5）整体合成：整体合成：整体系统有限元方程
（6）时间项处理：处理：显式或隐式
假设采用一阶向前差分（差分（欧拉格式），有），有
41
例5：求解下述问题的伽辽金有限元解
基本方程
伽辽金有限元积分方程
分部积分
42
整理重写
单元有限元方程
系数A NMQ ⋅u M ⋅u Q +A NM ⋅u M =F N
43
5 非线性问题有限元分析5 非线性问题有限元分析
整体有限元方程
整体有限元方程由于第一项的存在，的存在，是非线性的，是非线性的，对非线性项统称采用Newton-Raphson 法求解，求解，其基本思想是设法把非线性方程转化为某种线性方程求解已知非线性方程的一个根为
则Newton-Raphson 法的迭代公式为：式为：
为迭代的新的近似根！
44
多场流体FEM 求解（求解（V 、P 、Rho ）
45
流体力学FEM 与固体力学FEM 的区别
（1）数学基础—FEM for fluids一般基于Galerkin 方法，方法，FEM
for solids一般基于变分原理（基于变分原理（泛函极值）泛函极值）
（2）物理概念—FEM for solids反映虚功原理和最小势能原
理，FEM for fluids对应虚功率原理
（3）基本变量—FEM for solids（弹性力学）性力学）基于位移、基于位移、应变、应变、
应力三大变量，变量，对应几何方程、对应几何方程、物理方程及
平衡方程三大方程，方程，FEM for fluids基于基本
流体动力学变量（力学变量（速度、压力、密度、能
量）, 对应连续方程、对应连续方程、动量方程及能量方程。量方程及能量方程。
（4）空间描述—FEM for solids位移为基本变量，位移为基本变量，网格节点随
时间运动、调整，固FEM 网格为拉格朗日网
格，FEM for fluids网格一般为Euler 网格
46
FDM 、FVM 及FEM 的进一步讨论
有限差分：有限差分：
有限体积：有限体积：
有限元：将计算区域划分为差分网格，将偏微分方程的微商（导数）导数）用差商代替，得到各个离散点上含有限个未知数的差分（知数的差分（代数）方程组，方程组，是偏微分方程的数值近似解法--应用最早、最经典的CFD 方法。方法。将计算区域划分为一系列控制体积，体积，对待求解偏微分方程对各个控制体积积分得出离散方程组—需要对边界上的被求函数与导数的分布进行某种形式的假设，具有守恒性好、物理概念明确、明确、计算量相对较小，计算量相对较小，是目前CFD 应用最广的方法。的方法。了有限差分法中离散处理的内核，又采用了能量
极值原理与分段逼近思想，选择逼近试函数对区域进行积分的特点，行积分的特点，是目前CSD 应用最广的方法。的方法。
47吸收