转子三维反问题气动设计

收稿日期: 2001- 04- 10; 修订日期: 2001- 06- 04 基金项目: 国家自然科学基金资助项目 (59876041) 作者简介: 陈乃兴 (1933- ) , 男, 中国科学院工程热物理研究所研究员. 第 17 卷　第 1 期 2002 年 1 月 航空动力学报 Journa l of Aerospace Power V o l117 N o11 Jan. 　2002 文章编号: 100028055 (2002) 012023206 单转子风扇的三维反问题气动 陈乃兴1, 徐燕骥1, 黄伟光1, 陈俊杰2, 陈晓东2 (1. 中国科学院 工程热物理研究所, 北京 100080; 2. 奥克兰大学, 新西兰 奥克兰) 摘要: 采用准三维流函数反问题设计方法与三维N 2S 求解方法的相互迭代对单转子风扇进行气动设计。在 叶片初步计算中得到叶栅进出口气流角沿径向分布, 并将它作为本文的目标函数。采用准三维反问题求解方 法, 依次构造出各个 S 1 流面上的叶片几何形状和气流角分布。然后再采用N 2S 方程的求解方法, 对叶片进行 全三维流场的数值计算。通过N 2S 方程计算结果与目标函数的对比, 重新修正叶片出口气流角分布, 并作为 下一次反问题设计的目标函数。经过反问题与N 2S 方程求解的反复迭代, 最终得到满足设计要求的叶型。 关　键　词: 风扇设计; 反问题; 流函数; N 2S 求解; 叶轮机械气动热力学 中图分类号: V 23113　　　　　文献标识码: A An Inverse M ethod for 3-D Aerodynam ic D esign of Fan s and Pum ps CH EN N ai2x ing1, XU Yan2ji1, HU AN G W ei2guang1 John J J Chen2, CH EN X iao2dong2 (1. Ch inese A cadem y of Sciences,Beijing 100080, Ch ina; 2. T he A uck land U n iversity,A uck land P riva te Bag 92019,N ew Zealand) Abstract: T he m ethod fo r 32D design of fan s o r pum p s in the p resen t paper is p resen ted w h ich is based on the so lu t ion s of a quasi23D inverse m ethod and a 3D N 2S so lver. A t beg inn ing the flow ra te, the p ressu re ra t io and the geom etry of m erid ian pat tern w ere g iven. T he span2heigh t d ist ribu t ion s of in let and ex it p itch2w ise2averaged flow angles (o r b lade circu la t ion) of the b lade channel w ere ob ta ined by the p relim inary ca lcu la t ion s. T hey are served to be the ta rget param eters of the p resen t ca lcu la t ion. Successively u sing a quasi23D inverse m ethod and a N 2S so lver, a 32D b lade configu ra t ion can be ob ta ined. A s an exam p le, an ax ia l sing le ro to r fan is show n. Key words: fan design; inverse p rob lem ; stream funct ion; N 2S so lu t ion; aero therm odynam ic of tu rbom ach inery 1　前　言 　　随着计算机科学和计算技术的发展, 计算流 体力学在叶轮机械的产品设计中得到了广泛的应 用。众所周知, 一般的反问题和杂交问题的最大特 点, 是在给定压力面和吸力面上的压力分布或速 度分布条件下, 能直接得到叶片的几何形状, 其不 足之处在于很难给定优化的叶片压力分布, 无法 考虑机械强度的要求 (如叶片厚度分布等) , 并且 有时会得到一个非物理解。若在设计计算中引入 © 1994-2006 China Academic Journal Electronic Publishing House. All rights reserved. http://www.cnki.net 叶片负荷分布和厚度分布, 则可消除以上缺陷。 本文提出一个由流函数反问题与N 2S 算法 相结合的气动设计体系, 通过这两类不同算法的 相互迭代, 最终得到满足进出口气流角分布 (或环 量)、流量、总压比等设计参数要求的三维叶片。本 文计算体系中所采用的目标函数为叶片气流角 (或环量)。准三维流函数反问题造型的目标函数 是沿子午面坐标线的相对环量分布和叶片厚度, 在计算中, 通过适当调整叶片前后缘附近的环量 分布, 使得进出口气流角满足无冲角条件和 Ku tta2Joukow sk i 条件, 并且获得优化的设计叶 型。 由于篇幅所限, 本文只对某一单转子风扇进 行气动设计。计算结果表明, 本文所提出的方法可 按照设计者的不同设计要求, 如进出口气流角分 布或负荷分布, 自动构造出不同半径上的叶型, 并 且可获得详细的流场分布以及总压比、总温比、绝 热效率等气动性能参数。 2　计算方法 2. 1　基本思路 　　本文叶片设计的基本思路是: 即通过准三维 反问题和三维N 2S 求解方法的相互迭代, 来满足 叶片进出口沿径向气流角的分布规律, 其目标函 数为叶片气流角 (或环量)。计算过程中, 当然可采 用任何一种反问题和N 2S 解法。 2. 2　准三维流函数反问题求解 　　如上所述, 本文采用的设计方法是准三维反 问题与三维N 2S 解法的相互迭代的算法体系, 通 过反问题与N 2S 解法的反复迭代, 来获得满足叶 片负荷 (环量或气流角)的三维叶型。 本文所用的流函数反问题求解中将叶片环量 作为目标函数之一。环量表示叶片对流体所作的 功, 叶片总环量可表示成如下形式:# 0 = - (W rsinΒ) i∆Υ+ (W rsinΒ) e∆Υ+ (8 r2 sinΡ) i∆Υ- (8 r2 sinΡ) e∆Υ (1) 其中下标 i 和 e 表示叶片进出口值; W 为相对速 度, 与绝对速度的关系为W∼ = Vψ- 8ψ×rο; 8 为角速 度; Β和 Ρ分别表示气流角和子午面角; r 和 ∆Υ分 别为半径和一个叶片通道内的周向节距。 从上式不难看出, 对于给定的叶片几何参数 和转速, 存在着唯一的一个与进口气流角相对应 的总环量。叶片总环量也可由下式表示, 即: # 0 = 6t. e. l. e. d# = ∫t. e. l. e. d# d l d l (2) d# d l = - W sco sΒs + W p co sΒp (3) 式中 # 表示从叶片前缘点到计算站的当地环量; d# öd l 为当地环量的一阶导数; 下标 l. e. 和 t. e. 分别表示叶片的前缘点和尾缘点; s 和 p 分别表 示吸力面和压力面。 本文采用相对环量, 并作为目标参数分布, 其 表达式如下: 若 0 ≤ xθ ≤ a1, 则# ö# 0 = ca2xθ2 + ca3xθ3 (4a) 若 a1 ≤ xθ ≤ a2, 则## 0 = b1 + a1 (xθ - a2) b2 - b1a2 - a1 (4b) 　　在式 (4a)和 (4b)中 ca2 = g am 1 a1 + 3b1 a 2 1 和 ca3 = g am 1 a 2 1 - 2b1 a 3 1 若 a2 ≤ xθ ≤ 1, 则# ö# 0 = ca0 + ca1 (xθ - a2) + ca2 (xθ - a2) 2 + ca3 2 (xθ - a2) 3 (4c) 其中: ca0 = b2　, 　ca1 = g am 1 ca2 = 015 - g am 1 - 3ca3 (1 - a2) 21 - a2 和 ca3 = g am 1 (1 - a2) - 2 (1 - b2)(1 - a2) 3 (5) 式 (4a)～ (4b) 中的系数可采用L uu [ 1 ]的推荐值, 即: 　a1= 011　, 　b1= 0107　, 　a1= 019 b1 = 0193 , g am 1 = (b2 - b1) ö(a2 - a1) (6) 对于本文的算例, 通过流函数方程的求解, 我们可 以得到叶片的环量沿流向分布。叶片进出口气流 角变化量则可由环量求出。 本文所采用的反问题计算所采用的另一个目 标参数是沿叶片子午面方向和叶高方向的叶片厚 度分布, 这是考虑到叶片的强度要求。 2. 3　流函数方程 　　每一个回转流面上的叶型可由下面的流函数 方程求出, 即: A 1 52Υ57{ 2+ A 2 55l 5Υ57{ + A 3 52Υ5l2 + A 4 5Υ57{ + A 5 5Υ5l = A 6 (7) 式中的所有系数以及方程的求解过程详见文献 [4, 5 ]。只要将沿叶高方向所有坐标线上的任意回 转面进行求解完, 就可得到 1 个三维的叶片和一 42 航 空 动 力 学 报 第　17　卷 © 1994-2006 China Academic Journal Electronic Publishing House. All rights reserved. http://www.cnki.net 个为求解三维N 2S 方程用的三维计网。 2. 4　N-S 方程求解 　　本文应用压力修正法求解三维N 2S 方程, 湍 流模型为 k2Ε模型。为适应复杂叶片通道的几何 形状, 本文采用坐标变换, 将直角坐标的N 2S 方 程变换成任意曲线坐标形式, 并引入逆变速度分 量以简化方程的表述形式, 速度求解变量仍为直 角坐标系下的相对速度分量。本算法的详细推导 以及相应算例的验证, 请详见参考文献[8 ]。 任意曲线坐标系下的控制方程可表示如下:5(Θ g 5 )5t + 5(ΘW 15 )5x 1 + 5(ΘW 25 )5x 2 +5(ΘW 35 )5x 3 = 5 J # 5 g 11 555x 15x 1 +5 J # 5 g 22 555x 25x 2 + 5 J # 5 g 33 555x 35x 3 + S i (x 1, x 2, x 3) (8) 其中W 1,W 2,W 3 分别为相对逆变速度分量, 可表 示为: W 1 = J w x 5x 1dx + w y 5x 15y + w z 5x 15z 　 W 2 = J w x 5x 2dx + w y 5x 25y + w z 5x 25z W 3 = J w x 5x 3dx + w y 5x 35y + w z 5x 35z (9) 式中 J 为雅克比行列式; S i (x 1, x 2, x 3)为源项。任 意变量和有效扩散系数可写成:5 = 1,w x ,w y ,w z , I , k , Ε (10)# 5 = 0, Λef f , Λef f , Λef f , Κef f cp , ΛtΡk , ΛtΡΕ (11) 其中: w x , w y , w z , I , k , Ε分别为求解变量, 即直角 坐标系下的相对速度分量、转焓、湍流动能以及扩 散率; Λef f = Λl+ Λt 和 Κef f = Κl+ Κt 分别为有效粘性 系数和有效导热系数; 湍流粘性系数可由Λt = C ΛΘk 2öΕ求出; Ρk , ΡΕ 分别为经验系数, 取为 110 和 113; C 1, C 2, C Λ 分别为湍流模型的经验系 数, 取为 1147, 1192 和 0109。 本文计算中, 采用压力修正法求解 N 2S 方 程。控制方程的离散化时, 对流项采用迎风格式, 扩散项采用中心差分格式, 并应用线松弛迭代法 求解。 3　计算算例 　　为了验证本文提出方法的精度和可靠性, 本 文对一个低速轴流风扇进行了气动设计。该风扇 的流量为 19 kgös, 总压比为 11015, 转速为 1200 röm in。进口气流角可依据零攻角原则估算出:Β1 = - tan (8 r in letöV x , in let) (12) 其中 8 , V 分别表示转速和绝对速度; r 为半径; 下标 x 和 in let 分别表示轴向和进口参数。 图 1　目标相对气流角沿叶高的变化 图 2　反问题初次近似所的叶片形状 叶片进出口气流角的径向分布可由初步设计 计算获得, 并作为本文的目标函数, 如图 1 所示。 在同一半径处, 进出口气流角的变化量为 20°。图 2 为第一次反问题求解出的叶片根部、中部和顶 部叶型。所得的初步的三维网格如图 3 所示。气 流角分布可由N 2S 方程求出, 如图 4 所示。 3. 1　迭代过程 　　从图 4 可看出, 计算出的气流角与目标函数 存在着一定的差别, 因此需要对气流角分布进行 52第　1　期 陈乃兴等: 单转子风扇的三维反问题气动设计 © 1994-2006 China Academic Journal Electronic Publishing House. All rights reserved. http://www.cnki.net 图 3　由问题求解所得的三维计算网格 图 4　由N 2S 方程求解所得的相对气流角 图 5　逐次近似所得的相对气流角分布 重新修正。首先根据修正后的气流角, 应用反问题 求解方法重新生成三维叶型和网格; 然后采用N 2 S 算法对新生成的叶型重新计算。通过对计算值 与目标函数再次比较与分析, 确定下一次迭代的 修正值。经上述过程的反复迭代, 直至收敛, 并最 终得到满足设计要求的叶型。气流角分布的收敛 过程如图 5 所示。对于本文的算例, 整个迭代次数 仅需 3 次。迭代时采取类似数值求解的原则。 图 6　初次近似和最终解所得的根部叶型 图 7　初次近似和最终解所得的中部叶型 3. 2　计算结果 　　由于本文采用了N 2S 求解方法, 因此我们就 可得到详细的气动参数分布, 如周向平均的径向 分布, 不同回转面上的等参数等值线图以及总参 数分布等。图 6～ 图 8 分别为第一次和第三次迭 代后所生成的叶片根部、中部和顶部叶型。计算得 到的速度分量和绝热效率的径向分布如图 9 和图 62 航 空 动 力 学 报 第　17　卷 © 1994-2006 China Academic Journal Electronic Publishing House. All rights reserved. http://www.cnki.net 图 8　初次近似和最终解所得的顶部叶型 图 9　叶片前后的速度分量沿叶高的分布 图 10　风扇的绝热效率沿叶高的分布 图 11　叶中任意回转面的静压分布 图 12　压力面、吸力面和流道中心面上的静压分布 图 13　迭代的收敛历史 72第　1　期 陈乃兴等: 单转子风扇的三维反问题气动设计 © 1994-2006 China Academic Journal Electronic Publishing House. All rights reserved. http://www.cnki.net 10 所示。为节省篇幅文中只给出少量的细节, 如 图 11 的叶片中部回转面上的静压等值图, 图 12 的 S 2 面上的静压等值线图。图 13 为N 2S 方程求 解的收敛史。 3. 3　风扇特性 　　本文所设计的风扇特性曲线如图 14 和图 15 所示, 分别表示总压比、总温比以及绝热效率与流 量的关系。这对设备的正常运行非常重要。 图 14　总温比和总压比与流量之间的关系 图 15　绝热效率与流量之间的关系 4　结　论 　　本文提出一种反问题求解与N 2S 方程算法 相结合的气动设计方法。该方法将反问题与N 2S 正问题解有效地结合起来, 首先采用反问题求解 以构造出三维叶型和网格, 然后应用N 2S 算法计 算出叶片气动性能参数, 通过两种算法反复迭代 最终得到所需的叶型。本文设计方法具有很好的 实用价值, 特别是在反问题与N 2S 算法的选取上 具有很大的灵活性。本方法计算稳定性好、速度 快, 对于风扇算例, 仅需 3 次迭代, 即可获得满足 气动与强度设计要求的三维叶型, 因此本方法可 作为叶轮机械气动设计的一种有效工具。 致谢　本文工作得到国家重点基础研究项目 973 与国家自然科学基金项目的资助, 作者在此 表示深深的谢意。我们对于新西兰科学与技术研 究基金的部分资助和奥克兰大学给予第一作者的 计算机的使用, 也表示衷心的谢意。 参考文献: [ 1 ]　L uu T S, V iney B, Bencherif L. Inverse P roblem U sing S 1 and S 2 App roach fo r the Design of the Turbom achine w ith Sp litter B lades [ J ]. Revue F rancaise de M ecanique, 1992, (3). [ 2 ] 　 Chen N ai2xing, L i W ei2hong. A M ethod fo r So lving A erodynam ic H ybric P roblem of P rofile Cascade on S 1 Stream Surface of Revo lu tion by Emp loying Stream 2Function Equation Exp ressed w ith N on - O rthogonal Coo rdinate System [ A ]. P roc. In ternational Conference on Inverse D esign Concep ts in Engineering Sciences ( IC ID ES ) [ C ], Edited by G. S. D ulik ravich, 1984. [ 3 ]　Chen N ai2xing, Zhang Feng2xian, L i W ei2hong. A n Inverse (D esign) P roblem So lu tion M ethod fo r B lade Cascade F low on Stream Surface of Revo lu tion [ J ]. Journal of T urbom ach inery, T rans. A SM E, 1986, 108 (2) : 194- 199. [ 4 ]　Chen N ai2X ing, John J J , Chen X D. A SFC Inverse M ethod w ith the D istribu tion of T arget C ircu lation o r its F irst D erivative fo r T urbom ach ine D esign[A ]. P roceeding, the 7th In ternational Ro tational M ach inery Conference [ C ], Feb. , 1998: 1386- 1393. [ 5 ]　Chen N ai2xing. N um erical M ethods fo r Inverse So lu tion in A erodynam ic D esign of T urbom ach inery [ J ]. Recen t D evelopm ent of A erodynam ic D esign M ethodo logies2Inverse D esign and Op tim ization - , Edited by K. Fujii and G. S. D ulik ravich, F riedr. V iew eg and Sohn V erlagsgesellschaft m bH , B raunschw eigöW iesbaden, 1999: 113- 145. [6 ]　Chen N ai2xing, Xu Yan2ji. N um erical Computation of 32D T urbulen t F low in T urb ine Cascade [ J ]. Ch inese Sciences Bullet in, 1990, 35 (14). [ 7 ]　Chen N ai2xing, Xu Yan2ji. N um erical So lu tion of 22D and 32 D T urbulen t In ternal F low P roblem s [ J ]. In t. Journal of N um ericalM ethods in F lu ids, 1991, 13 (4). [ 8 ]　Chen N ai2xing, Xu Yan2ji, H uang W G, M a H W , J iang H K. N um erical Study on 32D T urbulen t F low in a L ow 2Speed A xial Single Ro to r Comp resso r and its Comparison w ith Experim ent [R ]. Contribu ted Paper, N o. 009, P roceeding of the Fourth In ternational Sympo sium on Experim ental and Computational A ero thermodynam ics of In ternal F low , A ug. 31- Sep t. 2, 1999. (责任编辑　杨再荣) 82 航 空 动 力 学 报 第　17　卷 © 1994-2006 China Academic Journal Electronic Publishing House. All rights reserved. http://www.cnki.net