基于多尺度算法的空蚀风险预报方法及装置

1.本发明属于工程计算流体力学领域，具体涉及基于多尺度算法的空蚀风险预报方法及装置。


背景技术：

2.空化是高速水动力学中的一个重要问题。在水中，空化往往会产生多尺度的空泡团，空泡团的生长和溃灭使得空化流动常常表现出非线性、强瞬态、高脉动等特点。同时，空化流动中气泡的溃灭是诱发空蚀的重要原因，使水泵、水轮机、螺旋桨的叶片发生磨损破坏，进一步影响流体机械的安全稳定运行。
3.目前，空蚀风险的评估大部分都是基于传统的欧拉空化流场进行计算的，这种方法的精确度非常依赖计算网格的分辨率，只能模拟大于网格尺度的空腔流动和其产生的空蚀破坏，无法评估流场中微尺度气泡溃灭对空蚀风险的贡献，给深入研究空蚀破坏的微观机理和提出有效的抗空蚀措施带来了巨大的挑战和困难。因此，有必要开发一种多尺度的空蚀风险预报方法，来深入研究流场中微观球形气泡运动对近壁面产生的影响，进一步揭示空蚀破坏的微观机理。


技术实现要素：

4.本发明是为了解决上述问题而进行的，目的在于提供基于多尺度算法的空蚀风险预报方法及装置，能够充分考虑空化流动中的宏观流场对微观球形气泡的影响，确定流场中球形气泡非对称坍塌对壁面的冲击和空蚀破坏风险。
5.本发明为了实现上述目的，采用了以下方案：
6.《方法》
7.本发明提供的基于多尺度算法的空蚀风险预报方法，其特征在于，包括以下步骤：
8.步骤1、采用欧拉方法计算宏观的空化流动；
9.将待预报区域作为计算区域，对计算区域划分网格单元，采用欧拉方法对每一个时刻网格单元进行计算求解，得到对应各时刻网格单元上的流场信息；
10.步骤2、基于每一个时刻网格单元上流场的信息，识别空化流动中的空腔并将小空腔转化为单个球形气泡；
11.步骤3、计算微观的球形气泡运动；
12.步骤4、根据多尺度流场信息利用球形气泡非对称坍塌模型计算球形气泡对壁面的冲击压力；
13.步骤5、根据球形气泡对壁面的冲击压力得到球形气泡对壁面的累计冲击能量进而确定流场的空蚀风险，实现空蚀风险预报。
14.优选地，本发明提供的基于多尺度算法的空蚀风险预报方法，在步骤1中，采用欧拉方法中的均质流模型计算宏观的空化流动，计算公式如下：
[0015][0016][0017]
式中，ρ为气液混合密度；uj是j方向上的速度；p是压力；μ是粘度；sb是源项。为了能够顺利求解，采用大涡模拟方法中的wall-adapting local eddy-viscosity(wale)亚格子模型和schnerr-sauer空化模型封闭计算方程组。同时在方程(2)中加入源项sb来体现微观球形气泡运动对宏观空化流场的影响，这种影响的计算公式如下：
[0018][0019]
式中，k表示一种高斯分布的调节因子；δub表示球形气泡速度的变化量；mb表示球形气泡的质量；vi表示第i个网格单元的体积；δt表示时间间隔。
[0020]
优选地，本发明提供的基于多尺度算法的空蚀风险预报方法，在步骤2中，任一时刻网格单元的空腔识别和转化过程具体包括如下子步骤：
[0021]
步骤2-1、对所有的网格单元依次进行编号，然后遍历每一个网格单元，读取各网格单元存储的流场含气率信息，若某个网格单元的含气率大于0，则将该网格单元的编号存储到列表中，得到存储有所有网格单元编号的网格编号表；
[0022]
步骤2-2、对网格编号表中的网格单元，按照编号顺序依次进行检查，采用空腔编号器为当前被检查的网格单元设置空腔编号n；
[0023]
步骤2-3、判断当前被检查的网格单元是否在相邻网格单元，若存在，则将所有相邻网格单元所隶属的空腔编号设置为当前被检查的网格单元的空腔编号，接着将当前被检查的网格单元和所有相邻网格单元的网格编号和空腔编号对应存储到空腔网格对应关系表中，然后进入步骤2-4；若不存在，则仅将当前被检测的网格单元的网格编号和空腔编号对应存储到空腔网格对应关系表中，返回步骤2-2检查下一个未被检查和判断过的网格单元，并令空腔编号器中n＝n+1；
[0024]
步骤2-4、进一步判断各相邻网格单元是否存也在未被检查和判断过的相邻网格单元，若存在，则将这些相邻网格单元所隶属的空腔编号设置为当前被检查的网格单元的空腔编号，并将这些相邻网格单元的网格编号和空腔编号对应存储到空腔网格对应关系表中，然后进入步骤2-5；若不存在，则认为这个空腔内包含的所有网格单元均已经检查完毕，返回步骤2-2检查下一个未被检查和判断过的网格单元，同时令空腔编号器中n＝n+1；
[0025]
步骤2-5、判断网格编号表中的所有网格单元均已经检查完毕，若是，进入步骤2-6，否则返回步骤2-2；
[0026]
步骤2-6、根据空腔网格对应关系表中空腔编号确定空腔数，根据各空腔编号所对应的网格编号确定各空腔所覆盖的网格单元数，进而计算得到各空腔的质心和体积；若某一空腔覆盖的网格单元数超出网格单元数第一阈值(一般可设置为8～30)，则在欧拉流场中删除该空腔，并在拉格朗日框架的原空腔质心位置处加入一个等体积的球形气泡；当拉格朗日框架中的球形气泡不断变大超过网格单元数第二阈值(略高于第一阈值，一般可设置为高于第一阈值5～10个网格单元数)时，将球形气泡删除，并将其影响以源项的形式加
入到欧拉框架下的气液两相输运方程中，重新计算新的大尺度空腔。
[0027]
优选地，本发明提供的基于多尺度算法的空蚀风险预报方法，在步骤3中，采用如下公式(1)计算微观流场中球形气泡的运动：
[0028][0029]
式中，mb表示气泡的质量；ub表示气泡运动速度；fa表示虚拟质量力；f
p
表示压力梯度力；f
buoy
表示浮力；f
drag
表示阻力；f
lift
表示升力，f
vol
表示体积变化力；
[0030]
采用显示-隐式耦合方法求解上述公式(4)：首先利用显式方法计算包含宏观欧拉场中液相影响的项，包括压力梯度力项、浮力项、部分虚拟质量力项，然后采用隐式方法预测并修正与气泡有关的项，包括升力项、阻力项和其余虚拟质量力项，得到气泡的运动轨迹；
[0031]
采用newton-paphson方法计算气泡的平衡半径，之后利用变时间步长的方法计算气泡的生长和溃灭过程，计算公式如下：
[0032][0033][0034]
式中，r表示球形气泡半径；表示球形气泡半径的变化速度；表示球形气泡半径变化的加速度；σ
st
表示球形气泡表面张力；pb表示球形气泡内的压力；p
ext
表示球形气泡中心处的环境压力；δ表示气体密度函数；c0表示纯水中的声速；p
gm
表示不可凝结气体的压力；k表示压缩常数；r
max
表示球形气泡的最大半径。
[0035]
优选地，本发明提供的基于多尺度算法的空蚀风险预报方法，在步骤4中，通过openfoam利用含相变的两相流动求解程序获得宏观欧拉空化流场，以此为基础设定欧拉场与拉格朗日场转化标准，然后利用多尺度程序计算空化流动，获得多尺度空化流场；采用如下公式计算球形气泡非对称坍塌产生的对壁面的冲击压力：
[0036]
p
imp,non-sphere
＝p
imp,sphere
(β1(γ)+β2(γ))
ꢀꢀꢀꢀꢀꢀ
(7)
[0037]
式中，p
imp,sphere
表示球形气泡对称坍塌产生的冲击压力，p0表示环境压力，r表示球形气泡中心到壁面的距离，f表示球形气泡的状态，可以表示为f’表示f的一阶导数；
[0038]
对于球形气泡非对称坍塌产生的一次和二次冲击压力，修正因子β的表达式为：
[0039]
[0040][0041]
式中，β1(γ)和β2(γ)表示球形气泡非对称一次坍塌和二次坍塌修正系数，γ表示球形气泡到壁面的距离与球形气泡最大半径的比。
[0042]
优选地，本发明提供的基于多尺度算法的空蚀风险预报方法，在步骤5中，统计一段时间内气泡对壁面的冲击压力，考虑球形气泡的冲击能量与球形气泡的冲击压力存在非线性关系，令球形气泡的冲击能量正比于球形气泡冲击压力的平方，利用计算得到的球形气泡对壁面的累计冲击能量进而确定流场的空蚀风险，实现空蚀风险预报。
[0043]
《装置》
[0044]
进一步，本发明还提供了能够自动实现上述《方法》的，其特征在于，包括：
[0045]
流场信息获取部，将待预报区域作为计算区域，对计算区域划分网格单元，采用欧拉方法对每一个时刻网格单元进行计算求解，得到对应各时刻网格单元上的流场信息；
[0046]
识别转化部，基于每一个时刻网格单元上流场的信息，识别空化流动中的空腔并将小空腔转化为单个球形气泡；
[0047]
气泡运动计算部，计算微观的球形气泡运动；
[0048]
冲击压力计算部，根据多尺度流场信息利用球形气泡非对称坍塌模型计算球形气泡对壁面的冲击压力；
[0049]
空蚀风险预报部，根据球形气泡对壁面的冲击压力得到球形气泡对壁面的累计冲击能量进而确定流场的空蚀风险，实现空蚀风险预报；以及
[0050]
控制部，与流场信息获取部、识别转化部、气泡运动计算部、冲击压力计算部、空蚀风险预报部均通信相连，控制它们的运行。
[0051]
优选地，本发明提供的基于多尺度算法的空蚀风险预报装置，还可以包括：输入显示部，与流场信息获取部、识别转化部、气泡运动计算部、冲击压力计算部、空蚀风险预报部、控制部均通信相连，根据用户输入的操作指令，显示相应的信息。
[0052]
优选地，本发明提供的基于多尺度算法的空蚀风险预报装置，还可以具有以下特征：输入显示部能够根据控制指令对流场信息获取部得到的流程信息进行显示，对识别转化部的识别和转化情况进行显示，对气泡运动计算部计算的气泡运动情况进行显示，对冲击压力计算部计算的冲击压力进行显示，将空蚀风险预报部确定的空蚀风险以数据图表方式进行显示，并能够将空蚀情况显示在待预报区域的空蚀风险强度预报图中，从而直观地反应空蚀情况。
[0053]
优选地，本发明提供的基于多尺度算法的空蚀风险预报方法，在流场信息获取部中，采用欧拉方法中的均质流模型计算宏观的空化流动，计算公式如下：
[0054][0055]
[0056]
式中，ρ为气液混合密度；uj是j方向上的速度；p是压力；μ是粘度；sb是源项。为了能够顺利求解，采用大涡模拟方法中的wall-adapting local eddy-viscosity(wale)亚格子模型和schnerr-sauer空化模型封闭计算方程组。同时在方程(2)中加入源项sb来体现微观球形气泡运动对宏观空化流场的影响，这种影响的计算公式如下：
[0057][0058]
式中，k表示高斯分布的调节因子；δub表示球形气泡速度的变化量；mb表示球形气泡的质量；vi表示第i个网格单元的体积；δt表示时间间隔。
[0059]
优选地，本发明提供的基于多尺度算法的空蚀风险预报装置，在识别转化部中，对任一时刻网格单元的空腔识别和转化处理具体包括：
[0060]
步骤2-1、对所有的网格单元依次进行编号，然后遍历每一个网格单元，读取各网格单元存储的流场含气率信息，若某个网格单元的含气率大于0，则将该网格单元的编号存储到列表中，得到存储有所有网格单元编号的网格编号表；
[0061]
步骤2-2、对网格编号表中的网格单元，按照编号顺序依次进行检查，采用空腔编号器为当前被检查的网格单元设置空腔编号n；
[0062]
步骤2-3、判断当前被检查的网格单元是否在相邻网格单元，若存在，则将所有相邻网格单元所隶属的空腔编号设置为当前被检查的网格单元的空腔编号，接着将当前被检查的网格单元和所有相邻网格单元的网格编号和空腔编号对应存储到空腔网格对应关系表中，然后进入步骤2-4；若不存在，则仅将当前被检测的网格单元的网格编号和空腔编号对应存储到空腔网格对应关系表中，返回步骤2-2检查下一个未被检查和判断过的网格单元，并令空腔编号器中n＝n+1；
[0063]
步骤2-4、进一步判断各相邻网格单元是否存也在未被检查和判断过的相邻网格单元，若存在，则将这些相邻网格单元所隶属的空腔编号设置为当前被检查的网格单元的空腔编号，并将这些相邻网格单元的网格编号和空腔编号对应存储到空腔网格对应关系表中，然后进入步骤2-5；若不存在，则认为这个空腔内包含的所有网格单元均已经检查完毕，返回步骤2-2检查下一个未被检查和判断过的网格单元，同时令空腔编号器中n＝n+1；
[0064]
步骤2-5、判断网格编号表中的所有网格单元均已经检查完毕，若是，进入步骤2-6，否则返回步骤2-2；
[0065]
步骤2-6、根据空腔网格对应关系表中空腔编号确定空腔数，根据各空腔编号所对应的网格编号确定各空腔所覆盖的网格单元数，进而计算得到各空腔的质心和体积；若某一空腔覆盖的网格单元数超出网格单元数第一阈值(一般可设置为8～30)，则在欧拉流场中删除该空腔，并在拉格朗日框架的原空腔质心位置处加入一个等体积的球形气泡；当拉格朗日框架中的球形气泡不断变大超过网格单元数第二阈值(略高于第一阈值，一般可设置为高于第一阈值5～10个网格单元数)时，将球形气泡删除，并将其影响以源项的形式加入到欧拉框架下的气液两相输运方程中，重新计算新的大尺度空腔。
[0066]
优选地，本发明提供的基于多尺度算法的空蚀风险预报装置，在气泡运动计算部中，采用如下公式(4)计算微观流场中球形气泡的运动：
[0067]
[0068]
式中，mb表示气泡的质量；ub表示气泡运动速度；fa表示虚拟质量力；f
p
表示压力梯度力；f
buoy
表示浮力；f
drag
表示阻力；f
lift
表示升力，f
vol
表示体积变化力；
[0069]
采用显示-隐式耦合方法求解上述公式(4)：首先利用显式方法计算包含宏观欧拉场中液相影响的项，包括压力梯度力项、浮力项、部分虚拟质量力项，然后采用隐式方法预测并修正与气泡有关的项，包括升力项、阻力项和其余虚拟质量力项，得到气泡的运动轨迹；
[0070]
采用newton-paphson方法计算气泡的平衡半径，之后利用变时间步长的方法计算气泡的生长和溃灭过程，计算公式如下：
[0071][0072][0073]
式中，r表示球形气泡半径；表示球形气泡半径的变化速度；表示球形气泡半径变化的加速度；σ
st
表示球形气泡表面张力；pb表示球形气泡内的压力；p
ext
表示球形气泡中心处的环境压力；δ表示气体密度函数；c0表示纯水中的声速；p
gm
表示不可凝结气体的压力；k表示压缩常数；r
max
表示球形气泡的最大半径。
[0074]
优选地，本发明提供的基于多尺度算法的空蚀风险预报装置，在冲击压力计算部中，通过openfoam利用含相变的两相流动求解程序获得宏观欧拉空化流场，以此为基础设定欧拉场与拉格朗日场转化标准，然后利用多尺度程序计算空化流动，获得多尺度空化流场；采用如下公式计算球形气泡非对称坍塌产生的对壁面的冲击压力：
[0075]
p
imp,non-sphere
＝p
imp,sphere
(β1(γ)+β2(γ))
ꢀꢀꢀꢀꢀꢀꢀ
(7)
[0076]
式中，p
imp,sphere
表示球形气泡对称坍塌产生的冲击压力，p0表示环境压力，r表示球形气泡中心到壁面的距离，f表示球形气泡的状态，可以表示为f’表示f的一阶导数；
[0077]
对于球形气泡非对称坍塌产生的一次和二次冲击压力，修正因子β的表达式为：
[0078][0079][0080]
式中，β1(γ)和β2(γ)表示球形气泡非对称一次坍塌和二次坍塌修正系数，γ表示球形气泡到壁面的距离与球形气泡最大半径的比。
[0081]
优选地，本发明提供的基于多尺度算法的空蚀风险预报装置，在空蚀风险预报部中，统计一段时间内气泡对壁面的冲击压力，考虑球形气泡的冲击能量与球形气泡的冲击压力存在非线性关系，令球形气泡的冲击能量正比于球形气泡冲击压力的平方，利用计算得到的球形气泡对壁面的累计冲击能量进而确定流场的空蚀风险，实现空蚀风险预报。
[0082]
发明的作用与效果
[0083]
本发明所提供的基于多尺度算法的空蚀风险预报方法及装置，采用多尺度的欧拉-拉格朗日耦合算法，通过球形气泡非对称坍塌模型来确定流场中球形气泡非对称坍塌对壁面的冲击和空蚀破坏风险，能够使计算更加符合实际流动情况，并准确得到空化流动中的跨尺度空泡团运动情况和微观气泡在壁面附近溃灭所产生的空蚀风险，为深入研究空蚀破坏的微观机理和提出有效的抗空蚀方法提供了可靠的技术支撑。
附图说明
[0084]
图1为本发明涉及的基于多尺度算法的空蚀风险预报方法的流程图；
[0085]
图2为本发明实施例中喷管空蚀预报的1/8个喷管空蚀几何模型示意图；
[0086]
图3为本发明实施例中喷管下底板空蚀风险强度预报图。
具体实施方式
[0087]
以下结合附图对本发明涉及的基于多尺度算法的空蚀风险预报方法及装置的具体实施方案进行详细地说明。
[0088]
《实施例》
[0089]
如图1所示，本实施所提供的基于多尺度算法的空蚀风险预报方法，包括如下步骤：
[0090]
步骤1、采用欧拉方法计算宏观的空化流动；
[0091]
如图2所示为本实施例中1/8个喷管空蚀预报算例的几何模型示意图，考虑到喷管是旋转结构，只选取其中的1/8部分喷管几何模型进行模拟计算。喷管垂直进口直径为16mm，水平段流道的高度为2.5mm，空化数σ＝0.9，上游来流流速v0＝6.25l/s。
[0092]
本实施例中，1/8个喷管为待预报空蚀情况的计算区域，对该计算区域划分网格单元，采用欧拉方法对每一个时刻网格单元进行计算求解，得到对应各时刻网格单元上的流场信息。
[0093]
采用欧拉方法中的均质流模型计算宏观的空化流动，计算公式如下：
[0094][0095][0096]
式中，ρ为气液混合密度；uj是j方向上的速度；p是压力；μ是粘度；sb是源项。为了能够顺利求解，采用大涡模拟方法中的wall-adapting local eddy-viscosity(wale)亚格子模型和schnerr-sauer空化模型封闭计算方程组。同时在方程(2)中加入源项sb来体现微观球形气泡运动对宏观空化流场的影响，这种影响的计算公式如下：
[0097][0098]
式中，k表示高斯分布的调节因子；δub表示球形气泡速度的变化量；mb表示球形气泡的质量；vi表示第i个网格单元的体积；δt表示时间间隔。
[0099]
步骤2、基于每一个时刻网格单元上流场的信息，识别空化流动中的空腔并将小空腔转化为单个球形气泡；任一时刻网格单元的空腔识别和转化过程具体包括如下子步骤：
[0100]
步骤2-1、对所有的网格单元以自然数的形式依次进行编号，然后从左到右、从上到下依次遍历每一个网格单元，读取每一个网格单元存储的流场含气率信息，若某个网格单元的含气率大于0，则将该网格单元的编号存储到列表中，得到存储有所有网格单元编号的网格编号表；
[0101]
步骤2-2、对网格编号表中的网格单元，按照编号从小到大顺序依次进行检查，采用空腔编号器为当前被检查的网格单元设置空腔编号n，初始值设为1；
[0102]
步骤2-3、判断当前被检查的网格单元是否在相邻网格单元，若存在，则将所有相邻网格单元所隶属的空腔编号设置为当前被检查的网格单元的空腔编号，接着将当前被检查的网格单元和所有相邻网格单元的网格编号和空腔编号对应存储到空腔网格对应关系表(该表包含两列，一列存储网格编号，另一列存储对应的空腔编号)中，然后进入步骤2-4；若不存在，则仅将当前被检查的网格单元的网格编号和空腔编号对应存储到空腔网格对应关系表中，返回步骤2-2按照网格单元的编号大小检查下一个未被检查和判断过(未被存储至空腔网格对应关系表中)的网格单元，并令空腔编号器中n＝n+1；
[0103]
步骤2-4、进一步判断各相邻网格单元是否存也在未被检查和判断过的相邻网格单元，若存在，则将这些相邻网格单元所隶属的空腔编号设置为当前被检查的网格单元的空腔编号，并将这些相邻网格单元的网格编号和空腔编号对应存储到空腔网格对应关系表中，然后进入步骤2-5；若不存在，则认为这个空腔内包含的所有网格单元均已经检查完毕，返回步骤2-2按照网格单元的编号大小检查下一个未被检查和判断过的网格单元，同时令空腔编号器中n＝n+1；
[0104]
步骤2-5、判断网格编号表中的所有网格单元均已经检查完毕，若是，表明流场中所有的空腔都识别完毕并编号，在空腔网格对应关系表中得到了每一个空腔的编号和其所包含的所有网格单元的编号，进入步骤2-6，否则返回步骤2-2；
[0105]
步骤2-6、根据空腔网格对应关系表中空腔编号确定空腔数，根据各空腔编号所对应的网格编号确定各空腔所覆盖的网格单元数，进而计算得到各空腔的质心和体积；若某一空腔覆盖的网格单元数超出使用者设置的网格单元数第一阈值(一般可设置为8～30)，则在欧拉流场中删除该空腔，并在拉格朗日框架的原空腔质心位置处加入一个等体积的球形气泡；当拉格朗日框架中的球形气泡不断变大超过使用者设置的网格单元数第二阈值(略高于第一阈值，一般可设置为高于第一阈值5～10个网格单元数)时，将球形气泡删除，并将其影响以源项的形式加入到欧拉框架下的气液两相输运方程中，重新计算新的大尺度空腔。
[0106]
步骤3、采用如下公式计算微观流场中球形气泡的运动：
[0107]
[0108]
式中，mb表示气泡的质量；ub表示气泡运动速度；fa表示虚拟质量力；f
p
表示压力梯度力；f
buoy
表示浮力；f
drag
表示阻力；f
lift
表示升力，f
vol
表示体积变化力；
[0109]
采用显示-隐式耦合方法求解上述公式(4)：首先利用显式方法计算包含宏观欧拉场中液相影响的项，包括压力梯度力项、浮力项、部分虚拟质量力项，然后采用隐式方法预测并修正与气泡有关的项，包括升力项、阻力项和其余虚拟质量力项，得到气泡的运动轨迹；
[0110]
采用newton-paphson方法计算气泡的平衡半径，之后利用变时间步长的方法计算气泡的生长和溃灭过程，计算公式如下：
[0111][0112][0113]
式中，r表示球形气泡半径；表示球形气泡半径的变化速度；表示球形气泡半径变化的加速度；σ
st
表示球形气泡表面张力；pb表示球形气泡内的压力；p
ext
表示球形气泡中心处的环境压力；δ表示气体密度函数，近似取为0.99882；c0表示纯水中的声速；p
gm
表示不可凝结气体的压力；k表示压缩常数，在气泡生长时取为1，在球形气泡溃灭时取为1.4；r
max
表示球形气泡的最大半径。
[0114]
步骤4、根据多尺度流场信息利用球形气泡非对称坍塌模型计算球形气泡对壁面的冲击压力，具体包括如下子步骤：
[0115]
步骤4-1、通过openfoam利用含相变的两相流动求解程序获得宏观欧拉空化流场；
[0116]
步骤4-2、在多尺度程序中设定欧拉场与拉格朗日场转化标准：当欧拉场中的宏观空腔所占据的网格单元数大于8时，空腔转化为拉格朗日场中的单个球形气泡；当拉格朗日场中的单个球形气泡所占据的网格单元数大于15时，单个球形气泡被替换为等体积的空腔，在欧拉场中进行计算；
[0117]
步骤4-3、在openfoam中以得到的宏观欧拉空化流场为基础利用多尺度程序计算空化流动，获得多尺度空化流场；
[0118]
步骤4-4、采用如下公式计算球形气泡非对称坍塌产生的对壁面的冲击压力：
[0119]
p
imp,non-sphere
＝p
imp,sphere
(β1(γ)+β2(γ))
ꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0120]
式中，p
imp,sphere
表示球形气泡对称坍塌产生的冲击压力，p0表示环境压力，r表示球形气泡中心到壁面的距离，f表示球形气泡的状态，可以表示为f’表示f的一阶导数；
[0121]
对于球形气泡非对称坍塌产生的一次和二次冲击压力，修正因子β的表达式为：
[0122][0123][0124]
式中，β1(γ)和β2(γ)表示球形气泡非对称一次坍塌和二次坍塌修正系数，γ表示球形气泡到壁面的距离与球形气泡最大半径的比。
[0125]
步骤5、根据球形气泡对壁面的冲击压力得到球形气泡对壁面的累计冲击能量进而确定流场的空蚀风险，实现空蚀风险预报。统计5～10个空化演变周期(对应的真实时间为0.015-0.03s)内气泡对壁面的冲击压力，考虑球形气泡的冲击能量与球形气泡的冲击压力存在非线性关系，令球形气泡的冲击能量正比于球形气泡冲击压力的平方，利用计算得到的球形气泡对壁面的累计冲击能量进而确定流场的空蚀风险，实现空蚀风险预报。
[0126]
图3为统计了一段时间后的喷管下底板的球形气泡冲击能量累计分布图，从图3中可以发现球形气泡的冲击能量大部分都分布在半径19～32mm内，将该预报情况与真实的喷管空蚀实验结果进行比较，发现吻合性非常好，证明了本发明所提出的基于多尺度算法的空蚀风险预报方法的有效性和可靠性。
[0127]
进一步，本实施例还提供了基于以上方法自动进行空蚀风险预报的装置，该装置包括流场信息获取部、识别转化部、气泡运动计算部、冲击压力计算部、空蚀风险预报部、输入显示部以及控制部。
[0128]
流场信息获取部执行上文步骤1描述的内容，将待预报区域作为计算区域，对计算区域划分网格单元，采用欧拉方法对每一个时刻网格单元进行计算求解，得到对应各时刻网格单元上的流场信息。
[0129]
识别转化部执行上文步骤2描述的内容，基于每一个时刻网格单元上流场的信息，识别空化流动中的空腔并将小空腔转化为单个球形气泡。
[0130]
气泡运动计算部执行上文步骤3描述的内容，计算微观的球形气泡运动。
[0131]
冲击压力计算部执行上文步骤4描述的内容，根据多尺度流场信息利用球形气泡非对称坍塌模型计算球形气泡对壁面的冲击压力。
[0132]
空蚀风险预报部执行上文步骤5描述的内容，根据球形气泡对壁面的冲击压力得到球形气泡对壁面的累计冲击能量进而确定流场的空蚀风险，实现空蚀风险预报。
[0133]
输入显示部根据用户输入的操作指令，显示相应的信息。例如，输入显示部能够根据控制指令对流场信息获取部得到的流程信息进行显示，对识别转化部的识别和转化情况进行显示，对气泡运动计算部计算的气泡运动情况进行显示，对冲击压力计算部计算的冲击压力进行显示，将空蚀风险预报部确定的空蚀风险以数据图表方式进行显示，并能够将空蚀情况显示在待预报区域的空蚀风险强度预报图中。
[0134]
控制部与流场信息获取部、识别转化部、气泡运动计算部、冲击压力计算部、空蚀风险预报部、输入显示部均通信相连，控制它们的运行。
[0135]
以上实施例仅仅是对本发明技术方案所做的举例说明。本发明所涉及的基于多尺
度算法的空蚀风险预报方法及装置并不仅仅限定于在以上实施例中所描述的内容，而是以权利要求所限定的范围为准。本发明所属领域技术人员在该实施例的基础上所做的任何修改或补充或等效替换，都在本发明的权利要求所要求保护的范围内。