叶片外换热计算程序STAN5收敛性改进

娄德仓

(中国燃气涡轮研究院，四川 成都 610500)

1 引言

涡轮叶片燃气边的换热系数，可通过数值求解流过型面的粘性流控制微分方程得到。这种解法的优点是：可得到壁面附近气流的速度和温度分布，且在此基础上求得的壁面摩阻系数和换热系数较为准确；能较方便地应用近年来发展的湍流模型；计算中用到的经验公式少，适应性较广。由于求解边界层方程的计算量较求解N-S方程的小得多，因而在工程计算中得到了较广泛的应用。Crawford等人于1975年采用Patankar-Spalding方法编制了边界层微分方程计算程序STAN5[1]，随后又对该程序进行了诸多改进。目前使用的TAXSTAN程序是其最新版本[2]，能考虑叶片表面气膜出流对外换热的影响。

上世纪80年代，STAN5程序被西工大刘松龄教授等引入国内，经改进后用于涡轮叶片外换热计算，结合气膜修正等经验公式，能用于带气膜冷却涡轮叶片的热分析。邓化愚等人[3]还对程序进行了改进，主要是加入了湍流度间隙因子，使程序能考虑来流湍流度对外换热的影响。蔡毅等人[4]通过微分法、积分法等的对比，发现STAN5程序在计算方法及精度上更具优势。但在应用STAN5程序进行外换热计算时，会出现不收敛等问题，这主要是由外边界的逆压梯度分布、计算网格及边界层卷吸量的控制等因素引起。因此，计算中经常要对涡轮气动数据进行人工干预，这会造成计算结果的人为误差。为此，本文对STAN5程序计算中存在的收敛性等问题进行分析，并采取相应措施进行改进。

2 边界层微分方程

描述轴对称可压缩流的边界层流动方程，包括连续性方程、动量方程和能量方程。其中连续性方程可表示为[5]：

式中：坐标x、y的方向分别与物体的型线相切和垂直，u和v分别为x、y方向的分速度，ρ为气流密度，r为径向坐标。

动量方程为：

式中：p为压力，μ为动力粘性系数，ρ----u′v′为湍流剪切应力，u′和v′分别为x、y方向的湍流脉动速度，z为单位质量流体上的彻体力。

总焓形式的能量方程为：

式中：I为气体总焓，λ为气体导热系数，cp为气体定压比热，i为气体静焓；I′为总焓脉动值，S为单位质量流体中的源项。

根据以上方程，用已知的边界条件和初始条件，便可解出ρ、u、v、I在边界层内的分布。对于湍流边界层，还须补充求湍流剪切力的关系式。

3 STAN5程序求解

3.1 初始条件

求解边界层微分方程过程中，常将驻点区域层流边界层作为计算初始条件。在STAN5程序中，根据Pohlhausen公式获得层流边界层内的速度分布：

式中：δ为边界层厚度，Λ为Pohlhausen参数，ue为边界层外气流速度。

边界层起始点温度分布(T)可通过温度分布与速度分布间的线性关系得到，关系式为：

式中：Tw为壁面温度，Te为边界层外气流温度。

3.2 湍流模型

现有的STAN5程序中采用的湍流模型有混合长度模型、单方程湍流模型，同时考虑了湍流度修正。其中较为常用的是Prandtl混合长度模型。

3.3 边界层网格变换

边界层微分方程求解方法很多，STAN5程序采用的是Patankar和Spalding提出的流函数坐标变换法[6]，首先从x、y坐标变换为x、ω(x,y)坐标，对边界层微分方程进行求解，然后再将结果返回到x、y坐标下。ω(x,y)及其与y的关系可用下式表示：

STAN5程序中的网格，其径向分布采用基于式(7)建立的速度分布与网格密度关系。为保证边界层底层网格有足够的密度，采用了滑移点控制方法。

式中：c1为常数，β为滑移点参数。

4 STAN5在计算中存在的主要问题及改进

以某涡轮叶片外换热计算过程为例。STAN5程序在计算叶片外部换热时，首先将叶片表面的等熵马赫数分布作为边界层计算中的外部速度边界。通过求解边界层微分方程，获得叶片表面换热系数分布。典型的涡轮叶片表面等熵马赫数分布如图1所示。从图中看，在涡轮叶片叶盆侧前缘扩张段和尾缘激波区存在明显的逆压梯度，而这些位置通常会造成外换热计算不收敛。因此在外换热计算过程中，通常要对外部势流速度进行人为调整，避免出现逆压梯度(如图1中虚框内所示)。由于边界层微分方程为椭圆型微分方程(即下游计算结果完全受上游的影响)，因此该调整会导致计算结果人为误差较大，影响外换热计算结果的精确性。

图1 涡轮叶片表面等熵马赫数分布Fig.1 Isentropic Mach number distribution of blade surface

4.1 收敛性问题及分析

根据边界层内部流动方程的量级分析，边界层内部的压力分布可简化为式(8)。其中沿流向的压力梯度决定于外部主流的速度分布。因此，在外部势流存在较强逆压梯度且边界层底部速度较小时，计算中就会出现速度负值。而边界层网格又与速度分布相关，最终导致计算发散。

通常，边界层内部出现速度负值时，STAN5程序自身会对速度负值进行调整(如式(9))，同时减少边界层外部的卷吸流量，重新求解方程。

如图2所示，修正后的速度分布与边界层内真实的速度分布有很大差别。若继续计算，仍会出现速度负值，导致计算发散。有时虽能得到收敛结果，但所采用的方法已改变了边界层外部的气流速度分布规律(图3)，并且这种收敛方式会对换热系数(α)产生较大影响(图4)。

4.2 收敛性的两点改进

通常在边界层方程计算中，若出现速度负值，表明边界层已破坏，这时边界层方程不再适用。程序自身的修正实质上是假设该区域内仍然为边界层型流动。在涡轮流道中，大部分区域属于顺压梯度的边界层流动，因此该局部区域内采用此假设仍然合理。本文在此基础上对逆压梯度下的速度分布u(y)及其网格密度进行修正。

4.2.1 逆压梯度下的速度分布

针对叶片外部速度场急速变化时产生的较强逆压梯度，可采取两种方法进行速度修正。

一是在当前计算站速度梯度出现负值时，对速度分布进行修正，如图5(a)所示。当前计算点xi的速度分布形式ui(y)，与前一计算点xi-1的速度分布形式ui-1(y)保持相似，即：

图5 边界层内部速度分布修正Fig.5 Modification of the velocity profile inside the boundary layer

结果表明该方法不能完全避免速度负值问题。

二是在当前计算站存在逆压梯度且速度二次梯度出现负值时，对速度分布进行修正，如图5(b)所示。当前计算点xi的速度分布形式ui(y)，与前一计算点xi-1的速度分布形式ui-1(y)保持相似，并保持外边界的速度值uei不变，即：

结果表明该方法能明显改善计算不收敛。

4.2.2 基于压力梯度变化的自适应网格

叶盆侧尾缘激波区马赫数变化剧烈，压力梯度迅速由顺压梯度变为逆压梯度，因此该区域对流向网格变化较敏感。虽采取了上述速度分布修正，但在有些马赫数分布下仍会出现计算不收敛。通常边界层在尾缘区域较厚，沿流向的空间步长较大(Δx=0.3δ)，这时会出现两种情况：一是当前节点位于最大(或较大)逆压梯度点，导致计算发散；二是当前节点跨过最大(或较大)逆压梯度点，计算稳定收敛。

为使边界层内部沿流向网格步长与主流压力梯度相联系，在逆压梯度下，采用式(12)计算网格步长：

式中：下标old表示上一步长的参数，Δx0为当前计算站下的初始步长，Δx1为经压力梯度修正后的计算步长。图6对比了压力梯度修正前后叶片尾缘附近的网格，可见修正后的网格明显变密。

经上述修正，STAN5程序计算中的不收敛现象得到解决。修正过程中，保持外边界速度不变，就能保证计算结果误差在合理范围内。

图6 修正前后的沿流向网格Fig.6 Comparison of the stream-wise mesh before and after modification

5 计算结果分析

图7 外换热计算结果分布Fig.7 External heat transfer coefficient of the blade

根据图1中的马赫数分布，改进后的STAN5程序外换热计算结果与人工修正马赫数分布计算结果间的对比如图7所示。从图中看出，两者间的外换热分布趋势一致，但在叶盆侧，两者的计算值相差约±4%，最大达20%。采用新型网格的计算结果中，对尾缘激波区域(如图中叶背尾缘)的换热系数分布更加详细。图8为边界层内部的速度分布云图。

图8 边界层内部的速度分布云图Fig.8 Contour of velocity inside the boundary layer

6 结论

本文分析了涡轮叶片外换热计算程序STAN5中的求解方法，包括边界层求解的微分方程、计算初始条件及网格变换等。针对程序在外流场存在较强逆压梯度时所出现的计算不收敛及错误收敛等现象，对该程序进行了两点改进：修正了逆压梯度下边界层内部速度分布，建立了边界层内部沿流向的网格分布与主流压力梯度间的联系。改进后的程序有效地避免了计算发散等问题。算例分析结果表明，本文提出的改进方法与以前的人工调整方式获得的结果相比，外换热系数分布趋势一致，局部误差约为4%，最大误差达20%。

[1]Crawford M E，Kays W M.STAN5-A Program for Numeri⁃cal Computation of Two-Dimensional Internal and Exter⁃nal Boundary Layer Flows[R].NASA CR-2742，1976.

[2]Kays W M，Crawford M E，Weigand B.Convective Heat and Mass Transfer[M].4thed.McGraw-Hill，2005.

[3]邓化愚，蔡 毅，刘玉芳.涡轮气冷叶片局部外换热系数计算方法的改进[J].燃气涡轮试验与研究，1989，2(2)：1—5.

[4]蔡 毅，刘松龄.涡轮叶片外换热系数计算方法和比较[J].燃气涡轮试验与研究，1995，8(1)：1—8.

[5]航空发动机手册总编委会.航空发动机手册：第16册——空气系统及传热分析[K].北京：航空工业出版社，2001.

[6]Patankar S V，Spalding D B.Heat and Mass Transfer in Boundary Layers:A General calculation Procedure[M].London：International Textbook Company Ltd.，1970.