CN112380626B - 一种零件间接触应力及其分布状态的计算方法 - Google Patents
一种零件间接触应力及其分布状态的计算方法 Download PDFInfo
- Publication number
- CN112380626B CN112380626B CN202011316830.7A CN202011316830A CN112380626B CN 112380626 B CN112380626 B CN 112380626B CN 202011316830 A CN202011316830 A CN 202011316830A CN 112380626 B CN112380626 B CN 112380626B
- Authority
- CN
- China
- Prior art keywords
- stress
- strain
- contact surface
- contact
- epsilon
- 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
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/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/20—Design optimisation, verification or simulation
- G06F30/25—Design optimisation, verification or simulation using particle-based methods
-
- 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
Abstract
本发明属于力学计算技术领域,公开了一种零件间接触应力及其分布状态的计算方法。为了计算出接触区域的接触应力和应力分布状态,通过实测得到接触区域相关周边密切非接触区位置的应变值,按照应力应变规律,结合有限元仿真逆向推算出接触区域的接触应力和应力分布状态。
Description
技术领域
本发明属于力学计算技术领域,特别涉及一种零件间接触应力及其分布状态的计算方法。
背景技术
飞机组件、部件产品的装配过程,几何精度保证了产品的正确可靠的装配位置关系,然而产品的装配应力却深刻影响着产品的性能和使用寿命,零件之间的挤压变形,强迫装配现象尤为严重,飞机上机构的某一状态位置对应的接触应力影响着产品的使用性能。因此分析零件的接触应力及其在接触面上的分布状态具有重要意义。
目前存在的测量接触应力的方法或装置有两类。
接触式的测量,此类测量方法需要将传感器塞到零件与零件之间的接触面上,依靠测力计的变形与应力的转换关系直接得出接触应力的大小。例如目前普遍使用的薄膜应力传感器。此种测量方法需要将其直接塞至接触面的位置,针对装配产品,零件与零件之间的接触面不可达的特性,接触式测量的适用范围极其有限。
非接触式的测量,非接触式测量可以满足零件接触应力不可达的特点。目前使用的方法是X衍射法和雷达超声测量法,然而对于非透明件或有一定厚度的零件,该测量方法无法实现接触应力的有效测量。
发明内容
本发明的目的是给出一种非接触式确定零件间接触应力大小分布状态的方法。为了计算出接触区域的接触应力和应力分布状态,通过实测得到接触区域相关周边密切非接触区位置的应变值,按照应力应变规律,结合有限元仿真逆向推算出接触区域的接触应力和应力分布状态。
为实现该目的,本发明采用了下面的技术方案。一种计算零件间接触应力及其分布状态的方法,步骤如下:
步骤1测量待测零件几何外形,关键特征,获取零件的真实几何参数;
步骤2根据实测数据建立有限元模型,对于不影响受力分析的部件,采用边界条件替代,在分析零件的接触面上施加均布的接触应力σ;后处理获取关键部位的应变值c=(c1,c2,c3...cm);
这里m指的是后处理时获取的应变源于多个位置和多个方向。
验证有限元模型的准确定性,在实际零件的实际接触面上施加均布的接触应力σ,在有限元后处理同样的位置上贴应变片测得实际应变值
∈′=(ε′1,ε′2,ε′3...ε′m)
步骤3根据实际经验,结合实际工况,按等梯度法,以固定大小δ为间隔,改变接触应力的均布载荷σ+k·δ(k=0,1,2...)的大小,得到多组σ-∈的关系数据。
步骤4装配前,贴附电阻应变片于待测零件上,应变片的贴附位置,布局的方向与步骤2中有限元仿真的后处理阶段测取的位置,方向保持一致。贴附完毕后,测得对应各处的应变初值
∈0=(ε01,ε02,ε03...ε0m)
上式中ε01,ε02,ε03...ε0m与步骤2所给应变ε1,ε2,ε3...εm对应。
步骤5装配完成后,再次读取应变仪上的示数,得到各处的应变值
∈1=(ε11,ε12,ε13...ε1m)
上式中ε11,ε12,ε13...ε1m与步骤2所给应变ε1,ε2,ε3...εm对应。
步骤6上述∈,∈0,∈1数据传入后台算法,得到接触面的平均应力大小,本步骤由计算机实现。
具体的:
步骤6.1处理步骤3获取的σ-∈的对应数据需进一步处理,对于每次模拟的应变值∈′=(ε′1,ε′2ε′3...ε′m)的应变分量ε′1,ε′2ε′3...ε′m,取
以σ为横坐标,以k组应变值ε′为纵坐标得到σ-ε′关系对应的散点图;通过最小二乘法拟合得到σ-ε′函数关系ε′=f(σ);
步骤6.2处理步骤4、5获取的数据,得到实测应变值。取
步骤7后台程序计算接触面上的接触应力分布状态,指的是接触应力在接触面上的大小随位置的不同,具体值是多少。本部分采用参数化建模,基于python语言对ABAQUS软件进行二次开发,建立零件接触状态有限元模型与后处理过程,以便不断改变单元的载荷,进行大量的循环计算,提高分析效率。同时通过MATLAB更新粒子的速度和位置,并将位置信息即施加的接触面应力分布写入Python的脚本语言中,以批处理方式调用ABAQUS分析求解,并读取结果文件,进行适应度评价。
具体的:
步骤7.1将待测接触面等分为n等份的矩形区域,分析实际工况,每个单元的受应力合力方向约定为固定的某一方向。
步骤7.2设置终止条件为最大迭代次数max,粒子规模为M,单个粒子为n维的向量
x=(σ1,σ2,σ3...σn)T
为接触面上应力分布的一个可行解,向量元素σ1,σ2,σ3...σn为接触面n个矩形区域的应力值。
设置粒子的最大速度Vmax。
步骤7.3设置个体的适应度函数。
fitness(x)=||∈1-∈0-∈i||
其中,∈j=(εi1,εi2,εj3...εim)表示第i次迭代求解的n个应变监测目标。∈1-∈0表示实测的m个应变监测值组成的向量。
步骤7.4初始化速度和位置。
步骤7.5传参,提交零件接触状态有限元模型分析计算,得到每个粒子初始设置下全局最优解Pg。
步骤7.6更新粒子的速度和位置
vid=ωVid+C1random(0,1)(Pid-Xid)+C2random(0,1)(Pgd-Xid)
Xid=Xid+Vid
其中,ω为惯性因子其值为非负,通过调ω的大小,可以对全局寻优性能和局部寻优性能进行调整。C1和C2为加速常数,前者为每个粒子的个体学习因子,后者为每个粒子的社会学习因子。Suganthan的实验表明:C1和C2为常数时可以得到较好的解,一般取C1=C2∈[0,4]。random(0,1)表示区间[0,1]上的随机数,Pid表示第i个变量的个体极值的第d维,Pgd表示全局最优解的第d维。
步骤7.7提交零件接触状态有限元模型计算。一次迭代结束后,根据适应度值得到粒子的个体最优解Pi,和全局最优解Pg。
Pi表示以此循环结束后第i个个体的历史最优解。Pg表示一次循环结束后,全体的最优解。
步骤7.8判断是否满足终止条件迭代次数达到max,若满足,则结束计算。否则跳到步骤7.6,继续循环计算。
步骤7.9输出最优解
best_x=(σ1,σ2,σ3...σn)T,即确定了接触面的应力分布状态。
有益效果:
1)本发明通过传感器得到实测变形,将变形结果作为输入参数根据仿真计算出目标区域的装配应力,解决了不可接触区域无法测量接触应力或测量误差较大的问题;
2)该计算方法通用性强,设备成本低;
3)测量接触应力便于操作员参数调整,增加装配产品的寿命。
附图说明
图1为零件间接触应力及其分布状态的计算方法流程图;
图2为待测试零件轴测图;
图3为待测试零件及其接触应力的位置示意图;
图4为有限元约束和载荷示意图;
图5为应变片布局方式示意图;
图6为接触面区域分块示意图;
图中:1-第一连接件,2-第二连接件,3-外部围框。
具体实施方式
下面结合附图和具体实例对本发明作进一步说明。
图1为本发明实例的确定零件间接触应力及其分布状态的流程图。
实例涉及求解的是图2所示的两零件之间的接触应力大小及其分布状态,第一连接件1和第一连接件2分别通过以四个螺栓连接固定在外部围框3上,根据本实际工况,由于零件制造的误差,第一第一连接件1和第二第一连接件2接触面之间存在相互挤压,作用在零件接触面上的接触应力的方向为X方向。如图3所示,第一连接件1与第一连接件2之间的接触面为单曲率的曲面,其接触面在YOZ平面上的投影所得面积大小为20mm×20mm。
本实例测量和计算操作步骤为:
S1扫描第一连接件1、第二连接件2的接触面,第一连接件1的长度、连接耳片等关键特征。得到真实几何参数。
S2通过得到真实数据,建立零件接触状态有限元模型。以边界约束代替外部围框3,如图4所示,本实例求解第一连接件1、第一连接件2之间的接触应力,建立第一连接件1的有限元模型,以作用在接触面上的边界条件即接触应力代替第一连接件2的有限元模型。
∈=(ε1,ε2,ε3...ε8)
∈′=(ε′1,ε′2,ε′3...ε′m)
S4采用参数化建模,基于python语言通过ABAQUS软件建立零件接触状态有限元模型,以便不断改变单元的载荷,进行大量的循环计算,改变均布接触应力分别为时进行有限元计算,同S3处理方式,得到20组数据。
S5装配前,贴附电阻应变片于待测第一连接件1上,应变片的贴附位置,布局的方向等按图示5方式,贴附完毕后,测得对应各处的应变初值
∈0=(ε01,ε02,ε03...ε08)
上式中ε01,ε02,ε03...ε08与S3所给应变ε1,ε2,ε3...ε8对应。
S6装配完成后,再次读取应变仪上的示数,得到各处的应变值
∈1=(ε11,ε12,ε13...ε18)
所以实测的应变值为∈1-∈0。
具体的,
以σ为横坐标,以ε为纵坐标将S4有限元模拟的数据画出,得到σ-ε关系对应的21个点构成的散点图。再通过最小二乘拟合得到σ-ε函数关系
ε=f(σ)
S7.2处理步骤S5、S6获取的数据,得到实测应变值。取
S8以第一连接件1接触面投影在YOZ面上的面积相等为原则,将接触面等分为10*10等份的矩形单元,如图6所示每个单元的受应力合力方向固定为X向。
S9参数初始化,设置终止条件为最大迭代次数max=50,粒子规模为M=20,单个粒子为n=100维的向量。
x=(σ1,σ2,σ3...σ100)T
为接触面上应力分布的一个可行解,向量元素σ1,σ2,σ3...c100为接触面n=100个矩形单元的应力值。
设置粒子的最大速度Vmax,随即初始化20个粒子的速度和位置;
设置适应度函数为
fitness(x)=||∈1-∈0-∈i||2
其中,∈i=(εi1,εi2,εi3...εi8)表示第i次迭代求解的8个应变监测目标。∈1-∈0表示实测的8个应变监测值。
S10传参,提交ABAQUS零件接触状态有限元模型分析计算,得到每个粒子初始设置下全局最优解Pg。
S11更新粒子的速度和位置
vid=ωVid+C1random(0,1)(Pid-Xid)+C2random(0,1)(Pgd-Xid)
Xid=Xid+Vid
其中,通过调ω的大小,可以对全局寻优性能和局部寻优性能进行调整,本实例设置ω=1;设置加速常数C1=C2=1.5。
S12提交ABAQUS零件接触状态有限元模型计算。一次迭代结束后,根据适应度值得到粒子的个体最优解Pi,和全局最优解Pg。
S13判断是否满足终止条件迭代次数达到max=50次,若满足,则结束计算。否则跳到S11,继续循环计算。
S14输出最优解
best_x-(σ1,σ2,σ3...σ100)T。
Claims (6)
1.一种零件间接触应力及其分布状态的确定方法,其特征在于:所述方法包括以下步骤:
步骤一:测量待测零件几何外形,得到待测零件的几何参数;
步骤二:根据待测零件的几何参数建立有限元模型;并在待测零件的接触面上施加均布的接触应力σ;基于有限元模型处理后得到待测零件接触面边缘上等间隔的m个位置点的一组应变值∈=(ε1,ε2,ε3...εm);
步骤三:按照等梯度法,以固定值δ为间隔,改变接触面上施加均布的接触应力的大小σ+k·δ(k=0,1,2…),得到k组σ-ε的关系数据;
步骤四:在待测零件装配前,在待测零件上贴附电阻应变片;电阻应变片数量为m个,贴附位置为所述步骤二中待测零件接触面边缘上等间隔的m个位置点,并测量各处的应变初值ε0=(ε01,ε02,ε03...ε0m);
步骤五:将待测零件装配完成后,再次读取电阻应变片的值,得到各处的应变值ε1=(ε11,ε12,ε13...ε1m);
步骤六:根据上述的∈、∈0和∈1计算接触面的平均应力大小;
步骤七:基于粒子群算法根据上述的∈、∈0、∈1和接触面的平均应力大小计算接触面的应力分布状态;具体过程如下:
7.1将待测接触面等分为n等份的矩形区域;
7.2设置终止条件为最大迭代次数max,粒子规模为M,单个粒子为n维的向量x=(σ1,σ2,σ3…σn)T;
向量x为接触面上应力分布的一个可行解,向量元素σ1,σ2,σ3…σn为接触面n个矩形区域的应力值;
基于接触面内平均应力设置位置空间
7.3设置粒子的最大速度Vmax;
设置个体的适应度函数fltness(x)=||∈1-∈0-∈j||;
其中,∈j=(εi1,εi2,εi3...εim)表示第i次迭代求解的m个应变监测目标,∈1-∈0表示实测的m个应变片实测应变值组成的向量;
7.4初始化粒子速度和位置;
7.5通过ABAQUS计算每个粒子初始设置下全局最优解Pg;
7.6更新粒子的速度和位置
vid=ωVid+C1random(0,1)(Pid-Xid)+C2random(0,1)(Pgd-Xid),
Xid=Xid+Vid;
其中,ω为惯性因子,其值为非负;
C1和C2为加速常数,C1为每个粒子的个体学习因子,C2为每个粒子的社会学习因子;random(0,1)表示区间[0,1]上的随机数,Pid表示第i个变量的个体极值的第d维,Pgd表示全局最优解的第d维;
7.7利用更新后的粒子速度和位置提交ABAQUS进行计算,根据适应度值得到粒子的个体最优解Pi,和全局最优解Pg;Pi表示以此循环结束后第i个个体的历史最优解;Pg表示一次循环结束后,全体的最优解;
7.8判断是否满足终止条件迭代次数达到max,若满足,则结束计算;否则重复步骤7.6-7.8直至循环结束;
7.9输出最优解best_x=(σ1,σ2,σ3…σn)T,最优解为接触面的应力分布状态。
3.根据权利要求2所述的一种零件间接触应力及其分布状态的确定方法,其特征在于:n等份的矩形区域内的受应力合力方向相同。
4.根据权利要求3所述的一种零件间接触应力及其分布状态的确定方法,其特征在于:适应度值为每次迭代后的应变监测目标组成的向量与实测应变片测量得到的应变值组成的向量之差构成的二范数。
5.根据权利要求4所述的一种零件间接触应力及其分布状态的确定方法,其特征在于:取C1=C2∈[0,4]。
6.根据权利要求5所述的一种零件间接触应力及其分布状态的确定方法,其特征在于:m个应变监测目标位置为m个应变片贴附位置。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202011316830.7A CN112380626B (zh) | 2020-11-20 | 2020-11-20 | 一种零件间接触应力及其分布状态的计算方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202011316830.7A CN112380626B (zh) | 2020-11-20 | 2020-11-20 | 一种零件间接触应力及其分布状态的计算方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN112380626A CN112380626A (zh) | 2021-02-19 |
CN112380626B true CN112380626B (zh) | 2022-09-06 |
Family
ID=74587544
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202011316830.7A Active CN112380626B (zh) | 2020-11-20 | 2020-11-20 | 一种零件间接触应力及其分布状态的计算方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN112380626B (zh) |
Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102184273A (zh) * | 2011-02-18 | 2011-09-14 | 洛阳轴研科技股份有限公司 | 斜撑式离合器楔块表面应力的有限元模型建立及修正方法 |
CN102759505A (zh) * | 2012-07-25 | 2012-10-31 | 上海交通大学 | 中厚板材料压缩试验辅助装置及流动应力曲线测定方法 |
CN106650034A (zh) * | 2016-11-30 | 2017-05-10 | 大工(青岛)新能源材料技术研究院有限公司 | 一种基于有限元分析超导股线接触仿真方法 |
CN106709142A (zh) * | 2016-11-18 | 2017-05-24 | 大连理工大学 | 一种获取螺栓连接结合面应力分布的方法 |
CN106907534A (zh) * | 2017-01-10 | 2017-06-30 | 海隆石油工业集团有限公司 | 一种多层非粘结柔性管抗拉性能快速评价方法 |
CN108256160A (zh) * | 2017-12-21 | 2018-07-06 | 中国石油天然气集团公司 | 一种热采工况特殊螺纹接头密封接触压应力的预测方法 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20160357893A1 (en) * | 2016-08-15 | 2016-12-08 | Xianwu Ling | Contact stiffness estimation based on structural frequency responses |
-
2020
- 2020-11-20 CN CN202011316830.7A patent/CN112380626B/zh active Active
Patent Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102184273A (zh) * | 2011-02-18 | 2011-09-14 | 洛阳轴研科技股份有限公司 | 斜撑式离合器楔块表面应力的有限元模型建立及修正方法 |
CN102759505A (zh) * | 2012-07-25 | 2012-10-31 | 上海交通大学 | 中厚板材料压缩试验辅助装置及流动应力曲线测定方法 |
CN106709142A (zh) * | 2016-11-18 | 2017-05-24 | 大连理工大学 | 一种获取螺栓连接结合面应力分布的方法 |
CN106650034A (zh) * | 2016-11-30 | 2017-05-10 | 大工(青岛)新能源材料技术研究院有限公司 | 一种基于有限元分析超导股线接触仿真方法 |
CN106907534A (zh) * | 2017-01-10 | 2017-06-30 | 海隆石油工业集团有限公司 | 一种多层非粘结柔性管抗拉性能快速评价方法 |
CN108256160A (zh) * | 2017-12-21 | 2018-07-06 | 中国石油天然气集团公司 | 一种热采工况特殊螺纹接头密封接触压应力的预测方法 |
Non-Patent Citations (4)
Title |
---|
"Contact stress analysis and calculation of stress concentration factors at the tool joint of a drill pipe";ShahaniA R等;《Materials & Design》;20091031;第30卷(第9期);第3615-3621页 * |
"Finite element analysis on the stress state of rope hoisting equipment based on the ABAQUS,";W.Du 等;《2017 8th International Conference on Mechanical and Intelligent Manufacturing Technologies (ICMIMT)》;20170630;第84-88页 * |
"滚滑轴承的力学特性及疲劳寿命分析";曾国文;《中国优秀硕士学位论文全文数据库工程科技Ⅱ辑》;20170215;第C029-216页 * |
"零件真实粗糙表面构建及微观接触性能分析";姜英杰 等;《机械设计与制造》;20180808(第8期);第8-10页 * |
Also Published As
Publication number | Publication date |
---|---|
CN112380626A (zh) | 2021-02-19 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN111546328B (zh) | 基于三维视觉测量的手眼标定方法 | |
CN112613118B (zh) | 火箭发动机内部不可测装配质量数字孪生建模与追溯方法 | |
Tian et al. | Determination of optimal samples for robot calibration based on error similarity | |
CN110285781B (zh) | 一种相对于基准面的平面平行度快速评定方法 | |
Stone et al. | Statistical performance evaluation of the S-model arm signature identification technique | |
CN112926152B (zh) | 一种数字孪生驱动的薄壁件装夹力精准控制与优化方法 | |
CN111540001B (zh) | 航空发动机涡轮叶片气膜孔轴线方向检测方法 | |
CN108846200B (zh) | 一种基于迭代法的准静态桥梁影响线识别方法 | |
CN111232239B (zh) | 曲面挠变位移场重构方法、装置及设备 | |
CN112380626B (zh) | 一种零件间接触应力及其分布状态的计算方法 | |
CN114611420A (zh) | 非定常气动力计算精度评估及修正方法 | |
Müller et al. | Utilization of single point uncertainties for geometry element regression analysis in dimensional X-ray computed tomography | |
CN110211189A (zh) | ToF相机深度误差建模校正方法及装置 | |
CN114417681B (zh) | 基于动态决策和神经网络的二维结构变形监测方法及装置 | |
CN114266776B (zh) | 一种应用复合裂纹位移场函数的数字图像相关方法 | |
CN107587955B (zh) | 基于深信度网络的火箭发动机推力偏移量的标定方法 | |
CN105387826A (zh) | 用于对尺寸偏差和工艺能力进行定量的方法和装置 | |
CN115752243A (zh) | 一种测量数据融合方法 | |
CN111210409B (zh) | 一种基于条件生成对抗网络的结构损伤识别方法 | |
CN109631813B (zh) | 一种大尺寸关节臂式坐标测量机的标定方法 | |
CN107607182B (zh) | 一种卡车称重系统以及称重方法 | |
CN105160069A (zh) | 一种基于改进的紧凑式教学优化算法的机械参数软测量方法 | |
CN111210877A (zh) | 一种推断物性参数的方法及装置 | |
CN113607083B (zh) | 一种基于深度学习的不规则形状光学元件面形的旋转平移绝对检测方法、设备及存储介质 | |
CN115546210B (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 |