CN116049997B - 飞行器表面烧蚀模拟混合处理方法、装置、设备及介质 - Google Patents

飞行器表面烧蚀模拟混合处理方法、装置、设备及介质 Download PDF

Info

Publication number
CN116049997B
CN116049997B CN202310319293.9A CN202310319293A CN116049997B CN 116049997 B CN116049997 B CN 116049997B CN 202310319293 A CN202310319293 A CN 202310319293A CN 116049997 B CN116049997 B CN 116049997B
Authority
CN
China
Prior art keywords
current
gas
mass
component
flow characteristic
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.)
Active
Application number
CN202310319293.9A
Other languages
English (en)
Other versions
CN116049997A (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.)
Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Original Assignee
Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
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 Computational Aerodynamics Institute of China Aerodynamics Research and Development Center filed Critical Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Priority to CN202310319293.9A priority Critical patent/CN116049997B/zh
Publication of CN116049997A publication Critical patent/CN116049997A/zh
Application granted granted Critical
Publication of CN116049997B publication Critical patent/CN116049997B/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/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • 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
    • 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/08—Thermal analysis or thermal optimisation
    • 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
    • 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)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • Mathematical Optimization (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Evolutionary Computation (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Computational Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Physics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本申请公开了一种飞行器表面烧蚀模拟混合处理方法、装置、设备及介质，涉及数值模拟计算和气动热防护领域，具体包括：获取第一流动特征参数；基于流动机理和烧蚀机制，计算第二流动特征参数；基于小扰动方法，计算相关偏导数；基于显隐混合方法，求解质量和能量平衡方程，得到所需目标参数。该方法从表面微元烧蚀机制及其质量和能量守恒的基本原理出发构建了表面烧蚀迭代模拟体系，依据各物理量在烧蚀体系中影响程度的差异，混合运用了显示和隐式相结合离散处理方案，实现了各气体组分质量平衡方程、能量平衡方程相对独立的隐式迭代计算，无需求解耦合矩阵，易于实现、稳定性好且计算效率高。

Description

飞行器表面烧蚀模拟混合处理方法、装置、设备及介质
技术领域
本发明涉及数值模拟计算和气动热防护领域，特别涉及一种飞行器表面烧蚀模拟混合处理方法、装置、设备及介质。
背景技术
高超声速飞行器再入飞行面临严酷的气动加热，飞行器表面常采用烧蚀材料进行热防护设计。烧蚀材料在高温下会发生热解、熔化、蒸发、侵蚀等复杂物理和化学变化, 其产物引射至表面边界层，与边界层内高温气体发生化学反应，从而改变飞行器绕流和尾迹的流动特性及其目标特性。由于其复杂性和重要性，高超声速飞行器表面烧蚀的数值模拟研究，一直受到人们的广泛关注。
高超声速飞行器表面烧蚀数值模拟的种类和形式有很多，既有较为单一的静态烧蚀过程模拟研究，也有考虑高超声速流动与烧蚀过程相互影响的耦合模拟研究。在这些研究中，通过求解表面微元质量和能量平衡方程获取烧蚀边界特征，是飞行器表面烧蚀数值模拟过程中最核心也是最常见的处理方法之一。
由于表面微元质量和能量平衡方程是一组非线性方程，无法直接解析得到计算结果，通常采用数值迭代模拟。其迭代方法大体可以分为两类：显示方法和隐式方法；其中，显示迭代方法实现较为简单，但稳定性相对较差，计算步长受到较大限制，迭代收敛较为缓慢，容易发散；而隐式方法虽然稳定性得到较大幅度的改善，但实现过程相对复杂，通常需要进行矩阵求逆或者矩阵耦合求解，且单步迭代耗时相对较长。
综上，如何提供一种稳定性较强的飞行器表面烧蚀模拟方案，并提高计算效率是目前有待解决的问题。
发明内容
有鉴于此，本发明的目的在于提供一种飞行器表面烧蚀模拟混合处理方法、装置、设备及介质，能够提供一种稳定性较强的飞行器表面烧蚀模拟方案，并提高计算效率。其具体方案如下：
第一方面，本申请公开了一种飞行器表面烧蚀模拟混合处理方法，包括：
在高超声速飞行器表面烧蚀过程中获取飞行器固壁面上与当前迭代步对应的第一流动特性参数；所述第一流动特性参数包括当前气体温度和当前气体各组分质量分数；
基于表面压力边界条件、烧蚀机制、气体状态方程、气体质量扩散和能力运输机制，并利用所述第一流动特性参数确定出与所述当前迭代步对应的第二流动特性参数；所述第二流动特性参数包括当前气体压强、当前气体各组分质量生成率、当前总质量生成率、当前总能量生成率、当前气体密度、当前气体组分扩散系数和当前热传导系数；
利用所述当前气体各组分质量分数和对应的组分质量扰动量确定出所述当前气体各组分质量生成率关于气体各组分质量分数的第一偏导数，并利用所述当前气体温度和对应的温度扰动量确定出所述当前总能量生成率关于气体温度的第二偏导数；
基于所述第一流动特性参数、所述第二流动特性参数、所述第一偏导数和所述第二偏导数，并按照预设的隐式处理和显示处理相结合的方法对质量平衡方程和能量平衡方程分别进行迭代求解，直到满足预设收敛条件，以得到迭代后的目标气体温度和目标气体各组分质量分数。
可选的，所述基于表面压力边界条件、烧蚀机制、气体状态方程、气体质量扩散和能力运输机制，并利用所述第一流动特性参数确定出与所述当前迭代步对应的第二流动特性参数，包括：
基于飞行器表面压力边界条件，并利用所述第一流动特性参数确定出与所述当前迭代步对应的当前气体压强；
基于飞行器表面材料的烧蚀机制，并利用所述第一流动特性参数和所述当前气体压强确定出与所述当前迭代步对应的当前气体各组分质量生成率、当前总质量生成率和当前总能量生成率；
基于气体状态方程，并利用所述第一流动特性参数和所述当前气体压强确定出与所述当前迭代步对应的当前气体密度；
基于气体质量扩散和能力运输机制，并利用所述第一流动特性参数、所述当前气体压强和所述当前气体密度确定出与所述当前迭代步对应的当前气体组分扩散系数和当前热传导系数；
基于所述当前气体压强、所述当前气体各组分质量生成率、所述当前总质量生成率、所述当前总能量生成率、所述当前气体密度、所述当前气体组分扩散系数和所述当前热传导系数构成第二流动特性参数。
可选的，所述基于飞行器表面材料的烧蚀机制，并利用所述第一流动特性参数和所述当前气体压强确定出与所述当前迭代步对应的当前气体各组分质量生成率、当前总质量生成率和当前总能量生成率，包括：
基于飞行器表面材料的烧蚀机制构建气体各组分质量生成率函数、总质量生成率函数和总能量生成率函数；
利用所述第一流动特性参数、所述当前气体压强、所述气体各组分质量生成率函数、所述总质量生成率函数和所述总能量生成率函数确定出与所述当前迭代步对应的当前气体各组分质量生成率、当前总质量生成率和当前总能量生成率。
可选的，所述利用所述当前气体各组分质量分数和对应的组分质量扰动量确定出所述当前气体各组分质量生成率关于气体各组分质量分数的第一偏导数，并利用所述当前气体温度和对应的温度扰动量确定出所述当前总能量生成率关于气体温度的第二偏导数，包括：
基于所述当前气体温度和预设偏导数调节系数确定出对应的温度扰动量，并基于所述当前气体各组分质量分数和所述预设偏导数调节系数确定出对应的组分质量扰动量；
利用所述当前气体各组分质量分数和所述组分质量扰动量确定出所述当前气体各组分质量生成率关于气体各组分质量分数的第一偏导数，并利用所述当前气体温度和所述温度扰动量确定出所述总能量生成率关于气体温度的第二偏导数。
可选的，所述在高超声速飞行器表面烧蚀过程中获取飞行器固壁面上与当前迭代步对应的第一流动特性参数的过程中，还包括：
获取流动界面上的输入气体温度和输入气体各组分质量分数。
可选的，所述基于所述第一流动特性参数、所述第二流动特性参数、所述第一偏导数和所述第二偏导数，并按照预设的隐式处理和显示处理相结合的方法对质量平衡方程和能量平衡方程分别进行迭代求解，包括：
确定出包括所述输入气体各组分质量分数、所述第一偏导数和表面网格微元法向厚度的质量平衡方程；
按照预设的隐式处理方法对所述质量平衡方程中的气体各组分质量分数和气体各组分质量生成率进行处理，并按照预设的显示处理方法利用所述第一流动特性参数和所述第二流动特性参数对所述质量平衡方程中其余参数进行处理，以对所述质量平衡方程进行迭代求解；
确定出包括所述输入气体温度、所述输入气体各组分质量分数、所述第二偏导数、所述表面网格微元法向厚度和气体各组分生成焓的能量平衡方程；
按照所述隐式处理方法对所述能量平衡方程中的气体温度、气体各组分质量分数和总能量生成率进行处理，并按照所述显示处理方法利用所述第一流动特性参数和所述第二流动特性参数对所述能量平衡方程中其余参数进行处理，以对所述能量平衡方程进行迭代求解。
可选的，所述按照预设的隐式处理和显示处理相结合的方法对质量平衡方程和能量平衡方程分别进行迭代求解，直到满足预设收敛条件，以得到迭代后的目标气体温度和目标气体各组分质量分数，包括：
按照预设的隐式处理和显示处理相结合的方法对质量平衡方程和能量平衡方程分别进行迭代求解，并判断迭代相对残差值是否小于预设残差阈值；
若所述迭代相对残差值小于所述预设残差阈值，则输出与下一迭代步对应的目标气体温度和目标气体各组分质量分数；
若所述迭代相对残差值大于所述预设残差阈值，则将所述下一迭代步作为当前迭代步，并重新跳转至所述在高超声速飞行器表面烧蚀过程中获取飞行器固壁面上与当前迭代步对应的第一流动特性参数的步骤，直到所述迭代相对残差值小于所述预设残差阈值，以得到迭代后的目标气体温度和目标气体各组分质量分数。
第二方面，本申请公开了一种飞行器表面烧蚀模拟混合处理装置，包括：
第一参数获取模块，用于在高超声速飞行器表面烧蚀过程中获取飞行器固壁面上与当前迭代步对应的第一流动特性参数；所述第一流动特性参数包括当前气体温度和当前气体各组分质量分数；
第二参数获取模块，用于基于表面压力边界条件、烧蚀机制、气体状态方程、气体质量扩散和能力运输机制，并利用所述第一流动特性参数确定出与所述当前迭代步对应的第二流动特性参数；所述第二流动特性参数包括当前气体压强、当前气体各组分质量生成率、当前总质量生成率、当前总能量生成率、当前气体密度、当前气体组分扩散系数和当前热传导系数；
偏导数确定模块，用于利用所述当前气体各组分质量分数和对应的组分质量扰动量确定出所述当前气体各组分质量生成率关于气体各组分质量分数的第一偏导数，并利用所述当前气体温度和对应的温度扰动量确定出所述当前总能量生成率关于气体温度的第二偏导数；
迭代求解模块，用于基于所述第一流动特性参数、所述第二流动特性参数、所述第一偏导数和所述第二偏导数，并按照预设的隐式处理和显示处理相结合的方法对质量平衡方程和能量平衡方程分别进行迭代求解，直到满足预设收敛条件，以得到迭代后的目标气体温度和目标气体各组分质量分数。
第三方面，本申请公开了一种电子设备，包括：
存储器，用于保存计算机程序；
处理器，用于执行所述计算机程序，以实现前述公开的飞行器表面烧蚀模拟混合处理方法的步骤。
第四方面，本申请公开了一种计算机可读存储介质，用于存储计算机程序；其中，所述计算机程序被处理器执行时实现前述公开的飞行器表面烧蚀模拟混合处理方法的步骤。
可见，本申请通过在高超声速飞行器表面烧蚀过程中获取飞行器固壁面上与当前迭代步对应的第一流动特性参数；所述第一流动特性参数包括当前气体温度和当前气体各组分质量分数；基于表面压力边界条件、烧蚀机制、气体状态方程、气体质量扩散和能力运输机制，并利用所述第一流动特性参数确定出与所述当前迭代步对应的第二流动特性参数；所述第二流动特性参数包括当前气体压强、当前气体各组分质量生成率、当前总质量生成率、当前总能量生成率、当前气体密度、当前气体组分扩散系数和当前热传导系数；利用所述当前气体各组分质量分数和对应的组分质量扰动量确定出所述当前气体各组分质量生成率关于气体各组分质量分数的第一偏导数，并利用所述当前气体温度和对应的温度扰动量确定出所述当前总能量生成率关于气体温度的第二偏导数；基于所述第一流动特性参数、所述第二流动特性参数、所述第一偏导数和所述第二偏导数，并按照预设的隐式处理和显示处理相结合的方法对质量平衡方程和能量平衡方程分别进行迭代求解，直到满足预设收敛条件，以得到迭代后的目标气体温度和目标气体各组分质量分数。由此可见，本申请在高超声速飞行器表面烧蚀过程中首先获取到飞行器固壁面上的流动特性参数，包括与当前迭代步对应的当前气体温度和当前气体各组分质量分数；再结合表面压力边界条件、烧蚀机制、气体状态方程、气体质量扩散和能力运输机制，利用当前气体温度和当前气体各组分质量分数确定出当前迭代步对应的当前气体压强、当前气体各组分质量生成率、当前总质量生成率、当前总能量生成率、当前气体密度、当前气体组分扩散系数和当前热传导系数等其他流动特性参数；然后在当前气体各组分质量分数和当前气体温度的基础上，引入各自对应的组分质量扰动量和温度扰动量，以分别确定出当前气体各组分质量生成率关于气体各组分质量分数的第一偏导数，以及当前总能量生成率关于气体温度的第二偏导数，通过引入小扰动实现偏导数的近似计算，不仅不需要推导偏导数具体的数学形式，实现简便，而且可广泛适用于各类材料的多种烧蚀机制；最后基于各流动特性参数、第一偏导数和第二偏导数，并按照预设的隐式处理和显示处理相结合的方法混合对质量平衡方程和能量平衡方程分别进行迭代求解，直到满足收敛条件时得到迭代后输出的固壁面上的目标气体温度和目标气体各组分质量分数。上述方案将隐式处理和显示处理相结合，混合运用显式和隐式离散处理方案，实现各气体组分质量平衡方程、能量平衡方程相对独立的隐式迭代计算，无需求解耦合矩阵，实现简便，稳定性好且计算效率高。
附图说明
为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。
图1为本申请公开的一种飞行器表面烧蚀模拟混合处理方法流程图；
图2为本申请公开的一种飞行器表面烧蚀模拟显隐混合计算处理流程图；
图3为本申请公开的一种具体的飞行器表面网格微元示意图；
图4为本申请公开的一种具体的飞行器表面烧蚀模拟混合处理方法流程图；
图5为本申请公开的一种驻点线主要气体组分质量分数分布示意图；
图6为本申请公开的一种现有文献中的驻点线主要气体组分质量分数分布示意图；
图7为本申请公开的一种碳-碳材料烧蚀条件下表面热流分布示意图；
图8为本申请公开的一种飞行器表面烧蚀模拟混合处理装置结构示意图；
图9为本申请公开的一种电子设备结构图。
具体实施方式
下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
当前数值迭代模拟的迭代方法大体可以分为两类：显示方法和隐式方法；其中，显示迭代方法实现较为简单，但稳定性相对较差，计算步长受到较大限制，迭代收敛较为缓慢，容易发散；而隐式方法虽然稳定性得到较大幅度的改善，但实现过程相对复杂，通常需要进行矩阵求逆或者矩阵耦合求解，且单步迭代耗时相对较长。为此，本申请实施例公开了一种飞行器表面烧蚀模拟混合处理方法、装置、设备及介质，能够提供一种稳定性较强的飞行器表面烧蚀模拟方案，并提高计算效率。
参见图1和图2所示，本申请实施例公开了一种飞行器表面烧蚀模拟混合处理方法，该方法包括：
步骤S11：在高超声速飞行器表面烧蚀过程中获取飞行器固壁面上与当前迭代步对应的第一流动特性参数；所述第一流动特性参数包括当前气体温度和当前气体各组分质量分数。
本实施例中，需要指出的是，上述在高超声速飞行器表面烧蚀过程中获取飞行器固壁面上与当前迭代步对应的第一流动特性参数的过程中，还包括：获取流动界面上的输入气体温度和输入气体各组分质量分数。可以理解的是，首先为取值步骤，导入数值计算网格，离散飞行器表面区域，表面网格微元上沿为流动界面，下沿为飞行器固壁面；表面网格微元法向厚度为h，h须足够小，使各参数在表面法向满足近似线性关系；流动界面上气体压强P up 、气体温度T up 和气体各组分质量分数
Figure SMS_1
由外部输入得到，在耦合流场模拟时，由高超声速飞行器流场数值模拟提供，这里/>
Figure SMS_2
，N s 为混合气体组分总个数。固壁面上的主要流动特性参数包括气体温度和气体各组分质量分数，因此在高超声速飞行器表面烧蚀过程中需要获取飞行器固壁面上与当前迭代步对应的第一流动特性参数，第一流动特性参数包括当前气体温度和当前气体各组分质量分数。令当前迭代步为第n步，则固壁面上的第一流动特性参数包括第n步迭代时表面气体温度T n 和气体各组分质量分数/>
Figure SMS_3
。n为迭代步数。当n=0，即迭代开始前，/>
Figure SMS_4
、/>
Figure SMS_5
，/>
Figure SMS_6
和/>
Figure SMS_7
为飞行器表面气体温度和气体各组分质量分数的初值。
可以理解的是，高超声速飞行器表面烧蚀数值模拟通常需要数值迭代求解“表面烧蚀能量平衡方程”和“表面烧蚀各气体组分质量平衡方程”。由于数值迭代过程需要在空间离散化的表面网格微元上进行，因此需要针对飞行器外形，采用公开的网格软件或公开的网格技术生成合适的表面网格，对飞行器表面进行离散，并将其数据化，为数值迭代做好准备。在飞行器流场与表面烧蚀效应耦合模拟时，上述表面网格通常是飞行器流场模拟采用网格的飞行器表面法向第一层网格。离散之后的网格微元及其参数如图3中所示，需要说明的是，图3为二维的简化示意图，仅为了便于理解，并不代表真实的网格。
由于微元法向厚度h足够小，满足“薄层近似”，在飞行器固壁面切向方向的能量和质量输运可以忽略，因此只需要关注表面网格微元上沿流动界面和下沿飞行器固壁面的能量、质量输运过程，数值模拟时仅需获取表面网格微元上沿和下沿的部分参数。表面网格微元上沿参数，一般由外部输入得到，在飞行器流场与表面烧蚀效应耦合模拟时，该参数可由流场模拟得到。飞行器流场与表面烧蚀效应耦合模拟，一般采用双重迭代进行：飞行器流场数值模拟进行父迭代，在每一步父迭代过程中，基于固壁面上的参数，流场模拟给出表面网格微元上沿流动界面上的参数；然后在表面网格微元上进行子迭代，子迭代求解能量平衡方程和质量平衡方程，得到飞行器固壁面上新的参数；最后基于固壁面上的新参数，进行流场的下一步父迭代，直至父迭代收敛。而本申请主要针对其中的子迭代过程，即基于网格微元上沿流动界面上的参数（即P up 、T up 和
Figure SMS_8
）和固壁面上温度、气体组分的旧参数（即初值/>
Figure SMS_9
和/>
Figure SMS_10
），通过高效的显式与隐式混合迭代求解表面烧蚀能量平衡方程和质量平衡方程，得到满足流动与表面材料烧蚀耦合机理的可用于下一步父迭代的固壁面温度和气体组分的新值。
此外，还需要指出的是，表面网格微元上沿参数和固壁面上气体温度、气体各组分质量分数为已知量，可以直接获取；并且，表面网格微元上沿参数，并不限定于气体压强 、温度 和气体各组分质量分数三者，还可以是其它参数形式，如气体密度、热流通量、质量通量、能量通量等，尽管参数可以有很多形式，但无论是哪种参数，最后可通过公开资料提供的变换方式，转化为气体压强P up 、温度T up 和气体各组分质量分数
Figure SMS_11
。
步骤S12：基于表面压力边界条件、烧蚀机制、气体状态方程、气体质量扩散和能力运输机制，并利用所述第一流动特性参数确定出与所述当前迭代步对应的第二流动特性参数；所述第二流动特性参数包括当前气体压强、当前气体各组分质量生成率、当前总质量生成率、当前总能量生成率、当前气体密度、当前气体组分扩散系数和当前热传导系数。
本实施例中，基于表面压力边界条件、烧蚀机制、气体状态方程、气体质量扩散和能力运输机制，以利用第一流动特性参数当前气体温度和当前气体各组分质量分数确定出与当前迭代步对应的第二流动特性参数，包括当前气体压强、当前气体各组分质量生成率、当前总质量生成率、当前总能量生成率、当前气体密度、当前气体组分扩散系数和当前热传导系数。也即，结合表面压力边界条件、烧蚀机制、气体状态方程、气体质量扩散和能量输运机制，计算第n步迭代所需的烧蚀参数和输运系数。上述飞行器具体的表面压力边界条件、烧蚀机制、气体状态方程、气体质量扩散和能量输运机制，其表达形式众多，可以由公开资料或者试验手段获得，本发明均能较好地适用，因此不作特殊限定。
步骤S13：利用所述当前气体各组分质量分数和对应的组分质量扰动量确定出所述当前气体各组分质量生成率关于气体各组分质量分数的第一偏导数，并利用所述当前气体温度和对应的温度扰动量确定出所述当前总能量生成率关于气体温度的第二偏导数。
本实施例中，在当前气体各组分质量分数T n 和当前气体温度
Figure SMS_12
的基础上，引入各自对应的组分质量扰动量/>
Figure SMS_13
和温度扰动量/>
Figure SMS_14
，以分别确定出当前气体各组分质量生成率关于气体各组分质量分数的第一偏导数/>
Figure SMS_15
，以及当前总能量生成率关于气体温度的第二偏导数/>
Figure SMS_16
。需要指出的是，表面材料烧蚀机制相关偏导数的直接解析计算较为困难：一方面由于常见表面材料烧蚀机制很复杂，包括热解、相变、氧化、催化、机械剥落等，其偏导数具体形式很复杂，推导过程相对困难；另一方面不同材料的烧蚀机制存在很大不同，物理模型千差万别，难以写出统一的解析形式的偏导数。而本申请实施例通过引入小扰动实现偏导数的近似计算，不仅不需要推导偏导数具体的数学形式，实现简便，而且可广泛适用于各类材料的多种烧蚀机制。
步骤S14：基于所述第一流动特性参数、所述第二流动特性参数、所述第一偏导数和所述第二偏导数，并按照预设的隐式处理和显示处理相结合的方法对质量平衡方程和能量平衡方程分别进行迭代求解，直到满足预设收敛条件，以得到迭代后的目标气体温度和目标气体各组分质量分数。
本实施例中，基于第一流动特性参数、第二流动特性参数流动特性参数、第一偏导数和第二偏导数，并按照预设的隐式处理和显示处理相结合的方法混合对质量平衡方程和能量平衡方程分别进行迭代求解，直到满足收敛条件时得到迭代后输出的固壁面上的目标气体温度和目标气体各组分质量分数。也即得到了表面微元固壁面上各气体组分质量分数和气体温度的新值：
Figure SMS_17
，/>
Figure SMS_18
，该值符合表面材料烧蚀机理，并与表面网格微元上沿的气体压强P up 、温度T up 和气体各组分质量分数/>
Figure SMS_19
相匹配，满足流动与烧蚀耦合模拟需求。本方案通过将隐式处理和显示处理相结合，混合运用显式和隐式离散处理方案，实现各气体组分质量平衡方程、能量平衡方程相对独立的隐式迭代计算，无需求解耦合矩阵，实现简便，稳定性好且计算效率高。
进一步的，上述按照预设的隐式处理和显示处理相结合的方法对质量平衡方程和能量平衡方程分别进行迭代求解，直到满足预设收敛条件，以得到迭代后的目标气体温度和目标气体各组分质量分数，包括：按照预设的隐式处理和显示处理相结合的方法对质量平衡方程和能量平衡方程分别进行迭代求解，并判断迭代相对残差值是否小于预设残差阈值；若所述迭代相对残差值小于所述预设残差阈值，则输出与下一迭代步对应的目标气体温度和目标气体各组分质量分数；若所述迭代相对残差值大于所述预设残差阈值，则将所述下一迭代步作为当前迭代步，并重新跳转至所述在高超声速飞行器表面烧蚀过程中获取飞行器固壁面上与当前迭代步对应的第一流动特性参数的步骤，直到所述迭代相对残差值小于所述预设残差阈值，以得到迭代后的目标气体温度和目标气体各组分质量分数。可以理解的是，在迭代求解过程中，需判断迭代是否收敛，本实施例中通过判断迭代相对残差值
Figure SMS_20
是否小于预设残差阈值/>
Figure SMS_21
，以确定是否收敛，/>
Figure SMS_22
的表达式如下，表示组分迭代最大相对残差：/>
Figure SMS_23
其中
Figure SMS_24
为气体温度的相对残差，/>
Figure SMS_25
为第s气体组分的迭代相对残差。由于迭代过程中/>
Figure SMS_26
可能接近或等于零，因此在分母中进行限制，防止分母为零。
当迭代相对残差值小于预设残差阈值，即
Figure SMS_27
时，迭代收敛，则输出与下一迭代步对应的目标气体温度和目标气体各组分质量分数；当迭代相对残差值大于预设残差阈值，即/>
Figure SMS_28
时，迭代未收敛，进行下一步迭代，即n=n+1，把下一迭代步作为当前迭代步，并重新跳转至在高超声速飞行器表面烧蚀过程中获取飞行器固壁面上与当前迭代步对应的第一流动特性参数的步骤，直到迭代相对残差值小于预设残差阈值，以得到迭代后的目标气体温度和目标气体各组分质量分数。其中，/>
Figure SMS_29
为精度指标，在流动控制方程求解初期，流场收敛程度低，为减少迭代次数，提升计算效率，/>
Figure SMS_30
取值相对较大，如/>
Figure SMS_31
；在流动控制方程求解后期，流场收敛程度高，为提升整体计算精度，/>
Figure SMS_32
取值相对较小，如
Figure SMS_33
。需要指出的是，判断收敛的方式有很多，本实施例仅以其中一种为例进行介绍。由此可见，本申请实施例设计了收敛精度指标的调节方案，可用流场与表面烧蚀耦合的数值求解，依据流场变化特征，在流场求解的初期和后期，设置不同精度指标，兼顾了计算效率和计算精准度。
可见，本申请通过在高超声速飞行器表面烧蚀过程中获取飞行器固壁面上与当前迭代步对应的第一流动特性参数；所述第一流动特性参数包括当前气体温度和当前气体各组分质量分数；基于表面压力边界条件、烧蚀机制、气体状态方程、气体质量扩散和能力运输机制，并利用所述第一流动特性参数确定出与所述当前迭代步对应的第二流动特性参数；所述第二流动特性参数包括当前气体压强、当前气体各组分质量生成率、当前总质量生成率、当前总能量生成率、当前气体密度、当前气体组分扩散系数和当前热传导系数；利用所述当前气体各组分质量分数和对应的组分质量扰动量确定出所述当前气体各组分质量生成率关于气体各组分质量分数的第一偏导数，并利用所述当前气体温度和对应的温度扰动量确定出所述当前总能量生成率关于气体温度的第二偏导数；基于所述第一流动特性参数、所述第二流动特性参数、所述第一偏导数和所述第二偏导数，并按照预设的隐式处理和显示处理相结合的方法对质量平衡方程和能量平衡方程分别进行迭代求解，直到满足预设收敛条件，以得到迭代后的目标气体温度和目标气体各组分质量分数。由此可见，本申请在高超声速飞行器表面烧蚀过程中首先获取到飞行器固壁面上的流动特性参数，包括与当前迭代步对应的当前气体温度和当前气体各组分质量分数；再结合表面压力边界条件、烧蚀机制、气体状态方程、气体质量扩散和能力运输机制，利用当前气体温度和当前气体各组分质量分数确定出当前迭代步对应的当前气体压强、当前气体各组分质量生成率、当前总质量生成率、当前总能量生成率、当前气体密度、当前气体组分扩散系数和当前热传导系数等其他流动特性参数；然后在当前气体各组分质量分数和当前气体温度的基础上，引入各自对应的组分质量扰动量和温度扰动量，以分别确定出当前气体各组分质量生成率关于气体各组分质量分数的第一偏导数，以及当前总能量生成率关于气体温度的第二偏导数，通过引入小扰动实现偏导数的近似计算，不仅不需要推导偏导数具体的数学形式，实现简便，而且可广泛适用于各类材料的多种烧蚀机制；最后基于各流动特性参数、第一偏导数和第二偏导数，并按照预设的隐式处理和显示处理相结合的方法混合对质量平衡方程和能量平衡方程分别进行迭代求解，直到满足收敛条件时得到迭代后输出的固壁面上的目标气体温度和目标气体各组分质量分数。上述方案将隐式处理和显示处理相结合，混合运用显式和隐式离散处理方案，实现各气体组分质量平衡方程、能量平衡方程相对独立的隐式迭代计算，无需求解耦合矩阵，实现简便，稳定性好且计算效率高。
参见图4所示，本申请实施例公开了一种具体的飞行器表面烧蚀模拟混合处理方法，相对于上一实施例，本实施例对技术方案作了进一步的说明和优化。具体包括：
步骤S21：在高超声速飞行器表面烧蚀过程中获取飞行器固壁面上与当前迭代步对应的第一流动特性参数；所述第一流动特性参数包括当前气体温度和当前气体各组分质量分数。
步骤S22：基于表面压力边界条件、烧蚀机制、气体状态方程、气体质量扩散和能力运输机制，并利用所述第一流动特性参数确定出与所述当前迭代步对应的第二流动特性参数；所述第二流动特性参数包括当前气体压强、当前气体各组分质量生成率、当前总质量生成率、当前总能量生成率、当前气体密度、当前气体组分扩散系数和当前热传导系数。
本实施例中，上述基于表面压力边界条件、烧蚀机制、气体状态方程、气体质量扩散和能力运输机制，并利用所述第一流动特性参数确定出与所述当前迭代步对应的第二流动特性参数，具体可以包括：基于飞行器表面压力边界条件，并利用所述第一流动特性参数确定出与所述当前迭代步对应的当前气体压强；基于飞行器表面材料的烧蚀机制，并利用所述第一流动特性参数和所述当前气体压强确定出与所述当前迭代步对应的当前气体各组分质量生成率、当前总质量生成率和当前总能量生成率；基于气体状态方程，并利用所述第一流动特性参数和所述当前气体压强确定出与所述当前迭代步对应的当前气体密度；基于气体质量扩散和能力运输机制，并利用所述第一流动特性参数、所述当前气体压强和所述当前气体密度确定出与所述当前迭代步对应的当前气体组分扩散系数和当前热传导系数；基于所述当前气体压强、所述当前气体各组分质量生成率、所述当前总质量生成率、所述当前总能量生成率、所述当前气体密度、所述当前气体组分扩散系数和所述当前热传导系数构成第二流动特性参数。也即，基于飞行器表面压力边界条件，由流动界面上气体压强P up 、当前气体温度T n 和当前气体各组分质量分数
Figure SMS_51
，计算得到第n步迭代的当前气体压强
Figure SMS_55
；基于飞行器表面材料的烧蚀机制，由当前气体温度T n 、当前气体各组分质量分数/>
Figure SMS_58
和当前气体压强/>
Figure SMS_36
计算烧蚀产生的第n步迭代的当前气体各组分质量生成率/>
Figure SMS_40
、当前总质量生成率/>
Figure SMS_44
和当前总能量生成率/>
Figure SMS_48
；基于气体状态方程，由当前气体温度T n 、当前气体各组分质量分数/>
Figure SMS_41
和当前气体压强/>
Figure SMS_45
计算第n步迭代的当前气体密度/>
Figure SMS_49
；基于气体质量扩散和能力运输机制，由当前气体温度T n 、当前气体密度/>
Figure SMS_53
、当前气体各组分质量分数/>
Figure SMS_52
和当前气体压强/>
Figure SMS_56
计算第n步迭代的当前气体组分扩散系数/>
Figure SMS_59
和热传导系数/>
Figure SMS_61
。上述的当前气体压强/>
Figure SMS_50
、当前气体各组分质量生成率/>
Figure SMS_54
、当前总质量生成率/>
Figure SMS_57
、当前总能量生成率
Figure SMS_60
、当前气体密度/>
Figure SMS_35
、当前气体组分扩散系数/>
Figure SMS_38
和当前热传导系数/>
Figure SMS_42
共同构成与当前迭代步，即第n步对应的第二流动特性参数。其中/>
Figure SMS_46
为气体各组分质量生成率函数
Figure SMS_37
的值，/>
Figure SMS_39
为总质量生成率函数/>
Figure SMS_43
的值，/>
Figure SMS_47
为总能量生成率/>
Figure SMS_34
的值。
进一步的，上述基于飞行器表面材料的烧蚀机制，并利用所述第一流动特性参数和所述当前气体压强确定出与所述当前迭代步对应的当前气体各组分质量生成率、当前总质量生成率和当前总能量生成率，包括：基于飞行器表面材料的烧蚀机制构建气体各组分质量生成率函数、总质量生成率函数和总能量生成率函数；利用所述第一流动特性参数、所述当前气体压强、所述气体各组分质量生成率函数、所述总质量生成率函数和所述总能量生成率函数确定出与所述当前迭代步对应的当前气体各组分质量生成率、当前总质量生成率和当前总能量生成率。
为了便于理解以上步骤，本实施例采用以下示例或函数进行补充说明，但并不限定于此：
飞行器表面最常用的表面压力边界条件为零压力梯度条件，采用一阶离散，可得：
Figure SMS_62
；飞行器表面烧蚀材料众多，尽管不同材料的烧蚀机制存在很大不同，物理模型差异很大，但均可以用以下函数抽象成统一的形式：
Figure SMS_63
其中
Figure SMS_64
和/>
Figure SMS_65
分别为烧蚀产生的各气体组分质量生成率函数、总质量生成率函数和能量生成率函数；/>
Figure SMS_66
为这些函数的自变量，分别为气体温度、气体各组分质量分数和气体压强。因此，烧蚀产生的第n步迭代气体各组分质量生成率、总质量生成率和能量生成率可写为：
Figure SMS_67
采用多组分气体的状态方程，由
Figure SMS_68
、/>
Figure SMS_69
和/>
Figure SMS_70
计算得到第n步迭代的表面气体密度
Figure SMS_71
:/>
Figure SMS_72
基于气体质量扩散和能量输运机制，由
Figure SMS_73
和/>
Figure SMS_74
，计算第n步迭代的气体组分扩散系数/>
Figure SMS_75
和热传导系数/>
Figure SMS_76
，其计算方法或模型很多，均能较好应用于本发明，这些模型可以从其他公开文献中查到或由试验得到，此处不作赘述。
步骤S23：基于所述当前气体温度和预设偏导数调节系数确定出对应的温度扰动量，并基于所述当前气体各组分质量分数和所述预设偏导数调节系数确定出对应的组分质量扰动量。
本实施例中，基于当前气体温度和预设偏导数调节系数确定出对应的温度扰动量，并基于当前气体各组分质量分数和预设偏导数调节系数确定出对应的组分质量扰动量。需要注意的是，温度扰动量
Figure SMS_77
和组分质量扰动量/>
Figure SMS_78
的取值很关键，既要足够小，以保证偏导数近似计算的合理性，又不能太小，以防止计算机存储数据的截断误差过大。可以按如下方式设置：
Figure SMS_79
其中，f为预设偏导数调节系数，对于单精度数据
Figure SMS_80
，对于双精度数据
Figure SMS_81
。
步骤S24：利用所述当前气体各组分质量分数和所述组分质量扰动量确定出所述当前气体各组分质量生成率关于气体各组分质量分数的第一偏导数，并利用所述当前气体温度和所述温度扰动量确定出所述总能量生成率关于气体温度的第二偏导数。
本实施例中，利用当前气体各组分质量分数
Figure SMS_82
和组分质量扰动量/>
Figure SMS_83
确定出当前气体各组分质量生成率/>
Figure SMS_84
关于气体各组分质量分数的第一偏导数/>
Figure SMS_85
：
Figure SMS_86
利用当前气体温度T n 和温度扰动量
Figure SMS_87
确定出总能量生成率/>
Figure SMS_88
关于气体温度的第二偏导数/>
Figure SMS_89
：
Figure SMS_90
其中，
Figure SMS_91
为函数/>
Figure SMS_92
的缩写，
Figure SMS_93
为函数/>
Figure SMS_94
的缩写。
步骤S25：确定出包括输入气体各组分质量分数、所述第一偏导数和表面网格微元法向厚度的质量平衡方程，按照预设的隐式处理方法对所述质量平衡方程中的气体各组分质量分数和气体各组分质量生成率进行处理，并按照预设的显示处理方法利用所述第一流动特性参数和所述第二流动特性参数对所述质量平衡方程中其余参数进行处理，以对所述质量平衡方程进行迭代求解。
本实施例中，混合离散求解各气体组分质量平衡方程：对于第s气体组分的质量平衡方程，隐式处理第s气体组分质量分数
Figure SMS_95
及其质量生成率/>
Figure SMS_96
，显式处理其余参数，结合一阶Taylor展开构建迭代关系，得到第n+1步迭代第s气体组分质量分数，即目标气体各组分质量分数/>
Figure SMS_97
：
Figure SMS_98
具体推导如下：首先确定出包括输入气体各组分质量分数、第一偏导数和表面网格微元法向厚度的质量平衡方程。可以理解的是，对于第s气体组分，在表面“无限薄”的微元内的满足以下质量平衡方程：
Figure SMS_99
其中
Figure SMS_100
为第s气体组分质量分数；I s 为第s气体组分质量分数的表面法向梯度，由于表面网格微元法向厚度h足够小，因此可由空间一阶离散得到/>
Figure SMS_101
；V为烧蚀引射法向速度，可由总质量生成率计算得到/>
Figure SMS_102
，因此上述质量平衡方程可写为：
Figure SMS_103
按照预设的隐式处理方法对质量平衡方程中的气体各组分质量分数和气体各组分质量生成率进行处理，并按照预设的显示处理方法利用第一流动特性参数和第二流动特性参数对质量平衡方程中其余参数进行处理，以对质量平衡方程进行迭代求解。具体的，其中第s气体组分质量分数
Figure SMS_104
及气体各组分质量生成率/>
Figure SMS_105
采用隐式处理，即采用/>
Figure SMS_106
代替方程中的/>
Figure SMS_107
和/>
Figure SMS_108
；方程中的其它量，采用显式处理近似，即直接采用第n步的值。这里
Figure SMS_109
为函数/>
Figure SMS_110
的缩写。因此上述质量平衡方程可写为：
Figure SMS_111
将
Figure SMS_112
在/>
Figure SMS_113
处一阶Taylor展开得/>
Figure SMS_114
，因此上述质量平衡方程可写为：
Figure SMS_115
整理上式得到迭代计算关系：
Figure SMS_116
该步骤中将隐式处理和显示处理相结合，吸收两者优点，避免了两者缺点，实现了各组分方程相对独立的隐式迭代计算，实现简便、迭代稳定、效率高：（1）对于第s气体组分，其烧蚀质量生成率是该组份变化最直接的主导因素，因此对该组分质量分数及其烧蚀质量生成率进行隐式处理，能有效提升迭代稳定性和效率；（2）
Figure SMS_117
等参量虽然受单个气体组分质量分数的影响，但其影响幅度相对较小，因此采用显示处理，能大幅简化其复杂耦合关系的求解，有效避免耦合矩阵及其求逆运算的复杂实现过程。
步骤S26：确定出包括输入气体温度、所述输入气体各组分质量分数、所述第二偏导数、所述表面网格微元法向厚度和气体各组分生成焓的能量平衡方程，按照所述隐式处理方法对所述能量平衡方程中的气体温度、气体各组分质量分数和总能量生成率进行处理，并按照所述显示处理方法利用所述第一流动特性参数和所述第二流动特性参数对所述能量平衡方程中其余参数进行处理，以对所述能量平衡方程进行迭代求解，直到满足预设收敛条件，以得到迭代后的目标气体温度和目标气体各组分质量分数。
本实施例中，混合离散求解表面能量平衡方程，方程中气体温度T、气体各组分质量分数
Figure SMS_118
和总能量生成率/>
Figure SMS_119
采用隐式处理，其它量采用显式处理，结合一阶Taylor展开构建迭代关系，计算得到第n+1步迭代表面气体温度，即目标气体温度T n+1 ：
Figure SMS_120
具体推导如下：首先确定出包括输入气体温度、输入气体各组分质量分数、第二偏导数、表面网格微元法向厚度和气体各组分生成焓的能量平衡方程。可以理解的是，在表面“无限薄”的微元内，满足能量平衡方程：
Figure SMS_121
其中h s 为第s组分气体的生成焓，
Figure SMS_122
为表面向外的辐射热流，/>
Figure SMS_123
为表面材料黑体辐射系数，/>
Figure SMS_124
为Stefan-Boltzman常数；/>
Figure SMS_125
为气体温度表面法向梯度，由于h足够小，因此可由空间一阶离散得到/>
Figure SMS_126
；/>
Figure SMS_127
为质量扩散热流，可由质量平衡方程变形得到/>
Figure SMS_128
。因此上述能量平衡方程可写为：
Figure SMS_129
按照所述隐式处理方法对该能量平衡方程中的气体温度、气体各组分质量分数和总能量生成率进行处理，并按照显示处理方法利用第一流动特性参数和第二流动特性参数对能量平衡方程中其余参数进行处理，以对能量平衡方程进行迭代求解。具体的，上述方程中T、
Figure SMS_131
和/>
Figure SMS_134
采用隐式处理，即采用T n+1、/>
Figure SMS_137
和/>
Figure SMS_132
代替T、/>
Figure SMS_135
和/>
Figure SMS_138
；对于/>
Figure SMS_140
，采用与组分质量平衡方程相同的离散处理方法，即用/>
Figure SMS_130
代替/>
Figure SMS_133
，保持能量平衡方程和组分质量平衡方程相容，也即在能量求解过程中，质量扩散热流采用质量平衡方程换算，保证了能量平衡方程和质量平衡方程的相容，避免了二者相对独立求解可能造成的数值振荡。上述能量方程中的其它参数采用显示处理，即直接采用第n步值，其中
Figure SMS_136
为函数/>
Figure SMS_139
的缩写。由此可得：/>
Figure SMS_141
其中
Figure SMS_142
和/>
Figure SMS_143
在T n 处一阶Taylor展开得：
Figure SMS_144
因此能量平衡方程可写为：
Figure SMS_145
整理得到迭代计算关系：
Figure SMS_146
该步骤将隐式处理和显示处理相结合，吸收两者优点，避免了两者缺点，实现简便、迭代稳定、效率高：（1）对于表面能量平衡过程，表面温度变化是最直接的表征，其烧蚀能量产生率、向外辐射热流是表面温度变化的主导因素，因此对表面温度及其烧蚀质量生成率、辐射热流进行隐式处理（辐射热流的隐式由温度隐式处理实现），能有效提升迭代稳定性和效率；（2）能量方程中的质量扩散热流
Figure SMS_147
，采用组分质量平衡方程进行换算，并采用与组分质量平衡方程相同的离散处理方法，既保证了能量平衡方程和质量平衡方程的相容，避免了二者相对独立求解可能造成的数值振荡，又避免了对于壁面/>
Figure SMS_148
的隐式处理难题：对于高超声速流动，表面气体密度受表面温度影响较为显著，采用显式处理会影响计算稳定性，而采用隐式处理则形式相对复杂；（3）其它/>
Figure SMS_149
等参量受温度影响，但其影响幅度相对较小，采用显示处理，能大幅简化其复杂耦合关系的求解，有效避免耦合矩阵及其求逆运算的复杂实现过程。
其中，关于上述步骤S21更加具体的处理过程可以参考前述实施例中公开的相应内容，在此不再进行赘述。
可见，本申请实施例从表面微元烧蚀机制及其质量和能量守恒的基本原理出发，构建表面烧蚀迭代模拟体系，依据各物理量在烧蚀体系中影响程度的差异，混合运用显式和隐式离散处理方案，实现各气体组分质量平衡方程、能量平衡方程相对独立的隐式迭代计算，无需求解耦合矩阵，易于实现、稳定性好且计算效率高。并且，在能量求解过程中，质量扩散热流采用质量平衡方程换算，保证了能量平衡方程和质量平衡方程的相容，避免了二者相对独立求解可能造成的数值振荡。
下面以球头碳-碳烧蚀模拟工况为例，对本发明的应用效果进行详细说明。
计算工况：采用球头模型，半径1.0m。计算飞行高度65km，来流速度8km/s。数值模拟考虑表面烧蚀的热化学非平衡流动。球头表面为碳-碳材料，烧蚀过程考虑碳的升华、氧化以及碳对氧原子的催化作用。采用本发明实现烧蚀过程中表面微元的气体组分质量平衡方程和能量平衡方程的求解，并将本发明的计算结果与非平衡流动求解相耦合。
图5为采用本发明计算得到的驻点线主要气体组分质量分数分布。图6为文献计算得到的驻点线主要气体组分质量分数分布。可以看出，本发明计算结果与文献结果（文献“Keenan J A, Candler G V. Simulation of Ablation in Earth Atmospheric Entry,AIAA-93-2789. ”的图6）符合较好，这说明本发明能较好地实现烧蚀过程各气体组分质量平衡方程的准确求解。
图7为采用本发明计算得到的表面热流分布。图中“Present”为本发明计算结果，“Δ”为文献结果（取值来源于文献 “Keenan J A, Candler G V. Simulation ofAblation in Earth Atmospheric Entry, AIAA-93-2789.”的图12）。可以看出，本发明计算结果与文献结果符合较好，这说明本发明能较好地实现烧蚀过程能量平衡方程的准确求解。
参见图8所示，本申请实施例公开了一种飞行器表面烧蚀模拟混合处理方法，该方法包括：
第一参数获取模块11，用于在高超声速飞行器表面烧蚀过程中获取飞行器固壁面上与当前迭代步对应的第一流动特性参数；所述第一流动特性参数包括当前气体温度和当前气体各组分质量分数；
第二参数获取模块12，用于基于表面压力边界条件、烧蚀机制、气体状态方程、气体质量扩散和能力运输机制，并利用所述第一流动特性参数确定出与所述当前迭代步对应的第二流动特性参数；所述第二流动特性参数包括当前气体压强、当前气体各组分质量生成率、当前总质量生成率、当前总能量生成率、当前气体密度、当前气体组分扩散系数和当前热传导系数；
偏导数确定模块13，用于利用所述当前气体各组分质量分数和对应的组分质量扰动量确定出所述当前气体各组分质量生成率关于气体各组分质量分数的第一偏导数，并利用所述当前气体温度和对应的温度扰动量确定出所述当前总能量生成率关于气体温度的第二偏导数；
迭代求解模块14，用于基于所述第一流动特性参数、所述第二流动特性参数、所述第一偏导数和所述第二偏导数，并按照预设的隐式处理和显示处理相结合的方法对质量平衡方程和能量平衡方程分别进行迭代求解，直到满足预设收敛条件，以得到迭代后的目标气体温度和目标气体各组分质量分数。
可见，本申请通过在高超声速飞行器表面烧蚀过程中获取飞行器固壁面上与当前迭代步对应的第一流动特性参数；所述第一流动特性参数包括当前气体温度和当前气体各组分质量分数；基于表面压力边界条件、烧蚀机制、气体状态方程、气体质量扩散和能力运输机制，并利用所述第一流动特性参数确定出与所述当前迭代步对应的第二流动特性参数；所述第二流动特性参数包括当前气体压强、当前气体各组分质量生成率、当前总质量生成率、当前总能量生成率、当前气体密度、当前气体组分扩散系数和当前热传导系数；利用所述当前气体各组分质量分数和对应的组分质量扰动量确定出所述当前气体各组分质量生成率关于气体各组分质量分数的第一偏导数，并利用所述当前气体温度和对应的温度扰动量确定出所述当前总能量生成率关于气体温度的第二偏导数；基于所述第一流动特性参数、所述第二流动特性参数、所述第一偏导数和所述第二偏导数，并按照预设的隐式处理和显示处理相结合的方法对质量平衡方程和能量平衡方程分别进行迭代求解，直到满足预设收敛条件，以得到迭代后的目标气体温度和目标气体各组分质量分数。由此可见，本申请在高超声速飞行器表面烧蚀过程中首先获取到飞行器固壁面上的流动特性参数，包括与当前迭代步对应的当前气体温度和当前气体各组分质量分数；再结合表面压力边界条件、烧蚀机制、气体状态方程、气体质量扩散和能力运输机制，利用当前气体温度和当前气体各组分质量分数确定出当前迭代步对应的当前气体压强、当前气体各组分质量生成率、当前总质量生成率、当前总能量生成率、当前气体密度、当前气体组分扩散系数和当前热传导系数等其他流动特性参数；然后在当前气体各组分质量分数和当前气体温度的基础上，引入各自对应的组分质量扰动量和温度扰动量，以分别确定出当前气体各组分质量生成率关于气体各组分质量分数的第一偏导数，以及当前总能量生成率关于气体温度的第二偏导数，通过引入小扰动实现偏导数的近似计算，不仅不需要推导偏导数具体的数学形式，实现简便，而且可广泛适用于各类材料的多种烧蚀机制；最后基于各流动特性参数、第一偏导数和第二偏导数，并按照预设的隐式处理和显示处理相结合的方法混合对质量平衡方程和能量平衡方程分别进行迭代求解，直到满足收敛条件时得到迭代后输出的固壁面上的目标气体温度和目标气体各组分质量分数。上述方案将隐式处理和显示处理相结合，混合运用显式和隐式离散处理方案，实现各气体组分质量平衡方程、能量平衡方程相对独立的隐式迭代计算，无需求解耦合矩阵，实现简便，稳定性好且计算效率高。
图9为本申请实施例提供的一种电子设备的结构示意图。具体可以包括：至少一个处理器21、至少一个存储器22、电源23、通信接口24、输入输出接口25和通信总线26。其中，所述存储器22用于存储计算机程序，所述计算机程序由所述处理器21加载并执行，以实现前述任一实施例公开的由电子设备执行的飞行器表面烧蚀模拟混合处理方法中的相关步骤。
本实施例中，电源23用于为电子设备20上的各硬件设备提供工作电压；通信接口24能够为电子设备20创建与外界设备之间的数据传输通道，其所遵循的通信协议是能够适用于本申请技术方案的任意通信协议，在此不对其进行具体限定；输入输出接口25，用于获取外界输入数据或向外界输出数据，其具体的接口类型可以根据具体应用需要进行选取，在此不进行具体限定。
其中，处理器21可以包括一个或多个处理核心，比如4核心处理器、8核心处理器等。处理器21可以采用DSP(Digital Signal Processing，数字信号处理)、FPGA(Field－Programmable Gate Array，现场可编程门阵列)、PLA(Programmable Logic Array，可编程逻辑阵列)中的至少一种硬件形式来实现。处理器21也可以包括主处理器和协处理器，主处理器是用于对在唤醒状态下的数据进行处理的处理器，也称CPU(Central ProcessingUnit，中央处理器)；协处理器是用于对在待机状态下的数据进行处理的低功耗处理器。在一些实施例中，处理器21可以在集成有GPU(Graphics Processing Unit，图像处理器)，GPU用于负责显示屏所需要显示的内容的渲染和绘制。一些实施例中，处理器21还可以包括AI(Artificial Intelligence，人工智能)处理器，该AI处理器用于处理有关机器学习的计算操作。
另外，存储器22作为资源存储的载体，可以是只读存储器、随机存储器、磁盘或者光盘等，其上所存储的资源包括操作系统221、计算机程序222及数据223等，存储方式可以是短暂存储或者永久存储。
其中，操作系统221用于管理与控制电子设备20上的各硬件设备以及计算机程序222，以实现处理器21对存储器22中海量数据223的运算与处理，其可以是Windows、Unix、Linux等。计算机程序222除了包括能够用于完成前述任一实施例公开的由电子设备20执行的飞行器表面烧蚀模拟混合处理方法的计算机程序之外，还可以进一步包括能够用于完成其他特定工作的计算机程序。数据223除了可以包括电子设备接收到的由外部设备传输进来的数据，也可以包括由自身输入输出接口25采集到的数据等。
进一步的，本申请实施例还公开了一种计算机可读存储介质，所述存储介质中存储有计算机程序，所述计算机程序被处理器加载并执行时，实现前述任一实施例公开的由飞行器表面烧蚀模拟混合处理过程中执行的方法步骤。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其它实施例的不同之处，各个实施例之间相同或相似部分互相参见即可。对于实施例公开的装置而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
专业人员还可以进一步意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本申请的范围。
结合本文中所公开的实施例描述的方法或算法的步骤可以直接用硬件、处理器执行的软件模块，或者二者的结合来实施。软件模块可以置于随机存储器（RAM）、内存、只读存储器（ROM）、电可编程ROM、电可擦除可编程ROM、寄存器、硬盘、可移动磁盘、CD-ROM、或技术领域内所公知的任意其它形式的存储介质中。
最后，还需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。
以上对本发明所提供的一种飞行器表面烧蚀模拟混合处理方法、装置、设备及存储介质进行了详细介绍，本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种飞行器表面烧蚀模拟混合处理方法，其特征在于，包括：
在高超声速飞行器表面烧蚀过程中获取飞行器固壁面上与当前迭代步对应的第一流动特性参数；所述第一流动特性参数包括当前气体温度和当前气体各组分质量分数；
基于表面压力边界条件、烧蚀机制、气体状态方程、气体质量扩散和能力运输机制，并利用所述第一流动特性参数确定出与所述当前迭代步对应的第二流动特性参数；所述第二流动特性参数包括当前气体压强、当前气体各组分质量生成率、当前总质量生成率、当前总能量生成率、当前气体密度、当前气体组分扩散系数和当前热传导系数；
利用所述当前气体各组分质量分数和对应的组分质量扰动量确定出所述当前气体各组分质量生成率关于气体各组分质量分数的第一偏导数，并利用所述当前气体温度和对应的温度扰动量确定出所述当前总能量生成率关于气体温度的第二偏导数；
基于所述第一流动特性参数、所述第二流动特性参数、所述第一偏导数和所述第二偏导数，并按照预设的隐式处理和显示处理相结合的方法对质量平衡方程和能量平衡方程分别进行迭代求解，直到满足预设收敛条件，以得到迭代后的目标气体温度和目标气体各组分质量分数。
2.根据权利要求1所述的飞行器表面烧蚀模拟混合处理方法，其特征在于，所述基于表面压力边界条件、烧蚀机制、气体状态方程、气体质量扩散和能力运输机制，并利用所述第一流动特性参数确定出与所述当前迭代步对应的第二流动特性参数，包括：
基于飞行器表面压力边界条件，并利用所述第一流动特性参数确定出与所述当前迭代步对应的当前气体压强；
基于飞行器表面材料的烧蚀机制，并利用所述第一流动特性参数和所述当前气体压强确定出与所述当前迭代步对应的当前气体各组分质量生成率、当前总质量生成率和当前总能量生成率；
基于气体状态方程，并利用所述第一流动特性参数和所述当前气体压强确定出与所述当前迭代步对应的当前气体密度；
基于气体质量扩散和能力运输机制，并利用所述第一流动特性参数、所述当前气体压强和所述当前气体密度确定出与所述当前迭代步对应的当前气体组分扩散系数和当前热传导系数；
基于所述当前气体压强、所述当前气体各组分质量生成率、所述当前总质量生成率、所述当前总能量生成率、所述当前气体密度、所述当前气体组分扩散系数和所述当前热传导系数构成第二流动特性参数。
3.根据权利要求2所述的飞行器表面烧蚀模拟混合处理方法，其特征在于，所述基于飞行器表面材料的烧蚀机制，并利用所述第一流动特性参数和所述当前气体压强确定出与所述当前迭代步对应的当前气体各组分质量生成率、当前总质量生成率和当前总能量生成率，包括：
基于飞行器表面材料的烧蚀机制构建气体各组分质量生成率函数、总质量生成率函数和总能量生成率函数；
利用所述第一流动特性参数、所述当前气体压强、所述气体各组分质量生成率函数、所述总质量生成率函数和所述总能量生成率函数确定出与所述当前迭代步对应的当前气体各组分质量生成率、当前总质量生成率和当前总能量生成率。
4.根据权利要求1所述的飞行器表面烧蚀模拟混合处理方法，其特征在于，所述利用所述当前气体各组分质量分数和对应的组分质量扰动量确定出所述当前气体各组分质量生成率关于气体各组分质量分数的第一偏导数，并利用所述当前气体温度和对应的温度扰动量确定出所述当前总能量生成率关于气体温度的第二偏导数，包括：
基于所述当前气体温度和预设偏导数调节系数确定出对应的温度扰动量，并基于所述当前气体各组分质量分数和所述预设偏导数调节系数确定出对应的组分质量扰动量；
利用所述当前气体各组分质量分数和所述组分质量扰动量确定出所述当前气体各组分质量生成率关于气体各组分质量分数的第一偏导数，并利用所述当前气体温度和所述温度扰动量确定出所述总能量生成率关于气体温度的第二偏导数。
5.根据权利要求1所述的飞行器表面烧蚀模拟混合处理方法，其特征在于，所述在高超声速飞行器表面烧蚀过程中获取飞行器固壁面上与当前迭代步对应的第一流动特性参数的过程中，还包括：
获取流动界面上的输入气体温度和输入气体各组分质量分数。
6.根据权利要求5所述的飞行器表面烧蚀模拟混合处理方法，其特征在于，所述基于所述第一流动特性参数、所述第二流动特性参数、所述第一偏导数和所述第二偏导数，并按照预设的隐式处理和显示处理相结合的方法对质量平衡方程和能量平衡方程分别进行迭代求解，包括：
确定出包括所述输入气体各组分质量分数、所述第一偏导数和表面网格微元法向厚度的质量平衡方程；
按照预设的隐式处理方法对所述质量平衡方程中的气体各组分质量分数和气体各组分质量生成率进行处理，并按照预设的显示处理方法利用所述第一流动特性参数和所述第二流动特性参数对所述质量平衡方程中其余参数进行处理，以对所述质量平衡方程进行迭代求解；
确定出包括所述输入气体温度、所述输入气体各组分质量分数、所述第二偏导数、所述表面网格微元法向厚度和气体各组分生成焓的能量平衡方程；
按照所述隐式处理方法对所述能量平衡方程中的气体温度、气体各组分质量分数和总能量生成率进行处理，并按照所述显示处理方法利用所述第一流动特性参数和所述第二流动特性参数对所述能量平衡方程中其余参数进行处理，以对所述能量平衡方程进行迭代求解。
7.根据权利要求1至6任一项所述的飞行器表面烧蚀模拟混合处理方法，其特征在于，所述按照预设的隐式处理和显示处理相结合的方法对质量平衡方程和能量平衡方程分别进行迭代求解，直到满足预设收敛条件，以得到迭代后的目标气体温度和目标气体各组分质量分数，包括：
按照预设的隐式处理和显示处理相结合的方法对质量平衡方程和能量平衡方程分别进行迭代求解，并判断迭代相对残差值是否小于预设残差阈值；
若所述迭代相对残差值小于所述预设残差阈值，则输出与下一迭代步对应的目标气体温度和目标气体各组分质量分数；
若所述迭代相对残差值大于所述预设残差阈值，则将所述下一迭代步作为当前迭代步，并重新跳转至所述在高超声速飞行器表面烧蚀过程中获取飞行器固壁面上与当前迭代步对应的第一流动特性参数的步骤，直到所述迭代相对残差值小于所述预设残差阈值，以得到迭代后的目标气体温度和目标气体各组分质量分数。
8.一种飞行器表面烧蚀模拟混合处理装置，其特征在于，包括：
第一参数获取模块，用于在高超声速飞行器表面烧蚀过程中获取飞行器固壁面上与当前迭代步对应的第一流动特性参数；所述第一流动特性参数包括当前气体温度和当前气体各组分质量分数；
第二参数获取模块，用于基于表面压力边界条件、烧蚀机制、气体状态方程、气体质量扩散和能力运输机制，并利用所述第一流动特性参数确定出与所述当前迭代步对应的第二流动特性参数；所述第二流动特性参数包括当前气体压强、当前气体各组分质量生成率、当前总质量生成率、当前总能量生成率、当前气体密度、当前气体组分扩散系数和当前热传导系数；
偏导数确定模块，用于利用所述当前气体各组分质量分数和对应的组分质量扰动量确定出所述当前气体各组分质量生成率关于气体各组分质量分数的第一偏导数，并利用所述当前气体温度和对应的温度扰动量确定出所述当前总能量生成率关于气体温度的第二偏导数；
迭代求解模块，用于基于所述第一流动特性参数、所述第二流动特性参数、所述第一偏导数和所述第二偏导数，并按照预设的隐式处理和显示处理相结合的方法对质量平衡方程和能量平衡方程分别进行迭代求解，直到满足预设收敛条件，以得到迭代后的目标气体温度和目标气体各组分质量分数。
9.一种电子设备，其特征在于，包括：
存储器，用于保存计算机程序；
处理器，用于执行所述计算机程序，以实现如权利要求1至7任一项所述的飞行器表面烧蚀模拟混合处理方法的步骤。
10.一种计算机可读存储介质，其特征在于，用于存储计算机程序；其中，所述计算机程序被处理器执行时实现如权利要求1至7任一项所述的飞行器表面烧蚀模拟混合处理方法的步骤。
CN202310319293.9A 2023-03-29 2023-03-29 飞行器表面烧蚀模拟混合处理方法、装置、设备及介质 Active CN116049997B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310319293.9A CN116049997B (zh) 2023-03-29 2023-03-29 飞行器表面烧蚀模拟混合处理方法、装置、设备及介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310319293.9A CN116049997B (zh) 2023-03-29 2023-03-29 飞行器表面烧蚀模拟混合处理方法、装置、设备及介质

Publications (2)

Publication Number Publication Date
CN116049997A CN116049997A (zh) 2023-05-02
CN116049997B true CN116049997B (zh) 2023-06-02

Family

ID=86124073

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310319293.9A Active CN116049997B (zh) 2023-03-29 2023-03-29 飞行器表面烧蚀模拟混合处理方法、装置、设备及介质

Country Status (1)

Country Link
CN (1) CN116049997B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117494400B (zh) * 2023-10-18 2024-04-19 中国空气动力研究与发展中心计算空气动力研究所 强对流条件下低温平板干模态结霜相似因素获得方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106844440A (zh) * 2016-12-15 2017-06-13 中国航天空气动力技术研究院 多组分防热材料热化学平衡烧蚀数据库构建方法
CN113065201A (zh) * 2021-05-08 2021-07-02 中国空气动力研究与发展中心计算空气动力研究所 一种考虑滑移修正的辐射平衡温度计算方法
CN114818149A (zh) * 2022-06-28 2022-07-29 中国空气动力研究与发展中心计算空气动力研究所 飞行器表面壁面热流分解方法、分析方法、介质及设备
CN115587551A (zh) * 2022-12-12 2023-01-10 北京航空航天大学 一种高超声速飞行器防热结构烧蚀行为多尺度预测方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6289977B2 (ja) * 2014-03-31 2018-03-07 三菱重工業株式会社 飛しょう体、及び、飛しょう体の動作方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106844440A (zh) * 2016-12-15 2017-06-13 中国航天空气动力技术研究院 多组分防热材料热化学平衡烧蚀数据库构建方法
CN113065201A (zh) * 2021-05-08 2021-07-02 中国空气动力研究与发展中心计算空气动力研究所 一种考虑滑移修正的辐射平衡温度计算方法
CN114818149A (zh) * 2022-06-28 2022-07-29 中国空气动力研究与发展中心计算空气动力研究所 飞行器表面壁面热流分解方法、分析方法、介质及设备
CN115587551A (zh) * 2022-12-12 2023-01-10 北京航空航天大学 一种高超声速飞行器防热结构烧蚀行为多尺度预测方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Supersonic turbulent channel flows over spanwise-oriented;JianQiang Chen等;《Physics of Fluids》;第34卷(第1期);016109 *
航天飞行器防热涂层烧蚀行为及机理研究;严旭;王洪波;范新中;杨东生;高超;贺晨;王毅飞;;航天器环境工程(第04期);369-376 *
高超声速再入体表面热解烧蚀效应数值模拟;丁明松等;《气体物理》;第8卷(第1期);58-67 *

Also Published As

Publication number Publication date
CN116049997A (zh) 2023-05-02

Similar Documents

Publication Publication Date Title
Zhu et al. Machine learning methods for turbulence modeling in subsonic flows around airfoils
Martins et al. A coupled-adjoint sensitivity analysis method for high-fidelity aero-structural design
Moura et al. Adaptive partial differential equation observer for battery state-of-charge/state-of-health estimation via an electrochemical model
Ceze et al. Anisotropic hp-adaptation framework for functional prediction
CN113051831B (zh) 机床热误差自学习预测模型建模方法及热误差控制方法
Miura An upwind-biased conservative advection scheme for spherical hexagonal–pentagonal grids
CN116049997B (zh) 飞行器表面烧蚀模拟混合处理方法、装置、设备及介质
Eyi et al. Aerothermodynamic design optimization of hypersonic vehicles
Turkel Improving the accuracy of central difference schemes
Biava et al. Fully implicit discrete-adjoint methods for rotorcraft applications
Ceze et al. Drag prediction using adaptive discontinuous finite elements
CN116738891B (zh) 一种增强飞行器流场模拟稳定性的lu-sgs改进方法
Rath et al. Symplectic Gaussian process regression of maps in Hamiltonian systems
Zhou et al. Reduced order unsteady aerodynamic model of a rigid aerofoil in gust encounters
Rao et al. A time-parallel approach to strong-constraint four-dimensional variational data assimilation
Liu et al. Three-dimensional high-order least square-based finite difference-finite volume method on unstructured grids
Xia et al. Particle swarm optimization of aerodynamic shapes with nonuniform shape parameter–based radial basis function
Xu et al. High order accurate and low dissipation method for unsteady compressible viscous flow computation on helicopter rotor in forward flight
Wang et al. A generalized minimal residual method-based immersed boundary-lattice Boltzmann flux solver coupled with finite element method for non-linear fluid-structure interaction problems
Yang et al. Improved automatic kernel construction for Gaussian process regression in small sample learning for predicting lift body aerodynamic performance
Yang et al. High-order gas kinetic flux solver for simulation of two dimensional incompressible flows
Li et al. Generalized linear‐constrained optimal power flow for distribution networks
Cai et al. Optimal guidance for hypersonic reentry using inversion and receding horizon control
García-Gutiérrez et al. Stochastic optimization of high-altitude airship envelopes based on kriging method
CN116341421B (zh) 高超声速流场数值模拟方法、系统、电子设备及存储介质

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