2-5有限元法在流体力学中的应用

第五章有限元法在流体力学中的应用
本章介绍有限元法在求解理想流体在粘性流体运动中的应用。

讨论了绕圆柱体、翼型和轴对称物体的势流，分析了求解粘性流动的流函数—涡度法流函数法和速度—压力法，同时导出粘性不可压流体的虚功原理。

§1 不可压无粘流动
真实流体是有粘性和可压缩的，理想不可压流体模型使数学问题简化，又能较好地反映许多流动现象。

1. 圆柱绕流
本节详细讨论有限无法的解题步骤。

考虑两平板间的圆柱绕流．如图5—1所示。

为了减小计算工作量，根据流动的对称性可取左上方的l／4流动区域作为计算区域。

选用流函数方法，则流函数 应满足以下Laplace方程和边界条件
203
204
22220(,)0(,)2(,)(,)0(,)x y x y x y aec x y bd y x y ab x y cd n ψψ
ψψ
⎧∂∂+=-∈Ω⎪∂∂⎪⎪-----∈⎧⎪⎪=-----∈⎨⎨⎪⎪-----∈⎩⎪⎪∂=-----∈⎪∂⎩
流线流线流线
流线 (5-1)
将计算区域划分成10个三角形单元。

单元序号、总体结点号和局部结点号都按规律编排．如图5—2所示。

从剖分图上所表示的总体结点号与单元结点号的关系，可以建立联缀表于下
表5-1
各结点的坐标值可在图5—2上读出。

如果要输入计算机运算必须列表。

本质边
205
界结点号与该点的流函数值列于下表
表5-2
选用平面线性三角形元素，插值函数为(3—15)式。

对二维Laplace 方程进行元素分析，得到了单元系数矩阵计算公式(3—19)和输入向量计算公式(3—20)。

现在对全部元素逐个计算系数矩阵。

例如元素1，其结点坐标为1x =0, 1y =2; 2x =0, 2y =1; 3x =2.5, 3y =2. 由(3—15)式可得
132 2.5a x x =-=; 213 2.5a x x =-=- 3210a x x =-=,
1231b y y =-=-; 2310b y y =-=;
3121b y y =-=; 0 1.25A =
从(3—19)式可计算出1K
1 1.45 1.250.21.2500.2K ⎛⎫
⎪
⎪
=
⎪ ⎪
⎝
⎭
--对称
依次可计算出全部子矩阵
20.20.201.45 1.251.25K ⎛⎫
⎪
⎪
= ⎪ ⎪
⎝
⎭
--
30.200.21.25 1.251.45K ⎛⎫
⎪
⎪
= ⎪ ⎪
⎝⎭
--
206
4 1.2
5 1.2501.450.20.2K ⎛⎫
⎪
⎪
= ⎪ ⎪
⎝⎭
--
50.50.5000.50.5K ⎛⎫
⎪
⎪
= ⎪ ⎪
⎝
⎭--
60.500.50.50.51K ⎛⎫
⎪
⎪
= ⎪ ⎪
⎝
⎭
--
70.50.5010.50.5K ⎛⎫
⎪
⎪
= ⎪ ⎪
⎝
⎭--
80.500.50.50.51K ⎛⎫
⎪
⎪
= ⎪ ⎪
⎝⎭--
90.500.50.50.51K ⎛⎫
⎪
⎪
= ⎪ ⎪
⎝⎭
--
1010.50.50.500.5K ⎛⎫
⎪
⎪
= ⎪ ⎪
⎝
⎭
--
根据联缀表把元素矩阵组合成总体系数矩阵 A=
1.450.20 1.25
2.4500 1.25 1.01000.50.52.90.400 1.254.9100 1.750.54.01000.52000.51.450.201.9501.5--⎡⎤
⎢⎥--⎢⎥⎢⎥--⎢⎥
--⎢⎥
⎢⎥---⎢⎥--⎢
⎥⎢⎥-⎢⎥
-⎢
⎥⎢⎥⎢⎥⎢⎥⎣⎦0对称
207
矩阵中零元素没有一一写出，下三角部分与上三角部分对称。

从(3—20)式计算元素输入向量，由于流函数满足齐次的自然边界条件
0n q n
ψ
∂==∂,所以输入向量为零，总体输入向量也为零，这样就得了总体有限元方程.
N A B ψ=
式中:
[]1210,,,T
N ψψψψ=
[]0,0,,0T
B =
用缩减方程的重新编号修正方法施加边界条件，本质边界结点的函数值是已知的。

把它们代入方程，修正右端项，再减去相应的方程，整理得
5674.910 2.914130121ψψψ-⎡⎤⎛⎫⎡⎤
⎪⎢⎥⎢⎥--= ⎪⎢⎥⎢⎥ ⎪⎢⎥⎢⎥-⎝⎭⎣⎦
⎣⎦ 解方程得到
5ψ＝0.845，6ψ=1.241，7ψ=1.121
这样求出了全部结点上的流函数。

为了求出每个单元形心处的速度，可以由单元的流函数近似表达式求导计算。

对元素e 来 说，有
T I
y
ψψ=Φ
[]112320
31,,2T I T I y u a a a A y A ψψψψψψ⎡⎤
∂⎢⎥==Φ==⎢⎥
∂⎢⎥⎣⎦
[]1231,,2T I
I T I x b b b B x A
ψνψψψ∂-=-
=-Φ==-∂ 例如单元1ψ=2, 2ψ=1, 3ψ=3，这样计算得到的速度为u=1,ν＝0。

二维绕圆柱流动还可以用势函数求解，
则定解问题可写成 22
220(,)001x y x
y cd aec bd n ab n
ωω
ωωω
⎧∂∂+=-∈Ω
⎪∂∂⎪⎪=------⎪⎨∂=-----⎪∂⎪⎪∂=------⎪∂⎩边界上
边界及上边界上
208
ω表示势函数，
为了使数值解唯一必须在部分边界上给定本质边界条件。

势函数边界同样标记在图5—l 上。

势因数满足Laplace 方程和相应的边界条件，与流函数不同仅在于有非齐次的自然边界条件。

采用与流函数方法完全一样的网格划分，可知计算得到的单元系数矩阵是完全一样的，总体矩阵也是完全一样的。

元素1和4具有非齐次自然边界条件．应该用(3—20)式计算输入向量。

元素l, 111,,022T
P ⎡⎤
=⎢⎥⎣⎦。

元素4, 411,,022T
P ⎡⎤=⎢⎥⎣⎦。

总体合 成得到1
1001000002
2T
B ⎡⎤=⎢⎥⎣⎦，这样就得到方程组
N A B ω=
巳知37100ωωω===，消去相应的三个方程得到一个7×7的 代数方程组，解得
1233.787, 1.204,0ωωω=== 4563.841, 1.261,0.616ωωω===
759100, 3.827, 1.491,0.ωωωω====
单元形心处的速度可以用下列公式计算
T I T I x
u B x
ωωω∂==Φ=∂ T I T I
y A y
ωνωω∂=
=Φ=∂ 式中I ω是单元的结点势函数向量[]123,,ωωω。

对于元素1来说，
1233.787, 3.841, 1.204ωωω===,这样计算得到u=1.033，v=-0.05。

这结果与流函
数方法得到的结果近似相等。

如果加密网格，就可以得到更好的结果。

2. 升力问题
考虑图5—3(a)所示的机翼绕流。

均匀来流u ∞平行于x 轴，机翼边界为1Γ，
209
后缘尖点为T ，流场外边界1Γ取在离机翼足够远处。

流函数ψ 满足以下方程和边界条件。

201020310a hu a yu a b ψψψψψ∞∞⎧∇=-----Ω⎪
=------Γ⎪⎪
=+---Γ⎨⎪=+---Γ⎪
⎪=-------Γ⎩在内在上在上在上
在上
(5-3) 其中a,b 是特定系数，h 是上下边界之间的距离。

机翼绕流的后驻点应位于后缘尖点处，在后缘T 点满足Kutta 条件
0u y
ψ
∂=
=∂；0v x ψ∂=-=∂； （5-4）
210
由于方程和边界条件是线性的,可用叠加原理求解，令
012a b ψ
ψψψ=++ （5-5）
其中0ψ，1ψ和2ψ：分别是下列问题的解
20010
000yu ψψψ∞⎧∇=-----Ω⎪
=-------Γ⎨⎪=------Γ⎩在内在上在上
21111
0001ψψψ⎧∇=-----Ω⎪
=-------Γ⎨⎪=--------Γ⎩在内
在上在上
22212
0010ψψψ⎧∇=-----Ω⎪
=-------Γ⎨⎪=--------Γ⎩在内
在上在上
用有限元方法分别解以上三个问题，得到各结点的0ψ、1ψ和2ψ，代入(5—5)式得到叠加解。

显然它满足问题(5—3)的全部方程和边界条件，特定常数a,b 可利用Kutta 条件(5—4)定出。

首先由流函数0ψ、1ψ和2ψ分别求出各个结点上的速度0,)u v 0(，1,)u v 1(和2,)u v 2(，然后在后缘点T 处利用Kutta 条件，应有
012
012u u au bu v v av bv =++⎧⎨=++⎩
解之可得到a 和b 。

图5—3(b)上给出了NACA4412具型以8 攻角置于均匀流场中所引起的流动图案，计算中采用了三角形单元。

与无升力体绕流一样，机具绕流也可以采用速度势函数求解. 3．轴对称问题
考虑圆管内绕轴对称物体的无旋流动，如图5—4(a)所示。

采用柱坐标系(r ，θ,z)，其势函数满足Laplace 方程。

211
12
2222221210S S r
r r z S q S n ϕϕϕϕϕϕ⎧∂∂∂⎪+--+=----Ω∂∂∂⎪⎪=--------------⎨⎪∂⎪=-------------⎪∂⎩在内
在上在上 (5-6)
写出与微分问题相应的伽辽金积分表达
2222221d r r r z
ϕϕϕδϕΩ
∂∂∂++Ω∂∂∂⎰（） =2S q ds n
ϕ
δϕ∂-∂⎰（） 分部积分上式的左边并整理得到弱解积分形式
2)rdrdz r r z z
ϕδϕϕδϕ
π∂∂∂∂+∂∂∂∂⎰⎰（
=0
2L
q rdl π
δϕ⎰
式中L 是元素的边长，L 绕轴旋转一周形成元素的边界面。

采用图5—4(b)所示轴对称的环形线性元素，它是将平面线性三角形元素绕对称轴旋转一周形成的环形体。

采用斜坐标系，那么插值函数可写成
{}123,,T ξξξΦ=
212
元素结点上势函数向量为
{}123(),,I T ϕϕϕϕ=
则逼近函数为
112233T I ϕϕξϕξϕξϕ=Φ=++
总体坐标和斜坐标系的关系为
T I
T I r r z z
⎧=Φ⎪⎨=Φ⎪⎩ 式中{}123(),,I T
r
r r r =。

{}123(),,I T z z z z =，是元素结点总体坐标向量。

将逼近函数表达式代入伽辽金公式，推导出元素有限元方程
I K P Φ=
式中影响系数矩阵和输入向量分别为
K=2)T T
r z rdrdz πΦΦ+ΦΦ⎰⎰r z （
P=0
2L
q dl πΦ⎰r
求出插值函数向量的偏导数r Φ和z Φ，代入上式得影响系数矩阵
K=22111212
131********
23230
2233()6a b a a b b a a b b r r r a b a a b b A a b π⎛⎫
+++++ ⎪
++
⎪ ⎪+⎝
⎭
（5-7） 式中 i k j a r r =-；i j k b z z =- i ＝1，2，3时J=2，3，1;k ＝3．1，2。

012212A b a b a =- ， 0A 三角形元素面积。

假设元素的“l 一2”边落在自然边界上且q 为常数，则可得转入向量计算公式
1212122230r r qL P r r π+⎡⎤
⎢⎥=+⎢
⎥⎢⎥⎣⎦
（5-8） 式中 12L ：是“l 一2”边的边长。

213
计算了各元素的K 和P ，然后总体合成，代入本质边界条件就可以解总体方程。

为了计算其它物理量，下面给出了相应的公式。

元素形心处的速度：
1122330
1122330
1
()21
()2z r U b b b A U a a a A ϕϕϕϕϕϕ=++=
++ （5-9）
附加质量m ：
m=11
2
21
()N
i p p ρ
ϕϕ
=+∑ (5-10)
式中，i=1,2,....N. N 是物体表面上所划的单元数。

p 是输入向量P 在元素结点上的值。

在文献[4]中，以圆球作为例子计算了三种状态。

球在无阻空间中运动，计算的附加质量系数λ＝0.4671，理论值是0.5。

计算值小于理论值是符合第二章2节例4中证明的附加质量极大值原理的.在圆管中球作匀速运动，计算的结果与T ．J ．Chung 在参考书(9)中给出的结果比较，虽然我们采用了较少的结点，但达到了相同的精度。

Chung 用流函数方法，采用轴对称四边形单元计算。

由于流函数满足Stokes 方程，是非自伴的，这样行成的影响系数矩阵是非对称的．不仅计算麻烦．而且不能利用半带宽存储。

四边形单元的短阵元素计算须用Guass 数值积分,计算量大且有误差。

而采用势函数方法和三角形元素恰好克服了以上两个缺点。

第三种状态计算了圆球在半盲管(一端封死)中的运动。

附加质量系数
λ＝0.897，这等于无限空间中附加质量系数的1.6倍。

这使我们想到，在计算水下管中发射弹道问题时，应考虑物体在管道中的附加质量系数。

轴对称不可压无粘流动也存在看流函数提法，流函数ψ应满足以Stokes 方程，而不是Laplace 方程。

2222
10r r r z ψψψ
∂∂∂-+=∂∂∂ （5-11） 应特别注意的是,Stokes 算子是非自伴的。

为了写出相应的迦辽金弱解积分表达式,先可将方程改写
22221120r r r z r r
ψψψψ
∂∂∂∂++-=∂∂∂∂
214
其伽辽金积分表达为
2222112)0d d r r r z r r ψψψψδψδψΩΩ∂∂∂∂++Ω-Ω=∂∂∂∂⎰⎰（）（
将上式第一项分部积分，并代入自然边界条件，得到
12(
)4()rdrdz drdz r r z z r r ψδψψδψψ
ππδψ∂∂∂∂∂++∂∂∂∂∂⎰⎰⎰⎰
=22S rdl n ψ
πδψ∂∂⎰ 假设近似解可表示成
T I ψψ=Φ ,I T δψδψ=Φ
那么
T I r r ψψ∂=Φ∂ ，T I
z z ψψ∂=Φ∂ ,I T r r δψδψ∂=Φ∂ ，,I T z z
δψ
δψ∂=Φ∂ 将以上各式代入弱积分表达式，得到单元有限元方程
I K P ψ=
式中影响系数矩阵为
K ＝2)4T T T r z r rdrdz drdz ππΦΦ+ΦΦ+ΦΦ⎰⎰⎰⎰r z （
式中第二项是非对称项，它使得K 成为非对称矩阵。

输入向量为
P=0
2L
r q dl π
Φ⎰
如果采用前面已用到的三角形洄旋环状体元素，则是K 和P 可以导出，得到相应的计算公式。

看来流函数方法不及势函数方法简便。

§2 不可压粘性流动
不可压粘性流体运动由速度散度为零的连续方程及Navier —Stokes 方程描述。

二维问题，引进流函数可导出流函数—涡量方程和四阶的流函数方程。

粘性流动中存在粘性应力，固壁边界上必须满足无滑条件，使得流体的运动一股是有旋的.
215
1．流函数问量法
以流函数ψ和涡量ω。

表示的粘性不可压流动方程和自然边界条件是
()
xx yy i y x x y xx yy ψψωωψωψωυωω+=-⎧⎪⎨+-=+⎪⎩ （5-12） s v S n g S n ψω
ψ
ω∂⎧=-------⎪⎪∂⎨
∂⎪=-------⎪∂⎩
在上在上
（5-13）
边界上还应给出本质边界条件，即ψ和ω的函数值。

根据具体边界不难绘出流函数的边界值。

而固壁上的涡量，则要通过区域内的流函数值，由涡量边界条件来确定。

流函数ψ、涡量ω和速度的关系为
,,v u
u v y x x y
ψψω∂∂∂∂=
=-=-∂∂∂∂ 分别写出流函数方程和涡量方程的伽辽金积分表达
2()d ψωδψΩ
∇+Ω⎰⎰
2()0i y x x y d ωψωψωυωδωΩ
+--∇Ω=⎰⎰
对以上两式中Laplace 算子进行分部积分并代入自然边界条件：得弱解
积分形式
(
)s S d v ds x x y zy
ψψδψψδψωδψδψΩ
∂∂∂∂+-Ω=-∂∂∂∂⎰⎰
⎰
()()i y x x y d x x y y ωδωωδωωψωψωδωυΩ⎡⎤
∂∂∂∂+-++Ω⎢⎥∂∂∂∂⎣⎦⎰⎰ =S vg ds ω
δω⎰
将求解区域剖分，单元中流函数ψ、涡量ω。

选用相同的插值展开形式，则有
T I ψψ=Φ ，T I ωω=Φ
216
其中 {}12(),(),()I I t t t ψψψψ= 。

{}12(),(),()I I t t t ωωωω= 是结点未知数向量。

()i t ψ和()i t ω分别是单元中i 结点在t 时刻流函数和涡量值。

这样就有
,I T δψδψ=Φ，,I T δωδω=Φ，T I i ωω=Φ
将流函数和涡量的插值表达式及其变分代入弱解积分形式，得
{}
,()0I T T T I T I x y s S d d v ds ψ
δψψωΩ
Ω
⎡⎤ΦΦ+ΦΦΩ-ΦΦΩ+Φ=⎣
⎦
⎰⎰⎰⎰⎰x y
{}
,(()()I T
T
I
T I
T T I T I T T I
y
x
x y x x y y S d d d vg ds
ω
δωωψψωυωΩ
Ω
Ω
ΦΦΩ+ΦΦΦ
-ΦΦΩ+ΦΦ+ΦΦΩ-Φ⎰⎰⎰⎰⎰⎰⎰=0
整理得到单元的有限元方程
I I
K M P ψψω-=
I I I M A K P ωωωυω++=
流函数的单元有限元方程是代数方程，而涡量的有限元方程是常微分方程组，必须在时间步进中求解。

式中各矩阵分别为
质量矩阵 T M d Ω
=ΦΦΩ⎰⎰
损耗矩阵 ()T T
x x y y K d Ω
=ΦΦ+ΦΦΩ⎰⎰
对流矩阵 ()T I T T I T y x x y A d ω
ψψ=ΦΦΦ-ΦΦΩ⎰⎰
输入向量为 s S P v d s
ψ
ψ=-Φ⎰
S P vg ds ω
ω=Φ⎰
I ψ是t 的函数，对流矩阵A 是与时间t 有关的非对称矩阵，因此每一时步
必须重新计算。

将所有单元进行总体合成，得到总体有限元方程
ϕϕ=+ΩK P M （5-14）
.
M A vK B ωΩ+Ω+Ω= （5-15）
式中ψ和Ω是总体结点未知数向量。

总体有限元方程是代数方程和非线性常微分方程的联立方程组。

将涡量
217
力程用显式差分貉式离散，可用交替迭代的方法求解耦合方程组。

用流函数涡量法，可以计算低雷诺数时绕障碍，如圆柱、圆球的流动，还可以模拟涡街的形成过程。

2．流函数方法
二维不可压粘性流动可用四阶的流函数方程描述
24()0d
dt
ψυψ∇-∇= （5-16） 写出相应的伽辽金积分表达
24()0d d dt ψυψδψΩ⎡⎤∇-∇Ω=⎢⎥⎣⎦
⎰⎰ 利用分部积分
2()d d dt ψδψΩ⎡⎤
∇Ω⎢⎥⎣⎦
⎰⎰ =()()()S d d d ds d dt n dt x x dt y y ψ
ψδψψδψδψΩ⎡⎤∂∂∂∂∂-+Ω⎢⎥∂∂∂∂∂⎣
⎦⎰
⎰⎰ 4d ψδψΩ⎡⎤∇Ω⎣⎦⎰⎰
=2222222
2222d x x x y x y y y ψδψψδψψδψ
Ω⎡⎤
∂∂∂∂∂∂++Ω⎢⎥∂∂∂∂∂∂∂∂⎣⎦
⎰⎰ -2222()()()()s s ds ds x n x y n y x n y n ψδψψδψψψδψ⎡⎤
⎡⎤∂∂∂∂∂∂∂∂∂∂+++⎢⎥⎢⎥∂∂∂∂∂∂∂∂∂∂⎣⎦⎣⎦
⎰⎰ 可得弱解积分表达式
2222222222()()(2)d d d dt x x dt y y x x x y x y y y ψδψψδψψδψψδψψδψυΩ⎡⎤
∂∂∂∂∂∂∂∂∂∂++++Ω⎢⎥∂∂∂∂∂∂∂∂∂∂∂∂⎣⎦
⎰⎰ =
2222()()()()()S d ds dt n x n x y n y x n y n ψψδψψδψψψδψυυδψ⎧⎫⎡⎤⎡⎤∂∂∂∂∂∂∂∂∂∂∂⎪⎪++-+⎨⎬⎢⎥⎢⎥∂∂∂∂∂∂∂∂∂∂∂⎪⎪⎣⎦⎣⎦⎩⎭
⎰ 假定在全部边界上给定本质边界条件，对于四阶算子，即给定
s ψψ=，
218
s S
v n
ψ∂=-∂，
0S
S
ψ∂=∂
那么在S 上有0δψ=及
x δψ∂∂和y
δψ
∂∂等于零，则伽辽金弱积分形式可简化为 2222222222()()(2)0d d d dt x x
dt y y x x x y x y y y ψδψψδψψδψψδψψδψυΩ⎡⎤
∂∂∂∂∂∂∂∂∂∂++++Ω=⎢⎥∂∂∂∂∂∂∂∂∂∂∂∂⎣⎦⎰⎰
设单元的近似函数为
T I ψψ=Φ 式中
{}12,,I φφφΦ=
{}12(),(),()I I t t t ψψψψ=
将近似函数代入弱解表达式，可以推导出单元有限元方程。

在推导中首先注意到全导数的运算
222
()()d dt x t x y x x x y ψψψψψψ
∂∂∂∂∂∂∂=+-∂∂∂∂∂∂∂∂ 222
()()d dt y t y y x y x y ψψψψψψ
∂∂∂∂∂∂∂=+-
∂∂∂∂∂∂∂∂ 然后把近似解代入积分表达式中．整理得到
I I I I +A()0M K ψψψψ+= （5-17）
式中
()T T
s x y y M d Ω
=ΦΦ+ΦΦΩ⎰⎰
{}I A()(()T I T T I T T I T T I T x y xx x xy y y xy x yy d ψψψψψΩ
=ΦΦΦ-ΦΦ+ΦΦΦ-ΦΦΩ⎰⎰
(2)T T T
xx xx xy yy K d υΩ
=ΦΦ+ΦΦ+ΦΦΩ⎰⎰xy yy
矩阵I A()ψ包含着未知向量 I ψ，这是对流项的非线性影响产生的。

如果写成下标表示的方程组，则可以看出方程的第二项包含ψ的二次项，所以得到的是
219
非线性的常微分方程组.
由于流函数方程是四阶偏微分方程，因此函数ψ及其一阶导数x ψ、y ψ，是本质变量，应取为结点未知数并保持其在元素内和元素之间的连续性。

为此，应选用函数及一阶导数都连续的Hermite 插值函数．要构造相容的Hermite 插值函数有时是困难的。

另外单元近似函数要用高阶插值，元案自由度大，这使得总体有限元方程的系数矩阵阶次很高．要占用大量的计算机内存.这种方法只具有一个未知函数，所以计算程序简单，也比较省时。

用流函数方法计算圆柱绕流，得到的表面压力分布曲线与有限差分解比较，符合得很好。

§3. 流体力学虚功率原理及速度压力法 不可压粘性流动基本方程如下
连续方程
0j j
v x ∂=∂ (5-18)
动量方程 jk k
k j
dv b dt x σρρ∂=+∂ k=1,2 (5-19) 本构方程 (
)j k
jk jk jk jk k
j
v v p p x x σδτδμ∂∂=-+=-++
∂∂ (5-20)
方程中重复指标暗指叠加运算。

本构方程给出了应力和应变率之间的关系。

其中应力张量jk σ是二阶张量，表示作用在流体微团小六面体的j 面上k 方向的应力分量。

由于流体微团动量矩平衡的原因．它是对称张量，即有
jk kj σσ=
jk τ是粘性应力张量。

10jk j k
j k δ---=⎧=⎨
---≠⎩
称为Kronecker δ k b 是质量力分
量。

流动区域边界可分成二类，一类是指定速度的，另一类是指定应力。

在指定速度的边界v S 上
220
j j v v = k=1,2
如图5—5所示．如果边界的单位外法向矢量为n ，其方向余弦cos(,)nj j a n x =
j=1,2. 那么其法向速度应满足 n j n j j n n v v a v a v === (5-21) 如固壁边界，流体质点粘附在固壁上，流体速度等于固壁运动速度，满足无滑条件。

若固壁静止，则其值为零。

在指定应力条件的边界P S ，上，有
n k n k
p p = nk p 是作用于边界上，法向为n 的单位微元面积上的应力矢量n P 在k 方向的分量．
见图5—5。

设液体微元j 面上的应力矢量为
j kj k i σσ=
式中k i 表示直角坐标系中的两个单位矢量。

那么n P 等于二个应力矢量在n 方向的投影之和，即
n j n j
n j j P a a i σσ== 那么应力边界条件可写成
n k n j
j k n k
P a P σ== （5-22）
221
如自由液面上，nk P 是法向为n 的自由面上结定的作用于液体的风应力，多大气静止常压时，其值为零.流体入口或出口边界，可给出速度或应力值。

不同液体分界面上，一船给出速度，压力连续性条件。

把本构方程(5—20)式代入动量方程(5—19)式，并利用连续方程，就可以导出Navier —Stokes 方程。

下面推导虚功率原理。

以压力的变分P δ为权函数．写出连续方程的伽辽金积分公式
0j j
v Pd x δΩ
∂Ω=∂⎰⎰
分部积分．得到弱解积分公式
()
j
n S j
P v d v Pds x υ
δδΩ
∂Ω=∂⎰⎰
⎰ （5-23）
以速度的变分为权函数，写出动量方程的伽辽金积分公式
()0jk k k k j dv b v d dt x σρρδΩ⎡⎤
∂-+Ω=⎢⎥∂⎢⎥⎣
⎦
⎰⎰ （5-24） 式中k v δ是重复指标,指叠加。

伽辽金积分公式是两个动量方程分别乘以各自的权函数，作内积后的和。

式中应力张量的导数
jk j
x σ∂∂包含着速度的二阶导数项。

对〔5—24)式中二所导数项进行分部积分，并代入应力边界条件(5—22)，则
p jk k
k nj jk k jk S j j v v d a v ds d x x σδδσδσΩΩ⎡⎤∂∂Ω=-Ω⎢⎥∂∂⎢⎥⎣⎦
⎰⎰⎰⎰⎰ =
p
k
nk k jk
S j
v p v ds d x δδσΩ
∂-Ω∂⎰
⎰⎰ 将上式代入(5—24)式
k
k k jk j dv v v d dt x δρδσΩ⎤∂⎡+Ω⎥⎢∂⎣⎥⎦
⎰⎰=p nk k k k S p v ds b v d δρδΩ+Ω⎰⎰⎰ 将(5—20)式代入上式，并考虑到全导数运算，上式可以写成
()()j k k k k j k jk j k j j v dv v v v v v p d dt x x x x δρδδμΩ⎧⎡⎤⎡⎤⎫∂∂∂∂⎪⎪++-++Ω⎢⎥⎢⎥⎨⎬∂∂∂∂⎢⎥⎢⎥⎪⎪⎣
⎦⎣⎦⎭⎩⎰⎰
222
=p
nk k k k S p v ds b v d δρδΩ
+Ω⎰⎰⎰ k=1,2 (5-25)
这就是不可压粘性流体的虚功率原理。

方程右边是由表面力和质量力所作的外部虚功率，左边是内应力的虚机械功率和惯性力的虚功率。

虚功原理是压分速度法有限元分析的基础。

设单元中压力和速度采用相同的插值展开式
T I p p =Φ , T I k k v v =Φ , T I
j j v v =Φ
结点未知数向量是时间t 的函数。

将近似解代入方程(6—23))注意到重复指标暗指累加，得到
1122I I
v C v C v F += （5-26）
式中1,12,2,T T
C d C d =
ΦΦΩ=ΦΦΩ⎰⎰⎰⎰
称为连续矩阵， v
v n S F v ds =Φ⎰,为边界流量向量。

将近似解代入(5—25)式．为了简便写成以下形式。

11111211I I I I
I Mv Av B v D v C p F +++-= （5-27） 22222221I I I I I Mv Av B v D v C p F +++-= （5-28）
式中
质量系数矩阵 T M d ρΩ
=ΦΦΩ⎰⎰
对流矩阵 1,12,2
()e
T I T T I T
d ρννΩ=ΦΦΦ+ΦΦΩ⎰⎰A 耗散矩阵 ,2,
2T
e
T d μΩ=ΦΦΩ⎰⎰1B 2,1,1T
e
T d μΩ=ΦΦΩ
⎰⎰B ,2,1T e
T d μΩ=ΦΦΩ
⎰⎰1D 2,1,2T e
T d μΩ=ΦΦΩ
⎰⎰D 压力矩阵 ,1T
e
T d Ω=ΦΦΩ⎰⎰1C
,2T
e
T d Ω=ΦΦΩ⎰⎰2C
外力向量 11e
p
n s b d p ds ρΩ=ΦΩ+Φ⎰⎰⎰1F
223 222e p
n s b d p ds ρΩ=ΦΩ+Φ⎰⎰⎰F 由于压力和速度采用了相同的插值形式，所以连续矩阵与压力矩阵相等。

对流矩阵中包含着未知的速度向量，方程(5—28)，(5—27)和(5—28)构成了非线性的常微分方程组。

速度和压力采取同一所次插值，速度可获得较精确的解,但压力会产生较大误差。

如果压力采用线性插值而速度采用二次插值，则二者都可获得好的结果。

将单元有限元方程总体合成，代入速度边界条件．得到总体有限元方程。

它仍然是非线性的常微分方程组。

一股可采用显式的差分格式求解，如欧拉法、中点法及Runge-Kutta 法等。

显式格式是条件稳定的,稳定性限制严格且难以用分析的方法得到稳定判据，往往要通过试算来选择适合的时间步长∆t.隐格式是无条件稳定的，虽可用较大的时间步长，但仍受到计算精度的限制，而且每一时步都要解非线性代数方程组．计算量很大，因此很少采用。

方程(5—27)和(5—28)中，如果去掉对时间的导数项，则方程表示了粘性不可压流体的定常运动。

这时需要用线性化迭代法解非线性代数方程组，如交替迭代法、Newton —RAphson 法和摄动法等，也可以用的间相关法求解非线性定常流动．。