CN106295045B - 识别拉杆转子界面局部脱开位置和程度的广义柔度矩阵法 - Google Patents

识别拉杆转子界面局部脱开位置和程度的广义柔度矩阵法 Download PDF

Info

Publication number
CN106295045B
CN106295045B CN201610687339.2A CN201610687339A CN106295045B CN 106295045 B CN106295045 B CN 106295045B CN 201610687339 A CN201610687339 A CN 201610687339A CN 106295045 B CN106295045 B CN 106295045B
Authority
CN
China
Prior art keywords
rod fastening
fastening rotor
interface
rotor
partial disengagement
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.)
Expired - Fee Related
Application number
CN201610687339.2A
Other languages
English (en)
Other versions
CN106295045A (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.)
Xian Jiaotong University
Original Assignee
Xian Jiaotong University
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 Xian Jiaotong University filed Critical Xian Jiaotong University
Priority to CN201610687339.2A priority Critical patent/CN106295045B/zh
Publication of CN106295045A publication Critical patent/CN106295045A/zh
Application granted granted Critical
Publication of CN106295045B publication Critical patent/CN106295045B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/23Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/10Geometric CAD
    • G06F30/17Mechanical parametric or variational design
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00Details relating to CAD techniques
    • G06F2111/04Constraint-based CAD

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)
  • Turbine Rotor Nozzle Sealing (AREA)

Abstract

本发明公开了一种识别拉杆转子界面局部脱开位置和程度的广义柔度矩阵法,通过对系统施加振型约束,得到能够使用任意单阶模态识别拉杆转子界面局部脱开位置和程度的改进广义柔度矩阵法。当燃气轮机拉杆转子界面发生局部脱开时,通过实测其临界转速和轴心涡动轨迹,结合转子静止时的固有频率和振型,即可监测出拉杆转子界面脱开的位置和程度。本发明能够使用任意单阶模态参数识别拉杆转子界面脱开的位置和程度,可为重型燃气轮机拉杆转子安全稳定运行及故障诊断提供指导。

Description

识别拉杆转子界面局部脱开位置和程度的广义柔度矩阵法
技术领域
本发明涉及重型燃气轮机、航空发动机、舰船用燃气轮机中的拉杆转子界面局部脱开识别方法,具体涉及一种使用改进的广义柔度矩阵法识别拉杆转子预紧松弛时界面局部脱开的方法。
背景技术
拉杆转子具有便于冷却通道设计、易于加工等优点,广泛应用于航空发动机和重型燃气轮机。与一般蒸汽轮机转子不同,拉杆转子不再是一个连续的整体,而是通过拉杆将多级轮盘压紧组合在一起。拉杆在长期交变载荷的作用下发生松弛、航空发动机空中停车都会导致拉杆预紧力的减小,此时拉杆转子的重力及振动可能引起轮盘间出现脱开、滑移等情况。拉杆转子界面的局部脱开可能导致振动幅值的进一步增大,甚至发生转子静子之间的碰摩,严重时导致停机甚至整个结构破坏。随着燃气轮机向大功率、高效率、高燃气温度方向发展,对拉杆转子的安全运行及故障诊断提出了更高要求。
拉杆转子各界面上弯矩随着转速改变而变化,通过不同阶次临界转速时各界面的接触状态可能不同,只能使用单阶模态参数对界面分离发生的位置和程度进行识别。而目前文献中的结构损伤识别方法由于需要多阶模态参数、未考虑陀螺效应,无法适用于识别拉杆转子界面的局部脱开。同时,文献中的结构损伤识别方法在实际应用中未考虑阻尼和测量误差的影响。
发明内容
针对目前损伤识别方法应用于燃气轮机拉杆转子界面局部分离所存在的问题,本发明提出了一种识别拉杆转子界面局部脱开位置和程度的广义柔度矩阵法。本发明为改进的广义柔度矩阵法,能够使用任意单阶模态识别燃气轮机拉杆转子界面局部脱开的位置和程度。
为达到以上目的,本发明采用了以下技术方案:
1)基于拉杆转子有限元模型建立完整无阻尼系统,对所述系统施加振型约束,得到约束系统;
2)经过步骤1)后,使用拉杆转子临界转速和涡动轨迹,结合约束系统静止时的固有频率和振型,计算界面局部脱开导致的约束系统的广义柔度矩阵变化量,利用所述广义柔度矩阵变化量以及拉杆转子各个相邻界面之间刚度损伤与所述广义柔度矩阵变化量的关系,对拉杆转子相邻界面之间的刚度减少程度进行求解,根据求解结果对拉杆转子界面脱开的位置和程度进行识别。
利用发生界面局部脱开前后转子的任意阶次临界转速变化量和所述固有频率,得到拉杆转子的耦合固有频率变化量;利用发生界面局部脱开前后转速为对应阶次临界转速时的转子轴心涡动轨迹变化量和所述振型,得到拉杆转子的耦合位移变化量;通过所述耦合固有频率变化量与耦合位移变化量并结合拉杆转子位移转换关系,计算所述广义柔度矩阵变化量。
所述广义柔度矩阵变化量具有以下两种表示形式:
以及
其中,n为拉杆转子界面数,r为小于n的自然数,di为单元刚度损伤系数,表示相邻界面之间刚度的减少程度,ωr+1为界面未发生局部脱开拉杆转子静止时的第r+1阶固有频率计算值;Φr+1为界面未发生局部脱开拉杆转子静止时的第r+1阶振型计算值;ωc,r+1为界面未发生局部脱开拉杆转子的第r+1阶临界转速计算值;为界面发生局部脱开拉杆转子的第r+1阶临界转速实测值;Zr+1为拉杆转子界面未发生局部脱开时通过r+1阶临界转速时的轴心涡动轨迹计算值;Zd,r+1为拉杆转子界面发生局部脱开后通过r+1阶临界转速时的轴心涡动轨迹实测值,D为位移向量约束转换矩阵。
所述拉杆转子位移转换关系的获取方法包括以下步骤:
a)假设完整无阻尼系统共有n=r+s阶模态,n为拉杆转子界面数,r为小于n的自然数,利用前r阶振型对系统施加约束,剔除系统的前r阶模态:
b)将系统振型矩阵Φ、质量矩阵M和位移向量x的分块形式代入步骤a)中公式后进行展开,得到位移向量xs与x之间的转换关系为:
x=Dxs
其中,I为r阶单位阵。
所述求解采用最小二乘法。
本发明的有益效果体现在:
与现有损伤识别方法相比,本发明方法的优点是:1)可以使用任意单阶模态参数识别拉杆转子界面局部分离的位置和程度;2)考虑了旋转导致的陀螺效应;3)适用于多处损伤的识别。此外,本发明消除了阻尼和转动对结构固有频率和振型的影响。并克服了随机误差对临界转速和轴心涡动轨迹进行界面局部脱开识别的影响。通过本发明方法,可以在燃气轮机拉杆转子界面发生局部分离的初期得到分离的位置及程度,减少停机时间,防止灾难性的事故。
附图说明
图1是拉杆转子节点划分示意图;图中:1为某界面,2为节点,3为梁单元。
图2是拉杆转子各轴段等效抗弯刚度示意图。
图3是基于改进广义柔度矩阵法的拉杆转子界面局部脱开识别方法流程图。
图4是本发明算例中界面局部脱开识别结果。
具体实施方式
下面结合附图和实施例对本发明作详细说明。
第一、建立拉杆转子有限元模型(完整无阻尼系统)
如图1所示,对某重型燃气轮机拉杆转子有限元模型进行节点划分。按照界面划分原则(轮盘接触面中心、轴径中心、轴截面未发生突变处以及轴端点处)设置n个节点(n=38),并按顺序编号。依据设置的节点将拉杆转子有限元模型共划分为n-1个梁单元(轴段)。
参见图2,对所述重型燃气轮机拉杆转子分段进行三维有限元计算,根据各节点在载荷作用下的位移,得出轴段两端面的相对转角θ,利用弯矩与转角关系分别计算拉杆转子各轴段等效抗弯刚度:
式中:Kbsn为轴段等效抗弯刚度;M为对轴段施加的弯矩;l为轴段(例如图2所示轮盘)沿轴向的长度,通过计算所得各轴段等效抗弯刚度得到梁单元刚度矩阵,再叠加成拉杆转子的整体刚度矩阵K。
第二、通过建立的拉杆转子有限元模型,结合施加振型约束前后的振型的转换关系,得出施加振型约束后系统的振动微分方程。
系统共有n=r+s阶模态,利用前r阶振型对系统施加约束,剔除系统的前r阶模态:
式中:Φi为拉杆转子的第i个振型;M为拉杆转子的质量矩阵;x为位移向量。系统振型矩阵Φ、质量矩阵M和位移向量x可以写成分块形式:
把式(2)代入式(1)可得到:
将式(3)展开,可得到xr=-(ΦrrMrrrsMsr)-1rrMrsrsMss)xs,即可得到位移向量xs与x之间的转换关系为:
x=Dxs (4)
式中:I为r阶单位矩阵。
将式(4)代入系统原始振动控制方程中,可得到施加振型约束后的振动控制方程:
式中:施加振型约束后拉杆转子的质量矩阵Mu=DTMD;施加振型约束后拉杆转子的整体刚度矩阵Ku=DTKD。
第三、计算施加振型约束后的拉杆转子有限元模型的固有频率和振型;同时考虑旋转导致的陀螺效应,计算未发生界面局部脱开故障(即无损)的拉杆转子各阶临界转速,及转速为各阶临界转速时的轴心涡动轨迹。
其中,施加振型约束后的系统和原始系统的模态参数之间的关系为:
式中:分别为施加振型约束后系统的第i个特征值(固有频率)和特征向量(振型)。λi+r和Φi+r分别为原始系统的第i+r个特征值和特征向量。由于D不是方阵,D-1为D的一个广义逆矩阵。至此在原始n自由度系统的基础上建立了一个自由度为n-r的约束系统。约束系统的首阶模态参数可以由原始系统的第r+1阶模态参数得到。
第四、计算约束系统的广义柔度矩阵变化量。
拉杆转子界面局部分离(即局部脱开)只会导致其刚度的减小,其质量属性不发生改变。定义界面局部分离参数为di,代表第i个轴段刚度的减少程度,即单元刚度损伤系数。结合施加振型约束前后质量阵和刚度阵之间的转换关系(Mu=DTMD和Ku=DTKD),可得到施加振型约束后柔度矩阵Fu对损伤(局部脱开)的敏感度,i=1,…,n-1:
Kai表示第i个梁单元的刚度矩阵。
施加振型约束后系统的广义柔度矩阵为:
fu=FuMuFuMuFu (8)
将式(8)对单元刚度损伤系数di求导,可得约束系统的广义柔度矩阵对损伤的敏感度:
结合式(7)和式(9),可以得到约束系统的广义柔度矩阵对单元刚度损伤的敏感度。利用泰勒展开并舍去高阶项,损伤导致的约束系统的广义柔度矩阵变化量为:
在式(10)的基础上,可以将约束系统的广义柔度矩阵变化量表示为:
式中:ωr+1为界面未发生分离拉杆转子静止时的第r+1阶固有频率计算值;Φr+1为界面未发生分离拉杆转子静止时的第r+1阶振型计算值;ωc,r+1为界面未发生分离拉杆转子的第r+1阶临界转速计算值;为界面发生分离拉杆转子的第r+1阶临界转速实测值;Zr+1为拉杆转子界面未发生分离时通过r+1阶临界转速时的轴心涡动轨迹计算值;Zd,r+1为转子界面发生分离后通过r+1阶临界转速时的轴心涡动轨迹实测值。
第五、监测实际拉杆转子的临界转速,得到拉杆转子发生界面局部脱开后的临界转速。在转子多处布置位移传感器,得到各界面轴心的涡动轨迹,由振型扩展方法得到相应阶次的完整振型。利用测得的拉杆转子临界转速和轴心涡动轨迹Zd,r+1,结合拉杆转子有限元模型的固有频率计算值ωr+1,临界转速计算值ωc,r+1和轴心涡动轨迹计算值,使用最小二乘法求解式(10)和式(11)组成的线性方程,即可得到各轴段对应的单元刚度损伤系数di。由于式中参数r为自然数(r<n),表明此方法可以仅利用任意单阶模态参数对结构损伤进行识别。
参见图3,本发明通过比较发生界面局部脱开前后拉杆转子的临界转速,得到临界转速变化量,利用发生界面局部脱开前后拉杆转子的临界转速变化量和有限元模型计算所得固有频率,得到转子的耦合固有频率变化量以任意单阶振型进行计算,比较发生界面局部脱开前后转速为相应阶次临界转速时的转子轴心涡动轨迹,得到转速为相应阶次临界转速时的转子轴心涡动轨迹变化量,利用转子轴心涡动轨迹变化量和未发生脱开时的振型(有限元模型计算所得),得到转子的耦合位移变化量(Φr+1+Zd,r+1-Zr+1)。通过耦合固有频率变化量与耦合位移变化量并结合转子位移转换关系(x=Dxs),计算出广义柔度矩阵变化量。最终基于转子的广义柔度矩阵变化量的两种表示形式(式(10)和式(11)),并使用最小二乘法求解得到拉杆转子各处界面局部脱开的程度。
通过数值算例验证本发明所提出方法。假设拉杆转子上三个界面出现局部分离,受影响的四个轴段编号分别为10、11、12和13,抗弯刚度减小量分别为原始刚度的5%、16%、12%和2%。考虑转动导致的陀螺效应等因素影响,计算得到界面发生分离后各阶临界转速和轴心涡动轨迹。图4为仅使用第二阶临界转速和轴心涡动轨迹得到的在拉杆预紧松弛时的界面局部脱开识别结果,轴段10到13识别出的损伤程度分别为0.053、0.224、0.144和0.008。数值算例结果表明,本发明提出方法能够利用拉杆转子单阶临界转速和相应轴心涡动轨迹识别出界面局部分离的程度和位置。

Claims (5)

1.一种识别拉杆转子界面局部脱开位置和程度的广义柔度矩阵法,其特征在于:包括下述步骤:
1)基于拉杆转子有限元模型建立完整无阻尼系统,对所述系统施加振型约束,得到约束系统;
2)经过步骤1)后,使用拉杆转子临界转速和涡动轨迹,结合约束系统静止时的固有频率和振型,计算界面局部脱开导致的约束系统的广义柔度矩阵变化量,利用所述广义柔度矩阵变化量以及拉杆转子各个相邻界面之间刚度损伤与所述广义柔度矩阵变化量的关系,对拉杆转子相邻界面之间的刚度减少程度进行求解,根据求解结果对拉杆转子界面脱开的位置和程度进行识别;
所述广义柔度矩阵变化量具有以下两种表示形式:
以及
其中,n为拉杆转子界面数,r为小于n的自然数,di为单元刚度损伤系数,表示相邻界面之间刚度的减少程度,ωr+1为界面未发生局部脱开拉杆转子静止时的第r+1阶固有频率计算值;Φr+1为界面未发生局部脱开拉杆转子静止时的第r+1阶振型计算值;ωc,r+1为界面未发生局部脱开拉杆转子的第r+1阶临界转速计算值;为界面发生局部脱开拉杆转子的第r+1阶临界转速实测值;Zr+1为拉杆转子界面未发生局部脱开时通过r+1阶临界转速时的轴心涡动轨迹计算值;Zd,r+1为拉杆转子界面发生局部脱开后通过r+1阶临界转速时的轴心涡动轨迹实测值,D为位移向量约束转换矩阵。
2.根据权利要求1所述一种识别拉杆转子界面局部脱开位置和程度的广义柔度矩阵法,其特征在于:利用发生界面局部脱开前后转子的任意阶次临界转速变化量和所述固有频率,得到拉杆转子的耦合固有频率变化量;利用发生界面局部脱开前后转速为对应阶次临界转速时的转子轴心涡动轨迹变化量和所述振型,得到拉杆转子的耦合位移变化量;通过所述耦合固有频率变化量与耦合位移变化量并结合拉杆转子位移转换关系,计算所述广义柔度矩阵变化量。
3.根据权利要求2所述一种识别拉杆转子界面局部脱开位置和程度的广义柔度矩阵法,其特征在于:所述拉杆转子位移转换关系的获取方法包括以下步骤:
a)假设完整无阻尼系统共有n=r+s阶模态,n为拉杆转子界面数,r为小于n的自然数,利用前r阶振型对系统施加约束,剔除系统的前r阶模态:
Φi为拉杆转子的第i个振型;
b)将系统振型矩阵Φ、质量矩阵M和位移向量x的分块形式代入步骤a)中公式后进行展开,得到位移向量xs与x之间的转换关系为:
x=Dxs
系统振型矩阵Φ、质量矩阵M和位移向量x的分块形式为:
4.根据权利要求3所述一种识别拉杆转子界面局部脱开位置和程度的广义柔度矩阵法,其特征在于:
其中,I为r阶单位阵。
5.根据权利要求1所述一种识别拉杆转子界面局部脱开位置和程度的广义柔度矩阵法,其特征在于:所述求解采用最小二乘法。
CN201610687339.2A 2016-08-18 2016-08-18 识别拉杆转子界面局部脱开位置和程度的广义柔度矩阵法 Expired - Fee Related CN106295045B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610687339.2A CN106295045B (zh) 2016-08-18 2016-08-18 识别拉杆转子界面局部脱开位置和程度的广义柔度矩阵法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610687339.2A CN106295045B (zh) 2016-08-18 2016-08-18 识别拉杆转子界面局部脱开位置和程度的广义柔度矩阵法

Publications (2)

Publication Number Publication Date
CN106295045A CN106295045A (zh) 2017-01-04
CN106295045B true CN106295045B (zh) 2019-07-23

Family

ID=57661401

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610687339.2A Expired - Fee Related CN106295045B (zh) 2016-08-18 2016-08-18 识别拉杆转子界面局部脱开位置和程度的广义柔度矩阵法

Country Status (1)

Country Link
CN (1) CN106295045B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112942450B (zh) * 2021-02-01 2022-06-03 湘电风能有限公司 基础环式风机基础疲劳损伤监测方法及系统
CN113358203B (zh) * 2021-06-02 2024-06-11 中国大唐集团科学技术研究院有限公司华东电力试验研究院 一种基于谐分量分解的转子固有频率识别方法和系统
CN113569349B (zh) * 2021-07-08 2023-10-24 中国科学院工程热物理研究所 一种考虑温度作用的拉杆转子动力学特性的分析方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102222123A (zh) * 2011-03-16 2011-10-19 西安交通大学 用于燃气轮机拉杆式转子扭振模态的计算验证方法
CN103970944A (zh) * 2014-04-29 2014-08-06 西安交通大学 轮盘间平面接触的燃气轮机拉杆转子预紧力设计校核方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102222123A (zh) * 2011-03-16 2011-10-19 西安交通大学 用于燃气轮机拉杆式转子扭振模态的计算验证方法
CN103970944A (zh) * 2014-04-29 2014-08-06 西安交通大学 轮盘间平面接触的燃气轮机拉杆转子预紧力设计校核方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Structural Damage Detection by Using Single Natural Frequency and the Corresponding Mode Shape;Bo Zhao等;《Shock and Vibration》;20160106;第20卷(第1期);第1-8页摘要、第2-3节,图7-图9
Structural Damage Detection Using Generalized Flexibility Matrix and Changes in Natural Frequencies;Jing Li等;《AIAA Journal》;20111202;第50卷(第5期);第1072-1078页
多盘转子系统动力学研究;胡四辈等;《甘肃科学学报》;20160228;第28卷(第1期);第110-114页

Also Published As

Publication number Publication date
CN106295045A (zh) 2017-01-04

Similar Documents

Publication Publication Date Title
CN106295045B (zh) 识别拉杆转子界面局部脱开位置和程度的广义柔度矩阵法
Hah et al. Tip clearance vortex oscillation, vortex shedding and rotating instabilities in an axial transonic compressor rotor
Provenza et al. Aeromechanical response of a distortion-tolerant boundary layer ingesting fan
Wang et al. An improved non-contact dynamic stress measurement method for turbomachinery rotating blades based on fundamental mistuning model
Dhini et al. Neural network based system for detecting and diagnosing faults in steam turbine of thermal power plant
Manwaring et al. Structures and Dynamics Committee Best Paper of 1996 Award: Inlet Distortion Generated Forced Response of a Low-Aspect-Ratio Transonic Fan
Storace et al. Unsteady flow and whirl-inducing forces in axial-flow compressors: part I—experiment
Kumar et al. A review: multiplicative faults and model-based condition monitoring strategies for fault diagnosis in rotary machines
Salhi et al. Identification of modal parameters and aeroelastic coefficients in bladed disk assemblies
Wang et al. Simulation analysis of casing vibration response and its verification under blade–casing rubbing fault
Wei et al. First-order approximate rigid-flexible coupled dynamics analysis of a simple aero-engine blade model with dynamic stiffening effect
Phibel et al. Comparison between a CFD code and a three-control-volume model for labyrinth seal flutter predictions
Eichner et al. Prediction of aerodynamically induced fan blade vibration due to boundary layer ingestion
AbdulBary et al. Fault Diagnosis in Rotating System Based on Vibration Analysis
Manwaring et al. Inlet distortion generated forced response of a low aspect ratio transonic fan
CN109883379A (zh) 基于模态振型的叶片位移应变测量方法
Miyakozawa et al. The effects of aerodynamic asymmetric perturbations on forced response of bladed disks
Rose et al. Experimentally observed unsteady work at inlet to and exit from an axial flow turbine rotor
Yang et al. Influence of Fir‐Tree Tenon Blade Parameters on Blade Dynamic Characteristics
Jacobs et al. Turbomachinery transient dynamics of radial rotor-stator contact occurrences with friction
Song et al. Numerical and Experimental Investigations on the Leakage Flow Characteristics and Temperature Distributions of Brush Seals
Guerin et al. Thermomechanical component mode synthesis for blade casing interaction prediction
Ma et al. Interval analysis method for rotordynamic with uncertain parameters
Tamai et al. Vibration analysis of shrouded turbine blades for a 30 MW gas turbine
Lian et al. Insight on aerodynamic damping of the civil transonic fan blade

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20190723