(19)中华 人民共和国 国家知识产权局 (12)发明 专利申请 (10)申请公布号 (43)申请公布日 (21)申请 号 202111636350.3 (22)申请日 2021.12.2 9 (71)申请人 武汉大学 地址 430072 湖北省武汉市武昌区珞珈山 武汉大学 (72)发明人 季斌　王紫阳　程怀玉　 (74)专利代理 机构 武汉科皓知识产权代理事务 所(特殊普通 合伙) 42222 代理人 俞琳娟 (51)Int.Cl. G06F 30/28(2020.01) G06F 113/08(2020.01) G06F 119/14(2020.01) (54)发明名称 基于多尺度算法的空蚀风险预报方法及装 置 (57)摘要 本发明提供基于多尺度算法的空蚀风险预 报方法及装置， 能够充分考虑空化流动中的宏观 流场对微观球形气泡的影 响， 确定流场中球形气 泡非对称 坍塌对壁面的冲击和空蚀破坏风险。 方 法包括： 步骤1、 将待预报区域作为计算区域， 对 计算区域划分网格单元， 采用欧拉方法对每一个 时刻网格单元进行计算求解， 得到对应各时刻网 格单元上的流场信息； 步骤2、 识别空化流动中的 空腔并将小空腔转化为单个球形气泡； 步骤3、 计 算微观的球形气泡运动； 步骤4、 根据多尺度流场 信息利用球形气泡非对称坍塌模型计算球形气 泡对壁面的冲击压力； 步骤5、 根据球形气泡对壁 面的冲击压力得到球形气泡对壁面的累计冲击 能量进而确定流场的空蚀风险， 实现 空蚀风险预 报。 权利要求书4页 说明书11页 附图3页 CN 114330165 A 2022.04.12 CN 114330165 A 1.基于多尺度算法的空蚀风险预报方法， 其特 征在于， 包括以下步骤： 步骤1、 采用欧拉方法计算宏观的空化 流动； 将待预报区域作为计算区域， 对计算区域划分网格单元， 采用欧拉方法对每一个时刻 网格单元进行计算 求解， 得到对应各时刻网格单 元上的流场信息； 步骤2、 基于每一个时刻网格单元上流场的信 息， 识别空化流动中的空腔并将小空腔转 化为单个球形气泡； 步骤3、 计算 微观的球形气泡 运动； 步骤4、 根据多尺度流场信息利用球形气泡非对称坍塌模型计算球形气泡对壁面的冲 击压力； 步骤5、 根据球形气泡对壁面的冲击压力得到球形气泡对壁面的累计冲击能量进而确 定流场的空蚀风险， 实现空蚀风险预报。 2.根据权利要求1所述的基于多尺度算法的空蚀风险预报方法， 其特 征在于： 其中， 在步骤2中， 任一时刻网格单 元的空腔识别和转 化过程具体包括如下子步骤： 步骤2‑1、 对所有的网格单元依次进行编号， 然后遍历每一个网格单元， 读取各网格单 元存储的流场含气率信息， 若某个网格单元的含气率大于0， 则将该网格单元的编号存储到 列表中， 得到存 储有所有网格单 元编号的网格编号表； 步骤2‑2、 对网格编号表中的网格单元， 按照编号顺序依次进行检查， 采用空腔编号器 为当前被 检查的网格单 元设置空腔编号 N； 步骤2‑3、 判断当前被检查的网格单元是否在相邻 网格单元， 若存在， 则将所有相邻 网 格单元所隶属的空腔编号设置为当前被检查的网格单元的空腔编号， 接着将当前被检查的 网格单元和所有相邻网格单元的网格编号和空腔编号对应存储到空腔网格对应关系表中， 然后进入步骤2 ‑4； 若不存在， 则仅 将当前被检查的网格单元的网格编号和空腔编号对应存 储到空腔网格对应关系表中， 返回步骤2 ‑2检查下一个未被检查和判断过的网格单元， 并令 空腔编号器中N ＝N+1； 步骤2‑4、 进一步判断各相邻网格单元是否存也在未被检查和判断过的相邻网格单元， 若存在， 则将这些相 邻网格单元所隶属的空腔编号设置为当前被检查的网格单元的空腔编 号， 并将这些相邻网格单元 的网格编号和空腔编号对应存储到空腔网格对应关系表中， 然 后进入步骤2 ‑5； 若不存在， 则认为这个空腔内包含的所有网格单元均已经检查完毕， 返回 步骤2‑2检查下一个未被 检查和判断过的网格单 元， 同时令空腔编号器中N ＝N+1； 步骤2‑5、 判断网格编号表中的所有网格单元均已经检查完毕， 若是， 进入步骤2 ‑6， 否 则返回步骤2 ‑2； 步骤2‑6、 根据空腔网格对应关系表中空腔编号确定空腔数， 根据各空腔编号所对应的 网格编号确定各空腔所覆盖的网格单元数， 进而计算得到各空腔的质心和体积； 若某一空 腔覆盖的网格单元数超出网格单元数第一阈值， 则在欧拉流场中删除该空腔， 并在拉格朗 日框架的原空腔质心 位置处加入一个等体积的球形气泡； 当拉格朗日框架中的球形气泡不 断变大超过网格单元数第二阈值时， 将球形气泡删除， 并将其影响以源项的形式加入到欧 拉框架下的气液两相输运方程中， 重新计算 新的大尺度空腔。 3.根据权利要求1所述的基于多尺度算法的空蚀风险预报方法， 其特 征在于： 其中， 在步骤3中， 采用如下公式计算 微观流场中球形气泡的运动：权　利　要　求　书 1/4 页 2 CN 114330165 A 2式中， mB表示气泡的质量； UB表示气泡运动速度； Fa表示虚拟质量力； Fp表示压力梯度力； Fbuoy表示浮力； Fdrag表示阻力； Flift表示升力， Fvol表示体积变化力； 采用显示 ‑隐式耦合方法求解上述公式(1)： 首先利用显式方法计算包含宏观欧拉场中 液相影响的项， 包括压力梯度力项、 浮力项、 部分虚拟质量力项， 然后采用隐式方法预测并 修正与气泡有关的项， 包括升力项、 阻力项和其 余虚拟质量力项， 得到气泡的运动轨 迹； 采用Newton ‑Paphson方法计算气泡的平衡半径， 之后利用变时间步长的方法计算气泡 的生长和溃灭过程， 计算公式如下： 式中， R表示球形气泡半径； 表示球形气泡半径的变化速度； 表示球形气泡半径变化 的加速度； σst表示球形气泡表面张力； pB表示球形气泡内的压力； pext表示球形气泡中心处 的环境压力； δ表示气体密度函数； c0表示纯水中的声速； pgm表示不可凝结气体的压力； k表 示压缩常数； Rmax表示球形气泡的最大半径。 4.根据权利要求1所述的基于多尺度算法的空蚀风险预报方法， 其特 征在于： 其中， 在步骤4中， 通过Op enFOAM利用含相变的两相流动求解程序获得宏观欧拉空化流 场， 以此为基础设定欧拉场与拉格朗日场转化标准， 然后利用多尺度程序计算空化流动， 获 得多尺度空化 流场； 采用如下公式计算球形气泡非对称坍塌产生的对壁 面的冲击 压力： pimp,non‑sphere＝pimp,sphere( β1(γ)+β2(γ))    (4) 式中， pimp,sphere表示球形气泡对称坍 塌产生的冲 击压力， p0表示环境压力， r表示球形气泡中心到壁面的距离， f表示球形气泡的 状态， f’表示f的一阶导数； 对于球形气泡非对称坍塌产生的一次和二次冲击 压力， 修正因子β 的表达式为： 式中， β1(γ)和β2(γ)表示球形气 泡非对称 一次坍塌和二次坍塌修正系数， γ表示球形权　利　要　求　书 2/4 页 3 CN 114330165 A 3