有限元法在流体力学中的应用.doc

《有限元法在流体力学中的应用.doc》由会员分享，可在线阅读，更多相关《有限元法在流体力学中的应用.doc（21页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第五章 有限元法在流体力学中的应用 本章介绍有限元法在求解理想流体在粘性流体运动中的应用。讨论了绕圆柱体、翼型和轴对称物体的势流，分析了求解粘性流动的流函数涡度法流函数法和速度压力法，同时导出粘性不可压流体的虚功原理。1 不可压无粘流动真实流体是有粘性和可压缩的，理想不可压流体模型使数学问题简化，又能较好地反映许多流动现象。1. 圆柱绕流本节详细讨论有限无法的解题步骤。考虑两平板间的圆柱绕流如图51所示。为了减小计算工作量，根据流动的对称性可取左上方的l4流动区域作为计算区域。 选用流函数方法，则流函数应满足以下Laplace方程和边界条件 (5-1)将计算区域划分成10个三角形单元。单元序号

2、、总体结点号和局部结点号都按规律编排如图52所示。从剖分图上所表示的总体结点号与单元结点号的关系，可以建立联缀表于下元素序号12345678910总体结点 号n11444226655n2459865710109n322593637810 表5-1各结点的坐标值可在图52上读出。如果要输入计算机运算必须列表。本质边界结点号与该点的流函数值列于下表边界结点号n12348910流函数2221000 表5-2选用平面线性三角形元素，插值函数为(315)式。对二维Laplace方程进行元素分析，得到了单元系数矩阵计算公式(319)和输入向量计算公式(320)。现在对全部元素逐个计算系数矩阵。例如元素1，

3、其结点坐标为=0, =2; =0, =1; =2.5, =2.由(315)式可得; ,; ; 从(319)式可计算出依次可计算出全部子矩阵根据联缀表把元素矩阵组合成总体系数矩阵A=矩阵中零元素没有一一写出，下三角部分与上三角部分对称。从(320)式计算元素输入向量，由于流函数满足齐次的自然边界条件,所以输入向量为零，总体输入向量也为零，这样就得了总体有限元方程.式中: 用缩减方程的重新编号修正方法施加边界条件，本质边界结点的函数值是已知的。把它们代入方程，修正右端项，再减去相应的方程，整理得 解方程得到 0.845，=1.241，=1.121 这样求出了全部结点上的流函数。为了求出每个单元形心

4、处的速度，可以由单元的流函数近似表达式求导计算。对元素e来说，有 例如单元=2, =1, =3，这样计算得到的速度为u=1,0。二维绕圆柱流动还可以用势函数求解，则定解问题可写成表示势函数，为了使数值解唯一必须在部分边界上给定本质边界条件。势函数边界同样标记在图5l上。势因数满足Laplace方程和相应的边界条件，与流函数不同仅在于有非齐次的自然边界条件。采用与流函数方法完全一样的网格划分，可知计算得到的单元系数矩阵是完全一样的，总体矩阵也是完全一样的。元素1和4具有非齐次自然边界条件应该用(320)式计算输入向量。元素l, 。元素4, 。总体合成得到，这样就得到方程组巳知，消去相应的三个方程

5、得到一个77的代数方程组，解得单元形心处的速度可以用下列公式计算式中是单元的结点势函数向量。对于元素1来说， ,这样计算得到u=1.033，v=-0.05。这结果与流函数方法得到的结果近似相等。如果加密网格，就可以得到更好的结果。2. 升力问题 考虑图53(a)所示的机翼绕流。均匀来流平行于x轴，机翼边界为，后缘尖点为，流场外边界取在离机翼足够远处。流函数 满足以下方程和边界条件。 (5-3)其中a,b是特定系数，h是上下边界之间的距离。机翼绕流的后驻点应位于后缘尖点处，在后缘T点满足Kutta条件 ； （5-4）由于方程和边界条件是线性的,可用叠加原理求解，令 （5-5） 其中，和：分别是下

6、列问题的解用有限元方法分别解以上三个问题，得到各结点的、和，代入(55)式得到叠加解。显然它满足问题(53)的全部方程和边界条件，特定常数a,b可利用Kutta条件(54)定出。首先由流函数、和分别求出各个结点上的速度，和，然后在后缘点T处利用Kutta条件，应有解之可得到a和b。 图53(b)上给出了NACA4412具型以攻角置于均匀流场中所引起的流动图案，计算中采用了三角形单元。与无升力体绕流一样，机具绕流也可以采用速度势函数求解.3轴对称问题 考虑圆管内绕轴对称物体的无旋流动，如图54(a)所示。采用柱坐标系(r，,z)，其势函数满足Laplace方程。 (5-6) 写出与微分问题相应的

7、伽辽金积分表达=分部积分上式的左边并整理得到弱解积分形式=式中L是元素的边长，L绕轴旋转一周形成元素的边界面。 采用图54(b)所示轴对称的环形线性元素，它是将平面线性三角形元素绕对称轴旋转一周形成的环形体。采用斜坐标系，那么插值函数可写成元素结点上势函数向量为则逼近函数为总体坐标和斜坐标系的关系为式中。，是元素结点总体坐标向量。 将逼近函数表达式代入伽辽金公式，推导出元素有限元方程式中影响系数矩阵和输入向量分别为K=P= 求出插值函数向量的偏导数和，代入上式得影响系数矩阵K= （5-7）式中 ； i1，2，3时J=2，3，1;k31，2。 ， 三角形元素面积。假设元素的“l一2”边落在自然边

8、界上且q为常数，则可得转入向量计算公式 （5-8）式中 ：是“l一2”边的边长。计算了各元素的K和P，然后总体合成，代入本质边界条件就可以解总体方程。为了计算其它物理量，下面给出了相应的公式。元素形心处的速度： （5-9）附加质量m：m= (5-10)式中，i=1,2,.N. N是物体表面上所划的单元数。p是输入向量P在元素结点上的值。在文献4中，以圆球作为例子计算了三种状态。球在无阻空间中运动，计算的附加质量系数0.4671，理论值是0.5。计算值小于理论值是符合第二章2节例4中证明的附加质量极大值原理的.在圆管中球作匀速运动，计算的结果与TJChung在参考书(9)中给出的结果比较，虽然我

9、们采用了较少的结点，但达到了相同的精度。Chung用流函数方法，采用轴对称四边形单元计算。由于流函数满足Stokes方程，是非自伴的，这样行成的影响系数矩阵是非对称的不仅计算麻烦而且不能利用半带宽存储。四边形单元的短阵元素计算须用Guass数值积分,计算量大且有误差。而采用势函数方法和三角形元素恰好克服了以上两个缺点。第三种状态计算了圆球在半盲管(一端封死)中的运动。附加质量系数0.897，这等于无限空间中附加质量系数的1.6倍。这使我们想到，在计算水下管中发射弹道问题时，应考虑物体在管道中的附加质量系数。轴对称不可压无粘流动也存在看流函数提法，流函数应满足以Stokes方程，而不是Lapla

10、ce方程。 （5-11） 应特别注意的是,Stokes算子是非自伴的。为了写出相应的迦辽金弱解积分表达式,先可将方程改写其伽辽金积分表达为将上式第一项分部积分，并代入自然边界条件，得到=假设近似解可表示成那么 ， ，将以上各式代入弱积分表达式，得到单元有限元方程式中影响系数矩阵为 K式中第二项是非对称项，它使得K成为非对称矩阵。输入向量为 P= 如果采用前面已用到的三角形洄旋环状体元素，则是K和P可以导出，得到相应的计算公式。看来流函数方法不及势函数方法简便。2 不可压粘性流动 不可压粘性流体运动由速度散度为零的连续方程及NavierStokes方程描述。二维问题，引进流函数可导出流函数涡量方

11、程和四阶的流函数方程。粘性流动中存在粘性应力，固壁边界上必须满足无滑条件，使得流体的运动一股是有旋的. 1流函数问量法以流函数和涡量。表示的粘性不可压流动方程和自然边界条件是 （5-12） （5-13）边界上还应给出本质边界条件，即和的函数值。根据具体边界不难绘出流函数的边界值。而固壁上的涡量，则要通过区域内的流函数值，由涡量边界条件来确定。流函数、涡量和速度的关系为分别写出流函数方程和涡量方程的伽辽金积分表达 对以上两式中Laplace算子进行分部积分并代入自然边界条件：得弱解积分形式=将求解区域剖分，单元中流函数、涡量。选用相同的插值展开形式，则有 ，其中 。是结点未知数向量。和分别是单元

12、中i结点在t时刻流函数和涡量值。这样就有，将流函数和涡量的插值表达式及其变分代入弱解积分形式，得=0整理得到单元的有限元方程 流函数的单元有限元方程是代数方程，而涡量的有限元方程是常微分方程组，必须在时间步进中求解。式中各矩阵分别为质量矩阵 损耗矩阵 对流矩阵 输入向量为 是t的函数，对流矩阵A是与时间t有关的非对称矩阵，因此每一时步必须重新计算。将所有单元进行总体合成，得到总体有限元方程 （5-14） （5-15）式中和是总体结点未知数向量。 总体有限元方程是代数方程和非线性常微分方程的联立方程组。将涡量力程用显式差分貉式离散，可用交替迭代的方法求解耦合方程组。 用流函数涡量法，可以计算低雷

13、诺数时绕障碍，如圆柱、圆球的流动，还可以模拟涡街的形成过程。 2流函数方法 二维不可压粘性流动可用四阶的流函数方程描述 （5-16）写出相应的伽辽金积分表达利用分部积分=-可得弱解积分表达式 = 假定在全部边界上给定本质边界条件，对于四阶算子，即给定，那么在S上有及和等于零，则伽辽金弱积分形式可简化为设单元的近似函数为 式中将近似函数代入弱解表达式，可以推导出单元有限元方程。在推导中首先注意到全导数的运算然后把近似解代入积分表达式中整理得到 （5-17）式中矩阵包含着未知向量 ，这是对流项的非线性影响产生的。如果写成下标表示的方程组，则可以看出方程的第二项包含的二次项，所以得到的是非线性的常微

14、分方程组.由于流函数方程是四阶偏微分方程，因此函数及其一阶导数、，是本质变量，应取为结点未知数并保持其在元素内和元素之间的连续性。为此，应选用函数及一阶导数都连续的Hermite插值函数要构造相容的Hermite插值函数有时是困难的。另外单元近似函数要用高阶插值，元案自由度大，这使得总体有限元方程的系数矩阵阶次很高要占用大量的计算机内存.这种方法只具有一个未知函数，所以计算程序简单，也比较省时。用流函数方法计算圆柱绕流，得到的表面压力分布曲线与有限差分解比较，符合得很好。3. 流体力学虚功率原理及速度压力法不可压粘性流动基本方程如下 连续方程 (5-18)动量方程 k=1,2 (5-19)本构

15、方程 (5-20)方程中重复指标暗指叠加运算。本构方程给出了应力和应变率之间的关系。其中应力张量是二阶张量，表示作用在流体微团小六面体的j面上k方向的应力分量。由于流体微团动量矩平衡的原因它是对称张量，即有是粘性应力张量。 称为Kronecker 是质量力分量。流动区域边界可分成二类，一类是指定速度的，另一类是指定应力。在指定速度的边界上 k=1,2如图55所示如果边界的单位外法向矢量为n，其方向余弦j=1,2. 那么其法向速度应满足 (5-21)如固壁边界，流体质点粘附在固壁上，流体速度等于固壁运动速度，满足无滑条件。若固壁静止，则其值为零。在指定应力条件的边界，上，有 是作用于边界上，法向

16、为n的单位微元面积上的应力矢量在k方向的分量见图55。 设液体微元j面上的应力矢量为 式中表示直角坐标系中的两个单位矢量。那么等于二个应力矢量在n方向的投影之和，即 那么应力边界条件可写成 （5-22）如自由液面上，是法向为n的自由面上结定的作用于液体的风应力，多大气静止常压时，其值为零.流体入口或出口边界，可给出速度或应力值。不同液体分界面上，一船给出速度，压力连续性条件。把本构方程(520)式代入动量方程(519)式，并利用连续方程，就可以导出NavierStokes方程。下面推导虚功率原理。 以压力的变分为权函数写出连续方程的伽辽金积分公式分部积分得到弱解积分公式 （5-23）以速度的变

17、分为权函数，写出动量方程的伽辽金积分公式 （5-24）式中是重复指标,指叠加。伽辽金积分公式是两个动量方程分别乘以各自的权函数，作内积后的和。式中应力张量的导数包含着速度的二阶导数项。对524)式中二所导数项进行分部积分，并代入应力边界条件(522)，则=将上式代入(524)式=将(520)式代入上式，并考虑到全导数运算，上式可以写成= k=1,2 (5-25)这就是不可压粘性流体的虚功率原理。方程右边是由表面力和质量力所作的外部虚功率，左边是内应力的虚机械功率和惯性力的虚功率。虚功原理是压分速度法有限元分析的基础。设单元中压力和速度采用相同的插值展开式 , , 结点未知数向量是时间t的函数。

18、将近似解代入方程(623)注意到重复指标暗指累加，得到 （5-26）式中 称为连续矩阵，,为边界流量向量。 将近似解代入(525)式为了简便写成以下形式。 （5-27） （5-28）式中质量系数矩阵 对流矩阵 耗散矩阵 压力矩阵 外力向量 由于压力和速度采用了相同的插值形式，所以连续矩阵与压力矩阵相等。对流矩阵中包含着未知的速度向量，方程(528)，(527)和(528)构成了非线性的常微分方程组。速度和压力采取同一所次插值，速度可获得较精确的解,但压力会产生较大误差。如果压力采用线性插值而速度采用二次插值，则二者都可获得好的结果。将单元有限元方程总体合成，代入速度边界条件得到总体有限元方程。它仍然是非线性的常微分方程组。一股可采用显式的差分格式求解，如欧拉法、中点法及Runge-Kutta法等。显式格式是条件稳定的,稳定性限制严格且难以用分析的方法得到稳定判据，往往要通过试算来选择适合的时间步长t.隐格式是无条件稳定的，虽可用较大的时间步长，但仍受到计算精度的限制，而且每一时步都要解非线性代数方程组计算量很大，因此很少采用。方程(527)和(528)中，如果去掉对时间的导数项，则方程表示了粘性不可压流体的定常运动。这时需要用线性化迭代法解非线性代数方程组，如交替迭代法、NewtonRAphson法和摄动法等，也可以用的间相关法求解非线性定常流动