（计算机应用技术专业论文）泥石流现象的真实感建模与绘制.pdf

浙江大学硕士学位论文 摘要 摘要 泥石流是介于流水与滑坡之间的一种地质作用，它经常发生在峡谷地区和地 震火山多发区，在暴雨期具有群发性，是一种世界范围的能给人类带来巨大生命、 财产损失的自然灾害现象。由于泥石流具有爆发突然，而且来势凶猛的特点，对 其进行场景拍摄相当困难，因此对泥石流现象的真实感建模与绘制工作是非常必 要的，其研究成果可应用于防灾救灾、影视特技、数字娱乐以及科普教育等方面。 本文第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 ) r n a r n i c s 平滑粒子流体动力学) 方法，为接下来本文对泥石流现象的模拟打下了技术基础。 本文第3 章针对泥石流的运动特性，提出了一种基于粒子的新的固液两相流 模型对泥石流的运动进行建模。首先分析了泥石流的构成及运动特性，并采用改 进的s p h 方法结合刚体模型来模拟泥石流中液相浆体和固相颗粒的运动，然后基 于粒子模型实现了固液两相问的相互作用，并采用g p g p u 技术加速了计算的过 程。 本文第4 章讨论了泥石流现象的真实感绘制。在对泥石流实施绘制时，对其 液相浆体部分采用隐式表面定义，在保证绘制效果的同时，大大加快了其表面提 取的速度：对固相颗粒则进行了精细的造型并依据其运动模型进行绘制。为了增 强场景的真实感，本文考虑并模拟了泥石流与场景的交互作用，着重展现了泥石 流冲毁、淹没房屋的效果。另外，通过推广第3 章的固液两相流模型，本文还实 现了洪水、山体滑坡等单相流自然灾害的动态模拟。 在论文的最后，作者对全文的研究工作进行了总结，并提出了进一步的研究 方向和任务。 关键词：自然灾害模拟，泥石流，真实感建模与绘制，基于物理的流体模拟，平 滑粒子流体动力学，固液两相流，通用图形处理器，洪水，山体滑坡 浙江大学硕士学位论文 a b s t r a c t a b s t r a c t d e b r i sf l o wi sak i n do fg e o l o g i c a la c t i o nb e t w e e nw a t e rf l o wa n dl a n d s l i d e i t f r e q u e n t l yo c c u r si nm o u n t a i na r e a sa f t e rh e a v yr a i n ( o t h e rc a u s e sm a yb ee a r t h q u a k e ， t h ee r u p t i o no fv o l c a n o t h em e l t i n gi c ea n ds n o w , e t c ) i tc a nb r e a ko u tw t hl i t t e r w a r n i n ga n dr u s hd o w n 丽廿lf a s ts p e e d ，s ot h a tt h ec a p t u r i n go fd e b r i sf l o wi sv e r y d i f f i c u l t t h er e a l i s t i cs i m u l a t i o nc a nv i v i d l ys h o w 血ed y n a m i cp r o c e s so fd e b r i sf l o w a n di sh e l p f u lf o rt h el o c a lp e o p l et os i n - r i v e i tc a na l s ob ea p p l i e di n 也ea r e a ss u c ha s d i s a s t e rr e s c u ef o rg o v e r n m e n t s ，s p e c i a le f f e c t si nm o v i e s ，d i 西t a le n t e r t a i n m e n t , a n d p o p u l a rs c i e n c ee d u c a t i o n , e t c t i l i sp a p e ri so r g a n i z e da sf o l l o w s c h a p t e r1 i n t r o d u c e st h e 删ch a z a r d sc a u s e db yd e b r i sf l o w , r e v i e w sp r e v i o u s r e l a t e dw o r ki nt h ea r e ao fg e o l o g ya n dc o m p u t e rg r a p h i c s ，a n db r i e f l yd e s c r i b e st h e m a i nr e s e a r c hw o r ko ft h i sp a p e r c h a p t e r2i n t r o d u c e st h et h e o r yb a s i so fp h y s i c s - b a s e df l u i ds i m u l a t i o n ，e s p e c i a l l y t h es 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 ) m e t h o d 田：1 i sc h a p t e ri st h eb a s eo ft h e f o l l o w i n gs i m u l a t i o no fd e b r i sf l o w c h a p t e r3p r o p o s e san e w u n i f i e dp a r t i c l e - b a s e ds o l i d 1 i q u i dt w o - p h a s ef l o wm o d e l a c c o r d i n gt ot h em o t i o nc h a r a c t e r i s t i c so fd e b r i sf l o w w bf i r s tu s et h es p hm e t h o dt o s i m u l a t et h em o v e m e n to fl i q u i dp h a s e ，a d o p tt h er i g i db o d ym o d e lo u s te x i s t s t r a n s l a t i o na n dr o t a t i o n ) t or e s t r i c tt h em o v e m e n to fs o l i dp h a s e ，t h e na c h i e v et h e i n t e r a c t i o nb e t w e e nt h ea b o v et w op h a s e sb a s e do np a r t i c l e s i na d d i t i o n , w eu s et h e g p g p ut e c h n o l o g i e st os p e e du pt h ec o m p u t i n g c h a p t e r4d i s c u s s e s 也er e n d e r i n gr e s u l t so fd e b r i sf l o w f o rt h el i q u i dp h a s e ，w e i m p r o v e 血eo r i g i n a li m p l i c i ts u r f a c ef u n c t i o n w h i c hs p e e d su pt h er e c o n s t r u c t i o no f s u r f a c e a n df o r t h es o l i dp h a s e w em o d e lt h es t o n e sc a r e f u l l y t h e nw ec o n s i d e ra n d s i m u l a t et h ei n t e r a c t i o n sb e t w e e nd e b r i sf l o wa n do t h e rs c e n eo b j e c t s ，e s p e c i a l l yt h e e f f e c to fd e b r i sf l o wd e s t r o y i n ga n ds u b m e r g i n gb u i l d i n g s a d d i t i o n a l l y , w ei m p l e m e n t t h ea n i m a t i o na n dr e n d e r i n go ff l o o da n dl a n d s l i d eb a s e do nt h es a m em o d e l t h es u i n n l s r ya n dc o n c l u s i o no ft h i sp a p e ra r eg i v e ni nc h a p t e r5 。嬲w e l la st h e d i s c u s s i o no f t h ef u l l e rr e s e a r c hw o r k k e y w o r d s ：n a t u r a ld i s a s t e rs i m u l a t i o n , d e b r i sf l o w , r e a l i s t i cm o d e l i n ga n dr e n d e r i n g ， p h y s i c s - b a s e df l u i ds i m u l a t i o n , 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 ，s o l i d l i q u i d t w o - p h a s ef l o w , g p g p u ，f l o o d ，l a n d s l i d e 浙江大学硕士学位论文图目录 图目录 2 0 0 8 年5 月1 2 日四川汶川大地震震后场面2 2 0 0 6 年2 月1 7 日菲律宾爆发严重的泥石流4 2 0 0 7 年7 月1 9 日云南省腾冲县泥石流掩埋现场5 汶川大地震造成泥石流掩埋村庄5 p e a c h e y 于1 9 8 6 年模拟的波浪效果1 0 二维的欧拉网格，图中各实心点即离散的速度场1 3 2 0 0 1 年f o s t e r 等绘制的子弹滑过水面的效果1 3 1 9 9 5 年o b f i e n 等绘制的浪花飞溅的效果1 4 m i e k y 于2 0 0 6 年采用s p h 方法实现的浪头效果1 5 核半径为1 时的默认核函数分布曲线1 7 核半径为1 时的压力核函数分布曲线1 8 核半径为1 时的粘滞核函数分布曲线1 8 泥石流运动的四个阶段2 3 g p g p u 中的纹理数据一3 0 2 0 0 5 年z h u 等采用的隐式表面定义函数中存在的错误3 3 几种石头的材质纹理3 4 固相石块的绘制效果3 5 基于高度图的山体建模3 6 各种树、草的b i l l b o a r d 纹理3 6 泥石流中的房屋建模效果图3 7 泥石流运动的全过程3 8 不同固相体积比的泥石流3 9 泥石流中固相石块的运动序列3 9 泥石流冲毁并淹没房屋的序列效果图4 0 山洪汹涌而下的序列效果图4 1 山体滑坡的全过程4 2 i i i 1 2 3 4 1 2 3 4 5 6 7 8 1 2 1 2 3 4 5 6 7 8 9 加n 记 1 l 1 1 2 2 2 2 2 2 2 2 3 3 4 4 4 4 4 4 4 4 4 4 4 4 图图图图图图图图图图图图图图图图图图图图图图图图图图 浙江大学硕士学位论文 表目录 表目录 表3 1不同容重天然泥石流对应的分界粒径2 3 表3 2泥浆的一些物理参数的取值2 4 i v 浙江大学硕士学位论文浙江大学研究生学位论文独创性声明 浙江大学研究生学位论文独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作及取得的 研究成果。除了文中特别加以标注和致谢的地方外，论文中不包含其他人已经发 表或撰写过的研究成果，也不包含为获得逝婆盘堂或其他教育机构的学位或 证书而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均已在论文 中作了明确的说明并表示谢意。 学位论文作者签名：扬岛乞签字日期：秒。2 年月g 日 学位论文版权使用授权书 本学位论文作者完全了解逝鎏盘鲎有权保留并向国家有关部门或机构 送交本论文的复印件和磁盘，允许论文被查阅和借阅。本人授权逝鎏盘堂可 以将学位论文的全部或部分内容编入有关数据库进行检索和传播，可以采用影 印、缩印或扫描等复制手段保存、汇编学位论文。 ( 保密的学位论文在解密后适用本授权书) 学位论文作者签名： 砌志) 乳 导师签名： 签字日期： 2 0 。g 年占月61 3 签字目期： 口加秽年痧月7 日 浙江大学硕士学位论文第l 章绪论 第1 章绪论 本章主要阐述了本论文课题选择的背景，介绍了泥石流现象给人类带来的巨 大灾难，回顾了地质学界对泥石流研究的历史和现状，以及计算机图形学在模拟 固液两相流体方面的研究进展，并简要描述了本文的工作。 1 1 课题背景 本节简要介绍了自然场景模拟以及自然灾害场景模拟的研究意义，并着重描 述了泥石流灾害带来人类的巨大灾难，由此引出本课题选题的背景和意义。 1 1 1 自然场景模拟 自然场景的真实感模拟，一直是计算机图形学研究的热点和难点，它可以广 泛应用于计算机动画、电脑游戏、影视特技、军事仿真、建筑景观设计以及虚拟 现实等方面。 自然场景，若按照内容来划分，大致可以分为以下几类【l 】： 地形地貌类：包括海洋、沙漠、湖泊、丘陵及山峰等； 植物植被类：包括各种树木、花草等； 大气现象类：包括日落日出、雨雪风暴、彩虹闪电等； 太空场景类：包括各种恒星、行星和卫星等。 前三类自然现象是我们日常生活中所常见的，而第四类自然现象也将随着人 类宇航及太空探索活动的日益频繁而越来越为我们所熟悉。 另一方面，依据人类的心理感知，自然景象还可以分为： 自然美景，如风、花、雪、月、落日彩虹、山川湖泊、碧海蓝天以及自 然界的各种奇观( 如：海市蜃楼、极光、峨嵋佛光) 等。这类场景可以 给人们美的享受，使人们在繁重的工作学习之余，获得身心的舒畅。 自然灾害，如地震、洪水、泥石流、山体滑坡、火山爆发、海啸、台风 等。这类场景则会给人们带来触目惊心的视觉冲击及心灵震撼，使人感 浙江大学硕士学位论文第l 章绪论 受到自然界不可抵御的巨大威力、悲惨的灾难场面及生命的可贵! 如最 近发生的我国四川汶川大地震( 2 0 0 8 年5 月1 2 日1 4 时2 8 分，震中位 于北纬3 1 0 度，东经1 0 3 4 度，地震震级为8 0 级) 造成了大约6 万人 死亡、2 万多人失踪的悲剧( 见图1 1 ) ，举国为之哀悼三日( 2 0 0 8 年5 月1 9 2 1 日) ! 图1 12 0 0 8 年5 月1 2 日四川汶川大地震震后场面( 源自搜狐网) 在计算机图形学领域，对许多美景型自然场景的真实感建模与绘制已取得了 较大的进展，如雪 2 - 4 1 、雨【5 j 、彩虹问、日落日出川、月夜【8 1 、云彩1 9 1 、花草 1 0 - 1 1 j 、 海市蜃楼n 2 】、极光1 3 1 和峨嵋佛光【1 4 】的真实感模拟。 但本文关注的则是自然灾害场景的真实感模拟工作，它较之风花雪月的模拟 更具实际意义。 1 1 2 自然灾害现象及其模拟 自然灾害，是人类赖以生存的自然界所发生的异常现象。其危害面广，破坏 性大，是对人民生命财产安全的最大威胁和对社会经济发展的重大制约因素。 我国幅员辽阔，南北跨越5 0 个纬度，又处于大陆性气候与海洋性气候的交 2 浙江大学硕士学位论文第l 章绪论 互地带，各大天气系统对我国都有影响，天气多变。而且我国地壳活动性大，地 质结构复杂，地形起伏大，地貌单元多，降雨量与植被分布不均匀，使地表与地 下径流变化很大。这些原因，使我国成为世界上自然灾害种类最多的国家。其中， 对我国影响最大的自然灾害有五大类： 气象灾害，如干旱、暴雨、冻雨、酸雨、雪灾、龙卷风等。 海洋灾害，如台风、海啸、赤潮等。 洪水灾害。 地质灾害，如地震、山体滑坡、泥石流等。 生物灾害，如农作物虫害、森林鼠害、蓝藻爆发成的绿潮等。 根据形成的过程，自然灾害可以分为两类： 突发性自然灾害，如火山爆发、地震、洪水、泥石流、山体滑波、台风 ( 飓风) 、沙尘暴、冰雹等。此类灾害当致灾因素的变化超过一定强度 时，就会在几天、几小时甚至几分、几秒钟内表现为灾害行为。 缓发性自然灾害，如土地沙漠化、水土流失、环境恶化等。此类灾害通 常需要几年或更长时间的发展。 由于突发性自然灾害不可预测，爆发突然，其场景拍摄往往比较困难，所以， 突发性自然灾害的场景模拟显得尤为重要，它有助于增强人们对自然灾害的认 知，而且可以指导人们在面临自然灾害时选择正确的方法逃生，并有助各级领导 采取正确的防灾救灾措施。 然而，在计算机图形学领域，有关自然灾害场景的真实感建模与绘制工作尚 不多见，但目前已有一些尝试工作，如龙卷风【1 5 】、火山云嗍和沙尘暴【1 7 】的真实 感模拟。 本文将重点讨论泥石流这一灾难性自然现象的真实感建模与绘制工作。 1 1 3 泥石流灾害现象 泥石流是山区沟谷中，由暴雨、冰雪融水等激发的，含有大量泥砂、石块的 特殊洪流。它经常发生在峡谷地区和地震火山多发区，在暴雨期具有群发性，是 山区最严重的自然灾害。 3 浙江大学硕士学位论文 第1 章绪论 泥石流的多发地带主要包括环太平洋褶皱带，阿尔卑斯喜马拉雅褶皱带， 以及欧亚大陆内部的一些褶皱山区。世界上有近5 0 多个国家存在泥石流的潜在 威胁，其中比较严重的有哥伦比亚、秘鲁、瑞士、中国和日本。我国有泥石流沟 l 万多条，大多分布在南方省市，如西藏、四川、贵州和云南等。 泥石流给人们带来的往往是毁灭性的灾难，历史上的一些比较严重的泥石流 灾难如： 1 9 7 0 年5 月，秘鲁发生7 8 级大地震，引发迄今世界上罕见的瓦斯卡兰 山特大泥石流，使秘鲁容加依城全部被毁，近7 万人丧生。 1 9 8 5 年1 1 月，哥伦比亚的鲁伊斯火山泥石流，以5 0 公里每小时的速度 冲击了近3 万平方公里的土地，其中包括城镇、农村、田地，哥伦比亚 的阿美罗城成为废墟，造成2 5 万人死亡，1 5 万家畜死亡，1 3 万人无家 可归，经济损失高达5 0 亿美元。 1 9 9 8 年5 月6 日，意大利南部那不勒斯等地突遭遇建国以来非常罕见的 泥石流灾难，造成1 0 0 多人死亡，2 0 0 0 多人无家可归。 2 0 0 6 年2 月1 7 日，菲律宾中东部莱特省因连日暴雨引发泥石流，造成 大约3 0 0 0 人失踪，见图1 2 。 图1 22 0 0 6 年2 月1 7 日菲律宾爆发严重的泥石流( 源自新华网) 2 0 0 7 年7 月1 9 日，我国云南省保山市腾冲县境内发生一起泥石流灾害， 造成2 7 人死亡，见图1 3 。 2 0 0 8 年5 月1 2 日，我国四川省汶川县发生8 0 级大地震，给当地人民 生命、财产带来巨大灾难。同时地震导致堆积在山坡上的大量泥、砂、 4 浙江大学硕士学位论文第l 章绪论 砾石等坍塌，又因连日降雨，震区泥石流灾害频发，见图1 4 。 图1 32 0 0 7 年7 月1 9 日云南省腾冲县泥石流掩埋现场( 源自新华网) 图1 4汶川大地震造成泥石流掩埋村庄( 2 0 0 8 年5 月1 9 日源自中国新闻网) 按照物质组成特征，泥石流可以分为以下三类： 砾石型泥石流。该类泥石流的固体物质中砾石含量大而泥砂含量少，其 泥砂含量占1 0 以下。 浙江大学硕士学位论文第l 章绪论 一般型泥石流。该类泥石流的固态颗粒分布很广，含有粘土、粉土、砂、 砾石甚至巨石，其固体物质中泥砂含量约在1 0 与5 0 之间。 泥流型泥石流。该类泥石流的固体物质中泥砂含量大而砾石含量少，其 泥砂含量占5 0 以上。 陡峻的地形地貌，松软的表面，以及短时间内大量积聚的水流，是发生泥石 流的基本条件。然而泥石流往往爆发突然，很难准确地预知其发生时间与发生地 点。 因此，基于物理的泥石流现象的真实感建模与绘制，可以展示泥石流的发生、 发展和成灾的动态过程，可广泛应用于防灾救灾、影视特技以及科普教育等方面。 1 2 国内外研究现状 在地质及物理学界，已有不少对泥石流这一自然灾害现象的研究工作，但其 关注重点是泥石流的成因、运动机制及数值模拟等，所涉及的计算模型由于过于 复杂并不适合于可视模拟。在计算机图形学领域，至今尚无针对此现象的可视化 研究工作。本节将综述地质及物理学界对泥石流现象的研究历史和现状，并介绍 计算机图形学领域可用于此现象模拟的相关技术。 1 2 1 地质学界对泥石流的研究 地质学界对泥石流的研究，始于b l a c k w e l d e r 。他于1 9 2 8 年发表的论文中指 出泥石流是自然界的一种地质作用，而非上帝的行动，从而开启了地质学界对泥 石流的科学研究的先河【培】。 早期的泥石流研究，侧重于泥石流发生过程的观察、地貌现象的描述和形成 环境的分析，以定性研究为主，机理研究的内容涉及较少【协2 1 1 。 1 9 7 0 年j o h n s o n 等首次提出了泥石流的运动模型。由于泥石流的运动速率与 应力之间呈非线性比例，不能应用牛顿流体模型，他们选用宾汉粘性流模型 ( b i n g h a mv i s c o u sf l u i dm o d e l ) 建立运动方程，并求得了泥石流的最大流速瞄1 。 t a k a h a s h i 则认为泥石流并不完全属于宾汉流体，他于1 9 8 0 年提出了泥石流拜格 6 浙江大掌硕士学位论文第1 覃绪论 诺膨胀流模型( b a g n o l dd i l a t a n tf l u i dm o d e l ) 瞄】。 一般认为，宾汉模型主要适用于粘性泥石流，而拜格诺模型主要适用于稀性 泥石流。c h e r t 于1 9 8 6 年提出了一个通用的泥石流粘塑流模型( g e n e r a l i z e d v i s c o p l a s t i cm o d e l ) ，并得出了这一模型的数值解【2 4 】。1 9 9 3 年o b r i e n 等对泥石 流宾汉模型和拜格诺模型的结合作了新的尝试，建立了膨胀塑流模型( d i l a t a n t p l a s t i cm o d e l ) 【2 5 】。i v e r s o n 于1 9 9 7 年借用连续混合流理论，根据能量守恒原理， 提出了复杂的泥石流混合流理论动量守恒方程( m i x t u r et h e o r ym o m u n t u m c o n s e r v a t i o ne q u a t i o n s ) ，用以分别计算泥石流体中的固体和液体的运动速度 2 6 1 。 从以上讨论可看出，地质及物理学界的研究侧重于泥石流的运动机制，他们 建立的运动方程，相对于可视化的真实感建模与绘制工作而言，过于复杂。本文 的主要工作是基于地质学界的研究成果，通过简化泥石流的运动方程，实现其真 实感建模与绘制。下节介绍计算机图形学领域已有的可应用于对此现象进行模拟 的相关工作。 1 2 2 图形学界对流体模拟的研究 泥石流是一种典型的固液两相流。而流体模拟在近几年一直是图形学界研究 的热点。 基于物理的流体模拟，一般都是通过求解n a v i e r - s t o k e s 方程来实现的，其大 致可以分为两类：基于网格的欧拉方法和基于粒子的拉格朗日方法。 基于网格的欧拉方法，始于f o s t e r 。1 9 9 6 年f o s t e r 等首次提出通过有限差分 法在均匀网格上求解n a v i e r - s t o k e s 方程，需要设置较短的时间步长 2 7 1 。s t a m 则 于1 9 9 9 年基于半拉格朗日法实现了n a v i e r - s t o k e s 方程的稳定求解，可以设置较 长的时间步长【2 引。2 0 0 1 年f o s t e r 等提出了粒子水平集( p a r t i c l el e v e ls e t ) 方法， 用来模拟液体的运动p 】。 近几年，基于粒子的拉格朗日方法，尤其是s p h ( s m o o t h e dp a r t i c l e h y d r o d y n a m i c s 平滑粒子流体动力学) 方法，逐渐成为欧拉方法的替代。r e e v e s 于 1 9 8 3 年首次将粒子系统引入计算机图形学领域【3 0 】。虽然粒子系统可以用来实现浪 7 塑鋈奎兰堡主堂垡笙茎星! 垩堕笙 花飞溅以及泡沫的效梨3 1 越】。不过由于粒子系统没有考虑粒子间的相互作用，所 以不能用来模拟液体的运动。而s p h 方法则成功的解决了这一问题。1 9 9 5 年s t a m 等首次采用s p h 方法成功地模拟了火焰及其他气体的扩散现象【3 3 1 。2 0 0 3 年m u l l e r 等采用s p h 方法实现了交互速度级的流体模拟【3 4 1 。 因此泥石流液相部分的模拟，可以通过s p h 方法来实现。 1 2 3 图形学界对固液交互的研究 泥石流存在固液两相，两相间的相互作用影响着泥石流整体的运动。计算机 图形学界在固液交互方面也有不少研究工作。 2 0 0 4 年c a r l s o n 等实现了剐性固体和液体间的双向交互作用，不过他们仅将 刚性固体作为液体运动时的边界【3 5 】。弹性固体和液体间的交互作用也取得了一定 的进展【3 刀。2 0 0 7 年s o l e n t h a l e r 等则提出了一个统一的粒子模型来模拟固体和液 体的交互【拥。 不过，在上述研究中，固相的运动相对于液相来说是辅助，往往作为液相运 动的边界，而且并非有着确定的“流动一方向。因此，至今并无有效的方法来模 拟固液两相流的运动。而本文的工作，正是要模拟固相流体的运动，并实现固液 两相间的相互作用。 1 3 本文的工作 泥石流是介于流水与滑坡之间的一种地质作用，是一种世界范围的能给人类 带来巨大生命财产损失的灾难现象。它经常发生在峡谷地区和地震火山多发区， 在暴雨期及地震期间具有群发性，是山区最严重的自然灾害。由于泥石流爆发突 然，而且来势凶猛，对其进行场景拍摄比较困难，故对泥石流现象的真实感建模 与绘制工作具有现实意义，其研究成果可以广泛应用于防灾救灾、影视特技以及 科普教育等方面。 本文提出了一种新的基于s p h ( s m o o t h e dp a r t i c l eh v d r o d y n a m i c s 平滑粒子 流体动力学) 的固液两相流模型来有效地模拟泥石流的运动，而且，通过调节固 浙江大学硕士学位论文第l 章绪论 相体积比及其他特征参数，该模型还能够扩展到洪水、山体滑坡等灾难现象的模 拟上。 本文的主要内容包括： 本文第2 章介绍了基于物理的流体模拟的理论基础，对s p h 方法作了详尽 的阐述，为接下来提出的泥石流运动模型打下了基础。 本文第3 章针对泥石流的运动特性，提出了一种新的基于s p h 的固液两相 流模型对泥石流的运动进行建模。首先分析了泥石流的结构及运动特性，将泥石 流的运动分为四个不同阶段，并根据分界粒径将泥石流中不同半径的颗粒进行固 液二相分离归类。然后采用采用s p h 方法对液相运动进行建模，运用刚体模型( 只 存在平动和转动) 对固相运动进行约束，通过粒子模型实现了固液两相间的相互 作用，并通过g p g p u 技术加速了计算的过程。 本文第4 章实现了泥石流现象的真实感动态绘制。通过改进原有的隐式表面 定义函数，在保证绘制效果的同时，加快了对泥石流液相浆体部分表面提取的速 度。为了增强场景的真实感，本文对泥石流场景物体进行精细建模并模拟泥石流 与周围环境的交互作用，实现了泥石流冲毁淹没房屋的效果。最后，通过对固液 两相流模型的扩展，本文还实现了洪水、山体滑坡等单相流灾难现象的绘制。 在论文的最后，作者对全文的研究工作进行了总结，并提出了进一步的研究 方向和任务。 9 浙江大学硕士学位论文第2 章基于物理的流体模拟理论基础 第2 章基于物理的流体模拟理论基础 2 1 引言 流体模拟在近几年一直是计算机图形学研究的热点。早期的流体模拟，由于 受计算能力所限，主要采用参数建模的方法。p e a c h e y 于1 9 8 6 年将波浪函数表示 成一系列线性波型的组合，并采用粒子系统来模拟浪花，如图2 1 所示f 3 9 1 ，从图 中可以看出，其结果的真实感还较差。t e s s e n d o r f 则于1 9 9 9 年采用基于统计的 f f t 经验模型来描述波幅较小的海平面【加j 。 图2 1 p e a c h e y 于1 9 8 6 年模拟的波浪效果【冽 但是上述模型，用户控制起来很困难，且难以模拟一些复杂的、细节更为丰 富的效果。为解决这个问题，很多研究者转向基于物理的模拟方法。 基于物理的流体模拟，大致可以分为两类：一是基于网格的欧拉方法，一是 基于粒子的拉格朗日方法。 本章首先介绍了n a v i e r - s t o k e s 方程，它是描述流体运动的基本方程，是实现 基于物理的流体真实感模拟的前提，再依次介绍了基于网格的欧拉方法和基于粒 子的拉格朗日方法，并着重介绍了作为拉格朗日方法代表的s p h 方法，它是本文 算法的基础，正是通过改进s p h 方法，我们才实现了对泥石流现象的模拟。 1 0 浙江大学硕士学位论文第2 章基于物理的流体模拟理论基础 2 2n a v i e r - s t o k e s 方程 n a v i e r - s t o k e s 方程，是描述粘性不可压缩流体动量守恒的运动方程，简称 n s 方程，因18 2 1 年由c l a u d e - l o u i sn a v i e r 和18 4 5 年由g e o r g ec r a b f i e ls t o k e s 分别导出而得名。 n a v i e r - s t o k e s 方程可以看成是牛顿第二定律在流体上的一个应用。它的基本 形式如下【4 1 l ： v 1 ，= 0( 2 1 ) 纵鼍+ 墨卜黧+ 黑+ 乏 q 2 v n s 扭a ”酬删朗枷栅 其中，式( 2 1 ) 为质量守恒方程，式( 2 2 ) 为动量守恒方程。式中，尸为 密度，1 ，为速度，f 为时间，p 为压强，为粘滞系数，厂为外力。 首先，式( 2 1 ) 定义速度的散度为0 ，以确保流体的体积守恒。而对于不 可压缩流体而言，密度是恒定的，所以体积守恒与质量守恒是等价的，式( 2 1 ) 本质为质量守恒方程。 其次，式( 2 2 ) 左边定义了流体的加速度，其中，霎为非稳态加速度，而y v v 项为流体的对流加速度，仅与空间相关，而与时间无关。具体推导如下： 丢v n 川”= 害+ 亲妄+ 言詈+ 笔砉 加i 加加加l 广知勿昆下 2 瓦+ 【瓦万瓦j 。ll j * 4 ，二o t 西j ( 2 3 ) 钟i 缸7 匆瑟l7 西i、 = o v + o t ( v i f f x ，器卜 ：竺+ 1 ，v y = 一+ 1 ，v y 对流( c o n v e c t i o n ) 是流体内分子最常见的运动方式，主要存在以下两种形 式： 浙江大学硕士学位论文第2 章基于物理的流体模拟理论基础 扩散( d i f f u s i o n ) 。即流体内粒子的随机布朗运动。 平流( a d v e c t i o n ) 。即流体内粒子的整体运动。如水面上漂浮着的树叶， 河流中掺杂着的泥沙，都是在平流作用下前进的。 式( 2 2 ) 右边三项均为作用力项，其中一即为压力项，胛2 v 为粘滞力项， 厂为外力项，主要包括重力，浮力( 针对气体) 和表面张力( 针对液体) 等。 式( 2 2 ) 本质为牛顿第二定律在不可压缩流体上的应用。正因为如此，不 可压缩的流体，又称为牛顿流体。 另外，n a v i e r - s t o k e s 方程还有其微观粒子形式，即拉格朗日形式： 尸譬= 一即+ 2 ，+ 厂 ( 2 4 ) 口l 拉格朗日形式下的流体，完全由粒子来表示。每个粒子保持有自己的质量， 所以并不需要额外的质量守恒方程：而且，我们在描述单个粒子的运动时，并不 需要考虑对流项( 因为对流是粒子运动的宏观表现) 。这就是为什么式( 2 4 ) 比 式( 2 2 ) 缺少对流项的原因。 2 3 基于网格的欧拉方法 本节将简要介绍基于网格的欧拉方法。 欧拉方法从研究流体所占据的空间中各个固定点处的运动着手，分析被运动 流体所充满的空间中每一个固定点上的流体的速度、压强、密度等参数随时间的 变化，以及研究由某个空间点转到另一空间点时这些参数的变化。 欧拉方法将n a v i e r - s t o k e s 方程，式( 2 1 ) 和式( 2 2 ) ，离散到网格上，然后 计算各个固定网格结点上状态量的变化，从而得到整个场。在进行网格化时，一 般采用交错网格，即将标量，如压强等，分布在网格单元的中心；而将速度之类 的矢量分布在单元的表面上，这样做的好处是容易保证守恒性条件，见图2 2 。 基于网格的欧拉方法，始于f o s t e r 。1 9 9 6 年f o s t e r 等首次提出通过有限差分 法在均匀网格上求解n a v i e r - s t o k e s 方程，但需要设置较短的时间步长【2 7 1 。s i a m 则于1 9 9 9 年基于半拉格朗日法实现了n a v i e r - s t o k e s 方程的稳定求解，可以设置 1 2 浙江大学硕士学位论文第2 章基于物理的流体模拟理论基础 较长的时间步长【2 引。2 0 0 1 年f o s t e r 等提出了粒子水平集( p a r t i c l el e v e ls e t ) 方 法，用来模拟液体的运动别，见图2 3 。2 0 0 4 年l o s a s s o 等采用八叉树结构，加 速了水以及烟等流体的模拟【4 2 1 。 。- 一：。0 一t ) 一 一)一一) 一 j 一二：二i 一。三一毒一0 0 。)t = )。- j 一卜一。卜一：- a - - - - - 一卜一 o000c 哆0 4 。 。 。 ； 。 c 一，)00)，)，) 图2 2二维的欧拉网格，图中各实心点即离散的速度场【叫 图2 3 2 0 0 1 年f o s t e r 等绘制的子弹滑过水面的效果c 2 9 l 欧拉方法能够很好地描述流体的密度、压强等属性，但细节的丰富与否，取 决于网格的精细程度。八叉树等自适应网格的提出，在一定程度上增加了流体的 运动细节，降低了算法的内存消耗，但仍无法达到交互的绘制速度。f o s t e r 等在 绘制泥浆时，采用了1 5 0 2 0 0 x 1 5 0 的网格，一帧需要3 分钟；而l o s a s s o 等在绘 制水面时，采用了5 1 2 5 1 2 5 1 2 的网格，一帧也需要4 至5 分钟。为提高计算和 绘制速度，近年来，不少研究者将更多的兴趣转向了基于粒子的拉格朗日方法。 1 3 浙江大学硕士学位论文第2 章基于物理的流体模拟理论基础 2 4 基于粒子的拉格朗日方法 本节将简要介绍基于粒子的拉格朗日方法。拉格朗日方法从分析流体各个微 团的运动着手，即研究流体中某一指定微团的速度、压强、密度等参数随时间的 变化，以及研究不同流体微团参数的变化，来研究整个流体的运动。 拉格朗日方法，采用n a v i e r - s t o k e s 方程的式( 2 4 ) 来模拟流体的运动。对 于单个粒子f ： 口，：堕：兰( 2 5 )口，= o = ol z ) j j 出p t 其中，q 为粒子j 的加速度，v 为粒子l 的速度，巧为粒子j 受到的合力，肛 为粒子f 的密度。显然，这种方法是对各个相对独立的粒子进行力的分析，然后 计算这些粒子下一个时刻的位置和其他状态量。 r e e v e s 于1 9 8 3 年首次将粒子系统引入计算机图形学领域【3 0 1 。粒子系统曾被 用来模拟浪花飞溅以及泡沫的效果 3 1 - 3 2 j ，见图2 4 。 图2 4 1 9 9 5 年o b r i e n 等绘制的浪花飞溅的效果【3 1 】 拉格朗日法的优点是容易表达，不需要对整个空间进行处理，而且容易保证 质量守恒，也易于控制。但拉格朗日法不易表现平滑的流体界面，其计算量会随 着粒子数的增多而不断增大。 2 5s p h 方法 19 9 5 年s t a m 等首次将s p h ( s m o o t h e dp a r t i c l eh y d r o d ) r n a m i c s 平滑粒子流 体动力学) 方法引入计算机图形学，用来模拟火和其他气体的扩散现象。近年 来，s p h 方法逐渐成为基于粒子的拉格朗日方法的主流方法。图2 5 为m i c k y 于 1 4 釜醚 浙江大学硕士学位论文第2 章基于物理的流体模拟理论基础 2 0 0 6 年采用s p h 方法实现的浪头效果f 4 3 】。 图2 5m i e k y 于2 0 0 6 年采用s p i - i 方法实现的浪头效果降】 本节将详细介绍s p h 方法的基本思想，运动方程以及实现细节，它是本文 算法的基础。 2 5 1s p h 方法的基本思想 在拉格朗日方法中，流体被离散为若干粒子，每个粒子有着自己的质量、速 度和位置等属性。而s p h 方法的基本思想是，通过周围粒子属性的插值来计算空 间q 任一位置厂处具有连续性的某属性a ( 如密度，作用力等) ： 么( ，) = j 。么( 厂，) 形驴一厂，h ) a r ( 2 6 ) 其中，w 是半径为h 的平滑核函数。 将式( 2 6 ) 由积分形式改写为求和形式，可得： 么( ，) = a j 巧w ( r - 5 ，办) ( 2 7 ) 其中，是对所有粒子的遍历( 其实因为核函数的存在，只需遍历与厂距离 小于办的粒子) ，巧2 詈，是粒子_ ，所占据的空间，4 则是。处的属性a 。将巧展 开，式( 2 7 ) 可转换为： 舯莩4 考卿吩j ；z ) ( 2 8 ) 浙江大学硕士学位论文第2 章基于物理的流体模拟理论基础 由于求解n a v i e r - s t o k e s 方程时会用到梯度算子和拉普拉斯算子，我们还需要 得到属性a 的梯度值和拉普拉斯值。 首先，属性a 对x 分量求偏导可得： 拿彳( ，) ：呈孚矽( ，一，：，办” o xo x 。p 。 。 = - 未( a j ，- r j , h ) + 4 考昙矿c ，一，：， ，。2 9 ) _ o 矿搠+ 4 詈去驴埘 = 4 孚熹形( ，吩办) pio x 由式( 2 9 ) 知，属性么对石分量求偏导，其实就是核函数矿对z 分量求偏导， 同理，属性么对y 分量和z 分量求偏导，即核函数矿对y 分量和z 分量求偏导。 所以，属性a 的梯度值为： 黝( ，) = y j 4m 竹jv 形( ，一删 ( 2 - 1 。) 而属性a 的拉普拉斯值为： v 2 ) 2 手4 詈俨形 ( 2 2 5 2s p h 方法的平滑核函数 由2 5 1 节可知，平滑核函数w 的选择对于s p h 方法是至关紧要的。 m o n a g h a n 于1 9 9 2 年指出，平滑核函数的选择需要满足以下两个条件删： j 。形驴，h ) a r = 1 ( 2 1 2 ) 1 i m w ( r ，