CN115563906A - 一种基于非定常欧拉两相流的多步长结冰计算方法和系统 - Google Patents

一种基于非定常欧拉两相流的多步长结冰计算方法和系统 Download PDF

Info

Publication number
CN115563906A
CN115563906A CN202211401890.8A CN202211401890A CN115563906A CN 115563906 A CN115563906 A CN 115563906A CN 202211401890 A CN202211401890 A CN 202211401890A CN 115563906 A CN115563906 A CN 115563906A
Authority
CN
China
Prior art keywords
icing
calculation
grid
field
water
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202211401890.8A
Other languages
English (en)
Other versions
CN115563906B (zh
Inventor
朱春玲
郭威
沈柯
邓思强
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University of Aeronautics and Astronautics
Original Assignee
Nanjing University of Aeronautics and Astronautics
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University of Aeronautics and Astronautics filed Critical Nanjing University of Aeronautics and Astronautics
Priority to CN202211401890.8A priority Critical patent/CN115563906B/zh
Publication of CN115563906A publication Critical patent/CN115563906A/zh
Application granted granted Critical
Publication of CN115563906B publication Critical patent/CN115563906B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Computing Systems (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Fluid Mechanics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明提供一种基于非定常欧拉两相流的多步长结冰计算方法和系统，其中方法包括获取空气流场初始条件和水滴场初始条件；构建基于任意物体几何模型生成的计算网格；将总结冰时间均匀划分为多个结冰时间步长；计算每个结冰时间步长内空气流场的非定常数值；计算每个结冰时间步长内水滴场的非定常数值；计算当前结冰时间步长的结冰外形；采用动网格技术重构计算网格；将计算网格重构前的空气流场和水滴场瞬态结果插值赋值给重构后的计算网格；总结冰时间结束，得到最终结冰计算结果。本发明符合结冰后的非定常效应，改进水滴场控制方程组且使用有保正性质的HLLC格式，改进了结冰计算流程无需假设初始水膜高度和结冰厚度，求解更稳定更准确。

Description

一种基于非定常欧拉两相流的多步长结冰计算方法和系统
技术领域
本发明属于飞机结冰数值模拟技术领域，尤其涉及一种基于非定常欧拉两相流的多步长结冰计算方法和系统。
背景技术
飞机结冰问题与飞行安全关联密切，一直是航空航天领域关注的重点。研究结冰问题的主要方法包括冰风洞试验、自然飞行试验和数值模拟等。数值模拟方法投入成本低，且可以研究更为丰富的结冰工况。
上世纪四十到五十年代，该时期研究的焦点是水滴运动轨迹的求解方法和结冰热力学模型。计算过冷水滴撞击特性的方法主要有 Lagrange 法和 Euler 法，前者是根据牛顿第二定律建立单个水滴的运动方程，在获得流场信息的基础上，对水滴运动方程进行求解，后者是将空气和水滴看作均相的两相流，在 Euler 坐标系下建立水滴运动微分方程，空气流场和水滴撞击特性经耦合计算一并获得。Taylor和 Langmuir先后根据牛顿第二定律推导出了水滴运动微分方程。Langmuir 提出了利用流场信息计算水滴运动轨迹的方法，由于当时获取流场信息的能力有限， 研究对象仅限于外形简单的物体。Bergrun发展了计算翼型前缘水滴撞击特性的方法，并以NACA0012 翼型为对象，进行了前缘水滴收集系数的计算；Brun 和 Mergler利用绕圆柱表面不可压流动流场信息，进行了水滴撞击特性计算。对于结冰热力学模型的研究，Messinger提出的模型反映出了结冰表面的主要传热项，并给出了相应的数值解法，为结冰的数值计算奠定了基础，他的研究得到了最为广泛的认可和应用。上世纪八十年代，该阶段结冰数值计算的研究工作主要有：基于流场数值解的水滴撞击特性计算、二维翼型前缘结冰冰形的计算、结冰数学模型的改进等。在该阶段，还出现了多尺寸分布水滴撞击特性的计算研究。Durst发展了 Euler法计算水滴撞击特性，并进行了对比验证。在 NASA的结冰研究风洞中，Olsen对翼型前缘结冰进行了大量的试验观测和研究， 给出了结冰数学模型的改进建议。MacArthur发展了模拟结冰的计算流程，将结冰的数值计算分成了空气流场计算、水滴撞击特性计算和结冰冰形计算三个主要步骤，结冰外形改变之后，重复以上步骤，直至算满给定的结冰时间。至今，几乎所有的结冰计算软件或代码仍然采用以上基本步骤。上世纪九十年代之后至今，飞机结冰问题的数值研究有了新的进展，体现为二维结冰数值模拟的深入研究、结冰模拟软件的兴起和验证、三维结冰数值模拟的出现、过冷大水滴（SLD）结冰研究的开展等。Myers等对Messinger模型做出重大改进，划分不同结冰类型，考虑水膜流动带来的非定常效应，提出了三维曲面上综合考虑水膜流动与冰层传热的Myers模型。目前已经开发出的结冰数值模拟代码包括：美国 NASA 结冰研究中心开发的 LEWICE及 LEWICE3D 软件；加拿大NTI公司研发的 FENSAP-ICE结冰软件，法国的 ONERA Icing Code结冰代码等等。
然而，传统的结冰数值模拟方法，大多数都无法实现结冰多步长自动化，需要手动再将结冰外形重新输入计算。且空气流场和水滴场大多都是定常计算，两相流结果都是稳态收敛解，无法体现结冰后复杂外形的空气流场和水滴场的非定常效应，这与实际结冰过程的非定常机理不符合。结冰方法大多基于Messinger模型或Myers模型，Messinger模型是准定常方法，与真实非定常结冰过程有差别，原始Myers模型是非定常求解过程，但有多个局限：（1）采用结冰厚度的作为结冰类型判据，网格重构后需要记录基于原曲面坐标系的结冰厚度，涉及复杂几何计算；（2）需要根据经验假设初始水膜和结冰厚度，这一点与实际物理过程不符，且假设的初值可能会引入误差。
发明内容
本发明针对现有技术中的不足，提供一种基于非定常欧拉两相流的多步长结冰计算方法和系统。
第一方面，本发明提供一种基于非定常欧拉两相流的多步长结冰计算方法，包括：
S1，获取空气流场初始条件和水滴场初始条件；所述空气流场初始条件包括来流压力、来流温度、来流速度和攻角；所述水滴场初始条件包括水滴平均容积直径和液态水含量；
S2，构建基于任意物体几何模型生成的计算网格；
S3，将总结冰时间均匀划分为多个结冰时间步长；
S4，根据空气流场初始条件和计算网格，计算每个结冰时间步长内空气流场的非定常数值；
S5，根据水滴场初始条件和计算网格，计算每个结冰时间步长内水滴场的非定常数值；
S6，根据空气流场和水滴场的非定常数值，计算当前结冰时间步长的结冰外形；
S7，根据当前结冰时间步长的结冰外形，采用动网格技术改动计算网格节点以重构计算网格，将重构后的计算网格用于下一个结冰时间步长的计算过程；
S8，将重构前的计算网格的空气流场和水滴场瞬态结果插值赋值给重构后的计算网格；
S9，重复步骤S2-S8，直至总结冰时间结束，得到最终结冰计算结果。
进一步地，所述根据空气流场初始条件和计算网格，计算每个结冰时间步长内空气流场的非定常数值，包括：
构建空气流场控制方程；所述空气流场控制方程为理想状态气体可压缩雷诺平均Navier-Stokes方程；空气流场控制方程空间离散采用HLLC格式，时间离散使用隐式LU-SGS双时间步迭代，外迭代为真实时间步Δt，内迭代伪时间步长利用CFL条件求解。
进一步地，所述根据水滴场初始条件和计算网格，计算每个结冰时间步长内水滴场的非定常数值，包括：
构建水滴场控制方程组的每个计算网格单元积分表达式：
Figure 581451DEST_PATH_IMAGE001
；
其中，t为计算时间；Ω为计算网格单元的体积；W d 为水滴场计算网格单元的体守恒通量；F d 为水滴场计算网格单元的面对流通量；S为计算网格单元其中一个面的面积；P和Q均为水滴场计算网格单元的源项；
Figure 892347DEST_PATH_IMAGE002
；
Figure 553135DEST_PATH_IMAGE003
；
Figure 101928DEST_PATH_IMAGE004
；
Figure 455549DEST_PATH_IMAGE005
；
其中，x、y和z分别为惯性坐标系的三个分量方向；ρ d 为水密度；α为水滴体积分数；ρ d α为表观水滴密度；u d 、v d 和w d 分别为水滴在x、y和z方向上的速度；V d 为水滴的逆变速度，V d =
Figure 362063DEST_PATH_IMAGE006
·
Figure 560963DEST_PATH_IMAGE007
；
Figure 229842DEST_PATH_IMAGE006
=（u d ，v d ，w d ）；计算单元的面单位法向量
Figure 754364DEST_PATH_IMAGE007
=（n x ，n y ，n z ）；n x ，n y 和n z 分别为面单位 法向量在x、y和z方向上的分量；g为重力加速度；d为水滴直径；K为空气和水滴作用因子；
Figure 508694DEST_PATH_IMAGE008
为质量守恒方程中的数值耗散项；θ为攻角；
根据以下公式计算水滴场计算网格单元的面对流通量：
Figure 511285DEST_PATH_IMAGE009
；
Figure 175615DEST_PATH_IMAGE010
；
Figure 605460DEST_PATH_IMAGE011
；
Figure 847085DEST_PATH_IMAGE012
；
Figure 653367DEST_PATH_IMAGE013
；
其中，下标为L的为面左侧单元的变量，下标为R的为面右侧单元的变量；V L 为面左 侧单元逆变速度；V L =
Figure 296838DEST_PATH_IMAGE014
·
Figure 897584DEST_PATH_IMAGE007
=
Figure 501871DEST_PATH_IMAGE015
；V R 为面右侧单元逆变速度；V R =
Figure 111844DEST_PATH_IMAGE016
·
Figure 875401DEST_PATH_IMAGE007
=
Figure 381469DEST_PATH_IMAGE017
；S M 为中间波速；S L 为左波速；S R 为右波速；
S M =
Figure 863266DEST_PATH_IMAGE018
；
S L =V L -q L
Figure 886716DEST_PATH_IMAGE019
；
S R =V R +q R
Figure 504780DEST_PATH_IMAGE019
；
Figure 181749DEST_PATH_IMAGE020
；
Figure 885262DEST_PATH_IMAGE021
；
Figure 837038DEST_PATH_IMAGE022
；
其中，q L 为面左侧单元密度比；q R 为面右侧单元密度比。
进一步地，所述根据空气流场和水滴场的非定常数值，计算当前结冰时间步长的结冰外形，包括：
获取结冰类型；
根据结冰类型采用不同的质能平衡方程计算结冰增长率；
将结冰增长率作为源项带入至水膜流动方程组进行计算；
循环计算结冰增长率并将结冰增长率作为源项带入至水膜流动方程组进行计算，直至结冰迭代时间步长结束，得到当前结冰时间步长的结冰外形。
第二方面，本发明提供一种基于非定常欧拉两相流的多步长结冰计算系统，包括：
获取模块，用于获取空气流场初始条件和水滴场初始条件；所述空气流场初始条件包括来流压力、来流温度、来流速度和攻角；所述水滴场初始条件包括水滴平均容积直径和液态水含量；
构建模块，用于构建基于任意物体几何模型生成的计算网格；
划分模块，用于将总结冰时间均匀划分为多个结冰时间步长；
第一计算模块，用于根据空气流场初始条件和计算网格，计算每个结冰时间步长内空气流场的非定常数值；
第二计算模块，用于根据水滴场初始条件和计算网格，计算每个结冰时间步长内水滴场的非定常数值；
第三计算模块，用于根据空气流场和水滴场的非定常数值，计算当前结冰时间步长的结冰外形；
网格重构模块，用于根据当前结冰时间步长的结冰外形，采用动网格技术改动计算网格节点以重构计算网格，将重构后的计算网格用于下一个结冰时间步长的计算过程；
插值赋值模块，用于将重构前的计算网格的空气流场和水滴场瞬态结果插值赋值给重构后的计算网格；
循环模块，用于重复执行构建模块至插值赋值模块的操作，直至总结冰时间结束，得到最终结冰计算结果。
进一步地，所述第一计算模块包括：
第一构建单元，用于构建空气流场控制方程；所述空气流场控制方程为理想状态气体可压缩雷诺平均Navier-Stokes方程；空气流场控制方程空间离散采用HLLC格式，时间离散使用隐式LU-SGS双时间步迭代，外迭代为真实时间步Δt，内迭代伪时间步长利用CFL条件求解。
进一步地，所述第二计算模块包括：
第二构建单元，用于构建水滴场控制方程组的每个计算网格单元积分表达式：
Figure 44028DEST_PATH_IMAGE001
；
其中，t为计算时间；Ω为计算网格单元的体积；W d 为水滴场计算网格单元的体守恒通量；F d 为水滴场计算网格单元的面对流通量；S为计算网格单元其中一个面的面积；P和Q均为水滴场计算网格单元的源项；
Figure 796958DEST_PATH_IMAGE002
；
Figure 987768DEST_PATH_IMAGE003
；
Figure 212076DEST_PATH_IMAGE004
；
Figure 804732DEST_PATH_IMAGE005
；
其中，x、y和z分别为惯性坐标系的三个分量方向；ρ d 为水密度；α为水滴体积分数；ρ d α为表观水滴密度；u d 、v d 和w d 分别为水滴在x、y和z方向上的速度；V d 为水滴的逆变速度，V d =
Figure 354662DEST_PATH_IMAGE006
·
Figure 501609DEST_PATH_IMAGE007
；
Figure 670553DEST_PATH_IMAGE023
=（u d ，v d ，w d ）；计算单元的面单位法向量
Figure 117715DEST_PATH_IMAGE007
=（n x ，n y ，n z ）；n x ，n y 和n z 分别为面单位 法向量在x、y和z方向上的分量；g为重力加速度；d为水滴直径；K为空气和水滴作用因子；
Figure 838547DEST_PATH_IMAGE008
为质量守恒方程中的数值耗散项；θ为攻角；
第一计算单元，用于根据以下公式计算水滴场计算网格单元的面对流通量：
Figure 738370DEST_PATH_IMAGE009
；
Figure 570059DEST_PATH_IMAGE010
；
Figure 747094DEST_PATH_IMAGE011
；
Figure 638826DEST_PATH_IMAGE012
；
Figure 760366DEST_PATH_IMAGE013
；
其中，下标为L的为面左侧单元的变量，下标为R的为面右侧单元的变量；V L 为面左 侧单元逆变速度；V L =
Figure 395747DEST_PATH_IMAGE014
·
Figure 817501DEST_PATH_IMAGE007
=
Figure 614556DEST_PATH_IMAGE015
；V R 为面右侧单元逆变速度；V R =
Figure 364337DEST_PATH_IMAGE016
·
Figure 537829DEST_PATH_IMAGE007
=
Figure 548511DEST_PATH_IMAGE017
；S M 为中间波速；S L 为左波速；S R 为右波速；
S M =
Figure 782046DEST_PATH_IMAGE018
；
S L =V L -q L
Figure 143757DEST_PATH_IMAGE019
；
S R =V R +q R
Figure 855361DEST_PATH_IMAGE019
；
Figure 625609DEST_PATH_IMAGE020
；
Figure 30045DEST_PATH_IMAGE021
；
Figure 613473DEST_PATH_IMAGE022
；
其中，q L 为面左侧单元密度比；q R 为面右侧单元密度比。
进一步地，所述第三计算模块包括：
获取单元，用于获取结冰类型；
第二计算单元，用于根据结冰类型采用不同的质能平衡方程计算结冰增长率；
第三计算单元，将结冰增长率作为源项带入至水膜流动方程组进行计算；
第四计算单元，用于循环计算结冰增长率并将结冰增长率作为源项带入至水膜流动方程组进行计算，直至结冰迭代时间步长结束，得到当前结冰时间步长的结冰外形。
第三方面，本发明提供一种计算机设备，包括处理器和存储器；其中，处理器执行存储器中保存的计算机程序时实现第一方面所述的基于非定常欧拉两相流的多步长结冰计算方法的步骤。
第四方面，本发明提供一种计算机可读存储介质，用于存储计算机程序；计算机程序被处理器执行时实现第一方面所述的基于非定常欧拉两相流的多步长结冰计算方法的步骤。
本发明提供一种基于非定常欧拉两相流的多步长结冰计算方法和系统，其中方法包括S1，获取空气流场初始条件和水滴场初始条件；所述空气流场初始条件包括来流压力、来流温度、来流速度和攻角；所述水滴场初始条件包括水滴平均容积直径和液态水含量；S2，构建基于任意物体几何模型生成的计算网格；S3，将总结冰时间均匀划分为多个结冰时间步长；S4，根据空气流场初始条件和计算网格，计算每个结冰时间步长内空气流场的非定常数值；S5，根据水滴场初始条件和计算网格，计算每个结冰时间步长内水滴场的非定常数值；S6，根据空气流场和水滴场的非定常数值，计算当前结冰时间步长的结冰外形；S7，根据当前结冰时间步长的结冰外形，采用动网格技术改动计算网格节点以重构计算网格，将重构后的计算网格用于下一个结冰时间步长的计算过程；S8，将重构前的计算网格的空气流场和水滴场瞬态结果插值赋值给重构后的计算网格；S9，重复步骤S2-S8，直至总结冰时间结束，得到最终结冰计算结果。
本发明的方法中欧拉两相流（空气流场和水滴场）的求解结果不再是稳态收敛结果，而是将其耦合求解过程变成了非定常过程，得出空气流场和水滴场的非定常积分结果，用于后续的结冰计算，更符合结冰后复杂外形附近空气流场和水滴场出现大分离现象或涡结构带来的非定常效应，水滴场求解时改进了水滴控制方程组且使用了适用于水滴场的具有保正性质的HLLC空间离散格式，使得结冰多步长计算时方程组求解更稳定不易发散。
本发明的方法中结冰计算方法基于Myers模型改进，将当前表面温度替换结冰厚度作为结冰类型的判据，避免了网格重构后较为复杂的几何运算，同时改进结冰计算流程，无需再根据经验假设初始水膜高度和结冰厚度，避免了假设初始值导致的误差。
本发明的方法将重构前计算网格的两相流结果插值赋值给重构后的新计算网格，无需初始化可直接对空气流场和水滴场进行非定常耦合求解，符合多步长计算的物理意义，实现真正的非定常多步长结冰过程。
附图说明
为了更清楚地说明本发明的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，对于本领域普通技术人员而言，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例提供的一种基于非定常欧拉两相流的多步长结冰计算方法的流程图；
图2为本发明实施例提供的结冰计算流程示意图；
图3为本发明实施例提供的控制体单元表面热量交换示意图；
图4为本发明实施例提供的三维曲面贴体坐标系图；
图5为本发明实施例提供的多步长结冰示意图；
图6为本发明实施例提供的DGRBF动网格方法流程示意图；
图7为本发明实施例提供的基于工况1得到的最终结冰结果图；
图8为本发明实施例提供的基于工况2得到的最终结冰结果图；
图9为本发明实施例提供的基于工况3得到的最终结冰结果图；
图10为本发明实施例提供的基于工况4得到的最终结冰结果图；
图11为本发明实施例提供的一种基于非定常欧拉两相流的多步长结冰计算系统的结构示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整的描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
在一实施例中，如图1所示，本发明实施例提供一种基于非定常欧拉两相流的多步长结冰计算方法，包括：
S1，获取空气流场初始条件和水滴场初始条件；所述空气流场初始条件包括来流压力、来流温度、来流速度和攻角；所述水滴场初始条件包括水滴平均容积直径和液态水含量。
结冰计算初始条件包括结冰计算总时间、结冰时间步长和结冰迭代时间步长。
S2，构建基于任意物体几何模型生成的计算网格。
计算的网格是基于物体几何信息由网格生成软件或自编程序生成的非结构化网格，网格文件格式为CGNS2.3形式；使用MPI和ParMETIS等库函数对计算网格进行并行分区，并对计算网格进行数据处理，包括对计算网格的所有单元进行编号，求出所有单元的体积、单元包含的每个面的面积及面左右单元间的梯度。
S3，将总结冰时间均匀划分为多个结冰时间步长。
假设结冰计算总时间为T total ，结冰时间步长为△T ice ，整个非定常结冰过程分为M个循环，即T total =M×△T ice 。
S4，根据空气流场初始条件和计算网格，计算每个结冰时间步长内空气流场的非定常数值。
在每个单结冰时间步长内，非定常耦合求解欧拉两相流（即空气流场和水滴场），采用隐式双时间步迭代求解（外迭代真实时间与结冰迭代时间步长相同），得出空气流场和水滴场的非定常积分结果。
结冰后物体表面形状变得极为复杂，明显的冰角特征会导致物面附近空气流场和水滴场出现大分离现象或涡结构，用定常求解两相流得到稳态收敛解是不准确的，为体现结冰后复杂外形导致空气流场和水滴场的非定常效应，本发明的方法使用隐式LU-SGS双时间步迭代非定常耦合求解空气流场和水滴场，假设每个单结冰时间步长为，隐式双时间步外迭代真实时间为△t，欧拉两相流计算循环次数为N，则△T ice =N×△t，得到N个循环空气流场和水滴场的非定常积分结果，包括物体表面附近空气流场速度梯度和温度梯度，表面局部水收集系数和收集水量。
构建空气流场控制方程；所述空气流场控制方程为理想状态气体可压缩雷诺平均Navier-Stokes方程；空气流场控制方程空间离散采用HLLC格式，时间离散使用隐式LU-SGS双时间步迭代，外迭代为真实时间步Δt，内迭代伪时间步长利用CFL条件求解。
S5，根据水滴场初始条件和计算网格，计算每个结冰时间步长内水滴场的非定常数值。
结冰后物体出现复杂外形（明显的冰角或冰瘤等），物面附近水滴场出现大分离现象或涡结构，这些非定常效应会导致表观水滴密度发生脉冲聚集，数值模拟易发散。原始水滴欧拉法方程组的非严格双曲性质易导致求解出现数值发散。将原始水滴欧拉法方程组分为适定双曲性质部分和源项，使方程求解稳定收敛，便于多步长非定常求解。
构建水滴场控制方程组的每个计算网格单元积分表达式：
Figure 128768DEST_PATH_IMAGE001
。
其中，t为计算时间；Ω为计算网格单元的体积；W d 为水滴场计算网格单元的体守恒通量；F d 为水滴场计算网格单元的面对流通量；S为计算网格单元其中一个面的面积；P和Q均为水滴场计算网格单元的源项。
Figure 379621DEST_PATH_IMAGE002
；
Figure 830325DEST_PATH_IMAGE003
；
Figure 901049DEST_PATH_IMAGE004
；
Figure 954456DEST_PATH_IMAGE005
。
其中，x、y和z分别为惯性坐标系的三个分量方向；ρ d 为水密度；α为水滴体积分数；ρ d α为表观水滴密度；u d 、v d 和w d 分别为水滴在x、y和z方向上的速度；V d 为水滴的逆变速度，V d =
Figure 59815DEST_PATH_IMAGE006
·
Figure 806054DEST_PATH_IMAGE007
；
Figure 364075DEST_PATH_IMAGE006
=（u d ，v d ，w d ）；计算单元的面单位法向量
Figure 96539DEST_PATH_IMAGE007
=（n x ，n y ，n z ）；n x ，n y 和n z 分别为面单位 法向量在x、y和z方向上的分量；g为重力加速度；d为水滴直径；K为空气和水滴作用因子；
Figure 321984DEST_PATH_IMAGE008
为质量守恒方程中的数值耗散项；θ为攻角。
水滴控制方程求解，将对流通量项右移，如下：
Figure 973545DEST_PATH_IMAGE024
。
其中，e表示计算网格单元有e个面；n k 表示计算网格单元的第k个面的单位法向量；ΔS k 表示计算网格单元的第k个面的面积。
根据以下公式计算水滴场计算网格单元的面对流通量：
Figure 18861DEST_PATH_IMAGE009
；
Figure 679650DEST_PATH_IMAGE010
；
Figure 634967DEST_PATH_IMAGE011
；
Figure 457430DEST_PATH_IMAGE012
；
Figure 990042DEST_PATH_IMAGE013
。
其中，F表示单元面对流通量；U表示单元体守恒通量；上标带*的变量是水滴场 HLLC格式构造的中间变量；下标为L的为面左侧单元的变量，下标为R的为面右侧单元的变 量；V L 为面左侧单元逆变速度；V L =
Figure 454522DEST_PATH_IMAGE014
·
Figure 123400DEST_PATH_IMAGE007
=
Figure 382343DEST_PATH_IMAGE015
；V R 为面右侧单元逆变速 度；V R =
Figure 776153DEST_PATH_IMAGE016
·
Figure 778744DEST_PATH_IMAGE007
=
Figure 567709DEST_PATH_IMAGE017
；S M 为中间波速；S L 为左波速；S R 为右波速；
S M =
Figure 997553DEST_PATH_IMAGE018
；
S L =V L -q L
Figure 239179DEST_PATH_IMAGE025
；
S R =V R +q R
Figure 45461DEST_PATH_IMAGE025
；
Figure 564298DEST_PATH_IMAGE020
；
Figure 165043DEST_PATH_IMAGE021
；
Figure 893965DEST_PATH_IMAGE022
。
其中，q L 为面左侧单元密度比；q R 为面右侧单元密度比。
上述空间离散格式具有良好的保正性质，满足表观水滴密度为正，使得方程组求解不易发散且便于多步长计算。
S6，根据空气流场和水滴场的非定常数值，计算当前结冰时间步长的结冰外形。
基于非定常水膜流动进行多个迭代时间步内水膜流动方程和结冰方程的计算，得出当前单结冰时间步长的结冰外形。结冰计算流程基于Myers结冰模型加以改进，以非定常水膜流动方程为核心，每一迭代时间步内，对所有表面控制体单元，均求解一次水膜流动方程组和结冰方程组，直到结冰迭代步结束得到当前单结冰时间步长的结冰外形，用于后续动网格模块进行网格重构。具体计算流程如图2所示。
将撞击过冷水滴全部冻结的状态定义为干结冰，部分冻结形成水膜的状态定义为湿结冰。对于霜冰工况，认为不出现湿结冰状态或占比很小，即无论是作用时间还是对应状态下的结冰量，干结冰状态都占绝对主导；明冰工况在初期也视作干结冰状态，但在很短的时间内就能转化为湿结冰状态，即湿结冰状态占主导；混合冰则介于两者之间。
无论是干结冰还是湿结冰状态，都涉及相变与传热问题。干结冰状态控制体和湿结冰状态控制体的表面热量交换如图3所示，图中T s表示壁面温度（等于冰层底部温度），T f表示控制体内最外层介质的表面温度，对于没有积冰和水膜存在的控制体，T f与T s等价，T i和T w分别表示冰层内部和水层内部的温度，T m为冰水界面温度。控制体的表面热量交换包括水滴撞击动能热流Q k，气动加热热流Q a，冻结潜热热流Q l，水滴显热热流Q d，对流换热热流Q c，蒸发热流Q evap，升华热流Q sub，辐射热流Q r等。
对水膜流动方程组的求解，基于欧拉两相流计算网格中的物体表面网格。三维曲面贴体坐标系如图4所示，贴体坐标系用(η,ζ,ξ)表示，η和ζ始终与曲面相切，ξ是曲面在该点的法向量，与η和ζ保持垂直，水膜在三维表面上沿η和ζ方向上流动，各方向上的单位宽度体积流量为F η 和F ζ 。方程组采用有限体积法进行离散，引入二阶迎风格式求解通量F η 和F ζ ；时间上采用三阶Runge-Kutta方法显式推进，计算时间步长取全局最小时间步长。
水膜流动方程组如下：
Figure 503938DEST_PATH_IMAGE026
。
其中，h w 表示水膜高度；LWC表示液态水含量；V 0表示水滴的远场速度，β表示局部水收集系数；m evap表示水膜蒸发质量；ρ i 表示冰密度；ρ w 表示水密度；p a 表示空气压力；V η 和V ζ 分别表示沿η和ζ方向的水流速度；τ η 和τ ζ 分别表示气流剪切力在沿η和ζ方向上的分量；vf η 和vf ζ 分别表示气流体积力在沿η和ζ方向上的分量。
结冰计算步骤具体如下：
步骤6.1，在每个迭代时间步内，首先判断结冰类型的依据由原Myers结冰模型的临界厚度改为当前表面温度T f ，T f ＜T m 表示干结冰状态（T m 为冰水混合物温度，通常为273.15K），反之则为湿结冰状态。
步骤6.2，根据结冰类型（干结冰状态或湿结冰状态）采用不同的质能平衡方程计算结冰增长率。
干结冰状态的控制体质量平衡方程（即结冰增长率计算式）如下：
Figure 1915DEST_PATH_IMAGE027
。
其中，t表示计算时间；h i 表示结冰厚度（即结冰增长率）；m imp 表示撞击水量；m in 表示上游溢流水量；m sub 表示冰升华质量。
干结冰状态表面能量平衡方程如下：
Figure 773562DEST_PATH_IMAGE028
。
其中，k i 表示冰的导热系数。
引入构造的中间变量E d 和e d ，则冰层温度分布表示为：
Figure 865146DEST_PATH_IMAGE029
；
其中，
Figure 13231DEST_PATH_IMAGE030
，
Figure 365715DEST_PATH_IMAGE031
，q d 表示水滴显 热系数；q c 表示对流换热系数；q sub 蒸发换热系数；q r 表示辐射换热系数。
湿结冰状态，控制体内表面热量平衡方程如下：
Figure 573842DEST_PATH_IMAGE032
；
其中，k w 表示冰的导热系数。
湿结冰状态的控制体质量平衡方程（结冰增长率由冰层和水层内部的导热决定）如下：
Figure 277356DEST_PATH_IMAGE033
；
其中，l f 表示结冰潜热。
引入构造的中间变量E w 和e w ，则水膜层温度分布表示为：
Figure 104498DEST_PATH_IMAGE034
；
其中，
Figure 311488DEST_PATH_IMAGE035
，
Figure 690517DEST_PATH_IMAGE036
；q evap 表示辐射换热 系数。
对于干结冰状态，则还需要用下式更新表面温度T f ；
表面温度T f 计算式为：
Figure 881327DEST_PATH_IMAGE037
。
步骤6.3，将结冰增长率作为源项加入上述的水膜流动方程组进行计算。
步骤6.4，每一单结冰的非定常迭代时间步长和两相流隐式双时间迭代求解的外迭代真实时间Δt相同，循环步骤6.1至步骤6.3，直到结冰迭代时间步长结束，得到当前单结冰时间步长的结冰外形，用于后续动网格模块进行网格重构。
改进的整体结冰计算流程不需要假设初始值（初始水膜高度和初始结冰厚度），迭代过程中能始终保持整个计算域冰水质量守恒。本发明结冰计算过程采用多步长方法，上一单结冰时间步长的结冰外形作为下一单结冰时间步长的干净翼型，上一单结冰时间步长的表面温度T f 用于初始化下一单结冰时间步长的壁面温度T s ，如图5所示。
S7，根据当前结冰时间步长的结冰外形，采用动网格技术改动计算网格节点以重构计算网格，将重构后的计算网格用于下一个结冰时间步长的计算过程。
动网格技术为基于Delaunay图映射和径向基函数的动网格算法，即DGRBF方法。
DGRBF方法以Delaunay图映射方法（即DGM方法）为基础，将最后求解映射关系的步骤改为通过位移或角度变化量（角变量）使用径向基函数（RBF）插值建立起Delaunay背景单元与网格节点之间的映射。DGRBF方法结合了RBF方法与DGM方法的优点，可以快速生成高质量的重构网格，且不改变网格节点数和网格单元数，便于后续多步长计算时插值赋值两相流结果。动网格方法流程示意图，如图6所示，具体步骤如下：
步骤7.1，从重构前的计算网格中，提取边界信息（包括远场边界、对称性边界、奇异性边界和周期性边界等），基于Delaunay准则使用Bowyer-Waston方法对边界点集进行划分，生成Delaunay背景图，此时的Delaunay背景图只用来记录计算网格节点之间的相对位置，保证在运动变形后能够正确地回插，而不必严格保持物面边界。
步骤7.2，使用阵面推进算法定位重构前的计算网格中各网格节点所属的背景图单元，在建立计算网格节点与Delaunay背景图之间的映射关系时，需要保证映射的唯一性，且所有节点都必须分类到背景图中的三角形（二维）或四面体（三维）单元内部。
步骤7.3，Delaunay背景图进行运动变形，基于运动变形产生的新Delaunay背景图和未产生运动变形时对应的初始Delaunay背景图，确定所有背景单元节点的位移。当几何体的运动或变形幅度过大时，可能会导致背景单元的节点位移尺度过大，从而在进行映射时产生负面积或负体积。此时，可将一次幅度较大的变形转化成多次幅度较小的变形，如果仍然出现负面积或负体积，则继续分成更多次幅度更小的变形，直到所有的变形步骤都没有产生问题。
步骤7.4，利用RBF插值建立起计算网格节点与Delaunay背景图之间的映射，Delaunay背景图节点变形前后的位移向量可以用RBF的加权和表示，即：
Figure 105635DEST_PATH_IMAGE038
；
其中，q为Delaunay背景单元的节点数，在二维情况下等于3，三维情况下等于4；α i 为系数，其值取决于具体的插值条件；γ i 是背景单元节点的坐标；φ是基函数，这里使用的是Wendland’s C2紧支型径向基函数，记为：φ（δ）=(1-δ)4(4δ+1)，其中δ是函数自变量。
将上述背景单元节点的坐标带入上式中，可得：
dt=Mt,tα；
其中dt=[s(γ1)，s(γ2)，┄，s(γq)]T表示Delaunay背景单元节点位移的列向量；α=[α1，α2，┄，αq]T是系数矩阵；Mt,t是一个q阶的矩阵；φ i,j =φ(||γ i -γ j ||)构成了矩阵元素，则有：
Figure 698290DEST_PATH_IMAGE039
；
得出系数矩阵α的计算公式：
Figure 356542DEST_PATH_IMAGE040
；
对于Delaunay图中任意单元中网格节点o的位移
Figure 34648DEST_PATH_IMAGE041
，由系 数矩阵α得到：
Figure 62647DEST_PATH_IMAGE042
。
其中，矩阵A=[φ o,1,φ o,2,┄,φ o,q ]；αx，αy和αz分别是α在笛卡尔坐标系中的沿坐标轴方向的三个分量。则对于网格节点o，其原始坐标为γo_1，则该点的新坐标γo_2为：
Figure 509809DEST_PATH_IMAGE043
；
最终，由所有网格节点新坐标得到重构后的计算网格。
S8，将重构前的计算网格的空气流场和水滴场瞬态结果插值赋值给重构后的计算网格。
采用多步长结冰方法时，因为动网格方法不改变网格节点数和网格单元数，将重构前网格的两相流结果插值赋值给重构后的新网格。插值方法由重构前后网格单元格心处两相流的体守恒变量W、网格单元格心处两相流的体守恒变量梯度∇W和限制器函数Φ计算完成，重构前网格单元的∇W和Φ都可以由邻居单元信息计算构造得到，都储存在网格单元的格心处。该插值方法从格心到格心，适用于任意网格单元类型，计算重构后新网格单元的变量的具体公式如下：
Figure 230640DEST_PATH_IMAGE044
。
其中，下标1的变量表示重构前计算网格单元的变量，下标2的变量表示重构后计算网格单元的变量；r 12表示网格单元1和网格单元2格心间的距离矢量。
随后，无需初始化直接对空气流场和水滴场进行非定常耦合求解，实现真正的非定常多步长结冰过程，符合结冰后复杂外形的空气流场和水滴场的非定常效应。
S9，重复步骤S2-S8，直至总结冰时间结束，得到最终结冰计算结果。
下面给出四个结冰工况算例作为本发明所公开方法的具体实施例，结冰具体工况如表1所示。
表1 四个结冰算例具体工况
Figure 864884DEST_PATH_IMAGE045
实施例一
NACA0012翼型典型霜冰结冰条件计算问题。以NACA0012翼型为研究对象，弦长为0.5334m，结冰工况如表1中工况1所示。工况1为典型霜冰结冰条件，以干结冰状态为主。单结冰时间步长取30s，隐式双时间求解欧拉两相流，外迭代真实时间取10-3s，内迭代CFL数取10，结冰模块迭代时间步长取10-3s。本发明方法得到的最终结冰结果图如图7所示，将本发明方法的多步长计算结果、单步长计算结果及LEWICE软件计算值与试验值对比，三个数值计算结果的结冰下极限位置都有偏差，本发明方法的多步长计算结果与试验值更为吻合，LEWICE软件计算结果的最大冰厚偏大，单步长计算结果总体结冰量偏大且在上翼面区域结冰量误差过大。该实施例说明了本发明在对称翼型典型霜冰结冰工况条件下的数值计算的有效性。
实施例二
NACA0012翼型典型明冰结冰条件计算问题。以NACA0012翼型为研究对象，弦长为0.5334m，结冰工况如表1中工况2所示。工况1为典型明冰结冰条件，以湿结冰状态为主。单结冰时间步长取30s，隐式双时间求解欧拉两相流，外迭代真实时间取10-3s，内迭代CFL数取10，结冰模块迭代时间步长取10-3s。本发明方法得到的最终结冰结果图如图8所示，将本发明方法的多步长计算结果、单步长计算结果及LEWICE软件计算值与试验值对比，从整体上看，本发明方法的多步长计算结果的结冰外形和最大冰厚与试验数据更接近，能够捕捉到冰角特征，上结冰极限位置略有误差，LEWICE软件计算结果在结冰上极限附近结冰量明显偏大，而单步长计算结果无法得出冰角特征且总体结冰量偏小。该实施例说明了本发明在对称翼型典型明冰结冰工况条件下的数值计算的有效性。
实施例三
非对称翼型正攻角明冰结冰条件计算问题。以一个非对称翼型为研究对象，弦长为0.5m，结冰工况如表1中工况3所示。工况3为明冰结冰条件，非对称翼型结冰以湿结冰状态为主。单结冰时间步长取30s，隐式双时间求解欧拉两相流，外迭代真实时间取10-3s，内迭代CFL数取10，结冰模块迭代时间步长取10-3s。本发明方法得到的最终结冰结果图如图9所示，展示了多步长计算总结冰时间第180s、540s和900s的计算结果，能够看出冰角特征的生长趋势和形成过程，正攻角条件下最大冰厚处靠下方，将多步长计算最终第900s的结果、单步长900s的计算结果与试验值相比，本发明方法的多步长计算结果的最大冰角位置和冰厚与试验值吻合度较高，结冰下极限位置与试验值极为接近，结冰上极限位置略有误差，而单步长计算结果的最大冰角位置偏下且总结冰量偏大，严重偏离试验值。该实施例说明了本发明在非对称翼型正攻角明冰结冰工况条件下的数值计算的有效性。
实施例四
非对称翼型负攻角明冰结冰条件计算问题。以一个非对称翼型为研究对象，弦长为0.5m，结冰工况如表1中工况4所示。工况4为明冰结冰条件，非对称翼型结冰以湿结冰状态为主。单结冰时间步长取30s，隐式双时间求解欧拉两相流，外迭代真实时间取10-3s，内迭代CFL数取10，结冰模块迭代时间步长取10-3s。本发明方法得到的最终结冰结果图如图10所示，展示了多步长计算总结冰时间第180s、540s和900s的计算结果，能够看出冰角特征的生长趋势和形成过程，负攻角条件下最大冰厚处靠上方，将多步长计算最终第900s的结果、单步长900s的计算结果与试验值相比，本发明方法的多步长计算结果的最大冰角位置和冰厚与试验值吻合度较高，结冰下极限位置与试验值极为接近，结冰上极限位置略有误差，而单步长计算结果的最大冰角位置偏下且总结冰量偏大，严重偏离试验值。该实施例说明了本发明在非对称翼型负攻角明冰结冰工况条件下的数值计算的有效性。
基于同一发明构思，本发明实施例还提供了一种基于非定常欧拉两相流的多步长结冰计算系统，由于该系统解决问题的原理与前述基于非定常欧拉两相流的多步长结冰计算方法相似，因此该系统的实施可以参见基于非定常欧拉两相流的多步长结冰计算方法的实施，重复之处不再赘述。
在另一实施例中，本发明一个实施例提供的基于非定常欧拉两相流的多步长结冰计算系统，如图11所示，包括：
获取模块10，用于获取空气流场初始条件和水滴场初始条件；所述空气流场初始条件包括来流压力、来流温度、来流速度和攻角；所述水滴场初始条件包括水滴平均容积直径和液态水含量。
构建模块20，用于构建基于任意物体几何模型生成的计算网格。
划分模块30，用于将总结冰时间均匀划分为多个结冰时间步长。
第一计算模块40，用于根据空气流场初始条件和计算网格，计算每个结冰时间步长内空气流场的非定常数值。
第二计算模块50，用于根据水滴场初始条件和计算网格，计算每个结冰时间步长内水滴场的非定常数值。
第三计算模块60，用于根据空气流场和水滴场的非定常数值，计算当前结冰时间步长的结冰外形。
网格重构模块70，用于根据当前结冰时间步长的结冰外形，采用动网格技术改动计算网格节点以重构计算网格，将重构后的计算网格用于下一个结冰时间步长的计算过程。
插值赋值模块80，用于将重构前的计算网格的空气流场和水滴场瞬态结果插值赋值给重构后的计算网格。
循环模块90，用于重复执行构建模块至插值赋值模块的操作，直至总结冰时间结束，得到最终结冰计算结果。
可选地，所述第一计算模块包括：
第一构建单元，用于构建空气流场控制方程；所述空气流场控制方程为理想状态气体可压缩雷诺平均Navier-Stokes方程；空气流场控制方程空间离散采用HLLC格式，时间离散使用隐式LU-SGS双时间步迭代，外迭代为真实时间步Δt，内迭代伪时间步长利用CFL条件求解。
可选地，所述第二计算模块包括：
第二构建单元，用于构建水滴场控制方程组的每个计算网格单元积分表达式：
Figure 571940DEST_PATH_IMAGE001
。
其中，t为计算时间；Ω为计算网格单元的体积；W d 为水滴场计算网格单元的体守恒通量；F d 为水滴场计算网格单元的面对流通量；S为计算网格单元其中一个面的面积；P和Q均为水滴场计算网格单元的源项；
Figure 139188DEST_PATH_IMAGE002
；
Figure 765341DEST_PATH_IMAGE003
；
Figure 152460DEST_PATH_IMAGE004
；
Figure 787841DEST_PATH_IMAGE005
。
其中，x、y和z分别为惯性坐标系的三个分量方向；ρ d 为水密度；α为水滴体积分数；ρ d α为表观水滴密度；u d 、v d 和w d 分别为水滴在x、y和z方向上的速度；V d 为水滴的逆变速度，V d =
Figure 944015DEST_PATH_IMAGE006
·
Figure 882016DEST_PATH_IMAGE007
；
Figure 756431DEST_PATH_IMAGE006
=（u d ，v d ，w d ）；计算单元的面单位法向量
Figure 664344DEST_PATH_IMAGE007
=（n x ，n y ，n z ）；n x ，n y 和n z 分别为面单位 法向量在x、y和z方向上的分量；g为重力加速度；d为水滴直径；K为空气和水滴作用因子；
Figure 940604DEST_PATH_IMAGE008
为质量守恒方程中的数值耗散项；θ为攻角。
第一计算单元，用于根据以下公式计算水滴场计算网格单元的面对流通量：
Figure 174140DEST_PATH_IMAGE009
；
Figure 411217DEST_PATH_IMAGE010
；
Figure 122821DEST_PATH_IMAGE011
；
Figure 253588DEST_PATH_IMAGE012
；
Figure 658025DEST_PATH_IMAGE013
。
其中，下标为L的为面左侧单元的变量，下标为R的为面右侧单元的变量；V L 为面左 侧单元逆变速度；V L =
Figure 241453DEST_PATH_IMAGE014
·
Figure 22327DEST_PATH_IMAGE007
=
Figure 381502DEST_PATH_IMAGE015
；V R 为面右侧单元逆变速度；V R =
Figure 956840DEST_PATH_IMAGE016
·
Figure 27564DEST_PATH_IMAGE007
=
Figure 346550DEST_PATH_IMAGE017
；S M 为中间波速；S L 为左波速；S R 为右波速；
S M =
Figure 451909DEST_PATH_IMAGE018
；
S L =V L -q L
Figure 932569DEST_PATH_IMAGE025
；
S R =V R +q R
Figure 631535DEST_PATH_IMAGE025
；
Figure 488632DEST_PATH_IMAGE020
；
Figure 448498DEST_PATH_IMAGE021
；
Figure 365638DEST_PATH_IMAGE022
。
其中，q L 为面左侧单元密度比；q R 为面右侧单元密度比。
可选地，所述第三计算模块包括：
获取单元，用于获取结冰类型。
第二计算单元，用于根据结冰类型采用不同的质能平衡方程计算结冰增长率。
第三计算单元，将结冰增长率作为源项带入至水膜流动方程组进行计算。
第四计算单元，用于循环计算结冰增长率并将结冰增长率作为源项带入至水膜流动方程组进行计算，直至结冰迭代时间步长结束，得到当前结冰时间步长的结冰外形。
关于上述各个模块更加具体的工作过程可以参考前述实施例公开的相应内容，在此不再进行赘述。
在另一实施例中，本发明提供一种计算机设备，包括处理器和存储器；其中，处理器执行存储器中保存的计算机程序时实现第一方面所述的基于非定常欧拉两相流的多步长结冰计算方法。
关于上述方法更加具体的过程可以参考前述实施例中公开的相应内容，在此不再进行赘述。
在另一实施例中，本发明提供一种计算机可读存储介质，用于存储计算机程序；计算机程序被处理器执行时实现第一方面所述的基于非定常欧拉两相流的多步长结冰计算方法。
关于上述方法更加具体的过程可以参考前述实施例中公开的相应内容，在此不再进行赘述。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其它实施例的不同之处，各个实施例之间相同或相似部分互相参见即可。对于实施例公开的系统、设备、存储介质而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
本领域的技术人员可以清楚地了解到本发明实施例中的技术可借助软件加必需的通用硬件平台的方式来实现。基于这样的理解，本发明实施例中的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备（可以是个人计算机，服务器，或者网络设备等）执行本发明各个实施例或者实施例的某些部分所述的方法。
以上结合具体实施方式和范例性实例对本发明进行了详细说明，不过这些说明并不能理解为对本发明的限制。本领域技术人员理解，在不偏离本发明精神和范围的情况下，可以对本发明技术方案及其实施方式进行多种等价替换、修饰或改进，这些均落入本发明的范围内。本发明的保护范围以所附权利要求为准。

Claims (10)

1.一种基于非定常欧拉两相流的多步长结冰计算方法，其特征在于，包括：
S1，获取空气流场初始条件和水滴场初始条件；所述空气流场初始条件包括来流压力、来流温度、来流速度和攻角；所述水滴场初始条件包括水滴平均容积直径和液态水含量；
S2，构建基于任意物体几何模型生成的计算网格；
S3，将总结冰时间均匀划分为多个结冰时间步长；
S4，根据空气流场初始条件和计算网格，计算每个结冰时间步长内空气流场的非定常数值；
S5，根据水滴场初始条件和计算网格，计算每个结冰时间步长内水滴场的非定常数值；
S6，根据空气流场和水滴场的非定常数值，计算当前结冰时间步长的结冰外形；
S7，根据当前结冰时间步长的结冰外形，采用动网格技术改动计算网格节点以重构计算网格，将重构后的计算网格用于下一个结冰时间步长的计算过程；
S8，将重构前的计算网格的空气流场和水滴场瞬态结果插值赋值给重构后的计算网格；
S9，重复步骤S2-S8，直至总结冰时间结束，得到最终结冰计算结果。
2.根据权利要求1所述的基于非定常欧拉两相流的多步长结冰计算方法，其特征在于，所述根据空气流场初始条件和计算网格，计算每个结冰时间步长内空气流场的非定常数值，包括：
构建空气流场控制方程；所述空气流场控制方程为理想状态气体可压缩雷诺平均Navier-Stokes方程；空气流场控制方程空间离散采用HLLC格式，时间离散使用隐式LU-SGS双时间步迭代，外迭代为真实时间步Δt，内迭代伪时间步长利用CFL条件求解。
3.根据权利要求1所述的基于非定常欧拉两相流的多步长结冰计算方法，其特征在于，所述根据水滴场初始条件和计算网格，计算每个结冰时间步长内水滴场的非定常数值，包括：
构建水滴场控制方程组的每个计算网格单元积分表达式：
Figure 198562DEST_PATH_IMAGE001
；
其中，t为计算时间；Ω为计算网格单元的体积；W d 为水滴场计算网格单元的体守恒通量；F d 为水滴场计算网格单元的面对流通量；S为计算网格单元其中一个面的面积；P和Q均为水滴场计算网格单元的源项；
Figure 714994DEST_PATH_IMAGE002
；
Figure 24753DEST_PATH_IMAGE003
；
Figure 626635DEST_PATH_IMAGE004
；
Figure 414463DEST_PATH_IMAGE005
；
其中，x、y和z分别为惯性坐标系的三个分量方向；ρ d 为水密度；α为水滴体积分数；ρ d α为 表观水滴密度；u d 、v d 和w d 分别为水滴在x、y和z方向上的速度；V d 为水滴的逆变速度，V d =
Figure 50981DEST_PATH_IMAGE006
·
Figure 266061DEST_PATH_IMAGE007
；
Figure 620819DEST_PATH_IMAGE006
=（u d ，v d ，w d ）；计算单元的面单位法向量
Figure 212337DEST_PATH_IMAGE007
=（n x ，n y ，n z ）；n x ，n y 和n z 分别为面单位法 向量在x、y和z方向上的分量；g为重力加速度；d为水滴直径；K为空气和水滴作用因子；
Figure 906624DEST_PATH_IMAGE008
为质量守恒方程中的数值耗散项；θ为攻角；
根据以下公式计算水滴场计算网格单元的面对流通量：
Figure 354923DEST_PATH_IMAGE009
；
Figure 869081DEST_PATH_IMAGE010
；
Figure 264290DEST_PATH_IMAGE011
；
Figure 875400DEST_PATH_IMAGE012
；
Figure 432283DEST_PATH_IMAGE013
；
其中，下标为L的为面左侧单元的变量，下标为R的为面右侧单元的变量；V L 为面左侧单 元逆变速度；V L =
Figure 496054DEST_PATH_IMAGE014
·
Figure 429375DEST_PATH_IMAGE007
=
Figure 98254DEST_PATH_IMAGE015
；V R 为面右侧单元逆变速度；V R =
Figure 888355DEST_PATH_IMAGE016
·
Figure 377106DEST_PATH_IMAGE007
=
Figure 176434DEST_PATH_IMAGE017
；S M 为中间波速；S L 为左波速；S R 为右波速；
S M =
Figure 434240DEST_PATH_IMAGE018
；
S L =V L -q L
Figure 598505DEST_PATH_IMAGE019
；
S R =V R +q R
Figure 636869DEST_PATH_IMAGE019
；
Figure 911992DEST_PATH_IMAGE020
；
Figure 86621DEST_PATH_IMAGE021
；
Figure 687367DEST_PATH_IMAGE022
；
其中，q L 为面左侧单元密度比；q R 为面右侧单元密度比。
4.根据权利要求1所述的基于非定常欧拉两相流的多步长结冰计算方法，其特征在于，所述根据空气流场和水滴场的非定常数值，计算当前结冰时间步长的结冰外形，包括：
获取结冰类型；
根据结冰类型采用不同的质能平衡方程计算结冰增长率；
将结冰增长率作为源项带入至水膜流动方程组进行计算；
循环计算结冰增长率并将结冰增长率作为源项带入至水膜流动方程组进行计算，直至结冰迭代时间步长结束，得到当前结冰时间步长的结冰外形。
5.一种基于非定常欧拉两相流的多步长结冰计算系统，其特征在于，包括：
获取模块，用于获取空气流场初始条件和水滴场初始条件；所述空气流场初始条件包括来流压力、来流温度、来流速度和攻角；所述水滴场初始条件包括水滴平均容积直径和液态水含量；
构建模块，用于构建基于任意物体几何模型生成的计算网格；
划分模块，用于将总结冰时间均匀划分为多个结冰时间步长；
第一计算模块，用于根据空气流场初始条件和计算网格，计算每个结冰时间步长内空气流场的非定常数值；
第二计算模块，用于根据水滴场初始条件和计算网格，计算每个结冰时间步长内水滴场的非定常数值；
第三计算模块，用于根据空气流场和水滴场的非定常数值，计算当前结冰时间步长的结冰外形；
网格重构模块，用于根据当前结冰时间步长的结冰外形，采用动网格技术改动计算网格节点以重构计算网格，将重构后的计算网格用于下一个结冰时间步长的计算过程；
插值赋值模块，用于将重构前的计算网格的空气流场和水滴场瞬态结果插值赋值给重构后的计算网格；
循环模块，用于重复执行构建模块至插值赋值模块的操作，直至总结冰时间结束，得到最终结冰计算结果。
6.根据权利要求5所述的基于非定常欧拉两相流的多步长结冰计算系统，其特征在于，所述第一计算模块包括：
第一构建单元，用于构建空气流场控制方程；所述空气流场控制方程为理想状态气体可压缩雷诺平均Navier-Stokes方程；空气流场控制方程空间离散采用HLLC格式，时间离散使用隐式LU-SGS双时间步迭代，外迭代为真实时间步Δt，内迭代伪时间步长利用CFL条件求解。
7.根据权利要求5所述的基于非定常欧拉两相流的多步长结冰计算系统，其特征在于，所述第二计算模块包括：
第二构建单元，用于构建水滴场控制方程组的每个计算网格单元积分表达式：
Figure 885130DEST_PATH_IMAGE001
；
其中，t为计算时间；Ω为计算网格单元的体积；W d 为水滴场计算网格单元的体守恒通量；F d 为水滴场计算网格单元的面对流通量；S为计算网格单元其中一个面的面积；P和Q均为水滴场计算网格单元的源项；
Figure 26262DEST_PATH_IMAGE002
；
Figure 993081DEST_PATH_IMAGE003
；
Figure 764727DEST_PATH_IMAGE004
；
Figure 512104DEST_PATH_IMAGE005
；
其中，x、y和z分别为惯性坐标系的三个分量方向；ρ d 为水密度；α为水滴体积分数；ρ d α为 表观水滴密度；u d 、v d 和w d 分别为水滴在x、y和z方向上的速度；V d 为水滴的逆变速度，V d =
Figure 129030DEST_PATH_IMAGE006
·
Figure 281181DEST_PATH_IMAGE007
；
Figure 958150DEST_PATH_IMAGE006
=（u d ，v d ，w d ）；计算单元的面单位法向量
Figure 192822DEST_PATH_IMAGE023
=（n x ，n y ，n z ）；n x ，n y 和n z 分别为面单位法向 量在x、y和z方向上的分量；g为重力加速度；d为水滴直径；K为空气和水滴作用因子；
Figure 879019DEST_PATH_IMAGE024
为质量守恒方程中的数值耗散项；θ为攻角；
第一计算单元，用于根据以下公式计算水滴场计算网格单元的面对流通量：
Figure 820430DEST_PATH_IMAGE009
；
Figure 730617DEST_PATH_IMAGE010
；
Figure 390268DEST_PATH_IMAGE011
；
Figure 614576DEST_PATH_IMAGE012
；
Figure 738390DEST_PATH_IMAGE013
；
其中，下标为L的为面左侧单元的变量，下标为R的为面右侧单元的变量；V L 为面左侧单 元逆变速度；V L =
Figure 491583DEST_PATH_IMAGE014
·
Figure 700847DEST_PATH_IMAGE007
=
Figure 728846DEST_PATH_IMAGE015
；V R 为面右侧单元逆变速度；V R =
Figure 644849DEST_PATH_IMAGE016
·
Figure 896839DEST_PATH_IMAGE007
=
Figure 265503DEST_PATH_IMAGE017
；S M 为中间波速；S L 为左波速；S R 为右波速；
S M =
Figure 628352DEST_PATH_IMAGE018
；
S L =V L -q L
Figure 930020DEST_PATH_IMAGE025
；
S R =V R +q R
Figure 290594DEST_PATH_IMAGE025
；
Figure 208872DEST_PATH_IMAGE020
；
Figure 313094DEST_PATH_IMAGE021
；
Figure 469269DEST_PATH_IMAGE022
；
其中，q L 为面左侧单元密度比；q R 为面右侧单元密度比。
8.根据权利要求5所述的基于非定常欧拉两相流的多步长结冰计算系统，其特征在于，所述第三计算模块包括：
获取单元，用于获取结冰类型；
第二计算单元，用于根据结冰类型采用不同的质能平衡方程计算结冰增长率；
第三计算单元，将结冰增长率作为源项带入至水膜流动方程组进行计算；
第四计算单元，用于循环计算结冰增长率并将结冰增长率作为源项带入至水膜流动方程组进行计算，直至结冰迭代时间步长结束，得到当前结冰时间步长的结冰外形。
9.一种计算机设备，其特征在于，包括处理器和存储器；其中，处理器执行存储器中保存的计算机程序时实现权利要求1-4任一项所述的基于非定常欧拉两相流的多步长结冰计算方法的步骤。
10.一种计算机可读存储介质，其特征在于，用于存储计算机程序；计算机程序被处理器执行时实现权利要求1-4任一项所述的基于非定常欧拉两相流的多步长结冰计算方法的步骤。
CN202211401890.8A 2022-11-10 2022-11-10 一种基于非定常欧拉两相流的多步长结冰计算方法和系统 Active CN115563906B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211401890.8A CN115563906B (zh) 2022-11-10 2022-11-10 一种基于非定常欧拉两相流的多步长结冰计算方法和系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211401890.8A CN115563906B (zh) 2022-11-10 2022-11-10 一种基于非定常欧拉两相流的多步长结冰计算方法和系统

Publications (2)

Publication Number Publication Date
CN115563906A true CN115563906A (zh) 2023-01-03
CN115563906B CN115563906B (zh) 2023-03-24

Family

ID=84769452

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211401890.8A Active CN115563906B (zh) 2022-11-10 2022-11-10 一种基于非定常欧拉两相流的多步长结冰计算方法和系统

Country Status (1)

Country Link
CN (1) CN115563906B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116562192A (zh) * 2023-07-06 2023-08-08 中国空气动力研究与发展中心计算空气动力研究所 一种飞机结冰冰形预测方法、装置、设备及存储介质
CN117875224A (zh) * 2024-03-12 2024-04-12 南京航空航天大学 一种非等分连续多粒径水滴耦合欧拉方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111396269A (zh) * 2020-06-08 2020-07-10 中国空气动力研究与发展中心低速空气动力研究所 一种多时间步非定常结冰计算方法、系统及存储介质
CN113434978A (zh) * 2021-06-28 2021-09-24 成都安世亚太科技有限公司 一种航空飞行器结冰可靠性智能评估方法
CN114462330A (zh) * 2022-01-19 2022-05-10 清华大学 飞机结冰冰形预测方法、装置、计算机设备和存储介质

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111396269A (zh) * 2020-06-08 2020-07-10 中国空气动力研究与发展中心低速空气动力研究所 一种多时间步非定常结冰计算方法、系统及存储介质
CN113434978A (zh) * 2021-06-28 2021-09-24 成都安世亚太科技有限公司 一种航空飞行器结冰可靠性智能评估方法
CN114462330A (zh) * 2022-01-19 2022-05-10 清华大学 飞机结冰冰形预测方法、装置、计算机设备和存储介质

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116562192A (zh) * 2023-07-06 2023-08-08 中国空气动力研究与发展中心计算空气动力研究所 一种飞机结冰冰形预测方法、装置、设备及存储介质
CN116562192B (zh) * 2023-07-06 2023-09-12 中国空气动力研究与发展中心计算空气动力研究所 一种飞机结冰冰形预测方法、装置、设备及存储介质
CN117875224A (zh) * 2024-03-12 2024-04-12 南京航空航天大学 一种非等分连续多粒径水滴耦合欧拉方法
CN117875224B (zh) * 2024-03-12 2024-06-04 南京航空航天大学 一种非等分连续多粒径水滴耦合欧拉方法

Also Published As

Publication number Publication date
CN115563906B (zh) 2023-03-24

Similar Documents

Publication Publication Date Title
CN115563906B (zh) 一种基于非定常欧拉两相流的多步长结冰计算方法和系统
CN109376403B (zh) 一种基于笛卡尔自适应重构技术的二维结冰数值模拟方法
Perić Analysis of pressure-velocity coupling on nonorthogonal grids
Demirdzic et al. A calculation procedure for turbulent flow in complex geometries
CN107220399A (zh) 基于埃尔米特插值基本加权无振荡格式的全流场模拟方法
CN114168796B (zh) 一种建立飞行器高空气动力数据库的方法
CN112818573B (zh) 一种用于非结构网格的获取边界层非当地变量信息的方法
CN115659517B (zh) 一种旋翼桨叶结冰准非定常数值模拟方法和系统
Wang et al. Multi-body separation simulation with an improved general mesh deformation method
CN109726433A (zh) 基于曲面边界条件的三维无粘低速绕流的数值模拟方法
CN112613243B (zh) 一种流体力学模拟的方法、装置及计算机可读存储介质
CN105760602A (zh) 一种有限体积加权基本无振荡格式的全流场数值模拟方法
CN109726465A (zh) 基于非结构曲边网格的三维无粘低速绕流的数值模拟方法
Hui The unified coordinate system in computational fluid dynamics
Ghoreishi et al. Vorticity-based polynomial adaptation for moving and deforming domains
CN115803744A (zh) 物理系统的计算分析
CN114385960A (zh) 一种基于能量平均温度的间壁式换热器性能计算方法
Thompson et al. Discrete surface evolution and mesh deformation for aircraft icing applications
Mathur Unsteady flow simulations using unstructured sliding meshes
CN109740281A (zh) 一种基于蒸汽流场计算液滴运动相变参数的方法
Bourgault-Cote et al. Two-dimensional/infinite swept wing ice accretion model
CN114611423A (zh) 一种三维多相可压缩流固耦合的快速计算方法
Byun et al. Wing-body aeroelasticity using finite-difference fluid/finite element structural equations on parallel computers
Koren Euler flow solutions for a transonic windtunnel section
Gudmundsson Performance evaluation of wet-cooling tower fills with computational fluid dynamics

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant