CN112800654B - 基于dscm-femu的岩石非均匀力学参数反演方法 - Google Patents

基于dscm-femu的岩石非均匀力学参数反演方法 Download PDF

Info

Publication number
CN112800654B
CN112800654B CN202110131501.3A CN202110131501A CN112800654B CN 112800654 B CN112800654 B CN 112800654B CN 202110131501 A CN202110131501 A CN 202110131501A CN 112800654 B CN112800654 B CN 112800654B
Authority
CN
China
Prior art keywords
unit
rock
finite element
elastic modulus
poisson ratio
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
CN202110131501.3A
Other languages
English (en)
Other versions
CN112800654A (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.)
China University of Mining and Technology Beijing CUMTB
Original Assignee
China University of Mining and Technology Beijing CUMTB
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 China University of Mining and Technology Beijing CUMTB filed Critical China University of Mining and Technology Beijing CUMTB
Priority to CN202110131501.3A priority Critical patent/CN112800654B/zh
Publication of CN112800654A publication Critical patent/CN112800654A/zh
Application granted granted Critical
Publication of CN112800654B publication Critical patent/CN112800654B/zh
Expired - Fee Related 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/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
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • 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

Landscapes

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

Abstract

本发明涉及一种基于DSCM‑FEMU的岩石非均匀力学参数反演方法。基于数字散斑相关方法观测岩石加载实验过程中试件的变形，计算应变场；基于有限元模型修正法，在ABAQUS中建立实验加载的数值计算模型，给定各单元初始的弹性参数，使用Python程序基于weibull分布对各单元的弹性参数赋值，进行有限元计算，通过Python程序自动输出数值模型的应力场；基于有限元理论，推导岩石非均匀力学参数反演目标函数，表示实验测量的应变场与数值计算应变场之间的差异，使用优化算法使目标函数最小化，输出各单元最匹配的弹性参数，若目标函数不满足迭代收敛条件，则通过Python程序自动修改各单元弹性参数，代入新一轮的有限元计算中，直到目标函数收敛，输出各单元的弹性参数。

Description

基于DSCM-FEMU的岩石非均匀力学参数反演方法
技术领域
本发明涉及一种基于DSCM-FEMU的岩石非均匀力学参数反演方法，属于光测力学、变形测量和数值模拟的技术领域。
背景技术
数字散斑相关方法(Digital Speckle Correlation Method,DSCM)是一种非接触式光学测量方法，其基于“灰度不变假设”，通过相关算法追踪随机散斑点在变形过程中的位置变化，具有全场性和高精度等特点。近年来，数字散斑相关方法在岩土工程多个领域的实验研究中被广泛应用，如在岩石材料实验加载过程中，通过数字散斑相关方法可以计算得到试件加载过程中的位移场和应变场，但是此方法不能对材料的应力场和材料参数进行计算测量，这一问题需要进行深入研究。
有限元模型修正法(Finite Element Model Updating,FEMU)是一种将实验手段和数值计算手段相结合进行参数计算的方法，在数值计算中给定模型力学参数代入计算，使用相关的算法不断更新迭代力学参数，使得数值计算结果数据与实验测量数据的差值最小。有限元模型修正法已经用于弹塑性、粘弹性和非线性剪切等条件下的材料参数识别，具有计算速度快、结果精度高等特点。
因此，将数字散斑相关方法和有限元模型修正方法相结合，构建DSCM-FEMU 法，对岩石的非均匀力学参数进行反演，可有效的解决数字散斑相关方法不能计算材料力学参数和实验应力场的问题。同时，基于数字散斑相关方法测量的非均匀应变场和有限元模型修正方法给定的非均匀力学参数，可有效解决由于岩石材料非均匀性导致的力学参数测量不精确的问题。
发明内容
针对现有技术的不足，本发明的目的是提供一种基于DSCM-FEMU的岩石非均匀力学参数反演方法，可计算岩石加载实验中试件弹性模量和泊松比的非均匀量值。将实验中测量的应变场、数值计算中的应力场和力学参数构建参数反演目标函数，使得实验测量得到的全场数据与数值计算的全场数据差异最小，得到实验中与岩石试件最匹配的力学参数。
为了实现上述目的，本发明所采用的技术方案是：
基于DSCM-FEMU的岩石非均匀力学参数反演方法，包括以下步骤：
(1)设计岩石加载实验方案，确定实验过程的边界条件；根据实验方案对岩石试件进行加载，同时，使用数字散斑相关方法观测试件在加载过程中的变形情况，加载完成后计算试件在加载过程中的应变场，在MATLAB软件中读取应变场数据1；
(2)根据步骤(1)的实验加载边界条件，在ABAQUS软件中建立实验的数值计算模型，编写Python程序，根据weibull分布函数给定每个单元的弹性模量E(i)和泊松比μ(i)，其中i为数值模型的单元编号，将各单元的弹性模量和泊松比导入 ABAQUS中，给每个单元赋予不同的力学参数属性，在ABAQUS中进行有限元计算；
(3)通过Python编写程序自动输出步骤(2)数值计算结果中的应力场数据，在MATLAB软件中读取不同单元的应力场数据1：
Figure BDA0002925479890000021
采用三次样条插值算法计算出不同单元与数值模型应力场坐标一致的数字散斑相关方法像素点应变场数据2：
Figure BDA0002925479890000022
(4)根据有限元理论，将步骤(3)中得到的应变场数据2与由不同单元的弹性模量、泊松比和应力场数据1计算得到的数值模拟应变场数据3的差异作为参数反演目标函数，构建岩石非均匀力学参数反演目标函数，将目标函数写入 MATLAB软件；
所述岩石非均匀力学参数反演目标函数为：
Figure BDA0002925479890000023
Figure BDA0002925479890000024
分别为数值计算中单元i的x方向应力、y方向应力、剪应力；
Figure BDA0002925479890000025
分别为在实验中通过数字散斑相关方法和在后处理中通过三次样条插值算法计算得到的单元i的x方向应变、y方向应变、剪应变；n为有限元模型中的单元数量；E(i)为单元i的岩石试件的弹性模量，μ(i)为单元i的泊松比；
(5)使用最小二乘法“trust-region-reflective”算法将目标函数最小化，给定收敛条件，若迭代满足收敛条件，输出此时有限元模型中各个单元的弹性模量和泊松比；若迭代不满足收敛条件，编写Python程序将优化算法给定更新后各单元的弹性模量和泊松比替代上一次迭代的弹性模量和泊松比，代入到有限元中进行新一轮计算；
(6)循环步骤(3)～步骤(5)，直到目标函数迭代收敛，得到与岩石试件最匹配的非均匀弹性模量和泊松比。
本发明有益效果：
本发明方法将数字散斑相关方法和有限元模型修正方法相结合，构建 DSCM-FEMU法，通过数字散斑相关方法观测岩石试件在加载过程中的变形，选取合理的物面分辨率和分析步长等参数，计算应变场。根据有限元模型修正法，给定岩石试件的初始力学参数，建立实验加载的数值计算模型，不同单元的力学参数量值不同，导出数值计算结果中的应力场。将实验中测量的应变场、数值计算中的应力场和各单元的力学参数构建参数反演目标函数，对岩石的非均匀力学参数进行反演，不仅可以通过数字散斑相关方法获得全场的变形数据作为边界条件代入到有限元数值计算中，使得有限元计算的边界条件更加准确，而且可以计算得到岩石的非均匀力学参数，进而计算加载过程中的应力场数据，对研究岩石材料的力学特性起到重要作用。
附图说明
图1为本发明方法流程框图。
图2为本发明中实施例的试件加载示意图。
图3为本发明中实施例的实验荷载-时间曲线。
图4为本发明中实施例的数值计算模型。
具体实施方式
以下结合实施例对本发明的具体实施方式作进一步详细说明。
实施例
本发明基于DSCM-FEMU的岩石非均匀力学参数反演方法的实现过程主要分为三部分(如图1所示)，第一部分，使用数字散斑相关方法观测岩石实验加载过程，计算应变场数据；第二部分，给定各单元初始的弹性参数，使用ABAQUS 有限元软件对实验过程进行数值模拟，导出应力场数据；第三部分，将实验的应变场数据和数值计算的应力场数据以及各单元弹性参数(弹性模量、泊松比)构建目标函数，通过优化算法使得实验测量数据与数值模拟数据差异最小化，输出岩石试件各单元的弹性参数。整个过程通过MATLAB调用ABAQUS和Python程序完成。
实施材料选取岩土实验中常用的花岗岩，计算测试软硬件环境为64位 Windows10操作系统、酷睿i7-7700HQ CPU(2.80GHZ)、8GB内存、1TB硬盘。数字散斑相关方法图像采集系统由CCD相机、冷光源和计算机组成，在实验过程中获取试件表面变形情况的散斑图像，转换为bmp格式进行存储，用于VIC-2D 软件进行后处理计算应变场。下面以花岗岩三点弯曲实验为例进行描述。
实施例应用本发明方法的基本过程分为以下10个步骤：
(1)将花岗岩材料加工成尺寸为400mm×100mm×50mm的长方体试件，打磨光滑，选取400mm×100mm的平面作为观测平面，在该平面上制作黑底白斑随机分布的人工散斑场，试件尺寸及支座位置如图2所示。
(2)实验加载装置为液压伺服控制试验机，加载方式为位移加载，加载速率为0.05mm/min。数字散斑相关方法图像采集系统采集速率为2帧/秒，物面分辨率为0.091mm/pixel，散斑图像的分辨率为1600pixel×1200pixel。
(3)实验开始前，对试件进行预加载，然后调整冷光源的位置和亮度，调整 CCD相机的位置、亮度和焦距，使计算机采集系统中的散斑场像素点清晰呈现，将加载系统与采集系统对时。
(4)实验开始，同时开启试验机的“加载”按钮和图像采集系统的“采集”按钮，当试件完成加载过程，点击“停止”按钮完成实验的加载和采集工作，保存加载数据和图像。实验加载的荷载-时间曲线如图3所示。
(5)使用VIC-2D软件对应变场数据进行分析计算。以图3中标识点0对应时刻的散斑图像作为参考图像，等时间间隔选取标识点1～4，计算对应时刻的应变场。在MATLAB中读取实验测量得到的不同单元的应变场数据1。
(6)在ABAQUS有限元软件中建立岩石三点弯曲实验的数值计算模型，尺寸为400mm×100mm×50mm，网格属性使用减缩积分单元，长度10mm，共2000 个单元，各单元弹性模量和泊松比通过weibull分布函数给定，其中给定均质度 m＝5，弹性模量均值E0＝36GPa，泊松比均值μ0＝0.18，使用自主编写的Python 程序对各个单元的力学参数自动赋值。上下端施加位移约束，以0.05mm/min的速率加载。具有非均匀初始力学参数的模型示意图如图4所示。
(7)计算完成后，使用自主编写的Python程序自动输出数值计算的应力场数据，在MATLAB中读取数值计算得到的不同单元的应力场数据1：
Figure BDA0002925479890000041
采用三次样条插值算法计算出不同单元与数值模型应力场数据坐标一致的数字散斑相关方法像素点应变场数据2：
Figure BDA0002925479890000051
其中i为数值模型的单元编号。
(8)根据有限元理论，将不同单元的应变场数据2与由不同单元的应力场数据1、弹性模量、泊松比计算得到的数值模拟应变场数据3的差异作为参数反演目标函数，构建岩石非均匀力学参数反演目标函数，将目标函数写入MATLAB软件，岩石非均匀力学参数反演目标函数为：
Figure BDA0002925479890000052
Figure BDA0002925479890000053
分别为数值计算中单元i的x方向应力、y方向应力、剪应力；
Figure BDA0002925479890000054
分别为在实验中通过数字散斑相关方法和在后处理中通过三次样条插值算法计算得到的单元i的x方向应变、y方向应变、剪应变；n为有限元模型中的单元数量；E(i)为单元i的岩石试件的弹性模量，μ(i)为单元i的泊松比。
(9)使用最小二乘法“trust-region-reflective”算法将目标函数最小化，收敛条件设置为ΔQ≤10-6，对于弹性模量E(i)和泊松比μ(i)，给定约束范围 10GPa≤E(i)≤100GPa，0.05<μ(i)<0.35。当迭代满足收敛条件时，输出各单元最优的弹性模量和泊松比。
(10)当迭代不满足收敛条件时，使用自行编写的Python程序，根据优化算法给定更新的弹性参数，修改数值模型中各个单元的弹性模量E(i)和泊松比μ(i)，重新导入ABAQUS进行新一轮的有限元计算，直到目标函数满足收敛条件时停止迭代，输出各单元的弹性模量和泊松比。
以图3中的标识点4为例，基于传统力学反演方法，计算区域得到的均匀力学参数为E＝35.2GPa，μ＝0.15，使用基于DSCM-FEMU的岩石非均匀力学参数反演方法计算结果见表1。
由表1可知，各个单元的弹性模量和泊松比各不相同，弹性模量为 12.69～51.30GPa，泊松比为0.05～0.32。岩石是一种具有非均匀性的天然材料，是由许多大小、形状和性质各不相同的矿物颗粒和胶结物组成的混合物，在细观上则表现为岩石的各细观结构物理力学性质的不均匀性，即岩石内部具有不同的弹性模量、泊松比、破坏强度等。细观上的非均匀性是造成岩石在外荷载作用下表现出宏观不均匀性和非线性特征的主要原因。
首先，通过本文给出的基于DSCM-FEMU的岩石非均匀力学参数反演方法，计算得到的非均匀力学参数相对于均匀的力学参数，更接近于岩石的真实非均匀特性；其次，基于DSCM-FEMU的岩石非均匀力学参数反演方法有效解决了由于岩石材料非均匀性导致的力学参数测量不精确的问题，岩石的非均匀力学参数可以作为岩土工程、采矿工程、水利水电工程等相关领域中理论推导和数值计算的数据基础，其精确性具有重要影响；最后，根据应力应变关系，得到更加精确的应力场数据，有效解决了数字散斑相关方法不能计算实验应力场的问题。
表1
Figure BDA0002925479890000071

Claims (1)

1.基于DSCM-FEMU的岩石非均匀力学参数反演方法，其特征在于，包括以下步骤：
(1)设计岩石加载实验方案，确定实验过程的边界条件；根据实验方案对岩石试件进行加载，同时，使用数字散斑相关方法观测试件在加载过程中的变形情况，加载完成后计算试件在加载过程中的应变场，在MATLAB软件中读取应变场数据1；
(2)根据步骤(1)的实验加载边界条件，在ABAQUS软件中建立实验的数值计算模型，编写Python程序，根据weibull分布函数给定每个单元的弹性模量E(i)和泊松比μ(i)，其中i为数值模型的单元编号，将各单元的弹性模量和泊松比导入ABAQUS中，给每个单元赋予不同的力学参数属性，在ABAQUS中进行有限元计算；
(3)通过Python编写程序自动输出步骤(2)数值计算结果中的应力场数据，在MATLAB软件中读取不同单元的应力场数据1：
Figure FDA0002925479880000011
采用三次样条插值算法计算出不同单元与数值模型应力场坐标一致的数字散斑相关方法像素点应变场数据2：
Figure FDA0002925479880000012
(4)根据有限元理论，将步骤(3)中得到的应变场数据2与由不同单元的弹性模量、泊松比和应力场数据1计算得到的数值模拟应变场数据3的差异作为参数反演目标函数，构建岩石非均匀力学参数反演目标函数，将目标函数写入MATLAB软件；
所述岩石非均匀力学参数反演目标函数为：
Figure FDA0002925479880000013
Figure FDA0002925479880000014
分别为数值计算中单元i的x方向应力、y方向应力、剪应力；
Figure FDA0002925479880000015
分别为在实验中通过数字散斑相关方法和在后处理中通过三次样条插值算法计算得到的单元i的x方向应变、y方向应变、剪应变；n为有限元模型中的单元数量；E(i)为单元i的岩石试件的弹性模量，μ(i)为单元i的泊松比；
(5)使用最小二乘法“trust-region-reflective”算法将目标函数最小化，给定收敛条件，若迭代满足收敛条件，输出此时有限元模型中各个单元的弹性模量和泊松比；若迭代不满足收敛条件，编写Python程序将优化算法给定更新后各单元的弹性模量和泊松比替代上一次迭代的弹性模量和泊松比，代入到有限元中进行新一轮计算；
(6)循环步骤(3)～步骤(5)，直到目标函数迭代收敛，得到与岩石试件最匹配的非均匀弹性模量和泊松比。
CN202110131501.3A 2021-01-30 2021-01-30 基于dscm-femu的岩石非均匀力学参数反演方法 Expired - Fee Related CN112800654B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110131501.3A CN112800654B (zh) 2021-01-30 2021-01-30 基于dscm-femu的岩石非均匀力学参数反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110131501.3A CN112800654B (zh) 2021-01-30 2021-01-30 基于dscm-femu的岩石非均匀力学参数反演方法

Publications (2)

Publication Number Publication Date
CN112800654A CN112800654A (zh) 2021-05-14
CN112800654B true CN112800654B (zh) 2021-07-30

Family

ID=75813141

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110131501.3A Expired - Fee Related CN112800654B (zh) 2021-01-30 2021-01-30 基于dscm-femu的岩石非均匀力学参数反演方法

Country Status (1)

Country Link
CN (1) CN112800654B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113340732B (zh) * 2021-05-31 2022-05-03 北京理工大学 基于图像自动分区的非均质材料多参数反演方法及装置
CN113962118A (zh) * 2021-06-28 2022-01-21 北方工业大学 一种岩土材料力学参数识别方法
CN113804248B (zh) * 2021-08-24 2023-09-22 中国石油大学(华东) 利用数字散斑和有限元技术的无损地应力测试装置及方法
CN116642750B (zh) * 2023-07-24 2023-10-20 长江三峡集团实业发展(北京)有限公司 一种岩石应变局部化起始时间的预测方法、装置及设备
CN117968992B (zh) * 2024-04-01 2024-06-04 泰富特钢悬架(济南)有限公司 一种钢板弹簧用钢板弹性检测装置

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108182335A (zh) * 2018-01-26 2018-06-19 山东科技大学 一种基于abaqus的岩石力学试验数值仿真方法
CN109870376A (zh) * 2019-02-03 2019-06-11 浙江大学 一种基于纳米压痕和数值模拟反演岩石矿物参数的方法
CN111414658A (zh) * 2020-03-17 2020-07-14 宜春学院 一种岩体力学参数反分析方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10013767B2 (en) * 2013-11-01 2018-07-03 The Research Foundation For The State University Of New York Method for measuring the interior three-dimensional movement, stress and strain of an object

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108182335A (zh) * 2018-01-26 2018-06-19 山东科技大学 一种基于abaqus的岩石力学试验数值仿真方法
CN109870376A (zh) * 2019-02-03 2019-06-11 浙江大学 一种基于纳米压痕和数值模拟反演岩石矿物参数的方法
CN111414658A (zh) * 2020-03-17 2020-07-14 宜春学院 一种岩体力学参数反分析方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Applications of digital correlation method to structure inspection;Junda Chen 等;《TSINGHUA SCIENCE AND TECHNOLOGY》;20070630;第12卷(第3期);273-243 *
基于数字散斑相关方法的岩石力学参数反演研究;吴佳宁 等;《北方工业大学学报》;20200430;第32卷(第2期);135-140 *
等幅循环加载岩石非均匀变形演化试验研究;杨小彬 等;《采矿与安全工程学报》;20190331;第36卷(第2期);388-395 *

Also Published As

Publication number Publication date
CN112800654A (zh) 2021-05-14

Similar Documents

Publication Publication Date Title
CN112800654B (zh) 基于dscm-femu的岩石非均匀力学参数反演方法
Zhang et al. Micromechanical modelling of deformation and fracture of hydrating cement paste using X-ray computed tomography characterisation
Zhao et al. The state of the art of two‐dimensional digital image correlation computational method
Gadelrab et al. Densification modeling of fused silica under nanoindentation
CN110532591B (zh) 基于dic-efg联合仿真分析裂纹尖端应变场的方法
Kalin et al. Comparing surface topography parameters of rough surfaces obtained with spectral moments and deterministic methods
CN111339700B (zh) 核电汽轮机叶片的疲劳损伤评估方法、装置和存储介质
CN104809362B (zh) 一种基于数值计算的包含非均匀变形散斑图制作方法
Huang et al. A digital volume correlation technique for 3-D deformation measurements of soft gels
CN103076347B (zh) 基于原位x射线断层照相的脆性材料力学损伤的测量方法
JP7356823B2 (ja) パラメータ推定装置、パラメータ推定方法及びプログラム
CN107590787B (zh) 一种扫描电子显微镜的图像畸变校正方法
Romani et al. Detection of crack onset in double cleavage drilled specimens of plaster under compression by digital image correlation–theoretical predictions based on a coupled criterion
Tal et al. Direct observations of damage during unconfined brittle failure of Carrara marble
Luo et al. Characterization of the stiffness distribution in two and three dimensions using boundary deformations: a preliminary study
Shi et al. Micromorphological characterization and random reconstruction of 3D particles based on spherical harmonic analysis
Garg et al. An integrated approach for evaluation of linear cohesive zone model’s performance in fracturing of rocks
Lucadamo et al. Characterization and simulation methods applied to the study of fretting wear in Zircaloy-4
CN112730102B (zh) 一种基于标准差椭圆的颗粒类材料剪切带演化测定方法
Mao et al. 3-D strain estimation in sandstone using improved digital volumetric speckle photography algorithm
Nakhodchi et al. The formation of fracture process zones in polygranular graphite as a precursor to fracture
Krishna et al. Numerical and experimental study on flexural behavior of reinforced concrete beams: Digital image correlation approach
CN114323950A (zh) 基于力学试验结果的分子动力学模型验证方法
Janíček et al. Accuracy testing of the DIC optical measuring method the Aramis system by GOM
Lee et al. A study of plastic zone formation by digital image processing

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
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20210730

Termination date: 20220130

CF01 Termination of patent right due to non-payment of annual fee