基于八叉树自适应网格技术的level set运动界面追踪方法-殷亚军.pdf

《基于八叉树自适应网格技术的level set运动界面追踪方法-殷亚军.pdf》由会员分享，可在线阅读，更多相关《基于八叉树自适应网格技术的level set运动界面追踪方法-殷亚军.pdf（10页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、2016 年 11 月 CIESC Journal 4732 November 2016第 67 卷 第 11 期 化 工 学 报 Vol.67 No.11DOI： 10.11949/j.issn.0438-1157.20160421 基于八叉树自适应网格技术的 Level Set 运动界面追踪方法 殷亚军，李阳东，涂志新，李文，沈旭，周建新 （华中科技大学材料成形与模具技术国家重点实验室，湖北 武汉 430074） 摘要 ： Level Set 方法因能有效地处理界面处复杂的拓扑结构变化以及大变形问题，广泛应用于界面追踪领域。在Level Set 方法追踪运动界面时引入八叉树网格技术，通过八

2、叉树网格的细化和粗化技术减少计算网格数量和计算内存并提高计算效率和计算精度。因为八叉树网格为非均匀网格，其相邻网格的层数值可能不相同，所以不能直接采用 WENO 格式离散 Level Set 函数得到网格处的函数值，进而提出八叉树网格离散模型解决这一问题，并提出基于八叉树网格距离场重新初始化方法减少 Level Set 方法的质量损失， 最后将基于八叉树网格技术的 Level Set方法应用于两个给定速度场的运动界面模拟算例以及基准件方腔的铸造充型过程的模拟。模拟结果表明该方法可以提高界面的精度，同时改善质量守恒性。 关键词 ：界面追踪；数值模拟；八叉树网格技术； Level Set 方法；网

3、格自适应；模型；优化 中图分类号 ： O 351.2 文献标志码 ： A 文章编号 ： 0438 1157（ 2016） 11 4732 10 Interface capturing method based on Level Set with octree adaptive mesh technology YIN Yajun, LI Yangdong, TU Zhixin, LI Wen, SHEN Xu, ZHOU Jianxin (State Key Laboratory of Materials Processing and Die & Mould Technology, Huazho

4、ng University of Science and Technology, Wuhan 430074, Hubei, China) Abstract: Level Set method is one of the dominant methods to calculate and track evolution of interfaces because it can efficiently handle complex topology change of interfaces and large deformation. This method often suffers from

5、large amount of numerical dissipation and high order accurate numerical discretization in time and space is used to improve mass conservation properties, which leads to an increase of computational overhead. When octree adaptive mesh technology is incorporated into Level Set to capture moving interf

6、aces, the amount of computational meshes and memory usage is reduced and the calculation precision and efficiency are enhanced by adaptive octree mesh refining or coarsening according to interface distance. Because octree mesh is non-uniform and layer numbers of neighbor grids may be different, the

7、WENO formula could not be directly applied to set values of octree grids. A discretization model and a technique of re-initializing distance field were proposed to calculate value of neighbor octree grids as well as to reduce mass loss in Level Set function. Numerical simulation on moving interface

8、in rotation flow field of Zalesak sphere, in shear flow field of sphere and in cast filling process of square cavity benchmark indicated that combination of Level Set method and octree adaptive mesh technology could successfully capture moving interface. The interface-tracking precision was enhanced

9、 by refinement of octree meshes and increase of maximum octree mesh layers. When the maximum octree mesh layer 2016-04-05 收到初稿， 2016-08-03 收到修改稿。 联系人 ：周建新。 第一作者 ：殷亚军（ 1985） ，男，博士，讲师。基金项目 ：国家自然科学基金项目（ 51475181， 51305149） ；中国博士后科学基金特别资助项目（ 2015T80795） 。 Received date: 2016-04-05. Corresponding author:

10、 Prof. ZHOU Jianxin, Foundation item: supported. by the National Natural Science Foundation of China(51475181, 51305149) and the China Postdoctoral Science Foundation Special Funding(2015T80795). 万方数据第 11 期 殷亚军等：基于八叉树自适应网格技术的 Level Set 运动界面追踪方法 4733was increased from 8 to 9 in the shear flow field,

11、mass conservation of the sphere was obviously improved and numerical dissipation caused by Level Set method was significantly reduced. Key words: interface capture; numerical simulation; octree mesh technology; Level Set method; adaptive mesh; model; optimization 引 言 许多物理问题会涉及气液自由界面。采用数值模拟计算自由界面时，界面

12、附近流体的物性会发生明显变化， 需要准确追踪自由界面的位置。 Level Set（水平集）方法是美国计算数学学者 Osher 等1-4提出的一种追踪界面的方法，该方法能够有效地解决复杂界面追踪问题，是目前界面追踪问题数值模拟的主流方法。 Level Set 方法可以克服一般界面前沿追踪方法必须构造出具体的界面而难以处理界面前沿发生的拓扑结构变化的弱点，不需要显式追踪运动界面，从而可以较容易地处理复杂的物质界面及其拓扑结构发生变化的情形，便于计算曲率及界面法向量，可以将间断问题转化为光滑问题加以解决。 与 VOF（ volume of fluid）法相比，该方法易于从二维问题推广到三维问题，追踪

13、界面精度较高。 但是，经典的 Level Set 方法守恒性差，初始时刻 Level Set 为距离函数， 经过有限的时间步长后不再保持距离函数的性质，需要重新初始化，重新初始化 Level Set 函数涉及的复杂数学过程极大地增加了计算量。 此外， 为了减少数值耗散， 在求解 Level Set 方程时通常采用高阶格式离散空间导数和时间导数，这些高阶计算格式会影响其计算效率。为了提高 Level Set 方法的计算速率与精度， 许多学者对此做了大量的研究。 Osher 等5对 Level Set 重新初始化方法进行了研究。 Sethian 等提出了 Fast Marching Method

14、(FMM，快速步进法 )6-7对 Level Set 函数进行重新初始化，提出窄带方法8使 Level Set 方法更加快速有效。 Peng 等9在全局 Level Set方法的基础上提出了一种快速局部 Level Set 方法，大大节省了计算时间。 Strain10首先引入树型结构的自适应网格，在交界面附近采用高分辨率网格，而在远离交界面的区域保持较粗的网格。 Sochnikov等11将 Fast Marching（快速步进）技术应用到四叉树型结构的自适应网格中， 使得计算量进一步减少。Costa 等12采用了自适应网格加密（ adaptive mesh refinement, AMR）技术

15、进行结构拓扑优化。宫翔飞 等13将 Level Set 方法与 AMR 技术融合用于数值模拟。黄筱云等提出了一种基于树型结构空间自适应的快速粒子 Level Set 方法14，并提出自适应四叉树网格下的 N-S 方程数值求解模型15。 王生辉等16将四叉树直角坐标网格技术与 Level Set 方法相结合进行运动界面的追踪。 Losasso 等17利用八叉树网格技术模拟了烟雾和水流的运动。朱雷等18采用八叉树网格模拟焊接温度场，对比了八叉树网格与均匀网格的计算效率，采用八叉树网格的计算时间只有均匀网格的 40%。 为了进一步提高 Level Set 方法的计算效率以及界面精度，本研究将八叉树自

16、适应网格技术应用于 Level Set 方法的求解，提出适用于 Level Set 方法追踪界面的八叉树自适应网格模型，并建立了基于八叉树自适应网格的 Level Set 函数的离散和重新初始化模型。采用八叉树网格容易进行网格的局部加密或粗化，减少网格数量以及计算量。同时，局部网格细化能够很好地提高界面追踪的准确性。本研究增加界面处网格的最大层数值，可以很大程度地提高界面追踪的精度以及减少 Level Set 方法引起的数值耗散，提高界面追踪的质量守恒性，这对于改善流场失真具有重大的意义。 1 八叉树网格技术 八叉树是一种用于描述三维空间的树状数据结构，其原理是将空间区域不断分解为 8 个同样

17、大小的子区域，分解的次数越多，子区域就越小，一直到同一区域的属性单一为止19。图 1 所示即为一 图 1 八叉树模型 Fig.1 Octree model 万方数据 化 工 学 报 第 67 卷 4734 个简单的八叉树节点产生过程。图中左侧为对应的空间区域划分，不同的区域用不同颜色的正方体表示；右侧为八叉树树型结构对应的网格节点关系，图中黑色节点表示子区域的父区域，叶子表示八叉树网格的叶子单元网格。八叉树网格以八叉树作为基本理论基础，网格组织方式与之类似，树中指代的区域即为数值模拟中对应的网格（网格单元统一为正方体单元） 。 2 八叉树网格自适应 Level Set 方法追踪运动界面时， 界

18、面精度和质量损失都与网格的分辨率有关，如果将网格全部细化，则会很大程度地增加计算时间和存储空间。需要自适应网格技术对界面附近的区域进行网格细化处理，对远离界面的区域进行网格粗化处理，这样既可以保证界面附近网格有较高的分辨率，提高界面区域的计算精度，又可以减少远离界面区域所占用的内存，减少计算时间和计算内存。 2.1 细化模型 采用八叉树自适应网格技术可以提高界面附近区域的网格分辨率，保证界面附近区域的计算精度。图 2 所示即为一个简单的网格细化实例，其中假设白色区域为运动界面，而如果网格如图左侧所示，则需将界面周围网格细化（右图中红色网格） ，转化为右侧图所示，才能满足计算分析要求。 图 2

19、网格细化实例 Fig.2 Example of adaptive octree grids refinement 2.2 粗化模型 八叉树自适应网格与均匀网格的优势就在于效率高、计算网格数少，许多局部计算区域可以用粗化网格代替众多细化网格。图 3 即为某一计算时刻网格粗化的例子，白色环形区域为运动界面，左边部分网格不太理想，因为中心红色显示部分网格没有必要采用细化网格，必须经过一定的粗化模型粗化，如图 3 右边部分，这样就无须在该红色网格 图 3 网格粗化实例 Fig.3 Example of adaptive octree grids coarsening 区域浪费大量的计算资源。 3 Le

20、vel Set 方程的离散 Level Set 方法的主要思想是将随时间运动的物质界面定义为一个函数的零等值面（线） ，即(x,t)=0。在任意某一时刻，只要求出 值，并求出其零等值面，就能够知道此时的活动界面位置。 在两相流或自由界面追踪问题中，速度等物理量的控制方程是 N-S 方程， 所对应的 Level Set 方程见式（ 1） ，其中 表示 Level Set 函数。 0uvwtxyz+=(1) 本研究采用五阶 WENO 格式离散 Level Set 函数对空间的偏导数20，采用三阶 TVD Runge-Kutta格式离散 Level Set 函数对时间的偏导数21。 4 八叉树网格离

21、散模型 对于均匀网格来说，采用上述离散方式直接离散 Level Set 函数。但是八叉树自适应网格为非均匀网格，每个网格的邻居网格层数值（层数值表征该网格的细化或粗化程度，其中树根节点网格的层数值为 0，层数值越大网格单元的划分越细）可能不同，尽管单元数据仍然存储在正方体网格中心，但相邻单元网格的大小不同，不能直接采用上述的离散方式，需要对非均匀网格重新进行处理，进而获得邻居网格的 Level Set 函数值进行导数离散。 八叉树网格与其相邻网格之间存在 3 种关系：当前网格与相邻网格均为叶子单元，如图 4（ a）所示； 当前网格为叶子单元， 相邻网格为子单元网格，如图 4（ b）所示；当前网

22、格为子单元网格，相邻网格为叶子单元，如图 4（ c）所示。其中，红色网格为当前网格，绿色网格为相邻网格。 对于图 4（ a） 所示情况， 计算当前网格 Level Set函数的导数可以采用均匀网格计算方法得到，直接万方数据第 11 期 殷亚军等：基于八叉树自适应网格技术的 Level Set 运动界面追踪方法 4735选取相邻网格的 值进行计算。对于图 4（ b）所示情况， 在计算当前网格 Level Set 函数的导数时首先计算相邻网格值，在其相邻网格的 8 个子单元中找出与当前计算网格面相邻的 4 个子单元，因为当前网格在该面的物理场变量的变化只与该面的流通量有关，故需要这 4 个单元才能

23、获得相邻网格值。本研究采用一种较为简单的方式获取相邻网格值，将4 个子单元对应的中心网格值取平均值计算相邻网格的 Level Set 函数值， 然后用相邻网格值计算当前网格 Level Set 函数的空间导数。 相邻网格值的计算如式（ 2）所示 () () () ()01451/8 1/8 1/8 1/814iiiii + += (2) 式中， i+1表示导数计算时邻居网格值，()01/8i+、()11/8i+、()41/8i+和()51/8i+分别表示与当前计算网格面接触的 4 个邻居子单元网格值。 对于图 4（ c）所示情况，计算子单元网格的Level Set 函数值，需要借助一定的插值方

24、法获得邻居网格值，将相邻网格构造出 8 个相同的虚单元网格 （图 5 中绿色网格） ， 利用插值公式计算相邻虚单元网格的函数值并作为邻居网格值，然后用该虚单元的值计算当前子单元网格的 Level Set 函数导数值。邻居网格的插值公式如式（ 3）所示 0 11111iiiiix yzxyz+=+ + + (3) 式中， 0i+1表示相邻虚单元网格值， i+1表示相邻叶子单元网格值， x、 y、 z 分别表示虚子单元网格中心到相邻叶子单元网格的距离。 图 5 相邻虚单元网格 Fig.5 Illustration of neighbor virtural subunit grid 5 Level

25、Set 函数重新初始化 求解 Level Set 函数时，保持 (x,t)为距离函数非常重要，但是 Level Set 方法守恒性差，经过几个计算时间步长后 (x,t)不再保持距离函数性质，需要对 t 时刻求解得到的 (x,t)函数进行修正，使函数 的近似值满足距离场的条件。 Sussman 等3提出了一种新的重新初始化方法 0()(| |1)0tS += (4) 222()|()Sx=+ (5) 其中， S(0)是函数符号。 Peng 等9采用式（ 5）所示格式可以得到更好的效果。本研究采用式（ 5）形式对 Level Set 函数进行重新初始化。 重新初始化 Level Set 函数值会发

26、生改变， 所以八叉树网格中的 Level Set 函数值也需要更新， 八叉树网格中距离场函数更新步骤如下：将界面处网格单元（图 6 中紫色网格）标记为 1，其余单元标记为 0， 若边界位置网格不是界面单元则将其 Level Set 函数值设置为一个极大值；从边界位置网格图 6（ a）中绿色网格 扫描遍历网格单元，当扫描到细化网格单元如图 6（ b）红色圆圈中的网格 图 6（ b）为图 6（ a）中红色圈出网格的局部放大图 ，本研究规定先扫描与当前扫描网格相邻的 4 个子单元网格，再扫描与这 4 个子单元网格相邻的子单元网格，在未扫描到界面单元之前这些网格单元的距离场函数值不更新；当扫描到标记为

27、 1 的界面网格单元时，直接跳过该单元，不更新该网格单元的距离场函数值，继续沿该方向扫描并更新后面标记为 0 的网格单元的 Level Set 函数值， 直至扫描到下一个边界位置；从已更新函数值的边界网格 图 6（ a）中绿色圆圈中的网格 的相反方向开始扫描，在扫描到界面网格之前不更新已扫描网格的函数值，当扫描到界面网格单元时直接跳过该网格，继续扫描网格单元并更新之后扫描网格单元的函数值，直至扫描到界面网格为止，如果当前扫描网格图 4 八叉树当前网格与其相邻网格的关系 Fig.4 Relationship between current grid and neighbor grid 万方数据

28、化 工 学 报 第 67 卷 4736 的邻居网格有多个界面网格则该网格的距离场函数值为 3 个方向邻居网格的最小值；重复步骤对其他网格进行重新初始化，直到所有八叉树网格单元距离场函数值全部重新初始化为止。 6 应 用 以基于八叉树网格技术的 Level Set 方法模拟不同网格细化程度的 Zalesak 圆球在旋转场中的运动情况以及圆球在剪切场中的运动情况为例，讨论基于八叉树网格技术的 Level Set 方法在三维空间的应用，并将该方法应用于追踪方腔铸件充型过程的界面前沿，模拟充型过程中液体的流动。 6.1 旋转流场 本算例模拟的是一个半径为 0.25 的 Zalesak 圆球界面在旋转流

29、场作用下的变化情况，初始位置为0.5,0.5,0.5， 切口宽度为 0.125， 切口深度为 0.3125，计算区域为 0,1 0,1 0,1 ，网格尺寸为128 128 128。网格的最大层数分别为 8 层和 9层，最小网格层数为 4 层。旋转流场的速度场为 00.50.5uvzwy=(6) 流场周期 T 为 2，时间步长为 0.05。 图 7 为 Zalesak 圆球在三维旋转流场中八叉树自适应网格示意图,图 8和图 9分别是网格最大层数为 8 层和 9 层时 Zalesak 圆球在旋转流场中的运动时间分别为 0、 T/4、 T/2、 3T/4 和 T 时刻的数值模拟结果。 图 10 为刘

30、长春22采用 Youngs-VOF 方法模拟缺口圆在旋转流场中的运动的结果， 图 11 为张明远23采用均匀网格的 Level Set 方法追踪 Zalesak 圆球运动界面的结果。 通过图 8 和图 9 的数值模拟结果可以得出：Zalesak 圆球经过一个计算周期后回到了初始位置，八叉树网格最大层数为 9 层时 Zalesak 圆球表面的精度远高于网格最大层数为 8 层时 Zalesak 圆球表面的精度，而且在 Zalesak 圆球缺口处仍然可以保持尖锐。比较图 9、图 10 和图 11 的模拟结果可以 图 6 八叉树网格值扫描示意图 Fig.6 Illustration of sweepi

31、ng octree grid value图 7 Zalesak 圆球在旋转流场中运动的八叉树自适应网格模型 Fig.7 Octree adaptive mesh model of Zalesak sphere moving in rotation flow filed 图 8 网格最大层数为 8 层时 Zalesak 圆球的数值模拟结果 Fig.8 Numerical simulation results of Zalesak sphere while maximum octree mesh layer is 8 万方数据第 11 期 殷亚军等：基于八叉树自适应网格技术的 Level Set

32、运动界面追踪方法 4737图 11 基于均匀网格 Level Set 方法追踪旋转流场 中的运动界面23 Fig.11 Level Set method based on uniform mesh captures moving interface in rotation flow field23 发现：图 10 Youngs-VOF 方法中红色圆圈中的尖锐部分有一些失真，本研究提出的方法在尖锐界面处的模拟精度略好一些； 图 11 中绿色圆圈中界面部分质量损失严重且界面精度也较低，其中的平面部分已经平滑成为了曲面，本研究图 9 的模拟结果则很好地保持了原来的界面形状，提高了界面追踪的精度和质量

33、守恒性。 6.2 剪切流场 本算例模拟的是一个半径为 0.125 的圆球在剪切流场作用下的变化情况，初始位置为0.5,0.5,0.25，计算区域为 0,1 0,1 0,1，网格尺寸为 128 128 128。网格的最大层数分别为 8层和 9 层，最小层数为 4 层。剪切流场的速度场为 图 9 网格最大层数为 9 层时 Zalesak 圆球的数值模拟结果 Fig.9 Numerical simulation results of Zalesak sphere while maximum octree mesh layer is 9 图 10 Youngs-VOF 方法追踪旋转流场中的运动界面22

34、 Fig.10 Youngs-VOF method captures moving interface in rotation flow field22万方数据 化 工 学 报 第 67 卷 4738 0cos( 0.5) sin( 0.5) sin( 0.5) cos( 0.5)uvy zwy= = (7) 流场周期 T 为 2，时间步长为 0.05。 图 12 为圆球在三维剪切流场中八叉树自适应网格模型,图 13 和图 14 分别是网格最大层数为 8层和 9 层时运动时间分别为 0、 T/4、 T/2、 3T/4 和 T时刻的数值模拟结果。图 15 为刘长春22采用Youngs-VOF 方

35、法模拟圆在二维剪切流场中的运动的结果， 图 16 为张明远23采用均匀网格的 Level Set方法追踪圆球在三维剪切流场中的运动界面的结果。 通过图 13 和图 14 的数值模拟结果可以得到：圆球在剪切流场中经过一个时间周期后又恢复到原来的形状，但是由于 Level Set 方法的质量守恒性差，圆球发生一些变形。比较图 13 和图 14，增加网格最大层数可以改善 Level Set 方法在剪切流场中的数值模拟结果， 提高 Level Set 方法的质量守恒性以及界面的精度。比较图 14、图 15 和图 16 可以发现， 图 15 中的界面追踪结果出现了破碎现象， 其结果很不理想，图 16 中

36、的均匀网格 Level Set 方法模拟结果的质量守恒性差，本研究提出的方法则可以较好地保持质量守恒性且界面形状保持较好。 6.3 方腔铸件充型过程数值模拟算例 为了进一步证明基于八叉树自适应网格技术的 Level Set 方法追踪界面运动的可行性， 本算例模拟方腔铸件充型流动过程中金属液体界面的运动过程，该过程中金属液体通过浇注系统流入方腔铸件并进行充型，其计算区域为 0,1 0,1 0,1，网格尺寸为 32 32 32。 数值模拟结果如图 17 所示，图中红色区域表示铸型，蓝色区域表示充型金属液体，紫色区域表示铸型中的空气即待充型区域。 由图 17 可以看出采用八叉树技术的 Level S

37、et 方法可以很好地模拟铸件充型过程中液体界面的运动位置，能够较为准确地描述铸造过程中金属液体的充型流动。 7 结 论 研究了基于八叉树网格技术的 Level Set 方法追 图 12 圆球在剪切流场中运动的八叉树自适应网格模型 Fig.12 Octree adaptive mesh model of sphere moving in shear flow filed 图 13 网格最大层数为 8 层时圆球的数值模拟结果 Fig.13 Numerical simulation results of sphere while maximum octree mesh layer is 8 图 14

38、 网格最大层数为 9 层时圆球的数值模拟结果 Fig.14 Numerical simulation results of sphere while maximum octree mesh layer is 9 万方数据第 11 期 殷亚军等：基于八叉树自适应网格技术的 Level Set 运动界面追踪方法 4739图 16 基于均匀网格 Level Set 方法追踪剪切流场 中的运动界面23 Fig.16 Level Set method based on uniform mesh captures moving interface in shear flow field23踪运动界面，引入

39、八叉树技术可以改善 Level Set方法模拟运动界面的一些不足，主要的研究工作总结如下。 （ 1）提出了适用于 Level Set 方法追踪界面的八叉树网格自适应模型（八叉树网格技术的细化与粗化模型） ， 通过对局部运动网格的细化与粗化可以减少数值模拟的计算网格与内存，提高计算的效率与精度。 （ 2）建立了一套基于八叉树自适应网格的Level Set 方法追踪运动界面求解方法，并提出了基于八叉树自适应网格 Level Set 函数值离散模型以及求解八叉树相邻网格距离场函数值的插值方法。 （ 3）提出了八叉树网格距离场重新初始化方法，通过快速扫描八叉树网格并判断网格是否为界面网格更新非界面网格

40、单元的 Level Set 函数值， 重新初始化八叉树网格中的距离场函数值。 （ 4）两个运动界面算例模拟结果表明：将八叉树技术引入 Level Set 方法可以成功地追踪运动界面，局部细化八叉树网格以及增加网格的最大层数可以提高 Level Set 方法追踪界面的精度以及质量守恒性，减少 Level Set 方法产生的数值耗散。并用该方法模拟铸件充型过程的金属液体流动，结果表明该方法可以较为准确地模拟铸件充型过程。 图 15 Youngs-VOF 方法追踪剪切流场中的运动界面22 Fig.15 Youngs-VOF method captures moving interface in sh

41、ear flow field22万方数据 化 工 学 报 第 67 卷 4740 符 号 说 明 S(0) 距离函数符号 T 运动周期， s t 时间， s t 时间步长， s u, v,w 分别为直角坐标 3 个方向的速度， s1x,y,z 三维直角坐标系的 3 个方向 x,y,z 分别为 3 个方向网格单元长度 Level Set 函数 References 1 OSHER S, SETHIAN J A. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formul

42、ationsJ. Journal of Computational Physics, 1988, 79(1): 12-49. 2 OSHER S, FEDKIW R P. Level set methods: an overview and some recent resultsJ. Journal of Computational Physics, 2001, 169(2): 463-502. 3 SUSSMAN M, SMEREKA P, OSHER S. A level set approach for computing solutions to incompressible two-

43、phase flowJ. Journal of Computational Physics, 1994, 114(1): 146-159. 4 WU K T, HAO L, WANG C, et al. Level set interface treatment and its application in Euler methodJ. Science China Physics Mechanics & Astronomy, 2010, 53(2):227-236. 5 CHEN S, MERRIMAN B, OSHER S, et al. A simple level set method

44、for solving Stefan problemsJ. Journal of Computational Physics, 1997, 135(1): 8-29. 6 SETHIAN J A. A fast marching level set method for monotonically advancing frontsJ. Proceedings of the National Academy of Sciences, 1996, 93(4): 1591-1595. 7 SETHIAN J A. Evolution, implementation, and application

45、of level set and fast marching methods for advancing frontsJ. Journal of Computational Physics, 2001, 169(2): 503-555. 8 ADALSTEINSSON D, SETHIAN J A. A fast level set method for propagating interfacesJ. Journal of Computational Physics, 1995, 118(2): 269-277. 9 PENG D, MERRIMAN B, OSHER S, et al. A

46、 PDE-based fast local level set methodJ. Journal of Computational Physics, 1999, 155(2): 410-438. 10 STRAIN J. Tree methods for moving interfacesJ. Journal of Computational Physics, 1999, 151(2): 616-648. 11 SOCHNIKOV V, EFRIMA S. Level set calculations of the evolution of boundaries on a dynamically adaptive gridJ. Internat