技术领域
[0001] 本
发明涉及一种诱导轮低温流体空化流固耦合数值预测方法，属于
叶轮机械技术领域。
背景技术
[0002]
涡轮泵是液体火箭
发动机系统的重要组件之一，其功能是有涡轮带动
离心泵告诉运转以提高火箭推进剂的压
力，从而提高为火箭产生的动力。在高速涡轮
泵叶轮中，由于转速高致使叶轮前缘吸力面压力较低容易产生空化，当
空泡扩大到堵塞
转子叶片流道时，泵的吸入性能就会急剧下降。而诱导轮具有稠度大，安装
角小等优点，所以常在
涡轮泵上游安装诱导轮以避免涡轮泵发生空化。
[0003] 而作为涡轮泵内主要的空化发生部位，诱导轮内部空化会使流场内产生强烈
压力脉动，由此会引发诱导轮流机振动问题。同时，液体
火箭发动机通常以液氢，液
氧的低温液体作为
燃料，而低温介质诱导轮的空化特性由于普通常温
水诱导轮不同。研究诱导轮低温流激振动问题，主要就是研究诱导轮在低温空化、旋转
失速以及动静部件间流场干涉等非定常流动引发的结构振动，建立流动特征与涡轮泵结构振动特性之间的联系。诱导轮低温空化流激振动是影响涡轮泵振动问题的关键因素，如何抑制诱导轮低温空化流激振动是解决涡轮泵减振的关键问题之一。
[0004] 实验是研究诱导轮流激特性最直接的方法，能够较准确地捕捉到流场局部测点的压力脉动和结构的振动特性，但受实验条件的影响很难将结构的振动与流场脉动很好的联系起来。此外，实验对于设备和场地有很高的要求，因此，通过基于流固耦合的流激振动数值预测是解决此类问题的重要方式。目前国内在水力机械的流固耦合数值计算研究相对较少，而对诱导轮低温空化流激振动的流固耦合研究则更少，缺乏系统完善的数值预测方法，因此发展一种诱导轮低温流固耦合数值预测方法有着重要意义。
发明内容
[0005] 针对现有对诱导轮低温空化流激振动的数值预测方法存在考虑不够全面、不能满足现有航天涡轮泵要求的问题，本发明公开的一种诱导轮低温流体空化流固耦合数值预测方法要解决的技术问题是：实现对诱导轮低温空化流固耦合数值预测，并能够提高预测
精度。本发明通过对诱导轮低温空化流固耦合数值预测，不仅能够获得准确的诱导轮叶片
载荷特性，更重要的是能够
对流场作用下诱导轮的叶片
变形情况进行评估，辅助实际进行诱导轮的设计与优化，节省实验成本和时间。
[0006] 本发明的目的是通过下述技术方案实现的。
[0007] 本发明公开的一种诱导轮低温流体空化流固耦合数值预测方法，根据实际工作条件，建立诱导轮模型及进出口管道模型并划分网格；设定诱导轮边界条件并确定流固耦合交界面；建立计算
流体力学模型，并对
湍流模型进行旋转修正，对空化模型进行
热力学修正；进行诱导轮空化定常初值计算；定常结果作为初值进行诱导轮空化非定常流固耦合数值计算；对计算结果进行后处理，获得不同时刻诱导轮的载荷特性以及叶片变形速率
云图，辅助实际进行诱导轮的设计与优化，节省实验成本和时间。
[0008] 本发明公开的一种诱导轮低温流体空化流固耦合数值预测方法，包括如下步骤：
[0009] 步骤一：建立诱导轮模型并划分网格。
[0010] 应用叶轮机械建模
软件建立诱导轮
几何模型，将诱导轮导入网格划分软件中，以实现结构化网格划分。
[0011] 作为优选，步骤一所述的机械建模软件优选Blade gen，网格划分软件优选ANSYS Turbo Grid、ANSYS ICEM。
[0012] 步骤二：建立进出口管道模型并划分网格。
[0013] 应用
三维建模软件建立诱导轮进出口管道几何模型，并将几何模型导入网格划分软件进行网格划分。步骤一和步骤二所划分的网格共同组成计算流域。
[0014] 作为优选，步骤二所述的三维建模软件优选Pro E、Solid Works，网格划分软件优选 ANSYS ICEM。步骤一与步骤二划分网格共同组成流域。
[0015] 步骤三：设定诱导轮边界条件。
[0016] 在有限元结构求解器中，设置诱导轮结构域材料；在诱导轮
转轴的
轴承支撑位置施加圆柱约束以固定约束位置轴向及径向位移；将导水锥、诱导轮
轮毂和诱导轮叶片设置为流固耦合交界面。
[0017] 在
计算流体动力学求解器中，设置进口为opening并给定总压，出口为
质量流量出口，流域所有壁面为无滑移壁面，进口段-诱导轮与诱导轮-出口段交界面类型采用Rotor Stater界面。
[0018] 步骤四：建立计算流体力学模型。
[0019] 计算流体力学模型包括流场控制方程和结构场控制方程。
[0020] 流场控制方程包括质量方程(1)和动量方程(2)：
[0021]
[0022]
[0023] 式中，表示函数对相应的变量求偏导数，ρ为流体的
密度，t为时间，ui、uj代表流体的速度分量，xi、xj代表流体的位置分量，p为流场入口处压强，μl和μt分别为流体的
层流和紊流粘性系数。
[0024] 结构场控制方程为：
[0025]
[0026] 式中， 为
节点加速度项， 为节点速度项，q(t)为节点位移项，M为结构体的质量矩阵，C为结构体阻尼矩阵，K为结构体
刚度矩阵，Q(t)为结构受到的载荷项。
[0027] 步骤五：建立考虑旋转修正的湍流模型。
[0028] 高速涡轮泵中诱导轮和叶轮具有很高的转速，流场中具有强烈的旋转
曲率效应，而传统的湍流模型很难捕捉到旋转曲率效应。为了捕捉到旋转曲率效应，采用的湍流模型是基于标准k-ε模型上对其中湍
动能生成项进行旋转曲率修正。标准的k-ε湍流模型的控制方程形式为：
[0029]
[0030]
[0031] 式中，k、ε分别为湍动能和湍流耗散率，t为时间，μ为湍流黏度，μt为湍流粘性系数，xi、 xj代表流体的位置分量，uj代表流体的速度分量，Pt为湍动能生成项，设定模型常数分别 Cε1、Cε2、σε、σk。
[0032] 作为优选，模型常数分别设定为：Cε1＝1.44，Cε2＝1.92，σε＝1.3，σk＝1.0。
[0033] 对其中湍动能生成项Pt进行旋转-曲率修正，乘以修正系数fr，即
[0034] Pt→Pt·fr (6)
[0035]
[0036] 其中，r*和 均与应变率张量和旋转率张量有关，设定模型常数分别Cr1，Cr2和Cr3。
[0037] 作为优选，模型常数分别设定为：Cr1＝1，Cr2＝2，Cr3＝1。
[0038] 式(4)、(5)、(6)、(7)为建立考虑旋转修正的湍流模型。
[0039] 步骤六：建立考虑低温介质热力学效应的空化模型。
[0040] 采用修正的Kubota空化模型以考虑
温度变化对质量传输过程的影响。Kubota模型是基于 Rayleigh-Plesset提出的空泡生长方程建立的，由R-P方程得到空化过程中相间质量传输速率，但Kubota空化模型是基于等温情况下得到的，计算得到的空化区域由当地
汽化压强决定。空化过程中液体汽化吸热，导致空泡附近液体温度降低，使得泡内和泡外形成温度差ΔT，空化区域温度下降的热力学效应对空泡的产生和发展具有重要影响。因此在Kubota空化模型的
蒸发和
凝结项中添加热力学效应项以考虑温度变化对质量传输过程的影响。
[0041] 任意t时刻时的热流量q为：
[0042]
[0043] 空泡界面上热平衡表示为：
[0044]
[0045] 其中，K＝αρlcp为热传导率，ΔT为时间，RB为空泡半径，ρv为气相密度，α为热扩散率， Lgv为蒸发
潜热。
[0046] 综合式(8)和式(9)可得：
[0047]
[0048] 基于式(10)对Kubota模型中蒸发项和凝结项进行修正，同时考虑液体中固有的气体含量，得到考虑热力学效应的质量传输蒸发项和凝结项分别为：
[0049]
[0050]
[0051] 其中：p为当地压力，pv为当地饱和蒸气压，Cp为液相
比热，Tl为当地温度，T∞为来流温度，ρl为液相密度，αv,αf为气相和杂质气体体积分数，Cε为蒸发系数，Cv为凝结系数。
[0052] 作为优选，蒸发和凝结系数Cε，Cv分别取1，0.0002。
[0053] 湍动能对空化产生重要的影响，在上述空化模型中，采用公式(13)计算湍动能k对当地汽化压强pv的影响：
[0054] pturb＝0.39ρmk (13)
[0055] 汽化压强采用下式计算：
[0056] pv＝(pv(Tl)+pturb/2) (14)
[0057] 式中pv(Tl)表示当地
饱和蒸汽压，ρm为平均密度。
[0058] 式(11)、(12)、(14)为建立考虑低温介质热力学效应的空化模型。
[0059] 步骤七：进行诱导轮流场空化定常初值计算。
[0060] 不考虑诱导轮内流场随时间的变化及诱导轮的变形，通过计算流体动力学求解器ANSYS CFX进行诱导轮流场定常初值计算，得到诱导轮流场内各参数初始值。
[0061] 步骤八：进行诱导轮空化非定常流固耦合数值计算。
[0062] 步骤八实现方法如下：
[0063] 步骤8.1，以步骤七所得计算结果作为初值，通过计算流体动力学求解器ANSYS CFX进行诱导轮空化非定常流动计算，得到流固耦合交界面的压力场；
[0064] 步骤8.2，将流固耦合交界面的压力场传递到结构场，在有限元结构求解器中对结构场控制方程进行离散求解，得到诱导轮转子系统的结构
应力与变形结果；
[0065] 步骤8.3，根据流固耦合交界面的变形结果，对网格进行更新，步骤8.1至步骤8.3为非定常流固耦合计算的一个
迭代步；
[0066] 步骤8.4，重复步骤8.1至步骤8.3直到迭代步计算结果收敛；
[0067] 步骤8.5，步骤8.1至步骤8.4为诱导轮非定常流固耦合计算的一个计算步，重复步骤8.1 至步骤8.4直到达到预定计算步数，即完成诱导轮流固耦合数值计算。
[0068] 步骤九：对步骤八所得诱导轮流固耦合数值计算结果进行后处理，提取不同时刻诱导轮的载荷特性以及叶片变形速率云图，得到空化对诱导轮结构变形作用随时间的变化关系，即实现对诱导轮低温空化流固耦合数值预测。
[0069] 还包括步骤十：利用步骤一至步骤九所述的预测方法预测低温流体诱导轮空化流固耦合现象，应用于相应工程问题，不仅能够获得准确的诱导轮叶片载荷特性，更重要的是能够对流场作用下诱导轮的叶片变形情况进行评估，辅助实际进行诱导轮的设计与优化，节省实验成本和时间。
[0070] 所述相应工程问题包括航空发动机涡轮泵流激振动问题、诱导轮空化不
稳定性问题等。
[0071] 有益效果：
[0072] 1、本发明公开的一种诱导轮低温流体空化流固耦合数值预测方法，通过对蒸发项以及凝结项进行热力学修正，充分考虑低温流体的空化热力学效应，能够准确地对流场空化情况进行预测。
[0073] 2、本发明公开的一种诱导轮低温流体空化流固耦合数值预测方法，通过对湍流模型中的湍动能生成项进行旋转修正，能够对准确地对旋转曲率效应进行预测。
[0074] 3、本发明公开的一种诱导轮低温流体空化流固耦合数值预测方法，通过对流场以及结构场的耦合计算，能够准确的对诱导轮空化过程中的结构变形进行预测。相比实验研究的方法，节约大量的时间和资源。
[0075] 4、本发明公开的一种诱导轮低温流体空化流固耦合数值预测方法，应用于航空发动机涡轮泵流激振动、诱导轮空化不稳定性问题领域，不仅能够获得准确的诱导轮叶片载荷特性，更重要的是能够对流场作用下诱导轮的叶片变形情况进行评估，辅助实际进行诱导轮的设计与优化，节省实验成本和时间。
附图说明
[0076] 图1为本发明公开一种诱导轮低温流体空化流固耦合数值预测方法
流程图；
[0077] 图2为本发明诱导轮网格划分结果；
[0078] 图3为本发明
实施例后处理结果，其中：图3a为t＝t0时刻叶片变形情况，图3b为 t＝t0+1/4T时刻叶片变形情况，图3c为t＝t0+2/4T时刻叶片变形情况，图3d为t＝t0+3/4T时刻叶片变形情况。
[0079] 图4为本发明流固耦合
算法具体实现流程图。
具体实施方式
[0080] 为了更好的说明本发明的目的和优点，下面结合附图和实例对发明内容做进一步说明。
[0081] 诱导轮在空化条件下内部流畅具有显著的压力脉动特性。在实际环境中这种压力脉动会对诱导轮结构产生一个交变载荷，使得诱导轮结构变形和振动，因此对诱导轮空化流动进行流固耦合计算十分必要。本实例通过一种诱导轮低温流体空化流固耦合数值计算方法，对诱导轮空化流固耦合特性进行预测，从而得到诱导轮的结构变形结果。
[0082] 如图1所示，本实施例公开的一种诱导轮低温流体空化流固耦合数值计算方法，具体实现步骤如下：
[0083] 步骤一：建立诱导轮模型并划分网格。
[0084] 应用叶轮机械建模软件Bladegen建立诱导轮几何模型，将诱导轮导入网格划分软件 ANSYS Turbo Grid中，以实现结构化网格划分。
[0085] 步骤二：建立进出口管道模型并划分网格。
[0086] 应用三维建模软件Solidworks建立诱导轮进出口管道几何模型，并将几何模型导入网格划分软件ANSYS ICEM进行网格划分。步骤一和步骤二所划分的网格共同组成计算流域。
[0087] 步骤三：设定诱导轮边界条件。
[0088] 在有限元结构求解器中，设置诱导轮结构域材料；在诱导轮转轴的轴承支撑位置施加圆柱约束以固定约束位置轴向及径向位移；将导水锥、诱导轮轮毂和诱导轮叶片设置为流固耦合交界面。
[0089] 在计算流体动力学求解器中，设置进口为opening并给定总压，出口为质量流量出口，流域所有壁面为无滑移壁面，进口段-诱导轮与诱导轮-出口段交界面类型采用Rotor Stater界面。
[0090] 步骤四：建立计算流体力学模型
[0091] 计算流体力学模型包括流场控制方程和结构场控制方程。
[0092] 流场控制方程包括质量方程(1)和动量方程(2)：
[0093]
[0094]
[0095] 式中，表示函数对相应的变量求偏导数，ρ为流体的密度，t为时间，ui、uj代表流体的速度分量，xi、xj代表流体的位置分量，p为流场入口处压强，μl和μt分别为流体的层流和紊流粘性系数。
[0096] 结构场控制方程为：
[0097]
[0098] 式中， 为节点加速度项， 为节点速度项，q(t)为节点位移项，M为结构体的质量矩阵，C为结构体阻尼矩阵，K为结构体刚度矩阵，Q(t)为结构受到的载荷项。
[0099] 步骤五：建立考虑旋转修正的湍流模型
[0100] 高速涡轮泵中诱导轮和叶轮具有很高的转速，流场中具有强烈的旋转曲率效应，而传统的湍流模型很难捕捉到旋转曲率效应。因此本方法采用的湍流模型是基于标准k-ε模型上，对其中湍动能生成项进行旋转曲率修正。标准的k-ε湍流模型的控制方程形式为：
[0101]
[0102]
[0103] 式中，k、ε分别为湍动能和湍流耗散率，t为时间，μ为湍流黏度，μt为湍流粘性系数，xi、 xj代表流体的位置分量，uj代表流体的速度分量，Pt为湍动能生成项，模型常数分别为： Cε1＝1.44，Cε2＝1.92，σε＝1.3，σk＝1.0。
[0104] 对其中湍动能生成项Pt进行旋转-曲率修正，乘以修正系数fr，即
[0105] Pt→Pt·fr (6)
[0106]
[0107] 其中，r*和 均与应变率张量和旋转率张量有关，模型常数分别设定为：Cr1＝1，Cr2＝2， Cr3＝1。
[0108] 步骤六：建立考虑热力学效应的空化模型
[0109] 本实施例采用修正的Kubota模型以考虑温度变化对质量传输过程的影响。Kubota模型是基于Rayleigh-Plesset提出的空泡生长方程建立的，由R-P方程可以得到空化过程中相间质量传输速率，但Kubota空化模型是基于等温情况下得到的，计算得到的空化区域由当地汽化压强决定。空化过程中液体汽化吸热，导致空泡附近液体温度降低，使得泡内和泡外形成温度差ΔT，空化区域温度下降的热力学效应对空泡的产生和发展具有重要影响。因此本文在 Kubota空化模型的蒸发和凝结项中添加了热力学效应项以考虑温度变化对质量传输过程的影响。
[0110] 任意t时刻时的热流量q为：
[0111]
[0112] 空泡界面上热平衡可以表示为：
[0113]
[0114] 其中，K＝αρlcp为热传导率，ΔT为时间，RB为空泡半径，ρv为气相密度，α为热扩散率， Lgv为蒸发潜热。
[0115] 综合式(8)和式(9)可得：
[0116]
[0117] 基于式(10)对Kubota模型中蒸发项和凝结项进行修正，同时考虑液体中固有的气体含量，得到考虑热力学效应的质量传输蒸发项和凝结项分别为：
[0118]
[0119]
[0120] 其中：p为当地压力，pv为当地饱和蒸气压，Cp为液相比热，Tl为当地温度，T∞为来流温度，ρl为液相密度，αv,αf为气相和杂质气体体积分数，Cε为蒸发系数，Cv为凝结系数。蒸发和凝结系数Cε，Cv分别取1，0.0002。
[0121] 湍动能对空化产生重要的影响，在上述空化模型中，采用公式(13)计算湍动能k对当地汽化压强pv的影响：
[0122] pturb＝0.39ρmk (13)
[0123] 汽化压强采用下式计算：
[0124] pv＝(pv(Tl)+pturb/2) (14)
[0125] 式中pv(Tl)表示当地饱和
蒸汽压，ρm为平均密度。
[0126] 步骤七：进行诱导轮流场空化定常初值计算
[0127] 不考虑诱导轮内流场随时间的变化及诱导轮的变形，通过计算流体动力学求解器ANSYS CFX进行诱导轮流场定常初值计算，得到诱导轮流场内各参数初始值。
[0128] 步骤八：进行诱导轮空化非定常流固耦合计算
[0129] 以步骤六所得计算结果作为初始条件进行诱导轮流场非定常的计算
[0130] 步骤8.1，在计算流体动力学求解器ANSYS CFX中设置计算时间步长为1.667*10-5s，总计算时间为0.5001s，以步骤六所得计算结果作为初值，进行诱导轮空化非定常流动计算，得到流固耦合交界面的压力场；
[0131] 步骤8.2，通过Workbench平台中的“System Coupling”模
块将流固耦合交界面的压力场传递到结构场，在有限元结构求解器Transient Structrual中设置计算时间步长为1.667*10-5s，总计算时间为0.5001s，与CFX设置保持一致。对结构场控制方程进行离散求解，得到诱导轮转子系统的结构应力与变形结果；
[0132] 步骤8.3，将Transient Structrual求解得到的结构变形结果传递到CFX中，对网格进行更新，得到新的流固耦合交界面的位置。步骤8.1至步骤8.3为非定常流固耦合计算的一个迭代步；
[0133] 步骤8.4，重复步骤步骤8.1至步骤8.3直到迭代步计算结果收敛；
[0134] 步骤8.5，步骤8.1至步骤8.4为诱导轮非定常流固耦合计算的一个计算步，重复步骤8.1 至步骤8.4直到计算步总时间达到0.5001s，即完成了诱导轮流固耦合数值计算。
[0135] 步骤九：用ANSYS CFD POST对步骤七计算结果进行后处理，提取四个不同时刻诱导轮沿转轴方向的变形，如图3所示，观察到诱导轮最大变形区域呈周期性变化，在t＝t0时刻，诱导轮变形最大位置位于叶片1和叶片2的叶顶前缘；而在t＝t0+1/4T时刻，最大变形位置则位于叶片1和叶片4前缘；当t＝t0+2/4T时刻，最大变形位置则位于叶片1和叶片2前缘；当 t＝t0+3/4T时刻，此时最大变形位置位于叶片1和叶片4前缘，如此往复。
[0136] 还包括步骤十：利用步骤一至步骤九所述的预测方法预测低温流体诱导轮空化流固耦合现象，应用于相应工程问题，不仅能够获得准确的诱导轮叶片载荷特性，更重要的是能够对流场作用下诱导轮的叶片变形情况进行评估，辅助实际进行诱导轮的设计与优化，节省实验成本和时间。
[0137] 所述相应工程问题包括航空发动机涡轮泵流激振动问题、诱导轮空化不稳定性问题等。
[0138] 本实施例通过一种诱导轮低温流体空化流固耦合数值计算方法，对低温介质下的诱导轮空化流固耦合特性进行了数值预测，得到了诱导轮叶片的载荷分布以及不同位置的变形情况，数值计算结果能够有效辅助诱导轮的设计工作，从而防止诱导轮叶片因空化流激振动而发生损坏。
[0139] 以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何
修改、等同替换、改进等，均应包含在本发明的保护范围之内。
