CN110378060B - 一种顶张式立管随机耦合振动的计算方法 - Google Patents

一种顶张式立管随机耦合振动的计算方法 Download PDF

Info

Publication number
CN110378060B
CN110378060B CN201910682572.5A CN201910682572A CN110378060B CN 110378060 B CN110378060 B CN 110378060B CN 201910682572 A CN201910682572 A CN 201910682572A CN 110378060 B CN110378060 B CN 110378060B
Authority
CN
China
Prior art keywords
formula
inner sleeve
oil pipe
centralizer
sleeve
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
CN201910682572.5A
Other languages
English (en)
Other versions
CN110378060A (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.)
Ocean University of China
Original Assignee
Ocean University of China
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 Ocean University of China filed Critical Ocean University of China
Priority to CN201910682572.5A priority Critical patent/CN110378060B/zh
Publication of CN110378060A publication Critical patent/CN110378060A/zh
Application granted granted Critical
Publication of CN110378060B publication Critical patent/CN110378060B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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
    • 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
    • G06F2119/00Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/06Power 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)
  • Earth Drilling (AREA)

Abstract

本发明涉及海洋深水立管的研究方法,尤其涉及一种顶张式立管随机耦合振动的计算方法。该方法采用三层立管组成的耦合模型进行计算分析,充分考虑了多层立管之间的耦合作用与立管管间间隙的影响。比现有的技术更真实的反映顶张式立管横向振动的运动和变形状态,且能够较精确地直接计算出外套管、内套管和油管的应力状态,而不必向之前一样采用截面荷载分配的方法二次计算外套管、内套管和油管的应力状态。

Description

一种顶张式立管随机耦合振动的计算方法
技术领域
本发明涉及海洋深水立管的研究方法,尤其涉及一种顶张式立管随机耦合振动的计算方法。
背景技术
深水顶张式立管是深水海洋立管的一种主要类型,它在深水油气开发中的作用是采油。深水顶张式立管主要结构形式为油管和同轴心的双套管组成的管中管结构。最内层的是油管,中间的和最外层的为套管。油管用于采油,内套管为气举线,外套管为隔水管,是承受海洋环境荷载的主要结构。为了防止相邻立管之间管体发生碰撞,油管与内套管之间以及内套管与外套管之间设置了扶正器。
由于采油的同时,伴有天然气产生,因此,内套管是作为气举线使用的,即它与油管之间形成的环形空间输送井下的天然气。这意味着,内套管在承受其它外部荷载的同时还承受天然气压力的作用,油管则同时受内压(原油压力)和外压(天然气压力)的作用。外套管主要起隔水和承载作用,其直径和壁厚均大于其它两根钢管,因此,承担了大部分的海洋环境荷载和张力荷载。
目前,在工程上顶张式立管的设计分析共有两种方法。
一是,采用等效管模型,采用弯曲刚度等效的原则将三根管等效为一根管,即等效管的弯曲刚度等于三根管的弯曲刚度之和:
EIeq=EItube+EIinner+EIouter (1)
式中:
EIeq——等效单层管的截面抗弯刚度;
EItube——油管的截面抗弯刚度;
EIinner——内层立管的截面抗弯刚度;
EIouter——外层立管的截面抗弯刚度。
为了保证等效管的模型与实际立管的海洋环境荷载相同,等效管的外径采用外层立管的外径。
二是,不考虑管间间隙对耦合响应的影响。且该分析方法是基于一定假设条件的,即:假设当相邻两根立管在扶正器处接触后,扶正器处两根立管的节点位移、速度、加速度是同步的。
现有的深水顶张式立管随机耦合振动的计算方法存在以下技术缺陷:
1、拉压刚度不等效
当等效管与三根管的弯曲刚度等效时,等效管的内径应等于:
Figure GDA0002815923920000021
而等效管与三根管的拉压刚度等效时,等效管的内径应等于:
Figure GDA0002815923920000022
公式(2)和公式(3)表明,按照弯曲刚度等效和拉压刚度等效计算得到的等效管内径是不同的,即按照弯曲刚度等效和拉压刚度等效计算得到的等效管的壁厚是不同的。因此,用等效管计算得到的轴向应力和变形是不同的,从而等效应力也不相同。所以,无论采用工作应力法或荷载抗力系数法校核结构强度和疲劳都将是不准确的。
2、压力不等效
采用等效管计算时,需要将油管的原油压力和内套管的天然气压力转换为等效管的内压。由于外套管不受油气压力的作用,因此,转换后,等效管的环向应力与三根管的环向应力均不相同,从而导致结构的应力状态与实际不符。所以,无论采用工作应力法或荷载抗力系数法校核结构强度和疲劳都将是不准确的。
3、动力响应不等效
除油气压力外,等效管结构受到的其它荷载与实际结构相同,但是,由于实际结构仅在离散的指定位置上有扶正器相互支撑,因此,等效管的静态力学性能与实际结构相似,但动态力学性能是不相似的。因为,三根管的质量和刚度均不相同,从而动力特性不同,即它们的固有频率和振型是不同的。
4、弯曲应力不等效
按照等效管理论并不能准确求出各层立管上的弯曲应力,按照截面弯曲刚度比值分配的外套管、内套管弯曲应力大于其实际应力,油管实际应力小于其实际应力,造成这种差别的原因是等效管模型不能正确反映各层立管的受力载荷形式以及扶正器在立管耦合运动中的作用。
5、没有考虑管间间隙对耦合响应的影响
已有的分析方法没有考虑管间间隙对耦合响应的影响,进而不能对扶正器高度(轴向位置)进行参数敏感性分析,不能准确模拟多层立管的实际模型。
发明内容
本发明提出一种解决深水顶张式立管随机耦合振动的分析方法,该分析方法计算时无需对立管进行等效处理,也可充分考虑管间间隙对耦合响应的影响。不仅能够解决上述五个问题,也可为扶正器的尺寸选型以及布置方案提供依据。
本发明的技术方案为:
本发明基于三层立管组成的耦合模型进行计算分析且充分考虑了管间间隙对耦合响应的影响,其具体计算方案如下:
(1)分别建立外套管、内套管与油管的运动方程:
Figure GDA0002815923920000041
Figure GDA0002815923920000042
Figure GDA0002815923920000043
式中:
m1、m2、m3分别为外套管、内套管与油管的单位长度质量,包括立管内部流体、气体质量及附加质量;
c1、c2、c3分别为外套管、内套管与油管的结构阻尼;
(EI)1、(EI)2、(EI)3分别为外套管、内套管与油管的抗弯刚度;
T1、T2、T3分别为外套管、内套管与油管的壁张力,q为波浪荷载;
p2为内套管扶正器与外套管接触时产生的相互作用力;
p3为油管扶正器与内套管接触时产生的相互作用力;
(2)对公式(1)~(3)进行有限元离散,建立各自的有限元方程:
Figure GDA0002815923920000044
Figure GDA0002815923920000045
Figure GDA0002815923920000046
式中:
[Mi]、[Ci]和[Ki]分别为外套管(i=1)、内套管(i=2)和油管(i=3)的质量矩阵、阻尼矩阵和刚度矩阵;
Figure GDA0002815923920000047
(Nd为节点总数)为外套管(i=1)、内套管(i=2)和油管(i=3)的节点位移向量(包括挠度和转角);
Figure GDA00028159239200000511
为内套管(i=2)和油管(i=3)扶正器分别与外套管和内套管的相互作用力;
{q}为节点波浪荷载;
(3)由于扶正器与相邻管柱之间有间隙,其接触状态是时断时续的,因此,方程(4)~(6)的求解需采用迭代法求解,为此,建立(4)~(6)式的增量迭代方程:
Figure GDA0002815923920000051
Figure GDA0002815923920000052
Figure GDA0002815923920000053
(4)对公式(7)~(9)进行耦合迭代分析
第一步:计算外套管当前时刻的响应
由式(7)的Newmark-β法公式
Figure GDA0002815923920000054
计算当前时刻的外套管位移响应增量
Figure GDA0002815923920000055
新增时间步第一次计算时,k=0。
上式中:
Figure GDA0002815923920000056
Figure GDA0002815923920000057
其中,
Figure GDA0002815923920000058
n为公式(7)与公式(8)的耦合迭代次数。
然后,由下式计算外套管的速度增量
Figure GDA0002815923920000059
和加速度增量
Figure GDA00028159239200000510
Figure GDA0002815923920000061
Figure GDA0002815923920000062
如果k=0,计算当前时刻的外套管位移
Figure GDA0002815923920000063
速度
Figure GDA0002815923920000064
和加速度
Figure GDA0002815923920000065
即:
Figure GDA0002815923920000066
Figure GDA0002815923920000067
Figure GDA0002815923920000068
然后,转至第二步。
如果k≠0,则由收敛条件:
Figure GDA0002815923920000069
或:
Figure GDA00028159239200000610
判断迭代是否满足收敛条件,不满足收敛条件则转至第二步的相应接触状态进行计算,否则,进行下一个时间步的计算或结束计算。
第二步:计算内套管和油管当前时刻的响应:
Figure GDA00028159239200000611
Figure GDA00028159239200000612
即时刻t内套管和油管均没有扶正器与外套管或内套管接触,则由下式分别计算内套管和油管的位移响应增量
Figure GDA00028159239200000613
Figure GDA00028159239200000614
Figure GDA00028159239200000615
其中:
Figure GDA00028159239200000616
Figure GDA00028159239200000617
式中,
Figure GDA0002815923920000071
l为公式(8)与公式(9)的耦合迭代次数。
然后,由下式分别计算内套管和油管的速度和加速度增量
Figure GDA0002815923920000072
Figure GDA0002815923920000073
Figure GDA0002815923920000074
及当前时刻的位移、速度和加速度:
Figure GDA0002815923920000075
然后,转至第三步判断当前时刻的接触状态。
Figure GDA0002815923920000076
Figure GDA0002815923920000077
即时刻t仅内套管有扶正器与外套管发生接触,则可由外套管的位移响应增量求出内套管与之接触的扶正器处的位移响应增量:
Figure GDA0002815923920000078
Figure GDA0002815923920000079
Figure GDA00028159239200000710
Figure GDA00028159239200000711
式中:
n(rm)为内套管扶正器rm的节点编号;
m′2为时刻t已接触的扶正器数量,m″2为t+t时刻发生接触的扶正器数量,
Figure GDA00028159239200000712
(k=0时,取
Figure GDA00028159239200000713
)为内套管扶正器与外套管的实时间隙,可按下式计算:
Figure GDA0002815923920000081
式中,δ2为内套管扶正器与外套管的静态间隙,δ2=d1/2-D2/2-h2;其中,d1为外套管内径,D2为内套管外径,h2为内套管扶正器的径向尺寸。式中的正负号由外套管的运动方向确定,向平衡位置运动时(|(a1,j)t+Δt|<|(a1,j)t|)取正号。
由于发生接触时,内套管的响应是由外套管的响应求出的,因此,公式(15)中的内套管响应只能采用前一次迭代的结果。为了避免误差累积,可采用迭代的方法对公式(15)的结果进行修正,即:
Figure GDA0002815923920000082
Figure GDA0002815923920000083
Figure GDA0002815923920000084
求出内套管接触点的挠度增量后,即可由下式计算接触点挠度的速度和加速度增量
Figure GDA0002815923920000085
并由式(8)的Newmark-β法公式
Figure GDA0002815923920000086
求出内套管其它未知位移增量
Figure GDA0002815923920000087
并由公式(11)(i=2)和公式(12)(i=2)求出相应的速度增量
Figure GDA0002815923920000088
和加速度增量
Figure GDA0002815923920000089
及位移
Figure GDA00028159239200000810
速度
Figure GDA00028159239200000811
和加速度
Figure GDA00028159239200000812
公式(19)中,m2为内套管扶正器与外套管发生接触的数量。
式(20)中,
Figure GDA0002815923920000091
为2Nd-m2个元素组成的向量,即不包括发生接触的内套管扶正器rm挠度的位移向量;等效刚度矩阵
Figure GDA0002815923920000092
和等效荷载向量
Figure GDA0002815923920000093
分别为:
Figure GDA0002815923920000094
Figure GDA0002815923920000095
上两式中的
Figure GDA0002815923920000096
Figure GDA0002815923920000097
为不包括j行和j列的质量矩阵、刚度矩阵和阻尼矩阵(j的定义同公式(19)),
Figure GDA0002815923920000098
是由
Figure GDA0002815923920000099
Figure GDA00028159239200000910
与j列刚度系数、阻尼系数和质量系数组成的“荷载”向量:
Figure GDA00028159239200000911
求出
Figure GDA00028159239200000912
Figure GDA00028159239200000913
后,即可由公式(8)计算
Figure GDA00028159239200000914
再返回第一步进行迭代计算。
油管的响应则采用公式(10)~式(12)计算,其中i=3
Figure GDA00028159239200000915
Figure GDA00028159239200000916
即时刻t仅油管有扶正器与内套管发生接触,则由式(10)i=2计算内套管的响应,如果k≠0,则根据收敛条件
Figure GDA00028159239200000917
判断迭代是否收敛。如果满足收敛条件,则转至第一步进行下一个时间步的计算。否则,基于接触条件
Figure GDA00028159239200000918
求出油管接触点的挠度增量,并由下式计算接触点挠度的速度和加速度增量
Figure GDA0002815923920000101
Figure GDA0002815923920000102
式中:n(sm)为油管扶正器sm的节点编号,m3为油管扶正器与内套管发生接触的数量。
然后,由式(9)的Newmark-β法公式
Figure GDA0002815923920000103
求出油管其它未知的位移增量
Figure GDA0002815923920000104
并由公式(11)(i=3)和公式(12)(i=3)求出相应的速度增量
Figure GDA0002815923920000105
和加速度增量
Figure GDA0002815923920000106
及位移
Figure GDA0002815923920000107
速度
Figure GDA0002815923920000108
和加速度
Figure GDA0002815923920000109
式(21)中,
Figure GDA00028159239200001010
为2Nd-m3个元素组成的向量,即不包括发生接触的油管扶正器sm挠度的位移向量;等效刚度矩阵
Figure GDA00028159239200001011
和等效荷载向量
Figure GDA00028159239200001012
分别为:
Figure GDA00028159239200001013
Figure GDA00028159239200001014
上两式中的
Figure GDA00028159239200001015
Figure GDA00028159239200001016
为不包括j行和j列的质量矩阵、刚度矩阵和阻尼矩阵(j=2n1(sm)-1,(m=1,2,…,m3),n1(sm)为发生接触的油管扶正器sm的节点编号),
Figure GDA00028159239200001017
是由
Figure GDA00028159239200001018
Figure GDA00028159239200001019
与j列刚度系数、阻尼系数和质量系数组成的“荷载”向量
Figure GDA00028159239200001020
求出
Figure GDA00028159239200001021
Figure GDA00028159239200001022
后,即可由公式(9)计算
Figure GDA00028159239200001023
再返回③进行迭代计算。
Figure GDA0002815923920000111
Figure GDA0002815923920000112
即时刻t内套管和油管均有扶正器与外套管和内套管发生接触,则根据外套管的运动条件首先由式(13)~式(20)计算内套管的响应,再将式中的下标1改为2、2改为3计算油管的响应,然后,由式(9)计算
Figure GDA0002815923920000113
并代入式(8)计算
Figure GDA0002815923920000114
并返回第一步进行迭代计算。
第三步:判断接触状态
如果
Figure GDA0002815923920000115
成立,则没有新增接触扶正器,转至第一步计算下一时刻的响应。否则,转至第二步的相应接触状态重新计算内套管和油管的响应并进行迭代。
上式中,δi+1为扶正器与相邻管柱的静态间隙,δi+1=di/2-Di+1/2-hi+1;其中,di为外套管(i=1)或内套管(i=2)的内径,Di+1为内套管(i=1)或油管(i=2)外径,hi+1为内套管(i=1)或油管(i=2)扶正器的径向尺寸。
当|(a1,j)t+Δt|>|(a1,j)t|,(j=2n(rm)-1)时,式(22)取正号。
本发明所达到的有益效果为:
本发明在深水顶张式立管随机耦合振动的分析中采用三层立管组成的耦合模型,且充分考虑了管间间隙对耦合响应的影响。计算过程真实反映顶张式立管横向振动的运动和变形状态。且能够较精确地直接计算出外套管、内套管和油管的应力状态及响应,而不必采用截面荷载分配的方法二次计算外套管、内套管和油管的应力状态。
附图说明
图1是本发明三层立管组成的耦合模型示意图。
图中,x轴代表宽度,z轴代表高度。
1、外套管;2、内套管;3、油管;4、扶正器。
具体实施方式
为便于本领域的技术人员理解本发明,下面结合附图说明本发明的具体实施方式。
如图1所示,本发明基于三层立管组成的耦合模型进行计算分析且充分考虑了管间间隙对耦合响应的影响,其具体计算方案如下:
(1)分别建立外套管1、内套管2与油管3的运动方程:
Figure GDA0002815923920000121
Figure GDA0002815923920000122
Figure GDA0002815923920000123
式中:
m1、m2、m3分别为外套管1、内套管2与油管3的单位长度质量,包括立管内部流体、气体质量及附加质量;
c1、c2、c3分别为外套管1、内套管2与油管3的结构阻尼;
(EI)1、(EI)2、(EI)3分别为外套管1、内套管2与油管3的抗弯刚度;
T1、T2、T3分别为外套管1、内套管2与油管3的壁张力,q为波浪荷载;p2为内套管扶正器4与外套管1接触时产生的相互作用力;
p3为油管扶正器4与内套管2接触时产生的相互作用力;
(2)对公式(1)~(3)进行有限元离散,建立各自的有限元方程:
Figure GDA0002815923920000124
Figure GDA0002815923920000125
Figure GDA0002815923920000131
式中:
[Mi]、[Ci]和[Ki]分别为外套管1(i=1)、内套管2(i=2)和油管3(i=3)的质量矩阵、阻尼矩阵和刚度矩阵;
Figure GDA0002815923920000132
(Nd为节点总数)为外套管1(i=1)、内套管2(i=2)和油管3(i=3)的节点位移向量(包括挠度和转角);
Figure GDA0002815923920000138
为内套管(i=2)和油管(i=3)扶正器4分别与外套管1和内套管2的相互作用力;
{q}为节点波浪荷载;
(3)由于扶正器4与相邻管柱之间有间隙,其接触状态是时断时续的,因此,方程(4)~(6)的求解需采用迭代法求解,为此,建立(4)~(6)式的增量迭代方程:
Figure GDA0002815923920000133
Figure GDA0002815923920000134
Figure GDA0002815923920000135
上述方程中考虑了结构的几何非线性及张力变化引起的几何刚度变化,因此,刚度矩阵随时间变化,而阻尼矩阵则是因瑞雷阻尼使然也随时间变化;
(4)对公式(7)~(9)进行耦合迭代分析
第一步:计算外套管1当前时刻的响应
由式(7)的Newmark-β法公式
Figure GDA0002815923920000136
计算当前时刻的外套管1位移响应增量
Figure GDA0002815923920000137
新增时间步第一次计算时,k=0。
上式中:
Figure GDA0002815923920000141
Figure GDA0002815923920000142
其中,
Figure GDA0002815923920000143
n为公式(7)与公式(8)的耦合迭代次数。
然后,由下式计算外套管1的速度增量
Figure GDA0002815923920000144
和加速度增量
Figure GDA0002815923920000145
Figure GDA0002815923920000146
Figure GDA0002815923920000147
如果k=0,计算当前时刻的外套管1位移
Figure GDA0002815923920000148
速度
Figure GDA0002815923920000149
和加速度
Figure GDA00028159239200001410
即:
Figure GDA00028159239200001411
Figure GDA00028159239200001412
Figure GDA00028159239200001413
然后,转至第二步。
如果k≠0,则由收敛条件:
Figure GDA00028159239200001414
或:
Figure GDA00028159239200001415
判断迭代是否满足收敛条件,不满足收敛条件则转至第二步的相应接触状态进行计算,否则,进行下一个时间步的计算或结束计算。
第二步:计算内套管2和油管3当前时刻的响应:
内套管2和油管3的响应取决于扶正器4的接触状态——内套管2和油管3均没有扶正器4与外套管1或内套管2接触
Figure GDA0002815923920000151
Figure GDA0002815923920000152
仅内套管2有扶正器4与外套管1接触
Figure GDA0002815923920000153
Figure GDA0002815923920000154
仅油管3有扶正器4与内套管2接触
Figure GDA0002815923920000155
Figure GDA0002815923920000156
和内套管2和油管3均有扶正器4与外套管1或内套管2接触
Figure GDA0002815923920000157
Figure GDA0002815923920000158
Figure GDA0002815923920000159
Figure GDA00028159239200001510
即时刻t内套管2和油管3均没有扶正器4与外套管1或内套管2接触,则由下式分别计算内套管2和油管3的位移响应增量
Figure GDA00028159239200001511
Figure GDA00028159239200001512
Figure GDA00028159239200001513
其中:
Figure GDA00028159239200001514
Figure GDA00028159239200001515
式中,
Figure GDA00028159239200001516
l为公式(8)与公式(9)的耦合迭代次数。
然后,由下式分别计算内套管2和油管3的速度和加速度增量
Figure GDA00028159239200001517
Figure GDA00028159239200001518
Figure GDA00028159239200001519
Figure GDA00028159239200001520
及当前时刻的位移、速度和加速度:
Figure GDA00028159239200001521
然后,转至第三步判断当前时刻的接触状态。
Figure GDA0002815923920000161
Figure GDA0002815923920000162
即时刻t仅内套管2有扶正器4与外套管1发生接触,则可由外套管1的位移响应增量求出内套管2与之接触的扶正器4处的位移响应增量:
Figure GDA0002815923920000163
Figure GDA0002815923920000164
Figure GDA0002815923920000165
Figure GDA0002815923920000166
式中:
n(rm)为内套管扶正器rm4的节点编号;
m′2为时刻t已接触的扶正器4数量,m″2为t+t时刻发生接触的扶正器4数量,
Figure GDA0002815923920000167
(k=0时,取
Figure GDA0002815923920000168
为内套管扶正器4与外套管1的实时间隙,可按下式计算:
Figure GDA0002815923920000169
式中,δ2为内套管扶正器4与外套管1的静态间隙,δ2=d1/2-D2/2-h2;其中,d1为外套管1内径,D2为内套管2外径,h2为内套管扶正器4的径向尺寸。式中的正负号由外套管1的运动方向确定,向平衡位置运动时
Figure GDA00028159239200001610
取正号。
由于发生接触时,内套管2的响应是由外套管1的响应求出的,因此,公式(15)中的内套管2响应只能采用前一次迭代的结果。为了避免误差累积,可采用迭代的方法对公式(15)的结果进行修正,即:
Figure GDA00028159239200001611
Figure GDA0002815923920000171
Figure GDA0002815923920000172
求出内套管2接触点的挠度增量后,即可由下式计算接触点挠度的速度和加速度增量
Figure GDA0002815923920000173
并由式(8)的Newmark-β法公式
Figure GDA0002815923920000174
求出内套管2其它未知位移增量
Figure GDA0002815923920000175
并由公式(11)(i=2)和公式(12)(i=2)求出相应的速度增量
Figure GDA0002815923920000176
和加速度增量
Figure GDA0002815923920000177
及位移
Figure GDA0002815923920000178
速度
Figure GDA0002815923920000179
和加速度
Figure GDA00028159239200001710
公式(19)中,m2为内套管扶正器4与外套管1发生接触的数量。
式(20)中,
Figure GDA00028159239200001711
为2Nd-m2个元素组成的向量,即不包括发生接触的内套管扶正器rm4挠度的位移向量;等效刚度矩阵
Figure GDA00028159239200001712
和等效荷载向量
Figure GDA00028159239200001713
分别为:
Figure GDA00028159239200001714
Figure GDA00028159239200001715
上两式中的
Figure GDA00028159239200001716
Figure GDA00028159239200001717
为不包括j行和j列的质量矩阵、刚度矩阵和阻尼矩阵(j的定义同公式(19)),
Figure GDA00028159239200001718
是由
Figure GDA00028159239200001719
Figure GDA00028159239200001720
与j列刚度系数、阻尼系数和质量系数组成的“荷载”向量:
Figure GDA0002815923920000181
求出
Figure GDA0002815923920000182
Figure GDA0002815923920000183
后,即可由公式(8)计算
Figure GDA0002815923920000184
再返回第一步进行迭代计算。
油管3的响应则采用公式(10)~式(12)计算,其中i=3
Figure GDA0002815923920000185
Figure GDA0002815923920000186
即时刻t仅油管3有扶正器4与内套管2发生接触,则由式(10)(i=2)计算内套管2的响应,如果k≠0,则根据收敛条件
Figure GDA0002815923920000187
判断迭代是否收敛。如果满足收敛条件,则转至第一步进行下一个时间步的计算。否则,基于接触条件
Figure GDA0002815923920000188
求出油管3接触点的挠度增量,并由下式计算接触点挠度的速度和加速度增量
Figure GDA0002815923920000189
Figure GDA00028159239200001810
式中:n(sm)为油管扶正器sm4的节点编号,m3为油管扶正器4与内套管2发生接触的数量。
然后,由式(9)的Newmark-β法公式
Figure GDA00028159239200001811
求出油管3其它未知的位移增量
Figure GDA00028159239200001812
并由公式(11)(i=3)和公式(12)(i=3)求出相应的速度增量
Figure GDA00028159239200001813
和加速度增量
Figure GDA00028159239200001814
及位移
Figure GDA00028159239200001815
速度
Figure GDA0002815923920000191
和加速度
Figure GDA0002815923920000192
式(21)中,
Figure GDA0002815923920000193
为2Nd-m3个元素组成的向量,即不包括发生接触的油管扶正器sm挠度的位移向量;等效刚度矩阵
Figure GDA0002815923920000194
和等效荷载向量
Figure GDA0002815923920000195
分别为:
Figure GDA0002815923920000196
Figure GDA0002815923920000197
上两式中的
Figure GDA0002815923920000198
Figure GDA0002815923920000199
为不包括j行和j列的质量矩阵、刚度矩阵和阻尼矩阵(j=2n1(sm)-1,(m=1,2,…,m3),n1(sm)为发生接触的油管扶正器sm的节点编号),
Figure GDA00028159239200001910
是由
Figure GDA00028159239200001911
Figure GDA00028159239200001912
与j列刚度系数、阻尼系数和质量系数组成的“荷载”向量
Figure GDA00028159239200001913
求出
Figure GDA00028159239200001914
Figure GDA00028159239200001915
后,即可由公式(9)计算
Figure GDA00028159239200001916
再返回③进行迭代计算。
Figure GDA00028159239200001917
Figure GDA00028159239200001918
即时刻t内套管2和油管3均有扶正器4与外套管1和内套管2发生接触,则根据外套管1的运动条件首先由式(13)~式(20)计算内套管2的响应,再将式中的下标1改为2、2改为3计算油管3的响应,然后,由式(9)计算
Figure GDA00028159239200001919
并代入式(8)计算
Figure GDA00028159239200001920
并返回第一步进行迭代计算。
需要指出的是,在当前时间步有新的扶正器4发生接触时,应确保接触点的过盈量尽可能小,即“正好”发生接触。这不仅是静态接触假定的需要,而且是计算收敛的前提。如果过盈量较大,可能导致死循环。因此,必须调整当前时间步长Δt,使扶正器4与相邻管柱处于非接触向接触过渡的临界点。正如材料非线性问题中,当有单元进入屈服时,必须调整荷载增量,以使进入屈服的单元处于屈服的临界点。
第三步:判断接触状态
如果
Figure GDA0002815923920000201
成立,则没有新增接触扶正器4,转至第一步计算下一时刻的响应。否则,转至第二步的相应接触状态重新计算内套管2和油管3的响应并进行迭代。
上式中,δi+1为扶正器4与相邻管柱的静态间隙,δi+1=di/2-Di+1/2-hi+1;其中,di为外套管1(i=1)或内套管2(i=2)的内径,Di+1为内套管2(i=1)或油管3(i=2)外径,hi+1为内套管2(i=1)或油管3(i=2)扶正器4的径向尺寸。当
Figure GDA0002815923920000202
时,式(22)取正号。

Claims (1)

1.一种顶张式立管随机耦合振动的计算方法,其特征在于:该方法基于三层立管组成的耦合模型进行计算分析且充分考虑了管间间隙对耦合响应的影响,其具体计算方案如下:
(1)分别建立外套管、内套管与油管的运动方程:
Figure FDA0002815923910000011
Figure FDA0002815923910000012
Figure FDA0002815923910000013
式中:
m1、m2、m3分别为外套管、内套管与油管的单位长度质量,包括立管内部流体、气体质量及附加质量;
c1、c2、c3分别为外套管、内套管与油管的结构阻尼;
(EI)1、(EI)2、(EI)3分别为外套管、内套管与油管的抗弯刚度;
T1、T2、T3分别为外套管、内套管与油管的壁张力,q为波浪荷载;
p2为内套管扶正器与外套管接触时产生的相互作用力;
p3为油管扶正器与内套管接触时产生的相互作用力;
(2)对公式(1)~(3)进行有限元离散,建立各自的有限元方程:
Figure FDA0002815923910000014
Figure FDA0002815923910000015
Figure FDA0002815923910000016
式中:
[Mi]、[Ci]和[Ki]分别为外套管i=1、内套管i=2和油管i=3的质量矩阵、阻尼矩阵和刚度矩阵;
Figure FDA0002815923910000021
为外套管i=1、内套管i=2和油管i=3的节点位移向量,所述节点位移向量包括挠度和转角,其中Nd为节点总数;
Figure FDA0002815923910000022
为内套管i=2和油管i=3扶正器分别与外套管和内套管的相互作用力;
{q}为节点波浪荷载;
(3)由于扶正器与相邻管柱之间有间隙,其接触状态是时断时续的,因此,方程(4)~(6)的求解需采用迭代法求解,为此,建立(4)~(6)式的增量迭代方程:
Figure FDA0002815923910000023
Figure FDA0002815923910000024
Figure FDA0002815923910000025
(4)对公式(7)~(9)进行耦合迭代分析
第一步:计算外套管当前时刻的响应
由式(7)的Newmark-β法公式
Figure FDA0002815923910000026
计算当前时刻的外套管位移响应增量
Figure FDA0002815923910000027
新增时间步第一次计算时,k=0;
上式中:
Figure FDA0002815923910000028
Figure FDA0002815923910000029
其中,
Figure FDA00028159239100000210
n为公式(7)与公式(8)的耦合迭代次数;
然后,由下式计算外套管的速度增量
Figure FDA0002815923910000031
和加速度增量
Figure FDA0002815923910000032
Figure FDA0002815923910000033
Figure FDA0002815923910000034
如果k=0,计算当前时刻的外套管位移
Figure FDA0002815923910000035
速度
Figure FDA0002815923910000036
和加速度
Figure FDA0002815923910000037
即:
Figure FDA0002815923910000038
Figure FDA0002815923910000039
Figure FDA00028159239100000310
然后,转至第二步;
如果k≠0,则由收敛条件:
Figure FDA00028159239100000311
或:
Figure FDA00028159239100000312
判断迭代是否满足收敛条件,不满足收敛条件则转至第二步的相应接触状态进行计算,否则,进行下一个时间步的计算或结束计算;
第二步:计算内套管和油管当前时刻的响应:
Figure FDA00028159239100000313
Figure FDA00028159239100000314
即时刻t内套管和油管均没有扶正器与外套管或内套管接触,则由下式分别计算内套管和油管的位移响应增量
Figure FDA00028159239100000315
Figure FDA00028159239100000316
Figure FDA00028159239100000317
其中:
Figure FDA00028159239100000318
Figure FDA0002815923910000041
式中,
Figure FDA0002815923910000042
l为公式(8)与公式(9)的耦合迭代次数;
然后,由下式分别计算内套管和油管的速度和加速度增量
Figure FDA0002815923910000043
Figure FDA0002815923910000044
Figure FDA0002815923910000045
及当前时刻的位移、速度和加速度:
Figure FDA0002815923910000046
然后,转至第三步判断当前时刻的接触状态;
Figure FDA0002815923910000047
Figure FDA0002815923910000048
即时刻t仅内套管有扶正器与外套管发生接触,则可由外套管的位移响应增量求出内套管与之接触的扶正器处的位移响应增量:
Figure FDA0002815923910000049
Figure FDA00028159239100000410
Figure FDA00028159239100000411
Figure FDA00028159239100000412
式中:
n(rm)为内套管扶正器rm的节点编号;
m′2为时刻t已接触的扶正器数量,m″2为t+△t时刻发生接触的扶正器数量,
Figure FDA00028159239100000413
为内套管扶正器与外套管的实时间隙,可按下式计算:
Figure FDA0002815923910000051
当k=0时,取
Figure FDA0002815923910000052
式中,δ2为内套管扶正器与外套管的静态间隙,δ2=d1/2-D2/2-h2;其中,d1为外套管内径,D2为内套管外径,h2为内套管扶正器的径向尺寸;式中的正负号由外套管的运动方向确定,向平衡位置运动,即|(a1,j)t+Δt|<|(a1,j)t|时,式中的正负号取正号;
由于发生接触时,内套管的响应是由外套管的响应求出的,因此,公式(15)中的内套管响应只能采用前一次迭代的结果;为了避免误差累积,可采用迭代的方法对公式(15)的结果进行修正,即:
Figure FDA0002815923910000053
Figure FDA0002815923910000054
Figure FDA0002815923910000055
求出内套管接触点的挠度增量后,即可由下式计算接触点挠度的速度和加速度增量
Figure FDA0002815923910000056
并由式(8)的Newmark-β法公式
Figure FDA0002815923910000057
求出内套管其它未知位移增量
Figure FDA0002815923910000058
并由公式(11)i=2和公式(12)i=2求出相应的速度增量
Figure FDA0002815923910000059
和加速度增量
Figure FDA00028159239100000510
及位移
Figure FDA00028159239100000511
速度
Figure FDA00028159239100000512
和加速度
Figure FDA0002815923910000061
公式(19)中,m2为内套管扶正器与外套管发生接触的数量;
式(20)中,
Figure FDA0002815923910000062
为2Nd-m2个元素组成的向量,即不包括发生接触的内套管扶正器rm挠度的位移向量;等效刚度矩阵
Figure FDA0002815923910000063
和等效荷载向量
Figure FDA0002815923910000064
分别为:
Figure FDA0002815923910000065
Figure FDA0002815923910000066
上两式中的
Figure FDA0002815923910000067
Figure FDA0002815923910000068
为不包括j行和j列的质量矩阵、刚度矩阵和阻尼矩阵,其中j的定义与公式(19)中j的定义相同,
Figure FDA0002815923910000069
是由
Figure FDA00028159239100000610
Figure FDA00028159239100000611
Figure FDA00028159239100000612
与j列刚度系数、阻尼系数和质量系数组成的“荷载”向量:
Figure FDA00028159239100000613
求出
Figure FDA00028159239100000614
Figure FDA00028159239100000615
后,即可由公式(8)计算
Figure FDA00028159239100000616
再返回第一步进行迭代计算;
油管的响应则采用公式(10)~式(12)计算,其中i=3
Figure FDA00028159239100000617
Figure FDA00028159239100000618
即时刻t仅油管有扶正器与内套管发生接触,则由公式(10)i=2计算内套管的响应,如果k≠0,则根据收敛条件
Figure FDA00028159239100000619
判断迭代是否收敛;如果满足收敛条件,则转至第一步进行下一个时间步的计算;否则,基于接触条件
Figure FDA00028159239100000620
求出油管接触点的挠度增量,并由下式计算接触点挠度的速度和加速度增量
Figure FDA0002815923910000071
Figure FDA0002815923910000072
式中:n(sm)为油管扶正器sm的节点编号,m3为油管扶正器与内套管发生接触的数量;
然后,由式(9)的Newmark-β法公式
Figure FDA0002815923910000073
求出油管其它未知的位移增量
Figure FDA0002815923910000074
并由公式(11)i=3和公式(12)i=3求出相应的速度增量
Figure FDA0002815923910000075
和加速度增量
Figure FDA0002815923910000076
及位移
Figure FDA0002815923910000077
速度
Figure FDA0002815923910000078
和加速度
Figure FDA0002815923910000079
式(21)中,
Figure FDA00028159239100000710
为2Nd-m3个元素组成的向量,即不包括发生接触的油管扶正器sm挠度的位移向量;等效刚度矩阵
Figure FDA00028159239100000711
和等效荷载向量
Figure FDA00028159239100000712
分别为:
Figure FDA00028159239100000713
Figure FDA00028159239100000714
上两式中的
Figure FDA00028159239100000715
Figure FDA00028159239100000716
为不包括j行和j列的质量矩阵、刚度矩阵和阻尼矩阵,其中j=2n1(sm)-1,(m=1,2,…,m3),n1(sm)为发生接触的油管扶正器sm的节点编号,
Figure FDA00028159239100000717
是由
Figure FDA00028159239100000718
Figure FDA00028159239100000719
与j列刚度系数、阻尼系数和质量系数组成的“荷载”向量
Figure FDA00028159239100000720
求出
Figure FDA0002815923910000081
Figure FDA0002815923910000082
后,即可由公式(9)计算
Figure FDA0002815923910000083
再返回③进行迭代计算;
Figure FDA0002815923910000084
Figure FDA0002815923910000085
即时刻t内套管和油管均有扶正器与外套管和内套管发生接触,则根据外套管的运动条件首先由式(13)~式(20)计算内套管的响应,再将式中的下标1改为2、2改为3计算油管的响应,然后,由式(9)计算
Figure FDA0002815923910000086
并代入式(8)计算
Figure FDA0002815923910000087
并返回第一步进行迭代计算;
第三步:判断接触状态
如果
Figure FDA0002815923910000088
成立,则没有新增接触扶正器,转至第一步计算下一时刻的响应;否则,转至第二步的相应接触状态重新计算内套管和油管的响应并进行迭代;
上式中,δi+1为扶正器与相邻管柱的静态间隙,δi+1=di/2-Di+1/2-hi+1;其中,di为外套管i=1或内套管i=2的内径,Di+1为内套管i=1或油管i=2外径,hi+1为内套管i=1或油管i=2扶正器的径向尺寸;
当|(a1,j)t+Δt|>|(a1,j)t|,(j=2n(rm)-1)时,式(22)取正号。
CN201910682572.5A 2019-07-26 2019-07-26 一种顶张式立管随机耦合振动的计算方法 Expired - Fee Related CN110378060B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910682572.5A CN110378060B (zh) 2019-07-26 2019-07-26 一种顶张式立管随机耦合振动的计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910682572.5A CN110378060B (zh) 2019-07-26 2019-07-26 一种顶张式立管随机耦合振动的计算方法

Publications (2)

Publication Number Publication Date
CN110378060A CN110378060A (zh) 2019-10-25
CN110378060B true CN110378060B (zh) 2021-02-09

Family

ID=68256450

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910682572.5A Expired - Fee Related CN110378060B (zh) 2019-07-26 2019-07-26 一种顶张式立管随机耦合振动的计算方法

Country Status (1)

Country Link
CN (1) CN110378060B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110795790B (zh) * 2019-10-31 2021-02-12 李鲁 一种复杂建筑结构非线性动力时程分析方法
TWI764312B (zh) 2020-10-08 2022-05-11 國立中央大學 基於等效節點割線質量逼近之結構體分析方法、裝置與電腦程式產品
CN112434472B (zh) * 2020-10-13 2024-04-05 华北电力大学 一种反应堆多层同轴筒体窄缝间隙附加质量计算方法
CN112834164B (zh) * 2020-12-31 2021-12-21 中国海洋大学 一种考虑约化速度和间距的尾流立管涡激升力确定方法
CN113094918B (zh) * 2021-04-22 2022-04-12 中国海洋大学 一种海洋深水湿式保温管设计外压的确定方法

Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101539477A (zh) * 2009-05-08 2009-09-23 中国海洋大学 一种深水顶张式立管涡激振动与疲劳分析的方法
CN102252897A (zh) * 2011-06-16 2011-11-23 中国海洋大学 一种深水顶张式立管弯曲振动分析方法
CN102353506A (zh) * 2011-06-16 2012-02-15 中国海洋大学 一种深水顶张式立管竖向振动分析方法
EP2886788A2 (en) * 2013-12-23 2015-06-24 2HOffshore, Inc. Riser fatigue monitoring
WO2016090217A1 (en) * 2014-12-05 2016-06-09 Schlumberger Canada Limited Monitoring tubing related equipment
CN106294927A (zh) * 2016-07-21 2017-01-04 中国海洋大学 一种载流管道动力特性的计算方法
CN106909710A (zh) * 2017-01-11 2017-06-30 中国海洋大学 深水顶张式立管全耦合动力分析方法
CN108491615A (zh) * 2018-03-17 2018-09-04 中国海洋大学 一种三层顶张式立管动力响应的有限元分析方法
WO2018215653A1 (en) * 2017-05-26 2018-11-29 Universitat Rovira I Virgili Device for passive suppression of vortex-induced vibrations (viv) in structures
CN109827734A (zh) * 2019-03-29 2019-05-31 西安建筑科技大学 一种评估内外流作用下深海立管涡激振动的方法
CN110008635A (zh) * 2019-04-19 2019-07-12 陕西新西商工程科技有限公司 利用Newmark精细积分结合法对弹塑性结构地震反应分析的方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102222168A (zh) * 2011-06-16 2011-10-19 中国海洋大学 一种深水钻井立管参激横向振动分析方法
CN109409006B (zh) * 2018-11-15 2022-12-20 中国地震局工程力学研究所 一种超高层结构动力时程分析方法

Patent Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101539477A (zh) * 2009-05-08 2009-09-23 中国海洋大学 一种深水顶张式立管涡激振动与疲劳分析的方法
CN102252897A (zh) * 2011-06-16 2011-11-23 中国海洋大学 一种深水顶张式立管弯曲振动分析方法
CN102353506A (zh) * 2011-06-16 2012-02-15 中国海洋大学 一种深水顶张式立管竖向振动分析方法
EP2886788A2 (en) * 2013-12-23 2015-06-24 2HOffshore, Inc. Riser fatigue monitoring
WO2016090217A1 (en) * 2014-12-05 2016-06-09 Schlumberger Canada Limited Monitoring tubing related equipment
CN106294927A (zh) * 2016-07-21 2017-01-04 中国海洋大学 一种载流管道动力特性的计算方法
CN106909710A (zh) * 2017-01-11 2017-06-30 中国海洋大学 深水顶张式立管全耦合动力分析方法
WO2018215653A1 (en) * 2017-05-26 2018-11-29 Universitat Rovira I Virgili Device for passive suppression of vortex-induced vibrations (viv) in structures
CN108491615A (zh) * 2018-03-17 2018-09-04 中国海洋大学 一种三层顶张式立管动力响应的有限元分析方法
CN109827734A (zh) * 2019-03-29 2019-05-31 西安建筑科技大学 一种评估内外流作用下深海立管涡激振动的方法
CN110008635A (zh) * 2019-04-19 2019-07-12 陕西新西商工程科技有限公司 利用Newmark精细积分结合法对弹塑性结构地震反应分析的方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
Riser-soil interaction model effects on the dynamic behavior of a steel catenary riser;Xinglan Bai et.al;《Marine Structures》;20150121;第41卷;53-76 *
基于惯性耦合的深水钢悬链线立管非线性分析方法研究;白兴兰;《中国博士学位论文全文数据库工程科技Ⅱ辑(月刊)》;20091115(第11期);C034-8 *
应用于干式半潜平台的TTR设计与结构分析;张嘉;《中国优秀硕士学位论文全文数据库工程科技Ⅱ辑(月刊)》;20180715(第07期);C036-161 *
深水顶张式立管全耦合动力响应及参数敏感性分析;罗坤洪;《海洋工程》;20190131;第37卷(第1期);108-116 *
深海顶张力立管参激_涡激耦合动力响应分析;邵卫东;《中国优秀硕士学位论文全文数据库基础科学辑(月刊)》;20120815(第08期);A010-7 *

Also Published As

Publication number Publication date
CN110378060A (zh) 2019-10-25

Similar Documents

Publication Publication Date Title
CN110378060B (zh) 一种顶张式立管随机耦合振动的计算方法
Arjomandi et al. Elastic buckling capacity of bonded and unbonded sandwich pipes under external hydrostatic pressure
CN103967428B (zh) 一种钻柱疲劳失效风险的评价方法
Xu et al. Collapse analyses of sandwich pipes under external pressure considering inter-layer adhesion behaviour
Arjomandi et al. Bending capacity of sandwich pipes
CN106909710B (zh) 深水顶张式立管全耦合动力分析方法
Sheng et al. Theory and model implementation for analyzing line structures subject to dynamic motions of large deformation and elongation using the absolute nodal coordinate formulation (ANCF) approach
Karamanos et al. Tubular members. I: Stability analysis and preliminary results
Chen et al. An extraction of the natural frequencies and mode shapes of marine risers by the method of differential transformation
Wang et al. A nonlinear dynamic model for 2D deepwater steel lazy-wave riser subjected to top–end imposed excitations
CN106503399B (zh) 垂直井悬挂管柱螺旋屈曲临界载荷的确定方法
Falser et al. Interaction between a compliant guide and a coiled tubing during sub-sea well intervention in deep water
Gu et al. Three-dimensional dynamic analysis of deep-water steel steep wave riser considering internal solitary wave
Chang et al. Axial-transverse coupled vortex-induced vibration characteristics of composite riser under gas-liquid two-phase internal flow
CN108491615B (zh) 一种三层顶张式立管动力响应的有限元分析方法
Chatjigeorgiou Second-order nonlinear dynamics of catenary pipelines: A frequency domain approach
Sun et al. The simulation model of sucker rod string transverse vibration under the space buckling deformation excitation and rod-tubing eccentric wear in vertical wells
Qin et al. Quasistatic nonlinear analysis of a drill pipe in subsea xmas tree installation
Yun et al. Improvement of the bending behavior of a flexible riser: Part I–nonlinear bending behavior considering the shear deformation of polymer layers
Chen et al. Modeling approach of hydropneumatic tensioner for top-tensioned riser
Luo et al. Numerical study on the effect of low Reynolds number flows in straight tube Coriolis flowmeters
CN110727977B (zh) 一种顶张式立管耦合运动的数值模拟方法
Zhao et al. Dynamic stability of a stepped drillstring conveying drilling fluid
Wan et al. Numerical model and program development of horizontal directional drilling for non-excavation technology
Wang et al. Experimental study of drilling riser and wellhead force by small scale testing

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: 20210209

Termination date: 20210726

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