CN109902376A - 一种基于连续介质力学的流固耦合高精度数值模拟方法 - Google Patents
一种基于连续介质力学的流固耦合高精度数值模拟方法 Download PDFInfo
- Publication number
- CN109902376A CN109902376A CN201910136050.5A CN201910136050A CN109902376A CN 109902376 A CN109902376 A CN 109902376A CN 201910136050 A CN201910136050 A CN 201910136050A CN 109902376 A CN109902376 A CN 109902376A
- Authority
- CN
- China
- Prior art keywords
- fluid
- calculation
- solid
- interface
- physical
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 110
- 238000004088 simulation Methods 0.000 title claims abstract description 46
- 239000012530 fluid Substances 0.000 title claims abstract description 44
- 238000010168 coupling process Methods 0.000 claims abstract description 52
- 230000008878 coupling Effects 0.000 claims abstract description 42
- 238000005859 coupling reaction Methods 0.000 claims abstract description 42
- 239000000463 material Substances 0.000 claims abstract description 36
- 230000035939 shock Effects 0.000 claims abstract description 25
- 230000001681 protective effect Effects 0.000 claims abstract description 5
- 238000011089 mechanical engineering Methods 0.000 claims abstract description 4
- 238000004364 calculation method Methods 0.000 claims description 85
- 239000007787 solid Substances 0.000 claims description 82
- 230000006870 function Effects 0.000 claims description 43
- 239000000126 substance Substances 0.000 claims description 16
- 239000006185 dispersion Substances 0.000 claims description 15
- 230000008569 process Effects 0.000 claims description 15
- 230000004907 flux Effects 0.000 claims description 12
- 230000005587 bubbling Effects 0.000 claims description 4
- 239000011159 matrix material Substances 0.000 claims description 4
- 230000035515 penetration Effects 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 claims description 2
- 230000006835 compression Effects 0.000 claims description 2
- 238000007906 compression Methods 0.000 claims description 2
- 238000012937 correction Methods 0.000 claims description 2
- 230000035945 sensitivity Effects 0.000 claims description 2
- 239000011343 solid material Substances 0.000 claims description 2
- 230000002123 temporal effect Effects 0.000 claims description 2
- 230000003993 interaction Effects 0.000 abstract description 11
- 238000004422 calculation algorithm Methods 0.000 abstract description 6
- 238000005516 engineering process Methods 0.000 abstract description 4
- 230000010355 oscillation Effects 0.000 abstract description 3
- 230000009514 concussion Effects 0.000 abstract 1
- 238000004880 explosion Methods 0.000 description 9
- 230000008901 benefit Effects 0.000 description 8
- 238000011160 research Methods 0.000 description 8
- 238000002474 experimental method Methods 0.000 description 7
- 238000012360 testing method Methods 0.000 description 4
- 238000004458 analytical method Methods 0.000 description 3
- 238000011161 development Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 230000008859 change Effects 0.000 description 2
- 238000009795 derivation Methods 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 230000007246 mechanism Effects 0.000 description 2
- 239000011365 complex material Substances 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000001066 destructive effect Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 239000002360 explosive Substances 0.000 description 1
- 239000007789 gas Substances 0.000 description 1
- 239000001307 helium Substances 0.000 description 1
- 229910052734 helium Inorganic materials 0.000 description 1
- SWQJXJOGLNCZEY-UHFFFAOYSA-N helium atom Chemical compound [He] SWQJXJOGLNCZEY-UHFFFAOYSA-N 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000005312 nonlinear dynamic Methods 0.000 description 1
- 238000012827 research and development Methods 0.000 description 1
- 238000000844 transformation Methods 0.000 description 1
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明提出一种基于连续介质力学的流固耦合的高精度数值模拟方法,属于流固耦合数值模拟技术领域。本发明实现方法如下:基于本质无震荡(WENO)有限差分方法,同时耦合水平集方法(Level‑Set)和真实虚拟流体方法(RGFM)的联合算法,实现连续介质力学的流固耦合的高精度数值模拟。WENO有限差分格式能通过选取空间模板格点实现空间的高精度离散,保证计算精度并显著降低冲击波的耗散;Level‑Set方法能有效处理多物质界面的复杂拓扑结构变化;RGFM方法能有效抑制冲击波与物质界面相互作用产生的非物理振荡问题。本发明能提高对流固耦合过程预测的准确性,进而解决流固耦合领域相关工程技术问题。流固耦合领域包括高速侵彻武器、防护装备、航空航天、机械工程领域。
Description
技术领域
本发明涉及一种基于连续介质力学的流固耦合的高精度数值模拟方法,属于流固耦合数值模拟技术领域。
背景技术
爆炸发生时会伴随着巨大能量的瞬时释放,产生的爆炸冲击波会对周围结构造成灾难性破坏。认识爆炸冲击波与结构目标相互作用机制,对于预防爆炸灾害的发生具有重要意义。目前,研究爆炸与冲击问题的主要手段包括理论分析、实验和数值模拟。由于爆炸与冲击问题具有强非线性动力学行为,给理论研究带来了巨大的挑战。实验研究虽然具有客观性等优点,但存在一定的危险性,且需要大量的科研经费作为支撑。而数值模拟在一定程度上比理论和实验对问题及机制的认识更加细致和深刻,可以重复地、不间断地呈现爆炸波传播及其与结构的相互作用,揭示实验很难观测到的发生在极短时间内的内部物理现象。综上,本发明提出的方法具有重要的技术背景。
传统的多物质相互作用研究方式主要包括理论分析和实验测试。
理论分析通过建立和求解包括流体、固体运动方程、物质状态方程、边界条件等在内的封闭方程组,得出相关物理状态的解析解。这种方法在早期研究中具有一定指导意义,但同时也存在很大缺陷:(1)理论模型较为简单,研究结果无法解决复杂多变的实际工程问题。(2)推导过程中引入了很多人为假设,这与实际工况存在很大差距。(3)理论求解涉及许多复杂的数学变换,研究难度大。
实验测试是最为直接和准确的研究方法,但由于爆炸冲击过程持续在ms量级,对实验测量技术提出了很高的要求。实验测试的不足在于:(1)非接触测量实验试件内部参数在ms量级的实验技术有待于进一步研究。(2)实验具有破坏性,相关测量仪器的损耗很大,实验成本较高。(3)实验中往往存在不确定因素,可重复性较差。
随着计算机的发展,数值模拟已成为研究爆炸与冲击问题的重要工具。该技术的优势在于:(1)解决实际工程问题的效率高,适应性强,操作简单。(2) 与实验研究相比,数值模拟成本低廉,无安全风险。从上世纪60年代起,以美国三大国家武器实验室为代表,进行了大量的爆炸与冲击数值计算程序的研发工作,实现了仿真软件如AUTODYN、LS-DYNA、AUTOREGAS等在世界范围内的绝对优势地位。这些软件功能非常强大且具有很好的人机交互界面,给国内用户带来了极大的便利,但是国外进口的大部分商业软件只是开放软件的部分功能,源代码是严格保密的,使得我们很难加入新的算法模块进行功能拓展,在应用上受制于人。国内外学者提出了无网格MPM、MMALE算法等。然而这些方法许多采用低阶精度格式,影响了计算结果的准确性,加大了冲击波的耗散,在高速侵彻武器及新型防护装备的研制的高速发展大背景下,以上两点不足亟待弥补。
发明内容
针对现有技术中流固耦合数值模拟方法采用低阶精度格式,存在影响计算结果的准确性、加大冲击波耗散的问题,本发明公开的一种基于连续介质力学的流固耦合高精度数值模拟方法要解决的技术问题是:实现基于连续介质力学的流固耦合高精度数值模拟,降低冲击波耗散,提高对流固耦合过程预测的准确性,进而解决流固耦合领域相关工程技术问题。
所述流固耦合领域包括高速侵彻武器、防护装备、航空航天、机械工程领域。
本发明的目的是通过下述技术方案实现的。
本发明提出一种基于连续介质力学的流固耦合的高精度数值模拟方法,实现方法如下:基于本质无震荡(WENO)有限差分方法,同时耦合水平集方法(Level-Set)和真实虚拟流体方法(RGFM)的联合算法,实现连续介质力学的流固耦合的高精度数值模拟。WENO有限差分格式能通过选取空间模板格点实现空间的高精度离散,保证计算精度并显著降低冲击波的耗散;Level-Set方法能有效处理多物质界面的复杂拓扑结构变化; RGFM方法能有效抑制冲击波与物质界面相互作用产生的非物理振荡问题。本发明能有效对流体、固体以及冲击波相互耦合的问题进行高精度模拟。
本发明公开的一种基于连续介质力学的流固耦合高精度数值模拟方法,包括如下步骤:
步骤1:确定流体、固体计算区域,建立直角坐标系,在直角坐标系将该区域划分为m*n个网格,其中,m表示x方向网格数量,n表示y方向网格数量。
步骤2:定义计算区域内初始状态下的用于区分界面位置的Level-Set 函数、各区域物理量、建立基于连续介质力学的流体、固体控制方程组、材料状态方程函数,根据用于区分界面位置的Level-Set函数划分的各区域,并为上述函数及方程组各未知物理量赋初始值。
步骤2.1:定义计算区域内初始状态下的用于区分界面位置的Level-Set 函数,根据实际工况,通过LS方法将待计算的流体、固体计算区域划分为多区域。
定义计算区域内初始状态下的用于区分界面位置的Level-Set函数,基于Level-Set函数求解多物质界面的方法为LS方法。LS方法采用符号距离函数实现在任何时刻对界面演化过程的追踪,界面始终处于符号距离函数为零的等值面位置:
在计算初始时刻,采用公式(2)定义符号距离函数的初值:
其中,x、y、z分别为网格节点的物理空间坐标,d表示网格节点(x,y,z)到界面的距离,Ω1、Ω2分别为两种物质界面两侧的区域。
根据实际工况,通过LS方法将待计算的流体、固体计算区域划分为多区域。
步骤2.2:定义步骤2.1划分的各区域初始状态物理量。
所述初始状态物理量包括初始状态下流体的密度ρ、速度u、熵S、变形梯度张量F,初始状态固流体的密度、速度、熵、变形梯度张量。其中,变形梯度张量用于衡量物质存在初始拉伸、压缩、扭转的初始变形状态。
步骤2.3:构建多区域连续介质力学控制方程组。
由于流体、固体在自然界中均表示为连续介质,所以,由连续介质力学表示的欧拉方程组为:
其中,E为单位质量材料的总能量,ε为比内能,×为旋度算符、▽为梯度算符。
步骤2.4:选取材料的状态方程函数,进而确定材料在变形、运动下的物理状态。
作为优选,不考虑化学反应的情况下,由于材料的内能由材料的变形和熵唯一决定,步骤2.4选取材料内能表示的状态方程函数如公式(4)所示,并确定内能形式状态方程函数的相关参数。
ε=ε(g,S)=ε(F,S)=ε(C,S)=ε(G,S) (4)
所述相关参数包括逆变形梯度张量g、变形梯度张量F,左柯西-格林应变张量C,右柯西-格林应变张量G。
作为进一步优选,步骤2.4所述内能形式状态方程函数优选用左柯西-格林应变张量来统一描述弹性固体以及流体的内能形式状态方程函数,如公式(5) 所示:
其中,γ为多方气体指数,κ为熵的函数,a、p∞为材料常数,χ为固体材料剪切弹性模量,对于流体材料χ=0。
步骤3:对步骤2划分的各区域设置边界条件。
针对每一个计算格点,由于不同的离散格式,需要周围不同范围点的物理信息,故对于边界点,需要将计算区域向外拓展N层。根据用于空间导数离散采用的格式,至少需要向外拓展3层网格。
对于自由边界,有:U-i=U1,i=1,2…N,即,计算区域边界外第一层、第二层、第三层虚拟网格上的物理量均等于计算区域边界网格上的物理量,对于更高阶的精度要求,需要向外拓展更多层的虚拟网格。
对于固定边界,有:un,-i=un,i,即计算区域边界外第一层法向速度等于计算区域边界内第一层,计算区域边界外第二层法向速度等于计算区域边界内第二层。其他物理量的边界设置方法与自由边界拓展方法相同。
步骤4:选取计算参数CFL计算时间步长。
CFL参数为介于0到1之间的常数,选取计算参数CFL计算时间步长:
Δx为计算区域中x方向的网格步长;Δy为计算区域中y方向的网格步长;u、 v为计算区域中当前网格在x方向的速度分量以及y方向的速度分量;cmax为计算区域中当前网格处的最大波速,包括P波、SV波以及SH波。
步骤5:基于步骤2各物理量的初始值赋值对计算区域初始化。
在t=0时刻,为步骤2各物理量的初始值赋值,密度、速度、以及公式(5) 中参量均为已知参量,通过所述已知参量利用公式(3)、公式(5)对总能量E、比内能ε进行求解,为计算区域及虚拟网格的参量赋值实现对计算区域初始化。
步骤6:根据Level-Set函数分别求解各区域多物质界面法向量。
根据Level-Set函数,分别各区域内各物质界面的法向量及曲率进行求解,求解公式为:
其中,为空间坐标点的梯度方向向量,k为该点的空间曲率。
步骤7:求解流体计算区域界面处的黎曼问题并确定流固耦合界面附近网格的物理状态。
采用RGFM方法,沿物质界面法向建立流体、固体两种物质之间的黎曼问题并求解出两种物质界面处的物理量,所述物理状态包括密度、速度、压力、熵、逆变形梯度张量。确定固体区域内紧邻流体界面的多层网格物理状态的具体方法为:
在计算流体时,将固体计算区域设置为虚拟流体网格,所述虚拟流体上的信息通过在界面上求解黎曼问题,然后将密度、熵、逆变形梯度张量通过延拓得到。延拓方程为:
It±N·▽I=0 (9)
其中I表示任意一个物理量,N为该点梯度方向矢量。
在迭代求解过程中,各物质区域的所有节点包括已经得到更新的边界节点的状态保持不变,从而得到固体区域内紧邻流体界面的多层网格物理状态。
经延拓之后虚拟区域的物理量为局部坐标系下的物理状态,需对以上迭代得到的所有多层网格物理状态施行逆向变换转变到全局坐标系下,然后通过后续步骤9采用高精度的WENO格式进行数值求解。
步骤8:在计算固体时,将流体计算区域设置为虚拟固体网格,按照步骤7 的方法求解固体计算区域界面处的黎曼问题并确定流固耦合界面附近网格的物理状态。
步骤9:采用高精度的WENO格式对单一介质计算区域进行空间离散。
采用高精度有限差分WENO格式对步骤2由连续介质力学表示的欧拉方程组空间导数项进行离散,具体空间离散方法如下:
步骤9.1:在各个节点上计算守恒变量Ui和通量F(Ui)的差商;
步骤9.2:在半节点xi+1/2,j,k处进行局部特征空间投影,然后利用WENO格式构造半节点处的数值通量具体实现方法如下:
步骤9.2a:通过Roe平均或简单的算术平均计算半节点上均值
步骤9.2b:计算在半节点处通量的雅可比矩阵
同时给出相应的左右特征向量LF=LF(Ui+1/2,j,k),RF=RF(Ui+1/2,j,k)和对角矩阵ΛF=ΛF(Ui+1/2,j,k);
步骤9.2c:实施局部特征投影,Vl,j,k=LUl,Wl,j,k=LF(Ul),l为i节点附近的节点;
步骤9.2d:采用Lax-Friedrichs通量分裂方法对Wl,j,k进行通量分裂,利用分裂之后的值以及WENO方法构造
步骤9.2e:特征空间返回物理空间
步骤9.3:构造数值通量
采用相似的过程可以得到数值通量和通过以下空间离散格式可以得到守恒变量关于时间导数的常微分方程:
采用中心差分格式对公式(11)右端项中的空间导数修正项进行离散。
步骤10:采用TVD Runge-Kutta格式对单一介质计算区域进行时间离散。
用TVD Runge-Kutta格式对单一介质计算区域进行时间离散。优选三阶TVDRunge-Kutta格式对公式(11)的时间导数项进行离散,时间三阶TVD Runge-Kutta 格式记公式(11)右端项为L,推进时间步长为Δt,三阶精度时间推进TVD Runge-Kutta格式的具体形式为:
步骤9-步骤10为对单一介质计算区域进行空间、时间离散。
步骤11:重复步骤9-步骤10,对各物质区域分别进行空间、时间离散,得到tn+1时刻各物质区域的物理状态。
对于t>0,将步骤2、3中的物理状态进行更新,作为计算区域内各网格参量在下一步计算时的初值。
步骤12:整合步骤11得到的不同区域流场物理状态,得到全流场物理状态,对Level-Set运动方程的空间导数项进行离散并进行求解,得到tn+1时刻各物质区域Level-Set函数。
对Level-Set运动方程的空间导数项进行离散并进行求解具体实现方法如下:
Level-Set运动方程如公式(13)所示,
采用HJ-WENO(Hamilton-Jocobi Weight Essentially Non-Oscillatory)数值方法对Level-Set运动方程的空间导数项进行高精度离散,采用显式TVD Runge-Kutta格式推进时间导数项。
空间导数采用如下方式进行逼近:
优选五阶HJ-WENO方法对和进行空间离散:
其中:
IS0=13(a-b)2+3(a-3b)2
IS1=13(b-c)2+3(b+c)2
IS2=13(c-d)2+3(3c-d)2
其中,ε为常数。对采用同样的方法进行离散。
通过步骤(10)中所述的显式TVD Runge-Kutta得到Level-Set的完全离散格式,进而求得出下一时刻计算区域内各网格单元的界面函数值。
步骤13:基于当前计算的Level-Set全场值,对步骤9-步骤12得到的固体、流体流场进行选取,整合全流场物理状态,得到n+1时刻全场的物理状态。
步骤14:判断当前计算时间tn+1与设定总时间ttotal的关系:若tn+1>ttotal,则结束数值模拟,输出计算区域内网格单元的压力值;若tn+1<ttotal,则返回步骤3,按照步骤13得到的物理状态设置边界条件,继续进行流固耦合的高精度数值模拟。
还包括步骤15:利用步骤1至步骤14进行基于连续介质力学的流固耦合高精度数值模拟,降低冲击波耗散,提高对流固耦合过程预测的准确性,进而解决流固耦合领域相关工程技术问题。
所述流固耦合领域包括高速侵彻武器、防护装备、航空航天航海、机械工程领域。
预测所述工程领域中的流固耦合过程,所述流固耦合过程包括激波打气泡、弹塑性块体撞击、聚能射流、Richtmyer–Meshkov不稳定性问题。
有益效果:
与传统的流体、固体相互作用模拟方法相比,具有以下优势:
1、本发明公开的一种基于连续介质力学的流固耦合高精度数值模拟方法,采用连续介质力学控制方程组,将流体、固体统一为一种形式,实现流体固体强耦合相互作用计算,及冲击波与界面相互作用的高精度捕捉,能准确预测所述工程领域中的流固耦合过程,所述流固耦合过程包括激波打气泡、弹塑性块体撞击、聚能射流、Richtmyer–Meshkov不稳定性问题。
2、本发明公开的一种基于连续介质力学的流固耦合高精度数值模拟方法,采用WENO格式对欧拉方程空间导数项进行离散,采用TVD Runge-Kutta格式对欧拉方程时间导数项进行离散,相比传统方法具有更高的精度,能够降低冲击波及稀疏波在材料内部传播的耗散。
3、本发明公开的一种基于连续介质力学的流固耦合高精度数值模拟方法,采用GFM方法中的RGFM方法处理冲击波与多物质界面的相互作用的强间断问题,具有物质界面高分辨率的间断特性;相比于经典的GFM方法,能有效地处理流体、固体复杂的状态方程,能有效抑制求解欧拉方程时界面处产生的非物理震荡现象,保障数值模拟过程的稳定运行。
4、本发明公开的一种基于连续介质力学的流固耦合高精度数值模拟方法,采用的Level-Set方法能够精确地捕捉运动界面,与传统的Lagrange界面追踪方法相比,Level-Set方法能够更加稳定有效的处理复杂物质界面及其拓扑结构变化,尤其是处理激波打气泡、弹塑性块体撞击、聚能射流、Richtmyer–Meshkov 不稳定性非线性大变形问题时,具有明显优势。
附图说明
图1为黎曼问题示意图;
图2为局部黎曼问题的构造示意图;
图3为RGFM方法构造示意图;
图4为算例的初始几何模型图;
图5为算例t=427μs时的全流场胞格图;
图6为算例在32μs、52μs、62μs、72μs、82μs、102μs、245μs、427μs、674μs 时刻实验纹影图与计算图对比;
图7为本实施例公开一种基于连续介质力学的流固耦合高精度数值模拟方法流程图。
具体实施方式
实施例1:
为了更好的说明本发明的目的和优点,下面结合附图与实施例对本发明作进一步说明。
本实施例在空气激波打氦气泡算例中的应用。
本实施例公开的算例中的初始几何模型在附图4中呈现,不同时刻,计算得到的He气泡模拟效果在附图6中呈现。
如图7所示,本实施例公开一种基于连续介质力学的流固耦合高精度数值模拟方法,具体实现方法为:
首先给出步骤1、2中所需的信息:如附图4所示,计算区域尺寸为 445mm*89mm,横向网格数2225,纵向网格数445,网格宽度0.5mm,He气泡中心点坐标为(225mm,44.5mm),半径r=25mm。设定总时间为983μs,为了算法稳定性,步骤4中CFL数取0.1。步骤1、2中He气泡及空气状态方程参数取如下所示,均为国际单位制:
表格1:
由此,完成步骤1、2。
进行步骤3,将上、下边界设置为反射边界条件,将左、右边界条件分别设置为出流和入流边界条件,保正物理模型的正确性。
然后进行步骤4对于计算区域进行计算时间步长的计算,通过不同的状态方程分别计算He气泡区域和空气所在区域最小的计算时间步长,然后取全流场计算时间步长的最小值,保正计算的稳定性。
接着进行步骤5至步骤14,得到空气中激波打He气泡的多物质及激波相互作用的数值模拟结果。现有数值模拟技术通过将空气和He气泡用统一的方程描述,方便得到界面两侧的物理状态,进而通过RGFM方法进行高精度数值仿真。
结果分析:
附图6为不同时刻计算得到的结果图,从图中可知,当冲击波从空气进入He气泡时,会产生入射激波,入射激波在He气泡内扫过He气泡边界,会产生波纹状的稀疏波,波纹状的稀疏波会汇聚到He气泡初始形状的1/3直径处,稀疏波的汇聚会在He气泡内形成低压、低密度区,导致右端空气-He界面射流由于不稳定产生发散现象。通过计算得到的图像与实验得到的图做对比,证明该本实施例公开一种基于连续介质力学的流固耦合高精度数值模拟方法的有效性。
以上所述的具体描述,对发明的目的、技术方案和有益效果进行了进一步详细说明,所应理解的是,以上所述仅为本发明的具体实施例而已,并不用于限定本发明的保护范围,凡在本发明的精神和原则之内,所做的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (10)
1.一种基于连续介质力学的流固耦合高精度数值模拟方法,其特征在于:包括如下步骤,
步骤1:确定流体、固体计算区域,建立直角坐标系,在直角坐标系将该区域划分为m*n个网格,其中,m表示x方向网格数量,n表示y方向网格数量;
步骤2:定义计算区域内初始状态下的用于区分界面位置的Level-Set函数、各区域物理量、建立基于连续介质力学的流体、固体控制方程组、材料状态方程函数,根据用于区分界面位置的Level-Set函数划分的各区域,并为上述函数及方程组各未知物理量赋初始值;
步骤3:对步骤2划分的各区域设置边界条件;
步骤4:选取计算参数CFL计算时间步长;
步骤5:基于步骤2各物理量的初始值赋值对计算区域初始化;
步骤6:根据Level-Set函数分别求解各区域多物质界面法向量;
步骤7:求解流体计算区域界面处的黎曼问题并确定流固耦合界面附近网格的物理状态;
步骤8:在计算固体时,将流体计算区域设置为虚拟固体网格,按照步骤7的方法求解固体计算区域界面处的黎曼问题并确定流固耦合界面附近网格的物理状态;
步骤9:采用高精度的WENO格式对单一介质计算区域进行空间离散;
步骤10:采用TVD Runge-Kutta格式对单一介质计算区域进行时间离散;
步骤11:重复步骤9-步骤10,对各物质区域分别进行空间、时间离散,得到tn+1时刻各物质区域的物理状态;
步骤12:整合步骤11得到的不同区域流场物理状态,得到全流场物理状态,对Level-Set运动方程的空间导数项进行离散并进行求解,得到tn+1时刻各物质区域Level-Set函数;
步骤13:基于当前计算的Level-Set全场值,对步骤9-步骤12得到的固体、流体流场进行选取,整合全流场物理状态,得到n+1时刻全场的物理状态;
步骤14:判断当前计算时间tn+1与设定总时间ttotal的关系:若tn+1>ttotal,则结束数值模拟,输出计算区域内网格单元的压力值;若tn+1<ttotal,则返回步骤3,按照步骤13得到的物理状态设置边界条件,继续进行流固耦合的高精度数值模拟。
2.如权利要求1所述的一种基于连续介质力学的流固耦合高精度数值模拟方法,其特征在于:还包括步骤15,利用步骤1至步骤14进行基于连续介质力学的流固耦合高精度数值模拟,降低冲击波耗散,提高对流固耦合过程预测的准确性,进而解决流固耦合领域相关工程技术问题。
3.如权利要求1或2所述的一种基于连续介质力学的流固耦合高精度数值模拟方法,其特征在于:步骤2实现方法为,
步骤2.1:定义计算区域内初始状态下的用于区分界面位置的level-Set函数,根据实际工况,通过LS方法将待计算的流体、固体计算区域划分为多区域;
定义计算区域内初始状态下的用于区分界面位置的Level-Set函数,基于Level-Set函数求解多物质界面的方法为LS方法;LS方法采用符号距离函数实现在任何时刻对界面演化过程的追踪,界面始终处于符号距离函数为零的等值面位置:
在计算初始时刻,采用公式(2)定义符号距离函数的初值:
其中,x、y、z分别为网格节点的物理空间坐标,d表示网格节点(x,y,z)到界面的距离,Ω1、Ω2分别为两种物质界面两侧的区域;
根据实际工况,通过LS方法将待计算的流体、固体计算区域划分为多区域;
步骤2.2:定义步骤2.1划分的各区域初始状态物理量;
所述初始状态物理量包括初始状态下流体的密度ρ、速度u、熵S、变形梯度张量F,初始状态固流体的密度、速度、熵、变形梯度张量;其中,变形梯度张量用于衡量物质存在初始拉伸、压缩、扭转的初始变形状态;
步骤2.3:构建多区域连续介质力学控制方程组;
由于流体、固体在自然界中均表示为连续介质,所以,由连续介质力学表示的欧拉方程组为:
其中,E为单位质量材料的总能量,ε为比内能,×为旋度算符、为梯度算符;
步骤2.4:选取材料的状态方程函数,进而确定材料在变形、运动下的物理状态;
不考虑化学反应的情况下,由于材料的内能由材料的变形和熵唯一决定,步骤2.4选取材料内能表示的状态方程函数如公式(4)所示,并确定内能形式状态方程函数的相关参数;
ε=ε(g,S)=ε(F,S)=ε(C,S)=ε(G,S) (4)
所述相关参数包括逆变形梯度张量g、变形梯度张量F,左柯西-格林应变张量C,右柯西-格林应变张量G;
步骤2.4所述内能形式状态方程函数优选用左柯西-格林应变张量来统一描述弹性固体以及流体的内能形式状态方程函数,如公式(5)所示:
其中,γ为多方气体指数,κ为熵的函数,a、p∞为材料常数,χ为固体材料剪切弹性模量,对于流体材料χ=0。
4.如权利要求3所述的一种基于连续介质力学的流固耦合高精度数值模拟方法,其特征在于:步骤3实现方法为,
针对每一个计算格点,由于不同的离散格式,需要周围不同范围点的物理信息,故对于边界点,需要将计算区域向外拓展N层;根据用于空间导数离散采用的格式,至少需要向外拓展3层网格;
对于自由边界,有:U-i=U1,i=1,2…N,即,计算区域边界外第一层、第二层、第三层虚拟网格上的物理量均等于计算区域边界网格上的物理量,对于更高阶的精度要求,需要向外拓展更多层的虚拟网格;
对于固定边界,有:un,-i=un,i,即计算区域边界外第一层法向速度等于计算区域边界内第一层,计算区域边界外第二层法向速度等于计算区域边界内第二层;其他物理量的边界设置方法与自由边界拓展方法相同。
5.如权利要求4所述的一种基于连续介质力学的流固耦合高精度数值模拟方法,其特征在于:步骤4实现方法为,
CFL参数为介于0到1之间的常数,选取计算参数CFL计算时间步长:
△x为计算区域中x方向的网格步长;△y为计算区域中y方向的网格步长;u、v为计算区域中当前网格在x方向的速度分量以及y方向的速度分量;cmax为计算区域中当前网格处的最大波速,包括P波、SV波以及SH波;
步骤5实现方法为,
在t=0时刻,为步骤2各物理量的初始值赋值,密度、速度、以及公式(5)中参量均为已知参量,通过所述已知参量利用公式(3)、公式(5)对总能量E、比内能ε进行求解,为计算区域及虚拟网格的参量赋值实现对计算区域初始化。
6.如权利要求5所述的一种基于连续介质力学的流固耦合高精度数值模拟方法,其特征在于:步骤6实现方法为,
根据Level-Set函数,分别各区域内各物质界面的法向量及曲率进行求解,求解公式为:
其中,为空间坐标点的梯度方向向量,k为该点的空间曲率;
步骤7实现方法为,
采用RGFM方法,沿物质界面法向建立流体、固体两种物质之间的黎曼问题并求解出两种物质界面处的物理量,所述物理状态包括密度、速度、压力、熵、逆变形梯度张量;确定固体区域内紧邻流体界面的多层网格物理状态的具体方法为:
在计算流体时,将固体计算区域设置为虚拟流体网格,所述虚拟流体上的信息通过在界面上求解黎曼问题,然后将密度、熵、逆变形梯度张量通过延拓得到;延拓方程为:
其中,I表示任意一个物理量,N为该点梯度方向矢量;
在迭代求解过程中,各物质区域的所有节点包括已经得到更新的边界节点的状态保持不变,从而得到固体区域内紧邻流体界面的多层网格物理状态;
经延拓之后虚拟区域的物理量为局部坐标系下的物理状态,需对以上迭代得到的所有多层网格物理状态施行逆向变换转变到全局坐标系下,然后通过后续步骤9采用高精度的WENO格式进行数值求解。
7.如权利要求6所述的一种基于连续介质力学的流固耦合高精度数值模拟方法,其特征在于:步骤9采用高精度有限差分WENO格式对步骤2由连续介质力学表示的欧拉方程组空间导数项进行离散,具体空间离散方法如下,
步骤9.1:在各个节点上计算守恒变量Ui和通量F(Ui)的差商;
步骤9.2:在半节点xi+1/2,j,k处进行局部特征空间投影,然后利用WENO格式构造半节点处的数值通量具体实现方法如下:
步骤9.2a:通过Roe平均或简单的算术平均计算半节点上均值
步骤9.2b:计算在半节点处通量的雅可比矩阵
同时给出相应的左右特征向量LF=LF(Ui+1/2,j,k),RF=RF(Ui+1/2,j,k)和对角矩阵ΛF=ΛF(Ui+1/2,j,k);
步骤9.2c:实施局部特征投影,Vl,j,k=LUl,Wl,j,k=LF(Ul),l为i节点附近的节点;
步骤9.2d:采用Lax-Friedrichs通量分裂方法对Wl,j,k进行通量分裂,利用分裂之后的值以及WENO方法构造
步骤9.2e:特征空间返回物理空间
步骤9.3:构造数值通量
采用相似的过程可以得到数值通量和通过以下空间离散格式可以得到守恒变量关于时间导数的常微分方程:
采用中心差分格式对公式(11)右端项中的空间导数修正项进行离散;
步骤10实现方法为,
用TVD Runge-Kutta格式对单一介质计算区域进行时间离散;优选三阶TVD Runge-Kutta格式对公式(11)的时间导数项进行离散,时间三阶TVD Runge-Kutta格式记公式(11)右端项为L,推进时间步长为△t,三阶精度时间推进TVD Runge-Kutta格式的具体形式为:
步骤9-步骤10为对单一介质计算区域进行空间、时间离散。
8.如权利要求7所述的一种基于连续介质力学的流固耦合高精度数值模拟方法,其特征在于:步骤11实现方法为,
对于t>0,将步骤2、3中的物理状态进行更新,作为计算区域内各网格参量在下一步计算时的初值。
9.如权利要求8所述的一种基于连续介质力学的流固耦合高精度数值模拟方法,其特征在于:步骤12实现方法为,
对Level-Set运动方程的空间导数项进行离散并进行求解具体实现方法如下:
Level-Set运动方程如公式(13)所示,
采用HJ-WENO(Hamilton-Jocobi Weight Essentially Non-Oscillatory)数值方法对Level-Set运动方程的空间导数项进行高精度离散,采用显式TVD Runge-Kutta格式推进时间导数项;
空间导数采用如下方式进行逼近:
优选五阶HJ-WENO方法对和进行空间离散:
其中:
IS0=13(a-b)2+3(a-3b)2
IS1=13(b-c)2+3(b+c)2
IS2=13(c-d)2+3(3c-d)2
其中,ε为常数;对采用同样的方法进行离散;
通过步骤(10)中所述的显式TVD Runge-Kutta得到Level-Set的完全离散格式,进而求得出下一时刻计算区域内各网格单元的界面函数值;
10.如权利要求9所述的一种基于连续介质力学的流固耦合高精度数值模拟方法,其特征在于:所述流固耦合领域包括高速侵彻武器、防护装备、航空航天航海、机械工程领域;
预测所述工程领域中的流固耦合过程,所述流固耦合过程包括激波打气泡、弹塑性块体撞击、聚能射流、Richtmyer-Meshkov不稳定性问题。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910136050.5A CN109902376B (zh) | 2019-02-25 | 2019-02-25 | 一种基于连续介质力学的流固耦合高精度数值模拟方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910136050.5A CN109902376B (zh) | 2019-02-25 | 2019-02-25 | 一种基于连续介质力学的流固耦合高精度数值模拟方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN109902376A true CN109902376A (zh) | 2019-06-18 |
CN109902376B CN109902376B (zh) | 2021-01-15 |
Family
ID=66945291
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910136050.5A Active CN109902376B (zh) | 2019-02-25 | 2019-02-25 | 一种基于连续介质力学的流固耦合高精度数值模拟方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN109902376B (zh) |
Cited By (17)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110598323A (zh) * | 2019-09-12 | 2019-12-20 | 中国水利水电科学研究院 | 一种渗透破坏离散元模拟方法 |
CN110852005A (zh) * | 2019-10-21 | 2020-02-28 | 北京理工大学 | 一种大规模并行计算的自适应扩大计算域的数值模拟方法 |
CN111159958A (zh) * | 2019-12-10 | 2020-05-15 | 北京航空航天大学 | 一种多介质耦合问题界面两侧状态的保物理特性获取方法 |
CN111339715A (zh) * | 2020-02-20 | 2020-06-26 | 中国石油大学(北京) | 液下可形变物质与固体表面的粘附力的确定方法及装置 |
CN111783277A (zh) * | 2020-06-04 | 2020-10-16 | 海仿(上海)科技有限公司 | 流体固体界面解耦算法、装置及设备 |
CN111783365A (zh) * | 2020-06-04 | 2020-10-16 | 三多(杭州)科技有限公司 | 应用于低压界面处理的虚拟介质方法、装置和设备 |
CN112100835A (zh) * | 2020-09-06 | 2020-12-18 | 西北工业大学 | 一种适用于复杂流动的高效高精度数值模拟方法 |
CN112307669A (zh) * | 2020-11-25 | 2021-02-02 | 北京理工大学 | 一种基于实际物理点迭代修正的多介质水平集方法 |
CN112861445A (zh) * | 2020-12-29 | 2021-05-28 | 中国船舶重工集团公司第七研究院 | 一种无网格数值模型的模拟方法 |
CN113378440A (zh) * | 2021-06-23 | 2021-09-10 | 四川大学 | 一种流固耦合数值模拟计算方法、装置及设备 |
CN113408168A (zh) * | 2021-06-18 | 2021-09-17 | 北京理工大学 | 一种基于黎曼问题精确求解的高精度数值模拟方法 |
CN113673185A (zh) * | 2021-08-25 | 2021-11-19 | 北京理工大学 | 一种加权双向映射的精确捕捉激波间断面方法 |
CN114169184A (zh) * | 2021-10-19 | 2022-03-11 | 北京理工大学 | 一种高精度自适应有限体积-有限差分耦合数值模拟方法 |
CN114444348A (zh) * | 2021-12-31 | 2022-05-06 | 海油发展珠海管道工程有限公司 | 一种螺旋列板涡激振动抑制装置的动力学设计方法 |
CN115329646A (zh) * | 2022-10-12 | 2022-11-11 | 国家超级计算天津中心 | 冲击波仿真方法、装置、设备及介质 |
CN118228640A (zh) * | 2024-05-23 | 2024-06-21 | 北京大学 | 磁流体数值模拟方法、装置及非易失性存储介质 |
CN118656990A (zh) * | 2024-08-19 | 2024-09-17 | 国家超级计算天津中心 | 一种流固耦合分析方法、装置、设备及存储介质 |
Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US7657415B2 (en) * | 2002-05-31 | 2010-02-02 | Schlumberger Technology Corporation | Subterranean formation treatment methods using a darcy scale and pore scale model |
CN102609582A (zh) * | 2012-02-06 | 2012-07-25 | 中国核电工程有限公司 | 一种评估砖墙在爆炸载荷作用下局部破坏结果的方法 |
CN103544342A (zh) * | 2013-09-30 | 2014-01-29 | 上海交通大学 | 基于混合模型的核电站防波堤越浪冲击模拟方法 |
CN104179514A (zh) * | 2014-08-18 | 2014-12-03 | 同济大学 | 水下隧道破碎围岩突水预测与渗流控制的方法 |
CN106354918A (zh) * | 2016-08-26 | 2017-01-25 | 中国科学院力学研究所 | 一种水力压裂中流固耦合问题数值模拟的构建方法 |
CN108491619A (zh) * | 2018-03-19 | 2018-09-04 | 浙江大学 | 基于物理与非物理混合的复杂场景流固耦合高效模拟方法 |
CN109214082A (zh) * | 2018-09-03 | 2019-01-15 | 北京理工大学 | 一种近场水下爆炸冲击波载荷的高精度数值模拟方法 |
-
2019
- 2019-02-25 CN CN201910136050.5A patent/CN109902376B/zh active Active
Patent Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US7657415B2 (en) * | 2002-05-31 | 2010-02-02 | Schlumberger Technology Corporation | Subterranean formation treatment methods using a darcy scale and pore scale model |
CN102609582A (zh) * | 2012-02-06 | 2012-07-25 | 中国核电工程有限公司 | 一种评估砖墙在爆炸载荷作用下局部破坏结果的方法 |
CN103544342A (zh) * | 2013-09-30 | 2014-01-29 | 上海交通大学 | 基于混合模型的核电站防波堤越浪冲击模拟方法 |
CN104179514A (zh) * | 2014-08-18 | 2014-12-03 | 同济大学 | 水下隧道破碎围岩突水预测与渗流控制的方法 |
CN106354918A (zh) * | 2016-08-26 | 2017-01-25 | 中国科学院力学研究所 | 一种水力压裂中流固耦合问题数值模拟的构建方法 |
CN108491619A (zh) * | 2018-03-19 | 2018-09-04 | 浙江大学 | 基于物理与非物理混合的复杂场景流固耦合高效模拟方法 |
CN109214082A (zh) * | 2018-09-03 | 2019-01-15 | 北京理工大学 | 一种近场水下爆炸冲击波载荷的高精度数值模拟方法 |
Non-Patent Citations (2)
Title |
---|
WU KAITENG等: "Level Set interface treatment and its application in Euler method", 《SCIENCE CHINA》 * |
郝莉 等: "水中障碍物爆炸毁伤效应的二维数值模拟", 《高压物理学报》 * |
Cited By (26)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110598323A (zh) * | 2019-09-12 | 2019-12-20 | 中国水利水电科学研究院 | 一种渗透破坏离散元模拟方法 |
CN110852005A (zh) * | 2019-10-21 | 2020-02-28 | 北京理工大学 | 一种大规模并行计算的自适应扩大计算域的数值模拟方法 |
CN110852005B (zh) * | 2019-10-21 | 2021-06-15 | 北京理工大学 | 一种大规模并行计算的自适应扩大计算域的数值模拟方法 |
CN111159958A (zh) * | 2019-12-10 | 2020-05-15 | 北京航空航天大学 | 一种多介质耦合问题界面两侧状态的保物理特性获取方法 |
CN111339715B (zh) * | 2020-02-20 | 2022-07-12 | 中国石油大学(北京) | 液下可形变物质与固体表面的粘附力的确定方法及装置 |
CN111339715A (zh) * | 2020-02-20 | 2020-06-26 | 中国石油大学(北京) | 液下可形变物质与固体表面的粘附力的确定方法及装置 |
CN111783277A (zh) * | 2020-06-04 | 2020-10-16 | 海仿(上海)科技有限公司 | 流体固体界面解耦算法、装置及设备 |
CN111783365A (zh) * | 2020-06-04 | 2020-10-16 | 三多(杭州)科技有限公司 | 应用于低压界面处理的虚拟介质方法、装置和设备 |
CN111783277B (zh) * | 2020-06-04 | 2024-04-12 | 海仿(上海)科技有限公司 | 流体固体界面解耦算法、装置及设备 |
CN111783365B (zh) * | 2020-06-04 | 2023-09-22 | 三多(杭州)科技有限公司 | 应用于低压界面处理的虚拟介质方法、装置和设备 |
CN112100835A (zh) * | 2020-09-06 | 2020-12-18 | 西北工业大学 | 一种适用于复杂流动的高效高精度数值模拟方法 |
CN112307669B (zh) * | 2020-11-25 | 2023-04-07 | 北京理工大学 | 一种基于实际物理点迭代修正的多介质水平集方法 |
CN112307669A (zh) * | 2020-11-25 | 2021-02-02 | 北京理工大学 | 一种基于实际物理点迭代修正的多介质水平集方法 |
CN112861445A (zh) * | 2020-12-29 | 2021-05-28 | 中国船舶重工集团公司第七研究院 | 一种无网格数值模型的模拟方法 |
CN113408168A (zh) * | 2021-06-18 | 2021-09-17 | 北京理工大学 | 一种基于黎曼问题精确求解的高精度数值模拟方法 |
CN113378440A (zh) * | 2021-06-23 | 2021-09-10 | 四川大学 | 一种流固耦合数值模拟计算方法、装置及设备 |
CN113673185A (zh) * | 2021-08-25 | 2021-11-19 | 北京理工大学 | 一种加权双向映射的精确捕捉激波间断面方法 |
CN113673185B (zh) * | 2021-08-25 | 2024-01-30 | 北京理工大学 | 一种加权双向映射的精确捕捉激波间断面方法 |
CN114169184A (zh) * | 2021-10-19 | 2022-03-11 | 北京理工大学 | 一种高精度自适应有限体积-有限差分耦合数值模拟方法 |
CN114444348A (zh) * | 2021-12-31 | 2022-05-06 | 海油发展珠海管道工程有限公司 | 一种螺旋列板涡激振动抑制装置的动力学设计方法 |
CN114444348B (zh) * | 2021-12-31 | 2024-04-16 | 海油发展珠海管道工程有限公司 | 一种螺旋列板涡激振动抑制装置的动力学设计方法 |
CN115329646A (zh) * | 2022-10-12 | 2022-11-11 | 国家超级计算天津中心 | 冲击波仿真方法、装置、设备及介质 |
CN115329646B (zh) * | 2022-10-12 | 2023-03-10 | 国家超级计算天津中心 | 冲击波仿真方法、装置、设备及介质 |
CN118228640A (zh) * | 2024-05-23 | 2024-06-21 | 北京大学 | 磁流体数值模拟方法、装置及非易失性存储介质 |
CN118228640B (zh) * | 2024-05-23 | 2024-09-10 | 北京大学 | 磁流体数值模拟方法、装置及非易失性存储介质 |
CN118656990A (zh) * | 2024-08-19 | 2024-09-17 | 国家超级计算天津中心 | 一种流固耦合分析方法、装置、设备及存储介质 |
Also Published As
Publication number | Publication date |
---|---|
CN109902376B (zh) | 2021-01-15 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN109902376B (zh) | 一种基于连续介质力学的流固耦合高精度数值模拟方法 | |
Filipovic et al. | An implicit algorithm within the arbitrary Lagrangian–Eulerian formulation for solving incompressible fluid flow with large boundary motions | |
CN108763683B (zh) | 一种三角函数框架下新weno格式构造方法 | |
Marcantoni et al. | High speed flow simulation using openfoam | |
Gupta et al. | Analysis and improvements of global–local enrichments for the generalized finite element method | |
Su et al. | Fast convergence and asymptotic preserving of the general synthetic iterative scheme | |
CN107220399A (zh) | 基于埃尔米特插值基本加权无振荡格式的全流场模拟方法 | |
CN109726465B (zh) | 基于非结构曲边网格的三维无粘低速绕流的数值模拟方法 | |
Lin et al. | Simulation of compressible two-phase flows with topology change of fluid–fluid interface by a robust cut-cell method | |
Pederson et al. | The Sedov blast wave as a radial piston verification test | |
Qiu et al. | An ellipsoidal Newton’s iteration method of nonlinear structural systems with uncertain-but-bounded parameters | |
Woodgate et al. | Fast prediction of transonic aeroelastic stability and limit cycles | |
He et al. | Characteristic-based and interface-sharpening algorithm for high-order simulations of immiscible compressible multi-material flows | |
CN111177965B (zh) | 基于定常问题求解的多重分辨weno格式定点快速扫描方法 | |
Hejranfar et al. | Chebyshev collocation spectral lattice Boltzmann method in generalized curvilinear coordinates | |
He et al. | New speedup algorithms for nonlinear dynamic time history analysis of supertall building structures under strong earthquakes | |
Vanharen et al. | Mesh adaptation for fluid-structure interaction problems | |
Korneev et al. | DiamondTorre GPU implementation algorithm of the RKDG solver for fluid dynamics and its using for the numerical simulation of the bubble-shock interaction problem | |
Farinatti Aymone | Mesh motion techniques for the ALE formulation in 3D large deformation problems | |
Ceze et al. | Development of a high-order space-time matrix-free adjoint solver | |
Valente et al. | Parameter identification and shape optimization: An integrated methodology in metal forming and structural applications | |
Martin et al. | An improved unsplit and convolutional perfectly matched layer absorbing technique for the navier-stokes equations using cut-off frequency shift | |
Hejranfar et al. | Arbitrary Lagrangian-Eulerian unstructured finite-volume lattice-Boltzmann method for computing two-dimensional compressible inviscid flows over moving bodies | |
Forti et al. | A parallel algorithm for the solution of large-scale nonconforming fluid-structure interaction problems in hemodynamics | |
Bubnovich et al. | Computation of transient natural convection in a square cavity by an implicit finite-difference scheme in terms of the stream function and temperature |
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 |