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所示。

从剖分图上所表示的总体结点号与单元结点号的关系，可以建立联缀表于下 元素序号 1 2 3 4 5 6 7 8 9 10 总体

结点 号 n1

1 4 4 4

2 2 6 6 5 5 n2 4 5 9 8 6 5 7 10 10 9 n3

2

2

5

9

3

6

3

7

8

10

表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

ωω

ωωω

⎧∂∂+=-∈Ω

⎪∂∂⎪

⎪=------⎪⎨∂=-----⎪∂⎪⎪∂=------⎪∂⎩边界上

边界及上边界上