流体力学有限元分析中的边界条件处理

流体力学有限元分析中的边界条件处理梁启国尹敏镐赵永凯高殿荣燕山大学大庆石化总厂摘要阐述了流体力学有限元分析中应用流函数-涡量法时典型边界条件的处理方法,并给出计算实例.关键词有限元涡-流函数法边界条件分类号O357.1引言求解不可压缩粘性流体二维流动问题的数值方法有速度-压力法、涡-流函数法和流函数法,其中涡-流函数法应用较广1～4.在应用涡-流函数法进行有限元分析时,数值边界条件的处理不但影响解的精度而且影响解的稳定性2.本文结合工程实际给出了在有限元分析中确定几种典型边界条件的方法.0边界条件的处理如图1所示的沿背部台阶的不可压缩粘性流体二维流动.C点为角点.在涡-流函数方程的求解过程中,不仅需要知道边界上的流函数和涡量值,有时还要对边界上的涡量值不断地进行修正.下面分别讨论各种边界上的流函数值和涡量值.设Δx和Δy分别为x方向和y方向的有限元划分网格间距,i和j分别为x方向和y方向的步长指标.1)壁面边界处理若壁面是不可渗透的,则沿壁面的流函数Ψ为常数.一阶精度的壁涡公式为:1图1沿背部台阶的边界条件2(Ψi,j-Ψi,j+1)Ωi,jcc+O(Δy)(1)=(Δy)2c二阶精度的壁涡公式为:3(Ψi,j-Ψi,j+1)Ωi,j+1+O((Δy)2)Ωi,jccc(2)=(Δy)22c图2用内点上的Ψi,j+1,Ωi,j+1表示边界上的涡量值.二阶壁涡公cc式在大网格雷诺数和变网格情况下往往引起数值不稳定和计算不收敛;一阶壁涡公式虽然精度低,但计算稳定、易于收敛.另外一种二阶壁涡公式5为:8Ψi,j+1-7Ψi,j-Ψi,j+2+O((Δy)2)Ωi,jccc(3)=2(Δy)2c用时间迭代法确定壁涡的方法6,7为:根据不可压缩粘性流体二维流动的流函数-涡量公式8～1092Ψ92Ψ=-Ω+9x29y292Ω92Ω9Ω9Ω9Ω1++u+v=9x29y29t9x9yRe对上述方程进行有限元离散,得总体矩阵方程KΨ=MΩ-qMΩQ+(K+A)Ω=s(4)(5)式中,ΩQ=9Ω/9t,K,M,A,q和s为系数矩阵.在求解式(4)和(5)时,为确定壁面涡量,首先由式(4)和t=t0时刻的涡量初值来计算t=t0+Δt时刻的流函数,有Ψt0+Δt=K-1(MΩt0-q)(6)然后再用式(4)和新的流函数值Ψt0+Δt求得壁涡Ωt0+Δt=M-1(KΨt0+Δt+q)(7)最后,从式(5)出发,并利用新的壁涡值可求t0+Δt时刻其内部各点的涡量值.以后每增加一个时间步长都按上述方法确定壁涡,所以也称壁面涡的时间迭代法.边界B6的处理和B1相同.2)对称边界处理由于流动的对称性,图1中的边界B5为对称边界;而对称线必为流线,因此在对称边界上的流函数ΨB为常数.5对称边界是滑移边界,显然99uvvB=0=0=09x9y5BB55由此可求对称边界上的涡量99uvΩB==0-9x9y5B53)上边界处理上边界B3一般为以下两种情形:其一是非滑移固体壁面,比如管道流动的上边界就属于这种情形.可采用对B1的处理方法来处理.其二是滑移面,如飞机在飞行中的无界流体的绕流问题:此时y方向的流场应是无限的,但数值计算时却要在有限区域内进行,于是在离开绕流物体足够远处确定一个上边界B3,其上的边界条件应与无穷远条件相匹配.滑移面B3可看作滑移固体壁面,通常可采用下述两种处理方法:第3期赵永凯等流体力学有限元分析中的边界条件处理191(i)在B3上规定为无穷远平行来流,u=V∞,v=0,流函数满足Ψi,J值由流量确定.为求B3上的涡量值,可推出类似的壁涡公式:一阶精度公式为=const.常数2(Ψi,J-1-Ψi,J+V∞Δy)Ωi,J+O(Δy)=(Δy)2二阶精度公式为3(Ψi,J-1-Ψi,J+V∞Δy)Ωi,j-1+O((Δy)2)Ωi,J=-(Δy)22上述边界条件,对于2Ψ=Ω是第一类Diricheit边界条件,易于求解.(ii)B3在上,与无穷远处平行流动完全一致,边界条件为Ωi,J=V∞.因为=0,ui,J99u9vv=0.可见该条件比Vi,J=0条件要宽些.但对2Ψ=Ωi,J所以=-9x9y9xi,Ji,JΩ是Neuman边界条件,可采用一阶精度的离散形式Ψi,J=Ψi,J+1+V∞Δy如果流量已知,也可以用Ψi,J=const.4)进流边界处理进流边界B2也称上游边界.其边界条件应按进流物理条件确定.有以下几种情况:(i)规定u0,j=V∞,v0,j=099u9vvu0,jV∞,Ω0,j=0.因为Ωi,j=(ii)规定,相当于规定=-=9x9y9x0,j0,j9v=0,这比规定v0=0要放松些.,j9x0,j92Ψ(iii)规定Ψ0,j,99v9uv=0.据此Ω0=也就给定了.显-=-,j9y29x9x9y0,j0,j0,j然,此时已不是均匀来流.通常Ψ值由边界层计算而得.92Ψ99u9u9vv(iv)规定Ψ0,j.而Ω0,j其中,=-=,=9y29x9y9y9x0,j0,j0,j0,j92Ψ92Ψ99v9v9v.采用假设,=0,可知==-.9x29x29x9x9x9x0,j0,j1,j0,j1,j最终得到92Ψ92ΨΩ0,j=-9x29y21,j0,jΨ0,j-2Ψ1,j+Ψ2,jΨ0,j+1-2Ψ0,j+Ψ0,j-1=--(Δx)2(Δy)25)出流边界条件在有限元分析中,出流边界条件的扰动会向上游传播,边界条件的处理将影响计算的稳定性和解的精度.出流边界条件的处理与进流边界条件的处理相似,也可按出流的物理条件确定.可有以下几种情况:Ω(i)规定vI,j=0,9=0.或ΨI,j=ΩI-1,j,ΩI,j=ΩI-1,j.9xI,j(ii)规定ΨI,j=ΨI-1,j,ΩI,j=0.9Ω=0,9v(iii)规定=0.9x9xI-1,jI-1,j22或ΩI,j=ΩI-1,j,ΨI,j=2ΨI-1,j-2ΨI-2,j.6)尖角点处理(i)凹尖角边界规定Ψi,0=0,Ωi,0=0.cc(ii)凸尖角边界规定流函数Ψc=const.对涡量Ωc的处理有以下几种方法:双值法如图4所示,由于C点即在壁面B1上又在壁面B6上,所以对C点可用壁涡公式求得两个壁面涡量:图3=2(Ψi,j-Ψi,j+1)/(Δy)2ΩC上ΩC下cccc=2(Ψi,j-Ψi+1,j)/(Δx)2cccc因为C是奇点,所以C点的壁面涡量可以不为单值.当ΩC分别用于上、下面结点(ic,jc+1)与(ic+1,jc)有限元离散时,分别取ΩC上与ΩC下.平均法图4Ωc=(ΩC上+ΩC下)/2取分离点法由于C点上游的流动,在C点没有分离;而C点下游的流动,在C点产生分离,所以在C点有两个涡量值.即=2(Ψi,j-Ψi,j+1)/(Δy)2ΩC上ΩC下=0cccc虚拟对称法将C点作为内点来处理,而把Ψ看成对于C点是对称的,其对称虚拟值Ψ3为Ψ33ic-1,jc=Ψic+1,jcΨi,j+1=Ψi,j-1cccc关于C点的Possion方程Ψi+1,j2Ψi,j+Ψ3-1,j3,j-1Ψi,j+1-2Ψi,j+Ψicccc-icccccccc=Ωi,jcc+(Δx)2(Δy)2第3期赵永凯等流体力学有限元分析中的边界条件处理193=0和虚拟值Ψ3的对称性,可得由于Ψi,jcc=2Ψi,j+1/(Δx)2+2Ψi+1,j/(Δy)2Ψi,jcccccc45°壁角法将尖角点C看成是一个法线方向为45°的小平面.可用法向壁涡公式Ωc=2Ψ(p)/(Δn)2式中,Ψ(p)由邻近内点Ψ值内插确定.计算实例在一平行板管道内有处于静止状态的流体,如图5所示.在t的给定速度为u=1,求充分发展后的稳态解8.解用有限元方法求解流函数-涡量方程时,将涉及到以下几种边界情况:(1)进流边界ab.此时进口为一单位均匀速度,按前述方法,流函数Ψab=y,涡量Ωab=0.(2)无滑移上边界bc.此时Ψbc=con2st,该常数值可以具体算出.涡量值的确定,可以用一阶、二阶壁涡公式,也可以用2=0时刻,入口处流体图5平行板管道流动时间迭代法.(3)对称边界ad.取Ψad=0,Ωad=0.9Ψ=0,9Ω(4)出流边界cd.=0.9x9xcdcd由于流动的对称性,取上半流场为计算对象.将流动区域进行有限元网格剖分,选用四结点等参四边形单元.计算的流动雷诺数Re=200.部分计算结果见表1～表3.表1采用一阶壁涡公式时的流函数涡量值流函数y涡量yx10.90.5010.90.5000.3251.0002.0003.3004.8008.00011.00021.1001111111110.90000.94530.97630.97940.98970.98240.98400.98430.98430.50000.52650.57670.61460.63620.65640.67410.67780.6778000000000010.93604.74244.12443.86853.52583.21133.14173.141002.35804.29333.72853.37133.08902.85352.80852.80900-0.4416-0.37420.24350.62030.97821.27801.33961.3392000000000表2采用二阶壁涡公式时的流函数涡量值流函数y涡量yx10.90.5010.90.5000.3251.0002.0003.3004.8008.00011.20021.1001111111110.90000.95490.97940.98000.98160.98310.98450.98480.98480.50000.53290.58610.61980.64000.65930.67570.67890.6887000000000012.00403.98624.17033.84793.51943.21803.15803.158403.04924.38643.66153.36323.08802.86272.82432.82560-0.5085-0.27940.32680.67241.01601.29611.34711.3457000000000表3固壁涡量采用时间迭代法时的流函数涡量值流函数y涡量yx10.90.5010.90.5000.3251.0002.0003.3004.8008.00011.20021.1001111111110.90000.95630.98090.98150.98320.98480.98630.98660.98660.50000.53400.58850.62270.64320.66270.67960.68300.6830000000000012.17804.00144.24603.91413.59003.28303.21143.209803.14284.45073.71693.41943.14882.91872.87122.87120-0.5195-0.26760.34400.69601.03851.32821.39111.3912000000000计算表明,一阶壁涡公式的计算精度低于二阶壁涡公式,二阶壁涡公式的计算精度低于时间迭代法.其原因是二阶壁涡公式有截断误差.结论给出的典型边界的处理方法和计算公式,对实际问题的分析和计算具有重要的参考价值和指导意义.对于复杂边界条件下的处理方法和计算公式,尚待进一步研究.3参考文献123456赵学端等.粘性流体力学.机械工业出版社,1993,376～390.析孝康等.计算流体力学.国防科技大学出版社,1989,251～269.章本照.流体力学中的有限元方法.机械工业出版社,1986,176～187.吴江航等.计算流体力学的理论、方法及应用.科学出版社,1988,172～199.RoachePJ.Computationalfluiddynamics,Hermasapublishers.1976,263～287.ManzanM,CominiG.Inflowandoutflowboundaryconditionsinthefiniteelements