（计算机软件与理论专业论文）流体模拟的混合模型研究.pdf

摘要 随着影视、游戏和虚拟现实行业的迅速发展，这些行业对新技术、 新效果有了更多的追求，好的效果能带给人们更多的身临其境的感 受。 在影视、游戏和虚拟现实领域，基于物理的模拟在特效制作中渐 渐成为主流。一方面，人们追求更高质量的特效，物理模型能够更好 的反应物理规律，呈现更自然的活动规律；另一方面，早期的特效制 作耗费成本很大。例如建立一个城堡模型，通常是建一个按比例缩放 的真实模型，然后采用一系列手段来产生真实的效果，效率低下。流 体模拟作为影视动漫特效的重要元素，具有重要意义。例如在怪物 史莱克，指环王，水世界等电影大片中，流体模拟震撼人心， 游戏行业同样也极为重视流体模拟的实时效果。 由于流体模拟的建模和模拟计算复杂，方法众多，本文在对前人 方法的总结和归纳的基础上提出了一种新的耦合的计算模型来模拟 两种流体的交互。文中比较欧拉方法和拉格朗日方法各自的优劣，对 流体受热运动进行了深入分析，把整个运动分为两部分：第一部分是 温度场自身的温度扩散和对流产生的速度，第二部分是流体自身所受 的压力、浮力、粘滞力、表面张力等。综合这两部分的计算特点，分 别使用稳定流方法和s p h 方法来进行模拟计算，然后建立数学模型对 这两种方法进行耦合。同时这种耦合的方式能修正n s 方程，以达到 使流体运动更自然的效果，然后使用m a c h i n gc u b e s 算法来进行表面 重构，最后绘制等值面，使用光照模型进行光照计算机，实现真实感 效果，能取得比l e v e ls e t 方法更优的时间效率。 本文第一章简述了流体模拟的发展历史和前人的研究成果，第二 章介绍了基于物理模型的流体模拟的理论基础，第三章针对流体受热 运动进行理论分析、建立物理模型，对欧拉和拉格朗日方法进行耦合。 第四章实现模拟过程，并针对热传播过程修正物理模型，给出实验结 果。第五章是对实验的分析和对理论发展的展望。 关键词稳定流；s p h ；欧拉法；拉格朗日法；外力修正；耦合； a b s t r a c t w i t ht h er a p i dd e v e l o p m e n to ff il m ，g a m e sa n dv i r t u a l r e a li t yi n d u s t r y ，t h e s ei n d u s t r i e sh a v e - m o r ed e m a n df o rn e w t e c h n o l o g i e sa n dn e we f f e c t s g o o de f f e c t sc a nb r i n gp e o p l e m o r ev i v i df e e l s i nt h e s ea r e a s ，p h y s i c s b a s e ds i m u l a t i o ni ns p e c i a le f f e c t s c r e a t i o ng r a d u a l l yb e c o m et h em a i n s t r e a m o nt h eo n eh a n d ， p e o p l eh a v eh i g h e r l e v e lp u r s u i tf o rt h eq u a li t yo fs p e c i a l e f f e c t s ，t h ep h y s i c a lm o d e lc a n b e t t e rr e s p o n dt op h y s i c a ll a w s a n dp r e s e n tam o r en a t u r a la c ti v itie s ：o nt h eo t h e rh a n d ，t h e e a r l ys p e c i a le f f e c t sc r e a t i o nc o s t sl a r g e l y f o re x a m p l e ，t h e e s t a b l i s h m e n to fac a s t l em o d e l ，w eo f t e nb u i i di tb ys c a l i n g ar e a lm o d ela n dt h e nu s eas e rie so fm e a s u r e st op r o d u c em o r e r e a l i s t i ce f f e c t sa n di n e f f i c i e n c y f l u i ds i m u l a t i o n ，a so n e o ft h ef i l me l e m e n t s ，p l a ya ni n c r e a s i n g l yi m p o r t a n tr o l e f o r e x a m p l e i nf il m s “s h r e k “l o r do ft h er i n g s “。w a t e rw o r l d 6 a n ds oo n ，p e o p l ew a ss h o c k e db yf l u i ds i m u l a t i o n a tt h es a m e t i m e ，t h eg a m ei n d u s t r yi sa l s oi n c r e a s i n g l yc o n c e r n e da b o u t t h ei m p a c to fe f f e c t so nt h eg a m e a st ot h ec o m p l i c a t i o no fc o m p u t a t i o ni nf l u i ds i m u l a t i o n m o d e l i n ga n ds i m u l a t i o na n dt h em e n t h o d sa r ee n o r m o u s ，t h i s p a p e rr a i s e san e wc o u p l e dc o m p u t a t i o n a lm o d e lt os i m u l a t et w o k i n d so ff l u id i n t e r a c ti o nb a s e do ns u m m a r izi n ga n d g e n e r a li z i n gt h ep r e d e c e s s o r sm e t h o d s i nt h i sp a p e r ，g o o da n d b a do ft h ee u l e rm e t h o da n dl a g r a n g em e t h o da r ea n a l y z e d r e s p e c t i v e1 y i nd e t a i1 ，t h em o v e m e n to fh e a t e df l u idis i i i a n a l y z e d ，t h ew h o l em o v e m e n ti sd i v i d e di n t ot w op a r t s ：t h e f i r s t p a r t i st h e v e l o c i t yg e n e r a t e db y t h e t e m p e r a t u r e d i f f u s i o no ft e m p e r a t u r ef i e l di t s e l fa n dc o n v e c t i o n ，t h e s e c o n dp a r tist h ep r e s s u r e ，b u o y a n c y ，v is c o u sf o r c e ，s u r f a c e t e n s i o no ff l u i di t s e l f i nv i e wo ft h e r e s p e c t i v e c h a r a c t e r i s t i c so ft h e s et w op a r t s ，t h ef l u i ds i m u l a t i o ni s c o m p u t e du s i n gt h e s t a b l ef l o wm e t h o da n d s p hm e t h o d r e s p e c t i v e l y t h e nam a t h e m a t i c a lm o d e li se s t a b l i s h e da n dt h e t w om e t h o d sa r ec o u p l e d a tt h es a m eti m et h isc o u p li n gm e t h o d c a nc o r r e c t i n gn s e q u a t i o n ，i no r d e rt om a k ef l u i dl o o k s n a t u r a l a n dt h e ns u r f a c er e c o n s t r u c t i o ni sc o m p l e t e db yu s i n g m a c h i n g c u b e sa lg o r it h m a tl a s t ，t h eis o s u r f a c eis r e n d e r e d ，w i t h1i g h t i n gm o d e lt h eli g h t i n gc o m p u t a t i o n sa r e ；a c h i e v e d ，t h er e a l i s t i ce f f e c ti sr e a l i z e du l t i m a t e l y ，t h et i m e e f f i c i e n c yi sb e t t e rt h a nl e v e ls e tm e t h o d t h i sf i r s tc h a p t e ro u t l i n e st h eh i s t o r yo ft h ed e v e l o p m e n t o ff l u i ds i m u l a t i o na n dp r e v i o u sr e s e a r c hr e s u l t s t h es e c o n d c h a p t e rd e s c r i b e st h et h e o r e ti c a lb a s isf o rf l u i ds i m u l a t i o n b a s e do nt h ep h y s i c a lm o d e l t h et h i r dc h a p t e ra n a l y z e st h ef l u i d h e a t e dm o v e m e n tt h e o r e t i c a l l ya n dc o u p l e de u l e rm e t h e dw i t h l a g r a n g em e t h o da f t e rs e l e c t i n gt h ep h y s i c a lm o d e l t h ef o u t h c h a p t e rr e a l i z e st h es i m u l a t i o np r o c e s sa n dc o r r e c t sp h y s i c a l m o d e li nh e a td i f f u s i o na n da c h i e v e se x p e r i m e n t a lr e s u l t s t h e f i f t h c h a p t e ra n a l y z e s t h ee x p e r i m e n t a g a i na n dg i v et h e e x p e c t a t i o n sf o rt h e o r yd e v e l o p m e n t k e yw o r d s ：s t a b l ef l u i d s ：s p h ；e u l e rm e t h o d ：l a g r a n g em e t h o d ： e x t e r n a lf o r c ec o r r e c ti o n ：c o u p li n g ： 湖南师范大学学位论文原创性声明 本人郑重声明：所呈交的学位论文，是本人在导师的指导下，独 立进行研究工作所取得的成果。除文中已经注明引用的内容外，本论 文不合任何其他个人或集体已经发表或撰写过的作品成果。对本文的 研究做出重要贡献的个人和集体，均已在文中以明确方式标明。本人 完全意识到本声明的法律结果由本人承担。 学位论文作者签名： 秘孵 2 0 o 年岁月弓日 湖南师范大学学位论文版权使用授权书 本学位论文作者完全了解学校有关保留、使用学位论文的规定， 研究生在校攻读学位期间论文工作的知识产权单位属湖南师范大学。 同意学校保留并向国家有关部门或机构送交论文的复印件和电子版， 允许论文被查阅和借阅。本人授权湖南师范大学可以将本学位论文的 全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或扫 描等复制手段保存和汇编本学位论文。 本学位论文属于 1 、保密口，在年解密后适用本授权书。 2 、不保密日。 ( 请在隧上相应方框内打“ ) 作者签名： 槲 日期：力加年岁 月； 日 剥醛各离蠢呜 日期肿年岁肌日 流体模拟的混合模型研究 1 绪论 1 1 引言 计算机图形学是流体力学发展的受益者之一，其主要是用电子计 算机和离散化的数值方法对流体力学问题进行数值模拟和分析。在早 期的流体模拟中，由于计算能力有限，主要采用参数建模的方法，具 有代表性的是1 9 8 6 年d a r w y nr p e a c h e y 博士发表的海面和波浪模拟 的文献n 1 ，文献中提出使用波形函数的组合来构造波浪函数，并处理波 的折射问题，然后采用粒子系统来模拟当波浪破碎或者遇到障碍物时 形成的浪花。随着计算机硬件的不断发展，计算能力不断提升，使的 基于物理模型的流体模拟迅速发展起来，本文主要探讨基于物理模型 的流体模拟，其中涉及到的现象包括自由表面流体、气泡等，这些现 象的模拟在计算流体力学领域是很热门的话题，而且工程软件多如牛 毛，其中商业化的软件如：c f x ，f l u e n t 等，有人会质疑为什么不直接 使用这些商业软件计算，然后使用绘制技术展现效果呢? ! 其实c f d 更强调对问题本身的精确求解，然后给出一个合理的分析；而计算机 图形学则是力图展现流体的视觉效果，需求不同，所以研究层次不同， 而且c f d 的研究者们根本就没有设计一个很好的计算模型来描述视觉 效果，这就需要图形学的研究者们去需求有效的解决方案了。在图形 学中好的效果和速度比精确性来的更重要，小的计算代价和较高的交 互性能为动画制作人员提供良好的控制方案，最终达到降低劳动成 本，提高制作效率的要求。目前市面上的三维图形软件产品有很多都 包含了流体模拟的模块，如3 d m a x 提供了g l u 3 d 流体插件；m a y a 提供 了f l u i de f f e c t s 模块；独立软件包如：r e a l f l o w 和r e a l w a v e 。这 些软件都表现的出类拔萃。 最后要提的是流体表面绘制的问题，对于自由边界的问题通常采 堡主堂垡笙壅 用m a r c h i n gc u b e s 方法来重新构造多边形表面，获得等值面的数据 后就可以采用几何体绘制的方法进行绘制，如果要产生真实感效果还 可以添加光照模型用于与表面交互。 1 2 流体模拟的发展简介晗3 早期的流体模拟，由于计算能力有限，主要采用参数化建模的方 法，如文献通过使用一系列简化了的线性波形函数的组合来拟合一 个波浪函数，从而合成水体表面的高度场，并处理波折射问题，并采 用粒子系统来模拟波浪破碎和碰撞到障碍物时形成的浪花。但由于文 献中的水粒子和网格只是在其初始位置附件运动，所以他们都无法表 示真正意义上的流动效果，也无法处理边界给水面带来的影响。如f f t 经验模型口1 可以很好的描述波幅较小的海面，但是对于变形较大，无 规则的水体，人们觉得控制困难，而且很多细节无法呈现，为了获得 更多更丰富的效果，研究者开始转向研究基于物理的模型。 基于物理的方法总共可以分为两种：第一种称为欧拉法，他是从 研究流体所占据的空间中各个固定点处运动着手的，分析空间中每一 个固定点上的流体的密度，速度，压强，温度等这些在每一时刻都变 化的参数，以及空间位置变化对网格带来的影响，这是种基于网格 的方法；第二种方法称为拉格朗日法，它与欧拉法刚好关注点相反， 它把流体看做由若干小微团组成，并跟踪分析这些微团的运动，分析 某一指定微团的速度，密度，压强，温度等随时间变化的参数，以及 在发生相变的时候研究一种微团转变为其他微团时候参数的变化，以 及整个流体的运动，这是一种基于粒子的方法。这两种方法各有优劣， 欧拉法的计算效率较高，而拉格朗日方便计算流体的自由边界，目前 还出现一些混合方法如：p a r t i c l el e v e ls e t 方法。 1 2 1 欧拉法 目前描述流体现象最完善的方程是n a v ie r s t o k e s 方程h 1 ，这个方 程是根据牛顿第二定律推导出来的，以下是n s 方程的欧拉形式， 流体模拟的混合模型研冤 v u ：0( 1 ) 质量方程( 不可压条件) i o u = 一( “v 扣+ 内2 “一跏+ f ( 2 ) 动量方程 d z 其中p 为密度，p 为压强，厂为体积力( 或外力) ，”为速度，y 为 运动粘性系数，欧拉法是将上述量离散到网格上去，然后计算各个网 格节点上的这些量的变化，从而得到整个场。目前流体的网格化方法 主要有两种，一种是把标量分布在网格中心，矢量放在网格表面，这 种方法比较容易保证守恒性条件，另一种是把所有的量都处于同一个 位置，这种方法不需要太多的插值运算，简单易于实现。 1 2 2 拉格朗日法 拉格朗曰法是一种基于粒子的方法，其n s 方程的形式和欧拉形 式略有不同， d u ：v v 2 “一跏p + 厂( 3 ) d t 由于拉格朗日法自身的定义，这种方法自然的保证了质量守恒， 所以相对欧拉形式少了一个平流项。如果将等式右边整理成一个力， 那么等式( 3 ) 就退化为牛顿第二定律： a t = d u i d t = l _ | p l 其中a ，是粒子i 的加速度，“，是粒子i 的速度，z 为该粒子受到的 合理，肛是粒子所在位置的密度。通过引进核函数，控制邻域半径， 积分计算粒子在下一时刻的位置和其他状态量。 拉格朗日法的优点是易于表达，无需处理整个空间而且自然的保 证质量守恒，最主要是比较容易控制。但是其缺点是由于依赖核函数， 在流体边界处领域半径内可能没有粒子，使得部分粒子的作用力失 效，边界出现弱效果。同时计算机随着粒子数的增多而加大。 硕士学位论文 1 3 国内外研究现状 在计算机图型学中，与流体相关的内容很多，为了满足人们观赏 需要，各种各样的流体模型被大量使用，5 1 ：1 - 烟雾、云彩、火焰、爆 炸、自由运动界面流体、混合流和多相流、甚至还有非真实感流体。 在早期的流体模拟中，学者们提出很多解决方案，如：文献 5 】中采 用隐式格式构造三对角方程组，快速稳定的求解水面高度场。文献【6 】 采用高度场来描述流体表面，然后针对物体与液体交互区域来施加作 用力，用粒子来模拟飞溅效果。 文献r 7 】首次提出在图形学中引入s p h ( s o o m t h e dp a r t i c l e h y d r o d y n a m i c s ) 方法来模拟流体，处理大变形物体。在文献【8 】中作者 通过s p h 来模拟自由表面流体，并考虑张力的作用，然后采用 m a r c h i n gc u d e s 算法来绘制表面。下面根据流体的两大类介绍国内外 最新研究成果。 1 3 1 欧拉法流体 如图1 1 ，是欧拉流体的网格结构，欧拉流体包含了单相流、自 由界面流等等。由于f o ss t a m 的稳定流【9 】的贡献，欧拉流体的数值稳 定性得到解决，由此半拉格朗日对流项与压力投影法在不可压流体领 域得到了广泛的应用，如火焰，烟雾，云等等的模拟。因为他们都不 存在自由运动界面的问题所以处理起来比较容易。而像液体则属于自 由界面流了，基于网格的方法中，目前已经得到广泛应用的自由界面 方法有：m a c ，l e v e ls e t ，v o f ，f e m ，f d m ，p a r t i c l el e v es e t 等。 目前最受关注的是l e v e ls e t 和p a r t i c l el e v e ls e t 两种方法。 流体模拟的混台模型研究 躺惑 一 ff 黼测参黪!鹫孬 霾 藏 糕 麟露懑 缫戮瓣 懑辫蕊 黛 隳漆鬻 嚣 霪 略“拳 隧一鬻 糕 j _ ? t 一。l g _ 图1 1 二维欧拉网格流体结构图 l e v e ls e t 方法是引入一个距离函数，定义为点到边界的有符号 距离，该函数在速度场的推动下移动。由于l e v e ls e t 方法的体积损失 很严重，所以文献【l 町提出在液体内部分布一些虚粒子来标示边界，在 曲率小的地方忽略粒子的影响，在曲率大的地方用来修正距离函数d 的值。由于这个修正启发了学者更进一步的研究，不久d o u g l a se n r i g h t 等在文献 1 1 e o 了提出了p a r t i c l e l e v e ls e t 方法，这种方法是在表面内 外很小距离的附近分布若干粒子，利益他们来检查界面的不精确度， 通过粒子局部的l e v e ls e t 函数。来修正，从而保证界面原来的高曲 率特征，以获得精确的视觉效果。t o l g ag o k t e k i n 等在文献【1 2 】中采用 p a r t i c l el e v e ls e t 来捕获自由运动界面，他们使用双倍采样来减轻体积 损失的问题。 v o l u m eo f f l u i d ( v o f ) 方法也是一种基于网格的方法，其特点是易 于满足体积守恒，它通过定义每个单元的流体体积比来捕获自由运动 界面，方程形式和l e v e ls e t 基本类似。但是通常这种方法并不单独使 用，如文献 1 3 】中结合了v o f 和l e v e ls e t 方法，同时这种方法的缺 点是无法处理边界的密度、压力和粘性等的不连续性，所以无法表现 表面张力。为了解决这个问题，有学者使用了g h o s t f l u i d 方法来避免 硕士学位论文 边界处的不连续，这种方法是利用函数矽将含有两种流体的流场分成 含有一种流体的两个流场，设边界内的流体分布为流体1 和流体2 ， 其中流体l 满足如0 ，流体2 满足矽0 。对于只含有流体2 的流场边 界内的点位真实点，具有真实的物理量，而边界外的点为虚拟点，其 为真实物理量在外面的延拓。随着影视业的发展，人们开始对真实感 有了更高的要求，人们越来越关注流体的细节，如在火焰模拟中，由 于存在数值耗散，需要加强对流效果，文献 1 4 1 提出了添加漩涡约束 来增强细节即( 国：v x u ) ，还有文献 1 5 中使用了拉格朗日示踪和欧拉 网格相结合的方法来增强对流效果，文献 1 6 贝j j 使用压强越阶条件和 速度梯度不连续性来对表面张力建模。又如液体中气泡的产生，文献 1 7 提出了气泡产生的物理模型，并且考虑了因流体热量变化进而产 生气泡的变化过程，并对气泡在流体中的运动做了深入的分析和研 究。不过到目前为止对于流体的细节的改进都只是在n s 方程的基础 上作出的改进，由于n s 方程本身是针对大规模流体简历的。所以很 多改进都不能很好的模拟动力学的特征，因此需要学者从建模的角度 来分析考虑流体细节的运动方式，构造合适的模拟方程，追求更好的 效果。欧拉流体在细节方面虽然考虑的越来越细致，但是时间消耗也 同时上升了不少，虽然前面有提到过半拉格朗日法和隐式迭代的引入 保证的数值的稳定，提高了迭代时间步长，提升了计算效率，但是这 只是解决了其中一个部分，随着l e v e ls e t 方法成为主流方式，由于 l e v e l s e t 方法自身的限制，每一时间步符号函数需要重新初始化，并 迭代至收敛，所以实时性难以满足要求，如文献 1 8 】 1 9 】中直接求解二 维n s 方程，然后结合一个中尺度的k o l m o g o r o v 速度场来近似逼近 三维情况，虽然速度有所提升，但是细节刻画就不满足要求并且仍然 不是真正意义上的实时。 迹堡量型塑堡量堡型塑塞 1 3 2 拉格朗日法流体 拉格朗日法是基于粒子系统的方法，它是对各个相对独立的粒子 进行受力分析，通过积分计算出这些粒子下一时刻的空间位置和速 度、密度、压力等状态量，图1 2 为拉格朗日法示意图。 图1 - 2 二维拉格朗日流体结构图 s p h ( s m o o t h e dp a r t i c l eh y d r o d y n a m i c s 】方法是粒子系统的一种插 值方法，通过构造光滑滤波函数来逼近狄拉克j 函数，来模拟粒子间 的作用力。也可以用来求解柔性物体的大变形。文献 1 8 】通过s p h 来 模拟自由表面流，并采用p o m ts p l a t t i n g 技术来绘制表面。文献1 1 9 1 采用m p s ( m o v i n g p a r t i c l es e m i i m p l i c i t ) 方法模拟多种类型的流体，该 方法将n s 方程转化为运动粒子相互作用的形式。基于拉格朗日方法 的流体表面捕获算法，如文献 1 8 】 2 0 】采用了等值面来构造多边形的方 法，这需要先确定流体的表面粒子。而文献2 l 】则是构造粒子系统下 的l e v e ls e t 方程，基于网格来重建表面。相比欧拉法来说拉格朗日法 更方便于处理自由表面流动问题，同时拉格朗日也有很多问题存在， 如s p h 存在数值稳定性方面的问题，边界弱化，同时还要求粒子规则 分布，当然随着c f d 的发展，越来越多的学者开始研究s p h 方法， 并提出了系列的优化和改进方法，如w k l i u 等人提出来一种再生 核粒子法瞰】来提高计算精度，s w e g l e 、d y k a 和c h e n 等提出了s p h 方 硕士学位论文 法不稳定的起因及稳定化方案【2 3 】和前面欧拉的g h o s tf l u i d 方法类似， 有人提出使用虚粒子来解决边界弱条件的问题，使得s p h 流体在边界 上能处理自然了。 欧拉法和拉格朗日法各有优缺点，为了更真实的模拟流体运动， 这两种方法往往都结合在一起使用，如文献 9 中的半拉格朗日法。总 之，两种方法的结合是今后流体模拟发展的趋势。 1 4 本文工作 本文建立了一种新的欧拉和拉格朗日耦合的模型，用来模拟流体 在受热情况下的运动情况，该模型可以进一步扩展，用来模拟液体沸 腾的效果。本文的主要内容包括： 第一章介绍流体模拟的国内外研究现状，和一些常见方法的优劣 对比。 第二章介绍了基于物理的流体理论基础，对稳定流方法和s p h 方 法进行详尽的阐述。 第三章针对液体受热运动的理论建立物理模型，对欧拉和拉格朗 日方法进行耦合。 第四章实现整个模拟过程，并针对热传播过程修正物理模型，追 求较好的视觉效果。 论文最后进行工作总结，提出下一步的研究方法和任务。 流体模拟的混合模型研究 流体模拟理论基础 2 1n a vie r - s t o k e s 方程 n a v i e r - s t o k e s 方程，是描述粘性不可压流体动量守恒的运动方 程，简称n _ s 方程，由c l a u d e - l o u i s 和g e o r g eg a b r i e ls t o k e s 分 别于1 8 2 1 和1 8 4 5 年提出。这个方程是针对大规模流体提出的最完善 的流体力学方程。其欧拉形式如下： 譬= 一( “v 咖一土跏+ v v 2 u + 厂 ( 4 ) v “= 0( 5 ) 其中式( 4 ) 是动量守恒方程，式子( 5 ) 是质量守恒方程。式子 中，“是速度矢量，户是流体密度，v 是动力粘度系数，表示作用于 流体的外力。首先，式( 5 ) 表示散度为0 ，以确保流体的体积守恒( 即 某表面的流进体积和流出体积一样多) 。由于不可压流体密度恒定， 那么体积和质量都是守恒的。其次，式( 4 ) 左边定义了流体加速度。 速度对时间的全导数是 d u ( r ( t ) ，f ) o u 缸o u 砂锄瑟抛 一= - i - 4 - - 一 d to xa t a ya t a za t a t = 罢+ 罢，考，警 妄，考，象丁 。6 ， 2 百+ l 瓦万西儿石亩瓦j r 烈 = 曼a 竺t + ( ia a x - ，专，鲁) “2 一“ i，瓦，瓦i 澎 ：丝+ 甜v “ 对流( c o n v e c t i o n ) 是流体内分子最常见的运动方式，一般分为 两种形式： 扩散( d i f f u s i o n ) 流体内分子随机运动的表现形式 平流( a d v e c t i o n ) 流体内部分子的整体运动，比较明显的如：在多 相流中，流体携带其他物质在水平方向上运动的表现形式。除去对流 硕士学位论文 项之后，式( 4 ) 右边三项分别为压力项一土即，粘性力项v v 2 “和外力 p 项厂，其中外力主要包括重力，浮力，表面张力等。因为式( 5 ) 的条 件，且又满足牛顿第二定律，所以不可压流体又称为牛顿流体。 除了欧拉形式，n _ s 方程还有拉格朗日形式： 皇竺：一土v p + 1 ，v 2 掰+ f( 7 ) d l p 由于拉格朗日形式下的流体是完全由粒子来构成的，每个粒子保 持自己的质量，所以不需要额外的质量守恒方程，而且因为要引入光 滑核函数，用来起到聚类的作用，反映了粒子与周围粒子的作用，所 以在描述单个粒子运动的时候也不需要对对流项进行计算，这就是 ( 7 ) 式和( 4 ) 式不同的原因。 2 2 稳定流方法 基于网格的欧拉方法始于f o s t e r 1 9 9 6 年他首次提出了通过使用 有限差分在均匀网格上求解n - s 方程，但是需要较短的时间步长来保 证数值稳定性。s t a m 在此基础上提出了著名的半拉格朗日方法p 1 来求 解n - s 方程，并且保证在大时间步情况下仍然能够保证数值稳定性。 一般求解n - s 方程时，用稳定流方法有以下几步：( 1 ) 计算外力项， ( 2 ) 计算平流项，( 3 ) 计算扩散项，( 4 ) 对前面的计算结果做投影。 对不同的流体来说，外力可以是浮力、重力、表面张力、漩涡力等等。 从( 4 ) 式我们可以知道，n - s 方程实际上转换为加速度和力的关 系式了，即： a = p r e $ s u 增+ 协c 0 。缈+ 础棚， ( 8 ) 首先来看外力项，由于本文的需要，主要考虑热力和漩涡限制力， 按照物理定律，浮力计算式为：k = c r ( t t o ) j - g 其中j 是垂直的方向矢量标志，r o 为给定环境温度，o r 为控制浮 力大小的系数，g 为重力。即，随着温度的增大，所受的浮力也增加。 流体模拟的混合模型研究 由于要模拟液体受热运动，所以考虑到液体是一种粘性流体，需要增 加漩涡限制来增加受热运动的表现力，同时防止受热运动不合物理规 律。漩涡限制的计算公式如下： 无= e ( 。i x w ) a x ( 9 ) 甲= 墨，r = vl 国i ，缈= v xu ( 1 0 ) i7 7 其中s 是比例系数，是旋度。再令厶删= 允，+ 丸，外力就计算 完成了。在求解n - s 方程的过程中，如果记要计算的状态量，如：密 度，速度，温度等等为w ( x ) ，则有： w ( 曲= w o ( x ) + 矿( 工，f ) ( 11 ) 在计算完外力项之后，就要计算对流项了，对流项是对流体携带 流体中的标量运动的描述，如：温度，密度等标量。对流项使得流体 携带这些标量沿着流体运动方向传输，用来弥补欧拉网格对空间量转 移无法计算的不足。如果描述成型如( 11 ) 式的形式，则有： w l ( x t + a t ) = w o ( x ) + u ( t ) a t ( 1 2 ) 但是这里有个问题，采用显示格式的积分，在较大时间步长中会 发生数值耗散，从而对一些参数做调整都会使计算结果剧烈振荡。这 里将采用隐式解法来求解对流项。如图2 - 1 所示，从每个网格单元逆 向追踪粒子轨迹，得到它前一时刻的位置，并通过线性差值来计算其 中的状态量。所以( 1 2 ) 式被描述成： ( x ) = w , ( x - u ( x ，t ) a t ) ( 1 3 ) 这样就能在较大时间步长内都能稳定求解，对流方程实际上可以 写成： 一o u ：0 ( 1 4 ) d t 然后对方程采用隐式解法来求解状态量，其具体可以描述如下：首先 可以想象在栅格中的一个点x 2 ，要得到它的速度值，记它的速度值为 帧士学位论文 甜g n + 1 ，在单位时间足够小的情况下，在缸时间里“；州= “；，从图1 可 以看出，计算方法就是计算栅格的前一时刻的速度。这里u 可以理解 为速度，温度等状态量。从已知公式拿：“出发，离散化为 d f = 毛一叱。其中是起点，是终点。即找到了起点，然后 求点的速度，速度的求法是使用与x p 相邻的点的速度进行线性插 值。 凡oo o卜o i ol o t 1 口( l d吨呸! 丝一 丑 。l o 旷 一q j b o 一o一ok o 图2 - 1 隐式求解对流示意图 接下来要计算扩散项，扩散项的p d e 为： 娑= 1 ，v 2 “( 1 5 ) d f 与求解平流一样需要使用一种稳定的解法，其隐式方程为： ( i - v a t v 2 ) w 3 ( x ) = w 2 ( x ) ( 1 6 ) 其中，是单位矩阵。 在计算完上述这些方程之后，得到的结果为一个散度非0 的新的 速度场，但是连续方程要求在每个时间段结束时得到散度为0 的速度， 利用h e l m h o l t z - h o d g e 胁3 分解知，矢量场w 能唯一的分解为如下形 式：w = 甜+ 跏，即一个无散度矢量场合一个标量场的梯度。使用 流体模拟的混合模型研究 h e l m h o l t z - h o d g e 分解定义一个投射算子p ，意义在于把一个矢量 场w 投射到它的无散度分量“上。对方程( 4 ) 两边使用投射算子，有： p 罢 = p ( 一( “v ) l l 一! 丫p + v v 2 u - i - f ) ( 1 7 ) 研p 因为“散度为0 ，左边的导数是无散度的，所以p 票：粤，并且 p ( 跏) = 0 得到： 譬= p ( 一( “v ) u + v v 2 甜+ ，) ( 1 8 ) 研 ( 1 8 ) 式的好处很明显，就是他把整个计算过程都封装到投射算子p 了。这就意味着我们在计算过程中不需要分部使用投射，只需要在计 算完所有状态量之后进行投射变换就可以得到无散度的速度场了。即 求解p o i s s o n 方程即可完成方程( 4 ) 的求解了。 2 3s p h 方法 s p h 方法于1 9 9 5 年被s t a m 等人首次引入计算机图形学中口1 ，用 来模拟火焰和其他气体的扩散现象。今年来s p h 方法异军突起，成为 新的c f d 计算方法，越来越多的学者开始研究s p h 方法的数值问题， 如常见的流体计算问题中的坝渍问题，目前用s p h 能取得满意的效果 和计算精度。以及一些学者提出s p h 方法的改进和修正方案，如：d ii t s 等提出的移动最小二乘粒子流体动力学方法眩5 2 引，目前s p h 方法可以 很好的模拟一般的流体力学问题。并逐渐成为基于粒子的拉格朗日方 法的主流方法。 如图2 - 2 ，为m i c k yk e l a g e r 于2 0 0 6 年采用s p h 方法实现的海浪 效果晗7 1 。 硕士学位论文 图2 - 2m i c k y 在2 0 0 6 年用s p h 实现的海浪效果 s p h 的核心思想归纳起来如下： ( 1 ) 将问题域以任意分布的粒子形式表示出来，粒子问不需要 连接。 ( 2 ) 用积分表示法来近似场函数，在s p h 方法中称为核近似法。 ( 3 ) 应用粒子来对核近似方程进步近似，也称为粒子近似。 该方法的过程是通过应用局部区域内的相邻粒子对应的值来叠加求 和取代场函数及其导数的积分。其中局部区域称为支持域。 ( 4 ) 在每一时间步内都要进行粒子近似，所使用的粒子取决于 当前局部分布的粒子。 ( 5 ) 将粒子近似法应用到所有的p d e s 的场函数中，得到一系列 只与时间相关的离散化形式的o d e s 。 ( 6 ) 使用显示积分法来求解o d e s 以获得最快的时间积分，并求 得与时间相关的状态量。 在s p h 方法中，函数的积分表达式的概念是由式( 1 9 ) 定 义的，其形式如下： ，( x ) = l f i x ) d 扛一x ) 如( 1 9 ) n 其中，为三维坐标向量z 的函数，d “一x ) 为狄拉克占函数，其性 质为： 流体模拟的混合模型研冗 万c x z ，) = 三三 c 2 0 ， 其中q 为包含x 的积分体积。只要在q 内( x ) 是已定义的和连续 的，则式( 1 9 ) 就是严格精确的。由于我们要使用离散的方法来进行 数值计算，所以不能直接使用狄拉克函数，所以需要使用光滑核函数 w ( x x ，五) 来取代狄拉克函数的核函数8 ( x 一工) ，于是厂( 工) 的积分表示式 为： 厂( x ) = l 厂( x 7 ) w ( x - x 7 ，矗) 出 ( 2 1 ) 矗 这就是积分表示的通式，为了方便反映我们关注的一些属性，这 里把流体中的状态量记为a ( 如速度，密度等等) ，则有： 彳( x ) = i a ( x ) w ( x - x 7 ， ) 出7 ( 2 2 ) 矗 将其积分形式转化为求和形式有： 彳( x ) = a s w ( x 一_ ，h ) a v j ( 2 3 ) 其中，是对所有粒子遍历，巧2 詈是粒子的体积，将_ 展开有： 砸) - z ，4 w ( 卜枷) 詈 ( 2 4 ) 这便是s p h 方法在描述粒子属性时采用的一般形式。由于在求解n - s 方程时会用到梯度算子和拉普拉斯算子，这里还需要得到状态量彳的 梯度计算形式和拉普拉斯计算形式。 从连续形式有： v 彳( x ) = 彳( x ) 矽 一，h ) d x ( 2 5 ) 矗 由于 i v a ( x ) w ( x - x ，j 1 1 ) = v 彳( x ) w ( x - x ，矗) 卜a ( x ) v w ( x - x ，j i l ) ( 2 6 ) 有： 硕士学位论文 v 彳( 工) = 么( x 7 ) 形( 工一工，h ) d x 一p ( x 7 ) v 矿( z z ，五) 出7 ( 2 7 ) 矗 矗 根据散度定理得到： v 4 ( x ) = p ( x ) 形( x 一工7 ，五) 褙一阻( 石) v w ( x 一工7 ，h ) d x ( 2 8 ) 矗矗 由于光滑函数是定义在紧支域内的，故当支持域位于问题域内，( 2 8 ) 式中的面积积分为零。如果支持域和问题域交错则光滑函数被边界截 断，面积积分不为零。当面积的积分不为零的时候就需要修正边界效 应的影响。这里就存在前面所说的边界弱效应。现在比较流行的解决 方法是使用虚粒子来加强边界，使得支持域始终处于问题域内部。这 样也就解决了离散求和的计算问题了，类似的对于拉普拉斯算子同样 成立。 于是有下式成立： v a ( 工) = z 4 m jv w ( x 一， ) ( 2 9 ) , , x j j rj v 2 a ( x ) = 4 _ m j v 2 一一，办) ( 3 0 ) j 卜l 解决了梯度和拉普拉斯算子的离散化之后，接着要解决光滑核函 数的选取，我们知道s p h 方法中光滑核函数( 光滑核函数也称为滤波 函数) 起着很重要的作用，因为它决定了函数表达式的精度和计算效 率，同时我们也可以根据需要构造自己的核函数，下面介绍一下核函 数需要满足那些条件： ( 1 ) 光滑核函数在支持域上要满足归一化条件，即： i w ( x - - x 7 ，h ) a x = 1 ( 3 1 ) q ( 2 ) 光滑核函数必须满足紧支性条件，即： w ( x x ) = o ，lx x l x h ( 3 2 ) 其中h 是光滑长度，鬈是比例因子，他们决定了紧支域的尺度。 ( 3 ) 在点x _ l 的粒子的支持域内任意一点处有： 煎堡塑垫盟塑鱼堡型堕塞 z ( x 门0 ( 非负性) ( 3 3 ) ( 4 ) 当粒子间的距离增加时，粒子的光滑函数值是单调递减的。 ( 5 ) 当光滑长度趋于0 时，光滑函数应该满足狄拉克j 函数条 件： 1 岫w (