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

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

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

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

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

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

选用流函数方法，则流函数 应满足以下Laplace方程和边界条件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 总体结点 号 n11 4 4 42 2 6 6 5 5 n2 4 5 9 8 6 5 7 10 10 9 n322593637810表5-1各结点的坐标值可在图5—2上读出。

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

本质边界结点号与该点的流函数值列于下表表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)式可计算出1K1 1.45 1.250.21.2500.2K ⎛⎫⎪⎪=⎪ ⎪⎝⎭--对称依次可计算出全部子矩阵20.20.201.45 1.251.25K ⎛⎫⎪⎪= ⎪⎪⎝⎭--30.200.21.25 1.251.45K⎛⎫⎪⎪=⎪⎪⎝⎭--41.25 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.252.4500 1.25 1.01000.50.52.90.400 1.254.9100 1.750.54.01000.52000.51.450.201.9501.5--⎡⎤⎢⎥--⎢⎥⎢⎥--⎢⎥--⎢⎥⎢⎥---⎢⎥--⎢⎥⎢⎥-⎢⎥-⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦0对称矩阵中零元素没有一一写出，下三角部分与上三角部分对称。

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

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

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

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

二维绕圆柱流动还可以用势函数求解，则定解问题可写成 ω表示势函数，为了使数值解唯一必须在部分边界上给定本质边界条件。

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

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

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

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

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

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

总体合成得到110010000022TB ⎡⎤=⎢⎥⎣⎦，这样就得到方程组22220(,)001x y xy cd aec bd n ab nωωωωω⎧∂∂+=-∈Ω⎪∂∂⎪⎪=------⎪⎨∂=-----⎪∂⎪⎪∂=------⎪∂⎩边界上边界及上边界上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 xu B xωωω∂==Φ=∂ T I T Iy A yωνωω∂==Φ=∂ 式中I ω是单元的结点势函数向量[]123,,ωωω。

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

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

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

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

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

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

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

机翼绕流的后驻点应位于后缘尖点处，在后缘T 点满足Kutta 条件0u y ψ∂==∂；0v xψ∂=-=∂； （5-4） 由于方程和边界条件是线性的,可用叠加原理求解，令012a b ψψψψ=++ （5-5）其中0ψ，1ψ和2ψ：分别是下列问题的解20010000yu ψψψ∞⎧∇=-----Ω⎪=-------Γ⎨⎪=------Γ⎩在内在上在上211110001ψψψ⎧∇=-----Ω⎪=-------Γ⎨⎪=--------Γ⎩在内在上在上222120010ψψψ⎧∇=-----Ω⎪=-------Γ⎨⎪=--------Γ⎩在内在上在上用有限元方法分别解以上三个问题，得到各结点的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 条件，应有012012u u au bu v v av bv =++⎧⎨=++⎩解之可得到a 和b 。

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

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

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

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

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

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

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

将逼近函数表达式代入伽辽金公式，推导出元素有限元方程I K P Φ=式中影响系数矩阵和输入向量分别为K=2)T T r z rdrdz πΦΦ+ΦΦ⎰⎰r z （P=02Lq dl πΦ⎰r求出插值函数向量的偏导数r Φ和z Φ，代入上式得影响系数矩阵K=22111212131********232302233()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”边的边长。