湍流工况下水导轴承承载力及其Jacobi矩阵的一维直接算法

湍流工况下水导轴承承载力及其Jacobi矩阵的一维直接算法 湍流工况下水导轴承承载力 及其 Ja c ob i矩阵的一维直接算法 闫志勇 , 陆 圯 , 郑铁生 ()复旦大学 力学与工程科学系 ,上海 200433 摘 要 : 基于变分原理 ,把湍流工况下水导轴承修正的 R eyno ld s方程转化为一个泛函极值问题 ,采用权函数 和线性插值函数乘积作为压力函数 ,提出计算承载力及其 J acob i矩阵的一维有限元直接求解算法 。算例表明 ,算 法不仅结果精度能得到保证 ,而且运算时间可以大为缩短 。 中图分类号 : O32; O353. 1; TM 312 文献标识码 : A DO I编码 : 10. 3969 / . j issn. 1006 - 1355. 2009. 06. 076 关键词 : 振动与波 ;水导轴承 ;湍流 ;变分原理 ;有限元 D irec t So lu t ion M e thod for L oa d Force an d Ja cob i M a tr ix of W a ter Gu ide Bea r in g in Turbu len t Reg im e YAN Z h i2yong, LU Y i, ZH EN G T ie2sheng ()D ep a rtm en t of M echan ic s and Enginee ring Sc ience, Fudan U n ive rsity, Shangha i 200433 , Ch ina A b stra c t: B a sed on the va ria tiona l theo ry, the am ended R eyno ld s equa tion of wa te r gu ide bea ring in tu rbu len t regim e is tran sfo rm ed to func tiona l extrem e va lue p rob lem. U sing the p roduc ts of we igh t func tion and linea r in te rpo la tion func tion a s p re ssu re func tion, one2d im en siona l fin ite e lem en t a lgo rithm fo r ca lcu la ting load fo rce and its J acob i m a trix is p ropo sed. The num e rica l re su lts show tha t th is m e thod can reduce the comp u ting tim e rem a rkab ly w ith the p rec ision gua ran teed. Key word s: vib ra tion and wave; wa te r gu ide bea ring; tu rbu len t; va ria tiona l p rinc ip le; fin ite e le2 m en t m e thod 水导轴承广泛地应用于大型立式水泵机组 、核 , 计 算量 巨 大 , 通常 采用 建 计算两个湍流润滑 系数 [ 3 ] 数据库方法 。戴学余等利用 N g2Pan理论 ,建立了 泵上 。由于制造 、安装误差 、来流等因素 ,运行时 低粘度非油润滑滑动轴承运动方程 ,并通过五点离 轮和主轴承受径向力作用 ,水导轴承起着稳定主 [ 4 ] 散差分格式对其进行求解 。 和叶轮的作用 。 综上所述 ,实际工程中水导轴承紊流计算模式 对水导轴承而言 ,由于采用水这种低粘度润滑 的求解 ,往往采取有限差分法 、有限体积法 、数据库 ,雷诺数一般较大 , 层流 R eyno ld s润滑方程往往 等方法 ,计算量巨大 。尤其是对水导轴承支撑的水 准确 ,需要采用紊流理论加以修正 。N g和 Pan根 轮机等机械系统进行非线性稳定性分析时 ,在求解 Pe icha rd t壁面定律提出了湍流润滑方程 ,适应于 系统的动力响应的每一步都必须重复计算承载力 。 [ 1 ] Coue tte流为主的高速轻载轴承 。张运清在考 惯郑铁生等运用变分原理 ,把流体润滑 R eyno ld s方程 性效应的情况下 ,提出了一种适用于高压密封 高转化为一个自由边值问题 ,并进而对此自由边值问 [ 2 ] 题等价变分原理进行了系统的归类和证明 ,为滑动 速重载轴承的复杂流场紊流 R eyno ld s方程 。 [ 5 ] 轴承非 线 性 油 膜 力 的 计 算 提 供 了 一 种 新 途 径 。 于需要采用修正平均法对惯性项 进行 处 理以 及 稿日期 : 2009 - 05 - 20; 修改日期 : 2009 - 07 - 28 本文基于自由边值理论和 N g2Pan紊流润滑理论 ,把 ( ) 者简介 : 闫志勇 1983 - ,男 ,河南南阳人 ,博士 ,主要从事转子动 湍流工况下水导轴 承满 足的 修 正 R eyno ld s方程 转 力学 、非线性动力学研究 。 化为一个泛函极值问题 ,选取权函数与线性插值函 E2m a il: 081029009 @ fudan. edu. cn 数乘积的形式作为压力函数 ,建立对应的一维有限 其中 9 99P 9P 元变分方程 ,并采用直接法求解承载力及其 Jacob i 33 += ( )L P G HG H θ z θζ9 9 θζ9 9 矩阵 。 ? ? 9H 9H ( ) θ ( ) θ = 2 x + y co s+ 2 y _x sinf = - - θ 99t 1 基本原理 1 1 (Ω ) (Ω ) ΩH= { q?Cq = 0 on 9} , K = { q? 取 0 1 湍流工况下 ,水导轴承满足如下修正的非稳态 (Ω) H, andq Ε 0 } , 根据自由边值理论 , 可得无量纲 0 [ 5 ] R eyno ld s方程的等价变分原理 为 R eyno ld s方程 3 3 9 Gh1 9 Gh θ 9p z 9p ( )( ) ( ) + = 4 J P = m inJ q2 θ9 q?K 9z μ θμ R 9 9z 其中 ω 9h 9h ( )1 - - 1 θ 2 99t ( ) ( )( )5 ( ) a q, q- f qJ q= 2 其中 G, G分别为无量纲湍流系数 。 θz 39p 9q 9q 9qθζ ( )h= dda p, q 如果水导轴承内以 Po iseu ille流为主 ,无量纲湍 + G G θ ζ ???θ θ ζ ζ9999 Ω 流润滑系数取为 ) ( θζ( ) f q= f?qdd4. 726 5 ???Ω > 12 344 R0. 428 6 ep R G= epp 2 一维变分方程求解 1 / 12 < 12 344 R ep 223 ρ9p 9p h 1 Ω 如图 1 所示 , 润滑区域 为一矩形 , 其中的水 + 。对于由剪切 其中 R= ep 2 φμ9 9z r 膜破裂边界曲线虽然微微弯曲 , 但其上的压力及压 引起的 Coue tte流 1 力导数皆为零 。因此将其设为直线 , 对总体承载力 = Gθ0. 90c 12 + 0. 0136R e 的影响微乎其微 。考虑 到 水导 轴承 在 边界 上压 力 1 (θζ) (ζ) (θ) (ζ) q ,= g r + a ( )6 G= zc0. 9612 + 0. 0043R e 其中 不为零 , 则压力函数 q可以分离变量为 ωρμ其中 R e =R c /。对于 Po iseu ille流和 Coue tte λ) ( λ) ( g = g - = 0 流两种流动同时存在的情况 ,取 P+ PP- P 1 2 1 2(ζ) ζ a =+ )( G= m in G, G λθ θ 2 p c 2( )2 (θ ζ) γλλ ×= [ 0,] ×[ - ,] )( G= m in G, G z p zc 满足压力边界条件 。令 λ γ 23 2 9r θ(ζ) ζ d= g d= Ghd c θ 11? ? θ9 λ -λ 0 γ 23 2 9g ζ(θ) θd = Gh r dc= d ζ 22? ? ζ9λ 0 -λ γ θ) θ((ζ) ζ= f? r dc= g d d 33 ? ? λ0 - γ 图 1 水导轴承单块瓦坐标系及侧面展开图 λ P- P 3 1 2Fig. 1 The coordinate system of single pad of water guide bearing (θ) θ h r d= Gd (ζ) ζζc= g ′d 44? λ ? 0 λ- 图 1 为水导轴承单块瓦坐标系及侧面展开图 , ( ) 把上式代入泛函 5 , 可得 π其中 O 为瓦块中心 , C 为轴颈中心 。根据 定理 , 1 ( ) ( ) J q= c d + c d + c d - c d( )7 1 1 2 2 4 4 3 3 取 R, C为 特 征 量 , 无 量 纲 后 的 修 正R eyno ld s 方 f r 2 θ)(程为 假定 r 已知 , 则 d, d, d, d已知 。根据 Eu2 1 2 3 4 + ( ) ΩL P = f in le r2L agrange方程可得 2 9P 0 9 g ( )Γ3 (ζ) P = 0, = 0 on d g - d= d 1 2 3 29n ( )8 ζ9 0 Ωλ) ( λ) (P = 0 in g = g - = 0 ( ) J r9d 1 由最小势能原理可知 , 可得令 k = ,可得 9U d 2 KU = b ( )16 (ζ) λ)ζ)( )( ( 9 g = co sh k- co sh k 其中 U Ε 0, b Ε 0。 ( ) 代入式 6 ,可得 根据 K为对称 、正定 、三对角矩阵的特点 , 可以 ζ) (θ)( - co sh k] r (θζ) λ) ( q ,= [ co sh kP+ [ 6 ] 1 利用修正的追赶法 进行求解 , 运算量小且简单易 P- P P2 21 ( )ζ 10 + + 行 。利用追赶法得到待定系数 a 后 , 即可得到迭代 λ2 2参数 k。进而可以得到瓦块局部坐标系下的瞬态油 压力分布函数的表达式 。 膜力 f, f及其 J acob i矩阵 x y( ) 由式 9 ,可得 T T fa bba x11 3 ( )17 = - c = - c 3λ)3( - sinh 2 kλ λλ)T T ( c= 2+co sh 2 k 1 f 2 ky b2 a a b 1 2 9f9f k xx9aT T 9a λ λ)( c= [ sinh 2 k- 2 k]2 bb 11 2 9y 9y 9x 9x J = = - C 13 9a9fT T 9a 2 9f yyλ)( λλ)bb( - sinh kc= 2co sh k 223 9x 9y k 9x 9y ( )18 9f 9f c= 0 T 9a 9aTxx4 b b 1 1 ?? ?? 9x 9y( ) 9x 泛函 6 可以简化为 9y J C= = - 23 9f9f 9a 9a1 T Tyy ( )) ) ( + c d - c d 11 J 9 q= c db b 1 1 2 2 3 3 2 2 ? ? ? ? 2 9x 9y 9 9 x y ? ? 由上可知 ,设定初始迭代值 k, 可以计算出 c, 1 ( ) 其中式 18 中 a 对 x, y, x , y 的偏导数可以由 ( )(θ) , c, 求解泛函 11 的极值问题 , 可以得到 r , 修 3 ( ) ( ) 通过求解 方 程 19 得 到 。方 程 19 可 以 和 方 程 参数 k。如此循环迭代至满意精度 , 从而求得压 ( )16 利用修正的追赶法一同求解 , 相关矩阵和向量 分布函数 q, 继而可以根据下式求得部分轴瓦的 只需要在正压力区内求解 。迭代过程中 , 由于充分 膜力分量 考虑了空穴位置 , R eyno ld s边界条件可以不需迭代 (θζ) θθζf= - q ,co sddx ???计算而直接满足 ,避免了大量的重复计算 。 Ω ( ) 9a 9 K 12 K = - Cb- a 3 2 9x 9x (θζ) θθζf= - q ,sinddy ???Ω 9a 9 K K= Cb-a 3 1 ( )( ) 9x9y压力分布函数表达式 10 代入式 12 , 整理可得 ( )19 γ 9a γθ) θθ (sin Cf= - cr co sd-λ x 3 = 2CbK 3 1? 0? 9x ( )13 γ 9a γ)(θ) θθ ( f= - cr sind-C1 - co sλ y 3 K= 2C b 0?3 2 ? 9y 其中 ( )λC= P+ P λ 1 2 3 计算算例 (θ)用分段线性插值 , 将 r 离散化为 T (作为比较 ,同 时采 用 本文 算法 周向 150 个 单 ( )θ) θ)14 ((r = U L T (θ)(θ) (θ) (θ) 其中 L { l, l, ?, l} 为总体 = 1 2 n )(元 和二维有限元算法 采用 2 ×3 等参单元 , 周向 T U = { u, u, ?, u} 为待定系数向 值函数向量 , 1 2 n )150 个单元 ,轴向 6 个单元 计算圆瓦水导轴承湍流 ( )。代入泛函 11 可得 1 T T 模型下承载力 ,其中圆瓦水导轴承具体参数如表 1 () ( )= U KU- U bJ U 15 2 22 法计算得到的承载力曲线 ,其中 f = f+ f, 轴颈中 其中 所示 。 x y γ γ 心位移和速度分别为 图 2为圆瓦水导轴承湍流模型利用上述两种方 3 T 3 T (θ) (θ) θ (θ) (θ) θ = cGh L ′L ′d+ cGh L Ldθ ζ 1 2 π?? ? ξu = - co s 0 0 u = 0 2 γ ? πv = 0 θ) (θ) θ(b = cf L d3 ξv = - sin ? 0 2 从图 2可以看出 , 即使在轴颈位移大扰动的情况下 , 图 3为利用本文算法计算得到的圆瓦水导轴承 本文算法可以保持比较高的计算精度 。此外 , 本文 主刚度系数 、交叉刚度系数 。从图 3 可以看出 , x 方 算法计 算 时 间 约 为 二 维 有 限 元 算 法 计 算 时 间 的 向交叉刚度系数大于主刚度系数 , y 方向也存在比 较大的交叉刚度系数 , 可能造成水导轴承 - 转子系 1 / 10。 表 1 圆瓦水导轴承具体计算参数统颤振失稳 。 Ta b. 1 pa ram e ter s of c ircu la r wa ter gu ide bea r in g 4 结语 名称 参数值 名称 参数值 2本文利用 变分 法来 求 解湍 流工 况 下水 导轴 承 入口压力 轴承长径比 1. 0 8 Kg / cm 2的承载力及其 J acob i矩阵 ,首先把水导轴承修正的 轴承上端压力 轴承半径 55 mm 4 Kg / cm 2R eyno ld s方程转化为一个自由边值问题 ,其次针对 轴承下端压力 轴承间隙 2. 0 mm 6 Kg / cm 3 31480 r /m in 水密度 自转角速度 水导轴承两端压力不为零的情况 ,给出了压力分布 10Kg /m - 3 150 水黏度 轴瓦包角 0. 45 10 Pa. S 函数 ,建立一 维 有限 元互 补 方程 。计算 结果 表 明 , 本文方法能有效地求解湍流工况下 水 导轴 承承 载 力及其 Jacob i矩 阵 , 而 且 可 以 保 证 相 当 高 的 计 算 精度 。 参考文献 : [ 1 ] N g C. W. , Pan C. H. T. A linea rized tu rbu len t lub rica2 tion theo ry [ J ]. Tran s A SM E Ser D , 1965, 87: 675 - 688. [ 2 ] 张运清 , 张直明 . 一种紊流润滑理论分析新方法 2复合 ( ) 型紊流模式理论 [ J ]. 摩擦学学报 , 1995 , 15 3 : 271 - 275. [ 3 ] 张新敏 , 等 . 一种湍流 润滑理论分析的 工 程 计 算 方 法 [ J ]. 润滑与密封 , 2002, 2: 4 - 7. [ 4 ] 戴学余 , 苗旭升 ,等 . 几种低粘度润滑介质下动静压轴 承的性能分析 [ J ]. 润滑与密封 , 2004 , 3: 9 - 13. [ 5 ] 郑铁生 , 张 文 , 马建敏 ,等 . 滑动轴承非线性油膜力 的几个理论问题 及应用 [ J ]. 水 动 力 学 研 究 与 进 展 A ( ) 辑 , 2003, 18 6: 761 - 768. [ 6 ] 刘大全 , 肖忠会 , 张 文 ,等 . 滑动轴承非线性油膜力 ( ) 的一位直接解法 [ J ]. 机械工程学报 , 2005, 41 2 : 51 - 56.