CN112765736A - 一种高超声速钝前缘绕流湍动能入口边界设置方法 - Google Patents

一种高超声速钝前缘绕流湍动能入口边界设置方法 Download PDF

Info

Publication number
CN112765736A
CN112765736A CN202110387351.2A CN202110387351A CN112765736A CN 112765736 A CN112765736 A CN 112765736A CN 202110387351 A CN202110387351 A CN 202110387351A CN 112765736 A CN112765736 A CN 112765736A
Authority
CN
China
Prior art keywords
turbulence
kinetic energy
disturbance
flow field
flow
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
CN202110387351.2A
Other languages
English (en)
Other versions
CN112765736B (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 CN202110387351.2A priority Critical patent/CN112765736B/zh
Publication of CN112765736A publication Critical patent/CN112765736A/zh
Application granted granted Critical
Publication of CN112765736B publication Critical patent/CN112765736B/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
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/10—Noise analysis or noise 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)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Computational Mathematics (AREA)
  • Aviation & Aerospace Engineering (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)
  • Aerodynamic Tests, Hydrodynamic Tests, Wind Tunnels, And Water Tanks (AREA)

Abstract

本发明公开了一种高超声速钝前缘绕流湍动能入口边界设置方法，包括如下步骤：S1，计算基础流场，获取定常基本流场变量
Figure DEST_PATH_IMAGE001
；S2，利用数值模拟系统进行扰动波直接模拟，获取非定常流场变量ϕ；S3，分析流场得到扰动场ϕ'，并获取湍动能分布特征；S4，设置湍动能入口边界；S5，计算湍流/转捩等；本发明不需对湍流模型中湍动能方程进行修改，通过入口边界条件的设置即可避免湍动能方程在激波下游的计算偏差，基于直接数值模拟结果设置更加准确的湍动能入口边界条件，具有计算结果合理、可实现性高的优势，同时可靠性高，可推广至三维流动情况等。

Description

一种高超声速钝前缘绕流湍动能入口边界设置方法
技术领域
本发明涉及钝前缘高超声速流动中边界层内初始扰动幅值估算领域，更为具体的，涉及一种高超声速钝前缘绕流湍动能入口边界设置方法。
背景技术
湍流/转捩是自然界广泛存在的物理现象，在航空/航天领域有着重要的科学研究意义和工程应用价值。在飞行器研制过程中，需要准确计算边界层湍流/转捩发生的位置、范围，因为相比层流状态，湍流/转捩会对边界层带来较大影响，主要表现壁面摩阻和热流，它们会对飞行器的气动特性，例如升阻比、航程、稳定性、操控特性，以及热防护和载荷效能等造成影响。目前关于湍流/转捩的数值计算方法主要有直接数值模拟（DNS）、大涡模拟（LES）和雷诺平均方程（RANS）等三种方法，前两种由于对计算网格要求较高，通常仅作流动机理方面研究，基于RANS方程的方法计算量小，可快速获得结果，具有较高工程实用价值。
目前工程实践中常用的湍流RANS方法包括：S-A模型、
Figure 249991DEST_PATH_IMAGE001
模型、
Figure 596658DEST_PATH_IMAGE002
模型、
Figure 730224DEST_PATH_IMAGE003
模 型、
Figure 430327DEST_PATH_IMAGE004
SST模型、KDO模型等。而常用的转捩RANS方法包括：
Figure 346330DEST_PATH_IMAGE005
模型、
Figure 254112DEST_PATH_IMAGE006
模型、
Figure 622777DEST_PATH_IMAGE007
模型等。这些湍流/转捩计算方法中，绝大多数都需要k方程，这是一个湍动能 方程，表征流场中湍流能量，在自由流中表示湍流强度。在流场数值模拟过程中，首先需要 对自由流边界（来流边界）赋值，然后通过RANS方法中的关于k的输运方程计算k在流场中的 分布，关键是在边界层中的值。那么在流场边界条件幅值时，就涉及k的设置方法。
目前在RANS方法计算中，在设置k的来流边界时，一般将风洞流场中自由流湍流强度Tu ∞%测量结果（例如脉动速度），根据湍流动能计算公式（式1），设为来流湍动能，或者按照基于脉动压力计算的湍流强度等价于基于速度的湍流强度的原则来设置自由流湍动能（高超声速流动中速度脉动较难测量，通常采用压力脉动表征湍流强度，也成为噪声水平，这里统一称为湍流强度）。例如，在对常规高超噪声风洞中流场的模拟中，湍流度约为Tu ∞=1~3%，而静音风洞或真实飞行环境中湍流度Tu ∞<0.1%。（下标“∞”表示自由流）
Figure 860991DEST_PATH_IMAGE008
（1）
对于流场中没有激波的低速流动（马赫数小于1），可以根据风洞湍流强度测量结果直接给来流边界设置k值，因为在自由流区湍流度变化不大。但是，在包含激波的高超声速流场中，特别是在钝前缘弓形激波情况下，强激波会显著改变自由流湍流强度，而这一过程不能被现有RANS方程中的k方程准确描述，由此会对边界层中初始湍流动能的评估造成较大偏差，进而影响湍流/转捩位置的准确计算。图1给出了Mach 6钝平板高超声速流场中RANS方程计算的湍流动能分布和DNS计算单频扰动波结果的对比，云图取值范围以自由流湍动能的20倍为界限。可以看到，两种方法计算结果相差较大，与DNS相比，RANS方程过高估计了激波后大部分流场中的湍动能，而低估了边界层中的湍动能，不能正确的反映对边界层湍流/转捩计算有影响的靠近壁面附近的湍动能分布。因此，在高超声速湍流/转捩的RANS数值模拟中，采用自由流湍流度设置来流边界条件的方法不再适用，不能正确估计激波后自由流湍动能，需要新的入口边界设置技术。
发明内容
本发明的目的在于克服现有技术的不足，提供一种高超声速钝前缘绕流湍动能入口边界设置方法，不需对湍流模型中湍动能方程进行修改，通过入口边界条件的设置即可避免湍动能方程在激波下游的计算偏差，基于直接数值模拟结果设置更加准确的湍动能入口边界条件，具有计算结果合理、可实现性高的优势，同时可靠性高，可推广至三维流动情况等。
本发明的目的是通过以下方案实现的：
一种高超声速钝前缘绕流湍动能入口边界设置方法，包括如下步骤：
S1，计算基础流场，获取定常基本流场变量
Figure 146348DEST_PATH_IMAGE009
；
S2，利用数值模拟系统进行扰动波直接模拟，获取非定常流场变量ϕ；
S3，分析流场得到扰动场ϕ'，并获取湍动能分布特征；
S4，设置湍动能入口边界；
S5，计算湍流/转捩。
进一步地，在步骤S1中，建立计算域，设定边界条件，采用数值计算平台，利用二维Navier-Stokes方程为控制方程，展开流动计算，计算得到没有施加扰动的基础流场解，获得流场空间区域的密度、速度、温度与压力信息。
进一步地，在步骤S2中，基于步骤S1获取的基础流场，开展自由流扰动波的非定常传播演化过程的直接数值模拟，获取流场中扰动波特征及湍动能分布。
进一步地，在步骤S3中，当扰动传播至整个流场形成周期解后，利用步骤S2获取的 瞬时非定常流场变量ϕ减去步骤S1获取的定常基本流场变量
Figure 975763DEST_PATH_IMAGE009
得到扰动场ϕ'，通过对脉 动量的傅里叶频谱分析，获取流场中扰动波演化和分布规律；根据脉动速度换算湍动能k， 获取湍动能分布规律
Figure 34986DEST_PATH_IMAGE010
Figure 388476DEST_PATH_IMAGE011
其中，k为湍动能，u'为x向速度脉动，v'为y向速度脉动。
进一步地，在步骤S4中，根据步骤S3中计算湍动能，选择钝前缘下游位置，沿壁面法向设置湍动能入口边界，按照步骤S3中结果提取指定湍动能入口边界的物理量，根据RANS和DNS计算中自由流湍动能的比值线性放大指定位置处的湍动能值，并设置RANS方程中k值，将其固定，不随流场的时间推进而改变；
Figure 747913DEST_PATH_IMAGE012
其中，
Figure 217072DEST_PATH_IMAGE013
为流体密度，
Figure 824245DEST_PATH_IMAGE014
为层流粘性系数，
Figure 466578DEST_PATH_IMAGE015
为湍动能，
Figure 149364DEST_PATH_IMAGE016
为模型系数，
Figure 101008DEST_PATH_IMAGE017
为湍流 粘性系数，
Figure 869244DEST_PATH_IMAGE018
为粘性应力，
Figure 315269DEST_PATH_IMAGE019
为应变率，
Figure 898566DEST_PATH_IMAGE020
为模型系数，
Figure 709527DEST_PATH_IMAGE021
为湍流耗散率，下标t表示湍流， 下标j表示计算维度，下标i表示空间维度，xj表示空间方向坐标，当j取1，2，3时，对应的xj表 示x,y,z三个空间方向坐标。
进一步地，在步骤S5中，基于RANS方法结合包含湍动能的湍流模型/转捩模型进行湍流/转捩计算，获得基于准确湍动能入口边界的湍流/转捩模型计算结果，所述计算结果包括全计算域的密度、温度、速度流场信息，获取扰流物体表面的Cf分布或热流分布。
进一步地，在步骤S1中，在无扰动、无体积力与外部热源的情况下，采用的控制方程具体表示如下：
Figure 27376DEST_PATH_IMAGE022
Figure 995201DEST_PATH_IMAGE023
Figure 652578DEST_PATH_IMAGE024
其中
Figure 431178DEST_PATH_IMAGE025
分别是密度、x方向速度、y方向速度、z方向速度和单位质量气体总 能量，
Figure 488520DEST_PATH_IMAGE026
分别为不同方向的粘性应力分量，
Figure 214031DEST_PATH_IMAGE027
分别为x和y方向的 热流，
Figure 522653DEST_PATH_IMAGE028
为压力，Q表示守恒变量，f和g分别表示x和y向矢通量，并采用不小于5阶精度的高 阶精度计算格式，以及能够分辨流场中扰动波尺度的网格，计算得到没有施加扰动的基础 流场。
进一步地，在步骤S2中，在对流程进行扰动传播模拟时选用声波扰动，考虑二维平面声波：
Figure 987001DEST_PATH_IMAGE029
声波振幅满足如下关系：
Figure 951546DEST_PATH_IMAGE030
其中，
Figure 277485DEST_PATH_IMAGE031
无量纲波数，
Figure 955460DEST_PATH_IMAGE032
是一小量，为扰动波振幅参数，
Figure 75863DEST_PATH_IMAGE033
为扰动波无量纲圆频 率；快、慢声波：
Figure 793283DEST_PATH_IMAGE034
，“+”为快声波，“−”为慢声波，
Figure 172180DEST_PATH_IMAGE035
为自由来流马赫数； 自由流声波参数：振幅
Figure 393077DEST_PATH_IMAGE036
，频率f = 100kHz；
Figure 215540DEST_PATH_IMAGE037
为x方向速度脉动幅值，
Figure 406874DEST_PATH_IMAGE038
为压力脉动幅值，
Figure 340195DEST_PATH_IMAGE039
为y方向脉动幅值；
通过在计算域入口边界引入声波，采用非定常计算方法直接数值模拟小扰动经过激波在边界层中引起扰动的过程。
与传统的湍动能入口边界设置技术相比，本发明的有益效果是：
（1）本发明实施例不需对湍流模型中湍动能方程进行修改，通过入口边界条件的设置即可避免湍动能方程在激波下游的计算偏差。
（2）本发明实施例基于直接数值模拟结果设置更加准确的湍动能入口边界条件，具有计算结果合理、可实现性高的优势，并且在钝平板的Mach6湍流计算中表现出比现有方法更高的可靠性。
（3）本发明实施例可简单推广至三维流动情况。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为RANS和DNS湍流动能计算对比图；其中，(a)为RANS方程计算湍流动能分布；（b）为DNS计算扰动波动能分布；
图2为湍动能入口边界位置；其中，(a)为传统入口条件设置，（b）为本发明实施例的入口条件设置；
图3为湍动能入口边界设置的技术路线图；
图4为基础流场图；
图5为瞬时扰动场图；
图6为湍动能分布；
图7为x=10mm 湍动能沿法向分布；
图8为湍动能分布；其中，(a)为采用传统方法的局部放大图，（b）为采用本发明实施例方法的局部放大图；（c）为采用传统方法的全场图，（d）为采用本发明实施例方法的全场图；
图9为湍流粘性系数分布；其中，(a)为采用传统方法的湍流粘性系数，（b）为采用新方法的湍流粘性系数；
图10为本发明实施例的方法步骤流程图。
具体实施方式
本说明书中所有实施例公开的所有特征，或隐含公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合和/或扩展、替换。
如图1~10所示，一种高超声速钝前缘绕流湍动能入口边界设置方法，包括如下步骤：
S1，计算基础流场，获取定常基本流场变量
Figure 664866DEST_PATH_IMAGE009
；
S2，利用数值模拟系统进行扰动波直接模拟，获取非定常流场变量ϕ；
S3，分析流场得到扰动场ϕ'，并获取湍动能分布特征；
S4，设置湍动能入口边界；
S5，计算湍流/转捩。
进一步地，在步骤S1中，建立计算域，设定边界条件，采用数值计算平台，利用二维Navier-Stokes方程为控制方程，展开流动计算，计算得到没有施加扰动的基础流场解，获得流场空间区域的密度、速度、温度与压力信息。
进一步地，在步骤S2中，基于步骤S1获取的基础流场，开展自由流扰动波的非定常传播演化过程的直接数值模拟，获取流场中扰动波特征及湍动能分布。
进一步地，在步骤S3中，当扰动传播至整个流场形成周期解后，利用步骤S2获取的 瞬时非定常流场变量ϕ减去步骤S1获取的定常基本流场变量
Figure 861492DEST_PATH_IMAGE009
得到扰动场ϕ'，通过对脉 动量的傅里叶频谱分析，获取流场中扰动波演化和分布规律；根据脉动速度换算湍动能k， 获取湍动能分布规律
Figure 537193DEST_PATH_IMAGE010
Figure 946309DEST_PATH_IMAGE011
其中，k为湍动能，u'为x向速度脉动，v'为y向速度脉动。
进一步地，在步骤S4中，根据步骤S3中计算湍动能，选择钝前缘下游位置，沿壁面法向设置湍动能入口边界，按照步骤S3中结果提取指定湍动能入口边界的物理量，根据RANS和DNS计算中自由流湍动能的比值线性放大指定位置处的湍动能值，并设置RANS方程中k值，将其固定，不随流场的时间推进而改变；
Figure 453382DEST_PATH_IMAGE012
其中，
Figure 555331DEST_PATH_IMAGE013
为流体密度，
Figure 721257DEST_PATH_IMAGE014
为层流粘性系数，
Figure 996381DEST_PATH_IMAGE015
为湍动能，
Figure 577535DEST_PATH_IMAGE016
为模型系数，
Figure 99652DEST_PATH_IMAGE017
为湍流 粘性系数，
Figure 766257DEST_PATH_IMAGE018
为粘性应力，
Figure 766443DEST_PATH_IMAGE019
为应变率，
Figure 202103DEST_PATH_IMAGE020
为模型系数，
Figure 895121DEST_PATH_IMAGE021
为湍流耗散率，下标t表示湍流， 下标j表示计算维度，下标i表示空间维度，xj表示空间方向坐标，当j取1，2，3时，对应的xj表 示x,y,z三个空间方向坐标。
进一步地，在步骤S5中，基于RANS方法结合包含湍动能的湍流模型/转捩模型进行湍流/转捩计算，获得基于准确湍动能入口边界的湍流/转捩模型计算结果，所述计算结果包括全计算域的密度、温度、速度流场信息，获取扰流物体表面的Cf分布或热流分布。
进一步地，在步骤S1中，在无扰动、无体积力与外部热源的情况下，采用的控制方程具体表示如下：
Figure 845760DEST_PATH_IMAGE022
Figure 906427DEST_PATH_IMAGE023
Figure 931015DEST_PATH_IMAGE024
其中
Figure 60514DEST_PATH_IMAGE025
分别是密度、x方向速度、y方向速度、z方向速度和单位质量气体总 能量，
Figure 436131DEST_PATH_IMAGE026
分别为不同方向的粘性应力分量，
Figure 43699DEST_PATH_IMAGE027
分别为x和y方向的 热流，
Figure 250689DEST_PATH_IMAGE028
为压力，Q表示守恒变量，f和g分别表示x和y向矢通量，并采用不小于5阶精度的高 阶精度计算格式，以及能够分辨流场中扰动波尺度的网格，计算得到没有施加扰动的基础 流场。
进一步地，在步骤S2中，在对流程进行扰动传播模拟时选用声波扰动，考虑二维平面声波：
Figure 36243DEST_PATH_IMAGE029
声波振幅满足如下关系：
Figure 882845DEST_PATH_IMAGE030
其中，
Figure 44836DEST_PATH_IMAGE031
无量纲波数，
Figure 561793DEST_PATH_IMAGE032
是一小量，为扰动波振幅参数，
Figure 580564DEST_PATH_IMAGE033
为扰动波无量纲圆频 率；快、慢声波：
Figure 914462DEST_PATH_IMAGE034
，“+”为快声波，“−”为慢声波，
Figure 411303DEST_PATH_IMAGE035
为自由来流马赫数； 自由流声波参数：振幅
Figure 530569DEST_PATH_IMAGE036
，频率f = 100kHz；
Figure 907192DEST_PATH_IMAGE037
为x方向速度脉动幅值，
Figure 541436DEST_PATH_IMAGE038
为压力脉动幅值，
Figure 28918DEST_PATH_IMAGE039
为y方向脉动幅值；
通过在计算域入口边界引入声波，采用非定常计算方法直接数值模拟小扰动经过激波在边界层中引起扰动的过程。
本发明实施例采用的新湍动能入口边界方法，可显著改善钝平板上游流场中湍动能的分布，避免了RANS方程计算中头部激波对湍动能方程的虚假影响，较好地计算了边界层中的湍动能和湍流粘性系数，更加合理反映了湍流效应在流场中的作用，相比传统方法而言更具优势。
本发明实施例解决的技术问题是高超声速钝前缘二维湍流RANS计算时的湍动能 入口边界的设置问题。在小扰动波直接数值模拟过程中发现，高超声速钝前缘平板上湍动 能分布和现有RANS方法计算的湍动能分布不同，现有入口边界直接设置自由流湍流度的方 法在RANS计算过程中出现对钝前缘区域湍动能计算的错误，现有上游入口边界湍动能的设 置方法对高超声速钝前缘流动情况不再适用。基于该技术问题，证明了现有在入口边界直 接设置湍流度的方法不在适用于高超声速钝前缘绕流（见现有技术部分），提出了适用于该 情况下的湍动能边界条件设置方法，并将该方法与基于Chant2.0高超声速计算平台的
Figure 799428DEST_PATH_IMAGE040
SST模型相耦合，获得了能够模拟高超声速钝前缘流动的湍流/转捩预测技术。该技术具有 低成本、易实现和高可靠性的特点；与现有入口湍动能设置方法相比，对于钝前缘几何外形 的湍流/转捩流动的计算更加合理。
利用直接数值模拟DNS的结果，对高超声速钝前缘流场的湍动能进行计算，提出新的自由流湍动能计算方法，将其用于RANS方程计算中，对高超声速钝前缘边界层湍流/转捩进行了数值模拟。
在本发明的实施例中，湍动能入口位置的选取，在所有CFD计算中都需要对计算域的边界设置流动变量的边界条件，基于RANS方程的湍流/转捩计算时需要设置入口或自由流湍动能，以满足控制方程的求解需求。入口边界（也称外边界）通常设置于满足自由流流动条件的位置上，采用指定变量的方式设置边界层条件，下游流场中的变量通过流动控制方程（RANS）计算获得。
由于RANS方程在高超声速钝前缘流动计算时的误差较大，不能真实反映钝前缘附近边界层中的湍动能分布，可以选择在钝前缘下游为湍动能入口边界，在此指定湍动能来避免RANS方程在前缘处的计算误差，以达到更加真实的计算效果。图2左侧图为现有的入口条件设置位置，右侧图为新的湍动能入口位置。
在本发明的实施例中，核心工作是基于直接数值模拟结果给出湍动能入口边界条件。
通过对钝前缘流场的扰动波直接数值模拟，获取扰动在流场中的分布，选择湍动能入口边界位置，参考该位置上的扰动能分布及与自由流湍动能的关系，根据该位置处扰动能设置RANS方程k的上游边界条件，进而完成对高超声速钝前缘绕流湍流/转捩的计算。
在本发明的实施中，对各步骤的解释如下：
步骤一，基础流场计算。
根据二维模型外形绘制计算网格（建立计算域），设定边界条件，采用Chant2.0数值计算平台，以二维Navier-Stokes方程为控制方程，展开流动计算。在无扰动、无体积力与外部热源的情况下，控制方程具体表示如下：
Figure 97685DEST_PATH_IMAGE041
其中
Figure 471422DEST_PATH_IMAGE042
分别是密度、x方向速度、y方向速度、z方向速度和单位质量气体总能 量，
Figure 247748DEST_PATH_IMAGE043
分别为不同方向的粘性应力分量，
Figure 138344DEST_PATH_IMAGE044
分别为x和y方向的热 流，
Figure 122349DEST_PATH_IMAGE045
为压力，Q表示守恒变量，f和g分别表示x和y向矢通量。采用可稳定捕捉激波间断和精 确计算边界层的高阶精度计算格式（不小于5阶精度）以及可以分辨流场中扰动波尺度的网 格，计算得到没有施加扰动的基础流场解
Figure 386977DEST_PATH_IMAGE046
，获得流场空间区域的密度、速度、温度与压力 等信息。
步骤二中，扰动波直接数值模拟。
基于步骤一种获取的定常基础流场，开展自由流扰动波的非定常传播演化过程的直接数值模拟，计算非定常流场解ϕ，获取扰动场ϕ'，同时提取流场中扰动波特征及湍动能分布。自由流中扰动包含声波、涡波和熵波等三种形式，在高超声速流动中，尤其是在地面风洞试验中，声波是扰动的主要成份，通常边界层触发转捩至湍流是由自由流中声波扰动引起。因此，在对流程进行扰动传播模拟时应该选用声波扰动。考虑二维平面声波：
Figure 685103DEST_PATH_IMAGE047
（6）
声波振幅满足如下色散关系：
Figure 633468DEST_PATH_IMAGE048
上式中，
Figure 525725DEST_PATH_IMAGE049
为无量纲波数，
Figure 356278DEST_PATH_IMAGE050
是一小量，为扰动波振幅参数。
Figure 271144DEST_PATH_IMAGE051
为扰动波无量 纲圆频率，快、慢声波：
Figure 323283DEST_PATH_IMAGE052
，“+”为快声波，“−”为慢声波。
Figure 930982DEST_PATH_IMAGE053
为自由来流 马赫数；自由流声波参数：振幅
Figure 248830DEST_PATH_IMAGE054
，频率f = 100kHz。
Figure 419918DEST_PATH_IMAGE055
为x方向速度脉动幅值，
Figure 139612DEST_PATH_IMAGE056
为压力脉动幅值，
Figure 121474DEST_PATH_IMAGE057
为y方向脉动幅值；
通过在计算域入口边界引入声波，采用非定常计算方法直接数值模拟（DNS）小扰动经过激波在边界层中引起扰动的过程。
步骤三，流场分析，获取湍动能分布特征。
当扰动传播至整个流场形成周期解后，利用步骤二获取的瞬时非定常流场变量ϕ减去步骤一获取的定常基本流场变量Ф得到扰动场ϕ'，通过对脉动量的傅里叶频谱分析，获取流场中扰动波演化和分布规律。根据脉动速度换算湍动能k，获取湍动能分布规律。
Figure 175887DEST_PATH_IMAGE058
其中，k为湍动能，u'为x向速度脉动，v'为y向速度脉动。
步骤四，设置湍动能入口边界。
根据步骤三中计算湍动能，选择钝前缘下游位置，沿壁面法向设置湍动能入口边界，按照步骤三中结果提取指定湍动能入口边界的物理量，根据RANS和DNS计算中自由流湍动能的比值线性放大指定位置处的湍动能值，并设置RANS方程中k值，将其固定，不随流场的时间推进而改变。
Figure 901397DEST_PATH_IMAGE059
其中，
Figure 475598DEST_PATH_IMAGE060
为流体密度，
Figure 892278DEST_PATH_IMAGE061
为湍动能，
Figure 122402DEST_PATH_IMAGE062
为模型系数，
Figure 166451DEST_PATH_IMAGE063
为湍流粘性系数，
Figure 798420DEST_PATH_IMAGE064
为粘性 应力，
Figure 433670DEST_PATH_IMAGE065
为应变率，
Figure 416669DEST_PATH_IMAGE066
为模型系数，
Figure 998829DEST_PATH_IMAGE067
为湍流耗散率，下标t表示湍流，下标j表示计算维 度。
步骤五，湍流/转捩计算。
基于RANS方法结合包含湍动能的湍流模型和转捩模型进行湍流/转捩计算，最终获得基于准确湍动能入口边界的湍流/转捩模型计算结果，包括全计算域的密度、温度、速度等流场信息，获取扰流物体表面的Cf分布或热流分布。
在本发明的其他实施例中，以头部半径5毫米圆弧为钝前缘的平板Mach 6 流动为例，采用本发明方法进行示范性的湍动能入口边界设置。平板长度为1米，流动条件为常规高超声速风洞流动，自由流马赫数为Ma=6，单位雷诺数为Re=1×107/m，来流静温为T=61K，风洞自由流噪声水平为3%，在实施例中具体包括如下步骤：
第一步，基础流场计算：以N-S方程为控制方程，采用5阶精度WCNS-E-5格式对无粘项进行离散，采用6阶中心格式对粘性项进行离散，采用LU-SGS格式进行时间推进。网格计算域包括激波，法向171个网格，流向1101个网格。图4为计算收敛后的基础流场解的马赫数，可以清晰地看到头部弓形激波，及马赫数的空间分布。
第二步，扰动波直接数值模拟：通过在基础流场的外边界引入平面声波来直接数值模拟扰动穿过激波在边界层中产生新扰动的过程。入射声波频率设为50kHz，无量纲幅值设为5e-4。采用三阶精度R-K时间推进方法模拟扰动的产生和演化过程。
第三步，获取湍动能分布特征：通过非定常扰动场和定常基础流场相减获得扰动场。图5为某时刻的瞬态动能脉动场。可以看到在远离壁面的自由流区和靠近壁面的边界层区都出现了扰动波。通过记录一个周期的脉动量，采用FFT变化得到扰动的幅值分布，包括湍动能幅值分布。图6为湍动能脉动幅值分布。
第四步，湍动能入口边界设置：选择平板前缘下游附近一个站位（x=10mm）的湍动能，获取沿法向的分布，见图7中曲线，以此作为RANS方法计算时的湍动能入口边界。RANS中无量纲自由流湍动能为1.35e-3，DNS模拟中自由流湍动能为6e-8，需要根据二者比例线性放大x=10mm处的湍动能。激波前的自由流区域湍动能直接赋值来流湍动能。
第五步，湍流/转捩RANS方程计算：将第四步中获得的湍动能分布作为x=10mm站位 处的k值，采用基于RANS方程的湍动能输运方程计算k在全流场边界层中的分布，借助
Figure 219726DEST_PATH_IMAGE068
SST模型计算湍流效应获得钝平板的湍流解。图8为传统湍动能来流边界赋值方法和改进后 的湍动能分布比较，图9为湍流粘性系数的比较。
从计算结果来看，本发明实施例采用的新湍动能入口边界方法，可显著改善钝平板上游流场中湍动能的分布，避免了RANS方程计算中头部激波对湍动能方程的虚假影响，较好地计算了边界层中的湍动能和湍流粘性系数，更加合理反映了湍流效应在流场中的作用，相比传统方法而言更具优势。
本发明实施例的创新点在于提出了新的湍动能入口边界条件设置方法，在基于RANS方程和湍动能k方程的湍流/转捩数值模拟中，新的湍动能入口边界设置方法相比传统方法而言，能更加合理、准确地反映湍动能在自由流区和边界层内的分布，为下游湍动能计算提供更加合理的上游边界值，更好地计算湍流粘性效应。
本发明功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，在一台计算机设备（可以是个人计算机，服务器，或者网络设备等）以及相应的软件中执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、或者光盘等各种可以存储程序代码的介质，进行测试或者实际的数据在程序实现中存在于只读存储器（Random Access Memory，RAM）、随机存取存储器（Random Access Memory，RAM）等。

Claims (8)

1.一种高超声速钝前缘绕流湍动能入口边界设置方法，其特征在于，包括如下步骤：
S1，计算基础流场，获取定常基本流场变量
Figure 873745DEST_PATH_IMAGE001
；
S2，利用数值模拟系统进行扰动波直接模拟，获取非定常流场变量ϕ；
S3，分析流场得到扰动场ϕ'，并获取湍动能分布特征；
S4，设置湍动能入口边界；
S5，计算湍流/转捩。
2.根据权利要求1所述的一种高超声速钝前缘绕流湍动能入口边界设置方法，其特征在于，在步骤S1中，建立计算域，设定边界条件，采用数值计算平台，利用二维Navier-Stokes方程为控制方程，展开流动计算，计算得到没有施加扰动的基础流场解，获得流场空间区域的密度、速度、温度与压力信息。
3.根据权利要求1所述的一种高超声速钝前缘绕流湍动能入口边界设置方法，其特征在于，在步骤S2中，基于步骤S1获取的基础流场，开展自由流扰动波的非定常传播演化过程的直接数值模拟，获取流场中扰动波特征及湍动能分布。
4.根据权利要求1所述的一种高超声速钝前缘绕流湍动能入口边界设置方法，其特征 在于，在步骤S3中，当扰动传播至整个流场形成周期解后，利用步骤S2获取的瞬时非定常流 场变量ϕ减去步骤S1获取的定常基本流场变量
Figure 980766DEST_PATH_IMAGE001
得到扰动场ϕ'，通过对脉动量的傅里叶 频谱分析，获取流场中扰动波演化和分布规律；根据脉动速度换算湍动能k，获取湍动能分 布规律
Figure 167028DEST_PATH_IMAGE002
Figure 827816DEST_PATH_IMAGE003
其中，k为湍动能，u'为x向速度脉动，v'为y向速度脉动。
5.根据权利要求1所述的一种高超声速钝前缘绕流湍动能入口边界设置方法，其特征在于，在步骤S4中，根据步骤S3中计算湍动能，选择钝前缘下游位置，沿壁面法向设置湍动能入口边界，按照步骤S3中结果提取指定湍动能入口边界的物理量，根据RANS和DNS计算中自由流湍动能的比值线性放大指定位置处的湍动能值，并设置RANS方程中k值，将其固定，不随流场的时间推进而改变；
Figure 297981DEST_PATH_IMAGE004
其中，
Figure 323706DEST_PATH_IMAGE005
为流体密度，
Figure 121897DEST_PATH_IMAGE006
为层流粘性系数，
Figure 711010DEST_PATH_IMAGE007
为湍动能，
Figure 911048DEST_PATH_IMAGE008
为模型系数，
Figure 107674DEST_PATH_IMAGE009
为湍流粘性 系数，
Figure 517795DEST_PATH_IMAGE010
为粘性应力，
Figure 520386DEST_PATH_IMAGE011
为应变率，
Figure 450296DEST_PATH_IMAGE012
为模型系数，
Figure 335600DEST_PATH_IMAGE013
为湍流耗散率，下标t表示湍流，下标j 表示计算维度，下标i表示空间维度，xj表示空间方向坐标，当j取1，2，3时，对应的xj表示x, y,z三个空间方向坐标。
6.根据权利要求1所述的一种高超声速钝前缘绕流湍动能入口边界设置方法，其特征在于，在步骤S5中，基于RANS方法结合包含湍动能的湍流模型/转捩模型进行湍流/转捩计算，获得基于准确湍动能入口边界的湍流/转捩模型计算结果，所述计算结果包括全计算域的密度、温度、速度流场信息，获取扰流物体表面的Cf分布或热流分布。
7.根据权利要求2所述的一种高超声速钝前缘绕流湍动能入口边界设置方法，其特征在于，在步骤S1中，在无扰动、无体积力与外部热源的情况下，采用的控制方程具体表示如下：
Figure 842805DEST_PATH_IMAGE014
Figure 524453DEST_PATH_IMAGE015
Figure 620454DEST_PATH_IMAGE016
其中
Figure 486779DEST_PATH_IMAGE017
分别是密度、x方向速度、y方向速度、z方向速度和单位质量气体总能量，
Figure 356646DEST_PATH_IMAGE018
分别为不同方向的粘性应力分量，
Figure 887990DEST_PATH_IMAGE019
分别为x和y方向的热流，
Figure 917126DEST_PATH_IMAGE020
为压力，Q表示守恒变量，f和g分别表示x和y向矢通量，并采用不小于5阶精度的高阶精度计 算格式，以及能够分辨流场中扰动波尺度的网格，计算得到没有施加扰动的基础流场。
8.根据权利要求3所述的一种高超声速钝前缘绕流湍动能入口边界设置方法，其特征在于，在步骤S2中，在对流程进行扰动传播模拟时选用声波扰动，考虑二维平面声波：
Figure 95298DEST_PATH_IMAGE021
声波振幅满足如下关系：
Figure 311515DEST_PATH_IMAGE022
其中，
Figure 646551DEST_PATH_IMAGE023
无量纲波数，
Figure 405559DEST_PATH_IMAGE024
是一小量，为扰动波振幅参数，
Figure 613687DEST_PATH_IMAGE025
为扰动波无量纲圆频率； 快、慢声波：
Figure 710343DEST_PATH_IMAGE026
，“+”为快声波，“−”为慢声波，
Figure 334223DEST_PATH_IMAGE027
为自由来流马赫数；自由 流声波参数：振幅
Figure 337951DEST_PATH_IMAGE028
，频率f = 100kHz；
Figure 841613DEST_PATH_IMAGE029
为x方向速度脉动幅值，
Figure 298002DEST_PATH_IMAGE030
为压 力脉动幅值，
Figure 725573DEST_PATH_IMAGE031
为y方向脉动幅值；
通过在计算域入口边界引入声波，采用非定常计算方法直接数值模拟小扰动经过激波在边界层中引起扰动的过程。
CN202110387351.2A 2021-04-12 2021-04-12 一种高超声速钝前缘绕流湍动能入口边界设置方法 Active CN112765736B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110387351.2A CN112765736B (zh) 2021-04-12 2021-04-12 一种高超声速钝前缘绕流湍动能入口边界设置方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110387351.2A CN112765736B (zh) 2021-04-12 2021-04-12 一种高超声速钝前缘绕流湍动能入口边界设置方法

Publications (2)

Publication Number Publication Date
CN112765736A true CN112765736A (zh) 2021-05-07
CN112765736B CN112765736B (zh) 2021-06-29

Family

ID=75691426

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110387351.2A Active CN112765736B (zh) 2021-04-12 2021-04-12 一种高超声速钝前缘绕流湍动能入口边界设置方法

Country Status (1)

Country Link
CN (1) CN112765736B (zh)

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113111608A (zh) * 2021-05-10 2021-07-13 中国空气动力研究与发展中心计算空气动力研究所 一种新型局部湍流脉动生成方法
CN113408218A (zh) * 2021-06-28 2021-09-17 西北工业大学 一种基于扰动方程的流动噪声模拟方法
CN113505543A (zh) * 2021-06-18 2021-10-15 中国空气动力研究与发展中心计算空气动力研究所 一种飞行器壁面热流分析方法
CN113946904A (zh) * 2021-08-31 2022-01-18 中国航天空气动力技术研究院 一种大尺寸低噪声喷管设计方法
CN113971379A (zh) * 2021-10-28 2022-01-25 中国人民解放军国防科技大学 超声速湍流火焰面/进度变量模型的温度求解简化方法
CN114154441A (zh) * 2022-02-10 2022-03-08 中国空气动力研究与发展中心计算空气动力研究所 一种飞行器的环境湍流场生成与模拟计算方法
CN115510694A (zh) * 2022-11-23 2022-12-23 中国航天三江集团有限公司 基于引射喷流的湍流边界层气动光学效应抑制方法
CN115859480A (zh) * 2023-02-08 2023-03-28 中国空气动力研究与发展中心计算空气动力研究所 基于确定发动机入口边界条件的气动分析方法及装置
CN116229021A (zh) * 2023-05-08 2023-06-06 中国空气动力研究与发展中心计算空气动力研究所 一种浸没边界虚拟网格嵌入方法、装置、设备及介质

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108388742A (zh) * 2018-03-05 2018-08-10 北京空间技术研制试验中心 一种空间飞行器带压分离的仿真分析方法
CN109766627A (zh) * 2019-01-08 2019-05-17 西南交通大学 一种基于滑板间距的受电弓非定常特性的分析方法
CN109969374A (zh) * 2019-04-09 2019-07-05 中国空气动力研究与发展中心计算空气动力研究所 用于高超声速边界层转捩研究的标模气动布局及设计方法
CN112052632A (zh) * 2020-07-27 2020-12-08 空气动力学国家重点实验室 一种高超声速流向转捩预测方法
CN112580272A (zh) * 2020-12-14 2021-03-30 中国市政工程华北设计研究总院有限公司 一种基于数值模拟的lng空温式气化器的优化设计方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108388742A (zh) * 2018-03-05 2018-08-10 北京空间技术研制试验中心 一种空间飞行器带压分离的仿真分析方法
CN109766627A (zh) * 2019-01-08 2019-05-17 西南交通大学 一种基于滑板间距的受电弓非定常特性的分析方法
CN109969374A (zh) * 2019-04-09 2019-07-05 中国空气动力研究与发展中心计算空气动力研究所 用于高超声速边界层转捩研究的标模气动布局及设计方法
CN112052632A (zh) * 2020-07-27 2020-12-08 空气动力学国家重点实验室 一种高超声速流向转捩预测方法
CN112580272A (zh) * 2020-12-14 2021-03-30 中国市政工程华北设计研究总院有限公司 一种基于数值模拟的lng空温式气化器的优化设计方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
STEVEN P. SCHNEIDER: ""Developing mechanism-based methods for estimating hypersonic boundary-layer transition in fight: The role of quiet tunnels"", 《PROGRESS IN AEROSPACE SCIENCES》 *
吴子牛 等: ""高超声速飞行器流动特征分析"", 《航空学报》 *
涂国华 等: ""MF-1钝锥边界层稳定性及转捩天地相关性研究"", 《中国科学：物理学 力学 天文学》 *
王义乾: ""基于边界层转捩直接数值模拟的湍流生成与维持机理研究"", 《中国博士学位论文全文数据库 基础科学辑（月刊）》 *
王雄 等: ""圆孔射流近场湍流特性DNS与RANS模拟的对比研究"", 《能源工程》 *
童福林 等: ""超声速膨胀角入射激波／湍流边界层干扰直接数值模拟"", 《航空学报》 *

Cited By (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113111608B (zh) * 2021-05-10 2022-06-28 中国空气动力研究与发展中心计算空气动力研究所 一种新型局部湍流脉动生成方法
CN113111608A (zh) * 2021-05-10 2021-07-13 中国空气动力研究与发展中心计算空气动力研究所 一种新型局部湍流脉动生成方法
CN113505543B (zh) * 2021-06-18 2023-05-26 中国空气动力研究与发展中心计算空气动力研究所 一种飞行器壁面热流分析方法
CN113505543A (zh) * 2021-06-18 2021-10-15 中国空气动力研究与发展中心计算空气动力研究所 一种飞行器壁面热流分析方法
CN113408218A (zh) * 2021-06-28 2021-09-17 西北工业大学 一种基于扰动方程的流动噪声模拟方法
CN113408218B (zh) * 2021-06-28 2022-09-16 西北工业大学 一种基于扰动方程的流动噪声模拟方法
CN113946904A (zh) * 2021-08-31 2022-01-18 中国航天空气动力技术研究院 一种大尺寸低噪声喷管设计方法
CN113946904B (zh) * 2021-08-31 2024-06-11 中国航天空气动力技术研究院 一种大尺寸低噪声喷管设计方法
CN113971379A (zh) * 2021-10-28 2022-01-25 中国人民解放军国防科技大学 超声速湍流火焰面/进度变量模型的温度求解简化方法
CN113971379B (zh) * 2021-10-28 2024-05-28 中国人民解放军国防科技大学 超声速湍流火焰面/进度变量模型的温度求解简化方法
CN114154441A (zh) * 2022-02-10 2022-03-08 中国空气动力研究与发展中心计算空气动力研究所 一种飞行器的环境湍流场生成与模拟计算方法
CN115510694A (zh) * 2022-11-23 2022-12-23 中国航天三江集团有限公司 基于引射喷流的湍流边界层气动光学效应抑制方法
CN115859480B (zh) * 2023-02-08 2023-05-02 中国空气动力研究与发展中心计算空气动力研究所 基于确定发动机入口边界条件的气动分析方法及装置
CN115859480A (zh) * 2023-02-08 2023-03-28 中国空气动力研究与发展中心计算空气动力研究所 基于确定发动机入口边界条件的气动分析方法及装置
CN116229021A (zh) * 2023-05-08 2023-06-06 中国空气动力研究与发展中心计算空气动力研究所 一种浸没边界虚拟网格嵌入方法、装置、设备及介质
CN116229021B (zh) * 2023-05-08 2023-08-25 中国空气动力研究与发展中心计算空气动力研究所 一种浸没边界虚拟网格嵌入方法、装置、设备及介质

Also Published As

Publication number Publication date
CN112765736B (zh) 2021-06-29

Similar Documents

Publication Publication Date Title
CN112765736B (zh) 一种高超声速钝前缘绕流湍动能入口边界设置方法
Hsiao et al. Scaling of tip vortex cavitation inception noise with a bubble dynamics model accounting for nuclei size distribution
Bakica et al. Accurate assessment of ship-propulsion characteristics using CFD
Sidebottom et al. A parametric study of turbulent flow past a circular cylinder using large eddy simulation
Zhang et al. Flow-induced noise analysis for 3D trash rack based on LES/Lighthill hybrid method
Zhang et al. Computation of vortical flow and flow induced noise by large eddy simulation with FW-H acoustic analogy and Powell vortex sound theory
Suzuki et al. Unsteady simulations of a fan/outlet-guide-vane system: tone–noise computation
Wang et al. Large-eddy-simulation prediction of an installed jet flow and noise with experimental validation
Liu et al. Design of directional probes for high-frequency turbine measurements
Hristov et al. Poststall hysteresis and flowfield unsteadiness on a NACA 0012 airfoil
Ziadé et al. Shear layer development, separation, and stability over a low-Reynolds number airfoil
Weckmüller et al. Cfd/caa coupling applied to dlr uhbr-fan: Comparison to experimantal data
Qin et al. Numerical simulation of hydrodynamic and noise characteristics for a blended-wing-body underwater glider
Zanon et al. Assessment of the broadband noise from an unducted axial fan including the effect of the inflow turbulence
Mueller Large eddy simulation of cross-flow around a square rod at incidence with application to tonal noise prediction
de Laborderie et al. Evaluation of a cascade-based acoustic model for fan tonal noise prediction
Koch et al. Numerical aeroacoustic analysis of a linear compressor cascade with tip gap
Liliedahl et al. Prediction of aeroacoustic resonance in cavities of hole-pattern stator seals
Ewert et al. A CAA based approach to tone haystacking
Zhang et al. Predictions of fan broadband noise using lifting surface method
Zhong et al. An efficient computation of cascade-gust interaction noise based on a hybrid analytical and boundary element method
Margha et al. Dynamic transition from regular to mach reflection over a moving wedge
Wang et al. A multiscale Euler–Lagrange model for high-frequency cavitation noise prediction
Vyas et al. Sidewall effects on exact Reynolds-stress budgets in an impinging shock wave/boundary layer interaction
Rendu et al. Investigation of shock-wave/boundary-layer interaction on aeroelastic stability: towards active control

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