《有限元法在流体力学中及应用》由会员分享，可在线阅读，更多相关《有限元法在流体力学中及应用（21页珍藏版）》请在金锄头文库上搜索。
1、0第五章 有限元法在流体力学中的应用本章介绍有限元法在求解理想流体在粘性流体运动中的应用。讨论了绕圆柱体、翼型和轴对称物体的势流，分析了求解粘性流动的流函数涡度法流函数法和速度压力法，同时导出粘性不可压流体的虚功原理。1 不可压无粘流动真实流体是有粘性和可压缩的，理想不可压流体模型使数学问题简化，又能较好地反映许多流动现象。1. 圆柱绕流本节详细讨论有限无法的解题步骤。考虑两平板间的圆柱绕流如图 51 所示。为了减小计算工作量，根据流动的对称性可取左上方的 l4 流动区域作为计算区域。选用流函数方法，则流函数 应满足以下 Laplace 方程和边界条件1(5-1)220(,)(,)0(,)xy
2、xyaecbdyxycn流 线流 线流 线流 线将计算区域划分成 10 个三角形单元。单元序号、总体结点号和局部结点号都按规律编排如图 52 所示。从剖分图上所表示的总体结点号与单元结点号的关系，可以建立联缀表于下元素序号 1 2 3 4 5 6 7 8 9 10n1 1 4 4 4 2 2 6 6 5 5n2 4 5 9 8 6 5 7 10 10 9总体结点 号 n3 2 2 5 9 3 6 3 7 8 10表 5-1各结点的坐标值可在图 52 上读出。如果要输入计算机运算必须列表。本质边2界结点号与该点的流函数值列于下表边界结点号n1 2 3 4 8 9 10流函数 2 2 2 1 0
3、0 0表 5-2选用平面线性三角形元素，插值函数为(315)式。对二维 Laplace 方程进行元素分析，得到了单元系数矩阵计算公式(319)和输入向量计算公式(320)。现在对全部元素逐个计算系数矩阵。例如元素 1，其结点坐标为 =0, =2; =0, =1; =2.5, =2.1xy2xy3x3y由(3 15)式可得; 132.5ax213.5a,10; ;123by2310by; 30.5A从(3 19)式可计算出 1K1.45.2对 称依次可计算出全部子矩阵 20.45.1K3.2341.250.K5.06.1K70.5.8.1K90.5.10.5K根据联缀表把元素矩阵组合成总体系数矩
4、阵A=1.450.21.250.1.05.9.41.25.70.5.0.1.45.2901.5 对 称4矩阵中零元素没有一一写出，下三角部分与上三角部分对称。从(3 20)式计算元素输入向量，由于流函数满足齐次的自然边界条件,所以输入向量为零，总体输入向量也为零，这样就得了总体有限元0nq方程. NAB式中: 1210,TNB用缩减方程的重新编号修正方法施加边界条件，本质边界结点的函数值是已知的。把它们代入方程，修正右端项，再减去相应的方程，整理得 5674.9102.931解方程得到0.845， =1.241， =1.121567这样求出了全部结点上的流函数。为了求出每个单元形心处的速度，可
5、以由单元的流函数近似表达式求导计算。对元素 e 来说，有 TIy112320,TI TIyuaAA123,TI ITIxbB例如单元 =2, =1, =3，这样计算得到的速度为 u=1, 0。123 二维绕圆柱流动还可以用势函数求解，则定解问题可写成220(,)1xyxycdaebn边 界 上边 界 及 上边 界 上5表示势函数，为了使数值解唯一必须在部分边界上给定本质边界条件。势函数边界同样标记在图 5l 上。势因数满足 Laplace 方程和相应的边界条件，与流函数不同仅在于有非齐次的自然边界条件。采用与流函数方法完全一样的网格划分，可知计算得到的单元系数矩阵是完全一样的，总体矩阵也是完全
6、一样的。元素 1 和 4 具有非齐次自然边界条件应该用(320)式计算输入向量。元素 l, 。元素 4, 。总体合1,02TP41,02TP成得到 ，这样就得到方程组1TBNAB巳知 ，消去相应的三个方程得到一个 77 的3710代数方程组，解得 123.8,.4,456310.1790,.7,.单元形心处的速度可以用下列公式计算 TITIxuBTITIyA式中 是单元的结点势函数向量 。对于元素 1 来说， I123,这样计算得到 u=1.033，v=-0.05。这结果与流函1233.78,.41,.20数方法得到的结果近似相等。如果加密网格，就可以得到更好的结果。2. 升力问题考虑图 53
7、(a)所示的机翼绕流。均匀来流 平行于 x 轴，机翼边界为u6，后缘尖点为 ，流场外边界 取在离机翼足够远处。流函数 满足以下1T1方程和边界条件。(5-3)20120310ahuyb在 内在 上在 上在 上在 上其中 a,b 是特定系数， h 是上下边界之间的距离。机翼绕流的后驻点应位于后缘尖点处，在后缘 T 点满足 Kutta 条件； ； （5-4）0uy0vx7由于方程和边界条件是线性的,可用叠加原理求解，令（5-5）012ab其中 ， 和 ：分别是下列问题的解0120100yu在 内在 上在 上21110在 内在 上在 上21200在 内在 上在 上用有限元方法分别解以上三个问题，得到
8、各结点的 、 和 ，代入(55)式012得到叠加解。显然它满足问题(53)的全部方程和边界条件，特定常数 a,b 可利用 Kutta 条件(5 4)定出。首先由流函数 、 和 分别求出各个结点上的速度 ， 和 ，012 0,)uv(1,)(2,)uv(然后在后缘点 T 处利用 Kutta 条件，应有012uabuvv解之可得到 a 和 b。图 53(b)上给出了 NACA4412 具型以 攻角置于均匀流场中所引起的流动8图案，计算中采用了三角形单元。与无升力体绕流一样，机具绕流也可以采用速度势函数求解.3轴对称问题考虑圆管内绕轴对称物体的无旋流动，如图 54(a)所示。采用柱坐标系(r ， ,
9、z)，其势函数满足 Laplace 方程。8(5-6) 12222120SrrzSqn 在 内在 上在 上写出与微分问题相应的伽辽金积分表达 2221drz（ ）= 2Sqsn（ ）分部积分上式的左边并整理得到弱解积分形式 )rdzrz（= 02Lqdl式中 L 是元素的边长，L 绕轴旋转一周形成元素的边界面。采用图 54(b)所示轴对称的环形线性元素，它是将平面线性三角形元素绕对称轴旋转一周形成的环形体。采用斜坐标系，那么插值函数可写成 123,T9元素结点上势函数向量为 123(),IT则逼近函数为 123TI总体坐标和斜坐标系的关系为 TIrz式中 。 ，是元素结点总体坐标向量。123(
10、),ITrr123(),ITzz将逼近函数表达式代入伽辽金公式，推导出元素有限元方程 IKP式中影响系数矩阵和输入向量分别为K= 2)Trzdr（P= 0Lqdlr求出插值函数向量的偏导数 和 ，代入上式得影响系数矩阵rzK= （5-7）211213123 20()6ababrA 式中 ； i1，2，3 时 J=2，3，1;k31，2。ikjarijkbz， 三角形元素面积。0121Ab0A假设元素的“l 一 2”边落在自然边界上且 q 为常数，则可得转入向量计算公式（5-8）12230rqLP式中 ：是“l 一 2”边的边长。12L10计算了各元素的 K 和 P，然后总体合成，代入本质边界条
11、件就可以解总体方程。为了计算其它物理量，下面给出了相应的公式。元素形心处的速度：（5-9）12301230()zrUbAaa附加质量 m：m= (5-10)12()Nip式中，i=1,2,.N. N 是物体表面上所划的单元数。p 是输入向量 P 在元素结点上的值。在文献4中，以圆球作为例子计算了三种状态。球在无阻空间中运动，计算的附加质量系数 0.4671，理论值是 0.5。计算值小于理论值是符合第二章2 节例 4 中证明的附加质量极大值原理的.在圆管中球作匀速运动，计算的结果与 T JChung 在参考书(9) 中给出的结果比较，虽然我们采用了较少的结点，但达到了相同的精度。Chung 用流
12、函数方法，采用轴对称四边形单元计算。由于流函数满足 Stokes 方程，是非自伴的，这样行成的影响系数矩阵是非对称的不仅计算麻烦而且不能利用半带宽存储。四边形单元的短阵元素计算须用 Guass 数值积分,计算量大且有误差。而采用势函数方法和三角形元素恰好克服了以上两个缺点。第三种状态计算了圆球在半盲管(一端封死)中的运动。附加质量系数 0.897，这等于无限空间中附加质量系数的 1.6 倍。这使我们想到，在计算水下管中发射弹道问题时，应考虑物体在管道中的附加质量系数。轴对称不可压无粘流动也存在看流函数提法，流函数 应满足以 Stokes 方程，而不是 Laplace 方程。（5-11）2210
13、rz应特别注意的是,Stokes 算子是非自伴的。为了写出相应的迦辽金弱解积分表达式,先可将方程改写 22110rzr11其伽辽金积分表达为 22112)0ddrzr （ ） （将上式第一项分部积分，并代入自然边界条件，得到2()4()rzrzrzr = 2Sdln假设近似解可表示成 TI,I那么，TIrTIz，,Ir,Iz将以上各式代入弱积分表达式，得到单元有限元方程 IKP式中影响系数矩阵为K 2)4TTTrz rrddz（式中第二项是非对称项，它使得 K 成为非对称矩阵。输入向量为 P= 02Lrql如果采用前面已用到的三角形洄旋环状体元素，则是 K 和 P 可以导出，得到相应的计算公式
14、。看来流函数方法不及势函数方法简便。2 不可压粘性流动不可压粘性流体运动由速度散度为零的连续方程及 NavierStokes 方程描述。二维问题，引进流函数可导出流函数涡量方程和四阶的流函数方程。粘性流动中存在粘性应力，固壁边界上必须满足无滑条件，使得流体的运动一股是有旋的.121流函数问量法以流函数 和涡量 。表示的粘性不可压流动方程和自然边界条件是（5-12）()xyixyxysvSng在 上在 上（5-13）边界上还应给出本质边界条件，即 和 的函数值。根据具体边界不难绘出流函数的边界值。而固壁上的涡量，则要通过区域内的流函数值，由涡量边界条件来确定。流函数 、涡量 和速度的关系为,vuuvyxy分别写出流函数方程和涡量方程的伽辽金积分表达 2()d2)0iyxyd对以上两式中 Laplace 算子进行分部积分并代入自然边界条件：得弱解积分形式 ( )sSdvdxyz ()()iyxyxy
