（计算机应用技术专业论文）基于格子boltzmann的流体景观建模与快速绘制.pdf

2 0 1 1m a s t e rd e g r e et h e s i su n i v e r s i t yc o d e ： 1 0 2 6 9 s t u d e n ti d ：5 1 0 8 1 5 0 0 0 2 7 eas tchinanorma _ univ er sit y q u i c km o d e n ga n dr e n d e r i n gf i u i d s c e n eu s i n gl ai c eb o l t z m a n n m e t h o d ( s u p p o r t e db yn a t i o n a in a t u r a is c i e n c ef o u n d a t i o no fc h i n a ，n o 6 0 6 0 3 0 7 6 ，n o 6 1 0 7 0 1 2 8 ) d e p a r t m e n t ： 三q ! 主巡坌 旦e 0 9 i n 金全 i 凸g ! d 主i 主些主金 m a j o r ：q m 巳竖主里 q 巳! i 坌主i q n ! 垒b n q ! q g y r e s e a r c ha r e a ： q 幽乜垡主皇【鱼 墨乜b i s u p e r v i s o r ： b 坌d g 垒q 巡墨d g 苎q i 坌! 金【q ! 金墨q a p p l l c a n t ： 圣b 型q 巳金d g 圣b 坌d g m a y2 0 1 1 2 华东师范大学学位论文原创性声明 郑重声明：本人呈交的学位论文基于格子b o l t z m a n n 的流体景观建模与伙速绘制， 是在华东师范大学攻读硬再博士( 请勾选) 学f 诊期间，在导帅的指导下进行的研究工作 及取得的研究成果。除文中已经注明引用的内容外，本论文不包含其他个人已经发表或 撰写过的研究成果。对本文的研究做出重要贡献的个人和集体，均已在文中作了明确说 明并表示谢意。 作者签名：邀盘鲴刍 日期：2 0 11 年5 月2 8 同 华东师范大学学位论文著作权使用声明 基于格子b o l t z m a n n 的流体景观建模与快速绘制系本人在华东师范大学攻读学 位期间在导师指导下完成的碌生博士( 请勾选) 学位论文，本论文的研究成果归华东师 范大学所有。本人同意华东师范大学根据相关规定保留和使用此学位论文，并向主管部 门和相关机构如国家图书馆、中信所和“知网 送交学位论文的印刷版和电子版；允许 学位论文进入华东师范大学图书馆及数据库被查阅、借阅；同意学校将学位论文加入全 国博士、硕士学位论文共建单位数据库进行检索，将学位论文的标题和摘要汇编出版， 采用影印、缩印或者其它方式合理复制学位论文。 本学位论文属于( 请勾选) ( ) 1 经华东师范大学相关部门审查核定的“内部”或“涉密”学位论文宰， 于年月同解密，解密后适用上述授权。 ( ) 2 不保密，适用上述授权。 导师签名本人签名缢剑刍 矽c 降f 月同 事“涉密”学位论文应是已经华尔师范人学学位评定委员会办公室或保密委员会审定过的学位 论文( 需附获批的华尔师范人学研究生中请学位论文“涉密”审批表方为有效) ，未经上 述部门审定的学位论文均为公开学位论文。此卢明栏不填写的，默认为公开学位论文，均适用 上述授权) 。 堂皇鹏硕士学位论文答辩委员会成员名单 姓名职称单位备注 谢步瀛教授，博导同济大学主席 何高奇副教授华东理工大学 全红艳副教授华东师范大学 摘 计算机自然景观模拟中诸多应用如电影特效，视频游戏中都需要对流体现象 忠实再现。然而简单的计算规则很难生成高度真实感的流体，使用物理模型对流 体现象合适的描述和建模成为计算机图形领域中的关键技术。另一方面物理计算 的复杂性往往占用了大量的计算机资源，限制当前应用，因此如何在视觉误差范 围内对其计算处理进行加速就显得尤为重要。 本文以基于物理的流体动画作为研究目标，并试图在物理模拟的真实感和效 率之间寻求一个良好的平衡点。格子波尔兹曼方法( l a t t i c eb o i t z m a n nm e t h o d ， l b m ) 作为流体力学计算的新方法，其计算快速，能适用于复杂边界障碍的特性 使它成为国内外基于物理的流体模拟研究中的热点。 本文将基于l b m 模型进行流体现象仿真，并在网格构建、动态环境交互、 多重介质混合模拟中寻找新突破，同时利用g p u 技术来进行l b m 计算的加速。 论文首先说明l b m 方法的历史与基本原理，并论述它并行性佳、边界处理灵活 等特点。然后讨论如何在g p u 编程模型中，实现格子b o i t z m a n n 方程计算，进 行实时的流体绕流模拟，并进一步采用g p u 更新的大规模粒子系统与实时光线 折射、r a y c a s t i n g 绘制技术，将物理计算与渲染完全在g p u 上加速实现；其次 提出了一种新的自由表面流体模型的自适应网格；另外根据刚体冲量运动学在 l b m 模型中模拟固体液体动态互相作用；最后使用优化的l b m 计算模型完成了 不同自由表面流体及流体型自然景观的真实感模拟。 关键词：l a t t i c eb o i t z m a n n 方法；流体模拟；图形处理器； s a b s t r a c t m a n ya p p l i c a t i o n s ，s u c h 舔 sm o v j e sa n dc o m p u t e rg a m e su s ec o m p u t e rg r a p h j c t e c h n o l o g yt 0c r e a t en a t u r es c e n e s r e a | s t i cf i u i ds i m u i a t i o ni sa l w a y sr e q u i r e d d u et ot h e c o m p i e x i t yo ft h ef i u i db e h a v i o d e s c r b i n ga n dm o d e l i n gf i u i db a s e0 np h y s i c si s ac r u c i a i t e c h n i q u ei nt h ef i e l do fc o m p u t e rg 陀p h i c w h i i et h es i m u | a t i o no fc o m p i e xf i u i ds y s t e m sb a s e o np h y s i c sr e q u i r e sm a s s i v e m p u t i n gr e s o u r c e 1 tl i m i t st h ep r o d u c t i o n0 fm a n ya p p i i c a t j o n s s o h o wt og e tap r e t t ys i m u i a t i o nr e s u i tw i t har e a s o n a b i ec o m p u t i n gs p e e di sv e r yi m p o r t a n t t h i sp a p e rw i t r yt 0b a i a n c et h ee f f e c ta n dt h ee 仟i c i e n c yo ff i u i ds i n l u i a t i o nb a s eo nt h e p h y s i c spr n c i p i e f o rt h en e wh o t s p o t ，l a t t i c eb o i t z m a n nm e t h o d ( l b m ) h a st h ea d v a n t a g ei n c o m p u t i n gs p e e d ，b o u n d a r yh a n d l i n g ，m a n yr e i a t i v ew o r i sa 伦b e i n ga d v a n c e db yw o r i d w i d e r e s e a r c h e r s t h ew o r k s0 ft h i sp a p e ra r ef o c u so ng r i dc o n s t r u c t i o n ，i n t e r a c t i o nw j t he n v i r o n m e n t a n dt o a c c e l e r a t et h el b mc a i c u i a t i o nb yt a k i n ga d v a n t a g eo ft h eg p u n l p u t i n gc a p t i v i t y h lm i sp a p t h eb a c k g r o u n da n df o u n d a t i o no fl b mw i l lb ee n 啪e r a t e d t h e n 、粑谢l lt a k en o t i c eo ft l l e c h a r a c t e ro fl b ms u c ha sh i 曲l ys i m u l t a n e o u sc a l c u l a t i o na n dn e x i b i l i t ) ro nb o u n d a r yp r o c e s s i n g 7 n l ef o l l o w i n gc h a p t e r 诵l l 昏v et h ei m p l 锄e n to fl b mo ng p up r o g r a 玑t h er e n d “n g 加d p h y s i c sc a l c u l a t i o no fn u i d 谢l lb o t hr u no ng p u t l l er e n e c t i o n r e 触c t i o n ，p a n i c l es ) r s t 锄a n d r a y - c a s tr e n d 耐n g 、) l ，i l lb eu s e do nn u i dv i s u a l i z a t i o n f u n h e m o r e ，a na d 印t i v e 鲥dm e t h o do f 丘s u m c ei si m p r o v e d f i n a 慨i m p u l s e 埘n c i p l ei n g i dm o t i 伽i su t i l i z c dt o0 p t i m i z et h e i n t e r a c t i o nb e m e 饥n u i da n dr i g i do b j e c t t h er e s tp a r t so fp a p e rs h o wh o wt o a p p l yt h e s e n l e t h o d so nr e a i i s t i cf i i j i ds i m u i a t i o n k e yw o r d ：l a t t i c eb o i t z m a n nm e t h o d ，f i u i ds i m u i a t i o n ，g p u 6 第1 章绪论 1 1 引言 流体模拟在诸多领域中是一项关键技术。石油、化工、气象预测、机械设 计等领域中需使用流体模拟进行数值分析。计算机自然景观模拟中流体模拟也有 较多应用，如电影特效、视频游戏中对流体现象的再现。电影如泰坦尼克号、未 来水世界等等都需要大量精致的水的特效，视频游戏中水花、血滴、烟雾等则需 要快速而令人信服的流体模拟。 利用物理原理来进行流体建模与模拟，可以使自然景观的仿真更加逼真。 然而流体运动蕴含着复杂规律使之模拟成为计算机图形学中难题之一。研究者们 不断地提出新的模拟方法并进行改进，推动着流体模拟课题的发展。 基于物理的流体模拟来自计算流体力学( c f d ) ，通过求解经典的n a v i e r s t o k e s 方程的方法，可以精确地计算出流体的运动状念，国内外也已经有了大量的研究 工作。但由于n s 方程式基于连续介质假设，对自由表面流体的离散和扩展性存 在困难、较难求解。近年来流行于计算流体领域的l a t t i c eb o i t a m a n nm o d e l ( l b m ) 方法则提供了另一种流体动力学表达方法，其主要思想是把大量时空离散的微观 粒子的迁移与碰撞抽象为分布函数的演化，利用这些微观粒子运动的统计结果来 反映流体现象的宏观特性，便于离散求解，具有较好地并行性。 1 1 1 研究背景 为了表现符合物理性质的流体，流体模拟的计算需考虑大量分子受力与相 互作用，即使是模拟一微小空间内的流体，都要求大量的计算，至使计算机在短 时间内无法完成。同时对于高动态流体越要求真实表现，就需要求解越复杂的流 体方程。另一方面，流体是充满动态变化，形状并无规则可循，高效地绘制也是 影响流体模拟应用性能的直观要素。因此，在尽量保证模拟效果的前提下，使用 一些近似解，在流体模拟效果与流体模拟效率之间取得平衡具有重要意义，这也 是国内外众多学者的努力方向。 流体景观建模与绘制的应用可分为离线渲染和实时绘制。离线渲染要求生 成高质量并且高度逼真的图像，所以常常会十分苛求模拟的精确性而使用大规模 的流场描述数据构建复杂场景，结果消耗大量的时问在每一帧的流体计算中。所 1 以我们需要对计算过程进行网格构建优化、物理计算优化或使用g p u 加速技术 以提高模拟效率，减少离线渲染时问。而同时对于虚拟现实与视频游戏的实时绘 制，常折中使用经验模型与随机生成算法，这些加速技术可以帮助计算机在限制 的时问内完成物理计算，为用户提供更真实的视觉效果。 l b m 方法能够方便地处理复杂边界与障碍物的流体【m e i9 9 】，也可以用来模 拟多相流体【z h u 0 6 1 【z h u 0 7 】。鉴于l b m 方法的并行特性，g p u 集群技术也用来 进行流体模拟时的加速【q u 0 4 】。但是以上l b m 方法多用来模拟空间充满的流体， 对于气液分离流体的表面特征表现较少。 具有气液分离的自由表面流体是一种最常见的流体，如溪流、淌流、水滴 等。这类自由表面的流体由于具有更多表面细节和物理机理，比充满空间的单一 流体更加难以逼真模拟。所以其探索与实现需要更多的物理与计算机知识，以及 科研时间投入。 1 2 国内外研究现状 在计算流体力学方面，一直以n a v i e r - s t o k e s 方程作为基本的物理模型，求 解n a v i e r s t o k e s 方程是传统的流体力学的计算方法，但它是存在较难求出解析解， 其求解在连续介质上，计算量大的缺点。n a v i e 卜s t o k e s 方程派生的求解主要分为 两种：基于网格的流体模拟与基于粒子的流体模拟。基于网格的方法称为欧拉方 法，而基于粒子的方法则称为拉格朗同方法。 基于网格的方法中常用的有基于m a c ( m a r k e r a n d c e l i ) 【n i c k9 6 】网格的方法。 1 9 9 9 年j o ss t a m 以m a c 网格技术为基础，提出了网格方法与粒子方法结合的半 拉格朗r ( s e m i l a g r a n g i a n ) 流体模拟方法【s t a 9 9 】。该方案是确保不可压缩性一项 具体化的步骤，并包含了处理粘性的方法，其流体模拟稳定性出色，模拟强流效 果突出。 基于粒子的方法中最引人瞩目的是s p h l 光滑粒子动力学s m 0 0 t hp a r t i c l e h y d r o d y n a m i c ) 方法，其最早于1 9 7 7 年提出【g i n g o i d7 7 】，是流体自由表面的一种 十分有效的模拟方案。流体飞溅、与场景问的交互都适合使用s p h 实现。由于粒 子本身携带质量，压强则可以通过和邻近粒子加权得出，所以整个计算式完成过 程十分简洁【m u e r 0 3 】【m u e r0 5 】。 网格的方法的缺陷是模拟区域由网格限定，流体无法在场景中任意处出现， 同时细节表现需要大量网格或是特殊的计算，表面需要跟踪构建。 粒子的方法则能让流体在场景中任意方位出现，且粒子自身就能表示流体， 相比网格的方法更能表现流体的剧烈形变。但是其遍历粒子可能会花费一定计算 量，为优化粒子查询需要建立粒子查找表以查询临近粒子。 气体在某些条件下( 如温度变化等) 会呈现可压缩性，而液体一般视为不 可压缩的流体。基于粒子的方法无法实现流体的绝对不可压性。相对的，基于网 格的方法可以表现不可压缩的流体。综上可知两种方法各有优劣，在实际运用中 则可以按需求做选择。 幽1 1s p h 模拟效果f l e n a e r t s 0 9 】 ： f i 9 1 1 s p hs i i n u l a t i o n si n 【l e n a e r t s 0 9 】 l b m 方法基于网格，同时借鉴了了统计学的思想。它宏观统计微观粒子的 运动，离散简化了流体模拟的求解计算，同时满足质量和动量的守恒，其计算都 产生在离散的网格节点上，故同时又具有计算并行的特性。 l b m 方法的一个优点是它能处理复杂边界与障碍物的流体，基于m e i 【m e i9 9 】 提出了曲面边界的处理方法，在2 0 0 2 年由w e i 等【w e i0 2 1 首先在图形处理器上 模拟了火焰，随后还模拟了空气中的烟气【w e i0 4 】。美国s t o n yb r o o k 大学【l 0 4 】 在图形硬件上实现的l b m 算法能在g p u 上处理复杂边晃的流体，其模拟的障碍 物还可以变形，比如在水中游动的水母。y ez h a o 【z h a o0 6 】使用格子b o i t z m a n n 模型模拟了数种自然现象如蜃景、融化、洪水蔓延等，在g p u 上实现了自适应 多分辨率的l b m 网格模型，还使用了g p u 集群模拟了烟雾在整个城市的蔓延 【q i u 0 4 1 。国内的l b m 研究也如火如荼，较多人关注方柱绕流模拟与流体的稳定 性【王0 3 】，使用不同的雷诺数模拟不同黏度的流体，另一些研究者在c u d a 中实 3 现了l b m 的计算【吴1 0 】。l b m 方法也适合模拟飘雪【w a n 9 0 6 1 等轻质粒子。朱红 斌等人【朱0 6 】在g p u 上引入了一种二元格子b o i t z m a n n 模型，实现了两种液体 组成的混合流的模拟，表现不同粘性的液体间的的扩散( 互溶) 和分离现象。 对于自由表面流体的模拟，k a s s 等【k a s s 9 0 】最早提出浅水方程的思想，其模 拟并非流体的截面，而将流体的高度假定为各个格点的宏观密度，再基于高度生 成水面。z h o u 等【z h o u 0 2 】在l b m 的d 2 q 9 网格中求解浅水方程，完成了渍堤模型 的数值模拟。t u b b s 等【t u b b s 0 9 】进一步提出了多层的浅水方程，将多个2 d 网格 层叠起来得到最后的格点高度，并对其计算性能和效果进行了验证。l i u 等 【l u 0 9 】2 0 0 9 年还提出了在多重网格下来实现浅水方程，即在需要更多细节的地方 人为划定密集的网格，从而在不大的运算量下能得到更好的效果。 流体表面的抽取是自由表面流体模拟的关键。e n r i g h t 等【e n r l g h t0 2 】进一步 提出p a r t i c i el e v e is e t 方法，为水的自由表面增加了更多的细节，大大改善了自 由表面模拟的真实感，然而其表面重构花费了较大的计算资源。t h u r e v 等 【t h u r e y 0 6 】使用控制粒子的思想完成对水体的控制的动画同时保留局部细节。基 于l b m 的自由表面方法提供了并行计算的动态边界处理，文献【t h u r e y 0 6 】中提出 了处理移动障碍和变化气液边界与l b m 的d 3 q 1 9 模型结合的数值模拟方法。海 浪模拟也是十分热门的一项图形学科目，l b m 模型加上g p u 技术模拟的海浪 【f n c h0 4 】可以同时满足大规模海浪模拟与效果需求。文献【t h u r e y 0 4 】中基于l b m 自由表面方法结合l e v e ls e t 方法来计算水体表面曲面和模拟较薄的水体，由于涉 及预测流体速度和保证质量守恒，难以实时绘制。文献【n i p u nk w a t r a1 0 】中描述 一种新颖的刚体与流体互动，其主要贡献在与完成骨骼( 关节) 动画与流体求解 的耦合，可以模拟人偶在水中游泳。 1 3 本文工作 本论文将分析l a t t i c eb o i t z m a n n 方法，利用l b m 方法建立流体模拟的物理模 型，优化计算模型并利用流体场景快速绘制的方法，模拟了几种典型的流体自然 景观；使用g p u 加速技术完成l b m 模型计算同时建立一风场的实时模拟框架； 另外针对自由表面流体的特点，基于l b m 方法结合浅水方程与网格标记法，来 进行自由表面流体的建模；并对复杂的流体表面进行抽取，最后实现了不同流体 的真实感仿真。 五 本文的主要贡献和创新点在于如下工作： ( 1 ) 利用g p u 来加速l b m 方法的代数运算，结合阴影图、g p u 粒子系统、r a y c a s t 绘制方法，使流体模拟和绘制都于g p u 中进行。 ( 2 ) 提出一种基于视点的非均匀网格优化方法，对自由表面流体的l b m 模型 进行了优化。 ( 3 ) 利用基于冲量的刚体运动学结合l 8 m 方法模拟固体物件与液体的交互作 用。 ( 4 ) 基于以上方法逼真模拟了不同情况下的自然景观场景，如风雪、烟、水滴、 溃堤、以及泥石流等现象。 本论文的整体结构安排如下： 第一部分：阐述文章研究背景和目的； 第二部分：简述l a t t i c eb o i t z m a n n 方法基本理论及应用； 第三部分：介绍基于l b m 方程驱动的大规模粒子系统场景模拟实现； 第四部分：介绍基于l b m 的自由表面流体模拟与刚体交互模拟； 第五部分：一些加速计算技术与大规模流体场景绘制实现； 第六部分：总结和展望。 第2 章l a t t i c eb o l t z m a n n 方法基本原理 近些年，i a t t i c eb o l t z m a n n 方法成为了流体物理建模与数值模拟的重要方法。 l b m 方法的一般思想是在网格节点上分布微观的粒子，让粒子在节点上发生碰 撞( c o l i s i o n ) ，并迁移或平流( p r o p a g a t i o no rs t r e a m i n g ) 运动到相邻节点上，在各个 时刻，微观粒子会有一个确定的分布。基于此分布，采用统计的方法，可以得到 物质的宏观物理量的分布，如密度场和速度场，从而实现了流体的运动模拟。 图2 1l b m 模型的发展 f i 9 2 1h i s t o r yo fl b mm e o 秽 2 1l a t t i c eb o l t z m a n n 方程 b o i t z m a n n 方程描述了流体的宏观特性，流体被描述为连续介质。在宏观体 系中，每个分子的微观运动都遵守力学规律，而b o i t z m a n n 方程的基本思想则是 求出每个分子在某个状态下的概率，通过统计的方法得出系统的宏观参数。 b o i t z m a n n 方程是统计力学中用来描述非平衡态分布函数演化规律的方程。连续 的b o i t z m a n n 方程( b o i t z m a n ne q u a t i o n ) 是一个比n a v i e r - s t o k e s 方程组更加复杂的 微分方程，它的精确求解是不现实的。 格子b o l t z m a n n 来源自格子气自动机，格子气自动机是元胞自动机( c e i i u l a r a u t o m a t a ) 在流体力学中的应用，元胞自动机是一个时问和空问的离散数据模型， 按简单的规则演化，可以用在细胞生长，城市交通流，舆论传播等复杂问题。1 9 8 8 年，最早的格子b o l t z m a n n 论文出现【m c n a m a r a 8 8 】，m c n a m a r a 和z a n e t t i 提出把 格子气自动的整数运算变成实数运算，粒子分布函数用实数表示，最早的格子 b o i t z m a n n 模型诞生。 格子b o i t z m a n n 方程通过c h a p m a n e n s k o g 展开，可以恢复为n a v i e r - s t o k e s 方 6 程。 2 2l b m 模型 格子b o l t z m a n n 方程是n a v i e 卜s t o k e s 方程组的近似求解方程，它的求解没 有n a v i e r s t o k e s 方程组复杂，但模拟效果与求解n a v i e r s t o k e s 方程相当。 2 2 1n a v i e r - s t o k e s 方程组 为了后面的话题展丌，先了解一下n a v i e r s t o k e s 方程组如何描述流体的物理 特性的，n a v i e r s t o k e s 方程组： v “= o 2 1 ) 鱼栅砚= j 即+ 舻“+ ， i j i p 2 2 ) 方程中如力为矢量速度场。舷力为标量压力场。p 是流体密度( 常数) ，y 是动力黏度。 罢描述液体加速度，一! v p 是压力项，闪：“代表扩散项，f ：饭，) 代表作i 描述液体加速度，一v p 是压力项，闪2 “代表扩散项，声= 执，) 代表作 成p 用于流体的外力。方程( 2 2 ) 满足时，同时方程( 2 1 ) 满足，即速度的散度v “ 为0 。 n a v i e 卜s t o k e s 方程组的求解涉及h e i m h o i t z h o d g e 分解推导；计算平流使用 隐式积分逆向求解e u i e r 方法；扩散的计算要求隐式求解p o i s s o n 方程( j a c o b i 迭 代) 。 2 2 2 格子b oi t z m a n n 方程 这一小节先介绍格子b o l t z m a n n 方程的理论。连续的b o l t z m a n n 方程形式： 和胪嘲 1 2 3 ) 厂代表分布函数( d i s t r b u t i o nf u n c t i o n s ( d f ) ) ，善是粒子速度，q ( 厂) 是碰撞的结 果。q 的选择至关重要，它保证了流体局部的质量与动量守恒。 根据【b h a t 响g a r 5 4 】的工作，在每个结点的包分布厂，存在一个局部平衡包分布尸， 7 尸线性依赖于结点的密度p 和速度“。由此，碰撞因子可以如下定义： ：三一尸) 2 4 ) f f 为流体的松弛时i 脚，决定了流体的粘性y ： y = 三p 一乏， 。2 5 ， 可以通过数值求解b o i t z m a n n 方程来模拟流体的宏观运动。然而求解 b o i t z m a n n 方程存在一难点：碰撞项是分布函数的非线性项，它与具体分子作用 力有关。因此b h a t n a g e r 、g r o s s 和k r o o k 提出b g k 模型，采用b g k 近似简化碰 撞项。 格子b o l t z m a n n 方程是b o i t z m a n n b g k 方程的一种特殊离散形式。包括速度 离散、时i 日j 离散和空问离散。 微观粒子无时无刻不在做着无规则的运动，而且速度的方向是任意连续的。 而在宏观运动方面我们可以把粒子的速度孝离散到有，q ，吃) 限维的速度， 其中6 表示速度的数目。同时，分布函数厂也相应的被离散为矾，彳，五) ，其 中厶= ，f ) ，g f _ o ，1 ，6 ，意味着某个格子节点薯在某一时刻f 在g f 方向的 分布函数。 单松弛模型中，c o l l s i o n ( 碰撞) 和s t r e a m ( 迁移) 是基本步骤。不含外力 项的格子b o i t z m a n n b g k 方程定义为： ，( x + ，f + 万f ) 一i ( x ，f ) = 一 ( ，( x ，f ) 一孑( p ，“) ) ( 2 6 ) 其中，( 五f ) 和劈( p ，“) 是粒子速度分布函数以及平衡分布函数，f 为无量纲松 弛时间。 此方程可分为两部分( 2 4 ) 与( 2 5 ) ，它将碰撞项完全线性化，以单松弛时 问逼近来表现碰撞效果，使得系统趋于平衡。 碰撞：元，( x ，f ) = ，( x ，f ) 一( 。( x ，f ) 一巧9 ( 尸，“) ) ( 2 7 ) 迁移：厶( x + 疣，f + 毋) = 乃( x ，f ) ( 2 8 ) 其中无，g f f ) 为碰撞后格点上的分布函数。 口= 4 ，+ 白“+ 去( 白甸2 一去( “功 c2 舢 8 其中后为方向上的平衡分布函数，铭，为对应粒子速度的权重，岛，为离散速度， p 为宏观密度， 为宏观速度，c 为格子声速，c = 夙所。缸、所分别为网格 步长和时间步长，通常x 和y 方向步长相同( 二维网格中) ，即万z = 砂。 格子b o i t z m a n n 方法对粒子分布函数进行统计计算，得到压力和流速等宏观 变量。 p ( x ，f ) = 如，f ) ， ( 2 1 0 ) g ，= 0 “( 圳) = 去心，f ) ( 2 1 1 ) y i = 0 处理单个 液体格了 “ 、l ；i ； 兰口j 髑_ * ；i ； 一 邻肼格子 d f 辽移 杉f l i ；迂移 进入的d f ( ：孽：渗i 蕊 ；旷网 计算密度速度、 d f 的碰撞 存储完成计并的 d f ，处坪f 格了 缓 图2 2 碰撞和迁移步骤流程图 f i g2 2s t 印so fs t r e a ma n d c o l l i s i o n 此外，时间和空间的离散并不是独立的，而是通过粒子的离散速度联系起来，有 舐= 乞，国。这一特性告诉我们如何把粒子在物理空间的运动分为碰撞和迁移， 即粒子在两个时间步之间由一个节点晗好运动到对应的相邻节点，并在该节点上 与其他粒子碰撞，这使格子b o l t z m a n n 方程具备了很好的并行特性和较强复杂边 界的处理能力。 流体模拟的格子b o i t z m a n n 方法有许多格子类型，即离散速度模型，q a n 等 人【q - a n 9 2 】提出的d d q m 系列模型是格子b o i t z m a n n 方法的基本模型。 d d 代表 了d 个维度，q m 代表了m 个离散速度。d 2 q 9 ，d 3 q 1 5 ，d 3 q 1 9 ，和d 3 c 眨7 模型 9 这些模型都有统一的分布函数公式( 2 9 ) 。 d 1 q 3 d 2 q 9 图2 3d 1 q 3 与d 2 q 9 模型构造 f i g2 3d 1 q 3 d 2 q 9m o d e i sg e o m e t r y 如图2 3 ，d 1 q 3 是一种一维的l b m 模型，d 2 q 9 则是最常用的二维的l b m 模型。 i 悻0 卢l 。sf 生， 一 。一_ i - 2 卜_ 1 ”一 ，1 r t 一一， 一， i 净3 t 二： ，i 图2 4d 3 q 1 5 与d 3 q 1 9 模型构造 f i g2 4d 3 q 1 9 & d 3 q 1 9m o d e ig e o m e t r y 图2 4 展示了l b m 的d 3 q 1 5 和d 3 q 1 9 模型。如图( a ) 离散速度分为四种：离 1 0 一，l，；_1， 一，置。一 散速度的位移大小i 1 分别为0 ，1 ，2 ，3 ，方向分别为指向自己( 静止) ， 立方体面中心上的6 个邻居，立方体棱中心上的1 2 个邻居和位于立方体顶点上 的8 个邻居。d 3 q 1 5 取用了大小为o ，l ，3 的离散速度组合，而d 3 q 1 9 取用了 大小为o ，1 ，2 的离散速度组合。 本论文实验采用的二维l b m 模型为d 2 q 9 。三维的则采用d 3 q 1 5 和d 3 q 1 9 ， d 3 q 1 9 在模拟的稳定性和模拟效果上优于d 3 q 1 5 ，而d 3 q 1 5 则可以在更少的存 储空i 日j 和计算量中实现。 2 3 边界处理 边界处理是在格子b o i t z m a n n 方法实施中非常关键的一项内容。根据应用 的不同，l b m 中有各种不同的边界条件处理方法。在边界附近迁移运算的求解 必须作特殊处理。 2 3 1 规整边界处理格式 ( a ) 周期性边界 ( b ) 标准反射( c ) 镜面反射 图2 s 常用的边界条件 f i 9 2 5t h ec o m m o n i yu s e db o u n d a r yc o n d i t i o n s 图2 5 是常用的三种边界条件处理：( a ) 为周期性边界( b ) 标准反射( c ) 镜面反射。其中周期性边界( p e r i o d i cb o u n d a r y ) 可以用于表现流场在空间呈现 周期性变化或在某个方向无穷大。当流体粒子从一侧边界沿着一个离散速度离开 流场，下一个时步就会从流场的另一侧边界沿同一速度向量进入流场。使其在 分布函数的计算的同时把两端边界连接，如图2 5 ( a ) 左右两边界上的节点成为 了邻居节点。 标准反射( b o u n c e - b a c k ) ，也称无滑移边界条件。对于静止固体边界上的粒 子作弹回处理。边界表面没有速度但是具有摩擦力，如图2 5 ( b ) 的情形，向着 边界表面移动的粒子会按原路弹回。 镜面反射( s p e c u i a rr e f i e c t i o n ) ，即滑移边界条件。假设流体与壁面没有摩擦 作用，如图2 ( c ) ，向边界运动的粒子，就会弹向原速度关于法线的反射方向。 在内部流场节点的分布函数计算完成后，根据反弹格式进行边界周围节点上 的分布函数更改。 图2 5 示意的2 维反弹情形可以类推扩展为3 维。图2 6 给出了另种标准 反射的分布函数计算方式，深色格子表示边界，下半部分格子是流体，边界上方 增加的一层虚拟节点( j = o ) 每一时间步分布函数都会被更换为( j = 2 ) 的弹回方 向的分布函数。 ；：i 业；：五 j o ；：譬；，名； 缪黟嬲黝缪嬲缈 磁貂缸端赫毋& 。级。钰。g 簇红磊# ；镕彩 84 7t 9 ，847 ：；：歹麓馨： 图2 6 另一标准反射的分布函数计算方式 f i 9 2 6d i s t r i b u t i o nf u n c t i o nc o m p u t a t i o no fs t a n d a r db o u n c e - i 培c kc o n d i t i o n 2 3 2 曲面边界处理 l b m 方法的一个优点是它能处理包含复杂边界与障碍物的流体。2 3 1 节介 绍的处理方式要求障碍物占满格点，然而在模拟自然场景中，往往存在更复杂的 物件对象，在论文【m e i9 9 】罩给出了曲面边界障碍的处理方法，障碍物边界允许 穿过网格节点，使网格节点部分被障碍物占据。 对任意的曲面边界，依据障碍物边界与网格离散速度向量相交的几何信息， 在交点上安排中间节点，在格点分布函数的计算后再根据中间节点的位置重新计 算边界两边格点的分布函数，这样完成的边界条件处理更为精确。 图2 7 在d 2 q 9 上布局的曲面边界 f i g2 7d 2 q 9c u r v e db o u n d a r yt r e a t m e n t 在图2 7 中，可以看到曲面边界的几何配置，曲线为边界，黑色圆点( ) 标记了粒子与边界的碰撞点位置，空心圆圈( 。) 和狄色圆点( 。) 分别标记了流 体和固体节点。曲面障碍物的边界( 在模拟空间为二维情况下是曲线，三维情况 下是网格面) 穿过格点间的链接( 这罩拿d 2 q 9 模型解释下链接：d 2 q 9 有8 个 离散速度指向周围8 个邻居格点，就是说每个d 2 q 9 的格点存在8 个连接邻居的 链接) 。 图2 7 中x w 是粒子与边界的碰撞点，也是障碍物边界与格点问的链接的交点， 交点两边的节点是和誓，是固体节点，号则是流体节点。当从出发的粒 子沿方向白运动，碰到边界x w 反弹，即沿返回( = 一乞，) 。这被称为 b o u n c e - b a c ko nt h ei i n k ( b b l ) 机铝0 【l a d d 9 4 】。 在【m e i9 9 】介绍的曲面边界处理中，用了虚构的分布函数模拟反弹的粒子。 是根据交点位置得到的流体部分的比例，大小为： ：譬型 1 2 1 2 ) l o 一i 碰撞步结束后，流体节点的分布函数已知，而来自固体节点的分布函数 乃( ，f ) 未知，为完成迁移步乃( _ = 吒+ 勺况，f + 挽) = 乃( ，f ) ，需要构造出此 分布函数，f j i i p p o v a 与h a n e l ( f h ) 提出以下线性差值方法【f i i i p p o v a 9 8 】： 无( ，f ) = ( 1 一z ) 五，( 。，f ) + z 臂( ，f ) 一2 ，p 砉( “。) ( 2 1 3 ) 其中兰“k ，f ) 是界面上的速度， z 是差值权重、控制五，( _ ，f ) 和 瑟( 讫，f ) 的线性差值。瑟( 毛，f ) 的计算方式如： 删w ) 魂嘶力 h 主撕+ 刍( ) 2 - 砉驴吩 ( 2 1 4 ) 吁三“( 号，f ) 是邻近界面的流体速度，的取值不唯一且与z 选取相关，以下是 f h 提出的选取方式，并证明了这种方式比单纯设置“6 ，= “，更能保持计算稳定： = 去( 一班p 是，a n 略掣，f o r 三； ( 2 1 s ) 旷即a n d z = 蹦，f o r 6 0 f p s d 3 q 1 56 4 6 4 + 6 45f p s 6 0 f p s 如表1 ，g p u 的l b m 计算实现运算速率比c p u 要高2 1 0 倍。所以当网格越 庞大，越能体现g p u 计算的性能优势。 3 2 风场中雪、烟现象的模拟 这罩我们设定模拟的流体对象为空气，空气的流动同时带动了一些轻质小物 体，如粉尘、烟雾、雪之类微粒。使用这章介绍的硬件加速计算技术，我们可以 2 6 o 实时模拟具有较高复杂度的空气流动 3 2 1 基于阴影图的障碍物标记 在l b m 流体求解之前，我们需要得到场景的几何信息，即在l b m 网格上标 注出障碍或地面。【e v e r i t t0 1 】中提出的基于g p u 的深度剥离体素化( v o e l i z a t i o n a sd e p t hp e e i n g ) 的技术可以处理任意复杂的边界。本文提出一种简化的障碍物生 成方法，使用阴影图( s h a d o wm a p ) 的思想，记录场景的深度信息到一张浮点 纹理中。 如图3 9 ，我们把摄像机安设在场景正上方并使用j 下交投影绘制场景。场景 中所有物体的y 值( 高度值) 会转换为相机空问的z 值( 深度值) 。 图3 9 深度图生成流程图。左边是在场景上方设置相机的示意；右边是物体被绘制剑相机空 间成为深度数据的示意。 f i g3 9c a m e r as e t t i n ga n dr e n d er n gd e p t hv a i u ep r o c e s s t h ei e f ts i d es h o w sc a m e r as e t0 nt h et o po fs c e n e t h er i g h ts i d es h o w sh o wa no b j e c tr e n d e r e di n t oc a m e r av i e wa n dd e p t hv a i u eg e n e r a t e d 按图3 9 所示设置相机，通过绘制入纹理可以生成图3 1 0 中的深度纹理。生 成的深度纹理将成为初始化l b m 网格节点的片段着色器的输入。将场景的几何 信息记录在二维纹理中也称为降维。经过这一步骤，场景的最大高度己记录，我 们可以比较l b m 网格的所在高度和深度纹理记录的场景高度，定义在场景表面 之下的网格点为固体节点，之后场景中的空气流动可以按l b m 方程求解。 该方法中深度纹理使用了一张纹理实现，也是使用g p u 通过绘制生成，所 以生成也很快速，只花费绘制一遍场景的时间，对于实时绘制的场景只使用少于 1 3 0 秒的时间即能绘制生成。 2 7 一一 图3 1 0 城市楼群和树林的深度纹理，离相机越近的部分越接近白色 f i g3 1 0d e p t ht e c t u