CN111985138A - 一种柔性结构横流与顺流方向涡激振动耦合响应预测方法 - Google Patents

一种柔性结构横流与顺流方向涡激振动耦合响应预测方法 Download PDF

Info

Publication number
CN111985138A
CN111985138A CN202010852497.5A CN202010852497A CN111985138A CN 111985138 A CN111985138 A CN 111985138A CN 202010852497 A CN202010852497 A CN 202010852497A CN 111985138 A CN111985138 A CN 111985138A
Authority
CN
China
Prior art keywords
vibration
coefficient
dimensionless
formula
equations
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
CN202010852497.5A
Other languages
English (en)
Other versions
CN111985138B (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.)
Harbin Institute of Technology Weihai
Original Assignee
Harbin Institute of Technology Weihai
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 Harbin Institute of Technology Weihai filed Critical Harbin Institute of Technology Weihai
Priority to CN202010852497.5A priority Critical patent/CN111985138B/zh
Publication of CN111985138A publication Critical patent/CN111985138A/zh
Application granted granted Critical
Publication of CN111985138B publication Critical patent/CN111985138B/zh
Active 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
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
    • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Business, Economics & Management (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Strategic Management (AREA)
  • Human Resources & Organizations (AREA)
  • Economics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Development Economics (AREA)
  • Fluid Mechanics (AREA)
  • Game Theory and Decision Science (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Mathematical Physics (AREA)
  • Marketing (AREA)
  • Operations Research (AREA)
  • Quality & Reliability (AREA)
  • Tourism & Hospitality (AREA)
  • General Business, Economics & Management (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
  • Bridges Or Land Bridges (AREA)

Abstract

一种柔性结构横流与顺流方向涡激振动耦合响应预测方法。目前缺少对柔性结构横流与顺流方向的涡激振动耦合响应预测措施，导致横流与顺流方向上的耦合VIV特性难以准确掌握，影响柔性圆柱相关研究结果准确性。本发明的涡激振动耦合响应预测方法为通过分别建立IL方向与CF方向振动方程，基于有限差分法对建立IL方向与CF方向振动方程进行求解，得出计算结果结合无量纲基本参数，评价IL方向VIV振动位移响应特性以及CF方向VIV振动位移响应特性。本发明用于海洋工程领域中。

Description

一种柔性结构横流与顺流方向涡激振动耦合响应预测方法
技术领域
本发明涉及一种涡激振动耦合响应预测方法，属于电数字数据处理技术领域。
背景技术
当柔性圆柱体结构在一定的来流作用下时，会在结构尾部形成交替脱落的漩涡，周期性的漩涡脱落会在结构上产生周期性的水动力载荷，从而诱发结构发生振动，称为“涡激振动”。VIV是一个典型的流固耦合问题，其结构振动是由旋涡脱落引起的，而振动结构会反过来影响结构后面的流场，从而导致施加在结构上的水动力发生变化。在海洋工程领域的研究过程中，VIV对采油立管、系泊缆绳和管道的动态特性的影响热点且重点的研究问题，由于VIV可能会增加施加在结构上的动载荷，并可能进一步导致结构的疲劳损伤，所以结构的VIV预测是结构设计过程中的关键问题之一。
近年来，许多学者对圆柱的VIV问题进行了研究。根据研究方法，这些研究可分为物理实验研究和数值研究两大类。采用实验方法的研究在内容和结果上大体是丰富和可靠的。但是，实验方法存在一定的局限性，如模型规模限制、成本昂贵等。数值方法可以克服其中的一些限制。例如，与实验方法相比，成本可以大大降低。两种基于确定作用在结构上的流体力的数值方法得到了广泛的应用。一种方法是计算流体力学(CFD)，它通过求解Navier-Stokes方程直接计算流体力。另一种方法是尾流振子法，它使用非线性振子通过经验来确定流体力。与CFD方法相比，尾流振子方法的明显优势在于其计算成本小。因此，尾流振子法在圆柱，特别是柔性长圆柱的VIV预测中得到了广泛的应用。
利用尾流振子模型对柔性圆柱进行的数值研究，根据结构振动方向可分为两类：一类是只考虑横流方向(CF方向)上的单向振动，另一类是同时考虑横流方向(CF方向)和顺流方向(IL方向)上的耦合振动。对于仅考虑横流方向(CF方向)的振动问题的研究相对较多，而对同时考虑横流方向(CF方向)和顺流方向(IL方向)两个方向的振动问题的研究较少，圆柱体结构在横流方向(CF方向)和顺流方向(IL方向)的两个方向上的耦合VIV特性仍不清楚，由于圆柱体结构在横流方向(CF方向)和顺流方向(IL方向)的两个方向上的耦合VIV特性未有明确结果，从而影响柔性圆柱相关研究结果准确性，在实际情况中遇到相关问题的处理方式不够准确全面。
发明内容：
针对上述问题，本发明公开了一种柔性结构横流与顺流方向涡激振动耦合响应预测方法。
本发明所采用的技术方案为：
一种柔性结构横流与顺流方向涡激振动耦合响应预测方法，所述涡激振动耦合响应预测方法为通过分别建立IL方向与CF方向振动方程，基于有限差分法对建立IL方向与CF方向振动方程进行求解，得出计算结果结合无量纲基本参数，评价IL方向VIV振动位移响应特性以及CF方向VIV振动位移响应特性。
本发明的有益效果为：
一、本发明专用于柔性结构横流与顺流方向的涡激振动耦合响应的计算分析过程，为研究柔性结构的横流方向与顺流方向耦合振动响应进行了研究，建立了完整的横流以及顺流方向的耦合振动模型，用于分析同时考虑横流方向以及顺流方向的耦合动力响应特性，该模型的选用能够模拟出结构的振动位移时间历程曲线、振动频率、振动轨迹以及振动响应位移包络线，结果全面准确且可靠。
二、本发明需解的方程是细长张力梁连续系统振动方程，本发明在计算过程中通过以单位长度平均拖曳力、单位长度振荡拖曳力、单位长度升力为中间推导计算指标进行合理推导，确保预测结论的准确性。
三、本发明在评价IL方向VIV振动位移响应特性以及CF方向VIV振动位移响应特性过程中，使用边界条件为圆柱体两端附近四个特定位置，四个特定位置分别为m＝0、1、m-1和m处对应的
Figure BDA0002645197110000031
和
Figure BDA0002645197110000032
的值，从而确保计算过程的全面可靠性。
四、本发明的研究分析对象柔性结构，具体为柔性圆柱体。本发明中振动位移响应不仅与时间t有关，还与轴线上的空间位置z有关，因此本发明中解的方程为偏微分方程。
五、本发明中的振动方程为4阶偏微分方程，结合二阶精度差分法解微分方程时，计算步骤复杂、合理且准确。
六、本发明采用张力和弯曲刚度联合的方向来考虑弹性恢复力，采取张力和弯曲刚度联合的方式使研究弹性恢复力的准确性更高。
附图说明：
为了易于说明，本发明由下述的具体实施及附图作以详细描述。
图1为均匀来流下两端铰接柔性圆柱体振动模型示意图；
图2为结构细长比L/D为100时的CF方向柔性圆柱体振动位移包络线以及振动位移响应变化云图；
图3为结构细长比L/D为200时的CF方向柔性圆柱体振动位移包络线以及振动位移响应变化云图；
图4为结构细长比L/D为500时的CF方向柔性圆柱体振动位移包络线以及振动位移响应变化云图；
图5为结构细长比L/D为1000时的CF方向柔性圆柱体振动位移包络线以及振动位移响应变化云图；
图6为结构细长比L/D为100时的IL方向柔性圆柱体振动位移包络线以及振动位移响应变化云图；
图7为结构细长比L/D为200时的IL方向柔性圆柱体振动位移包络线以及振动位移响应变化云图；
图8为结构细长比L/D为500时的IL方向柔性圆柱体振动位移包络线以及振动位移响应变化云图；
图9为结构细长比L/D为1000时的IL方向柔性圆柱体振动位移包络线以及振动位移响应变化云图。
具体实施方式：
为使本发明的目的、技术方案和优点更加清楚明了，以振动台混合试验原理为基础，说明应用本发明方法开展振动台混合试验的基本原理，但是应该理解，这些描述只是示例性的，而并非要限制本发明的范围。此外，在以下说明中，省略了对公知结构和技术的描述，以避免不必要地混淆本发明的概念。
具体实施方式一：结合图1、图2、图3、图4、图5、图6、图7、图8和图9说明本实施方式，本实施方式中所述涡激振动耦合响应预测方法为通过分别建立IL方向与CF方向振动方程，基于有限差分法对建立IL方向与CF方向振动方程进行求解，得出计算结果结合无量纲基本参数，评价IL方向的VIV振动位移响应特性以及CF方向的VIV振动位移响应特性。
如图1所示，取一长度为L、直径为D的柔性圆柱体在均匀来流U作用下引起的CF方向以及IL方向互为耦合的涡激振动响应问题，柔性圆柱体两端采用铰接边界条件，坐标系的原点O位于圆柱体的底端，其中X方向为IL方向，Y方向为CF方向，Z方向为铅直方向，柔性圆柱体上的张力记为Θ，柔性圆柱体上的弯曲刚度记为EI。
本实施方式中将柔性圆柱体作为张力梁模型，分别建立张力梁模型的IL方向振动方程以及CF方向振动方程：
Figure BDA0002645197110000041
Figure BDA0002645197110000042
在公式(1)和公式(2)中，m为振动系统单位长度的质量，R为阻尼系数，T为时间，Fx(Z,T)为X方向尾流动力学引起的单位长度的外部水动力激励力，Fy(Z,T)为Y方向尾流动力学引起的单位长度的外部水动力激励力；阻尼系数R包括流体阻尼系数Rf，流体阻尼系数Rf＝γΩfρD2＝(2πStU/D)γρD2；上式中Ωf为漩涡脱落频率，St为斯脱哈尔数，γ为黏滞力系数，斯脱哈尔数St、黏滞力系数γ与流体阻力系数
Figure BDA0002645197110000051
的关系式为：
Figure BDA0002645197110000052
结合公式(1)、公式(2)、流体阻尼系数Rf以及斯脱哈尔数St、黏滞力系数γ与流体阻力系数
Figure BDA0002645197110000053
的关系式得到Fx(Z,T)和Fy(Z,T)表达式为：
Figure BDA0002645197110000054
Figure BDA0002645197110000055
公式(3)和公式(4)中
Figure BDA0002645197110000056
为单位长度平均拖曳力，FD(Z,T)为单位长度振荡拖曳力，FL(Z,T)为单位长度升力，单位长度平均拖曳力
Figure BDA0002645197110000057
单位长度振荡拖曳力FD(Z,T)和单位长度升力FL(Z,T)分别表示为:
Figure BDA0002645197110000058
公式(5)中
Figure BDA0002645197110000059
为平均拖曳力系数且其为常数，CD(Z,T)为振荡拖曳力系数，CL(Z,T)为升力系数；
振荡拖曳力系数CD(Z,T)的表达式为CD(Z,T)＝CD0·p(Z,T)/2；
升力系数CL(Z,T)表示为CL(Z,T)＝CL0·q(Z,T)/2；
上式中CD0为静止圆柱体的振荡拖曳力系数，CL0为静止圆柱体的升力系数，p(Z,T)为与结构上的振荡拖曳力系数有关的无量纲尾流变量，q(Z,T)为与结构上的升力系数有关的无量纲尾流变量；
采用改进的Van der pol方程来满足尾流振子的非线性特性，表达式为：
Figure BDA00026451971100000510
Figure BDA00026451971100000511
公式(6)和公式(7)中εx、εy、Ax以及Ay分别为第一经验参数、第二经验参数、第三经验参数和第四经验参数，x和y分别是IL和CF方向上的无量纲振动振幅，将公式(1)、(2)、(6)和(7)转化为无量纲形式，配合使用的表达式为：
x＝X/D,y＝Y/D,z＝Z/D,t＝T·Ωf (8)
公式(8)中t和z分别是无量纲时间和无量纲空间位置，x和y分别是IL和CF方向上的无量纲振动振幅，将公式(8)带入公式(1)、(2)、(6)和(7)中，整理得到IL方向以及CF方向结构和尾流振子的无量纲方程为为：
Figure BDA0002645197110000061
Figure BDA0002645197110000062
Figure BDA0002645197110000063
Figure BDA0002645197110000064
公式(9)、(10)、(11)和(12)中，质量比μ，无量纲系统质量参数
Figure BDA0002645197110000065
MD和ML，无量纲张力c和无量纲弯曲刚度b的表达式分别为：
Figure BDA0002645197110000066
本实施方式中基于有限差分法对建立IL方向与CF方向振动方程进行求解的过程为在时间和空间上采用标准二阶精度中心差分格式对公式(9)～(12)进行先离散后迭代求解过程，L为柔性圆柱体的长度，D为柔性圆柱体的直径，将结构无量纲总长度L/D划分为M段；将无量纲总时间ttotal划分为N段，从而数值计算时空间步长Δz＝L/(D×M)，时间步长Δt＝ttotal/N；
被划分后的M+1个空间点记为：z＝zi(i＝0,1,2,…,M)；
被划分后的N+1时间点记为：t＝tj(j＝0,1,2,..,N)；
当tn时刻zm位置处参数x、y、p和q表示为
Figure BDA0002645197110000067
以及
Figure BDA0002645197110000068
时，公式(9)～(12)中各偏导数项的二阶精度差分格式表达式分别为：
Figure BDA0002645197110000071
Figure BDA0002645197110000072
Figure BDA0002645197110000073
Figure BDA0002645197110000074
将公式(14)～(17)代入公式(9)～(12)整理得到：
Figure BDA0002645197110000075
Figure BDA0002645197110000076
Figure BDA0002645197110000077
Figure BDA0002645197110000078
x、y的初始条件设为在整个轴线上柔性圆柱体振动位移以及速度均为0，即：
Figure BDA0002645197110000079
p和q的初始条件设为：p和q均为一微小的振幅以及
Figure BDA00026451971100000710
将
Figure BDA00026451971100000711
和
Figure BDA00026451971100000712
的表达式代入公式(14)得到：
Figure BDA00026451971100000713
将公式(22)代入公式(18)～(21)中得到t1时刻x、y、p以及q的值，表达式分别为：
Figure BDA0002645197110000081
至此，得出n＝0以及n＝1时刻的
Figure BDA0002645197110000082
以及
Figure BDA0002645197110000083
的值，由公式(18)～(21)可知，当n≥2时，无需边界条件即可直接获得2≤m≤M-2位置处的
Figure BDA0002645197110000084
的值，同时必须使用边界条件确定柔性圆柱体两端附近四个特定位置，四个特定位置分别为m＝0、1、m-1和m处的
Figure BDA0002645197110000085
和
Figure BDA0002645197110000086
的值。
本实施方式中作为张力梁模型的柔性圆柱体的两端均采用铰接边界条件，即x和y方向的位移和弯矩始终为零，表达式为：
Figure BDA0002645197110000087
当m＝0以及m＝M时，则需结合位移边界条件对其加以求解，由两端处位移为0得到：
Figure BDA0002645197110000088
当m＝1以及m＝M-1时，则需结合弯矩为0边界条件，由两端弯矩为0得到：
Figure BDA0002645197110000089
Figure BDA00026451971100000810
将公式(27)代入公式(20)中得到m＝1以及m＝M-1时y的表达式如下：
Figure BDA00026451971100000811
Figure BDA0002645197110000091
将公式(26)代入公式(18)再结合式(28)和(29)中求得的
Figure BDA0002645197110000092
以及
Figure BDA0002645197110000093
得到m＝1以及m＝M-1时x的表达式如下：
Figure BDA0002645197110000094
Figure BDA0002645197110000095
至此得到n+1时刻整个轴线上(0≤m≤M)所有位置的振动位移
Figure BDA0002645197110000096
和
Figure BDA0002645197110000097
将计算得到的
Figure BDA0002645197110000098
和
Figure BDA0002645197110000099
代入公式(19)以及(21)即可得到n+1时刻的
Figure BDA00026451971100000910
和
Figure BDA00026451971100000911
的值，依次类推对公式(18)～(21)进行反复迭代求解便可得到整个计算时间域内的x、y、p以及q的值。
本实施方式中计算结果结合无量纲基本参数，评价IL方向VIV振动位移响应特性以及CF方向VIV振动位移响应特性的过程为基于无量纲基本参数数据，对计算得出
Figure BDA00026451971100000912
和
Figure BDA00026451971100000913
值进行分析计算的过程。
具体实施方式二：本实施方式为具体实施方式一的进一步限定，阻尼系数R还包括结构阻尼系数Rs，即阻尼系数R包括结构阻尼系数Rs和流体阻尼系数Rf，流体阻尼系数Rf表示为：Rf＝γΩfρD2＝(2πStU/D)γρD2，其中Ωf为漩涡脱落频率；St为斯脱哈尔数；γ为黏滞力系数，与流体阻力系数
Figure BDA00026451971100000914
的关系为：
Figure BDA00026451971100000915
对于海洋工程中的低质量阻尼比介质，而海洋工程中的低质量阻尼比介质具体为水，通常结构阻尼与流体阻尼比起是小量，因此本发明忽略了结构阻尼的影响，只考虑了流体阻尼。
具体实施方式三：本实施方式为具体实施方式一或二的进一步限定，本发明基于分析数据，对实例进行计算分析过程中，无量纲基本参数的具体数值如下表一：
表一 无量纲基本参数
Figure BDA0002645197110000101
本实施方式中柔性圆柱体处于静止状态时称为静止圆柱体。
本实施方式中将柔性圆柱体作为张力梁模型，分别建立张力梁模型的IL方向振动方程以及CF方向振动方程：
Figure BDA0002645197110000102
Figure BDA0002645197110000103
在公式(1)和公式(2)中，m为振动系统单位长度的质量，R为阻尼系数，T为时间，Fx(Z,T)为X方向尾流动力学引起的单位长度的外部水动力激励力，Fy(Z,T)为Y方向尾流动力学引起的单位长度的外部水动力激励力；阻尼系数R包括流体阻尼系数Rf，流体阻尼系数Rf＝γΩfρD2＝(2πStU/D)γρD2；上式中Ωf为漩涡脱落频率，St为斯脱哈尔数，γ为黏滞力系数，斯脱哈尔数St＝0.2、黏滞力系数γ与流体阻力系数
Figure BDA0002645197110000111
的关系式为：
Figure BDA0002645197110000112
结合公式(1)、公式(2)、流体阻尼系数Rf以及斯脱哈尔数St、黏滞力系数γ与流体阻力系数
Figure BDA0002645197110000113
的关系式得到Fx(Z,T)和Fy(Z,T)表达式为：
Figure BDA0002645197110000114
Figure BDA0002645197110000115
公式(3)和公式(4)中
Figure BDA0002645197110000116
为单位长度平均拖曳力，FD(Z,T)为单位长度振荡拖曳力，FL(Z,T)为单位长度升力，单位长度平均拖曳力
Figure BDA0002645197110000117
单位长度振荡拖曳力FD(Z,T)和单位长度升力FL(Z,T)分别表示为:
Figure BDA0002645197110000118
公式(5)中
Figure BDA0002645197110000119
为平均拖曳力系数且其为常数，CD(Z,T)为振荡拖曳力系数，CL(Z,T)为升力系数；
振荡拖曳力系数CD(Z,T)的表达式为CD(Z,T)＝CD0·p(Z,T)/2；
升力系数CL(Z,T)表示为CL(Z,T)＝CL0·q(Z,T)/2；
上式中CD0为静止圆柱体的振荡拖曳力系数，CL0为静止圆柱体的升力系数，p(Z,T)为与结构上的振荡拖曳力系数有关的无量纲尾流变量，q(Z,T)为与结构上的升力系数有关的无量纲尾流变量，根据表一可知，
Figure BDA00026451971100001110
CD0＝0.2，CL0＝0.3；
采用改进的Van der pol方程来满足尾流振子的非线性特性，表达式为：
Figure BDA00026451971100001111
Figure BDA00026451971100001112
公式(6)和公式(7)中εx、εy、Ax以及Ay分别为第一经验参数、第二经验参数、第三经验参数和第四经验参数，根据表一可知εx＝εy＝0.3，Ax＝Ay＝12，将公式(1)、(2)、(6)和(7)转化为无量纲形式，配合使用的表达式为：
x＝X/D,y＝Y/D,z＝Z/D,t＝T·Ωf (8)
公式(8)中t和z分别是无量纲时间和无量纲空间位置，x和y分别是IL和CF方向上的无量纲振动振幅，将公式(8)带入公式(1)、(2)、(6)和(7)中，整理得到IL方向以及CF方向结构和尾流振子的无量纲方程为：
Figure BDA0002645197110000121
Figure BDA0002645197110000122
Figure BDA0002645197110000123
公式(9)、(10)、(11)和(12)中，根据表一可知，质量比μ＝2.785，无量纲系统质量参数
Figure BDA0002645197110000124
MD和ML，根据表一可知，无量纲张力c＝23.6，无量纲弯曲刚度b＝303，无量纲张力c和无量纲弯曲刚度b的表达式分别为：
Figure BDA0002645197110000125
本实施方式中基于有限差分法对建立IL方向与CF方向振动方程进行求解的过程为在时间和空间上采用标准二阶精度中心差分格式对公式(9)～(12)进行先离散后迭代求解过程，将结构无量纲总长度L/D划分为M段；将无量纲总时间ttotal划分为N段，从而数值计算时空间步长Δz＝L/(D×M)，时间步长Δt＝ttotal/N；
本实施方式中结构无量纲总长度L/D即为结构无量纲总长度，根据表一可知，L/D的取值分别为100、200、500或1000。
如图3所示，图中表示结构细长比L/D为200时的CF方向柔性圆柱体振动位移包络线，图中标号1指向位置为振动位移包络线中第一结点位置。
如图4所示，图中表示结构细长比L/D为500时的CF方向柔性圆柱体振动位移包络线，图中标号2指向位置为振动位移包络线中第二结点位置。
如图5所示，图中表示结构细长比L/D为1000时的CF方向柔性圆柱体振动位移包络线，图中标号3指向位置为振动位移包络线中第三结点位置。
如图6所示，图中表示结构细长比L/D为100时的IL方向柔性圆柱体振动位移包络线，图中标号4指向位置为振动位移包络线中第四结点位置。
如图7所示，图中表示结构细长比L/D为200时的IL方向柔性圆柱体振动位移包络线，图中标号5指向位置为振动位移包络线中第五结点位置。
如图8所示，图中表示结构细长比L/D为500时的IL方向柔性圆柱体振动位移包络线，图中标号6指向位置为振动位移包络线中第六结点位置。
如图9所示，图中表示结构细长比L/D为1000时的IL方向柔性圆柱体振动位移包络线，图中标号7指向位置为振动位移包络线中第七结点位置。
根据表一可知，ttotal＝3000；Δt＝0.001，Δz＝1；
被划分后的M+1个空间点记为：z＝zi(i＝0,1,2,…,M)；
被划分后的N+1时间点记为：t＝tj(j＝0,1,2,..,N)；
当tn时刻zm位置处参数x、y、p和q表示为
Figure BDA0002645197110000131
以及
Figure BDA0002645197110000132
时，公式(9)～(12)中各偏导数项的二阶精度差分格式表达式分别为：
Figure BDA0002645197110000133
Figure BDA0002645197110000134
Figure BDA0002645197110000135
Figure BDA0002645197110000136
将公式(14)～(17)代入公式(9)～(12)整理得到：
Figure BDA0002645197110000137
Figure BDA0002645197110000138
Figure BDA0002645197110000141
Figure BDA0002645197110000142
x、y的初始条件设为在整个轴线上圆柱体振动位移以及速度均为0，即：
Figure BDA0002645197110000143
p和q的初始条件设为：p和q均为一微小的振幅以及
Figure BDA0002645197110000144
将
Figure BDA0002645197110000145
和
Figure BDA0002645197110000146
的表达式代入公式(14)得到：
Figure BDA0002645197110000147
将公式(22)代入公式(18)～(21)中得到t1时刻x、y、p以及q的值，表达式分别为：
Figure BDA0002645197110000148
至此，得出n＝0以及n＝1时刻的
Figure BDA0002645197110000149
以及
Figure BDA00026451971100001410
的值，由公式(18)～(21)可知，当n≥2时，无需边界条件即可直接获得2≤m≤M-2位置处的
Figure BDA00026451971100001411
的值，同时必须使用边界条件确定圆柱体两端附近四个特定位置，四个特定位置分别为m＝0、1、m-1和m处的
Figure BDA00026451971100001412
和
Figure BDA00026451971100001413
的值。
本实施方式中作为张力梁模型的柔性圆柱体的两端均采用铰接边界条件，即x和y方向的位移和弯矩始终为零，表达式为：
Figure BDA00026451971100001414
当m＝0以及m＝M时，则需结合位移边界条件对其加以求解，由两端处位移为0得到：
Figure BDA00026451971100001415
当m＝1以及m＝M-1时，则需结合弯矩为0边界条件，由两端弯矩为0得到：
Figure BDA0002645197110000151
Figure BDA0002645197110000152
将公式(27)代入公式(20)中得到m＝1以及m＝M-1时y的表达式如下：
Figure BDA0002645197110000153
Figure BDA0002645197110000154
将公式(26)代入公式(18)再结合式(28)和(29)中求得的
Figure BDA0002645197110000155
以及
Figure BDA0002645197110000156
得到m＝1以及m＝M-1时x的表达式如下：
Figure BDA0002645197110000157
Figure BDA0002645197110000158
至此得到n+1时刻整个轴线上(0≤m≤M)所有位置的振动位移
Figure BDA0002645197110000159
和
Figure BDA00026451971100001510
将计算得到的
Figure BDA00026451971100001511
和
Figure BDA00026451971100001512
代入公式(19)以及(21)即可得到n+1时刻的
Figure BDA00026451971100001513
和
Figure BDA00026451971100001514
的值，依次类推对公式(18)～(21)进行反复迭代求解便可得到整个计算时间域内的x、y、p以及q的值。
具体实施方式四：结合图2至图9结合具体数据说明本发明的张力梁模型的选用能够模拟出结构的振动位移时间历程曲线、振动频率、振动轨迹以及振动响应位移包络线的过程。
以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (5)

1.一种柔性结构横流与顺流方向涡激振动耦合响应预测方法，其特征在于：所述涡激振动耦合响应预测方法为通过分别建立IL方向与CF方向振动方程，基于有限差分法对建立IL方向与CF方向振动方程进行求解，得出计算结果结合无量纲基本参数，评价IL方向VIV振动位移响应特性以及CF方向VIV振动位移响应特性。
2.根据权利要求1所述的一种柔性结构横流与顺流方向涡激振动耦合响应预测方法，其特征在于：取一长度为L、直径为D的柔性圆柱体作为张力梁模型，在均匀来流U作用下引起的CF方向以及IL方向互为耦合的涡激振动响应，柔性圆柱体两端采用铰接边界条件，坐标系的原点O位于柔性圆柱体的底端，其中X方向为IL方向，Y方向为CF方向，Z方向为铅直方向，柔性圆柱体上的张力为Θ，柔性圆柱体上的弯曲刚度为EI，分别建立张力梁模型的IL方向振动方程以及CF方向振动方程：
Figure FDA0002645197100000011
Figure FDA0002645197100000012
在公式(1)和公式(2)中，m为振动系统单位长度的质量，R为阻尼系数，T为时间，Fx(Z,T)为X方向尾流动力学引起的单位长度的外部水动力激励力，Fy(Z,T)为Y方向尾流动力学引起的单位长度的外部水动力激励力；阻尼系数R包括流体阻尼系数Rf，流体阻尼系数Rf＝γΩfρD2＝(2πStU/D)γρD2；上式中Ωf为漩涡脱落频率，St为斯脱哈尔数，γ为黏滞力系数，斯脱哈尔数St、黏滞力系数γ与流体阻力系数
Figure FDA0002645197100000013
的关系式为：
Figure FDA0002645197100000014
结合公式(1)、公式(2)、流体阻尼系数Rf以及斯脱哈尔数St、黏滞力系数γ与流体阻力系数
Figure FDA0002645197100000015
的关系式得到Fx(Z,T)和Fy(Z,T)表达式为：
Figure FDA0002645197100000016
Figure FDA0002645197100000017
公式(3)和公式(4)中U为均匀来流流速，
Figure FDA0002645197100000018
为单位长度平均拖曳力，FD(Z,T)为单位长度振荡拖曳力，FL(Z,T)为单位长度升力，单位长度平均拖曳力
Figure FDA0002645197100000021
单位长度振荡拖曳力FD(Z,T)和单位长度升力FL(Z,T)分别表示为:
Figure FDA0002645197100000022
公式(5)中
Figure FDA0002645197100000023
为平均拖曳力系数且其为常数，CD(Z,T)为振荡拖曳力系数，CL(Z,T)为升力系数；
振荡拖曳力系数CD(Z,T)的表达式为CD(Z,T)＝CD0·p(Z,T)/2；
升力系数CL(Z,T)表示为CL(Z,T)＝CL0·q(Z,T)/2；
上式中L为柔性圆柱体的长度，D为柔性圆柱体的直径，CD0为柔性圆柱体处于静止状态下的振荡拖曳力系数，CL0为柔性圆柱体处于静止状态下的升力系数，p(Z,T)为与柔性圆柱体上的振荡拖曳力系数有关的无量纲尾流变量，q(Z,T)为与柔性圆柱体上的升力系数有关的无量纲尾流变量；
采用改进的Van der pol方程来满足尾流振子的非线性特性，表达式为：
Figure FDA0002645197100000024
Figure FDA0002645197100000025
公式(6)和公式(7)中εx、εy、Ax以及Ay分别为第一经验参数、第二经验参数、第三经验参数和第四经验参数，将公式(1)、(2)、(6)和(7)转化为无量纲形式，配合使用的表达式为：
x＝X/D,y＝Y/D,z＝Z/D,t＝T·Ωf (8)
公式(8)中t和z分别是无量纲时间和无量纲空间位置，x和y分别是IL和CF方向上的无量纲振动振幅，将公式(8)带入公式(1)、(2)、(6)和(7)中，整理得到IL方向以及CF方向结构和尾流振子的无量纲方程为：
Figure FDA0002645197100000026
Figure FDA0002645197100000027
Figure FDA0002645197100000028
Figure FDA0002645197100000031
公式(9)、(10)、(11)和(12)中，质量比μ，无量纲系统质量参数
Figure FDA0002645197100000032
MD和ML，无量纲张力c和无量纲弯曲刚度b的表达式分别为：
Figure FDA0002645197100000033
3.根据权利要求1或2所述的一种柔性结构横流与顺流方向涡激振动耦合响应预测方法，其特征在于：基于有限差分法对建立IL方向与CF方向振动方程进行求解的过程为：
在时间和空间上采用标准二阶精度中心差分格式对公式(9)～(12)进行先离散后迭代求解过程，将结构无量纲总长度L/D划分为M段；将无量纲总时间ttotal划分为N段，从而数值计算时空间步长Δz＝L/(D×M)，时间步长Δt＝ttotal/N；
被划分后的M+1个空间点记为：z＝zi(i＝0,1,2,…,M)；
被划分后的N+1时间点记为：t＝tj(j＝0,1,2,..,N)；
当tn时刻zm位置处参数x、y、p和q表示为
Figure FDA0002645197100000034
以及
Figure FDA0002645197100000035
时，公式(9)～(12)中各偏导数项的二阶精度差分格式表达式分别为：
Figure FDA0002645197100000036
Figure FDA0002645197100000037
Figure FDA0002645197100000038
Figure FDA0002645197100000039
将公式(14)～(17)代入公式(9)～(12)整理得到：
Figure FDA0002645197100000041
Figure FDA0002645197100000042
Figure FDA0002645197100000043
Figure FDA0002645197100000044
x、y的初始条件设为在整个轴线上柔性圆柱体振动位移以及速度均为0，即：
Figure FDA0002645197100000045
p和q的初始条件设为：p和q均为一微小的振幅以及
Figure FDA0002645197100000046
将
Figure FDA0002645197100000047
和
Figure FDA0002645197100000048
的表达式代入公式(14)得到：
Figure FDA0002645197100000049
将公式(22)代入公式(18)～(21)中得到t1时刻x、y、p以及q的值，表达式分别为：
Figure FDA00026451971000000410
至此，得出n＝0以及n＝1时刻的
Figure FDA00026451971000000411
以及
Figure FDA00026451971000000412
的值，由公式(18)～(21)可知，当n≥2时，无需边界条件即可直接获得2≤m≤M-2位置处的
Figure FDA00026451971000000413
的值，同时必须使用边界条件确定柔性圆柱体两端附近四个特定位置，四个特定位置分别为m＝0、1、m-1和m处的
Figure FDA00026451971000000414
和
Figure FDA00026451971000000415
的值。
4.根据权利要求3所述的一种柔性结构横流与顺流方向涡激振动耦合响应预测方法，其特征在于：作为张力梁模型的柔性圆柱体的两端均采用铰接边界条件，即x和y方向的位移和弯矩始终为零，表达式为：
Figure FDA0002645197100000051
当m＝0以及m＝M时，则需结合位移边界条件对其加以求解，由两端处位移为0得到：
Figure FDA0002645197100000052
当m＝1以及m＝M-1时，则需结合弯矩为0边界条件，由两端弯矩为0得到：
Figure FDA0002645197100000053
Figure FDA0002645197100000054
将公式(27)代入公式(20)中得到m＝1以及m＝M-1时y的表达式如下：
Figure FDA0002645197100000055
Figure FDA0002645197100000056
将公式(26)代入公式(18)再结合式(28)和(29)中求得的
Figure FDA0002645197100000057
以及
Figure FDA0002645197100000058
得到m＝1以及m＝M-1时x的表达式如下：
Figure FDA0002645197100000059
Figure FDA0002645197100000061
至此得到n+1时刻整个轴线上(0≤m≤M)所有位置的振动位移
Figure FDA0002645197100000062
和
Figure FDA0002645197100000063
将计算得到的
Figure FDA0002645197100000064
和
Figure FDA0002645197100000065
代入公式(19)以及(21)即得到n+1时刻的
Figure FDA0002645197100000066
和
Figure FDA0002645197100000067
的值，依次类推对公式(18)～(21)进行反复迭代求解便可得到整个计算时间域内的x、y、p以及q的值。
5.根据权利要求4所述的一种柔性结构横流与顺流方向涡激振动耦合响应预测方法，其特征在于：计算结果结合无量纲基本参数，评价IL方向VIV振动位移响应特性以及CF方向VIV振动位移响应特性的过程为基于无量纲基本参数数据，对计算得出
Figure FDA0002645197100000068
和
Figure FDA0002645197100000069
值进行分析计算的过程。
CN202010852497.5A 2020-08-21 2020-08-21 一种柔性结构横流与顺流方向涡激振动耦合响应预测方法 Active CN111985138B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010852497.5A CN111985138B (zh) 2020-08-21 2020-08-21 一种柔性结构横流与顺流方向涡激振动耦合响应预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010852497.5A CN111985138B (zh) 2020-08-21 2020-08-21 一种柔性结构横流与顺流方向涡激振动耦合响应预测方法

Publications (2)

Publication Number Publication Date
CN111985138A true CN111985138A (zh) 2020-11-24
CN111985138B CN111985138B (zh) 2023-12-19

Family

ID=73442961

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010852497.5A Active CN111985138B (zh) 2020-08-21 2020-08-21 一种柔性结构横流与顺流方向涡激振动耦合响应预测方法

Country Status (1)

Country Link
CN (1) CN111985138B (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113111420A (zh) * 2021-04-19 2021-07-13 哈尔滨工业大学（威海） 一种边界激励细长张力梁不稳定区间的快速预测方法
CN113572306A (zh) * 2021-08-17 2021-10-29 哈尔滨工业大学（威海） 基于变质量的宽频viv能量收集装置及其效率验证方法
CN113791541A (zh) * 2021-09-17 2021-12-14 江西理工大学 一种理论预测耦合非全同振子系统中同步包络的参数的方法
CN115391881A (zh) * 2022-08-09 2022-11-25 哈尔滨工业大学 一种桥塔尾流区吊索风致振动数值预测方法
CN115859748A (zh) * 2023-02-15 2023-03-28 山东科技大学 一种拖曳式温盐深测量仪柔性缆绳涡激振动分析方法

Citations (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20030097209A1 (en) * 2001-11-16 2003-05-22 Cedric Le Cunff System and method for limiting vortex-induced vibrations on an offshore production riser
US20090045289A1 (en) * 2007-01-10 2009-02-19 Continuum Dynamics, Inc. Flow-driven oscillating acoustic attenuator
CN101539477A (zh) * 2009-05-08 2009-09-23 中国海洋大学 一种深水顶张式立管涡激振动与疲劳分析的方法
CN102445318A (zh) * 2011-09-30 2012-05-09 中国海洋大学 一种顶张式立管顺流向振动分析方法
CN102507084A (zh) * 2011-09-30 2012-06-20 中国海洋大学 一种尾流立管的时域涡激升力确定方法
CN105184030A (zh) * 2015-11-02 2015-12-23 北京理工大学 基于弯扭复合弹簧质点模型柔性线缆位姿模拟方法及装置
CN108229043A (zh) * 2018-01-12 2018-06-29 中国海洋大学 考虑涡激效应的深海spar型浮式风机疲劳损伤分析方法
CN108414191A (zh) * 2018-02-09 2018-08-17 山东交通学院 一种考虑长径比影响的浮式圆柱涡激运动分析系统及方法
CN110046451A (zh) * 2019-04-25 2019-07-23 西南石油大学 变张力细长柔性圆柱体涡激振动响应预测方法
CN110110408A (zh) * 2019-04-25 2019-08-09 西南石油大学 刚性圆柱体横流与顺流方向涡激振动耦合响应预测方法
CN110598337A (zh) * 2019-09-17 2019-12-20 中国海洋大学 一种圆柱体涡激振动的流固耦合时域分析方法
CN110826277A (zh) * 2019-11-06 2020-02-21 中国石油大学(华东) 一种预测柔性或钢制悬链线型立管与海床土体相互作用所形成海沟长度与位置的计算方法
CN110889178A (zh) * 2019-11-26 2020-03-17 北京工业大学 一种预测谐波减速器柔轮寿命的方法
CN111177926A (zh) * 2019-12-30 2020-05-19 清华大学深圳国际研究生院 铺管过程中的涡激振动预报方法

Patent Citations (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20030097209A1 (en) * 2001-11-16 2003-05-22 Cedric Le Cunff System and method for limiting vortex-induced vibrations on an offshore production riser
US20090045289A1 (en) * 2007-01-10 2009-02-19 Continuum Dynamics, Inc. Flow-driven oscillating acoustic attenuator
CN101539477A (zh) * 2009-05-08 2009-09-23 中国海洋大学 一种深水顶张式立管涡激振动与疲劳分析的方法
CN102445318A (zh) * 2011-09-30 2012-05-09 中国海洋大学 一种顶张式立管顺流向振动分析方法
CN102507084A (zh) * 2011-09-30 2012-06-20 中国海洋大学 一种尾流立管的时域涡激升力确定方法
CN105184030A (zh) * 2015-11-02 2015-12-23 北京理工大学 基于弯扭复合弹簧质点模型柔性线缆位姿模拟方法及装置
CN108229043A (zh) * 2018-01-12 2018-06-29 中国海洋大学 考虑涡激效应的深海spar型浮式风机疲劳损伤分析方法
CN108414191A (zh) * 2018-02-09 2018-08-17 山东交通学院 一种考虑长径比影响的浮式圆柱涡激运动分析系统及方法
CN110046451A (zh) * 2019-04-25 2019-07-23 西南石油大学 变张力细长柔性圆柱体涡激振动响应预测方法
CN110110408A (zh) * 2019-04-25 2019-08-09 西南石油大学 刚性圆柱体横流与顺流方向涡激振动耦合响应预测方法
CN110598337A (zh) * 2019-09-17 2019-12-20 中国海洋大学 一种圆柱体涡激振动的流固耦合时域分析方法
CN110826277A (zh) * 2019-11-06 2020-02-21 中国石油大学(华东) 一种预测柔性或钢制悬链线型立管与海床土体相互作用所形成海沟长度与位置的计算方法
CN110889178A (zh) * 2019-11-26 2020-03-17 北京工业大学 一种预测谐波减速器柔轮寿命的方法
CN111177926A (zh) * 2019-12-30 2020-05-19 清华大学深圳国际研究生院 铺管过程中的涡激振动预报方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
冯绍军;熊友明;高云;: "刚性圆柱体涡激振动响应模态特性研究", 石油工程建设, no. 02 *
高云;邹丽;宗智;: "两端铰接的细长柔性圆柱体涡激振动响应特性数值研究", 力学学报, no. 01 *

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113111420A (zh) * 2021-04-19 2021-07-13 哈尔滨工业大学（威海） 一种边界激励细长张力梁不稳定区间的快速预测方法
CN113111420B (zh) * 2021-04-19 2022-04-22 哈尔滨工业大学（威海） 一种边界激励细长张力梁不稳定区间的快速预测方法
CN113572306A (zh) * 2021-08-17 2021-10-29 哈尔滨工业大学（威海） 基于变质量的宽频viv能量收集装置及其效率验证方法
CN113572306B (zh) * 2021-08-17 2023-07-04 哈尔滨工业大学（威海） 基于变质量的宽频viv能量收集装置及其效率验证方法
CN113791541A (zh) * 2021-09-17 2021-12-14 江西理工大学 一种理论预测耦合非全同振子系统中同步包络的参数的方法
CN113791541B (zh) * 2021-09-17 2023-08-18 江西理工大学 一种理论预测耦合非全同振子系统中同步包络的参数的方法
CN115391881A (zh) * 2022-08-09 2022-11-25 哈尔滨工业大学 一种桥塔尾流区吊索风致振动数值预测方法
CN115391881B (zh) * 2022-08-09 2023-04-18 哈尔滨工业大学 一种桥塔尾流区吊索风致振动数值预测方法
CN115859748A (zh) * 2023-02-15 2023-03-28 山东科技大学 一种拖曳式温盐深测量仪柔性缆绳涡激振动分析方法

Also Published As

Publication number Publication date
CN111985138B (zh) 2023-12-19

Similar Documents

Publication Publication Date Title
CN111985138A (zh) 一种柔性结构横流与顺流方向涡激振动耦合响应预测方法
Mittal et al. Flow-induced oscillations of two cylinders in tandem and staggered arrangements
Guilmineau et al. Numerical simulation of vortex-induced vibration of a circular cylinder with low mass-damping in a turbulent flow
Wu et al. Experimental and numerical investigation of hydroelastic response of a flexible hydrofoil in cavitating flow
Ge et al. Flow-induced vibrations of long circular cylinders modeled by coupled nonlinear oscillators
CN110110408B (zh) 刚性圆柱体横流与顺流方向涡激振动耦合响应预测方法
Munir et al. Numerical investigation of the effect of plane boundary on two-degree-of-freedom of vortex-induced vibration of a circular cylinder in oscillatory flow
Tutar et al. Large eddy simulation of a smooth circular cylinder oscillating normal to a uniform flow
Doan et al. Modeling of fluid–structure interaction for simulating vortex-induced vibration of flexible riser: finite difference method combined with wake oscillator model
Bossio et al. Numerical study on the influence of acoustic natural frequencies on the dynamic behaviour of submerged and confined disk-like structures
Li et al. Wave-induced accelerations of a fish-farm elastic floater: experimental and numerical studies
Gao et al. Time-domain prediction of the coupled cross-flow and in-line vortex-induced vibrations of a flexible cylinder using a wake oscillator model
Qu et al. Modelling of coupled cross-flow and in-line vortex-induced vibrations of flexible cylindrical structures. Part I: model description and validation
Fitt et al. The unsteady motion of two-dimensional flags with bending stiffness
Chen et al. An unsteady flow theory for vortex-induced vibration
He et al. Vortex-induced vibrations of a pipe subjected to unsynchronized support motions
Tang et al. A fully nonlinear BEM-beam coupled solver for fluid–structure interactions of flexible ships in waves
Bahmani et al. Response characteristics of a vortex-excited circular cylinder in laminar flow
Bossio V et al. Numerical modeling of the dynamical interaction between slug flow and vortex induced vibration in horizontal submarine pipelines
Bose et al. Dynamical stability analysis of a fluid structure interaction system using a high fidelity Navier-Stokes solver
CN113111420B (zh) 一种边界激励细长张力梁不稳定区间的快速预测方法
Yin et al. Wave effects on vortex-induced vibrations of a top-tensioned riser
Wang et al. Analysis of the hydrodynamic damping characteristics on a symmetrical hydrofoil
Wan et al. Overset-rans computations of two surface ships moving in viscous fluids
He et al. Numerical simulation of interactions among air, water, and rigid/flexible solid bodies

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