蜗壳内湍流流场的单向流线上风有限元法数值计算.docx

蜗壳内湍流流场的单向流线上风有限元法数值计算王仁人周建华(山东轻工业学院机电工程系济南250100) (上海交通大学上海200030)摘要 本文采用单向流线上风 Galer kin 有限元法 ,并采用压力速度同次插值函数 ,求解了形状比较复杂的涡轮蜗壳内二维可压缩湍流流动 。关键词蜗壳 ,湍流流场 ,单向流线上风 ,有限元中图法分类号 T K41118 ,O35715涡轮蜗壳出口的速度系数提高 0101 , 相应地整个涡轮的效率亦提高 1 %左右 。因此研究蜗壳内部流动性能对提高涡轮效率起着重要的作用 。为全面了解蜗壳内的流动 , 以便设计出 性能良好的蜗壳 , 对蜗壳内流动进行详细的测量和分析计算是很有必要的 。蜗壳内的流动是很复杂的 , 其计算方法主要有有限差分法和有限元法 。有限差分法应用最广泛 , 但其对复杂边界条件的适应性差 , 虽然可采用曲线贴体坐标的方法来解决 , 但仍不方 便 。与之相比有限元法主要优点之一是 :对于求解区域的单元剖分没有特别的限制 , 这对处理 具有复杂边界区域的工程实际问题 , 格外方便 , 而且完全可以根据问题的物理特点 , 在求解区域中安排单元网格的疏密 。作为“车用发动机变截面涡轮增压技术的系列研究 3 ”课题的一部分 , 本文采用单向流线 上风 Galer kin 有限元法 , 并采用压力速度同次插值函数 , 求解形状比较复杂的涡轮蜗壳内二维 可压缩湍流流动 , 为变截面涡轮蜗壳的优化设计提供技术支持 。1 二维粘性可压缩湍流流动的控制方程定义参考长度 L re , 参考速度 U re后 , 各变量的无量纲形式为 :xyuvX =,Y =,U =,V =L re p - p0L reU reU reU reL reU reL reP =Ree =Ret =,U 2ef ftre对二维粘性可压缩湍流流动 , 无量纲化后的各方程可统一用下列方程表示 :U 9 + V9 = 9 99 9+ S+( 1)9X 9 Y 9X 9X 9 Y 9 Y3 此课题得到山东省优秀中青年科学家奖励基金的资助 收稿日期 :1998 - 05 - 01 30 山东轻工业学院学报第 13 卷表 1 列出了当 分别为速度 、湍流动能 K 和耗散率时 , 在直角坐标系下与之相对应的 、S 。这里 , 对流扩散方程采用了非守恒形式 , 因蜗壳内流动为亚音速流动 , 在密度变化不 很大的情况下 , 不会产生太大的误差 。表 1无量纲化后的各方程S方程连续性方程1001 9U( 2) 1 9V9P +9 9X 方向动量方程1 /R ee-+( 3)UR ee 9X1 9VR ee 9X9X 9X9 Y 1 9U9P +9 91 /R eeY 方向动量方程-+( 4)V9 Y 9 YR ee 9 Y9XR ee 9 Y+ t /kL re湍流动量方程KG - U( 5)U Lre rere+ t / 2 L re湍流动能耗散率方程C1 G - C2( 6)U LKKU rere re2229U9V9U + 9V U re其中 :G = 2+L re t9X9 Y9X 9 Y湍流粘性系数 t 由湍流动能 K 和其耗散率得出 :t = CK2 /此 K - 模型中的经验系数为 : C1 = 1144 , C2 = 1192 , C = 0109 ,k = 110 , = 1130( 7)在靠近固体壁面的粘性支层中 , 上述 K - 高 Re 数模型不能适用 。本文采用壁面函数法把高 Re 数 K - 模型加以修正 , 使它可以一直用到壁面附近的粘性支层内 。在有限元法中 ,具体做法是将壁面到旺盛湍流区域之间取一层边界单元 , 对此单元 , 壁面上的节点 Kp = 0 ,p= 0 , 而非壁面上的节点 , Kp 从 K 方程中求得 ,p 由代数法由得 。即 ,+t = Y pl /u p( 8)1 /4 1 /2其中 : Y + = Y/, u + = B + (ln Y + ) /k , k = 01400142 , B = 510515C K pp ppp2控制方程的有限元离散由 ( 1) 式可得 :1 91 9U 9 + V 9999XW d A +SW d AW d A =( 9)+Ree 9XRee 9 Y9X9 Y9 Y( e)( e)( e)AAA这里 , W 为权函数 , A ( e) 为单元面积 。方程中的对流项 、扩散项及源项具有不同的物理特性 , 因此应采用不同的离散方法 。离散后得到的代数方程的一般形式为 : A i j ij = D i j ij + S i211扩散项的离散( 10)方程 ( 9) 右边第一项为扩散项 , 可直接采用 Galer kin 加权余量法进行离散 : 1 9W 9 9W9R9X 9Xd A + W 9n d S( 1)D =+9 Yee( e)( e)AS假设传输量在边界上的外法向梯度为零 , 积分域用图 ( 1) 所示线性四边形等参数单元离99 Y31第 2 期王仁人等 : 蜗壳内湍流流场的单向流线上风有限元法数值计算散 ,此扩散项离散后得到系数为 9 W i 9N j 9W i 1 R 9X 9X9 Y 9 YD ij=d A+ee( e)A( 12)式中 : N j 为单元内插值函数 , W i 为权函数 。212对流项的离散采用流线上风法离散对流项 。在无扩散项 或源项时 , 沿流线 值为常数 。在流线坐标系 中 , 在一个单元内可假定 :U s 9= co nstant( 13)9S设 p 、i 点分 别 为 上 、下 游 点 , 用 差 分 表 示上式偏导数项 , 对流项积分可表示为 :图 1线性四边形单元U sA =(i - p ) W d A( 14)s( e)A2 2( X i - X p ) 2 + ( Y i - Y p ) 2式中 : s =U s =U + Vi i上游结点对单元的传输量可由单元内线性插值得到 :N pjjp( 15)=其中 N pj 为单元插值函数 。将 ( 15) 式代入 ( 14) 式即得到对流项离散后的系数 :U s= - s N pj W i d A( j i )( 16 a)A ij( e)AU s= - s W i d A( 16 b)A ii( e)A某单元的下游结点定义为在该结点处其速度向量的负方向指向该单元内 。不失一般性考虑结点 1 , 在此结点处的速度向量为 V? 1 , 如果满足以下条件 , 则此结点为下游结点 :( | V? 1 | 0) ( V? 1 ?n 1 0) ( V? 1 ?n 4 0)( 17)单元被上游和下游结点间的流线分为两部分 , 由于流线上无传输量穿过 , 流线两边单元内的通流量应相等以满足连续性条件 。根据此特点即可确定上游结点 p ( X p , Y p ) 的位置 。上游 结点 p ( X p , Y p ) 的坐标 X p 、Y p , 可根据图 ( 1) 所示的插值因子计算 , 而插值因子根据通过单元边的通量比求得 , 即 :Fp = max min F1 /F2 , 1 , 0 , Fn = max min F4 /F3 , 1 , 0 X p = ( 1 - Fp ) X 2 + ( 1 - Fn ) X 4 + Fp Fn X 3Y p = ( 1 - Fp ) Y 2 + ( 1 - Fn ) Y 4 + Fp Fn Y 3( 18)( 19 a) ( 19 b)式中 , F1 、F2 、F3 、F4 为各单元边上的通流量 箭头所示方向为正 ( 参见图 1) 。同时 , 上游点的插值函数 N pj 可表示为 :N p1 = 1 , N p2 = 1 - Fp , N p3 = Fp Fn , N p4 = 1 - Fn( 19)这样 , 由 ( 16 a) 、( 16 b) 式即可求得到对流 - 扩散方程的对流项离散后的单元系数矩阵的元素 。需要指出的是 , 这里是用上 、下游结点间的直线段代替流线弧长的 , 这样的简化不会对9N j 32 山东轻工业学院学报第 13 卷对流项的计算精度产生大的影响 。213源项的离散为了加速迭代及解的稳定性 , 将各方程中的源项分为两部分 , 一部分与该方程中的变量有 关 , 另一部分与方程中的变量无关 , 分别记为 S p 、S c 。对其分别进行离散 , S p 离散后放在方程的左边 , S c 离散后放在方程的右边 。= S pijj + S ciS i参照表 1 可很容易地写出源项的离散形式 : 1 9W i 9N j 1 9V 9 P 9U 方程中 : S UUci- +d A , Sd A( 20)= -=pijRee 9XRee9X 9X9X9 Y( e)( e)AA 1 9 W i 9N j1 9U9 P9V 方程中 S V = -Vci= - +d A , Sd A( 21)pi j9 Y9 Y9 Y9XRee 9XRee( )e( e)AA L reK 方程中 : S KKciGW d AW d A , S( 22)= -=pijiiK U re( e)( e)AAL re方程中 : S pi j = - C2= C1 G K W i d AW i d A , S ci( 23)KUre( e)( e)AA214压力方程的推导及速度修正为了推出压力方程 , 首先从离散的动量方程得到速度和压力的关系 , 然后将其代入守恒型 连续性方程代替方程中的速度 , 从而得到求解压力的代数方程 。由离散的总体动量方程可以获得压力和速度之间的需要的关系式 。动量方程离散后可表示为 : 1 9VW i 9 9 P9 Y= - ai j U j -d AW i d A( 24 a)aii U i+Ree 9X9Xj( e)( e)AAW i 9 9 1 9U P9X( e)= - ai j V j -d AW i d A( 24 b)aii V i+Ree 9 Y9 Yj( e)AA速度可以表示为 :U i = i -m p 9 9 P P+ n , V =V - m+n( 25)u iip 9 Yv9Xai j U jai j V j其中 , i = - a a( 26), V i= -iji jjjm p = W i d A /aii( 27)( e)A 1 9V 1 9U9 99XW i d A/ai i , n vW i d A/aii ( 28)n u =9 YR ee 9XRee 9 Y( e)( e)AA式 ( 25) 就是速度和压力间的关系式 , 在结点处速度和压力间的关系不是严格精确的 , 而仅仅是一逼近 , 对单元完成上述积分 , 可得到下列形式的压力单元方程 : 9W i 9W i9 P9 P9XN j j -i =d A( 29)m j 9X+N j V j -m j 9 Y9 Y( e)A上述积分可用 2 2 高斯积分法求得 , 可得到下列形式的单元矩阵 : A P= F 。以常规方式将单元压力方程装配起来就可得到总体压力方程 , 其矩阵是正定的和对称的 。求解不会遇到特33第 2 期王仁人等 : 蜗壳内湍流流场的单向流线上风有限元法数值计算殊的问题 。215密度的修正方法密度由能量方程 ( 总焓不变) 求得 , 无量纲化后则有 :p+ 1 ( U 2 + V 2 )H0 =( 30) -1U 22re设流场中某点压力 、速度和密度分别为 : P0 、U 0 、V 0 、0 , 可得到流场中任一点的密度为p0 - 1 2( U 2 - U 2 ) + ( V 2 - V 2 ) = P / 0 -U re( 31)0023求解迭代方法及步骤首先假定一密度场 , 按湍流不可压缩流动迭代到基本收敛后 , 再用式 ( 31) 对密度场进行修正更新 , 然后继续迭代直至密度场收敛 。为了检验上述计算方法 , 首先计算了图 ( 2) 所示的后掠台阶流动 ( H /h = 115 , H 为流 通道的高度) , Reh = 42000 。计算区域划分为 405 个单元 , 448 个结点 。计算结果如图 ( 3) 、( 4) 所示 。可见 , 在台阶后存在明显的回流区 ; 沿 X 方向 , 速度的分布逐渐趋平坦 , 在 X /h = 10 截面上速度分布规律已非常接近直管内流动规律 ; 在台阶附近压力分布呈二维特性 , 在 X /h= 415 截面上压力分布已恢复均匀 。在整个计算过程中 , 迭代稳定 , 收敛较快 , 没有出现压力 振荡 , 这说明上述计算方法正确可行 。图 2 后掠台阶有限元网格划分图 3 后掠台阶流动流场矢量图 34 山东轻工业学院学报第 13 卷图 4 后掠台阶流动不同截面上 X 方向速度分布4蜗壳内流动计算结果及分析蜗壳内流动计算有限元网格的划分如图 ( 5) 所示 , 共划分为 611 个节点 , 564 个单元 , 图( 6) 为蜗壳截面 E 、H 和 L 上的通流速度分布 。可以看出 , 通流速度在蜗壳出口处最大 , 在同 一个截面上随着半径 R 的增大逐渐减小 , 其趋势是基本符合角动量守恒规律 ; 由截面 E 到截面 L , 同一半径处降低 。计算结果表明 , 可以用上述方法研究蜗壳几何形状对蜗壳内流动的影 响 。在相同的出口和进口速度分布下 , 改变蜗壳流道截面的收缩率 , 可以得到不同的通流速度 分布和出口气流角分布 , 根据设计要求对蜗壳的各种形状进行优选 。这对蜗壳设计工作是很 有用的 。图 5 蜗壳内流动有限元网格划分5 结论511采用单向流线上风 Galer kin 有限元法 , 并采用压力速度同次插值函数 , 求解二维可压缩湍流流 , 数值检验结果表明 , 在整个计算过程中 , 迭代稳定 、收敛较快 , 没有出现压力振荡 , 且具 有较高的精度 , 这说明上述计算方法正确可行 。512求解了形状比较复杂的变截面涡轮蜗壳内的二维可压缩湍流流动 , 计算结果表明 , 蜗壳 截面上的通流速度是基本符合角动量守恒规律 。513可以用上述方法研究蜗壳几何形状对蜗壳内流动的影响 , 从而可根据设计要求对蜗壳的35第 2 期王仁人等 :蜗壳内湍流流场的单向流线上风有限元法数值计算图 6 截面 E 、H 、L 上的通流速度分布各种形状进行优选 。参考文献1234王仁人. 变截面涡轮蜗壳内气体流动与增压器性能的研究.博士学位论文,上海交通大学 1996 年 11 月.周建华 ,二维三维流场的有限元模拟及一维非定常流计算模型的研究 ,博士学位论文,上海交通大学 ,1996 年 1 月. S. V . Patankar ,Numerical Heat Transfer and Fluid Flow , Hemisp here ,Mc Graw - Hill New Yo r k ,1980 .J Heirich O C Zienkiewicz. The Finite Element Met ho d and Up winding Techniques in t he Numerical Solutio n of Co nvectio n Do minated Flow Pro blems. In : T I R hughes , Ed. , Finite Element Met ho ds fo r Co nvectio n Do minated Flows ,AMD Vol . 34 ( A SM E ,New Yo r k ,1979) .T J Rice , R J Schnip ke . A Mo noto ne St ream line Up wind Finite Element Met ho d fo r Co nvectio n Do m