（计算数学专业论文）基于谱方法的实时流体模拟.pdf

摘要 用计算机真实地再现水、云、烟为代表的流体在自然界中的运 动过程和运动规律，即流体模拟是计算机图形学中的一个重要研究 领域。流体模拟技术广泛应用于流体力学研究、航空航天、气象研 究、海洋工程、汽车设计以及计算机动画等众多领域，是有关基础 科学理论研究和应用工程设计研究的重要手段。因此，流体模拟的 研究在理论研究和工程实践上都有着重要的意义。 本文从计算机图形学的角度，依据流体力学的基本理论，用谱 方法模拟不可压均匀低速流体的运动规律。基于谱方法的流体模拟 是有别于目前通用的半拉格朗日稳定流体模拟方法的另一种尝试。 该方法将流体运动从通常的实空间转换到谱空间进行运算，对纳维一 斯托克斯方程采用谱配置法进行求解，从而直接消去了压力项，并 将原方程中的扩散、平流、外力项做为一个整体用简单的代数运算 来进行处理，避免了把纳维一斯托克斯方程拆成几项分别求解带来的 误差和用迭代法求解泊松方程的繁复运算，大大降低了实时流体模 拟过程中的运算量。在计算机的实现上，本方法充分利用图形处理 器( g p u ) 强大的并行运算能力，采用渲染到纹理的计算方式，帧缓存 扩展的纹理管理方法，实现了对流体的实时模拟。通过实验比对可 以看到，基于谱方法的流体模拟在模拟的实时性和流体细节精度等 方面较半拉格朗日稳定流体模拟方法都有着不同程度的提高。 关键词：实时流体模拟，纳维一斯托克斯方程，谱方法，图形处理器 a b s t r a c t b u i i d i n ga n i m a t i o nt o o l sf o rf l u i dm o t i o n sb yc o m p u t e ro n m a n yn a t u r a lp h e n o m e n as u c ha sw a t e r ，c l o u da n ds m o k e ，w h ic h iso n eo ft h em o s ti m p o r t a n tp r o b l e m si nc o m p u t e rg r a p h i c s f l u i ds i m u l a t i o nh a sl o t so fa p p l i c a t i o n si nm a n ya r e a ss u c h a sf l u i dm e c h a n i c s ，a e r o n u t i c s ，a s t r o n a u t i c s ，m e t e o r o l o g y ， o c e a ne n g i n e e r i n g ，a u t od e s i g n ，c o m p u t e ra n i m a t i o na n ds oo n s o ，t h er e s e a r c h e so ff l u i ds i m u l a t i o na r ev e r ys i g n i f i c a n t i nt h e o r yr e s e a r c h e sa n de n g i n e e r i n gp r a c ti a c e i nt h i sp a p e r ，a c c o r d i n gt ot h eb a s i ck n o w l e d g eo nf l u i d m e c h a n i c s ，w eu s es p e c t r a lm e t h o d st os i m u l a t et h ei n c o m p r e s s i b l ef l u i dw h o s ed e n s i t ya n dt e m p e r a t u r ea r en e a r l yc o n s t a n t i t i sad i f f e r e n tw a yf r o mt h es t a b l ef l u i d ss i m u l a t i o n w e s o l v et h i sp r o b l e mi ns p e c t r a ls p a c ew h i c hi sd i f f e r e n tf r o m t h eu s u a lm e t h o di nr e a ls p a c e b yu s i n gs p e c t r a lc o l l e c t i o n m e t h o d so nt h en a v i e r s t o k e se q u a t i o n ，w ec a no m i tt h ep r e s s u r e i t e ma n dc o n s i d e rt h ed i f f u s i o n ，a d v e c t i o na n da d d i t i o n a l f o r c ea sa ni n t e g e rt od e a lw i t h t h i si sas i m p l ea l g e b r a o p e r a t i o na n di t i sv e r ye a s yt oi m p l e m e n t b yt h i sw a y ，w e n e e d td e c o m p o s et h en a v i e r s t o k e se q u a ti o na n ds o l v et h e p o i s s i o ne q u a t i o n i tr e d u c et h ec a l c u l a t i o n so nr e a l t i m ef l u i d si m u l a ti o n i nc o m p u t e ri m p le m e n t a t i o n ， w ea ls ou s et h ew a y i i o fr e n d e r i n gt ot e x t u r ea n df r a m e b u f f e re x t e n s i o n w em a n a g e t op u ta l lt h et e x t u r e sn e e d e di ng r a p h i ch a r d w a r em e m o r yw h i c h i n c r e a s e st h er a t eo fu t i l i z a t i o no fg p u t h et e s tr e s u l t ss h o w t h a to u rm e t h o dc a ni m p r o v et h ea c c u r a c ya n de f f i c i e n c y c o m p a r e dw i t ht h es t a b l ef l u i d s k e yw o r d s ：r e a l t i m ef l u i ds i m u l a t i o n ，n a v i e r s t o k e se q u a t i o n ， s p e c t r a lm e t h o d ，g p u ( g r a p h i c sp r o c e s s i n gu n i t ) i 湖南师范大学学位论文原创性声明 本人郑重声明：所呈交的学位论文，是本人在导师的指导下，独 立进行研究工作所取得的成果除文中已经注明引用的内容外，本论 文不合任何其他个人或集体已经发表或撰写过的作品成果。对本文的 研究做出重要贡献的个人和集体，均已在文中以明确方式标明。本人 完全意识到本声明的法律结果由本人承担。 学位论文作者签名：陇弁怨 一。7 年月中日 湖南师范大学学位论文版权使用授权书 本学位论文作者完全了解学校有关保留、使用学位论文的规定， 研究生在校攻读学位期间论文工作的知识产权单位属湖南师范大学 同意学校保留并向国家有关部门或机构送交论文的复印件和电子版， 允许论文被查阅和借阅。本人授权湖南师范大学可以将本学位论文的 全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或扫 描等复制手段保存和汇编本学位论文。 本学位论文属于 1 、保密影在年解密后适用本授权书。 2 、不保密瓯 ( 请在以上相应方框内打“ ) 作者签名： 发辱碍日期：p d 陴6 月够日 导师签名： 乞幻侈价 日期： 呻年 6 7 月v 日“ 基于谱方法的实时流体模拟 1 1 引言 1 绪论 随着交互式计算机图形学的不断发展和计算机硬件的不断更 新，越来越多的领域如计算机辅助设计、科学仿真计算、电视广告、 计算机动画等方面，计算机图形学都起到了至关重要的作用。 在计算机图形学中，研究人员总是努力试图利用计算机来再现 周围的真实环境，其中以水、云、烟为代表的流体在自然界中的运 动过程和运动规律，即流体模拟已成为当前的一个研究热点、新兴 领域，也是未来计算机图形学发展的一个重要方向。 流体模拟在当代社会有着非常广泛的应用，如在流体力学的基 础研究中，时常需要形象地展示流体的运动形态便于研究者对其进 行观察分析；在航空航天领域常常需要事先模拟航天器在大气中飞 行时周围的气流环境，以检测某些大气条件对航天器飞行的影响； 在汽车外形设计中流体模拟研究已占主导地位，合理的外形能使汽 车具有更好的动力学性能和更低的耗油率；在建筑领域也常常要对 暴雨、台风等恶劣气候条件进行模拟，从电脑模型上检验建筑物抵 御这些恶劣灾害环境的能力；在计算机动画、电影特技制作中云彩 的变化、河水的流动、大气的扩散、爆炸等大量特效场面都在不同 程度上运用了流体模拟技术；甚至在气象、海洋工程等方面都会涉 及到流体模拟的应用。流体模拟不仅是广大计算机图形学研究者关 硕士学位论文 注的焦点，更是相关基础科学理论研究和应用工程设计研究的重要 手段。因此，流体模拟的研究在理论研究和工程实践上都有着重要 的意义。 计算机图形学中的流体模拟有别于流体工程计算和流体理论研 究，本文提及的实时流体模拟都是基于计算机图形学的角度即人眼 观察角度，将不可压均匀低速流体做为研究对象对其运动进行模拟， 该运动过程并不需要达到流体工程计算或流体理论研究所要求的精 度和准确度，只需流畅、逼真的表现流体的运动形态即可。一个有 效的实时流体模拟( r e a l t i m ef l u i ds i m u l a t i o n ) 通常具有两个特 征：一是模拟流体运动形态的逼真性，二是模拟流体运动形态的实 时性，而这两个特征从某种意义上来说又是相互矛盾的。要提高流 体运动形态的逼真度势必提高流体计算的精度，就会增加计算量， 导致实时性的降低；另一方面，要增强流体模拟的实时性往往会忽 视流体运动形态上的某些细节丰富程度，从而导致模拟逼真度的下 降。因此，如何又快又好地实现实时流体模拟既是实时流体模拟的 两个重要标准又是实时流体模拟研究中的一项难题。 1 2 实时流体模拟的研究发展及现状 早期的流体模拟由于计算机的计算能力有限，往往采用参数建模 的方法。文献 1 ， 2 根据波动理论用三角函数来模拟波浪的效果， 但是该方法所表示水的粒子只是在其起初始位置附近运动，无法真实 地表现流体的流动效果，并且这类模型不能模拟复杂的、细节更为丰 基于谱方法的实时流体模拟 富的流体效果，因此，大量的图形学工作者转向基于物理模型的流体 模拟研究。 流体力学( f l u i dd y n a m i c s ) 中常用一些复杂的数学物理方程 来精确描述流体的运动规律，其中最基本、最常用的是纳维一斯托克 斯方程组( n a v ie r s t o k e se q u a ti o n s ，简称n s e s ) ，它是描述不可 压缩流体的物理模型。文献 3 采用二维n s e s 求解流体表面速度场， 然后对流体表面压强伯努利方程求解高度场，模拟运动物体或其它 障碍物对流体运动的影响，但效果不理想；真正用三维n s e s 来模拟 流体运动是始自文献 4 用m a c ( m a r k e ra n dc e l l s ) 求解n s e s ，该 方法采用显式格式，必须满足极小的时间步长才能使整个计算收敛， 而在较大时间步长时则容易导致n s e s 求解产生很严重的耗散，使得 该方法只适用于较小时间步长，计算量大，实时性不高；直到1 9 9 9 年，s t a r e 5 在计算机图形领域中采用半拉格朗日稳定流体方法，才 使得实时流体模拟成为可能。该方法允许采用较大的时间步长，并 且能够较好地获得流体运动的基本特征，较之前的各种流体模拟方 法在逼真性和实时性方面都有大幅度的提高，奠定了现代计算机图 形学实时流体模拟的基础，但是半拉格朗日稳定流体方法在数值耗 散方面还是比较严重的，部分流体的旋涡细节被忽视了；f e d k i w 在 文献 6 中引入旋涡约束因子( v o r t i c i t yc o n f i n e m e n t ) 以弥补半拉 格朗日稳定流体方法带来的数值耗散，恢复由高旋度区域指向低旋 度区域的向心力，产生一个局部区域中的旋涡，从而更加逼真地模 拟了烟雾；文献 7 中引入热力学参量参与流体方程的求解，加入了 硕士学位论文 浮力和温度对n s e s 的影响，成功地模拟了云的生成和飘动。 随着g p u 在运算性能上和可编程度上的飞速提高，原先在c p u 上实 现的研究成果有相当一部分被移植倒g p u 上实现。从2 0 0 4 年开始，利 用g p u 进行通用计算的研究( g e n e r a lp u r p o s eo ng p u ，简称g p g p u ) 8 得到开展和普及，一些基于g p g p u 的线性代数和偏微分方程的数值解 法得以实现，这些方法被证明比在c p u 上实现有成倍提高。其中比较 著名的包括：文献 9 利用g p u 实现了大型稀疏线性方程组的高斯一赛 德耳迭代( 简称g - s 迭代) 法和共轭梯度法；文献 1 0 利用g p u 进行基于 有限元方法的图像处理；文献 11 在g p u 上应用多重网格算法求解边 界值问题；文献 1 2 实现了基于像素编程的稀疏非结构化矩阵的共轭 梯度法和多重网格法。 g p u 在数值计算方面的研究发展，使得大量的计算机图形学研究 者开始利用g p u 的高速计算和渲染技术来进行流体模拟方面的研究， 提高模拟的实时性。文献 7 ， 1 3 把流体模型和热力学模型，以及水 的蒸发和冷凝的模拟结合起来在g p u 上实现了云的模拟；文献 1 4 在 g p u 上利用四层网格依靠渲染到纹理的计算方式实现了多重网格法， 改进了二维实时流体模拟。 目前，计算机图形学中的实时流体模拟大多都是采用半拉格朗日 稳定流体模拟方法，并在其基础上采取某些细节优化策略，弥补半拉 格朗日方法产生的数值耗散，提高模拟的精度和真实性，或优化g p u 加速策略，使流体模拟实时性不断提高，力求计算机实时流体模拟达 到逼真、自然的效果，在满足人们的视觉效果的同时实现实时交互的 基于谱方法的实时流体模拟 功能。 近几年，实时流体模拟在对上述方面进行优化研究的同时，还逐 渐形成了两个新的研究发展方向：一类是使流体照指定的规则进行流 动形变，以控制流体的运动形态。文献 1 5 3 通过对外力场的控制和改 进扩散方程进而形成了烟雾聚集的初步理论控制模型；文献 16 采用 类似的方法并加入对烟雾的密度控制实现了烟雾的聚集；文献 1 7 采用关键帧技术将球状的烟雾逐步形变成一匹马的形态。另一类是在 复杂边界上或较复杂几何体表面实现流体模拟。文献 1 8 在g p u 上采 用离散格子玻尔兹曼方法( l b m ) 实现不同障碍物的流场模拟；文献 1 9 中在g p u 上利用模板缓冲对三维场景进行剪切造作，将三维计算 域离散成一组切片，并将这组切片平铺为一幅二维纹理，实现了带有 复杂边界的三维实时流体模拟。文献e 2 0 ， 2 1 则对复杂的模拟区域 采用三角形或四面体剖分，在剖分网格表面实现了流体模拟，提高了 模拟的场景复杂度。 1 3 主要内容及基本框架 本文在对目前常用的半拉格朗日稳定流体模拟方法进行研究基 础之上，针对其不足，提出了基于谱方法的实时流体模拟并在g p u 上 实现其应用。 本文的主要研究内容是： ( 1 ) 了解流体模拟的发展及现状。 ( 2 ) 立足于流体理论知识和流体物理模型n s e s 的研究，分析基 硕士学位论文 于半拉格朗日稳定流体模拟方法的优点与不足。 ( 3 ) 提出一种基于谱方法的实时流体模拟方法。 ( 4 ) 研究谱方法在g p u 上的实现问题。 根据研究内容，本文的基本框架如下： ( 1 ) 第一章是绪论，简要介绍流体模拟的研究发展过程，流体 模拟的现状和本文的研究内容、创新点。 ( 2 ) 第二章简要介绍流体基本理论和流体模型，并对半拉格朗 日稳定流体方法进行较为全面的分析。 ( 3 ) 第三章提出一种基于谱方法的实时流体模拟方法，并在g p u 上实现该方法的应用，最后给出实验结果并与半拉格朗日稳定流体 方法和加入旋涡约束因子的稳定流体方法进行实验分析比较。 ( 4 ) 第四章是结束语，对全文进行一个总结和回顾以及进一步 工作展望。 1 4 本文的创新点 针对目前计算机图形学流体模拟中普遍采用的半拉格朗日稳定 流体模拟方法的不足，本文将计算数学中的谱方法思想引入到计算 机图形学流体模拟的算法设计中，提出了基于谱方法的实时流体模 拟方法，用谱方法的思想处理流体模拟问题是一种大胆地尝试。它 可能带来的好处有：( 1 ) 谱方法是一种对偏微分方程进行高精度求 解的算法。用谱方法处理流体问题不仅有可能提高流体模拟的精度， 也丰富了计算机图形学流体模拟的理论和方法。( 2 ) 本方法将流体 基于谱方法的实时流体模拟 运动求解从通常的实空间转换到谱空间下进行，对n s e s 用谱配置法 在谱空间下消去了n s e s 中的压力项，避免了在实空间下对n s e s 拆 分成几个偏微方程带来的误差以及拆分后对压力进行繁杂的迭代求 解时所产生的数值耗散。( 3 ) 本方法的流体模拟求解的实现过程将 n s e s 中的扩散、平流、外力项做为一个整体，对实空间与谱空间下 的速度场进行f f l 运算后，即可用简单的代数运算完成n s e s 求解， 避免了逐项计算扩散、平流、外力项的繁杂过程，简化了流体模型 n s e s 的求解，大大降低了实时流体运动模拟过程中的运算量。( 4 ) 在计算机的实现上，充分利用图形处理器( g p u ) 的强大的并行运算 能力，采用渲染到纹理的计算方式，帧缓存扩展的纹理管理方法，实 现了对流体的实时模拟，提高了流体模拟效率。 硕士学位论文 流体模拟的基本理论 为了描述流体的运动规律，流体力学领域的众多研究者已经做 了大量的工作，形成了一整套较为完善的理论体系。计算机图形学 目前主要是针对不可压均匀低速流体运动现象，研究不可压流体的 物理模型纳维一斯托克斯方程组( n a v i e r - s t o k e se q u a t i o n s ，简称 n s e s ) 完成一些简单流体，如水、云、烟的实时模拟。本章首先简 单介绍流体力学的一些基本理论和流体的物理模型n s e s ，最后详细 阐述半拉格朗日稳定流体模拟方法的基本思路，并分析该方法的优 点与不足。 2 1 有关流体的基础知识 流体就是受到任何微小剪切力的作用都能够导致其连续变形的 物质，它是由大量的、不断地作热运动而且无固定平衡位置的分 子构成的。简单来说，流体就是可以流动的物质。在现实生活中， 随处都有运动着的流体的影子。如沿岸流淌的河水、升腾成云的水 汽、香烟燃烧扩散的青烟等等，所有这些全部都是运动的流体。 2 1 1 流体的基本特征 流体通常分为液体和气体，具有易变形性、粘性和压缩性三 个基本特征，其中易变形性是流体的决定性特征。 1 ) 易变形性 基于谱方法的实时流体模拟 流体与固体的分子结构存在差异，流体具有易变形性，它 没有固定的形状且具有流动性。液体分子间具有凝聚力，有一定 的体积，但没有确定的形状。注入任何容器的液体，都将在容器 内填满一个与液体自身体积相同的空间；气体则是由大量运动着 的分子组成，分子之间相互碰撞，不断扩散。气体没有确定的体 积和形状，且会充满它所注入的任何容器。 流体的易变形性决定了流体与固体有着许多不同的运动特 征。当受到剪切力持续作用时，固体的变形是有限的，流体能产 生无限大变形或流动；固体内的切应力由变形量决定，流体内的 切应力由变形速率决定；当剪切力停止作用时，固体形变能恢复 或者部分恢复，而流体则不做任何恢复。 2 ) 粘性 流体在运动过程中，由流体内部摩擦力产生阻碍其本身运动 的性质称为流体的粘性。如酒精和蜂蜜这两种不同的液体，我们 会明显地发现蜂蜜比酒精更难改变自身的运动状态，也就是通常 所说的蜂蜜比酒精粘稠，即粘度是运动流体对剪切力反抗的度 量。用剪切力f 即单位接触面积上的力表示流体内摩擦力。实验 表明，剪切应力r 的大小同速度梯度华成正比，这就是著名的牛 砂 顿粘性定律，即：f = 华，其中表示粘度或粘滞系数( v is c o s it y 掣 o ra b s o l u t ev is c o s i t y ) ，它反映了流体粘度的动力性质，所以 又称动力粘度：与动力粘度对应的是运动粘度) ( k i n e m a t ic 硕士学位论文 v is o s i t y ) ，二者的关系是u ：丝，d 表示在单位速度梯度作用下 p 切应力对单位体积质量作用产生的阻力加速度。 流体的粘性与温度有着密切的关系，液体的粘度随着温度升 高而减小，气体则相反。因为，液体分子之间的内聚力对液体的 粘性起主要作用，当温度升高时，分子间距变大，内聚力相应变 小，因而粘度下降。而对气体而言，相邻流层之间分子动量的交 换对气体粘性起主要作用，当温度升高时，气体的热运动加强， 动量交换增多，各层之间的制动加大，因而粘度增大。 3 ) 压缩性 压缩性就是在一定温度下，压强变化引起流体体积或密度变 化的性质。液体的可压缩性很小，其密度随温度和压强的变化也 很小。而气体的密度随着温度和压强的改变将发生显著的变化， 因此气体具有明显的可压缩性。 当流体的粘性和可压缩性都很小时，即为无粘性不可压缩 时，该流体可近似看作是理想流体，它是为研究流体的运动和状 态而引入的一个理想模型。 2 1 2 流体模型的分类 通常，根据流体的粘性大小可将其分为无粘性流体与粘性流 体，根据流体是否压缩分为可压缩流体和不可压缩流体。 1 ) 无粘性流体与粘性流体 无粘性流体：由于空气和水的粘度很小，当流层间的速度梯 - l o - 基丁谱方法的实时“【体模拟 度不大时，流体的粘性应切力可忽略不计可建立无粘性流体模 型( u = 0 ) 。建立在无粘性流体模型基础l 的伯努利方程、环量 理论和表面波坪论等，在解释水和卒气流动中的机械能守恒、机 翼升力和水波运动等方面取得了成功，但无粘性流体模拟在解释 管道和渠道流动压强损失及绕流物体阻力方面却无能为力。 粘性流体：同无粘性流体模型相比，粘性流体模型更接近真 实的流体。粘性流体流动存在层流和湍流两种形态。层流足指流 体的流动是分层的或分片的。湍流则与之相反，在湍流中，速度 分量在其平均值上还叠加有随机的湍流脉动( 见图2 - 1 ) 。 图2 1湍流 对于特定的流体可根据流体流动的速度和通道的形状或大 小来确定该流体的流动是层流还是湍流。点燃的香烟升起的烟 柱，这就是由层流到湍流过渡的例子，在开始的一段时间里，烟 以平缓的层流方式升起，然后是突然地烟柱开始混合，从而转变 成湍流。 2 ) 可压缩流体与不| j 压缩流体 一般情况f ，液体的压缩性可以忽略，建立不町压流体模型 ( p = 常数) 。常温常抵气体作低速流动时，气体密度的相对变化 硕士学位论文 较小，可以按不可压缩流体处理。当气体作高速流动时，考虑到 其密度变化带来的影响，把它看成可压缩流体。 牛麓毓俸p l - 0 撇) 薯阡籁流俸( p 常数) 正珏流体p - p0 p ) 耕压流俸( 9 - p ( p , t ， 图2 - 2 流体的大致分类 流体模型的大致分类可以参见图2 2 。当然，还可根据不同 需要建立不同类型的流体模型给分析带来方便，如重力流体、均 质流体、正压流体、斜压流体等等。 2 1 3 流体的连续介质假设与流体运动的描述方法 刚体力学中，常考虑质点的空间位置如何随时间变化。从这 个信息出发，考虑质点的空间位失，速度，加速度。然而，流体 力学并不是处理单独的质点，而是依据欧拉提出的流体连续介质 假设，研究连续介质。 流体的连续介质假设就是不考虑分子间存在的间隙，而把流 体看作是由无数连续分布的流体微团或质点组成的连续介质。连 续介质模型的假设使表征流体属性的物理量，如密度、压力( 压 强) 、温度、速度、粘度、应力等都可以看成是连续分布的，可 以用时间、空间的单值连续函数来表示，从而就可以用微分方程 来描述流体的平衡、运动规律了。流体连续介质假设是对实际流 1 2 基于谱方法的实时流体模拟 体的数学抽象，就像几何学对自然物图形的抽象一样。只有在稀 薄气体和激波层内等少数情况下，连续介质假设才不适用。 事实上，我们不必留意流体中各个质点或某些小团流体的运 动轨迹，而是关心在空间中某个相对于任何一个固定坐标系的位 置上，流体的速度、加速度和热力学参量随时间变化的关系。在 时间过程中，空间点上的流体不断地为流过来的流体进行更新。 因此，不考虑各个单独流体质点的轨迹，也不管某个特定时间是 哪些流体质点在这个位置上，只关心空间某些点上的历史过程， 这种描述流体的方法称为欧拉描述方法。与之相反的则是拉格朗 日描述方法，它像刚体运动力学中那样，记录各个流体质点的轨 迹。 2 2 n a vie r - s t o k e s 方程组 在科学理论研究中，对复杂现象构建理论模型时要对该现象进 行一定的简化，流体模拟也不例外。如果某流体的任何局部区域的 体积不随时间而改变，则该流体为不可压缩的；如果该流体的密度p 在空间中是常数，则流体是均匀的。这是流体力学中普遍的假设， 也适用于水和空气等真实流体的模拟。本文所研究的流体是不可压 缩均匀低速流体。 在其它热力学参量确定的情况下，流体模拟的关键就是正确地、 及时地确定当前流体的速度场。速度不仅决定流体本身的运动形态， 而且还决定了掺杂在流体中的其他物质如何运动。流体的速度既随 硕士学位论文 时间又随空间而改变，因此把它表示为一个矢量场，矢量场是矢量 函数在一个参数化空间的映射。流体的速度矢量场是对于每个位置 工= ( 工，y ) g f f 定义的，也存在与时间t 相关的速度，a ( x ，f ) = ( 甜( 工，f ) ，v ( x ，f ) ) 。 用空间坐标工= ( 工，y ) 和时间变量t ，在一般笛卡尔坐标系上模拟 流体，流体由它的速度矢量场a ( x ，f ) 和标量压力场p ( x ，f ) 表示，这些 场既随时间又随空间而改变。对不可压缩流体，其初始速度和初始 压力在t ：0 时已知，那么流体随时间的状态变化可以用n a v i e r - s t o k e s 方程组( n s e s ) 来描述，这是研究不可压流体的基本物理模 型。 整个n s e s 主要由两个部分构成：一个是动量方程，另一个是不 可压条件。 詈邓v 矿吉即+ u v 2 f f + 于( 2 - 1 ) v 厅：0( 2 - 2 ) 其中，厅表示流体的速度矢量，p 表示压强，歹表示外部施加的作用 力，d 表示运动粘性系数( 运动粘度) ，p 表示流体密度，u 和p 均为 常数。v 表示矢量的梯度算子，如在二维空间中v ：( 昙，晏) ，在三维 甜卯 空间中v 表示为( 昙，晏，昙) 。v 是散度算子，在二维空间下v 一昙+ 昙。 廊卯0 2o x 卯 v 2 是拉普拉斯算子，v 2 ：v v 。 n s e s 中，露代表二维或三维空间矢量，因此方程( 2 1 ) 是由具 有二维或三维分量的速度场构成的空间矢量方程。它是牛顿第二定 律户：所五应用于流体运动中的特殊形式，也称( 2 - 1 ) 式为动量方程。 基于谱方法的实时流体模拟 考虑到研究的流体是不可压缩的，即体积不变，用积分方法可以得 到流体的不可压缩条件，即( 2 - 2 ) 式，用于保证流体速度场厅的散 度为零。本文对流体物理模型n s e s 不予进行推导和证明，它是根据 流体力学力学中的四大基本规律( 质量守恒、牛顿第二运动定理、 能量守恒定律、热力学第二定理) 通过微积分方法得到的。 流体物理模型n s e s 看起来很简单，但是由于方程( 2 - 1 ) 等号右 边存在非线性项，使得方程求解非常困难。事实上，n s e s 求解的难题 代表了流体模拟乃至流体力学研究中遇到的主要挑战之一。我们逐项 观察n s e s ，了解该方程中影响流体运动的因素，更好的理解流体的物 理模型。 1 ) 平流项：( 2 - 1 ) 式右端第一项表示流体的速度场自身的传送， 称之为平流项。流体具有速度使得流体沿着流动方向传送物质、密度 等其他的量，可以将其想像成将染料注入流动的液体中，染料沿着液 体的速度场被传送。事实上，正如液体的速度场能够携带染料一样， 它还携带着液体本身。 2 ) 压力项：( 2 1 ) 式右端第二项称为压力项。流体在运动过程 中，其流体分子也在流动中相互“挤压 和“碰撞”，当外力作用在 流体上时，外力并不是立即地传递到整个流体，而是由离外力较近的 流体分子将该流体分子周围离外力相对远一点的流体分子推开，流体 内部产生了压力。 3 ) 扩散项：( 2 - 1 ) 式右端第三项是扩散项，该项表明粘度对流 体扩散阻碍的衡量。观察各种液体发现，有些液体比其他液体更“粘 硕士学位论文 稠 ，例如蜂蜜流动很慢而医用酒精却流动的很快。 ( 4 ) 外力项：( 2 - 1 ) 式右端第四项是外部施加到流体上的作用 力，这些外力可以是整体施加于流体的力，也可以是局部施加在流体 特定区域的力。 本文在上述两节中，大致介绍了流体力学的一些基础理论和流体 的物理模型n s e s ，这是进行流体实时模拟的基本理论依据，更多 有关流体力学基本知识和流体物理模型的详细内容可以参考相关书 籍 2 2 ， 2 3 ， 2 4 ， 2 5 。 2 3 基于半拉格朗日稳定流体模拟方法 实时流体模拟就是在计算机上实现流体运动模型求解，并以动 画的形式将流体的运动形态形象实时地表现出来。目前，计算机图 形学中实时流体模拟大多采用s i a m 提出的半拉格朗日稳定流体模拟 方法 5 ，在其基础加入旋涡约束因子，或改进迭代格式提高求解精 度等方法增强流体的某些细节，以实现流体的实时模拟。因此，研 究半拉格朗日稳定流体模拟方法，对掌握目前计算机图形学实时流 体模拟有着重要的理论和实际意义。 首先，s t a r e 引入了h e l m h o l t z h o d g e 分解，大大简化了n s e s 的求解。 2 3 1 h e i m h o l t z - h o d g e 分解 h e l m h o l t z h o d g e 分解定理：区域d 上的任意一个矢量场访可以 基于谱方法的实时流体模拟 唯一的被分解为一个零散度矢量场订与一个标量场p 的梯度的和。 访= 厅+ v p ( 2 3 ) 这里面是无散度矢量，即v 厅= 0 ，且与劬平行，即在d 上露，l = 0 ， 这是一个简化n s e s 的有效方法。 1 ) 在每个时间段内通过逐次计算平流、扩散、外力三项后可得 到一个散度非零的新速度场形，但n s e s ( 2 - 2 ) 式要求每个时间段结 束后得到的速度场的散度为零。根据该分解定理，只需对话减去压 力场的梯度即可得到散度为零的速度场历= 访一跏。 2 ) 该定理同时引出了计算压力场的方法。对( 2 - 3 ) 式两边同 时作用散度算子，得到流体压力的p o i s s o n 方程。 v 访= v 2 p ( 2 4 ) 这也就是说在得到散度速度力之后，求解方程( 2 - 4 ) 得到压力p ， 最后将帚和p 带入方程( 2 3 ) 即可得到新的零散度速度场露。 2 3 2 投影算子p 根据h e l m h o l t z h o d g e 分解，可定义一个投影算子p ，它将一个 矢量场茹投射到它的零散度分量厅上。将尸作用于( 2 - 1 ) 式两端， 由于筋= 厅，p v p = 0 则有： 粤= p ( 一( 厅v 够+ u v 2 订+ 歹) ( 2 5 ) 优 这样，投影算子p 就封装了流体模拟的整个算法。首先计算 ( 2 - 5 ) 式右端括号里面的平流、扩散和外力项，得到速度场谚，然 后对它应用投影算子p 得到一个新的零散度的速度场厅。因此整个 硕士学位论文 n s e s 求解过程可以表示为： s ( u ) = p d a f ( 露) ( 2 6 ) 其中各分步算子定义如下：投影p 、扩散d 、平流彳和外力f ，每一 个算子都是一个步骤，它取一个场作为输入，产生新的场作为输出， 在实现时依次从右向左进行运算，如图2 - 3 所示。 外力f乎漉a扩散d投影p w o ( x ) 叫w l ( x ) _ w 2 ( x ) 一w 3 ( x ) 一w 4 ( x ) 图2 3 半拉格朗日稳定流体模拟方法的n s e s 求解步骤 2 3 3 平流、扩散项的半拉格朗日处理方法 计算平流量，就必须更新每个坐标点上的量。假设每个网格单 元由一个粒子代表，每个粒子在所时间段内沿着速度场向前移动的 距离为，- ，由e u l e r 方法显然有厂( f + 田) = ，( f ) + a ( t ) s t ，但在较大的 时间步长下该表示法对n s e s 求解是不稳定的。故s t a m 提出了半拉 格朗日方法，它不是通过计算粒子在当前时间移动到什么位置来得 到平流量，而是对每个流体粒子逆时追踪该粒子的运动轨迹，得到 它在万f 时刻之前的位置，并复制该位置的物理量作为其新的物理 量，半拉格朗日方法对于任意时间步长都是稳定的并且容易实现。 若更新的物理量为q ( 它可能是速度、密度、温度或被液体携带的任 何量) 使用方程： q ( 工，t + a t ) = q ( 工一a ( x ，t ) s t ，f ) ( 2 7 ) 来隐式地表平流量，图2 - 4 描绘了如何实现平流量的计算。隐式的 基r 谱方法的蛮时流体模拟 、f 流量计算是对速度场进行逆向追踪，以确定是如何向前携带量。 为了计算以双固标记的单元的平流，逆时追踪速度场到绿色的x ，该 t 意黄的量就是救固位置的平流量。而点物理量的汁算可对最靠进 点的4 个网格值进行双线性插值得到。 、+ 一 罔2 - 4 流体甲流量的半拉格朗日算法 对于扩散项，考虑粘性扩散的偏微方程为娑= o vz i ，与平流项类 叭 似地使用显示方法计算也会导致对较大时问步琵的求解不稳定。 s t a r a 再次采用半拉格朗日方法对粘性扩散方程采用隐式表达： ( ，一o s t v2 ) i ( ，r + 曲) = i ( z ，f )( 2 - 8 ) 其中，是单位矩阵。这样( 28 ) 式对于任何时问步长都是稳定的。 接f 来，只需完成两个p o i s s o n 方程的求解，一个是压力方程 ( 2 - 4 ) ，另一个是扩散方程( 2 - 8 ) 。p o i s s o r l 方程是形如a ，= b 的矩 阵方程式，这里x 足所需求解的数值矢量( p 或i ) ，6 是一个矢量常 数，a 是一个矩阵，采用j a c o b i 迭代法进行求解。依照图2 - 3 的流 程完成整个的n s e s 的求解。然后，考虑流体边界条件以及流场的可 视化实现、仁托格朗口稳定流体模拟方法。 箱 硕士学位论文 研究发现烟雾、空气或其它低粘度流体运动时，会引起一些大 小不一的旋流，这种旋流是旋涡状态的。半拉格朗日稳定流体模拟 方法由于某种原因在求解n s e s 时产生比较严重的数值耗散，减弱了 局部旋涡的细节。因此，f e d k i w 6 在半拉格朗日稳定流体模拟方法 的基础上提出在模拟中人工地加入旋涡状态，用来恢复这些细节的 运动。首先，计算旋涡缈；v ”，其中“旋度 算子vx ，可以通过 梯度的定义和叉积算子推导出来。然后规范旋涡矢量场少：鲁，这里 吲 ，7 = v l 叫，该矢量场中的矢量，从低旋涡状态区域指向较高旋涡状态 区域，人工计算一个虚拟的力厶= 占舻缈) 面，用它近似恢复耗散掉 的旋涡特征，以增强半拉格朗日稳定流体模拟中的旋涡细节。 目前，在计算机图形学界已广泛采用基于半拉格朗日稳定流体 模拟方法。 2 4 半拉格朗日稳定流体模拟方法的不足 通过对半拉格朗日稳定流体模拟方法的研究，我们发现该方法 较早期的参数建模方法在模拟的真实性方面有着显著的提高，并且 解决了c f d 中较大的时间步长下数值耗散的问题，解决了模拟的实 时性。因此目前计算机图形学实时流体模拟通常都是采用这种稳定 流体技术，但是该方法仍然存在着几点不足。 首先，从算法上看，半拉格朗日稳定流体模拟方法是通过将n s e s 方程拆分为几个偏微分方程，用迭代法求解，导出更新速度场，完 成n s e s 的求解。这种拆分的方法本身就不可避免的带来了误差。 2 n 基于谱方法的实时流体模拟 其次，该方法在求解n s e s 方程中并没有彻底解决掉压力项部分 的难题，并且在对拆分后的扩散方程与压力方程都是采用j a c o b i 迭 代法求解的。这种迭代求解本身也是一个相对较为耗时且精度不高 的过程。因此，大量的研究者针对这点进行迭代格式的改进以提高求 解精度，或通过加入旋涡约束因子 6 的方法来弥补迭代产生的数值 耗散。 针对半拉格朗日稳定流体模拟方法的两点不足，本文提出一种 基于谱方法的实时流体模拟方法，它是有别于计算机图形学中现行 流体模拟方法的另一种大胆尝试。该方法将谱方法引入到实时流体 模拟中，把n s e s 做为一个整体将流体从通常的实空间转换到谱空间 进行计算，避免了拆分n s e s 所产生的误差，并在谱空间下消去压力 项，避免了多次迭代运算产生的数值耗散，提高了求解精度。 硕士学位论文 基于谱方法的实时流体模拟 3 1 谱方法的基本思想 谱方法 2 6 ， 2 7 ( s p e c t r a lm e t h o d s ) 是一种求解偏微分方程 的数值方法。早在1 8 2 0 年，n a v i e r 运用双重三角级数求解弹性薄板 问。然而，由于该方法计算量大，使其长期没有得到广泛的应用。 直到1 9 6 5 年，c o o l e y 和t u k e y 提出快速傅立叶变换( f a s tf o u r i e r t r a n s f o r m ，简称f f t ) 2 8 才使得谱方法的应用出现转机。现在， 谱方法的数值分析理论不断完善，并广泛地运用于物理、大气、海 洋等领域的数值计算。谱方法已和有限差分法、有限元法一起成为 偏微分方程数值解的三种基本方法。 众所周知，各种数值方法都要考虑一个近似解的有限维子空间 b 和一个到巩的投影算子昂。在谱方法中，b 是选一个正交规范 完备函数系 形( x ) ) ，( 万= 0 , 1 2 ) 的前n 个基函数所张的空间；而投 影算子目选取为把函数“从某个泛函空间投影到风中的算子，即 日“= 艺口州( f ) ( 工) ，这里口州o ) 是以适当的方式确定的系数，使 艘肛一b “i i ：o ，则”笔口叫( f ) 嵋( 石) 。选定近似解的形式后，按照某 种意义下：在谱方法中试探函数取为无穷可微的整体函数，根据检 验函数的不同选取，谱方法可分为g a l e r k i n 谱方法、t a u 方法和配 置法。g a l e r k i n 谱方法有时就称为谱方法，配置法有时也称为拟谱 方法。在g a l e r k i n 谱方法中检验函数与试探函数属于同一个空间， 基于谱方法的实时流体模拟 并要求满足边界条件：t a u 方法类似于g a l e r k i n 谱方法，但不要求 检验函数满足边界条件，而是利用边界条件再补充一些方程，最后 得到一个封闭的方程组；配置法则是取检验函数为以那些配置点为 中心的d i r a c - 艿函数，使得微分方程在这些配点上精确成立。 谱方法有其自身的优点“无穷阶收敛性 ，即如果原问题的解充 分光滑，那么谱方法的收敛阶将是无穷阶的；此外，快速算法的使 用可以大大减少计算量。当然谱方法也有其不足，其一是要求原问 题解的正则性较好；其二是要求解区域比较规则，一般是典型乘积 型区域。 本文所采用的是谱配置法，因此，这里主要介绍谱配置的基本 思想。考虑微分方程的初边值问题 三( “) = ，工q ，0 t r ( 3 1 ) “i 铀= “量，“i r = o = “， 设近似解的形式为 甜。( 五，) - t 4 0 ( 五f ) + 口o 蟛( j ) ( 3 2 ) ，互l 其中九是已知的解析函数，通常称为试探函数( t r i a lf u n c t i o n ) ， 这里选取试探函数为妒，( x ) = e 弘，u o ( x ，f ) 被选择来满足初始和边界条 件。这里口) 为待求系数 令残量 r ( a i ，0 2 ，jo $ 口，z ，y ) = l ( u 。) - f ( 3 3 ) 为了获得口，的方程，令带权残量的内积为0 ， ( r ，吼( x ) ) = o ， k = 1 , 2 ，n( 3 - 4 ) 硕士学位论文 其中( 力称为权函数或检验函数( t e s tf u n c t i o n ) ，这里取权函数 为d i r a c - 8 函数 吸(