（应用数学专业论文）vof运动界面重构的流体体积分数保持法.pdf

摘要 在矩形网格上，基于v o f 方法的基本思想，针对运动界面的重构问题，本文提出 了流体体积分数保持法，并对该方法进行了一些研究。 该方法是在单个网格内用斜线段近似运动界面，并要求相邻网格内的斜线段在公 共边上的交点重合，通过保持网格内的流体体积分数不变，从而在界面网格上建立非 线性代数方程组，通过求解该非线性代数方程组，确定运动界面的边界与网格线的交 点位置，从而重构出运动界面的形状和位置。 本文对所建立的非线性代数方程组给出了最速下降法和直接法两种求解方法。最 速下降法得到的最优解，使运动界面的重构不但能够达到较高精度，而且可应用于大 型网格计算。在对流体体积分数保持法得到的非线性代数方程组的研究基础上，得到 方程组统一的表达式，给出了一种直接求解方法，该求解方法在解存在的情况下，所得 到的解精确度很高，而且在处理大型网格计算时运算速度很快。并在直接求解方法研 究的基础上，进一步讨论了流体体积分数保持法解的适定性问题。通过对流体体积分数 保持法的改进，将其推广应用到含多个闭合型的界面重构和非闭合型的界面重构问题 。 关键词：运动界面重构，v o f 方法，流体体积分数保持法，界面网格 a b s t r a c t i nr e c t a n g u l a rm e s h e s ，t h ek e e p i n gv o l u m ef r a c t i o no ff l u i d si n v a r i a b l em e t h o di s p r e s c n t a df o rr e c o n s 打u c f i o no f m o v i n g - i n t e r f a c e s , b a s e do nt h ev o fm e t h o d , a n di ss t u d i e d a sf o f l o w s 。 t h i sm e t h o d 咖a p p r o x i m a t em o v i n gi n t e r f a c e sb yad i a g o n a ll i n ei nas i n g l eg r i 正a n d p o i n t so ft h ei n t e r s e c t i o no fd i a g o n a ll i n e sa n d t h ep u b l i cs i d e sb e t w e e nt h ea d j a c e n t 舒d s c o i n c i d e a c c o r d i n gt ok e e pt h ev o l u m ef r a c t i o ni n v a r i a b l ei ne a c hm e s h , t h en o n l i n e a r a l g e b r a i ce q u a t i o n s a l eb u i l ti ns u r f a c em e s h e s b ys o l v i n gt h ee q u a t i o n s ，t h ep o s i t i o no f t h e i n t e r s e c t i o no f m o v i n g - i n t e r f a c e sa n dt h eg r i dl i n e sa r cd e t e r m i n e d , 8 0t h es h a p ea n dl o c a t i o n o f m o v i n gi n t e r f a c e sa r es i m u l a t e d i tp r o v i d e st w om e t h o d s ，t h es t e e p e s td e s c e n ta l g o r i t h ma n dt h ed i r e c tm e t h o d ，t os o l v et h e e s t a b l i s h e de q u a t i o n s t h es t e e p e s td e s c e n ta l g o r i t h mc a l la c c u r a t e l ys i m u l a t em o v i n g i n t e r f a c e sa n db ea p p l i e dt ot h el a r g eg r i d sc o m p u t i n g b a s e do nt h er e s e a r c ho ft h e e s t a b l i s h e dn o n l i n e a ra l g e b r a i ce q u a t i o n si nt h ek e e p i n gv o l u m ef r a c t i o no ff l u i d si n v a r i a b l e m e t h o d , t h eu n i f i e de x p r e s s i o no f e q u a t i o n si si n f e r r e da n dt h ed i r e c tm e t h o di sp r o p o s e d o n t h ec o n d i t i o no fs o l u t i o no ft h ed i r e c tm e t h o de x i s t e d ，t h es o l u t i o ni sv e r ya c c u r a t e ，a n dt h e c o m p u t i n gs p e e do f t h e d i r e c tm e t h o di sv e r yf a s tt oh a n d l el a r g e g r i d sc o m p u t i n g a n di ti s f u r t h e rt od i s c u s sw e l l - p o s e d n e s so ft h es o l u t i o nb yt h ed i r e c tm e t h o d t oi m p r o v et h e k e e p i n gv o l u m ef r a c t i o no ff l u i d si n v a r i a b l em e t h o d ，i ti sp r o m o t e dt os o l v et h es i m u l a t i o n o f m a n yd o s e di n t e r f a c e sa n do p e i li n t e r f a c e s k e y w o r d s ：r e c o n s t r u c t i o n so fm o v i n gi n t e r f a c e s ；v o l u m eo ff l u i d sm e t h o d ；k e e p i n gv o l u m e f r a c t i o no f f l u i d sm e t h o d ；s u l f a t em a s h e s 学位论文独创性声明： 本人所呈交的学位论文是我个人在导师指导下进行的研究工作及取 得的研究成果。尽我所知，除了文中特别加以标注和致谢的地方外，论文 中不包含其他人已经发表或撰写过的研究成果。与我一同工作的同事对本 研究所做的任何贡献均已在论文中作了明确的说明并表示了谢意。如不 实，本人负全部责任。 一4 一t ：- 一 论文作者( 签名) ：t 登喙垦圣叫年，月乙7 日 学位论文使用授权说明 河海大学、中国科学技术信息研究所、国家图书馆、中国学术期刊( 光 盘版) 电子杂志社有权保留本人所送交学位论文的复印件或电子文档，可 以采用影印、缩印或其他复制手段保存论文。本人电子文档的内容和纸质 论文的内容相一致。除在保密期内的保密论文外，允许论文被查阅和借阅。 论文全部或部分内容的公布( 包括刊登) 授权河海大学研究生院办理。 论文作者( 签名) - 蓝军斗硼年 - 卜月习 日 第一章绪论 第一章绪论 1 1 问题的提出及研究意义 运动界面的流体流动现象广泛存在于机械、能源、电子、动力、水利和建筑等多 种工程领域中。如海面的波动、海浪翻滚，漫滩、涌浪，溢流坝【1 l ，金属的铸造嘲， 气泡的运动1 3 1 ，石油的加工、运输，空间飞行器的设计，化学反应及扩散过程，电子 器件的制作，等等。对于这人们最为关注的是运动界面的位置等特征。因而，模拟和 追踪它们的轨迹、特征和发展，是了解这类现象的关键。在理论上这种界面流动也是 流体力学中一个具有挑战性的研究工作，因为在许多间断问题的数值计算和数值模拟 中，对它们的数值模拟要比强间断的激波问题更加困难。许多的运动界面恰恰多为这 类界面。不仅如此，由于它是运动的，常常体现其活动边界的某种不适定性，更增加 了研究的难度，因此长期以来受到各国学者的重视。 运动界面是指流体在这些界面上的一个和多个物理量发生间断或不连续的变化。 它的特殊性体现在两个方面，一方面它是流场的边界，是流场求解的必要条件；另一 方面它的位置是未知的，如何去确定它的位置，本身就作为问题解的一个重要部分。 事实上，运动界面问题的提出和讨论可以追溯到一个多世纪以前。1 8 8 9 年，s t e f e n 最 早研究了所谓“活动边界”问题的一个典型实例一冰水融化问题，因此，该类问题通 常又称为s t e f e n 问题。在生活中经常遇到的诸如石蜡熔化、水结冰、泄露在海面的石 油的扩散等边界随时间变化的问题，都是这种活动边界问题【4 1 。此外，“自由面问题”， 诸如溃坝，波形在斜坡上爬升、变形、以及破碎，液体绕障碍物产生的波形的扰动【5 1 等， 与我们的生活密切相关。另外还有一类问题，如多相、多组分的气体和液体等的内部 速度、压力、密度等所产生的间断，或者内部产生的各种不稳定现象。而这些间断或 不稳定状态，会依据一定的物理或化学的规律移动变化。 从对自由面问题的数值计算研究，发展到对运动界面追踪的数值模拟研究。这个 问题具有很深刻的现实意义和应用价值。例如，h i r t 等人提出的v o f 方法的运动界面 追踪技术，现在已发展和形成了很多实用软件，如v o f - 3 d ，n a s a - v o f 【6 7 1 ，s o l a v 0 f 嗍等，并且已经成功地应用于半导体工业、喷墨技术等近代高新技术产业中。 从理论上讲，它一般是介质面问题，不同密度液体间的不稳定问题，固化和熔化界面 问题，等等。这些问题是通常所谓的接触间断问题，对它们的数值模拟比较困难。特 河海大学硕士论文 别是要得到精细而锐利的高分辨率图像更加不易，因而对它的数值研究和分析日益受 到关注。它是最现实、最直接地调查，分析和把握许多微观未知过程的途径，同时也 是理论方面，特别是数值方法和计算科学发展的必须。 1 2 运动界面追踪的研究方法 运动界面是一个非常广义的概念，但就具体的实际问题而言，各种运动界面、介 质面和间断面，有着不同的数学表达形式。在当前的追踪方法中，针对不同问题提出 了各种不同的追踪技术，并进行了较为广泛的研究，取得了较丰富的成果。归结起来， 界面追踪技术可以分为三类：拉格朗日法( l a g r a n g i a n ) 、欧拉法( e u l e r i a n ) 和任意拉 格朗日一欧拉法( a r b i t r a r y l a g r a n g i a n - e u l e r i a n ) 。对于拉格朗日法，由于物质点与网格 点不存在相对运动，所以在控制方程中不出现对流项，大大简化求解过程。另外，网 格能自动跟踪流体域的移动边界，可以准确地跟踪运动界面，方便处理边界条件。但 是，当流体大变形时网格会发生畸形，甚至扭曲，所以拉格朗日法仅能应用于小变形 情况；欧拉法在流体运动中网格固定在空间不动，易于处理复杂的运动界面情况，但 必须使用其它方法捕捉运动界面，同时会出现数值的震荡和耗散问题，目前对这些问 题已取得较好的解决方法，其中的v o f ( v o l u m eo f f l u i d s ) 方法以其存储小、应用简单、 追踪的界面锐利性好和精细，而深受广大c f d 工作者的喜爱；任意拉格朗日一欧拉法 吸收前面两种方法的优点，在网格取得合理情况下，既使网格重分计算量大大降低于 拉格朗日法，又可以消除数值震荡和过大的耗散，提高模拟精度。但是取得合理的网 格运动并不容易，特别是当运动界面上出现翻卷、吞并和飞溅等现象。 在解决运动界面问题方面，出现了很多行之有效的方法。最早处理自由面问题的 有效方法，是h a r l o w 和w e l c h 等人提出的格子类( c e l l - t y p e ) 方法1 9 j 们。这类方法在 e u l e r 网格上进行差分离散，在网格中布置若干标记点( m a k e r ) 或者流体质点( p a r t i c l e ) 来标记流体。它是通过这些标记点或质点的追踪，来模拟自由面的大体位置。网格是 e u l e r i a n 型的，固定不动，而标记点却随流体运动，是l a g r a n g i a i l 型的。在格子类方法 中，最著名的p i c t l ”( p a r t i c l e i nc e l l ) 、m a c t l 2 蚓( m a r k 盯a n dc e l l ) 和f l i c t “1 ( f l m d i n c e l l ) 等。这类方法能给出自由面的粗糙近似，但不能给出更加细致的描述，更主要的 是计算时仍需要较大的储存量，尤其是对于三维的情况。 以下介绍几种常见的现代运动界面跟踪的数值方法： ( 一) v o f ( v o l u m eo f f l u i d s ) 方法 在当前的运动界面追踪方法中，v o f ( v o l u m eo f f l u i d s ) 方法是最为重要的方法 之一，它的特点是将运动界面在空间网格内定义成一种流体体积函数( 这种流体体积 函数既适应于单介质问题，也适应于双介质或多介质、多相问题) ，并构造这种流体体 2 第一章绪论 积函数的发展方程，从而随着主场的模拟过程，通过流体运输，精细地确定该运动界 面的位置、形状和变形方向，达到追踪的效果。v o f 方法的发展在以下三个方面【1 5 】： ( 1 ) 基本场模型的改进。引入运动界面表面张力的效应 fv 矿= o 1 警+ v ( 痧) ：琊+ 去v 什专藤+ f c n 忍d p = c p l + o - c o p , ， = c h + ( i - c ) a 著冉( 网= o ( 2 ) 主场求解方法逐步采用了高分辨率的各种数值方法，如r o e “】，o o d u n o v ”， t v d 1 引，e n o 1 ”，w e n o ，f v l 4 等。 ( 3 ) 重构运动界面采用了更为精细的构造，网格内界面迁移计算越来越巧妙。如 y o u n g s 2 ”，g u e y f f i e r ，l i j i e 2 1 1 等的论文。 ( - 3 波前追踪方法( f r o n tt r a c k i n gm e t h o d ) 波前追踪方法，是将间断面、物质界面、自由面等看作活动界面，先计算它们的 位置等特征，再结合边界条件、跳跃条件等，求解物理量的控制方程，追踪活动界面， 并在活动边界处做高精度、高分辨率的处理，以求得逼真的物理图像。 活动界面追踪方法的历史可以追溯较远，早在2 0 世纪4 0 - 5 0 年代，由于战争的激 励，强激波的追踪计算发展的非常迅速，提出了多种优秀的激波装配追踪的方法。 到1 9 8 0 年，c h e n 和g l i m m ”等人正式地将波前追踪方法作为一种方法提出来， 处理界面相交、变形与边界的相互作用，而且避免弱界面的消失和虚假界面的产生， 比较成功地计算了激波、反射，特别是k e l v i n - h e l m h o l t z 的不稳定性。 1 9 8 6 年，m a r s h a l l ，c h a r d e 1 和t e s s i e r e s t “1 分别对运动界面和守恒律型的控制方 程提出了波前追踪方法，前者采用插值外推方法，后者采用r i e m a n 问题方法的g o d u n o v 格式。 1 9 9 5 年，l e v e q u e 和s h y u e 提出了一种波前追踪算法，用于解决一维守恒问题。 这种方法在处理激波、接触间断问题时方便有效。它是基于r o e 近似解将解看作元素 波的叠加，每个元素波以各自的速度传播。计算得到下一时刻的元素波后，再叠加起 来得到物理解。在产生间断的地方增加追踪点，来达到波前追踪的目的。 但是，波前追踪方法的间断仍离不开r i e m a n n 间断解的特征和装配技术，对运动 的物理界面很难达到高效、高逼真的效果。对于高维问题实施难度也比较大。 ( 三) l e v e l s e t 方法 河海大学硕士论文 1 9 8 8 年0 s h c r 2 6 1 等人提出了一种零等值面方法。这种方法在许多复杂的界面追踪 问题中得到成功，引起了广泛的兴趣和研究。以o s h e r 为首的l e v e ls e t 研究群体在这 方面发表了大量有指导意义的文章【2 7 ”。 l e v e ls e t 方法的主要思想是将运动界面定义为一个函数的零等值面( 线) ，即 r ( f ) = x q 1 + q 2 ：伊( 而f ) = o ，其中，q 1 、q 2 分别为两种不同的流体区域。在流体运 动的过程中，然后始终保持它是零等值面，而且在界面附近保持单调。运动界面的初 值定义为： i + d i s t ( x , r ( o ) ) x q 1 妒( 工，0 ) = - d i s t ( x ，r ( o ) ) x q 2( 1 2 2 ) 【0 xr ( o ) d i s t ( x ，r ( o ) ) 表示x 到r ( o ) 的距离。当伊以适当的速度移动时，在任意时刻，只要求出矿 的值，就可以确定运动界面的位置，再由此求解物理量的控制方程。 ( 四) p h a s e - f i e l d 方法 1 9 9 5 年，a n t a n n o v s k i i 2 9 1 和j a c q m i n 【”1 提出了一种所谓的相域方法，对复杂的三 维自由界面问题的相域模型进行了讨论和模拟。1 9 9 9 年，j a c q m i n 提出了d i f f u s e i n t e r f a c em o d e l s 方法【3 l i ，考虑了扩散、界面效应的效果。不同于v o f 的重构、l e v e l s e t 的的重构或是f r o n tt r a c k i n g 的装配，这种方法允许运动界面有一定的扩散。它 对多相、多组分和多维问题的适应性也很有启发性。 1 3 本文工作 运动界面的重构是v o f 方法追踪运动界面的研究课题之一，本文以v o f 方法为基 本理论背景，是对界面重构的深入研究。主要工作： ( 1 ) 介绍了v o f 方法的原理和发展应用，总结了v o f 方法中各种界面重构的方法 和捕捉界面的方法，这些方法的提出丰富了v o f 方法的内容，拓展了v o f 方法的应用 范围，奠定了v o f 方法在现代运动界面追踪的地位。 ( 2 ) 在矩形计算网格上，本文提出了一种的v o f 运动界面重构的算法一流体体积 分数保持法，阐述这种方法的基本思想。并通过优化方法进行计算，结合实例进一步 分析了这种方法的可行性及应用性。 ( 3 ) 应用直接法求解在流体体积分数保持法中得到的非线性代数方程组。通过对 三种不同类型界面网格对应方程的分析，得到了一个统一的方程。在此基础上，运用 直接法把这个方程组的求解问题转化成一个一元至多二次方程的求解问题，把问题简 4 第一章绪论 单化，既而求解。通过算例表明，直接法不仅能够使求解精度很高，而且运算速度很 快，特别是对于大型的网格计算，计算时间将比优化法节省很多。 通过直接法的算法过程分析了一些方程组的解的适定性问题，讨论了几个方面的 情况：首先讨论了无穷多个解的情况，推导出流体体积分数在满足什么条件时方程组 的解是有无穷多个的；接着讨论了不合要求的解、两个解和无解的情况，通过给出的 反例，我们推测出现这些情况的原因是，流体体积分数保持法的前提要求，即相邻两 个界面网格内的界面在公共边上相交于一点，提出的过于严格，这就需要我们在以后 的工作中继续对这种界面重构方法进行探讨，找到更加合理的前提条件。 ( 4 ) 通过对流体体积分数保持法的改进，将它推广到含多个闭合型的界面重构， 和非闭合型的界面重构。通过算例表明，这种方法能够较好地模拟出多个闭合型界面 的重构；对非闭合型的界面重构，关键是将非闭合型界面重构问题，转化为闭合型界 面重构问题。非闭合型界面的非闭合部分，通过用直线型界面( 或者其它的曲线型界 面) 补充成闭合型界面，最后，将模拟出的界面形状去掉补上的部分，从而达到非闭 合型界面重构的目的。并通过对比，这种由非闭合型转化为闭合型界面的方法，不会 对界面的重构带来不利影响。 5 河海大学硕士论文 第二章v o f 方法的理论基础及发展应用 1 9 8 1 年，h i r t 和n i c h o l se 3 2 1 在j c p 杂志上首先提出了著名的v o f ( v o l u m eo f f l u i d ) 方法。对于运动界面追踪问题的数值研究开创了新思路，通过巧妙地加入了一 个流体体积函数的方程，避开了工程浩大的m a k e r 点方案，成功地进行了运动界面的 数值模拟。 2 1 v o f 方法的基本原理 2 1 1 流体体积方程的定义m 1 v o f 方法是在整个流场中定义一个函数f ，在每个网格中，这个函数定义为一种流 体( 目标流体) 的体积与网格体积的比值。只要求得每个网格上的，值，就可以构造出 运动界面，然后求解物理方程时，可以在界面附近作特殊的精细处理，提高分辨率和 精度。由于每个网格只需保存一个，值，可以大大减少储存量。 设计算区域为q ，流体a 所在区域记为q 1 ，流体b 所在区域记为q 2 ，首先定义 一个函数口( ；，f ) 制= 佗；： 亿， 对于两种不相容的流体组成的流场，口b f ) 满足拉格朗日流体体积的特性： 丝+ “丝+ v 丝：o( 2 1 2 ) a t徼 o y 其中矿= ( “，y ) 是流体的速度场，为单个网格的体积。在每个网格毛上定义乃为口( ；，f ) 在网格上的积分， 弓2 壶黔) q l3 称之为v o f 函数，同样它满足： 望+ “篓+ v 娑：0(210 4 ) + “+ v =l z 4 ， o t缸加 6 第二章v o f 方法的理论基础及发展应用 称之为v o f ( v o l u m eo f f l u i d ) 方程。 由上看出每个网格上的流体体积函数实际上是 一网格中的流体体积 网格体积 显然，f = 1 的网格充满流体a ，我们称之为流体网格( f l u i d ) ；当f = 0 时为空网格 ( v o i d ) ；而0 ， l 的网格，则是含有流体界面的网格，称为界面网格( s u r f a c e ) 。 v o f 方法是通过求解流体体积函数，实现对运动界面的追踪。由于流体体积函数是 单个网格中一种流体体积与网格体积之比，是具有物理意义的。但在求解v o f 方程时， 要保持它的这一物理意义是困难的。首先要求求解格式是守恒的，而且每一时间步得 到的函数值在跨过界面时，即在含有界面的网格上，v o f 函数要0 f 1 ，在其它网格 上必须是1 或0 。 2 1 2 流体体积方程的求解 h i r t 和n i c h o l s 针对界面追踪的一般方法以及通过装配格子内界面的斜率来确定 运动界面的斜率和位置，提出了v o f 方法。在这以后，a s h g r i z 和p o o ，k i m 等人 通过对格子内界面斜率更加精确的构造得到了其它不同的重构方法，尽管他们的精确 性和形式是不同的。r i d e r 和k o t h e ”】，g u e y f f i e r 和l i j i e ，y o u n g s 等人利用了运 动界面的法向矢量近似值来构造界面。r u d m a n t ”却1 ，v i n c e n t 和c a l t a g i r o n e d 8 ， u b b i n k 和i s s a ”1 等人采用了更精细的方法求解流体体积函数方程。 1 h i r t 和n i c h o l s 的直线近似4 0 i 黝 钐笏 钐：厂f ， 产l y - - t j t - - ! f件i 图2 1 1h i r t - n i c h o l s 的重构模板 1 9 8 1 年，h i r t 和n i c h o l s 在提出v o f 方法时，提出的用直线近似的界面重构方法。 7 河海大学硕士论文 它是在一个网格内，用水平和竖直两种直线近似该网格内的流体界面。采用旌主一受主 ( d o n o r - a c c e p t o r ) 的差分逼近格式计算运动界面的流体体积函数方程。 界面的重构如下：将流体的自由表面看做局部的单值函数r ( x ) 和x ( x ) ，采用9 个 网格的模板( 如图2 1 1 所示) ，计算i - 1 ，j ，i + 1 网格列的巧值和j 一1 ，j ，j + l 网格列的蜀 值，估算出每个网格上界面的斜率值d y d x 和d x 痧，然后根据流体体积函数和斜率 的大小确定网格( f ，_ ，) 上的自由面的位置和方向。 耳= f _ | s y , ，l = i - i ，i ，i + 1 f 坐1 ：! f 当二圣! ! d x ) i6 x l n + 2 8 x f + 6 x x t = f x p l = j - 1 ，j ，j + l ( 2 1 5 ) ( 2 1 6 ) ( 2 1 7 ) f ，_ d x l ：上唾掣 ( 2 1 8 ) 、由，) l8 y j + i + 2 6 y j + 6 y i i 如果 削矧 则定义网格内的界面是水平的：否则，网格内的界面是竖直的。 以上，即是h i n 和n i c h o l s 的界面重构的方法。尽管它只是“零”阶自由面格式， 却为以后的高精度界面重构方法开拓了道路。 2 f l a i r 界面重构技术, 3 3 1 1 9 9 1 年，a s h g r i z 和p o o 的f l a i r ( f l u xl i n e s e g m e n tm o d e lf o ra d v e c ti o na n d i n t e r f a c er e c o n s t r u c t i o n ) 方法在任意相邻两个网格内用斜线段近似界面。然后通过 计算一个时间步内流过该网格边界的流体体积量，作为修改流体体积函数的数值流通 量。由于界面构造涉及到两个网格，它分成多种情况分别进行计算。首先利用网格边 界的流体速度，确定施主单元和受主单元，分别用屹，e 表示其流体体积函数，然后 按以下几种情况计算： ( 1 ) 施主网格是满网格，即，n = 1 0 ； ( 2 ) 施主网格是空网格，即昂= 0 0 ； 第二章v o f 方法的理论基础及发展应用 ( 3 ) 施主网格和受主网格都是半网格，即o o 昂 1 0 ，0 0 e 1 0 ； ( 4 ) 施主网格既是半网格，受主网格为空或者满网格，即0 0 ： 亿， 溉：= 瓮貉u i + i 2 - - - ： 汜s ， 对于二维和三维情况，采用方向分裂算法。对x 方向和对y 方向的扫描应顺序交替， 以免引入系统误差。对于方向分裂引起的流体体积分数大于1 的问题，r u d m a n 3 7 1 采用 这样方法处理：在每一时间步开始时，对每个方向计算网格内的流体体积，代替计算 体积分数，利用l a g r a n g e 观点计算半时间后的网格体积，再相除得到体积分数。 2 2v o f 方法的应用 v o f 方法，在每个网格内只需保存一个流体体积函数，大大降低计算储存量，并且 通过巧妙的界面重构，可以得到精细的运动界面，因此受到普遍的关注。它广泛应用 于实际工程中，并得到较好的结果。尤学一等【撕1 采用v o f 法追踪了重力流液液和气液 相界面的迁移，发现v o f 法可很好地追踪强非线性、大形变的相界面随时间的变化； 李维仲等【4 7 1 采用了v o f 方法跟踪界面，并采用p l i c 进行界面重构，模拟了竖直通道 宽度对气泡在液体中的非定常运动、变形以及传热特性的影响，得到了一系列有价值 的结果，并与实验结果比较，表明数值模拟结果与实验结果吻合的较好；沙海飞等1 利用v o f 方法处理自由水面对多孔溢洪道泄流进行三维数值模拟，给出了泄流能力、 水面线和压力分布。数值模拟结果与模型试验结果对比表明，使用的模型能精确的模 1 2 第二章v o f 方法的理论基础及发展应用 拟溢洪道的过坝水流：袁丽蓉等采用标准双方程k 一占紊流模型和v o f 方法建立二维 数值模型，模拟有横流情况下底部排污问题。将回流区的分析结果与试验进行了对比， 吻合较好：端木玉等啪1 采用有限差分法离散流体运动方程，v o f 法追踪自由表面，通 度系数概念处理液舱内的障碍物，计算分析了不同结构形式液舱的防晃效果；t s e n g 等归”采用p l i c - v o f 方法模拟观管网系统中微观井的充液过程，指出井的形状液体的表 面张力对充液结果影响很大；李梅娥等通过s o l a - v o f 方法对金属铸造成型过程中 的流体运动规律进行数值模拟，得出改变浇铸系统浇口和改变液体温度应力分布可以 提高液面稳定、改善充型形态的结论；h i e up h u n gd a n g 等例采用h f - v o f 方法对斜坡 海岸上的波浪翻滚现象进行模拟，指出采用n _ s 双方程耦合求解可改进自由面单元中 的压力和速度插值；付小莉等运用水动力学的n - s 方程和k 一占紊流数值模型，结合 计算水气两相流的方法v o f ，对葛洲坝下游中华鲟产卵河段的流场进行了计算，并结合 中华鲟的生态特性，分析其产卵繁殖对流场特性的要求，期望在三峡水库下游河道整 治中注意保护其产卵场地；鹿子娟等1 用v o f 方法模拟环氧树脂在毛细管中的浸润过 程，得到较好的模拟结果，为纤维束集合体中的轴向毛细浸润模拟奠定了基础。 v o f 方法经过不断地研究发展，界面重构越来越精细，流体输运方程由二维拓展 至三维，差分格式复杂多样，计算网格从二维网格发展到三维网格、正交网格、非正交 自适应网格、无结构网格等等。v o f 方法广泛应用于各种领域，所解决的问题涉及科 学工程计算、水利水电、化学和生物工程、材料加工、国防建设，等等方面，具有广 阔的发展前景。v o f 方法作为运动界面追踪的重要方法，将在未来发挥越来越大的作 用。 2 3 本章小结 本章主要介绍了v o f 方法的原理和发展应用，总结了v o f 方法中，各种界面重构 的方法和捕捉界面的方法，对比了它们的优缺点。由h i r t 和n i c h o l s 首先提出了通过 界面重构的方案，解决运动界面的追踪问题。从此以后，对运动界面追踪的研究，由 开始的捕捉类方法，捕捉类方法加半构造的性方法，发展到完全装配和构造性的重构 型方法。在这个过程中，v o f 方法得到不断的创新和改进。 各种v o f 型运动界面追踪的方法的提出，丰富了v o f 方法的内容，拓展了v o f 方 法的应用范围，奠定了v o f 方法在现代运动界面追踪的地位。使得v o f 方法能够在各 种领域得到广泛应用，具有广阔的发展前景。 河海大学硕士论文 第三章流体体积分数保持法及其优化算法 在运动界面追踪问题中，采用数值模拟的方法“捕捉”界面时，由于数值方法的 数值耗散和数值色散性、非线性效应等，可能所捕捉到的只是一个模糊的轮廓，所以 要采用某种技术处理或者几何的装配方法，重新将运动界面锐利化甚至曲面化。这就 是所谓的界面重构。事实上，近来对于运动界面的重构，有许多学者不再用“捕捉” 类方法，即使是高分辨率方法，如) , j s c l ，e n o ，w e n o 等也不用，而是直接采用装配， 或者称为再构造的方法，甚至每一时间步都采用运动界面的“重构”技术，去实现界 面追踪。譬如y o u n g s ，u b b i n k 和i s s a 的方法。这种更强调网格内界面的运动和重构 的方法，在当今尤其显示出它们的优势，日益成为v o f 型运动界面追踪和重构方法的 主流。 v o f 型运动界面重构的一般步骤如下【”】： ( 1 ) 给定初始运动界面位置，即流体体积函数在某个时刻f = 厶的分布为r , 3 ； ( 2 ) 通过流体体积控制方程，计算f = o 。时刻的辱岁； ( 3 ) 采用直接或者几何的装配方法，重构f = 时刻的e ，； ( 4 ) 检验# 了1 是否在【o ，1 】的取值范围内，再进行必要的技术处理。 3 1 基本概念 针对以上阐述的运动界面重构的过程中的第三、四步，本文提出了一种在矩形网 格上的运动界面重构方法一流体体积分数保持法。在二维的矩形网格上，该方法在 单个网格内用斜线段近似运动界面，并要求相邻网格内的斜线段在公共边上的交点重 合，通过保持网格内的流体体积分数不变，从而在界面网格上建立非线性代数方程组， 最后通过求解该非线性代数方程组，确定运动界面的边界与网格线的交点位置，从而 重构出运动界面的形状和位置。为叙述方便，先引进如下一些概念。 假设一个网格中最多只存在一条可用直线段近似的边界。与某网格有一条公共边 的网格，称之为某网格的邻网格。我们把0 f 1 的网格或f = 1 但至少有一个邻网格 满足f = 0 的网格称为界面网格。如果该网格的f = 1 ，该网格的四个邻网格的，值都 大于零，则称该网格为内网格。用1 - 3 6 记单元编号，1 1 h 2 0 陋型点编号。如图3 i 1 所示，7 ，1 3 ，1 9 ，2 5 ，2 6 ，3 2 ，3 3 ，3 4 ，3 5 ，2 9 ，3 0 ，2 4 ，1 8 ，1 2 ，l l ，5 ，4 ，3 ，2 ， 8 号网格是界面网格。如果边界所围内部为流场，则9 ，1 0 ，1 4 ，1 5 ，1 6 ，1 7 ，2 0 ，2 1 ， 1 4 第三章基于矩形网格的v o f 运动界面重构的流体体积分数保持法 l 乡 14习 石 1 19 i 1 1 目硼 露， 门2 口 小目2 旬萨 91 口 闷 一毛 1 41 5 161 7瞄 7 习 ：_ j 。9 7 2 02 12 22 32 4 j丁习 甥 2 6 2 72 8 2 哑 幽 3 d t， 3 l 固习。3 3 目3 4 目 掣3 6 两q 码z l飘码 图3 1 1 网格及型点示意图 2 2 ，2 3 ，2 7 ，2 8 号网格为内网格。把近似界面的直线段与网格线的交点日q 做型点，且 假设两个邻网格公共边上的型点重合( 如图3 1 1 所示的点) 。一般情况下，在含有 近似边界的网格内，其网格的边界线上有且只有两个型点。由所有界面网格按正方向 ( 规定沿此方向行走，流场位于左手边) 顺次相邻连接组成的区域称为网格界面，其 包含实际边界和近似边界。如果网格界面能围成一个闭合回路，则称该网格界面为闭 路网格界面( 图3 1 1 所示的界面网格构成闭路网格界面) ，否则称之为非闭路网格界 面。本文只对闭路网格界面情况进行研究。对网格界面上的网格按界面正向进行顺序 连接并给予编号，如图3 1 1 所示，7 ，1 3 ，1 9 ，2 5 ，2 6 ，3 2 ，3 3 ，3 4 ，3 5 ，2 9 ，3 0 ， 2 4 ，1 8 ，1 2 ，1 1 ，5 ，4 ，3 ，2 ，8 号网格，可分别编号为1 ，2 ，2 0 ，则对应各网 格上的型点编号记为纠，臣0 | ，要求型点1 1 陋于界面网格的最后一个网格和1 号 网格的公共边上，型点p 陋于界面网格的1 号网格和2 号网格的公共边上，其余依次 类推。 沿着界面网格的正方向行走，按每个界面网格与其前后相邻的界面网格的位置关 系，一般可把界面网格划分为3 种类型，而每种类型又可分为4 种情况( 图3 1 2 ) 。 扎 h 河海大学硕士论文 片 ( c ) b b b c a ) 类型1 ：直线型 ( d ) 类型2 ： 类型3 ：内凹型 图3 1 2 界面网格分为3 类共1 2 种情况 c a ) 外凸型 图3 1 2 中的是沿着界面网格正向行走的次序编号，从走到，再走到 。“i ”记当前界面网格在网格界面中的网格序号。类型l 的特征是，与当前网格前 后相邻的界面网格的中心点在同一条直线上，或平行工轴或平行y 轴，称为“直线型”； 类型2 的特征是，与当前网格前后相邻的网格的中心点不在同一条直线上，但流场的 边界往外凸，称为“外凸型”；类型3 的特征是，与当前网格前后相邻的网格的中心点 不在同一条直线上，但流场的边界往内凹，称为“内凹型”。每一种类型的界面网格， 根据界面走向与工轴正负向和y 轴正负向的平行关系，又可分为如图3 1 2 所示的( a ) ， ( b ) ，( c ) ，( d ) 四种情况。各种情况图中所示的q ，为界面型点，0 为紧邻q ，的位 于流场一侧的网格节点，作为将来确定型点坐标的参照点。引进型点位置变量 f ，= 1 只q ，i 缸，( 只9 平行于x 轴时) 或0 = 1 只q f i 缈，( 只q ，平行于y 轴时) ，其中缸， a y , 分别为第f 个网格的网格边界在x 向和y 向的长度。记录与每个型点q ，相关的信息， 包括型点编号，界面网格编号等于型点编号，界面网格类型，以及参照点p ，的坐标 ( x 弓，y 厅) 。 如果第，个界面网格上的流体体积分数为正，则可根据界面网格的类型，建立保 1 6 图 瞄b 鼠酱庇 第三章基于矩形网格的v o f 运动界面重构的流体体积分数保持法 持流体体积分数的对应方程如下( 见表3 1 1 ) 。 表3 1 1 界面网格类型与保持流体体积分数 第，个界面网格类型保持流体体积分数的方程 类型1 ：直线型 圭( 。) = 五 类型2 ：外凸型 1 i 。= ，f 类型3 ：内凹型 1 一三( 1 一t d o - - t l + 1 ) = 石 假设网格界面由n 个界面网格组成，则当1 = 1 , 2 ，n 时，可形成含个方程的非 线性代数方程组，此时f 。= t l 。求解该方程组，得到f ，= 1 , 2 ，n 。根据记录的型 点q 的相关信息，最终可以计算出型点的坐标( 乜，y 。) ( 计算公式见表3 1 2 ) 。 表3 1 2 型点q ，的坐( x 岛，y 岛) 计算公式 第，个界面网格四 型点g 种情况 x q iy e , ( a ) y 一t i ( b ) 石器+ f j y 矗 ( c ) 七t | ( d ) x 一t | yp f 3 2 流体体积分数保持法的最速下降法 本文介绍的运动界面重构算法一流体体积保持体积分数法对界面网格划分为三种 不同的类型，通过巧妙的几何构造，找出界面与网格线的交点，从而刻画界面形状。 基于上节的基本概念的介绍，本节将具体地介绍该方法，基本步骤如下： 第一步： 在整个计算区域中，搜索界面网格，并通过矩阵彳来记录搜索情况。界面网格， 即体积分数y , a ，v f ，_ ，满足：o 1 ，或曩。，= l 但与一个空网格相邻。根据其定 义确定矩阵a 中的各元素如下： 1 7 河海大学硕士论文 若位于第m 行行列的网格为界面网格，则记以，= - 1 ；否则，以，= o 统计界面 网格的个数，并计下界面网格总数为。 第二步： 将搜索到的界面网格，即以。= - 1 的网格进行编号用k 累计界面网格的个数。 编号原则：( 图3 2 1 中位于第m 行刀列的网格为中心网格) m l m r e , n - 1 m ，1 1m n + l m + l , n 图3 2 1 中心网格 ( 一) 首先确定第一个