2021年《计算传热学的近代进展》-第四章.pdf

计算传热学的近代进展 第四章 捕捉界面的VOSET方法 主讲 陶文铨 西安交通大学能源与动力工程学院 热流科学与工程教育部重点实验室 2021年3月15日, 西安 1/125 第四章 捕捉界面的VOSET方法 4.1 VOF与Level Set方法简介 4.2 捕捉界面的VOSET方法 4.3 VOSET方法求解有相变的气液两相流问题 4.4 VOSET向三维的推广 2/125 4.1 VOF与Level Set方法简介 4.1.1 引言 4.1.2 VOF方法简介 4.1.3 Level Set 方法简介 3/125 4.1 VOF 与 LEVEL SET 方法简介 4.1.1 引言 两相流是一类非常普遍的自然现象，广泛引用于多个工程 领域。例如：液滴/气泡流动，凝结，沸腾，晶体生长。两相流 中的运动界面的存在对直接数值模拟提出了新的要求。 4/125 气液两相流数值模拟方法的分类 气液两相流的模拟由于存在两相而使数值方法的复杂性 大大增加，至今连如何分类也还没有统一的意见；在连续介 质的范围内就气液两相流的数值模拟而言，可以大别为以 下两大类： 1. 精细扑捉相界面的方法：随着流体的运动相界面的位置移动 能够精细地予以描述， 如VOF； 2. 非精细扑捉相界面的方法：求解流场及容积含气率等相分 布的变化，如均相模型，双流体模型等；即使分辨不同的相， 也只是假定分散相是均匀的颗粒，并不精细刻画其界面的变化， 如颗粒轨道模型，等。 5/125 目前常用的几种的相界面捕捉方法： (1) Volume of fluid (VOF)，流体体积方法 (2) Level Set (LS)，水平集方法 (3) Front tracking，前沿跟踪方法 Unverdi S.O., Tryggvason G. , A front tracking method for voscous incompressible multifluid flow. J. Compt. Physics, 1992, 100:25-34 (4) Phase-field，相场法 Jacqmin D., A Calculation of two-phase Navier-Stokes flows using phase-field modeling. J. Compt. Physics, 1999, 155:96-127 (5) VOSET，流体体积和水平集耦合方法 6/125 4.1.2 VOF方法简介 1981年 由 Hirt 和 Nichols 提出。 1.流体体积函数 1) 定义 VOF方法用流体体积函数 (volume function) C 来标识两 相的分布。一个计算网格的流体体积函数指主相(reference phase)流体在整个网格中所占的体积比例。 根据流体体积函数的定义： C=0 表示网格不含主相； C = 1 表示网格完全被主相占据； 00 有以下 四种情况 y x 通过对界面方向进行翻转，可均化归为 nx > 0, ny > 0 的问题， 可能的界面形状只有以下四种类型。 21/125 相界面 (a) 类型1 (b ) 类型2 (c) 类型3 (d) 类型4 对于每一类型的界面，都可以根据以下条件用几何的方 法确定其位置： (1) 相界面垂直于法线 (nx, ny); (2) 阴影部分面积比等于流体体积函数； 据此，可求出线段的两个端点，也就完成了界面重构。 22/125 界面法线矢量的两个分量可用相邻九点的C 值确 定: nix, j  (Ci 1, j 1  2Ci 1, j  Ci 1, j 1  Ci 1, j 1  2Ci 1, j  Ci 1, j 1 ) / 8 x niy, j  (Ci 1, j 1  2Ci , j 1  Ci 1, j 1  Ci 1, j 1  2Ci , j 1  Ci 1, j 1 ) / 8 y       i, j        23/98 2) 界面推进 考虑网格(i, j)的右边界： Fi+1/2 如果ui+1/2 > 0，红色虚线部分流体从右边界流过界面； 如果ui+1/2 < 0，在网格(i+1, j)中计算界面 i+1/2上流入的流量。 用同样的方法计算网格上、下、左边界主相流体的流出量： Cit,jt  Cit, j  Fi 1/ 2, j  Fi 1/ 2, j  Fi , j 1/ 2  Fi , j 1/ 2  V 24/125 这是一种方向分裂 (split)的 推进算法，其不足之处在于一 部分体积被重复计算。 虚线所示面积被重复计算。 G. Tryggvason, R. Scardovelli, S. Zaleski, Direct Numerical Simulations of Gas-Liquid Multiphase Flows, Cambridge University Press, New York, 2011 25/125 目前已发展出基于许多PLIC的非分裂(unsplit)推进算法， 例如： J. López, J. Hernández, P. Gómez, F. Faura, A volume of fluid method based on multidimensional advection and spline interface reconstruction, J. Comput. Phys, 195, 718-742, 2004 J. Hernández, J. López, P. Gómez, F. Faura, A new volume of fluid method in three dimensions. Part I: Multidimensional advection method with facematched flux polyhedra, Int. J. Numer. Meths. Fluids, 58, 897-921,2008 下列论文的作者也公布了他们的计算程序: VOFTools, a package of FORTRAN subroutines with analytical and geometrical tools for 2D/3D VOF methods in general grids http://www.dimf.upct.es/personal/lrj/voftools.html 26/125 4.1.3 Level Set 方法简介 1.符号距离函数 Level Set方法用一个连续函数 ϕ (Level set函数)的零等值 面表示相界面。用其正负表示位于哪一相。 ϕ<0 表示位于主相中 ϕ>0 表示位于另一相中 ϕ=0 表示位于相界面上 相界面 符号距离函数 (Signed Distance Function)就是一种Level Set函数。对于某一点的符号距离函数，其绝对值等于这个点 到相界面的最近距离，符号则取决于这个点位于哪一相。 S. Osher, J.A. Sethian. Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations. J. Comput. Phys, 79 (1988)12-49 27/125 用两个例子来说明符号距离函数 (正负值随定义而异）。 (1) 一维区域，主相占据 [-1, 1] 区间; (2) 二维区域，主相占据以原点为圆心，半径为 1 的圆内部。则 符号距离函数为一组同心圆。 f 在主相中点 f =1 f =1− x outside y f  1  x 2 y 2 在相界面上 f =0 inside inside -1 outside 1 x -0.6 -0.4 -0.2 0 0.2 0.4 inside outside 主相 表示符号距离函数的直线(x大于0） 28/125 x 符号距离函数的特点： （1）在相界面附近连续光滑； （2） f  1 距离对距离求导数！ 2.符号距离函数在两相流模拟中的应用 1) 计算界面方向与曲率 符号距离函数可以用于计算界面方向与曲率： n  f f    f 由于符号距离函数在相界面附近连续光滑的特点，用其计算 界面曲率与方向与曲率比用流体体积函数要准确的多。 29/125 2) 计算光顺化的Heaviside函数 用符号距离函数可以定义光顺化的Heaviside函数： 0 1 f 1  H  (f )   [1   sin(f /  )]   2 1 (f   ) ( | f |  ) (f   ) ε 表示界面处光顺化的宽度，一般取 1.5Δ (Δ为网格尺寸) 3) 设置两相流的物性分布 利用光顺化的Heaviside函数H和两相的物性可以设置对计 算区域的物性分布：    g (1  H  (f ))  l H  (f )    g (1  H  (f ))  l H  (f ) 以这种物性来算流场以考虑存在两相的事实。主相是液相 还是气相取决于模拟什么问题。 30/125 4) 计算表面张力 气液两相流的模拟还需要考虑相界面上的表面张力。CSF (Continuum surface force) 模型可以把表面张力转换为相界面附 近控制容积中的体积力，并用符号距离函数计算： Fsv   (f )H (f ) σ 为表面张力系数；  (f ) 为界面曲率，通过符号距离函数计算； H (f ) 为界面H函数的梯度。 符号距离函数直接表示了节点到相界面的距离。因此，符 号距离函数为两相流的直接模拟提供了很多方便，是首选的 Level Set函数。 J. U. Brackbill, D. B. Kothe, C. Zemach. A continuum method for modeling surface tension. J Comput Phys. 100(1992) 335-354 31/125 3. 用Level Set函数追踪相界面 Level Set函数的推进也可表示为纯对流方程： f  u f  0 t 如果用上述纯对流方程来推进Level Set函数时，一般设置 了初始符号距离函数，但是在推进过程中距离函数的特征会改 变。因此，在每一个时层完成推进后，一般还需做进一步数学 处理，重新初始化，以将推进后的Level Set 函数具有符号距离 函数的特点。 符号距离函数应该满足： f  1 S. Osher, R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, Springer, 2002 32/125 重新初始化方程的求解目标是使Level Set函数在零等值面不 改变的情况下满足 f  1 ，恢复为符号距离函数。 Level Set 推进方程和重新初始化方程的求解一般采用高阶的 ENO(Essential Non-Oscillatory) 格式。重新初始化的方程也是 一个纯对流方程！ Level Set 方法的一个优势是简单方便，从二维到三维的推 广也十分容易。因此在更复杂的两相流 (例如相变) 模拟应用广 泛。 Level Set方法的缺点是：除了计算上比较费时外，它仅追 踪 ϕ = 0 的等值面，难以保证相的总体质量守恒性。 A Harten , B Engquist, S. Osher, et al., Uniformly high order essentially nonoscillatory schemes , J. Comput Physics, 1987, 71:231-303 33/125 4.2 捕捉界面的VOSET方法 4.2.1 VOSET方法的基本思想 4.2.2 VOSET方法中距离函数的计算 D.L. Sun, W.Q. Tao, A coupled volume-offluid and level set method (VOSET) for computing incompressible two-phase flows, Int. J. Heat Mass Transfer, 53, 645-655, 2010 推荐阅读(13) 孙东亮 34/125 4.2 捕捉界面的VOSET 方法 4.2.1 VOSET方法的基本思想 VOF方法和Level Set方法各有优缺点，而且正好互补。 如果能够同时拥有VOF中的流体体积函数和Level Set中 的符号距离函数，就可以让它们发挥各自的优势，执行各自 “擅长”的任务。但不宜简单的叠加，否则计算工作量会显 著增加，文献中曾提出过这样的方法。 35/125 流体体积函数负责界面重构与界面推进，保证总体质量守 恒性。 符号距离函数负责计算界面方向与曲率，计算表面张力，设 置物性。 VOSET用流体体积函数追踪相界面，同时根据流体体积函 数，用几何方法生成一个与之相对应的符号距离函数。这样 VOSET就能同时具备VOF和Level Set方法的功能。 VOSET不是第一个将VOF及Level set 结合起来的方法， 此前有Sussman （2000）提出的CLSVOF方法及其以后的各种 改进；但是在该系列方法中同时求解VOF及Level set 函数两个 微分方程，使得计算工作量及复杂程度显著增加。VOSET方 法中只求解流体体积函数，但用几何方法生成距离函数。 36/125 4.2.2 VOSET方法中距离函数的计算 1.迭代计算的思路 1） 为什么要迭代 准确的界面方向在界面捕捉中有重要意义，也与符号距离 函数有密切的关系，先用已知的流体体积函数初步计算的界面 方向需要改进： nix, j  (Ci 1, j 1  2Ci 1, j  Ci 1, j 1  Ci 1, j 1  2Ci 1, j  Ci 1, j 1 ) / 8 x niy, j  (Ci 1, j 1  2Ci , j 1  Ci 1, j 1  Ci 1, j 1  2Ci , j 1  Ci 1, j 1 ) / 8 y (1) 知道了初步的界面方向，就可以用PLIC方法重构相界面； (2) 获得重构界面后，再用下面介绍的几何方法计算符号距离 函数，而通过距离函数又可以求出更准确的界面方向。 PLIC界面重构和计算符号距离函数的准确性都要求有准确 的界面方向。因此需要用迭代的方法计算符号距离函数。 37/125 2） 如何迭代 在没有符号距离函数的时候，首先用流体体积函数计算一 个初始的、粗略的界面方向； 在此基础上，重复以下过程： (1) 用PLIC重构相界面：知道 nx , n y 及 C 即可重构； (2) 根据重构界面，用几何方法计算符号距离函数； (3) 用符号距离函数计算界面方向： n  f    f f 这些过程一般需要迭代三次，就可以得到准确的界面方向 和符号距离函数。 前面已介绍过第 (1)步 和 第(3)步，下面详细介绍第 (2)步： 计算符号距离函数的几何方法。 38/125 2. 计算符号距离函数的几何方法 假定已用PLIC方法获得重构界面。 步骤 1：在整个计算区域内为符号距离函数赋初始值 相界面 M -M M  M fi , j    M 0 if Ci , j  0.5 if Ci , j  0.5 M为计算域内两个方向的最大长度，这样可以保证界 面附近的符号距离函数都包含在设定的初值范围之内。 39/125 步骤 2：标记界面附近的计算网格 标记包含相界面一定距离范围（例如：三个网格长度） 之内的网格。这一步又可显著减少计算量。因为实际计算 中只需要相界面附近的单元的节点计算符号距离函数。 相界面 标示区域 40/125 步骤 3：计算被标记的单元中心与相界面的距离 假设A (i, j)是一个被标记的计算单元，在A的附近搜寻带有 相界面的网格。设此时其附近的相界面已重构为一个线段(BC)。 于是可利用一个三角形(ΔABC)计算给定点到这个线段的最短距 离,无非三种情况：   1  90  1  90 , 2  90 1  90 , 2  90 41/125 最短距离 相界面 i, j 从节点(i, j ）到附近7×7网格 内的相界面的所有最小距离 在以 (i,j)为中心的7X7个网格内进行界面搜索；在完成了 节点(i, j)附近相界面的搜寻、并求出它与这些线段之间的最短 距离后，从这些最短距离中选出最小值。这个值就是点(i, j)到 相界面的最短距离d，也就是该节点的符号距离函数的绝对值。 42/125 步骤4：设置距离函数的符号 一个点的符号距离函数的正负取决于这个点位于哪一 相中。 符号可以用流体体积函数来判定。如果一个网格主相 流体占据的体积超过一半，就可以认为这个网格的中心位 于主相内，反之则位于另一相中。因此： d  fi , j  0  d  if Ci , j  0.5 if Ci , j  0.5 if Ci , j  0.5 43/125 3.符号距离函数计算示例 Frame 001  08 May 2007  Frame 001  08 May 2007  以下是迭代几何方法得到的不同形状相界面附近的符号距 离函数。 符号距离函数 对于无相变 的两相流问题， 如空气泡与水的 流动，获得了符 相界面 (a) (b) 号近距离函数后 完成了界面重构， X X 可以进行下一时 间步长的的流场 求解及界面的推 进。 (c) (d) 0.04 0.03 0.03 Y Y 0.04 0.02 0.02 Frame 001  08 May 2007  Frame 001  08 May 2007  0.01 0.01 0.04 0.01 0.02 0.03 0.03 0.04 0.04 0.01 0.05 0.02 0.03 0.04 0.05 0.04 0.05 0.01 0.02 0.03 0.04 0.05 Y Y 0.03 0.02 0.02 0.01 0.01 0.01 0.02 0.03 44/125 4.3 采用VOSET方法求解气液两相流问题 4.3.1 无相变过程的控制方程 4.3.2 无相变过程的求解步骤 4.3.3 有相变过程的控制方程 4.3.4 有相变过程能量方程的求解 4.3.5 数值计算实例 45/125 4.3 采用VOSET方法求解气液两相流问题 4.3.1 无相变过程的控制方程 流体体积函数方程 连续性方程 C     u C   0 t  u  0 u T    u  u   p      u   u   g  H   动量方程   t 表面张力 重力 每个控制容积的密度、粘性、表面张力均通过Level Set函 数计算，把两相流当做变物性的单相流计算，确定流场。 Level Set函数用几何方法根据流体体积函数直接求得，因 46/125 此不需要数值求解其函数的偏微分方程，大大节省计算资源。 4.3.2 无相变过程的求解步骤 给定初始时刻流体体积函数C0，初始时刻速度u0,v 0 ； 在第 n 个时层的推进步骤如下： (1) 用 PLIC 方法，根据流体体积函数Cn重构界面，并用几何方 法计算相应的Level Set函数； (2) 用基于Level Set函数的Heaviside函数计算物性分布和表面张 力； (3) 用某种流场求解算法（SIMPLER/IDEAL/Projection算法） 求解下一时层的速度场 un+1, vn+1和压力场 pn+1 ； (4) 据计算得到的速度及界面计算下一时层各控制容积的流体体 积函数Cn+1 。 47/125 4.3.3 有相变过程的控制方程(液体沸腾，C蒸汽体积分数) C m (1) 流体体积函数方程:    (Cu )  t v m 为液体蒸发率  1 1 (2) 连续方程：  u = m      v l  表面张 u (3) 动量方程：  u (u )  力作用 t 1 {p    (f )[(u )  (u ) T ]   (f ) g + (f )H }  (f ) (4) 能量方程 ： T  u  T   l  2T 液体区域 t 分区求解 T  u T   v 2T 汽体区域 t 48/125 4.3.4 有相变过程能量方程的求解及液体蒸发率的确定 利用相界面上温度为饱和温度的条件，两个区域各自求解。 1.边界温度确定：对包含界面的单元 A, B,节点A,B 就是边界 节点。A,B点的温度需要通过插值获得:通过该两点做垂直于界 面的线段至另一个单元，得A’及B’： fA A点的温度：  TA'  Tintf fA  d TA  Tintf d 为A’与A的距离, f 为符号距离函数 fB  B 点的温度： TB'  Tintf fB  d TB  Tintf A’及B’本身的温度用插值获得；本次迭代计算获得的A 及 B点的温度，作为下一次迭代计算两个区域已知的边界温度。 49/125  2.液体蒸发率 m 计算： 4 取决于界面的热流密度。采 用探测点法来确定。 在相界面上找一点 I， 沿界面法线向两侧等间距 2 确定点1,2,3,4；利用该四 点周围已知温度插值得出四 点的温度，然后按照Fourier 定律 线性插值: T1  Ti T  n l h 二次插值: T2  4T1  3Ti T  n l 2 h  3 1 T T q  l  v n l n v Ti  T3 T  n v h Ti  4T3 +3T4 T  n v 2 h 50/125 4.3.5 数值计算实例 1.用不同的方法计算圆的曲率 给定一个半径为 1 的静止的圆，分别用 VOF, Level Set, VOSET 在不同网格尺寸下计算圆的曲率。同时VOSET在迭代 求解符号距离函数中选择了不同的迭代次数：N = 1, 2, 3。比较 这些方法计算的曲率误差的大小。 圆的曲率与半径满足：  exact R  1 M L2误差的定义： L2   ( R   i 1 M exact R) 2 其中M为圆形相 界面涉及到的计 算单元数 κ 为计算得到的曲率。 51/125 (1) VOSET 计算曲率的准确性远优于 VOF； (2) 迭代次数增加，VOSET计算曲率的准确性会提高；迭 代三次偏差已经下降到5%或以下，已经可以接受。 52/125 Frame 001  14 Jul 2000  Frame 001  14 Jul 2000  (b) (a) VOF Level Set (c) VOSET 53/125 2. 椭圆形气泡的震荡 无重力的条件下，在静止的液体中放置一个初始的椭 圆气泡。 椭圆界面上的曲率并不均匀，因此表面张力不均匀。 不均匀的表面张力会导致气泡的震荡。 液体的粘性使震荡逐渐减弱，最终气泡会变成圆形 并保持静止。 54/125 Frame 001  02 Jan 2008  Frame 001  03 Jan 2008  t = 0, 压力 Frame 001  02 Jan 2008  t = 3, 速度 Frame 001  02 Jan 2008  Frame 001  02 Jan 2008  Frame 001  02 Jan 2008  Frame 001  02 Jan 2008  Frame 001  02 Jan 2008  Frame 001  02 Jan 2008  t = 5, 相界面 VOF Level Set VOSET VOSET继承了Level Set计算表面张力的准确性，优于VOF。 55/125 VOF/ VOSET Level Set 计算时刻主相流体的质量 质量比率（Ratio）= 初始时刻主相流体的质量 VOSET 继承了 VOF 的总体质量守恒性，优于 Level Set。 56/125 3.单个气泡的上升 在静止的液体放置一个气泡，气泡在浮力的作用下上 升。 气泡在上升的过程中受到浮力、粘性力和表面张力的 共同作用，最终形成的稳定的气泡形状由无量纲数 Morton number和 Eotovs number 决定。 M  g / l 4 l 3 Eo  gd ( l   g ) /  2 e 表征粘性力与表面张力的对 比关系 表征浮力与表面张力的对比 关系 57/125 Case 1: Eo = 1.0, M = 0.001，会形成球圆形气泡 58/125 Case 2: Eo = 10.0, M = 0.1，会形成椭圆形气泡 59/125 Case 3: Eo = 100.0, M = 1000, 会形成球帽形气泡 60/125 4. 浅液层沸腾过程 液面 Critical liquid height hc 气泡 加热壁面 Normal BHT 浅液层中的沸腾 When the height of liquid over a heating plate is less than a certain value boiling heat transfer can be significantly enhanced. 61/125 Shallow liquid layer 4 5x10 4 热流密度/ W·m-2 4x10 4 3x10 4 2x10 深水条件 浅水条件 4 1x10 Deep liquid layer 0 0.0 0.2 0.4 0.6 时间/ s 0.8 1.0 62/125 5.膜态沸腾模拟 Experiment by Reimann and Grigull (1975) 郭东之 3 Nu 2 Klimenko correlation 1 低过热度 高过热度 0 数值模拟结果 Klimenko t(s) 0 1 2 3 4 63/125 4.4 VOSET向三维的推广 4.4.1 VOSET从二维推广到三维需要解决的问题 4.4.2 VOSET从二维到三维的简单推广---沿用二维 问题的方法 4.4.3 VOSET从二维到三维的智慧推广---利用体积 分数与界面位置之间的内在联系 4.4.4 三维VOSET应用例子 4.4.5 国内外学术界的评价及推广 64/125 4.4 VOSET向三维的推广 4.4.1 VOSET从二维推广到三维需要解决的问题 VOSET中重构相界面和计算符号距离函数都采用了几何方 法，因此在二维中实施VOSET需要解决的是一个平面几何问题。 在三维坐标系中实施VOSET需要解决的是一个立体几何问 题。因此，VOSET向三维推广的关键问题在于如何解决更复杂 的几何问题。 从二维到三维要解决三个问题： （1）如何推进三维体积函数； （2）如何计算符号距离函数； （3）如何重构三维相界面？ 65/125 二维 三维 相界面的推进比较简单，三个方向分别做推进，不是讨论 重点； (1) 计算点到相界面的最短距离 在二维中，通过解一个三角形就可以计算出一个点到 相界的最短距离，同样的问题在三维中应该如何计算？ (2) 相界面的描述与重构 二维问题中用线段两个端点的坐标就可以描述重构相 界面，三维中的重构相界面应该如何描述？ 66/125 4.4.2 VOSET从二维到三维的简单推广---沿用二维问题 的方法 一种简单的推广方法就是完全沿用二维的方法，只是增加 了一个维度。 1. 流体体积函数的推进 (同二维情形) 三维空间含相界面的一个控制容积在一个时间步长后体积 函数的变化为： (Cin, j ,1k  Cin, j ,k )xyz  Fi n1/ 2, j ,k  Fi n1/ 2, j ,k  Fi ,nj 1/ 2,k  Fi ,nj 1/ 2,k  Fi ,nj ,k 1/ 2  Fi ,nj ,k 1/ 2 u i+1/2 Δt 对于存在相界面的控制容积： Fi 1/2  0 ui 1/2 t  (1  C )x ui 1/2 t  (1  C )x ui 1/2 t  (1  C )x i-1 Ci Δx ui+1/2 i i+1 (1  C )x 67/125 2. 符号距离函数的迭代计算 (基本同二维) 在节点(i,j,k) 的三个坐标方向取7x7x7个控制容积， 寻找该点到位于多个控制容积内相界面距离，并从中选出 最短距离作为该点的符号函数。 F点在BCDE相界面面内 F点在BCDE相界面面外 68/125 为确定垂足F的坐标，首先需要计算垂线AF的距离： | AB  n | AF  |n| 如果垂足F位于界面BCDE内， 垂线AF就是节点A到界面BCDE的 最小距离；如果垂足F位于界面BCDE外，首先需要分别 计算节点A到界面各个边（线段BC，CD，DE，EB）的 最小距离（计算方法如二维问题），然后比较这些最小 距离，得出节点A到界面BCDE的最小距离d。 d  fi , j  0  d if Ci , j  0.5 if Ci , j  0.5 if Ci , j  0.5 69/125 3. PLIC界面重构 (二维简单推广) 当 (i,j,k) 控制容积内是： 0