CN109522663A - 一种基于有限元模型的汽轮机转子蠕变损伤实时计算方法 - Google Patents

一种基于有限元模型的汽轮机转子蠕变损伤实时计算方法 Download PDF

Info

Publication number
CN109522663A
CN109522663A CN201811424356.2A CN201811424356A CN109522663A CN 109522663 A CN109522663 A CN 109522663A CN 201811424356 A CN201811424356 A CN 201811424356A CN 109522663 A CN109522663 A CN 109522663A
Authority
CN
China
Prior art keywords
rotor
stress
real
temperature
creep
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
CN201811424356.2A
Other languages
English (en)
Other versions
CN109522663B (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.)
University of Jinan
Original Assignee
University of Jinan
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 University of Jinan filed Critical University of Jinan
Priority to CN201811424356.2A priority Critical patent/CN109522663B/zh
Publication of CN109522663A publication Critical patent/CN109522663A/zh
Application granted granted Critical
Publication of CN109522663B publication Critical patent/CN109522663B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design
    • 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/06—Power analysis or power optimisation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)

Abstract

本发明提出了一种基于有限元模型的汽轮机转子蠕变损伤实时计算方法，包括：通过有限元分析软件计算转子的温度场和应力场；提取关键部位的温度和应力数据；将蠕变损伤计算相关的温度、压力、转速以及应力数据进行归一化处理；采用多元回归方法建立蠕变应力实时计算模型；解析并拟合应力‑损伤函数关系，实时计算蠕变损伤；建立蠕变损伤实时计算系统结构。本发明可以解决现有技术中转子蠕变损伤损伤实时计算的问题，并且具有较高的精度。

Description

一种基于有限元模型的汽轮机转子蠕变损伤实时计算方法
技术领域
本发明涉及汽轮机转子蠕变损伤计算与评估领域，特别是涉及有限元数值计算方法和多元回归实时计算方法。
背景技术
转子是汽轮机机组的核心部件，通常运行在高温高压的恶劣工作条件，容易产生裂纹，金属材料会出现损伤。转子的损伤主要包括低周疲劳损伤和高温蠕变损伤，低周疲劳损伤和高温蠕变损伤对转子产生的作用是耦合的，在机组实际运行中有必要考虑到两类损伤的耦合作用。金属材料的蠕变损伤主要是指长期高温环境下，塑性应变逐渐变大的现象。转子金属材料一般是低碳合金钢，当汽轮机机组进入平稳状态后，各类热工参数逐渐稳定，热应力变化较小，其损伤形式由低周疲劳损伤过渡到高温蠕变损伤。危害主要体现在，由于金属蠕变损伤的累积，转子部件发生过量的塑性变形而无法使用；或者蠕变进入到加速发展阶段，材料会发生蠕变破裂，导致严重的安全事故。因此，对于汽轮机机组，尤其是最昂贵的转子部件，需要进行连续的、实时的、高精度的蠕变损伤计算。
现有技术中，既有单独计算高温蠕变损伤的方法，也有考虑疲劳-蠕变耦合情况的方法，例如非线性力学方法、数值计算法等，这类方法理论基础好，计算精度高，但是算法复杂，计算速度慢。另外一类是解析法，这类算法强调计算速度，对损伤计算过程和模型做了一定的简化，而牺牲了计算精度。这两类方法对热能装备的损伤计算都具有积极的指导作用。但是，如何将二者相结合，既能满足热能装备的蠕变损伤实时计算，又可以保证计算精度，是现在急需解决的一个问题。
发明内容
鉴于以上所述现有技术的缺点，本发明的目的在于提供一种基于有限元模型的汽轮机转子蠕变损伤实时计算方法，用于解决现有技术对汽轮机转子蠕变损伤计算中实时性与高精度难以共融的问题。
为实现上述目的及其他相关目的，提供一种基于有限元模型的汽轮机转子蠕变损伤实时计算方法，其特征在于，所述方法包括：通过有限元分析软件计算转子的温度场和应力场；提取关键部位的温度和应力数据；将蠕变损伤计算相关的温度和应力数据进行归一化处理；采用多元回归方法建立蠕变应力实时计算模型；解析并拟合应力-损伤函数关系，实时计算蠕变损伤；建立蠕变损伤实时计算系统结构。
优选地，通过有限元分析软件计算汽轮机转子的温度场和应力场。
(1)温度场计算可采用二维轴对称模型，考虑温度函数为T(z，r，t)，则有：
转子外叶轮外表面与蒸汽对流换热，为第三类边界条件：
其中，r为转子径向参量；r0为表面半径；h为蒸汽与转子换热系数；Tf为蒸汽温度。转子与轴承润滑油接触部位可认为是第一类边界条件，即温度与润滑油温相等:
转子两端忽略转子热流密度。边界条件表示如下；
其中，Γ表示转子中心及两端边界，q为热流密度。
(2)应力场计算时对转子固定点部位进行零位移约束具体为：
u|z＝0＝0
优选地，提取关键部位的温度和应力数据。一般情况下转子调节级根部是应力较大的部位，可提取调节级根部应力最大处的温度和应力曲线，再通过数据导出模块将所需的温度数据和应力数据分别导出保存。
优选地，将蠕变损伤计算相关的温度、压力、转速以及应力数据进行归一化处理。由于不同评价指标往往具有不同的量纲和量纲单位，这样的情况会影响到数据分析的结果。为了消除指标之间的量纲影响，数据集需要进行数据标准化处理，具体为：
其中max为原始样本数据的最大值，min为原始样本数据的最小值，x*为归一化后的数值，x为归一化前的数值。
优选地，采用多元回归方法建立蠕变应力实时计算模型。多元线性回归模型：
y＝b0+b1x1+L+bmxm+ε ε～N(0,σ2)
m元线性回归方程估计值为：
目标函数为最小化残差平方和：
其中x为自变量；y为因变量；b0～bm为回归系数；ε为随机变量；N为满足期望为0，方差为σ2的正态分布；X为自变量x组成的矩阵；Y为n个因变量y组成的矩阵；XT为X的转置矩阵；为b0～bm的估计值；为y的估计值；为yi的估计值；σ2为方差；为σ2的无偏估计值；Qe为残差平方和。
优选地，解析并拟合应力-损伤函数关系，实时计算蠕变损伤。已知最大应力，根据以下三式计算高温蠕变损伤。
(1)循环应力-应变关系式为：
其中E是弹性模量，K'是循环强度系数，n'是循环应变硬化指数。
(2)蠕变损伤发生在平稳状态，此时主蒸汽温度为540℃的时候，应变与循环致裂周次Nf关系为：
ε＝0.00332(2Nf)-0.0697+0.6264(2Nf)-0.7553，
(3)损伤d和应变ε的关系为：
其中的多项式拟合系数P1＝9.7904×1012，P2＝-3.8943×1011， P3＝6.011×109，P4＝-4.6037×107，P5＝1.8292×105，P6＝-3.3086×102， P7＝0.27029，P8＝0。
优选地，建立蠕变损伤实时计算系统架构。将训练温度输入到ADINA有限元软件中来求得热应力数据，再将压力和分析得到的离心力和求得的热应力一起输入到多元回归模型中，同时将训练的温度、压力、转速以及实测的温度、压力、转速也输入到多元回归模型中进行拟合对比分析，将拟合后得到的应力进行应力-应变的关系式得到应变，再通过应变-损伤的关系模型进行多项式拟合，计算得到蠕变损伤。
附图说明
图1显示为一种基于有限元模型的汽轮机转子蠕变损伤实时计算方法流程示意图。
图2显示为某国产300MW汽轮机转子冷态启动温度曲线。
图3显示为冷态启动转子蠕变应力曲线。
图4显示为冷态启动蠕变损伤应力-应变曲线。
图5显示为蠕变应力计算训练数据曲线。
图6显示为蠕变应力计算测试数据曲线。
图7显示为高温蠕变损伤实时计算曲线图。
图8显示为一种基于有限元模型的汽轮机转子蠕变损伤实时计算系统架构。
具体实施方式
以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易的了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。
请参阅图1至图8。需要说明的是，本实施例中所提供的图示仅以示意方式说明本发明的基本构想，遂图示中仅显示与本发明中有关的组件而非按照实际实施时的组件数目、形状及尺寸绘制，其实际实施时各组件的形态、数量及比例可为一种随意的改变，且其组件布局形态也可能更为复杂。
在汽轮机转子高温蠕变损伤计算领域，计算方法主要可以分为理论分析法和解析法。前者理论基础好，计算精度高，但是算法复杂，计算速度慢。后者强调计算速度，对损伤计算过程和模型做了简化，牺牲了计算精度。目前需要寻求一种可以具备一定精度的转子高温蠕变损伤实时计算方法，本发明就是基于这些理念而形成的。
本发明的目的在于提供一种转子蠕变损伤实时计算方法，用以解决现有方法中转子蠕变损伤实时计算和高精度计算难以共融的问题。以下将详细描述本发明的一种基于有限元模型的汽轮机转子蠕变损伤实时计算方法的原理和实施方式，使本领域技术人员不需要创造性劳动即可理解本发明的一种基于有限元模型的汽轮机转子蠕变损伤实时计算方法。
如图1所示，本发明一种基于有限元模型的汽轮机转子蠕变损伤实时计算方法，所述方法步骤包括：
S1，通过有限元分析软件计算转子的温度场和应力场；
S2，提取关键部位的温度和应力数据；
S3，将蠕变损伤计算相关的温度、压力、转速以及应力数据进行归一化处理；
S4，采用多元回归方法建立蠕变应力实时计算模型；
S5，解析并拟合应力-损伤函数关系，实时计算蠕变损伤；
S6，建立蠕变损伤实时计算系统架构。
下面结合具体实施方式对本发明进行详细说明。本发明实施在ADINA 8.5 软件和MATLAB R2016a环境下完成。具体方法如下：使用的汽轮机转子是以某电厂的1号机组为原型，该汽轮机是单轴、凝汽式汽轮机。转子的材料为30Cr1Mo1V。高压转子的长度约为4800mm，具有单列调节级和11级的压力级。通过ADINA建立该转子二维模型，根据《国产300MW汽轮机运行导则》，在调节级和高压级施加工况模拟。
首先执行步骤S1，使用有限元分析软件计算汽轮机转子温度场和应力场。转子呈对称形状，建立适当的几何模型，同时施加荷载、确定边界条件及材料参数，通过有限元分析软件进行分析求解。
在步骤S2中，提取关键部位的温度和应力数据。根据汽轮机转子的温度场和应力场云图可以发现，转子前几级温度较大，且温度梯度较为明显；从应力场云图也可看到，转子调节级根部是应力较大的部位。提取调节级根部应力最大处的温度曲线和应力曲线以及温度和应力的等值线，再通过数据导出模块将所需的温度数据、疲劳损伤实测应力数据、蠕变损伤实测应力数据和蠕变损伤应力- 应变数据导出为Excel文档并保存。某国产300MW汽轮机转子冷态启动温度曲线如图2所示，冷态启动转子蠕变应力曲线如图3所示，冷态启动蠕变损伤应力-应变曲线如图4所示。
在步骤S3中，将蠕变损伤计算相关的温度、压力、转速以及应力数据进行归一化处理。原始数据经过数据归一化处理后，各指标处于同一数量级，本发明采用Min-Max标准化对原始数据的线性变换，使结果值映射到[min,max]之间。转换函数为：在MATLAB中调用mapminmax函数，对温度、压力、转速以及应力数据进行归一化处理。
在步骤S4中，采用多元回归方法建立蠕变应力实时计算模型，根据多元回归模型y＝b0+b1x1+L+bmxm+ε ε～N(0,σ2)，采用最小二乘法和求偏导数，整理后得到方程组，引入矩阵：
则方程组可写成：XTXB＝XTY，于是可解得：
的估计值：
方程为m元线性回归方程。在MATLAB软件中调用regress多元回归函数、mse均方误差函数对预测数据和实测数据进行拟合分析。蠕变应力计算训练数据图像如图5所示，蠕变应力计算测试数据图像如图6所示。
在步骤S5中，解析并拟合应力-损伤函数关系，实时计算蠕变损伤。如果已知最大应力值，那么可以通过过程量的数值计算得到高温蠕变损伤。循环应力-应变关系式为：其中E是弹性模量，K'是循环强度系数，n'是循环应变硬化指数。最大应力发生在主蒸汽温度达到540℃的时候，此时汽轮机转子高温蠕变损伤可以通过下式计算得到：ε＝0.00332(2Nf)-0.0697+0.6264(2Nf)-0.7553，循环致裂次数Nf与应变ε的关系可以通过多项式拟合得到，并可以进一步得到以金属材料应变作为输入，以循环致裂次数Nf作为输出的关系式。因此，通过多项式拟合重新建立起应力-损伤关系模型，解决了损伤的实时计算问题。其中损伤d和应变ε的关系可通过多项式拟合得到，高温蠕变损伤实时计算曲线如图7所示。
在步骤S6中，建立蠕变损伤实时计算系统架构。将训练温度输入到ADINA有限元软件中来求得热应力数据，再将压力和分析得到的离心力和求得的热应力一起输入到多元回归模型中，同时将训练的温度、压力、转速以及实测的温度、压力、转速也输入到多元回归模型中进行拟合对比分析，将拟合后得到的应力进行应力-应变的关系式得到应变，再通过应变-损伤的关系模型进行多项式拟合，从中即可得到蠕变损伤，如图8所示。

Claims (6)

1.一种基于有限元模型的汽轮机转子蠕变损伤实时计算方法，其特征在于，所述方法包括：
通过有限元分析软件计算转子的温度场和应力场；
提取关键部位的温度和应力数据；
将蠕变损伤计算相关的温度、压力、转速以及应力数据进行归一化处理；
采用多元回归方法建立蠕变应力实时计算模型；
解析并拟合应力-损伤函数关系，实时计算蠕变损伤；
建立蠕变损伤实时计算系统结构。
2.根据权利要求1所述的转子蠕变损伤实时计算方法，其特征在于，通过有限元分析软件计算汽轮机转子的温度场和应力场，边界条件具体为：
(1)转子外叶轮外表面与蒸汽对流换热，为第三类边界条件：
其中，r为转子径向参量；r0为表面半径；h为蒸汽与转子换热系数；Tf为蒸汽温度。
(2)转子与轴承润滑油接触部位可认为是第一类边界条件，即温度与润滑油温相等:
(3)转子两端忽略转子热流密度，边界条件表示如下：
其中，Γ表示转子中心及两端边界，q为热流密度。
(4)应力场计算时对转子固定点部位进行零位移约束具体为：
u|z＝0＝0。
3.根据权利要求1所述的转子蠕变损伤实时计算方法，其特征在于，将蠕变损伤计算相关的温度、压力、转速以及应力数据进行归一化处理具体为：
其中max为原始样本数据的最大值，min为原始样本数据的最小值，x*为归一化后的数值，x为归一化前的数值。
4.根据权利要求1所述的转子蠕变损伤实时计算方法，其特征在于，采用多元回归方法建立蠕变应力实时计算模型。
(1)多元线性回归模型具体为：
y＝b0+b1x1+L+bmxm+ε ε～N(0,σ2)，
(2)m元线性回归方程估计值为：
(3)目标函数为最小化残差平方和：
其中x为自变量；y为因变量；b0～bm为回归系数；ε为随机变量；N为满足期望为0，方差为σ2的正态分布；X为自变量x组成的矩阵；Y为n个因变量y组成的矩阵；XT为X的转置矩阵；为b0～bm的估计值；为y的估计值；为yi的估计值；σ2为方差；为σ2的无偏估计值；Qe为残差平方和。
5.根据权利要求1所述的转子蠕变损伤实时计算方法，其特征在于，解析并拟合应力-损伤函数关系，实时计算蠕变损伤。
(1)循环应力-应变关系式为：
其中E是弹性模量，K'是循环强度系数，n'是循环应变硬化指数。
(2)主蒸汽温度为540℃的时候，应变与循环致裂周次Nf关系为：
ε＝0.00332(2Nf)-0.0697+0.6264(2Nf)-0.7553，
(3)损伤d和应变ε的关系为：
其中的多项式拟合系数P1＝9.7904×1012，P2＝-3.8943×1011，P3＝6.011×109，P4＝-4.6037×107，P5＝1.8292×105，P6＝-3.3086×102，P7＝0.27029，P8＝0。
6.根据权利要求1所述的转子蠕变损伤实时计算方法，其特征在于，建立蠕变损伤实时计算系统架构，具体如图8所示。
CN201811424356.2A 2018-11-27 2018-11-27 一种基于有限元模型的汽轮机转子蠕变损伤实时计算方法 Active CN109522663B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811424356.2A CN109522663B (zh) 2018-11-27 2018-11-27 一种基于有限元模型的汽轮机转子蠕变损伤实时计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811424356.2A CN109522663B (zh) 2018-11-27 2018-11-27 一种基于有限元模型的汽轮机转子蠕变损伤实时计算方法

Publications (2)

Publication Number Publication Date
CN109522663A true CN109522663A (zh) 2019-03-26
CN109522663B CN109522663B (zh) 2023-05-02

Family

ID=65793127

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811424356.2A Active CN109522663B (zh) 2018-11-27 2018-11-27 一种基于有限元模型的汽轮机转子蠕变损伤实时计算方法

Country Status (1)

Country Link
CN (1) CN109522663B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110096769A (zh) * 2019-04-16 2019-08-06 西北工业大学 一种发动机单晶叶片热机械疲劳寿命快速评估方法
CN112613148A (zh) * 2020-12-30 2021-04-06 一重集团大连工程技术有限公司 一种基于数值分析变形数据的核电装备设计方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2713160A2 (en) * 2012-10-01 2014-04-02 Hitachi Ltd. Method and system for evaluating creep damage of high temperature component
CN103940662A (zh) * 2013-01-18 2014-07-23 安阳工学院 高温材料应力松弛剩余应力和损伤的预测方法
CN107145687A (zh) * 2017-06-07 2017-09-08 上海电力学院 汽轮机转子启动曲线优化及蠕变‑疲劳寿命评估的方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2713160A2 (en) * 2012-10-01 2014-04-02 Hitachi Ltd. Method and system for evaluating creep damage of high temperature component
CN103940662A (zh) * 2013-01-18 2014-07-23 安阳工学院 高温材料应力松弛剩余应力和损伤的预测方法
CN107145687A (zh) * 2017-06-07 2017-09-08 上海电力学院 汽轮机转子启动曲线优化及蠕变‑疲劳寿命评估的方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
王立滨等: "损伤力学方法在汽轮机高温蠕变损伤分析中的应用", 《汽轮机技术》 *
荆建平等: "汽轮机转子疲劳-蠕变损伤的非线性损伤力学分析", 《中国电机工程学报》 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110096769A (zh) * 2019-04-16 2019-08-06 西北工业大学 一种发动机单晶叶片热机械疲劳寿命快速评估方法
CN112613148A (zh) * 2020-12-30 2021-04-06 一重集团大连工程技术有限公司 一种基于数值分析变形数据的核电装备设计方法
CN112613148B (zh) * 2020-12-30 2024-05-28 一重集团大连工程技术有限公司 一种基于数值分析变形数据的核电装备设计方法

Also Published As

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

Similar Documents

Publication Publication Date Title
Marconcini et al. Large eddy simulation and RANS analysis of the end-wall flow in a linear low-pressure-turbine cascade—Part II: loss generation
CN102914427B (zh) 一种多轴随机载荷下疲劳损伤评估方法
CN101750216B (zh) 次同步振荡引起的汽轮发电机轴系疲劳损伤在线分析方法
KR101474102B1 (ko) 임펠러 형상 최적설계방법
CN109522663A (zh) 一种基于有限元模型的汽轮机转子蠕变损伤实时计算方法
Trampert et al. Thermomechanical fatigue life prediction of cylinder heads in combustion engines
CN103984813A (zh) 一种离心压缩机裂纹叶轮结构的振动建模与分析方法
CN108629092B (zh) 一种基于尺寸效应修正的涡轮盘分区可靠性分析方法
CN107729658A (zh) 一种可靠性优化设计的模糊智能分解协调极值响应面法
CN101201307A (zh) 一种自动绘制材料热加工图的方法
CN113128047A (zh) 一种航空发动机推力计算方法、存储介质及计算设备
Yao et al. Vibration characteristics of pre-twisted rotating Ti-SiC composite airfoil blade
Mu et al. Study on reconstruction and prediction methods of pressure field on blade surfaces for oil-filling process in a hydrodynamic retarder
Prasad et al. Static and dynamic analysis of a centrifugal pump impeller
Xu et al. Random buckling bearing capacity of super-large cooling towers considering stochastic material properties and wind loads
CN109558669A (zh) 一种基于有限元模型的汽轮机转子疲劳损伤在线计算方法
Ogbereyivwe et al. Families of Means-Based Modified Newtons Method for Solving Nonlinear Models
CN118036389A (zh) 预测航空发动机叶片剩余寿命的方法及航空发动机
GB2555010B (en) Determining loads on a wind turbine
CN117669254A (zh) 一种分析离心泵叶片误差影响流场的方法
Wang et al. Real-time variable geometry triaxial gas turbine model for hardware-in-the-loop simulation experiments
Xiao et al. Prediction of separation flows around a 6: 1 prolate spheroid using RANS/LES hybrid approaches
CN106021692B (zh) 一种蒸汽轮机性能与可靠性一体化设计方法
CN106548029A (zh) 韦伯燃烧规则经验系数的确定方法
Shlyannikov et al. Method for assessment of the residual life of turbine disks

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