CN111259325A - 一种基于局部曲率自适应修正的改进水平集方法 - Google Patents
一种基于局部曲率自适应修正的改进水平集方法 Download PDFInfo
- Publication number
- CN111259325A CN111259325A CN202010091056.8A CN202010091056A CN111259325A CN 111259325 A CN111259325 A CN 111259325A CN 202010091056 A CN202010091056 A CN 202010091056A CN 111259325 A CN111259325 A CN 111259325A
- Authority
- CN
- China
- Prior art keywords
- level set
- particles
- interface
- set function
- grid
- 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.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 82
- 238000012937 correction Methods 0.000 title claims abstract description 58
- 230000003044 adaptive effect Effects 0.000 title claims abstract description 35
- 239000002245 particle Substances 0.000 claims abstract description 184
- 239000000700 radioactive tracer Substances 0.000 claims abstract description 63
- 238000004364 calculation method Methods 0.000 claims abstract description 58
- 238000004088 simulation Methods 0.000 claims abstract description 37
- 230000003993 interaction Effects 0.000 claims abstract description 27
- 239000006185 dispersion Substances 0.000 claims abstract description 15
- 230000006870 function Effects 0.000 claims description 135
- 230000008569 process Effects 0.000 claims description 14
- 238000004422 calculation algorithm Methods 0.000 claims description 12
- 230000001681 protective effect Effects 0.000 claims description 11
- 239000000126 substance Substances 0.000 claims description 9
- 230000035515 penetration Effects 0.000 claims description 8
- 238000011089 mechanical engineering Methods 0.000 claims description 7
- 239000013598 vector Substances 0.000 claims description 7
- 238000004519 manufacturing process Methods 0.000 claims description 6
- 230000035939 shock Effects 0.000 claims description 6
- 230000005587 bubbling Effects 0.000 claims description 3
- 238000013461 design Methods 0.000 claims description 3
- 238000009826 distribution Methods 0.000 claims description 3
- 230000010354 integration Effects 0.000 claims description 3
- 238000011056 performance test Methods 0.000 claims description 3
- 230000009467 reduction Effects 0.000 claims description 2
- 101150049349 setA gene Proteins 0.000 claims description 2
- 230000008859 change Effects 0.000 abstract description 4
- 238000004880 explosion Methods 0.000 description 8
- 230000008901 benefit Effects 0.000 description 4
- 238000012545 processing Methods 0.000 description 4
- 238000011161 development Methods 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 238000007796 conventional method Methods 0.000 description 2
- 238000010168 coupling process Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 239000002360 explosive Substances 0.000 description 2
- 239000012530 fluid Substances 0.000 description 2
- 239000007788 liquid Substances 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 238000003672 processing method Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 230000001052 transient effect Effects 0.000 description 2
- 230000008878 coupling Effects 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000007246 mechanism 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
- 238000004445 quantitative analysis Methods 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 230000007480 spreading Effects 0.000 description 1
- 238000003892 spreading Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/15—Correlation function computation including computation of convolution operations
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Computing Systems (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开的一种基于局部曲率自适应修正的改进水平集方法，属于多介质相互作用数值模拟技术领域。本发明实现方法如下：基于水平集Level Set方法，通过WENO有限差分格式和TVD Runge‑Kutta格式对水平集函数进行离散；基于水平集函数的局部曲率自适应播撒示踪粒子，通过示踪粒子对水平集函数进行修正；WENO有限差分格式和TVD Runge‑Kutta格式能实现水平集函数空间与时间项的高精度离散，保证界面水平集函数的计算精度；示踪粒子能在界面曲率变化较大的位置对水平集函数进行修正，确保多介质界面相互作用的数值仿真预测精度与效率。本发明有效的对复杂的工程问题多介质界面的运动进行高精度的数值模拟，解决多介质相互作用技术领域相应工程技术问题。
Description
技术领域
本发明涉及一种基于局部曲率自适应修正的改进水平集方法，属于多介质相互作用数值模拟技术领域。
背景技术
多介质相互作用广泛的存在于自然界，并在工业、军事中有着重要的应用。认识多介质相互作用机理，对于爆炸冲击、武器弹药领域的发展具有重要意义。目前，研究多介质相互作用问题的主要手段包括理论分析、实验和数值模拟。由于爆炸冲击瞬态作用下，多介质相互作用的强非线性动力学行为，界面拓扑性质会发生复杂的改变，理论研究仅限于流致震动等界面几何性质简单的流固耦合问题。实验研究难以捕捉材料瞬态变化的历史过程，目前仍然缺乏针对一些特定问题的测试手段，且需要耗费大量科研经费。而数值模拟可以细致、深刻的真实再现爆炸冲击下炸药、壳体、靶板结构多介质相互作用的全物理过程，从多维度定点、定量进行分析，揭示实验难以观测到的物理现象。多介质界面处理是爆炸冲击高精度数值模拟的关键组成部分，综上，本发明提出的方法具有重要的工程应用价值和技术背景。
由于爆炸冲击下材料会产生拉伸、扭曲等大变形，只能通过欧拉算法进行高精度计算。爆炸冲击多介质相互作用数值模拟欧拉算法中的界面处理分为“界面追踪”和“界面捕捉”两种类型。“界面追踪”基本为拉格朗日方法，在数值模拟中，通过显式的推进界面运动来区分不同种物质。界面追踪方法表示的界面位置清晰、直观，但是当界面附近发生较大变形甚至是拓扑性质改变时，界面追踪方法会失效，产生较大误差，最终导致计算的不稳定。“界面捕捉”基本为欧拉方法，在计算中隐式地求解网格点上的体积分数、符号距离函数等信息，再通过求解得到的物理信息重构得到界面位置。由于界面信息是通过“重构”得到的，界面捕捉方法能够自然的处理多介质相互作用过程中界面拓扑性质发生显著改变的现象。
随着数值模拟在爆炸冲击、武器弹药领域中的广泛应用，复杂工况对多介质相互作用界面处理方法提出了更高的要求，单一的“界面追踪”或“界面捕捉”已经不能满足实际问题的高效率、高分辨率的界面处理。国内外学者提出并改进了水平集(Level Set)、流体体积分数(Volume of Fluid)等算法。然而这些方法许多在提高精度的同时，耗费了大量的计算量，计算效率低下，或者在牺牲精度的情况下提高计算效率，在高速侵彻武器及新型防护装备的研制及航空、机械领域的高速发展大背景下，以上不足影响了多介质相互作用数值模拟在实际工况下的有效应用，亟待弥补。
发明内容
针对现有技术中多介质相互作用界面处理方法精度低、计算效率低，影响实际工况下界面处理不稳定、不精细的问题，本发明公开的一种基于局部曲率自适应修正的改进水平集方法要解决的技术问题是：实现界面追踪的高精度数值模拟，大幅提高计算效率，提高多介质耦合过程预测的准确性，进而进行多介质相互作用技术领域工程技术问题高精度数值模拟与预测，解决多介质相互作用技术领域相应工程技术问题。
所述多介质相互作用技术领域包括爆炸冲击、高速侵彻武器、防护装备、航空航天、机械工程领域。
本发明的目的是通过下述技术方案实现的。
本发明公开的一种基于局部曲率自适应修正的改进水平集方法，实现方法如下：基于水平集Level Set方法，通过WENO有限差分格式和TVD Runge-Kutta格式对水平集函数进行离散。基于水平集函数的局部曲率自适应播撒示踪粒子，通过示踪粒子对水平集函数进行修正。WENO有限差分格式和TVD Runge-Kutta格式能实现水平集函数空间与时间项的高精度离散，保证界面水平集函数的计算精度；示踪粒子能在界面曲率变化较大的位置对水平集函数进行修正，确保多介质界面相互作用的数值仿真预测精度与效率。本发明有效的对复杂的工程问题多介质界面的运动进行高精度的数值模拟。
本发明公开的一种基于局部曲率自适应修正的改进水平集方法，包括如下步骤：
步骤1：确定界面运动计算区域，建立欧拉直角坐标系，在直角坐标系内进行网格划分，并设置边界条件及计算终止条件。
步骤1实现方法为：确定界面运动计算区域，建立欧拉直角坐标系，在直角坐标系内将计算区域划分为m*n个网格，其中，m表示x方向网格数量，n表示y方向网格数量。设置边界条件及计算终止条件。
步骤2：定义水平集Level Set函数，通过水平集函数确定界面位置。
步骤2实现方法为：定义计算区域内初始状态下的用于区分界面位置的水平集函数如公式(1)所示，水平集函数满足空间欧拉点距离界面最短的线段长度，正负号表示水平集函数所区分的不同的物质。
在计算初始时刻，采用公式(1)定义水平集函数的初值：
其中，x、y分别为网格节点的物理空间坐标，d表示网格节点(x,y)到界面的距离，Ω1、Ω2分别为两种物质界面两侧的区域，Γ(0)为初始时刻两种物质组成的界面。
水平集函数采用符号距离实现在任何时刻对界面演化过程的追踪，界面始终处于符号距离函数为零的零等值面位置：
其中，Γ(t)为两种物质组成的运动界面，t为运动时间。
步骤3：根据步骤2定义的水平集函数分别求解不同区域法向量与曲率。
步骤3实现方法为：根据步骤2定义的水平集函数分别求解不同区域法向量与曲率，求解公式分别为公式(3)和(4)：
其中，
为空间坐标点的梯度方向法向量，k为该点的空间曲率。
步骤4：根据步骤2定义的水平集函数确定各区域，并根据需要定义不同区域对应的背景速度场。
步骤4实现方法为：根据步骤2定义的水平集函数确定各区域，并根据需要定义不同区域对应的背景速度场u(u,v)，其中u为x方向的速度分量,v为y方向的速度分量。
步骤5：根据步骤3计算的求解区域内各网格点曲率，分别得到全场曲率与局部曲率。在界面附近三个网格宽度区域内，根据网格点的全场曲率与局部曲率初始化自适应播撒示踪粒子，实现在曲率相对较大的区域布置较多的示踪粒子，在曲率相对较小的区域布置较少的示踪粒子。定义示踪粒子半径范围，通过插值确定示踪粒子半径及正负，通过随机算法，将粒子在界面法向方向进行均匀排布。确保界面附近更多的粒子与界面相切，增加界面重构精度。进而可对曲率相对较大的易产生质量损失的区域产生有效的修正，提高模拟预测精度和效率。
步骤5具体实现方法为：根据步骤3计算求解区域的全场曲率与局部曲率，即基于步骤3计算得出全流场曲率绝对值的平均值与当前网格曲率的绝对值，对于界面附近三个网格宽度内的网格，单一计算网格内需布置的粒子数为：
其中，
为全流场曲率绝对值的平均值，|κ(i,j)|为当前格点曲率的绝对值，i，j分别表示计算区域内x，y方向的坐标索引，floor为向下取整函数，即求不大于要求值的最大整数值函数。p,q为整数。
为了提高数值模拟精度，作为优选，公式(5)中，p为3，q为4，为：
为了避免奇异点处曲率的奇异性导致过多示踪粒子的引入，采用如下限制器控制示踪粒子数量：
Nneed≤24。 (7)
对于示踪粒子，其半径rp的最小值与最大值分别为：
rmin＝0.1min(Δx,Δy)， (8)
rmax＝0.5min(Δx,Δy)。 (9)
示踪粒子初始均匀的布置在水平集函数绝对值小于3max(Δx,Δy)的欧拉网格内，其中，Δx,Δy分别为x，y方向的网格宽度。界面两侧粒子分布层的最小与最大厚度分别为：
bmin＝rmin， (10)
bmax＝3·max(Δx,Δy)。 (11)
在保证界面切向粒子近似均匀分布的同时，有必要通过随机算法，将粒子在界面法向方向进行均匀排布，最终将有效粒子所处位置插值出的水平集函数的绝对值
作为该粒子的半径rp，通过
标记粒子的正负sp，确保界面附近更多的粒子与界面相切，增加界面重构精度，即示踪粒子半径表示为：
在示踪粒子初始化完成后，将其标记为初始化粒子。
步骤6：根据步骤1确定的网格宽度及步骤4确定的速度参数，选取CFL参数，得到计算时间步长。
选取计算参数CFL计算时间步长，CFL参数为介于0到1之间的常数，选取计算参数CFL计算时间步长：
步骤7：对水平集函数的空间项进行离散，对单一介质计算区域进行时间离散，从而推进水平集函数，得到下一个时间步的水平集函数。
为了提高数值模拟精度，作为优选，步骤7中采用WENO有限差分格式对水平集函数的空间项进行离散，采用TVD Runge-Kutta格式对单一介质计算区域进行时间离散，从而推进水平集函数，得到下一个时间步的水平集函数。
步骤7实现方法为：
步骤7.1：对步骤2定义的水平集函数进行推进求解，具体推进的对流方程如下：
步骤7.2：对公式(14)的空间导数采用如下方式进行离散：
其中，
为水平集函数在具体计算节点i上的x方向空间导数，
为计算节点i上的x方向向前差分空间导数，
为计算节点i上的x方向向后差分空间导数，ui为计算节点i上x方向速度分量。
优选五阶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
其中，ε为常数。对
采用同样的方法进行离散。
步骤7.3：采用TVD Runge-Kutta格式对计算区域进行时间离散。优选三阶TVDRunge-Kutta格式对公式(14)的时间导数项进行离散，时间三阶TVD Runge-Kutta得到水平集函数的完全离散格式，进而求得出下一时间步计算区域内各网格单元的界面函数值。
其中，
分别为水平集函数的对应不同中间变量，
分别为初始时间步与下一时间步的水平集函数，L为关于水平集函数
的函数。
步骤8：对步骤4定义的速度场参数插值得到示踪粒子速度，结合步骤5定义的示踪粒子其他参数，推进示踪粒子的运动。
步骤8实现方法为：通过步骤4定义的速度场插值得到示踪粒子速度，结合步骤5定义的示踪粒子参数，通过式(21)直接推进布置在界面附近的示踪粒子的运动，修正网格点上的水平集函数。
其中，xp为示踪粒子的位置，u(xp)为示踪粒子的速度。
步骤9：根据步骤7与步骤8分别推进水平集函数与示踪粒子，根据新的水平集函数对示踪粒子进行插值，插值结果与步骤8得到的示踪粒子进行对比得到“逃逸粒子”，进而对逃逸粒子附近的水平集函数值进行修正。将正逃逸粒子与负逃逸粒子修正得到水平集函数值进行整合，增加鲁棒性。
在推进水平集函数及示踪粒子之后，通过新的水平集函数对示踪粒子进行插值，对于插值后示踪粒子对应正负符号
与初始示踪粒子正负符号Sp不同的粒子，标记为“逃逸粒子”，通过公式(22)对逃逸粒子附近点进行修正：
其中，
分别为正负逃逸粒子周围欧拉网格结点对应的水平集函数值。
对于每个粒子都重复以上步骤，并最终将正逃逸粒子与负逃逸粒子修正得到水平集函数值进行整合，为了鲁棒性，采用公式(23)：
其中，
分别为正负逃逸粒子修正过后的全场水平集函数，
为整合的全场水平集函数。
步骤10：在步骤7求解过程中水平集函数将失去其原有性质。通过求解重新初始化函数，使水平集函数保持其符号距离性质，从而提高计算精度。
在步骤7求解过程中，由于数值误差，水平集函数难以保持原有性质，界面附近点的误差会导致界面计算精度的大幅降低，甚至导致错误的结果，因此，通过重新初始化将全流场水平集函数进行修正：
其中，为每一步初始水平集函数，通过时间步τ进行虚拟推进。
为每一步初始水平集函数对应的正负号符号。
为了提高精度，作为优选，步骤10中时间步τ采用0.1dx进行虚拟推进。
步骤11：示踪粒子再次修正水平集函数。步骤10进行过程中，需保持示踪粒子静止，再用逃逸的示踪粒子去判断重新初始化产生的误差，并采用步骤9中的方法修正水平集函数。
步骤12：对示踪粒子重置。对步骤11修正后的水平集函数进行插值，采用修正后的粒子的水平集函数对应的插值对示踪粒子的半径进行重新赋值，重新赋值后的粒子半径应继续保持与界面相切。
步骤13：根据步骤12中的水平集函数增加或删除动态修正粒子，并对增加的动态修正粒子按照步骤5确定半径及正负信息，从而确保界面持续稳定的追踪。
在界面发生过度拉伸的位置，通过局部曲率计算得到该网格内需要存在示踪粒子的总数Nneed(i,j)，通过与当前网格内存在粒子总数Ntotal(i,j)对比，确定额外引入的修正粒子数量Ndynamic(i,j)，随机布置在网格内，后续添加的示踪粒子标记为动态修正粒子。当界面发生缩短扭曲等情况时，若通过公式(6)计算所得到的当前网格示踪粒子所需数量Nneed(i,j)小于当前网格内初始化粒子Ninitial(i,j)与动态修正粒子数量Ndynamic(i,j)的总和，则通过随机算法删除网格内部分或全部动态修正粒子，最终保证网格内粒子数量维持在预设较低水平。
在添加动态修正粒子之后，及时通过当前粒子插值得到的水平集函数值
初始化动态修正粒子半径rp及粒子正负符号Sp等信息，确保界面持续稳定的追踪。
步骤14：重复迭代步骤7至步骤13，直至达到最终计算步数或最终计算时间，实现不同介质界面追踪的高精度数值模拟预测。
还包括步骤15：根据步骤1至步骤14所述的一种基于局部曲率自适应修正的改进水平集方法进行多介质相互作用技术领域工程技术问题高精度数值模拟与预测，解决多介质相互作用技术领域相应工程技术问题。
步骤15所述多介质相互作用技术领域包括爆炸冲击、高速侵彻武器、防护装备、航空航天、机械工程领域。
所述高速侵彻武器工程技术问题包括激波打气泡、弹塑性块体撞击、聚能射流问题。
所述防护装备工程技术问题包括防护装备设计与制造、性能测试问题。
所述航空航天工程技术问题包括机翼绕流、激波模拟、紊流问题。
所述机械工程技术问题包括精确制造及其相关模拟问题。
作为优选，所述一种基于局部曲率自适应修正的改进水平集方法中，通过增加水平集函数维度，将二维x、y拓展至三维x、y、z，对于多介质界面，通过多个水平集函数确定，实现三维多介质数值模拟与预测。
有益效果：
1、本发明公开的一种基于局部曲率自适应修正的改进水平集方法，采用基于局部曲率的粒子重新播撒及修正方法，在提高精度的同时，大幅提升计算效率。基于界面曲率，将粒子连续的布置在界面附近，在角点等曲率较大的位置布置更多的粒子以达到更高的修正效率。
2、本发明公开的一种基于局部曲率自适应修正的改进水平集方法，能够有效的处理界面拉伸等大变形问题，动态修正粒子的引入确保界面在过度拉伸的情况下依然能够维持稳定的界面精度。
3、本发明公开的一种基于局部曲率自适应修正的改进水平集方法，由于初始化粒子携带的初始界面信息未被人为修正，避免信息的丢失，确保界面修正的准确性。具有更高的鲁棒性，并在保证界面精度的同时，提高计算效率。
4、本发明公开的一种基于局部曲率自适应修正的改进水平集方法，能够进行多介质相互作用技术领域工程技术问题高精度数值模拟与预测，解决多介质相互作用技术领域相应工程技术问题，所述多介质相互作用技术领域包括高速侵彻武器、防护装备、航空航天、机械工程领域。所述高速侵彻武器工程技术问题包括激波打气泡、弹塑性块体撞击、聚能射流问题。所述防护装备工程技术问题包括防护装备设计与制造、性能测试问题。所述航空航天工程技术问题包括机翼绕流、紊流、激波模拟问题。所述机械工程技术问题包括精确制造及其相关模拟问题。
附图说明
图1是本发明的一种基于局部曲率自适应修正的改进水平集方法流程图；
图2是本发明步骤2五角星形液滴界面示意图；
图3是原始方法粒子初始化云图；
图4是本发明步骤5自适应修正方法粒子初始化云图；
图5是本发明自适应修正方法与原始方法最终时间步界面对比示意图。
具体实施方式
实施例1：
为了更好的说明本发明的目的和优点，下面结合附图与实施例对本发明作进一步说明。
本实施例为基于局部曲率自适应修正的改进水平集方法对于五角星形液滴界面的数值模拟算例中的应用，五角星形液滴算例由于界面曲率为连续分布，是很好的检验自适应修正粒子算法性能的算例。
本实施例中所涉及单位制均为国际单位制。
如图1所示，本实施例公开一种基于局部曲率自适应修正的改进水平集方法，具体实现方法为：
步骤1：确定多介质界面运动计算区域，建立直角坐标系，在直角坐标系内进行网格划分，并设置边界条件及计算终止条件。
计算区域为[0,1]×[0,1]，五角星形液滴中心坐标为[0.5,0.5]，网格数为200×200。同时将上、下、左、右边界条件均设置为出流边界条件，保正物理模型的正确性。设置总时长为16s，在时间为8s时速度场反向，在时间为16s时，水平集函数零等值线回到原位置。
步骤2：定义水平集Level Set函数。初始水平集函数在极坐标下为：
r(θ)＝r+acos(bθ)， (25)
其中，r＝0.3，a＝0.1，b＝5，初始界面函数如附图2所示。
步骤3：根据步骤2定义的水平集函数分别求解不同区域法向量与曲率。
步骤4：定义背景速度场。背景速度流场为旋转速度场：
步骤5：根据步骤3计算求解区域内各网格点的曲率，得到全场曲率与局部曲率。在界面附近三个网格宽度区域内，根据网格点的全场曲率与局部曲率初始化自适应播撒示踪粒子。
在进行步骤5时，原始方法与基于局部曲率自适应修正的改进水平集方法的粒子初始化云图如图3、图4所示。
步骤6：选取CFL参数，得到计算时间步长。为对比传统方法和改进的粒子方法，采用统一的时间步0.001s，其满足CFL条件。
步骤7：对水平集函数的空间项进行离散，对单一介质计算区域进行时间离散，从而推进水平集函数，得到下一个时间步的水平集函数。
步骤8：对步骤4定义的速度场参数插值得到示踪粒子速度，结合步骤5定义的示踪粒子其他参数，推进示踪粒子的运动。
步骤9：根据步骤7与步骤8分别推进水平集函数与示踪粒子，根据新的水平集函数对示踪粒子进行插值，插值结果与步骤8得到的示踪粒子进行对比得到“逃逸粒子”，进而对逃逸粒子附近的水平集函数值进行修正。将正逃逸粒子与负逃逸粒子修正得到符号距离函数值进行整合，增加鲁棒性。
步骤10：在步骤7求解过程中水平集函数将失去其原有性质。通过求解重新初始化函数，使水平集函数保持其符号距离性质，从而提高计算精度。
步骤11：示踪粒子再次修正符号距离函数。
步骤12：对示踪粒子重置。
步骤13：根据步骤12中的水平集函数增加或删除动态修正粒子，并对增加的动态修正粒子按照步骤5确定半径及正负信息，从而确保界面持续稳定的追踪。
步骤14：重复迭代步骤7至步骤13，直至达到最终计算步数或最终计算时间，实现不同介质界面追踪的高精度数值模拟预测
完成步骤14后，原始的传统的带示踪粒子的水平集方法与基于局部曲率自适应修正的改进水平集方法的计算时间与质量损失见表1。为方便对比，不同计算方法所需示踪粒子见表1。不同计算方法在计算结束后，其界面对比如图5所示。
表1不同方法五角星形液滴算例对比
仿真预测结果分析：
图4为自适应修正粒子算法初始粒子分布图。在距离中心更近的位置处，由于曲率最大，粒子数相较于传统方法(图3)明显变多，而在接近直线的曲边附近，粒子保持在一个较低的水平。图5为初始曲线与原始方法、自适应修正粒子算法计算结果对比图，在曲率大的位置，自适应修正粒子方法达到比原始方法更好的计算结果，效率提高约50％。
以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (10)
1.一种基于局部曲率自适应修正的改进水平集方法，其特征在于：包括如下步骤，
步骤1：确定界面运动计算区域，建立欧拉直角坐标系，在直角坐标系内进行网格划分，并设置边界条件及计算终止条件；
步骤2：定义水平集Level Set函数，通过水平集函数确定界面位置；
步骤3：根据步骤2定义的水平集函数分别求解不同区域法向量与曲率；
步骤4：根据步骤2定义的水平集函数确定各区域，并根据需要定义不同区域对应的背景速度场；
步骤5：根据步骤3计算的求解区域内各网格点曲率，分别得到全场曲率与局部曲率；在界面附近三个网格宽度区域内，根据网格点的全场曲率与局部曲率初始化自适应播撒示踪粒子，实现在曲率相对较大的区域布置较多的示踪粒子，在曲率相对较小的区域布置较少的示踪粒子；定义示踪粒子半径范围，通过插值确定示踪粒子半径及正负，通过随机算法，将粒子在界面法向方向进行均匀排布；确保界面附近更多的粒子与界面相切，增加界面重构精度；进而可对曲率相对较大的易产生质量损失的区域产生有效的修正，提高模拟预测精度和效率；
步骤6：根据步骤1确定的网格宽度及步骤4确定的速度参数，选取CFL参数，得到计算时间步长；
步骤7：对水平集函数的空间项进行离散，对单一介质计算区域进行时间离散，从而推进水平集函数，得到下一个时间步的水平集函数；
步骤8：对步骤4定义的速度场参数插值得到示踪粒子速度，结合步骤5定义的示踪粒子其他参数，推进示踪粒子的运动；
步骤9：根据步骤7与步骤8分别推进水平集函数与示踪粒子，根据新的水平集函数对示踪粒子进行插值，插值结果与步骤8得到的示踪粒子进行对比得到“逃逸粒子”，进而对逃逸粒子附近的水平集函数值进行修正；将正逃逸粒子与负逃逸粒子修正得到水平集函数值进行整合，增加鲁棒性；
步骤10：在步骤7求解过程中水平集函数将失去其原有性质；通过求解重新初始化函数，使水平集函数保持其符号距离性质，从而提高计算精度；
步骤11：示踪粒子再次修正水平集函数；步骤10进行过程中，需保持示踪粒子静止，再用逃逸的示踪粒子去判断重新初始化产生的误差，并采用步骤9中的方法修正水平集函数；
步骤12：对示踪粒子重置；对步骤11修正后的水平集函数进行插值，采用修正后的粒子的水平集函数对应的插值对示踪粒子的半径进行重新赋值，重新赋值后的粒子半径应继续保持与界面相切；
步骤13：根据步骤12中的水平集函数增加或删除动态修正粒子，并对增加的动态修正粒子按照步骤5确定半径及正负信息，从而确保界面持续稳定的追踪；
步骤14：重复迭代步骤7至步骤13，直至达到最终计算步数或最终计算时间，实现不同介质界面追踪的高精度数值模拟预测。
2.如权利要求1所述的一种基于局部曲率自适应修正的改进水平集方法，其特征在于：还包括步骤15：根据步骤1至步骤14所述的一种基于局部曲率自适应修正的改进水平集方法进行多介质相互作用技术领域工程技术问题高精度数值模拟与预测，解决多介质相互作用技术领域相应工程技术问题。
3.如权利要求2所述的一种基于局部曲率自适应修正的改进水平集方法，其特征在于：步骤15所述多介质相互作用技术领域包括高速侵彻武器、防护装备、航空航天、机械工程领域；
所述高速侵彻武器工程技术问题包括激波打气泡、弹塑性块体撞击、聚能射流问题；
所述防护装备工程技术问题包括防护装备设计与制造、性能测试问题；
所述航空航天工程技术问题包括机翼绕流、激波模拟、紊流问题；
所述机械工程技术问题包括精确制造及其相关模拟问题。
4.如权利要求3所述的一种基于局部曲率自适应修正的改进水平集方法，其特征在于：通过增加水平集函数维度，将二维x、y拓展至三维x、y、z，对于多介质界面，通过多个水平集函数确定，实现三维多介质数值模拟与预测。
5.如权利要求1、2或3所述的一种基于局部曲率自适应修正的改进水平集方法，其特征在于：
步骤1实现方法为，确定界面运动计算区域，建立欧拉直角坐标系，在直角坐标系内将计算区域划分为m*n个网格，其中，m表示x方向网格数量，n表示y方向网格数量；设置边界条件及计算终止条件。
步骤2实现方法为，定义计算区域内初始状态下的用于区分界面位置的水平集函数如公式(1)所示，水平集函数满足空间欧拉点距离界面最短的线段长度，正负号表示水平集函数所区分的不同的物质；
在计算初始时刻，采用公式(1)定义水平集函数的初值：
其中，x、y分别为网格节点的物理空间坐标，d表示网格节点(x,y)到界面的距离，Ω1、Ω2分别为两种物质界面两侧的区域，Γ(0)为初始时刻两种物质组成的界面；
水平集函数采用符号距离实现在任何时刻对界面演化过程的追踪，界面始终处于符号距离函数为零的零等值面位置：
其中，Γ(t)为两种物质组成的运动界面，t为运动时间；
步骤3实现方法为，根据步骤2定义的水平集函数分别求解不同区域法向量与曲率，求解公式分别为公式(3)和(4)：
其中，
为空间坐标点的梯度方向法向量，k为该点的空间曲率；
步骤4实现方法为，根据步骤2定义的水平集函数确定各区域，并根据需要定义不同区域对应的背景速度场u(u,v)，其中u为x方向的速度分量,v为y方向的速度分量。
6.如权利要求5所述的一种基于局部曲率自适应修正的改进水平集方法，其特征在于：步骤5具体实现方法为，根据步骤3计算求解区域的全场曲率与局部曲率，即基于步骤3计算得出全流场曲率绝对值的平均值与当前网格曲率的绝对值，对于界面附近三个网格宽度内的网格，单一计算网格内需布置的粒子数为：
其中，
为全流场曲率绝对值的平均值，|κ(i,j)|为当前格点曲率的绝对值，i，j分别表示计算区域内x，y方向的坐标索引，floor为向下取整函数，即求不大于要求值的最大整数值函数；p,q为整数；
为了提高数值模拟精度，作为优选，公式(5)中，p为3，q为4，为：
为了避免奇异点处曲率的奇异性导致过多示踪粒子的引入，采用如下限制器控制示踪粒子数量：
Nneed≤24； (7)
对于示踪粒子，其半径rp的最小值与最大值分别为：
rmin＝0.1min(Δx,Δy)， (8)
rmax＝0.5min(Δx,Δy)； (9)
示踪粒子初始均匀的布置在水平集函数绝对值小于3max(Δx,Δy)的欧拉网格内，其中，Δx,Δy分别为x，y方向的网格宽度；界面两侧粒子分布层的最小与最大厚度分别为：
bmin＝rmin， (10)
bmax＝3·max(Δx,Δy)； (11)
在保证界面切向粒子近似均匀分布的同时，有必要通过随机算法，将粒子在界面法向方向进行均匀排布，最终将有效粒子所处位置插值出的水平集函数的绝对值
作为该粒子的半径rp，通过
标记粒子的正负sp，确保界面附近更多的粒子与界面相切，增加界面重构精度，即示踪粒子半径表示为：
在示踪粒子初始化完成后，将其标记为初始化粒子。
7.如权利要求6所述的一种基于局部曲率自适应修正的改进水平集方法，其特征在于：步骤六实现方法为，
选取计算参数CFL计算时间步长，CFL参数为介于0到1之间的常数，选取计算参数CFL计算时间步长：
8.如权利要求7所述的一种基于局部曲率自适应修正的改进水平集方法，其特征在于：为了提高数值模拟精度，步骤7中采用WENO有限差分格式对水平集函数的空间项进行离散，采用TVD Runge-Kutta格式对单一介质计算区域进行时间离散，从而推进水平集函数，得到下一个时间步的水平集函数；
步骤7实现方法为，
步骤7.1：对步骤2定义的水平集函数进行推进求解，具体推进的对流方程如下：
步骤7.2：对公式(14)的空间导数采用如下方式进行离散：
其中，
为水平集函数在具体计算节点i上的x方向空间导数，
为计算节点i上的x方向向前差分空间导数，
为计算节点i上的x方向向后差分空间导数，ui为计算节点i上x方向速度分量；
选五阶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
其中，ε为常数；对
采用同样的方法进行离散；
步骤7.3：采用TVD Runge-Kutta格式对计算区域进行时间离散；优选三阶TVD Runge-Kutta格式对公式(14)的时间导数项进行离散，时间三阶TVD Runge-Kutta得到水平集函数的完全离散格式，进而求得出下一时间步计算区域内各网格单元的界面函数值；
其中，
分别为水平集函数的对应不同中间变量，
分别为初始时间步与下一时间步的水平集函数，L为关于水平集函数
的函数。
9.如权利要求8所述的一种基于局部曲率自适应修正的改进水平集方法，其特征在于：
步骤8实现方法为，通过步骤4定义的速度场插值得到示踪粒子速度，结合步骤5定义的示踪粒子参数，通过式(21)直接推进布置在界面附近的示踪粒子的运动，修正网格点上的水平集函数；
其中，xp为示踪粒子的位置，u(xp)为示踪粒子的速度；
在推进水平集函数及示踪粒子之后，通过新的水平集函数对示踪粒子进行插值，对于插值后示踪粒子对应正负符号
与初始示踪粒子正负符号Sp不同的粒子，标记为“逃逸粒子”，通过公式(22)对逃逸粒子附近点进行修正：
其中，
分别为正负逃逸粒子周围欧拉网格结点对应的水平集函数值；
对于每个粒子都重复以上步骤，并最终将正逃逸粒子与负逃逸粒子修正得到水平集函数值进行整合，为了鲁棒性，采用公式(23)：
其中，
分别为正负逃逸粒子修正过后的全场水平集函数，
为整合的全场水平集函数。
10.如权利要求9所述的一种基于局部曲率自适应修正的改进水平集方法，其特征在于：
步骤10实现方法为，
在步骤7求解过程中，由于数值误差，水平集函数难以保持原有性质，界面附近点的误差会导致界面计算精度的大幅降低，甚至导致错误的结果，因此，通过重新初始化将全流场水平集函数进行修正：
其中，
为每一步初始水平集函数，通过时间步τ进行虚拟推进；
为每一步初始水平集函数对应的正负号符号；
步骤13实现方法为，在界面发生过度拉伸的位置，通过局部曲率计算得到该网格内需要存在示踪粒子的总数Nneed(i,j)，通过与当前网格内存在粒子总数Ntotal(i,j)对比，确定额外引入的修正粒子数量Ndynamic(i,j)，随机布置在网格内，后续添加的示踪粒子标记为动态修正粒子；当界面发生缩短扭曲等情况时，若通过公式(6)计算所得到的当前网格示踪粒子所需数量Nneed(i,j)小于当前网格内初始化粒子Ninitial(i,j)与动态修正粒子数量Ndynamic(i,j)的总和，则通过随机算法删除网格内部分或全部动态修正粒子，最终保证网格内粒子数量维持在预设较低水平；
在添加动态修正粒子之后，及时通过当前粒子插值得到的水平集函数值
初始化动态修正粒子半径rp及粒子正负符号Sp等信息，确保界面持续稳定的追踪。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010091056.8A CN111259325A (zh) | 2020-02-13 | 2020-02-13 | 一种基于局部曲率自适应修正的改进水平集方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010091056.8A CN111259325A (zh) | 2020-02-13 | 2020-02-13 | 一种基于局部曲率自适应修正的改进水平集方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN111259325A true CN111259325A (zh) | 2020-06-09 |
Family
ID=70951284
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010091056.8A Pending CN111259325A (zh) | 2020-02-13 | 2020-02-13 | 一种基于局部曲率自适应修正的改进水平集方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111259325A (zh) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112307669A (zh) * | 2020-11-25 | 2021-02-02 | 北京理工大学 | 一种基于实际物理点迭代修正的多介质水平集方法 |
| CN112613251A (zh) * | 2020-12-29 | 2021-04-06 | 中国航天空气动力技术研究院 | 一种介质类型分步标记方法 |
| CN113378440A (zh) * | 2021-06-23 | 2021-09-10 | 四川大学 | 一种流固耦合数值模拟计算方法、装置及设备 |
| CN115438604A (zh) * | 2022-11-08 | 2022-12-06 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种基于质数体系的网格标识方法 |
-
2020
- 2020-02-13 CN CN202010091056.8A patent/CN111259325A/zh active Pending
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112307669A (zh) * | 2020-11-25 | 2021-02-02 | 北京理工大学 | 一种基于实际物理点迭代修正的多介质水平集方法 |
| CN112307669B (zh) * | 2020-11-25 | 2023-04-07 | 北京理工大学 | 一种基于实际物理点迭代修正的多介质水平集方法 |
| CN112613251A (zh) * | 2020-12-29 | 2021-04-06 | 中国航天空气动力技术研究院 | 一种介质类型分步标记方法 |
| CN112613251B (zh) * | 2020-12-29 | 2023-07-28 | 中国航天空气动力技术研究院 | 一种介质类型分步标记方法 |
| CN113378440A (zh) * | 2021-06-23 | 2021-09-10 | 四川大学 | 一种流固耦合数值模拟计算方法、装置及设备 |
| CN115438604A (zh) * | 2022-11-08 | 2022-12-06 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种基于质数体系的网格标识方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111259325A (zh) | 一种基于局部曲率自适应修正的改进水平集方法 | |
| Wang | Adaptive high-order methods in computational fluid dynamics | |
| CN111241742B (zh) | 一种多相流计算方法 | |
| CN109726465A (zh) | 基于非结构曲边网格的三维无粘低速绕流的数值模拟方法 | |
| CN112613243B (zh) | 一种流体力学模拟的方法、装置及计算机可读存储介质 | |
| CN112947534A (zh) | 一种高超声速飞行器下压段自适应伪谱法轨迹优化方法 | |
| CN112001109A (zh) | 再生核粒子算法实现结构冲击动力学仿真方法 | |
| CN109657408A (zh) | 一种再生核粒子算法实现结构线性静力学仿真方法 | |
| CN109492234B (zh) | 一种改进的流固耦合插值方法 | |
| CN113536461A (zh) | 用于高超声速强激波流场气动热预测的湍流模型修正方法 | |
| CN112307669B (zh) | 一种基于实际物理点迭代修正的多介质水平集方法 | |
| CN107944113A (zh) | 一种计算三维高速平动目标电磁散射场的方法 | |
| Hui | The unified coordinate system in computational fluid dynamics | |
| CN105718619A (zh) | 一种基于有限元法的飞行器燃油质量特性确定方法 | |
| CN109726496B (zh) | 一种基于iisph提高不可压缩水体模拟效率的实现方法 | |
| CN105987695A (zh) | 一种高速旋转弹姿态解算用四区间拉格朗日方法 | |
| CN116522660A (zh) | 一种用于爆炸问题的数值模拟方法及系统 | |
| CN112307668B (zh) | 一种基于粒子的粘液类效果模拟方法 | |
| CN115730533A (zh) | 一种侧向喷流干扰流场模拟计算的网格自适应混合判据 | |
| Kutler | Application of selected finite difference techniques to the solution of conical flow problems | |
| CN112560326B (zh) | 压力场的确定方法及装置 | |
| Gao et al. | A local frame based hexahedral mesh optimization | |
| CN109684766A (zh) | 含转角的大变形柔性梁单元建模方法 | |
| CN113886964B (zh) | 基于边界替代混合模型的飞行器气动分析方法 | |
| Kim et al. | Three-dimensional building-cube method for inviscid compressible flow computations |
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 | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20200609 |
|
| RJ01 | Rejection of invention patent application after publication |