（计算机应用技术专业论文）基于gpu的实时大理石纹理仿真.pdf

摘要 大理石纹釉是一项长久以来受全世界人民所喜爱的手工艺技术。其制作 流程的基本特征是，对悬浮在基底液表层的颜料的流动方式进行操纵，然后 将得到的图案转移到纸或其他具有吸纳性质的材料上进行保存。对传统的：【 艺手段来说，设计过程和最终结果的质量都在很大程度上受到各种物理条件 的限制。 本文对大理石纹釉的制作原理进行了深入探索，并详细地解析了如何运 用计算机来自动模拟这门传统艺术。仿真的主体部分是建立在真实大理石纹 釉工艺的物理模型的基础上，即二维流体动力学仿真。在多栅解法器的辅助 下用数值方法解n a v i e 卜s t o k e s 方程，得到了一个二维流场，而色彩在该流 场中的水平对流就产生了错综复杂的图案。附加的图像锐化方案则使得最终 的艺术作品更加明亮和精致。由于大量利用了现代图形硬件g p u 所包含的 优秀并行性能，本文实现了一个真正实时的交互式设计系统，为用户提供了 一种能够自由定制制图工具，并在反复尝试中寻求最佳设计的虚拟环境。由 本系统生成的大理石纹理图案具有无缝拼贴的属性，因此能够直接被用于装 饰三维对象以及绘制场景。 关键词：大理石纹釉，非真实感绘制，流体仿真，g p u 计算，纹理合成 a b s t r a c t m a r b n gi sat r a d i t i o n a ia r t1 0 v e db yp e o p l ea o v e rt h ew o r l df o r c e n t u r j e s t h ep r o c e s so fm a r b j n gj sc h a r a c t e r j z e db yl h em a n j p u i a t j o no f f i o a t l n gp i g m e n t ，a n dt h et r a r i s f e r0 ft h er e s u i t i n gd e s i g no np a p e r o rs o m e o t h e rm a t er | a 1 w 往hl h et r a d i t i o n a lm a r b i i n g ，t h ed e s i g n i n gp r o c e s sa n d t h eq u a t yo fr e s u 代i n gd e s i g n sa r el a r g e i yc o n s t r a i n e db yv a n o u sp h y s i c a l c o n d i 铂o n s 1 nt h i sp a p e r ，w ep r o v i c i ea ne x p i o r a t i o ni n t om a m i l n 9p r o c e s sa n d s h o w sh o wt h i st r a d i t i o n a la r tc a nb es i m u i a t e da i g o m h m i c a y t h em a j o r p a r ta fo u rs i m u i a t i o ni sb u i l to nt h ep h y s i c a im o d e io ft h er e a im a r b l i n g ，le a2 d 钉u i dd y n a m i c ss i m u i a t i o n t h i se n a b i e su st 0g e n e r a t er e a i i s t i c m a r b l i n gd e s i g n s i n t r i c a t ed e s i g n sa r eg e n e r a t e da st h er e s u i to fc o l o r a d v e c t i o ni nt h e2 df i a wf i e i d s0 b t a i n e db yn u m e r i c a ys o l v i n gt h e n a v i e 卜s t o k e se q u a t i o n so nt h eg p uw i t ham u 眭i g r i ds o i v e r a na d d i t i o n a i i m a g es h a r p e n i n gs c h e m ei sa d o p t e dt om a k e t h er e s u i t i n ga r t w o r km u c h m o r eb 训i a n t t h a n k st ot h ei a r g ea m o u n to fp a f a e l j s m - 几m o d e m g f a p h i c sh a r d w a r e ，w ea r ea b i et oi m p i e m e n ta ni n t e r a c t i v ec o m p u t e r s y s t e mw h i c hp r o v i d e sa r e a l _ t i m ee n v - r o n m e n tf o ru s e r st od e s i g nt h e i r o w nv i r t u a ip a t t e r n i n gt o o i sa n de x p i o i tt h eb e s td e s i g n st h r o u g ht r i a ia n d e r r o r p r o c e d u r a l i yg e n e r a t e dm a r b i i n gt e ) ( t u r e sc a na i s ob eu s e df o r r e n d e r n 9o b j e c t sa n ds c e n e sd i r e c t l y k e y w o r d s ：m a r b l i n g ，n o np h o 幻r e a s t i cr e n d e r i n g ，f i u i ds i m u i a t i o n ，g p u c o m p u t i n g ，t e ) ( t u r es y n t h e s i s 第1 章绪论 1 1 大理石纹釉工艺 大理石纹釉工艺( m a r b l i n g ) ，是一种在纸或织物上制作彩色的流体状图案的传 统艺术。作为一项历史悠久的手工艺技术，大理石纹釉工艺几个世纪以来为全世界 人民所喜爱。这项艺术起源于11 世纪的西南亚地区( 土耳其或是波斯) ，其最初的 应用是起到装饰性的目的，同时用于官方公文的背景图案以防l e 删改和伪造。直到 1 6 世纪，这项技术被当年的十字军战士带到了西欧，成为了当时书籍装订业必不可 少的组成部分【2 】。在制作精良的书本封面内侧印上错综复杂的大理石纹釉图案，既 能够遮盖装订处的折痕、细绳以及胶水的痕迹，又可以在深色的皮革封面和内部的 白色页面中形成完美的过渡。一种古老的防盗版技术即是将大理石纹釉图案印制于 书的侧面，由于撕掉其中任何一页都会使图案看上去不连贯，从而有效地预防了书 籍的人为破坏。当时的大理石纹釉手：】：艺者甚至成立了独立于书籍装订业的专门行 会，以保证这些美妙图案的配方和制作工艺秘不外传。 随着人类历史的发展，大理石纹釉工艺在书籍制作领域的功能已经逐渐被其他 先进的装订和批量印刷技术所取代，但它的装饰作用仍然保留了下来，并且在我们 的r 常生活中随处可见，比如像框，便签本，拼贴画，墙纸，灯罩，以及礼品外包 装等等。在今天，人理石纹釉工艺已经不再是少数人赖以生存的技艺，也正是因为 如此，一些世世代代流传下来的制作技法得到了公开。作为一项简单易学，而又不 乏趣味性和创造性的艺术活动，它在近千年后的今天获得了广大的追随者和爱好者。 世界各地喜爱这门艺术的人们，正利用现代更加先进的技术和广泛的资源将其不断 地发扬光大。 1 1 1 手工制作过程 大理石纹釉工艺的制作流程概括说来就是对悬浮在基底液表层的颜料的流动方 式进行操纵，并将得到的图样印到纸或其他具有吸纳性质的材料上以保存下来。传 统的制作工序大致包含以下三个步骤【1 9 ，【2 0 】： 浙江人学硕l j 学位论文 第1 章绪论 i 准备一个浅的平底碟( 一般是方形) 。容器底部预先倒入层浓稠的基底 液，通常使用水，但是必须加入一些添加剂使其浓度增大。在液体表层喷撒各 种色彩的颜料，就得到了一幅基本图样。眼药水滴瓶，移液管，还有竹刷，都 可以作为喷散颜料的工具。竹刷的妙用就在于可以一下甩出许多大小不一的液 滴，产生类似雪花的效果。承载这些颜料的基底液与颜料间必须互不相溶，而 且要足够的浓稠，以保证颜料飘浮在其表丽。图1 1 所展示的就是基本图样生 成过程巾的一个操作。 图1 1 ：基本图样的制作 2 】。 2 用本棒，梳子或其他各种工具来制作图案。梳子在碟中的来回搅动使得彩 色的液滴被拉长变形；沿着水平面轻轻吹动则能产生水波纹的近似效果。类似 的方法还有许多。一种复杂的图案往往需要对这些液体迸行反复的操作。齿的 形状和空间分布不同的梳子能够生成截然不同的效果，同样的，工具的操纵方 式也对结果有非常大的影响。 3 当图案制作完成后，将一张纸( 或者织 丁二平底碟中，吸收飘浮在表面的色彩，然后 程序制作出来的 的被制作出来。 设计作品，每一幅都是独 物，木板，皮革等等) 小心地放置 掀起，漂洗，并且晾干。以这样的 无二的，不可能再有一幅一模一样 图1 2 展示了人们在工作室中制作图案和晾制半成品的场景。 浙江大学硕 学位论文第1 章绪论 幽1 2 ：利片j 长形的梳f 搅动液体。后面悬挂着的是处丁二晾制过程中 的二r 成品【2 0 】。 “大理石纹釉”这个名称其实带有一定的误导性。事实卜，只有极少数用最自 然的色彩制作出的简单图案比较具有大理石花纹的风格；其余的绝犬多数作品则完 全没有这种特征。艺术家们为他们独特的设计赋予了于奇百怪的名字，如孔雀，喷 泉，龟甲，箭尾，贝壳等等，形象地描述了他们所希望表达的艺术思想。在下面的 四幅作品中，我们感受到了艺术家l r i s n e v i n s 的独具匠心。 ( a ) 蓝色的太刚黑子( b 1 u es u n s p o t ) ( b ) 海扇壳( c o c “e ) ( c ) 斑马( z e b r a ) ( d ) 古式斑点( a n t i q u es p o t ) 幽l - 3 ：一些大埋石纹釉作品。图例及命名全部来自丁二【2 2 】。 浙江夫学硕上学位论文第1 章绪论 1 1 2 传统工艺的局限- 眭 虽然基本步骤看上去十分简单，一幅完美的作品并不是一朝一夕能够练就的。 在传统的大理石纹釉工艺中，人为形成的各种物理条件，比如颜料的浮动性，基底 液的粘性，甚至于周遭环境的湿度，都在很人程度上限制了设计的过程以及成品的 质量。如果在制作过程中发生了误操作，就必须重新开始整个设计过程，之前的所 有努力等于是前功尽弃。这对丁二初学者来说应该是时常经历的。而另一方面，用于 承载基底液和颜料的乎底碟的大小，决定了所完成的设计作品的大小，从而直接限 制了成品的应用范围。 1 2 本文的研究内容与意义 本文对大理石纹釉工艺进行了深入的探索，并致力于研究如何以算法模拟这门 传统艺术。接着，文章展示了个能够交互式地生成大理石纹理( 即可直接用于无 缝拼贴的大理石纹釉图案) 的计算机系统。本文提供的设计系统总共有以下四个方 面的贡献： 仿真系统将艺术， 作者从物理的局限性中释放了出来。人们能够以改变粘 滞系数的方式柬控制基底液的属性，随意搭配色彩方案，设计他们自己的绘制 工具，甚至还可以撤销或重复设计过程中的每一步，通过反复的试验发掘出最 完美的设计方案。 本系统的核心原理是基于对真实大理石纹釉工艺的物理仿真，即二维流体 的动力学仿真( 2 dn u i dd y n 姗i c ss i m u l a t i o n ) 。根据s t 锄在文章【2 8 】和【2 9 】中 提出的稳态流体算法( s t a b l e 玎u i d s ) ，采用数值方法解n a v i e r - s t o k e s 方程得到 一个二维流场( 2 dn o wf i e l d ) 。色彩在该流场中的水平对流( a d v e c t i o n ) 形成 了复杂的纹样。所有这些都使得本系统生成的设计作品更加逼真。 由于近几年图形硬件并行性能的显著提升，新的图形处理单元( g r 印h i c s p r o c e s s i n gu n i t ，g p u ) 在数值应用领域已经表现出大大高于c p u 的运算性能。 本文利用了这一优势，并采用一种运行于g p u 上的多栅( m u l t i 鲥d ) 解算器将 其近涉提升，使得解n a v i e r s t o k e s 方程的过程更加快捷。本系统能够对用户 的各种操作给出实时响应，就如同真实的大理石纹釉形成过程一样。 浙江大学顾十学位论立第l 帝绪论 最后，本文所提供的系统使大理石纹釉的设计过程集成到现代的c a d c a m 系统巾成为可能。程式化地生成的大珲石纹理能够直接被用于装饰三维对象以 及绘制场景。 1 2 1 图形硬件的优势 与c p u 的单一可编程处理器不同，g p u 拥有至少两个可编程处理器：顶点处 理器( v e n e xp r o c e s s o r ) 与片段处理器( 疔a g m e mp f o c e s s o r ) ，以及其他写定的硬件 单元。处理器与硬件单元，以及应用程序之| 1 ；| j 通过数据流( d a t af l o w ) 相互链接起 来。图1 4 是g p u 流水线的图解。 剀1 _ 4 ：g p u 硬件流水线【2 4 j 。 g p u 具有极高的并行处理性能。纹理( t e x t u r e ) 是g p u 上用于并行处理的基 本数据类型，相当于c p u 上的数组。纹理数据包含四个颜色通道，因此可以直接 储存有四个以下分量的矢量数据，或者在单张纹理中存入多组标量场。c p u 上要处 理数组的每一个元素时通常使用循环迭代的方式，而g p u 并不支持对纹理像素点 的迭代。事实上，g p u 的片段处理器能够在所有片段上同时实施完全相同的操作， 即单指令，多数据( s i n 西ei n s t n l c t i o nm l l l t i p l ed a t a ，s i m d ) 运算。 很长一段时间以来，对图形流水线的编程都停留在底层的汇编语言阶段。虽然 对底层接口的直接操作提供了最大的灵活性，但是却严重影响程序的可读性和开发 浙江人学硕上学位论文 第l 章绪论 效率。n v i d i a 公司针对这种情况开发了c g 语言，一种类似于c 的g p u 编程语 言，从而大大简化了图形硬件程序的开发过程。本文将要提供的计算机仿真系统， 其核心部分就是通过c g 语占来实现的。 1 3 学术背景 1 3 1 五花八门的艺术仿真系统 非真实感绘制技术是近年来计算机图形学研究领域的一大热门话题，尤其是针 对某些特殊的艺术形式所进行的仿真。目前运用的技术手段可大致归结为三类：( 1 ) 直接在网格或其他三维模型上进行绘制；( 2 ) 输入原始图像或照片，并将其转换成 为所要模拟的绘制风格；f 3 ) 由用户来定制输入的交互式系统。有时也将几种不同 的手段相互结合起来，从而在各自的优势和劣势中寻求一个平衡点。下面分别回顾 近几年来针对多种不同类型艺术形式的仿真所进行的研究。 硬笔绘画 硬笔绘画中的一大类就是铅笔和钢笔绘画。s o u s a 等人1 9 9 9 年的文章【2 7 】对铅 笔绘画仿真的整体构架进行建模，并做了逐层深入的剖析。首先是底层的制图材料 如笔，绘图纸，橡皮擦等；建立于其上的是笔触和标记等绘制图元；接着研究绘制 技法，如轮廓线，明暗，阴影等的生成算法；位于最高层的是各个绘图步骤间的排 序和组合过程。该仿真系统的独特之处在于提供这样的功能，将一个_ 三维模型的传 统绘制过程，从勾勒外轮廓，到阴影填充，一步步地展现出来。图1 5 是其中的一 组仿真效果。 图l - 5 ：s i e pl ，h b 铅笔，描画一致轮廓；s t e p2 ，3 h ，b 铅笔，人面积轻度明暗处理；s t 印 3 ，2 b ，3 b 铅笔，加大力度，铅笔在阴影处倾斜着色，重复三遍【2 7 】。 浙江人学硕。i 学位论文第1 章绪论 2 0 0 4 年关于铅笔绘画的研究有两个不同的侧重点。w i l s o n 等人提出当绘制对 象的几何结构十分复杂，拥有大量微小的、相i 重叠的细节时，过去的算法无法将 其区分开束。因此他们在三维绘制管线中集成了二维的图像处理机制来解决这个问 题【3 l 】。y 锄a m o t o 等人的研究【3 2 】则区别于前两者：并非基于三维模型，而是将现 有的数字图像经过一系列的分片，边界剥离，笔触生成，过滤，以及合成处理，最 终转换成为彩色铅笔画。用户可以调整区域范围和定制绘图颜色，但是该系统所能 够提供的交互程度也仅限于此。 木炭画是另一种近似于铅笔画的艺术门类。d l l r a n d 等人2 0 0 1 年讨论了如何在 用户输入笔触的辅助下将照片转化为铅笔或木炭画的风格【6 。m a j u m d e r 等人2 0 0 2 年的木炭绘制【1 8 】主要贡献在于，将图形硬件引入到了三维非真实感绘制之中。一 方而，利用图形硬件的p h o n gs h a d i n g 能力实现列比度增强的绘制效果；另一方面， 经图形硬件加速过的系统达到了实时绘制的标准。 r u d o l f 等人2 0 0 3 年的蜡笔画仿真用到了前面提过的第三种技术手段：拾取用 户的笔触，并依据一定的物理原理将其合成为图画。他们以高度场来表示绘图纸， 而将蜡笔接触纸的部分表示为二维掩码。蜡的沉积量根据蜡笔剖面，接触力和摩擦 力计算得到。之前沉积下来的蜡在蜡笔运动的过程中被拖长，拖长程度的和蜡的软 硬度等信息相关。蜡迹的绘制则是基于一个光的透射和散射的简化近似模型【2 6 】。 图1 6 是他们实现的仿真效果。 蚓l - 6 ：用户定义的笔触( 左) 所生成的蜡笔绘制效果( 右) 【2 6 】。 浙 l 大学硕上学位论文 筘l 章绪论 水彩画和印象画 c u n i s 等人的1 9 9 7 年所作的水彩画仿真【5 】，首先基于一个描述浅水的流体仿真 和光学特性的物理模型，独立生成一组半透明色料绘制效果的纹理序列。在此基础 上，他们的系统提供了非常全面的功能：交互式的水彩绘画，自动图像水彩化，以 及三维场景的非真实感绘制。 关于印象画的仿真我们分别考查了l i t 晰n o w i c z 和h e n z m a n n 的研究。这两者 的共同之处在于都选择了图像梯度的法线作为移动画刷的方向，其中莳者还特别强 调了在梯度值接近于零的区域运用一种离散的数据插值方法进行处理f 1 7 】。另一方 面，前者还着眼于视频剪辑的印象化操作，提出用光学流场( o p t i c a in o wf i e l d ) 来 增强动画各帧之间的时问连贯性。后者致力于以多种笔刷尺度和长弧形的笔触轨迹 来传达图像的细节等级( 1 e v e lo f d e 协i 1 ) ：通过改变各种参数阈值可以限制或夸大笔 触的曲率，也可以使绘图效果更加粗略或精确【1 3 】。 其他艺术形式 图1 7 是h a u s n e r 于2 0 0 1 年发表的马赛克装饰画效果图【1 2 】。依据已知图像和 用户制定的边界特征，该算法既重现了图像的本来色彩，又同时在视觉上增强了边 界效果。在下图中我们可以看到，当限制马赛克瓷片为相同大小时，得到的图案均 匀且抽象；而可变大小的装饰效果则突出了整幅图案的重心，并且在细节部位要显 得清晰许多。 ( a ) 2 0 0 0 片等大瓷片 图l - 7 ：米开朗基罗作品l y b i a i ls b y l 8 ( b ) 2 0 0 0 片_ 二种大小的瓷片 的马赛克装饰画效果c 1 2 1 。 浙江大学硕 学位论文 第1 章绪论 凯尔特编结工艺( c e l t i ck n o t w o r k ) 是凯尔特民族一l j 古老的艺术。k a p l a j l 等 人2 0 0 3 年对这项工艺技术进行了研究 1 4 】。他们的仿真系统可以设计出各种式样的 复杂绳结，并在绳结的生成过程中交织和连接其他图案，见图1 8 。是更为神奇的是， 该系统还能够自动生成包裹三维模型表面的绳结。这种绳结的逻辑结构完全正确， 也就是说，依据计算机给出的方法制作编结工艺品绝对是可行的。 国1 8 ：计算机生成的凯尔特编结作品【】4 】。 1 3 _ 2 大理石纹理仿真的研究 为本次论文设计带来灵感的是2 0 0 3 年茅晓阳等人发表的大理石纹理生成系统 a t e l i e 订【1 9 】。a t e l i e r m 的基础是传统大理石纹釉工艺的物理模型。图案制作的过 程被看作一个二维流体动力学问题，并用数值方法进行解算。虽然该系统提供了交 互接口，a t e l i e t m 所能够提供的响应速度仅仅是准实时的。用户在进行每一步操作 之前都必须预先制定好整个运行轨迹，等待系统在十几秒之后给出结果一幅静 态的效果图，而后再制定下一步的操作。与此不同的是，本文在研究过程中采纳了 【1 1 】的方法，将流体的动态仿真移植到了g p u 上，实现了一个真正实时的设计系统。 在本系统提供的虚拟环境中，用户能够体验到在液体中移动制作工具，看着图案随 着他们的动作而变化的整个过程。使用者还可以同时进行多种操作，并见到它们之 问的相互作用。 塑塑型塑型二一 兰! 兰堕堡 本文后继的篇幅将会做如下安排：第二章主要介绍大理石纹理仿真的物理基础 二维流体动力学的n a v j e 卜s t o k e s 方程的解算过程，以及针对g p u 的特性所进 行的优化；第i 章讲述了大理石纹釉图案的彤成，包括多图层上的色彩平流，颜料 的添加，以及鼹示效果的锐化处理；第四章介绍本系统提供的各项功能，并展示了 若干设计作品；最后，第五章中在总结本次设计经验的基石 | 上，对本项研究的未来 发展方向进行了讨论。 o 第2 章二维流体动力学仿真 在大理石纹釉的手工制作过程中，无论是向基底液中添加新的颜料，还是用木 棒或其他工具去搅动液体，都会引起液体的流动，并带动飘浮在其表面的各种色彩， 形成各式各样的图案。因此，只要计算出基底液在用户操作卜获得的速度场，就可 以根据颜料在流场中的水平对流，生成大理石纹釉图案。山于所有颜料的运动范围 都只局限于基底液的表面，本文实际上需要解决的就是一个二维流体的动力学仿真 问题。 如何进行高效的流体动力学演算，长久以来+ 直是计算机图形学领域的热门话 题。为了探讨这个问题曾经发表了多篇文章，其中比较具有代表性的有k a s s 等人 1 9 9 0 年发表的【1 5 】，f e d k i w 等人2 0 0 1 年发表的【7 】，以及s t a m 发表的【2 8 】和【2 9 】。 本文所采用的技术主要参考了s t a n l 于1 9 9 9 年提出的“稳态流体( s t 曲l en u i d s ) ” 方法。然而，本次设计在仿真系统的实现过程中用g p u 代替了c p u ，以期获得更 高的时间效能。 2 1 流体状态的数学表达 为了在g p u 上进行流体仿真，必须为任意时刻的流体状态选择一个合理的数 学表达方式。首先是运算域的表示。大理石纹釉的制作是在一个浅的方形平底碟中 进行的，平底碟的内部平面就是它的“运算域”，我们将其映射到一个长为m 宽为 的方形栅格上，栅格内部由m 个正方形单元格构成，其空间坐标为x = ( x ，y 1 。 流体的速度随着空间和时间的不同而变化，因而将其表示为一个二维矢量场( v e c t o r f i e l d ) ，即速度场u f x ，r 1 ，其中r 为时间变量。同理，流体内部的压力可以表示为一 个标量场( s c a l a r f i e l d ) ，即压力场p ( x ，f 1 。为了方便计算，我们在空间上将这两个 场离散化，对应于前面所述的运算域，定义各个单元格内的速度和压力各自位于它 们的中心，如下页图2 1 所示。 本文在1 2 1 节中曾经提到，g p u 的纹理数据包含( 工，y ，z ，w ) 四个颜色通道。在 图形硬件的浮点运算单元中，这四个通道是作为一个整体参与计算的。也就是说， 浙江人学硕士学位论文第2 章二维溢c 体动力学仿真 单独计算其中任何一个通道的值，和同时对四个通道进行计算，在时间上没有任何 区别。利用这一点，我们将运算域划分为四个象限，分别将它彳f j 安置到一张纹理的 四个通道中，而这张纹理在两个维度上都只有原本运算域大小的一半。运用了这样 的压缩方法之后，既从一定程度上节约了图形硬件的内存空间，更为重要的是，通 过后期在此基础上的g p u 片段程序优化，省去了大量冗余计算，从而极大地提高 了算法效率。 烈2 - l ：任意时刻流体状态的数学表达1 1 1 。箭头表示速度的方向。 2 2 不可压缩流的n a v j e 卜s t o k e s 方程 本文对于二维流体的仿真是在假设其为不可压缩( i n c o m p r e s s i b l e ) 且密度均匀 ( h o m o g e n e o u s ) 的前提条件下进行的。这两个条件的含义是，流体的密度不随时间 和时间的变化而改变。这两个条件在流体仿真中是非常常见的，尤其当我们研究的 对象为水或其他液体时，增加此项限制并不影响仿真的物理适用性。 在上一节已经定义了流体的速度场u ( x ，f ) 和压力场p ( x ，) 。如果已知在初始时 间f _ o 时的速度和压力，那么二维不可压缩流体在后继时刻的运动状态可以下面的 n a v i e r s t o k e s 方程来表示： 浙江人学硼上学位论文 第2 章一维流体动力学仿真 害。( u + v ) u 一吉即+ 心2 u + f ( s 1 ) vr u = o ( 5 2 1 其中p 表示流体的密度，p 为粘度，而f 指的是作用于流体的外力。公式( 5 1 ) 右侧 的四个项都是加速度，它们的物理意义从左至右对应如下： 水平对流( a d v e c t i o n ) 流体的速度使其将置于其巾的各种物质，如颜料等，搬运到它的其他部位。实 际上，速度场本身也是被“搬运”的物质之。 压力( p r e s s u r e ) 压力的产生是由于流体分子的相互挤压。对流体施加的外力并非瞬间传播到所 有区间，而是由靠近外力的分子推挤远离外力的分子形成加速度，从而将作用力传 播出去。 扩散( d i f 如s i o n ) 不同的流体具有不同的粘稠度。粘度高的流体，流动时受到的阻力也比较大。 而这种阻力同时也引起了动力的扩散。 外部作用力( e x t e m a lf o r c e s ) 在大理石纹理仿真中，外部作用力来自于添加颜料或搅动液体等动作。外力的 其它来源一般还包括重力和浮力等等。 n a v i e r - s t o k e s 方程的第二个等式称为连续性方程( c o n t i n u i t ye q u a t i o n ) ，其中 v u 表示速度场的散度。v - u = o 的物理含义是，在流体内部任意一块子区域的边界 上速度变化的总和为零。也就是说，公式( 5 2 ) 满足了我们之前关于流体不可压缩的 假设。 在少数特定的简单物理设定下，n a v i e r - s t o k e s 方程可以解析地计算。然而本文 关注的是软件所能达到的视觉效果，并不需要达到工程学中要求的精确程度，所以 一般的数值解法已经能够符合本研究的要求。解方程的基本流程为，随着时间的增 长，逐步更新和累加速度场。每一次对速度场的更新需要经过四步操作，接下来将 会作详细的阐述。 浙江大学硕上学位论史 第2 章一维流体动力学仿真 2 _ 2 1 水平对流 速度场的水平对流是流体速度场传送其自身的过程。为了度量这一过程，我们 将每个栅格单元看作一个粒子，并对粒子的运动速度进行更新。比较直观的做法是， 直接将每个粒子从原来所在的位置r ，沿着它在一段时问研内运动的轨迹移动到新 的位置： r ( ，+ 占，) = r ( f ) + u ( f ) 研 ( 5 3 ) 这种计算方法存在两个方面的问题：首先，在速度u ( f ) 或步长研较大的情况下， 显式方法容易导致不稳定因素：其次，受g p u 的硬件条件所限，片段程序不允许改 变进行写操作的片段对象的位置。s t 锄在 2 8 】中提出的隐式方法很好地解决了这两 个问题，其采用的运算式如下： u ( x ，f + 研) = u ( x u ( x ，f ) 新，f ) ( 5 4 ) 圈圈圜圜圈圈 阑 随一l 图2 - 2 ：速度场的水平对流。邻近x 点的四个单元的速度经过插值后成 为x 点在f + 研时刻的新速度值。 式( 5 4 ) 的意义如图2 2 所示。为了得到点x 在f + 西时刻的速度值，沿着它的运 动轨迹回溯到在研时间之前的位置x ，并将点x 处的速度拷贝到x 。x 的准确速度 通过对最接近它的四个单元格进行插值得到，即图中方框标出的四个点。显而易见， 在使用这种方法之后，无论步长研取多大的值，更新后的速度都不会超过整个速度 场在前一步时的速度上限，从而保证了算法的稳定性。下面是用c g 语言实现速度 场水平对流的代码。根据2 1 节中提到的纹理压缩方法，速度场u f x ，f 1 被分解成为 浙江大学预l 学位论文 第2 章一维流体动力学仿真 代表其水平方向分量的“和代表其竖直方向分量的v ，两张l 4 分辨率大小的纹理。 代码中的f 4 t e x r e c t ( ) 是c g 的内置函数，其功能是一次性读取纹理某个坐标点上 四个通道的值。我们还另外定义了一个通用函数f l t e x r e c t c o m p r e s s b i l e r p ( ) ，用于 从压缩后的速度纹理中自动提取最适合的四个片段并进行插值。 v o i da d v e c t v e l o c i t v ( f l o a t 2c o o r d s：w p o s ， o u tf l o a t 4w n e w：c o l o r u n i f o 堋f l o a tt i m e s t e d u n jf o r mf l o a t 2w j n d o w d i m s ， u n jf o r msa i l 】口l e r r e c tu ， u u n i f o r msa i l l 口l e r r e c tv ， v u n if o r 珊sa i l l p l e r r e c tw ) v e l o c i t y f 0 1 1 0 wt h ev e l o c i t yf i e l db a c ki nt i m o c o o r d s 木= 2 ： f l o a t 4u u = c o o r d s x x x x + f l o a t 4 ( 一0 5 0 5 c o o r d s ) ： f l 。a t 4v v = c o 。r d s y y y y + f l o a t 4 ( 一o 5 ，一o 5 o5 ，0 5 ) 一t i m e s t e p 术f 4 t e x r e c t ( v c o o r d s ) ： i n t e r p 0 1 a t ea n dw r it e t ot h eo u t p u tf r a g m e n t w n e w x = f l t e x r e c t c o m p r e s s b i l e r p ( w ，f l o a t 2 ( u u x w n e w y = f l t e x r e c t c o m p r e s s b i l e r p ( w ，f l o a t 2 ( u u ，y w n e w z = f l t e x r e c t c o m p r e s s b i l e r p ( w ， f l o a t 2 ( u u z w n e w w = f l t e x r e c t c o m p r e s s b i l e r p ( w ，f l o a t 2 ( u u w w i n d o w d i m s ) w i n d o w d i m s ) w i n d o w d j m s ) w i n d o w d j m s ) 2 2 2 粘性扩散 公式( 5 1 ) 中的粘性扩散项胛2 u 可以单独表示成为一个偏微分方程等= 内2 u 。类 似于水平对流，该方程也有显式与隐式两种解法： u ( x ，f + m ) = u ( x ，) + 曲v 2 u ( x ，f ) ( 5 5 ) ( 1 一l 谚f v 2 ) u ( x ，+ 占，) = u ( x ，f ) ( 5 ，6 ) 其中i 是单位矩阵。同上一节阐述的道理一样，显式解法在这里仍然是不稳定的， 因此必须采用式( 5 6 ) 的隐式方法。这就涉及到了解泊松方程的问题。 泊松方程指的是具有a x = b 形式的矩阵方程，其中x 为未知数矢量，b 为常数 矢量，而a 为矩阵。在本例中，该矩阵是以式( 5 6 ) 左侧的( 1 一沁，v 2 ) 形式出现的， 其中v 2 为拉普拉斯算子( l a p l a c i a no p e r a t o r ) ，其定义为： 浙江人学硕i 学位论文 第2 章二维流体动力学仿真 v z u ：粤+ 宴 。靠一c 少 ( 5 7 ) 在本文的数学模型中，正方格栅格单元在两个维度上是等长的。当指定该长度 为l 时，拉普拉斯算子可离散化为下而的表达形式： v 2 u 业裔监+ f 融) = u h l ，+ u 卜l ，+ u j + l + u ，j ( 渺) 2 ( 5 8 ) 一4 u ， 其中下标f 和，分别表示栅格的行、列坐标。用式( 5 8 ) 代入方程( 5 6 ) ，经过重新组合 得到下面的式了： 尘坐坐塑盟型 甑， 其中6 = u 螂，口= l + 4 衙，= 一坩f 。式( 5 9 ) 一般称作雅各比迭代法( j a c o b i i t e m t i o n ) ，即从最初对解的猜测u ( o ) 丌始，在每一步算出一个改良了的解u ( “，逐渐 地逼近正确的结果。j a c o b i 迭代的c g 代码如下，其中h 4 t e x r e c t h e i 曲b o r s ( ) 是我 们定义的用于读取上，下，左，右，以及当前片段值的通用函数。 c o o r d s：w p o s x n e w：c o l o r r a l p h a ， b e t a w i n d o w d j m s h u n i f o r ms a m p l e r r e c tx ，xf i e l d u n i f o r ms a m p l e r r e c tb ) bf i e l d h a l f 4x c ，x l ，x r ，x b ，x t ： h 4t e x r e c t n e jg h b o r s ( x ， c o o r d s ， j a c o b i 迭代实现起来非常简单，但它的收敛速度较慢。而其他一些复杂的方法， 如共轭梯度( c o n j u g a t eg r a d i e n t ) 法和多栅( m u l t i g r i d ) 法( 见 3 】， 8 【9 】和 1 6 】) ， 则一般能取得更高的效率。本文在实现过程中就引入了g p u 上的多栅解法，关丁 这一点将会在后面的2 3 节中有详细的讲述。 1 6 f f 陀 a a a胁叭 砣h c 堇c 墨胁眦唑罾 = 罾伽d岫 浙江人学硕士学位论史 第2 章_ 二维流体动力学仿真 2 _ 2 3 施加外力 施加外力这一步没有什么复杂的算法只要直接把外力叠加到当前的速度上 就可以了。添加外力的方向为搅动液体时经过路径的切线方向，而力的大小取决于 用户预先设定的定值。唯一需要特殊处理的情况是向基底液中加入颜料的操作：虽 然用户并没有“搅动”液体，但先滴入的液滴会被后滴下的液滴挤到一边，也可能 后者刚好滴在前者的内部而迫使其向四周蔓延。为了真实地再现这些效果，我们在 每次滴入液体时，赋给这滴液体内部范围内的所有粒子一个从中心指向外的初始作 用力。这个力的大小为一个常量，在当前实现的系统中该常量值为1 。 2 2 4 无发散投影 经过上述的三个步骤，得到了一个新的速度场w 。这个速度场是不符合连续性 疗程的，也就是说，w 是发散的。求解速度场的最后一步操作，就是将w 转化成为 一个无发散( d i v e r g e n c e 一行e e ) 的速度场u 。 c h o r i n 等人在他们的h e l m h o l t z 坩o d g e 分解定理中指出，任何一个矢量场都可 以被分解成为另外两个矢量场的和 4 】，其中一个为无发散的矢量场，另一个为标量 场的梯度( 芦a d i e n t ) ，即： w = u + 即( 5 1 0 ) 其中的p 正好对应于n a v i e r s t o k e s 方程中的压力场。根据式( 5 1 0 ) ，只要求得压力 场p ，从w 中减去其梯度，就能得到我们需要的无发散速度场。这就是所谓的“无 发敝投影”。 在式( 5 1 0 ) 的左右两侧同时应用发散度算子得到： v w = v - ( u + 丫p ) = v u + v 2 p ( 5 1 1 ) 由于公式( 5 2 ) 要求v - u = o ，式( 5 1 1 ) 可以简化为： v 2 p = v w( 5 1 2 ) 我们非常高兴地发现，只要求解这个新的泊松压力方程，就能够得到压力场p ，并 最终求得u 。 浙江大学硕十学位论文 第2 章二维流体动力学仿真 在2 2 2 节中推导出的求解泊松方程的j a c o b i 迭代法在这里仍然适用，仅仪需 要改动式( 5 9 ) 的一些参数。除了将未知数由u 改为p 之外还需要改变的有：6 = v w ， 口= 一4 ，声= 1 。另外，在程序实现过程中使用到的散度算子v 的离散形式如下： 相应的c g 代码： v 。：塑+ 堡： 缸咖 q 十f ，j 一“卜1 ， 2 ( 5 1 3 ) v o i dd jv e r g e n c e ( h a l f 2c o o r d s： w p o s ， o u th a l f 4d i v：c o l o r u n if o r mh a lf 2w jn d o w d i m s i j u n i f o r ms a m p le r r e c t u ， u u n i f o ms a 皿p 1e r r e c tv ) v h a l f 4w ，w l ，w r ，w b ，w t ： w = h 4 t e x r e c t ( u ，c o o r d s ) ： w l = h 4 t e x r e c r “，h 2 n o r 衄j p o s ( c o o r d s h a j f 2 ( 1 ，0 ) ，w i n d o w d i m s m ) w r = h 4 t e x r e c t ( u ，h 2 n o r m a l p o s ( c o o r d s + h a l f 2 ( 1 ，o ) ，w i n d o w d i m s h ) ) w l= h a l f 4 ( w l y ，w x ，w l w ，w z ) ： w r= h 8 l f 4 ( w y ，w r x ，w w ，w r z ) ； = h 4 乞e x r e c t ( v c o o r d s ) ： w b = h 4 t e x r e c t ( v ，h 2 n o r 佃a l p o s ( c o o i d s h a l f 2 ( o ，1 ) ，w i n d o w d i m s l i ) ) w t = h 4 t e x r e c t ( v ，h 2 n o r m a l p o s ( c 0 0 r d s + h a l f 2 ( o ，1 ) ，w i n d o w d i m s h ) ) w b = h a l f 4 ( w b z w ，w x y ) ： w t = h a l f 4 ( w z w ，w t x y ) ： d j v = o 5 丰( w rw l + w t w b ) ： 梯度算子v 的离散形式以及求解最终速度场u 的水平方向分量的c g 代码： 勖= ( 塞，豺( 学，学 ( 5 1 4 ) v o j ds u b t r a c t g r a d i e n t u ( h a l f 2 c o o r d s ：w p o s o u th a l f 4u w：c o l o r u n i f o r mh a l f 2w i n d o w d i m s h u n i f o r ms a m p l e r r e c t p ，p r e s s u r o u n i f o r msa i i l p l e r r e c tu ) u h a l f 4p c ，p l ，p r ：