闫清东，刘博深，魏巍（.北京理工大学机械与车辆学院，北京0008；2.北京理工大学车辆传动国家重点实验室，北京0008）



液力变矩器流体-固体耦合压力脉动分析

闫清东1，2，刘博深1，魏巍1，2
（1.北京理工大学机械与车辆学院，北京100081；2.北京理工大学车辆传动国家重点实验室，北京100081）

摘要：在冲焊型高功率密度液力变矩器的设计过程中，需要考虑在油液非定常流动下叶片所受压力载荷脉动，以及在载荷脉动激励下结构的振动响应。采用基于动网格的流体-固体耦合方法，沿叶片入口至出口方向设定监测点，分析对应位置压力载荷脉动与叶轮振动时域特性，对载荷脉动进行频率转换并对叶轮模态进行频域分析。分析表明:涡轮叶片所受压力载荷脉动幅值最大处位于叶片入口与外环连接处；压力载荷脉动与叶片振动的幅值沿叶片入口到出口逐渐减弱，且随着速比升高载荷脉动幅值与叶片振动响应明显减弱；涡轮脉动峰值频率在叶轮第2阶与第3阶模态之间，随速比升高，压力载荷脉动频域幅值明显减弱。

关键词：兵器科学与技术；液力变矩器；流体-固体耦合；压力脉动；结构振动

0 引言

液力变矩器是一种广泛应用于传动的液力元件，其叶轮内部流场压力脉动是引起叶轮工作不稳定的主要原因之一。变矩器内流动复杂，尤其是在起动工况时泵轮、涡轮较大的转速差引起的非定常流动现象较为明显。

王玉川等［1］采用基于RNG和滑移网格的方法对不同工况下离心泵进行瞬态流场分析，分析了叶片压力脉动主频与叶轮转频的关系，认为叶轮出口处附近随时间变化的漩涡是内部流动不均匀的主要原因，这使得压力脉动强度增大。田辉等［2］分析了单级单吸蜗壳式离心泵的叶轮蜗壳的非定常流动特性，认为蜗壳的几何不对称影响了动静干涉作用在泵内部的传播。王秀礼等［3］采用数值模拟对核主泵3种不同进口直径短叶片进行对比，分析研究了泵内流量和叶片叶型对压力脉动的影响。Marathe等［4］采用五孔探针试验的试验方法测量了液力变矩器导轮叶片压力分布，认为涡轮与导轮间转速差造成油液流动不稳定是导轮叶片压力波动的主要原因。Schlienger等［5］研究了轴流式涡轮中涡流-尾流与叶片非定常流动机制的内在联系，认为油液流动不稳定的主要原因是动静干涉造成的流道漩涡和二次流。Tsukamoto等［6］分析了离心泵内流动特性，确定了泵入出口压力脉动与流体非常流动的关系。

上述文献主要是对研究对象内流体的流动状态进行了分析，但目前尚未结合流体-固体耦合方法对压力脉动与叶轮振动的内在联系进行分析研究。

本文采用流体-固体耦合计算方法计算单位时间步内流场特性与结构特性，考虑了冲焊型液力变矩器叶轮振动变形对流场的影响。采用动态网格保证叶轮振动变形时结构与流场交界面网格数据传递。通过该方法能够分析研究叶轮在更接近实际工作情况下时域与频域载荷脉动激励与叶片振动响应的内在联系，为进一步分析油液非定常流动状态和疲劳失效提供理论依据［7］。

1 数值方法

1. 1 流体-固体耦合理论

流体-固体耦合分析是在每一个时间迭代步长中同时分别进行流场分析与结构分析，通过流体-固体耦合交界面传递流场压力与结构位移，采用动网格对流场网格进行更新，完成单位时间步内的流体-固体耦合交互计算。

流体连续性方程和Navier-Stokes方程矢量形式为

式中：U为单元流体速度；f为流体单元所受质量力；ρ为流体密度，在计算中假设油液为不可压缩，其密度为常量；p为作用在微元上的表面压力；μ为流体运动黏度；t为时间；▽为哈密顿算子；▽2为拉普拉斯算子。

其结构方程为

式中：M为结构质量矩阵；C为阻尼矩阵；K为结构刚度矩阵；L为结构变形位移；F为载荷矩阵；A为分布载荷作用的表面积；S为形函数；p为节点所受压力。

流体-固体耦合同时遵循最基本的守恒原则。在流体-固体耦合面处如果满足流体与固体包括应力、位移等变量相等或守恒，即可满足方程

式中：τf为流体节点上的应力；τs为结构节点上的应力；n为流体-固体耦合面法向矢量；df为流体-固体耦合面上流体的位移；ds为流体-固体耦合面上结构的位移。

流体-固体耦合中的数据传递是指将流体计算结果和固体结构计算结果通过流体-固体耦合面相互传递。流体与结构的求解网格不同，需要进行数据传递插值运算。数据传递包括单元匹配和插值计算。本文采用桶式搜索匹配算法与保形插值法（见图1）.利用桶式算法划分源网格到一个桶空间中，桶空间就是网格中的一组单元，并以此计算源网格与目标网格的匹配权重系数。

图1 保形插值法Fig. 1 Shape-preserving interpolation

如图1所示，在保形插值法中，有目标网格节点Ti与源网格节点Si，Si在源网格面进行插值。传递数据为

式中：n为源网格桶内的节点数；wi为匹配权重系数；fTi为目标网格节点接受到的变量值；fSi为源网格节点上变量值。此耦合求解方法通过分别求解流场与结构方程，可以同时计算流场分析与结构变形，提高计算效率。

1. 2 分析模型

流场计算中湍流模型采用SST应力运输模型［8］，扩散项离散采用2阶中央差分格式，速度和压力耦合采用SIMPLE算法。为模拟计算叶轮之间相互作用，采用全叶轮流道滑移网格方法进行流场分析，流场网格选择为四面体非结构网格，总网格数为183万个单元［9］。为保证流体与结构协同变形，流体区域采用动网格技术保证网格更新时计算收敛，选取弹簧光顺模型进行网格更新，保证网格运动与计算收敛。

为了充分计算叶轮间不同转速造成的相互影响，计算步长应保证叶轮之间充分的相对转动。算例中泵轮叶片数为41，叶片间隔角度为8. 78°；涡轮叶片数为37，叶片间隔角度为9. 73°；导轮叶片数为19，叶片间隔角度为18. 95°.以泵轮与涡轮之间相对转动为基准，间隔Δθ为单个计算步长的Δθ=2°，取相对旋转3圈为总计算时长。

式中：ts为计算步长；Δθ为间隔角度；ωP为泵轮转速；ωT为涡轮转速；ttot为总计算时长；n为计算步数。泵轮转速为1 200 r/ min，计算时取液力变矩器起动工况和高效工况，速比分别为i = 0和i = 0. 7.速比为i = 0时，计算步长为ts= 0. 3×10-3s，总计算时长ttot=0. 18 s；i =0. 7时，计算步长为ts=0. 9× 10-3s，总计算时长ttot=0. 54 s.

表1 叶轮材料参数Tab. 1 Material parameters of turbine

冲焊型液力变矩器泵轮与涡轮是钢板冲压而成的。导轮是用铸铝制造的，刚度较高，变形较小，在计算中视为刚体。表1给出了涡轮的材料参数。泵轮与涡轮的约束形式如图2所示。泵轮约束1为泵轮外环外缘与视为刚体的罩壳连接，泵轮约束2为泵轮内缘与视为刚体的泵轮毂焊接固连。涡轮约束为涡轮内缘与涡轮毂铆接，涡轮毂视为刚体［10 -11］。

图2 叶轮约束Fig. 2 Structural constraint

1. 3 流体-固体耦合计算中的压力脉动监测点

液力变矩器在实际工作中叶片所受载荷是压力面与吸力面的压力差：

式中：pp为叶片压力面监测点所受压力；ps为吸力面监测点所受压力；Δp为该监测点叶片载荷。

为了分析叶片压力载荷脉动，仿真计算中在叶片表面设置监测压力载荷值与振动幅值的监测点，根据监测位置的不同监测点（见图3）主要分为两部分。第一部分为叶片监测点，点1～点6为6组监测点，每一组两个监测点分别位于叶片压力面与吸力面，且两点共处于同一表面法向量。两点压力差为该点压力面沿其曲面法向作用的压力载荷，下文中简称点来代表监测点组。上述监测点可以得到计算总时长中沿叶片方向从入口到出口压力载荷的变化，以此分析其压力载荷波动及对应位置叶片的振动。第二部分为叶轮圆周监测点，叶片压力面上监测点A与C在点2两侧均布，监测点B与点2重合。此监测点可以得到在计算总时长中叶片入口处压力载荷沿叶轮圆周方向的变化，以此分析叶轮周向压力载荷波动［3］。

图3 叶片压力载荷监测点Fig. 3 Pairs of sampling points on blade

2 计算结果分析

2. 1 变矩器原始特性

变矩器试验台主要由冲焊型液力变矩器样机及试验包箱、电机、测功机、转速、转矩传感器和泵站组成。调节电机保持泵轮转速恒定，通过测功机对涡轮增加载荷，直至加载到最大负荷即速比为i = 0.每次加载稳定后，同时测得液力变矩器的泵轮和涡轮的转速和转矩。

以非耦合稳态流场分析结果作为计算初值，分别采用瞬态流场计算与流体-固体耦合计算两种方法计算变矩器原始特性。对比仿真计算与试验所得的原始特性如图4所示。图4中：Ke为试验变矩比，Kfsi为流体-固体固耦合计算变矩比，Kcfd为流场分析计算变矩比；λe为试验泵轮转矩系数，λfsi为流体-固体耦合计算泵轮转矩系数，λcfd为流场分析计算泵轮转矩系数；ηe为试验效率，ηfsi为流体-固体耦合计算效率，ηcfd为流场分析计算效率。流体-固体耦合计算原始特性与相较非耦合流程计算结果更接近试验数据，验证了仿真计算中油液载荷计算结果的正确性［12］。

图4 原始特性对比Fig. 4 Comparison of simulation and experimental results

2. 2 流体-固体耦合叶轮变形计算结果

为表征叶轮整体变形趋势，图5示出了不同工况速比下总计算时间内最大变形量的均方根值（RMS）.图5中，TRMS为涡轮最大变形量均方根值，PRMS为泵轮最大变形量均方根值。从图5可见，随速比升高即涡轮转速升高，涡轮的变形量明显下降。分析可知泵轮变形量较小，故本文主要以涡轮为分析对象。

图5 流体-固体耦合叶轮最大变形量Fig. 5 Deformation of pump and turbine

3 时域分析

3. 1 叶轮周向压力载荷脉动

对比启动工况与高效率工况下叶轮周向监测点压力载荷时域脉动幅值，启动工况时速比i =0，叶轮负荷较高；高效工况i =0. 7，样机效率达到最高。速比i =0. 8时，变矩比为1，液力变矩器进入偶合器工况。故选取启动工况（i =0）与高效工况（i =0. 7）作为分析对象工况，以研究非闭锁工况下压力脉动随涡轮转速即工况变化的趋势与特点，如图6所示。

图6 不同工况速比叶轮周向压力载荷脉动Fig. 6 Pressure load fluctuation along the tangential direction of turbine in different working conditions

在起动工况下i = 0时，3个监测点压力载荷脉动幅值基本相同，在0. 235～0. 245 MPa之间波动，且波动峰值相位基本相同。在高效工况下i = 0.7时，3个监测点压力载荷脉动幅值明显下降，在0.06～0. 066 MPa之间波动，波动曲线整体呈波动趋势，且相比于起动工况波动周期增加。

由以上分析可知，在固定工况速比下，各监测点压力载荷脉动幅值与峰值相位基本相同；随工况速比升高，泵轮与涡轮之间转速差缩小，涡轮转速受泵轮出口油液冲击减弱，压力载荷脉动幅值减小。周向各监测点载荷脉动幅值与载荷相位均差别较小，故可根据单叶片来表征全部叶片的载荷脉动与振动响应。

3. 2 叶片压力载荷脉动

分析叶片压力载荷时域特性［13 -14］。如图7所示，监测点1与点2处于涡轮叶片入口处载荷波动较高，点2为叶片与外环连接处所受冲击最大，载荷峰值为0. 567 MPa，最大波动幅值为0. 012 MPa；点3与点4位于叶片中段，受油液冲击较小；点5与点6位于涡轮叶片出口处，载荷与波动幅值均较小。

图7 工况i =0时的叶片压力载荷Fig.7 Time series of pulsing pressure load of blade for i =0

以载荷脉动为激励的叶片振动响应如图8所示。点1与点2处振动幅值较高，点2变形响应峰值可达0. 029 mm，最高振幅为1. 57×10-3mm.其他位置振动幅值沿叶片入口到出口逐渐衰减。出口处监测点6振动幅值已减弱到0. 29×10-3mm，下降约81. 5%.振动响应变化趋势与叶片载荷变化趋势相同，振动幅值与载荷脉动幅值沿叶片入口到出口逐渐衰减。

如图9所示，速比i =0. 7工况下压力载荷脉动明显减弱，入口处监控点2载荷峰值为0. 044 MPa，最大载荷脉动幅值为6. 98×10-3MPa.

如图10所示，速比i =0. 7工况下，各点变形响应随载荷脉动的减小而减小，点2变形响应峰值为9. 4×10-3mm，最高振幅为0. 15×10-3mm.

图11表示时域下工况变化对监测点2的载荷脉动与叶轮振动的影响，pmax为压力载荷峰值，Apmax为最大载荷脉动幅值，Smax为叶轮振动峰值，ASmax为叶轮振幅。

图8 工况i =0时的叶片变形响应Fig. 8 Vibration of blade for i =0

图9 工况i =0. 7时的叶片压力载荷Fig 9 Time series of pulsing pressure load of blade for i =0. 7

图10 工况i =0. 7时的叶片变形响应Fig. 10 Vibration of blade for i =0. 7

根据以上分析可知，由于油液的冲击，涡轮叶片所受压力载荷脉动幅值最高位置为叶片入口与外环连接处，压力载荷脉动沿叶片入口到出口逐渐减弱，同时以载荷作为激励叶片振动幅值随着随之减弱。随着工况变化，速比升高，泵轮与涡轮间转速差逐渐缩小，载荷脉动幅值与叶片变形响应明显减弱。

4 频域分析

为了分析压力载荷脉动与涡轮的振动关系，研究二者在频域内的激励与响应的关系，需要计算叶轮振动模态。表2为与流体-固体耦合计算采用相同结构设定计算的涡轮前5阶模态。

图11 工况变化对监测点2的载荷脉动与叶轮振动时域影响Fig. 11 Time series of pulsing pressure loads on sampling point 2 at different speed ratios

表2 涡轮前5阶段模态Tab. 2 The first five orders of turbine modal

对压力脉动载荷采用快速傅里叶变换方法进行频域转换，得到叶轮载荷脉动的频率特性，以及各点在不同工况时的压力载荷脉动频域特性［15 -17］，如图12和图13所示。

图12 工况i =0时各监测点压力载荷脉动频域Fig. 12 Frequency domain of pressure load fluctuation on blade for i =0

从图12中可以看出，起动工况i = 0下各个监测点频率和相位一致，约为820 Hz，最高脉动幅值频率820 Hz在叶轮第2阶与第3阶模态之间。其中监测点2载荷脉动频域幅值最高为9. 96 kPa，监测点1的幅值为3. 89 kPa，其余各点脉动频域幅值较低，频域幅值变化趋势与时域相同。从图13中可以看出，i =0. 7时的各个监测点载荷脉动峰值相位虽较接近，入口处点1与点2频率相位前移至740 Hz，靠近叶轮第2阶模态，但由于叶轮间转速差减小，监测点1与点2的脉动峰值明显减弱。由频域分析可知，油液载荷峰值振幅频域在叶轮第2阶与第3阶模态之间是叶轮振动的主要原因，随涡轮转速升高，涡轮载荷峰值频率降低且载荷幅值下降。

图13 工况i =0. 7时各监测点压力载荷脉动频域Fig. 13 Frequency domain of pressure load fluctuation on blade in speed ratio i =0. 7

5 结论

1）采用动网格技术实现冲焊型液力变矩器流体-固体耦合计算，通过流场动网格变形保证了计算收敛。对比流体-固体耦合计算原始特性与试验数据，验证了流体-固体耦合分析在流场计算结果的合理性。

2）在固定速比工况下液力变矩器周向各涡轮叶片压力载荷脉动幅值与峰值相位基本相同；单涡轮叶片所受载荷脉动与振动幅值沿入口到出口方向减弱，压力载荷脉动幅值最大位置位于叶片入口与外环连接处。

3）涡轮脉动峰值频率在叶轮第2阶与第3阶模态之间，随速比升高，泵轮与涡轮间相互作用减弱，压力载荷脉动频域幅值明显下降。

参考文献（References）

［1］ 王玉川，谭磊，曹树良，等.离心泵叶轮区瞬态流动及压力脉动特性［J］.机械工程学报，2014，50（10）：163 -169. WANG Yu-chuan，TAN Lei，CAO Shu-liang，et al. Characteristics of transient flow and pressure fluctuation in impeller for centrifugal pump［J］. Journal of Mechanical Engineering，2014，50（10）：163 -169.（in Chinese）

［2］ 田辉，郭涛，孙秀玲，等.离心泵内部动静干涉作用的数值模拟［J］.农业机械学报，2014，50（10）：92 -95. TIAN Hui，GUO Tao，SUN Xiu-ling，et al. Numerical simulation of unsteady flow in a centrifugal pump［J］. Transactions of the Chinese Society for Agricultural Machinery，2014，50（10）：92 - 95.（in Chinese）

［3］ 王秀礼，袁寿其，朱荣生，等.基于CFD数值模拟的复合叶轮核主泵压力脉动特性研究［J］.原子能科学技术，2014，48（1）：99 -105. WANG Xiu-li，YUAN Shou-qi，ZHU Rong-sheng，et al. Numerical simulation on pressure fluctuation of reactor coolant pump with complex impeller based on CFD technique［J］. Atomic Energy Science and Technology，2014，48（1）：99 -105.（in Chinese）

［4］ Marathe B V，Lakshminarayana B. Experimental investigation of steady and unsteady flow field downstream of an automotive torque converter turbine and inside the stator：part II—unsteady pressure on the stator blade surface［J］. Journal of Turbomachinery，1997，119（3）：634 -645.

［5］ Schlienger J，Kalfas A I，Abhari R S. Vortex-wake- blade interaction in a shrouded axial turbine［J］. Journal of Turbomachinery，2005，127（4）：699 -707.

［6］ Wang H，Tsukamoto H. Experimental and numerical study of unsteady flow in a diffuser pump at off-design conditions［J］. Journal of Fluids Engineering，2003，125（5）：767 -778.

［7］ 魏巍，闫清东，朱颜.液力变矩器叶片流固耦合强度分析［J］.兵工学报，2008，29（10）：1158 -1162. WEI Wei，YAN Qing-dong，ZHU Yan. Strength analysis of fluid solid interaction field of hydrodynamic torque converter vanes［J］. Acta Armamentarii，2008，29（10）：1158 -1162.（in Chinese）

［8］ Kamakoti R，Wei S. Fluid-structure interaction for aeroelastic applications［J］. Progress in Aerospace Sciences，2004，40（8）：535 -558.

［9］ Schmucker H，Flemming F，Coulson S. Two-way coupled fluid structure interaction simulation of a propeller turbine［C］∥25th IAHR Symposium on Hydraulic Machinery and systems. Timisoara，Romania：IAHR，2010.

［10］ Yan Q D，Wei W. Numeric simulation of single passage ternary turbulence model in hydraulic torque converter［J］. Journal of Beijing Institute of Technology，2003，12（2）：172 -175.

［11］ Saeed R A，Galybin A N，Popov V. 3D fluid-structure modelling and vibration analysis for fault diagnosis of Francis turbine using multiple ANN and multiple ANFIS［J］. IMES Discussion Paper，2013，34（34）：259 -27.

［12］ Pei J，Yuan S Q，Yuan J P，et al. The influence of the flow rate on periodic flow unsteadiness behaviors in a sewage centrifugal pump［J］. Journal of Hydrodynamics，2013，25（5）：702 -709.

［13］ Pei J，Dohmen H J，Yuan S Q，et al. Investigation of unsteady flow-induced impeller oscillations of a single-blade pump under off-design conditions［J］. Journal of Fluids & Structures，2012，35（35）：89 -104.

［14］ Liu C，Alexandrina U，Houston G. W，et al. Parametric analysis and optimization of inlet deflection angle in torque converters［J］. Journal of Fluids Engineering. 2015，137（3）：031101.

［15］ Flack R，Brun K. Fundamental analysis of the secondary flows and jet-wake in a torque converter pump—partⅠ：model and flow in a rotating passage［J］. Journal of Fluids Engineering. 2005，127（1）：66 -74.

［16］ Macphee D，Beyene A. Fluid-structure interaction of a morphing symmetrical wind turbine blade subjected to variable load†［J］. International Journal of Energy Research，2013，37（1）：69 -79.

［17］ Wilkosz B，Zimmermann M，Schwarz P，et al. Numerical investigation of the unsteady iInteraction within a close-coupled centrifugal compressor used in an aero engine［J］. China Pharmacy，2013，136（4）：179 -182.

Pressure Load Fluctuation Analysis of Torque Converter Based on Fluid-structure Interaction

YAN Qing-dong1，2，LIU Bo-shen1，WEI Wei1，2
（1. School of Mechanical Engineering，Beijing Institute of Technology，Beijing 100081，China；2. National Key Lab of Vehicular Transmission，Beijing Institute of Technology，Beijing 100081，China）

Abstract：The load fluctuation and vibration of blade should be analyzed in the design of high power density and compact stamping-welded torque converter. The simulation is based on fluid structure interaction and dynamic mesh. The load fluctuation and vibration in time domain are extracted by setting the monitoring points along the blade from inlet to outlet. Fast Fourier transform and modal analysis are used in frequency domain analysis. The analysis results indicate that the max load fluctuation and blade vibration are at the inlet of blade because of the hydraulic impact，and the amplitudes of load fluctuation and blade vibration decay along inlet to outlet at a certain speed ratio. As the speed of turbine is higher，the amplitudes also decay. The peak frequency of load fluctuation is between the second-order and third-order modals，and also the amplitude of load fluctuation decreases at higher speed ratio.

Key words：ordnance science and technology；hydrodynamic torque converter；fluid-structure interaction；pressure fluctuation；structural vibration

中图分类号：TH137. 332

文献标志码：A

文章编号：1000-1093（2016）04-0577-07

DOI：10. 3969/ j. issn. 1000-1093. 2016. 04. 001

收稿日期：2015-08-05

基金项目：国家自然科学基金青年科学基金项目（51475041）；国家部委预先研究项目（40402050202）；国家部委基础产品创新科研基金项目（VTDP-2104）

作者简介：闫清东（1964—），男，教授，博士生导师。E-mail：yanqd@ bit. email. cn