低耗散且精确嵌入边界条件的流体模拟算法

第２３卷第７期　２Ｏ１１年７月　计算机辅助设计与图形学学报　Ｃｏｍｐｕｔｅｒ—Ａｉｄｅｄ　Ｄｅｓｉｇｎ＆Ｃｏｍｐｕｔｅｒ　Ｇｒａｐｈｉｃｓ　Ｖｏ１．２３　Ｎｏ．７　Ｊｕｌｙ　２０１１　低耗散且精确嵌入边界条件的流体模拟算法　杨　猛　。　，黄海明　，刘金月４　。　（中国科学院计算技术研究所智能信息处理重点实验室北京１００１９０）　（首都师范大学计算机科合研究院。（中国科学院研究生院北京北京　１０００４８｝　１０００４９）　（ｍａｘｐａｎｅ２＠１　６３．ｃｏｒｎ）　摘　要：针对欧拉流体模型中的一些弊端，提出一种基于物理的流体模拟算法．该算法采用一种改进的网格配置及　插值策略减小数值耗散，并引入一种实用的摩擦力ＶＬ￥０减少流体在固体墙处的“过度攀爬”效应；求解压力时通过离　散化等价模型隐式地嵌入Ｎｅｕｍａｎｎ边界条件，同时修改梯度算子嵌入Ｄｉｒｉｃｈｌｅｔ边界条件，从而消除了速度场的失真　现象．最后，通过对比实验展示了算法改进之处，证明了其高效、精确及稳定性．　关键词：流体模拟；数值耗散；边界条件；体积百分比；摩擦力　中图法分类号：ＴＰ３９１．９　Ｆｌｕｉｄ　Ｓｉｍｕｌａｔｉｏｎ　Ａｌｇｏｒｉｔｈｍ　ｗｉｔｈ　Ｌｏｗ　Ｄｉｓｓｉｐａｔｉｏｎ　ａｎｄ　Ａｃｃｕｒａｔｅｌｙ　Ｅｍｂｅｄｄｅｄ　Ｂｏｕｎｄａｒｙ　Ｃｏｎｄｉｔｉｏｎｓ　Ｙａｎｇ　Ｍｅｎｇ　’”．Ｈｕａｎｇ　Ｈａｉｍｉｎｇ　’。　（Ｋｅｙ　Ｌａｂｏｒａｔｏｒｙ　ｏ／Ｉｎｔｅｌｌｉｇｅｎｔ　Ｉｎｆｏｒｍａｔｉｏｎ　１ＯＯ１　９Ｏ）　Ｊｉｎｇａｎｇ　’。　Ｉｎｓｔｉｔｕｔｅ　ｏｆ　Ｃｏｍｐｕｔｉｎｇ　Ｔｅｃｈｎｏｌｏｇｙ，Ｃｈｉｎｅｓｅ　Ａｃａｄｅｍｙ　ｏｆ　Ｓｃｉｅｎｃ　ｅｓ，Ｂｅｉｊｉｎｇ　（Ｊｏｉｎ　Ｆａｃｕｌｔｙ　ｏｆ　Ｃｏｍｐｕｔｅｒ　Ｓｃｉｅｎｔｉｆｉｆ　Ｒｅｓｅａｒｃｈ，Ｃａｐｉｔａｌ　Ｎｏｒｍａｌ　Ｕｎｉｖｅｒｓｉｔｙ，Ｂｅｉｊｉｎｇ　１００１　９０）　。（Ｇｒａｄｕａｔｅ　Ｕｎｉｖｅｒｓｉｔｙ　ｏｆ　Ｃｈｉｎｅｓｅ　Ａｃａｄｅｍｙ　ｏｆ　Ｓｃｉｅｎｃｅ，Ｂｅｉｊｉｎｇ　１０００４９）　，Ａｂｓｔｒａｃｔ：Ｔｏ　ａｄｄｒｅｓｓ　ｓｏｍｅ　ｄｒａｗｂａｃｋ　ｉｎ　ｔｈｅ　Ｅｕｌｅｒ　ｂａｓｅｄ　ｆｌｕｉｄ　ｍｏｄｅｌｓ，ｔｈｉｓ　ｐａｐｅｒ　ｐｒｅｓｅｎｔｓ　ａ　ｐｈｙｓｉｃａｌｌｙ　ｂａｓｅｄ　ｆｌｕｉｄ　ｓｉｍｕｌａｔｉｏｎ　ａｌｇｏｒｉｔｈｍ．Ｔｈｅ　ｐｒｏｐｏｓｅｄ　ａｌｇｏｒｉｔｈｍ　ｒｅｄｕｃｅｓ　ｔｈｅ　ｎｕｍｅｒｉｃａｌ　ｄｉｓｓｉｐａｔｉｏｎ　ｔｈｒｏｕｇｈ　ａｎ　ｉｍｐｒｏｖｅｄ　ｇｒｉｄ　ｃｏｎｆｉｇｕｒａｔｉｏｎ　ａｎｄ　ｉｎｔｅｒｐｏｌａｔｉｏｎ　ｓｔｒａｔｅｇｙ，ａｎｄ　ｉｎｔｒｏｄｕｃｅｓ　ａ　ｐｒａｃｔｉｃａｌ　ｆｒｉｃｔｉｏｎ　ｍｅｃｈａｎｉｓｍ　ｔｏ　ｒｅｐｅｌ　ｔｈｅ“Ｅｘｃｅｓｓｉｖｅｌｙ　Ｃｌｉｍｂｉｎｇ　ｕｐ”ｅｆｆｅｃｔ　ｏｆ　ｔｈｅ　ｆｌｕｉｄ　ｎｅａｒ　ｔｈｅ　ｓｏｌｉｄ　ｗａｌｌｓ．Ｍｅａｎｗｈｉｌｅ，ｔｈｉｓ　ａｌｇｏｒｉｔｈｍ　ｄｉｓｃｒｅｔｉｚｅｓ　ａｎ　ｅｑｕｉｖａｌｅｎｔ　ｍｏｄｅｌ　ｔｏ　ｉｍｐｌｉｃｉｔｌｙ　ｅｍｂｅｄ　ｔｈｅ　Ｎｅｕｍａｎｎ　ｂｏｕｎｄａｒｙ　ｃｏｎｄｉｔｉｏｎｓ　ｆｏｒ　ｓｏｖｌｉｎｇ　ｔｈｅ　ｐｒｅｓｓｕｒｅ，ａｎｄ　ｍｏｄｉｆｉｅｓ　ｔｈｅ　ｇｒａｄｉｅｎｔ　ｏｐｅｒａｔｏｒ　ｔｏ　ｅｍｂｅｄ　ｔｈｅ　Ｄｉｒｉｃｈｌｅｔ　ｂｏｕｎｄａｒｙ　ｃｏｎｄｉｔｉｏｎ，ｗｈｉｃｈ　ｅｌｉｍｉｎａｔｅｓ　ｔｈｅ　ａｒｔｉｆａｃｔｓ　ｉｎ　ｔｈｅ　ｖｅｌｏｃｉｔｙ　ｆ　ｉｅｌｄ．Ａｔ　ｌａｓｔ，ｔｈｅ　ｅｘｐｅｒｉｍｅｎｔａｌ　ｒｅｓｕｌｔｓ　ｗｉｔｈ　ｃｏｍｐａｒｉｓｏｎｓ　ａｒｅ　ｐｒｅｓｅｎｔｅｄ　ｔｏ　ｄｅｍｏｎｓｔｒａｔｅ　ｔｈｅ　ｉｍｐｒｏｖｅｍｅｎｔｓ　ｉｎ　ｅｆｆｉｃｉｅｎｃｙ，ａｃｃｕｒａｃｙ　ａｎｄ　ｓｔａｂｉｌｉｔｙ　ｏｆ　ｔｈｅ　ｐｒｏｐｏｓｅｄ　ａｌｇｏｒｉｔｈｍ．　Ｋｅｙ　ｗｏｒｄｓ：ｆｌｕｉｄ　ｓｉｍｕｌａｔｉｏｎ；ｎｕｍｅｒｉｃａｌ　ｄｉｓｓｉｐａｔｉｏｎ；ｂｏｕｎｄａｒｙ　ｃｏｎｄｉｔｉｏｎ；ｖｏｌｕｍｅ　ｆｒａｃｔｉｏｎ；ｆｒｉｃｔｉｏｎ　基于物理的流体模拟一直是图形学领域研究的　的核心是较精确地数值求解流体控制方程——纳　维一斯托克斯方程（Ｎａｖｉｅｒ—Ｓｔｏｋｅｓ，Ｎ—Ｓ），但图形学　热点和难点，其中多数的经验与知识都源于计算流　体力学（ｃｏｍｐｕｔａｔｉｏｎａｌ　ｆｌｕｉｄ　ｄｙｎａｍｉｃｓ，ＣＦＤ）．ＣＦＤ　研究追求的则是兼顾效率与精度的优化折中．虽然　收稿日期：２Ｏ１Ｏ—Ｏ９　Ｏ８；修回日期：２０１１　Ｏ４　Ｏ７．基金项目：国家自然科学基金（６０２７３０１９，６０３７３０４２，６０４９６３２６，６０５７３０６３，６０５７３０６４）．杨猛　（１９８３　），男，博士研究生，主要研究方向为计算机图形学、基于物理的模拟；黄海明（１　９７２　），男，博士，主要研究方向为计算机图形学、计算机　动画；刘金刚（１９６３），男，博士，教授，博士生导师，主要研究方向为虚拟现实、操作系统设计．　计算机辅助设计与图形学学报　第２３卷　流体模拟的最终目的是视觉效果，但其核心及重点　仍是集中在流体建模方面，即高效、稳定地求解每个　时刻流场物理属性（如速度、压力等）的数值，而不是　集中在用于流体可视化的渲染技术上．本文给出一　种基于欧拉模型的流体模拟算法，通过修改网格配　置、修改离散算子等技术，对以往的一些研究成果进　行了改进，获得了良好的效果并提高了计算效率．　１　相关工作　基于物理的流体模拟的数值方法主要分为２　类：一种是欧拉法，即将空间划分为网格，在网格的　一些特殊位置上计算各种流场变量；另一种是拉格　朗日法，即用粒子对流体建模，追踪每个粒子的受力　状况，再根据牛顿定律来决定粒子的运动．这２类方　法优劣互补，如粒子法易于求解对流项，而欧拉法适　合求解压力项．　国内外研究人员发表了大量的对这２类方法文　献．其中Ｓｔａｒｅ＿１　引入了半拉格朗日平流法，从而大　幅提高了求解对流项的效率；Ｆｏｓｔｅｒ等　在此基础　上扩展并引入Ｌｅｖｅｌ　Ｓｅｔ方法；Ｅｎｒｉｇｈｔ等　用ＰＬＳ　（ｐａｒｔｉｃｌｅ　ｌｅｖｅｌ　ｓｅｔ）方法来精确地追踪流体表面．为　模拟小规模水流，Ｍｉｉｌｌｅｒ等ｌ＿４　用ＳＰＨ（ｓｍｏｏｔｈｅｄ　ｐａｒｔｉｃｌｅ　ｈｙｄｒｏｄｙｎａｍｉｃｓ）方法来体现倾泻、飞溅等细　节，Ｆｅｄｋｉｗ等　通过添加漩涡约束来抑制速度场旋　度的耗散，Ｌｏｓａｓｓｏ等　使用八叉树模型来增强边　界细节，Ｙａｎ等ｌ＿７　使用适应性ＳＰＨ实现实时流体模　拟，而周世哲等　８　基于多重网格法实现了实时模拟．　在流固耦合方面，Ｃａｒｌｓｏｎ等　使用了分布式拉　格朗日乘数法，Ｐｅｓｋｉｎｌ１。。引入了浸入式边界法，　Ｇｕｅｎｄｅｌｍａｎ等口　模拟了流体与轻薄物体的耦合；　而Ｋｌｉｎｇｎｅｒ等ｌ＿１　则利用非结构化网格、Ｂａｔｔｙ等ｌ＿１胡　利用动能最小化方法实现了流体一刚体的双向耦合．　Ｌｉ等　提示一种从微观角度求解Ｎ—ｓ方程的　方法——Ｌａｔｔｉｃｅ　Ｂｏｌｔｚｍａｎｎ　Ｍｅｔｈｏｄ．有关流体模拟　的综述文章详见文献Ｅ１５—１６］．　本文算法基于规则正交网格的欧拉模型，通过　抑制对流项求解产生的数值耗散，并保证对流项的　求解效率；在正交网格上处理复杂的边界条件，减小　正交网格上因几何估算导致的精度损失；解决了流　体沿固体墙的过度攀爬问题．　２流体控制方程及数值方法　不可压无粘流体的控制方程形式为　３Ｕｉ～Ｕ　ＶＵ一　ｖｐ＋ｆ　（１）　Ｄ　Ｖ・Ｕ一０　（２）　其中，ｕ表示流体速度，规定其沿　，Ｙ，　轴正方向　的分量分别为“，　，叫；』Ｄ为密度，　为压力，－厂表示体　积力（如重力、摩擦力）加速度；ｖ表示梯度算子　（Ｏ／ａｚ，ａ／ａｙ，ａ／ａ　），Ｖ表示散度算子．本文仍然采用　算子分割＿１　的思想来求解Ｎ－Ｓ方程，即先忽略式（１）　中的压力项，求解对流项和体积力项来获得一个中　间速度ｕ　，再结合体积守恒条件式（２）求解压力　场，最后利用压力场更新ｕ　，得到下一时刻的速度　场ｕ　，　一Ｕｎ．ｖ　＋ｒ，　凸　ｖ．（　）一Ｖ．Ｕ＊／　ｕ　一Ｕ＊一　ｖ　（３）　１０　ｌ０　其中，上标，ｚ表示第　时刻；Ａｔ为时间步长．式（３）　中求解压力时需要将边界条件植入线性系统（泊松　方程）中，流场的边界条件主要有２个：　一０　（４）　【，”　・”一Ｕｓ・”　（５）　式（４）规定流体自由界面处的压力值　。为０，即　Ｄｉｒｉｃｈｌｅｔ边界条件．ｕ　表示固体界面运动的速度；　式（５）规定流体不可流入固体内部，即Ｎｅｕｍａｎｎ边界　条件．如果场景中没有运动的固体，则Ｕ　Ｏ．　本文算法用改进的网格配置和插值方法来抑制　数值耗散；利用等价模型隐式嵌入Ｎｅｕｍａｎｎ边界　条件，并修改梯度算子来嵌入Ｄｉｒｉｃｈｌｅｔ边界条件，　以避免非结构化网格的开销；同时设计了有效的摩　擦力机制来抑制攀爬现象．　３本文算法模型及技术细节　３．１算法框架　本文算法定义了２种距离场：一种表示到自由　界面的有向距离值，记为　，在单元中心处计算；另　一种表示到固体界面的有向距离值，记为　。，在顶　点上采样．在每个时间步△　内，流场解算过程如下：　Ｓｔｅｐ１．平流粒子，用ＰＬＳ算法来更新距离场．　Ｓｔｅｐ２．求解对流项，并计算体积力的影响，得到ｕ　．　Ｓｔｅｐ３．将【，　从单元顶点插值到面的中心．　Ｓｔｅｐ４．求解压力场，并根据式（３）更新【，　得到　；再　根据　计算对固体边界的摩擦力．　Ｓｔｅｐ５．将【，　外推至空气域及固体域＿３］．　边界流体粒子的生成及距离场的计算均由　Ｍｏｋｂｅｒｉ提供的ＰＩ　Ｓ库［”　完成．　第７期　杨猛，等：低耗散且精确嵌入边界条件的流体模拟算法　１１７５　３．２改进的网格配置及对耗散的抑制　在经典的ＭＡＣ网格中　］，速度分量值在单元　面心处采样，如图１　ａ所示；而本文算法则是在单元　化．因此，仅对顶点处速度的增量进行滤波，就相当于更多地　保留了原来顶点速度的尖锐特征，这就是本文算法能抑制数　值耗散的理论依据．　的顶点处对整个速度向量【，采样，如图１　ｂ所示．在　顶点上计算速度场有利于累计某些源项（如漩涡力、　摩擦力等）的贡献，减少离散的粒度．　一　Ｄ　０一　●　／　，＼　ｚ＼ｌ／　／　｝／Ｂ　ａ　ＭＡＣ单元配置　ｂ改进的网格配置示例　图１面采样与顶点采样　具体步骤如下：　Ｓｔｅｐ１．平流每个粒子，粒子的速度由邻近网格点的速　度场插值计算求得．考虑到插值精度，本文不再使用线性插　值，而是使用Ｃａｔｍｕｌｌ—Ｒｏｍ三次样条插值［５］，这种：话值方法　既达到高阶精度又保证了数值的稳定性．平流所有粒子并处　理碰撞　后，用ＰＬＳ算法　来更新距离场　．　Ｓｔｅｐ２．用经典的半拉格朗日平流法　求解嬲　流项．为　提高效率，本文算法先用线性插值的方法来求“回溯”的速　度，然后在回溯到的位置采用Ｃａｔｍｕｌｌ—Ｒｏｍ插值．接着计算　，的影响【，　十一△　．　Ｓｔｅｐ３．均值计算非常直接，例如图１　ｂ中，如果用Ｘ　代　表点ｘ处中间速度的“分量，则有Ｏ　：（Ａ　＋堍＋Ｃ　＋　Ｄ　）／４．　Ｓｔｅｐ４．首先求解压力场，算法细节将在第３．３节详细论　述．接着，理论上应该先根据式（３）更新面心处的速度场；然　后通过均值运算计算顶点处的速度向量，但是反复的均值运　算会带来大量的数值耗散．为解决这个问题，本文算法将对　各个面心上的值一Ａｔ　Ｖｐ　进行均值运算，将结果作为顶点　上的对应的速度分量的增量．图２　ａ中，以Ｏ点的　分量的计　算方式为例，有　（　＋一～ＡｔＥ（ｖｐ／￣）　＋（ｖｐ／ｐ）ｒ＋（ｖｐ／ｐ）Ｇ＋（　ｏ／ｐ）Ｈ　３／４．　均值操作就像一个滤波器，速度场的尖锐特征会被平滑　ａ增量插值法示例　ｂ体积权重的计算网格示例　图２增量插值与体积权重示例　Ｓｔｅｐ５．将速度场外推至附近的空气域及固体域中，方　法同文献ＥａＪ第３．３．１节所述．　３．３　改进的压力求解及Ｎｅｕｍａｎｎ边界条件的嵌入　将Ｎｅｕｍａｎｎ边界条件植入泊松方程时，经典　的做法是以固体边界的速度为依据直接在边界处的　单元面上显式指定Ｎｅｕｍａｎｎ边界条件［２　，但是　这种方式极容易导致速度场的“阶梯化”现象．因此，　本文借助于一个等价模型隐式地处理了边界条件　Ｅａ＇ｔＰ，Ｕ　ｒｅｍｅ…Ｊ　Ｊ　Ｊ　ｎＦ　［１０　一ｕ　ｌＩ。＋　ｒ　ｒ　２Ａｔ　ｖｐ・【，Ｊ一２Ａｔ｝Ｉ　ｕｓ・ｐｎ；　Ｊ　Ｊ　ａｎＦ　其中，Ｅｏｃｔｒｅｒｎｅ表示求极值，积分域ｎ　是非固体域　的体积，ａ力　表示Ｑ　的表面，ｎ表示指向固体域的表　面单位法向量，　表示求向量的模．等价模型的证　明详见附录．本文将等价模型的面积分项通过散度　定理转为２个体积分项的和，矩阵形式的等价模型　离散式可表示为　（ｕ－ｕ　）　Ｖｕ（Ｕ－ｄ　）＋　（６）　２Ａｔ（ＶｕＧｐ＋Ｐ　Ｇ　Ｕｓ＋Ｕ　ｒＶＰＧｐ）　其中，Ｐ表示由所有离散的待求压力值组成的解向　量；Ｖ　和’，Ｐ都是对角矩阵，　每行对角线元素的　值为对应的速度分量的权重，而　每行对角线元素　的值为对应压力的权重．本文算法中，在某一点上的　流场变量值的权重定义为：以该点为几何中心、边平　行于网格方向的［０，△　］。立方体中流体所占的体积　百分比，这里的“流体”包括液体与空气．例如图２　ｂ　中，点ＰＩ处的权重为左边实线表示的立方体中流体　的体积百分比，而点“…，　处的权重为虚线表示的立　方体中流体的体积百分比；矩阵Ｇ表示标准差　分离散形式的梯度算子．对式（６）求Ｕ，Ｐ的极值并　令结果为零向量，则得到的线性系统为　ｆｌ　　／Ｇ　Ｖｕ△　Ｖ　ｕＧ　０　Ｊ１　『Ｌｐ－【，］一ｆｊ　ｌ　Ｇ　Ｖ【，　一　Ｇ　ｕ【，　／△　１（Ｊ　７）　求解式（７）所示的线性系统，方法需要谨慎选　择．由于系数矩阵是对称的，但不是对角占优，因此　使用超松弛迭代法不能保证解的收敛．然而，考虑到　系数矩阵的左上角子矩阵　是对角矩阵并且可逆，　所以可以使用高斯块消去法来消去【，，得到　ＡｔＧ　Ｇｐ：Ｇ　’／ＵＵ　十　Ｇ　Ｕｓ—Ｇ　Ｕｓ（８）　式（８）所示的线性系统中，系数矩阵是对称且　计算机辅助设计与图形学学报　第２３卷　正定的，因此可使用高效的二次型矩阵求解器ＰＣＧ　（ｐｒｅｃｏｎｄｉｔｉｏｎｅｄ　ｃｏｎｊｕｇａｔｅ　ｇｒａｄｉｅｎｔ）Ｅ１ｓ２来求解ｐ，再　根据式（３）来计算【，．　３．４　与’／Ｐ的构造　本文设计了一种递归算法来求某一个Ｅｏ，△　］。　单元中的流体体积百分比，记为　．步骤如下：　Ｓｔｅｐ１．求该偏移单元中心处的　值，若　２≥０．８６Ａｘ　（即到固体边界的有向距离大于一个Ｅ０，△　］。单元的外接球　半径），则返回　一１；否则，进行第一层分割（该单元８等分，　２Ｄ时为４等分），即深度为１的递归．　Ｓｔｅｐ２．对每个子单元，用父单元顶点处的　值插值计　算该子单元中心处的　值，如果该值大于该子单元外接球　半径，则　＋一１／８；否则，对该子单元再次分割，进行深度为　２的递归，如图３　ａ所示．对于深度为”的递归，子单元外接　球半径为（Ｏ．８６Ａｘ）／（８。』　Ｊ），而子调用返回时Ｖ的相应增量为　１／８　．如果递归深度到达最大值ｈ时，如果　＞Ｏ，则Ｖ＋一　１／８　．　该递归算法的误差量级约为０（１／８　）．　３．５　Ｄｉｒｉｃｈｌｅｔ边界条件的嵌入　式（７）所示的线性系统隐式地处理了Ｎｅｕｍａｎｎ　边界，使求解器不必再考虑固体边界处的法向量．下　面考虑Ｄｉｒｉｃｈｌｅｔ条件的嵌入：２个相邻压力点分别　为Ｌ与Ｒ，如果Ｒ处于液体域，Ｌ处于空气域，定义　一　ｌＲ　ｌ／（１　Ｒ　Ｉ＋ｌ　Ｉ）．　根据文献［１９１，自由界面处的压力可通过Ｌ和　Ｒ的值线性插值求得，且仍能达到二阶精度，即　户Ｒ（１一　）＋ｐＩ　一０　Ｉ　一ｐＲ（　一１）／０．　则相应的梯度离散式中，Ｐ　可被消去，有　一（　Ｒ—ＰＩ　）／△　—ｐＲ／ＯＡＸ．　因此，梯度算子Ｇ被修改．将修改后的Ｇ代入式（７）　中并不影响矩阵的特征．　在式（８）的线性系统中，由于Ｎｅｕｍａｎｎ条件是　隐式处理的，因此固体内部权值大于０的压力点可　能也会作为未知量进行求解．但如果固体深处权值　为０的压力点在矩阵中形成零行，则ＰＣＧ的工作效　率会受到严重影响，而且这些压力对ｕ并无影响，　因此必须防止它们成为待求未知量．本文算法针对　此问题进行了优化，即对位于固体内部的压力点，如　果它的　值处于［一０．５Ａｘ，Ｏ］之间，则算法会将它　的ｐ　值置为一０．５Ａｘ，并且在构造矩阵时进行判　断：只有　＜Ｏ的压力点有效，而对于　。＜一０．５Ａｘ　的压力点，本文用上述的Ｄｉｒｉｃｈｌｅｔ嵌入法来将它们　消去．这种方法有效地消除了矩阵的零行，保证了　ＰＣＧ的求解效率．　３．６流体与固体间的作用　如果需要计算液体对固体的作用，本文仍然应　用单向耦合ｕ列的思想，即用固体的速度作为流场的　边界条件，而流体的压力提供固体的加速．此时需要　一个线性转换算子Ｒ，Ｒｐ将固体边界附近（Ａｘ距离　内）的压力值转化为对固体的合力与力矩　ＦＭ一　０］　Ｒ—Ｑ　ＫＱ—Ｑ　ｌ　Ｄ　Ｅ一－ＪＱ，　ｒＤ　］　Ｑ—ｌ二　ｌＬ。　２　Ｊ　Ｑ３．　其中，表示固体质量的对角阵Ｍ与惯性张量Ｅ　组成了６×６质量矩阵Ｋ；Ｑ　是３×Ｎ矩阵（Ｎ为待　求压力点的个数），它的每一行分别表示３２，Ｙ，　方　向上的离散梯度；Ｑ。也是３×』＼『矩阵，它是固体质　心到所有压力点的向量的离散旋度算子；Ｑ。一Ｉ一　，　为单位矩阵．算子Ｒ和文献［１２］第３．４节所定　义的转换算子在本质上是相同的，只是鉴于网格的　结构不同，本文算法的Ｒ相当于用正交网格中的固　体体积权重代替了四面体网格的面积权重．如果流　体对固体做功，则ＡｔＫ　Ｒｐ则作为固体线速度与角　动量的增量　．　现在计算下一时刻的边界摩擦力．由于摩擦力　只影响边界处的速度场，恰好ＰＬＳ库仅在边界附近　生成粒子而固体的模型表示是三角形网格，因此本　文利用流体粒子与分布边界处的每个三角形内的３　个高斯积分点之间的作用力来为摩擦力建模，对于　某一三角形的一个高斯点Ｄ和附近的一个粒子Ｐ，　两者间的作用力为ｆ（Ｄ，Ｐ）一ｖ（Ｕ　一Ｕ　）Ｗ（ｒ，ｈ）；　其中，ｕ。，ｕ　表示Ｄ，Ｐ点下一时刻的速度，７７为摩　擦系数，ｒ为Ｄ和Ｐ两点间的距离，ｈ为影响半径．　线性衰减函数Ｈ　为　Ｗ　Ｉ（＾一ｒ），０≤ｒ≤＾　（ｒ，ｈ）一＜丁　１ｃｎ　．　ｌ０，其他　本文算法计算摩擦力场的流程如下：对每个流　体粒子，累计其附近的三角形中所有处于影响半径　内的高斯点对该粒子的合力；然后以不同的权重将　该合力的贡献分散累加至粒子所在的网格单元顶点　上（如图３　ｂ的２Ｄ示意图）．本文算法中对单元顶点　分配权重时使用的是三线性权值，由此便能得到在　顶点处的摩擦力场，该力场作为体积力，的一部分　对以后的速度场产生影响．因此本文算法在顶点上计　算速度场的优势便是防止额外力场被离散得过粗．　第7期杨猛等：低耗散且精确嵌入边界条件的流体模拟算法．子填满流域无阶梯化现象且内部和外推的速度场，，更加平滑．千广“、、。i．．，-、●●●．‘‘’，，，，，，，，●，(00)，(10)，a适应性分割的2D示例b固体边界处摩擦力机制图示图32D体积权重计算示例及摩擦力模型图62种雁力求解万式对比4实验结果及讨论本文的实验环境为Windo，图ws7显示了椭球型固体与水流相向运动的场，XP；平台配置为，景如果采用文献[12]的算法取四面体数为．10。，Intel2．83GHz4CoreCPU2GB内存NVIDIA105．则速度为0．022帧，s；本文再取1．80××32来匹配GeForce9800GT4显卡渲染工具为PBRP．“J版之前的四面体数约．0．s本文算法大幅提高效帧／图50展示了矩形水体在网格分辨率，1。(]××率的原因主要在于避免了非结构化网格的重新网格化过程造成的大量计算资源的消耗且非结构化网，的容器内坍塌翻涌的场景每个时间步的计算耗2～时约3．5s。渲染方式采用体绘制图．4中左边．格求解对流项时所用的插值方法更加复杂．应用经典的网格配置¨，一；右边是应用本文算法的网，格配置与插值方法可见其中波浪的尖锐程度更高说明本文算法实现了更低的数值耗散黍鬻滋豢鬻0．一图7流体与阎体的相向运动图。12种网格配置的数值耗散对比(第一100帧)5结i吾图54展示了．个相同分辨率的场景每个时间，．步约2～3S，渲染方式采用面绘制其中右图的场6本文分析了经典欧拉流体模型中的不足之处、，景应用了第23．节中的摩擦力机制有效地抑制了，．通过改进网格配置和插值策略利州等价模型隐式处理Neumann左图左上边界处的攀爬现象图Ax．5中取印一0．9，h边界条件达到了在正交网格上保，一持低耗散且精确处理不规则边界的目的并提出了，边界摩擦力的建模方案我们未来的工作是利用．GPU的高并行性实现加速计算并研究更精良的表，．面追踪机制以及多相流体与固体的耦合方案参考文献(R图5eferences)：添加摩擦力的效果对比(第125帧)[1]SAiam图为6显示了在倾斜固体边界下(粗黑线所示)的I。J．Stanblefluidsfe，CoEC]／／ermputerGraphicsP．roceedingYors，速度场(箭头所示)实时()penG50×50．绘制网格分辨率．nnualCoressreneeSies．ACMSIGGRAPHNewk：ACMP1999e：121128图，6中左图是显式设定，Neumann“边界”[2]FCosterNu．FGdkiwpRs．Practicals，animationofliqnuids[qS／／ies，条件的效果可以看到右边界处粒子滞留在阶梯，ompterrahicProceedingorAnnualCoress．ferenceer上；右图应用本文算法的压力求解方法可以看到粒ACMSIGGRAPH．NewYk：ACMP2001：2330１１７８　计算机辅助设计与图形学学报　第２３卷　［３］　Ｅｎｒｉｇｈｔ　Ｄ，Ｍａｒｓｃｈｎｅｒ　Ｓ，Ｆｅｄｋｉｗ　Ｒ．　Ａｎｉｍａｔｉｏｎ　ａｎｄ　ｒｅｎｄｅｒｉｎｇ　ｏｆ　ｃｏｍＣｅｘ　ｗａｔｅｒ　ｓｕｒｆａｃｅｓ［Ｊ］．ＡＣＭ　Ｔｒａｎｓａｃｔｉｏｎｓ　ｏｎ　Ｇｒａｐｈｉｃｓ，２００２，２１（３）：７３６—７４４　Ｅ４］　Ｍａｉｌｅｒ　Ｍ，Ｃｈａｒｙｐａｒ　Ｄ，Ｇｒｏｓｓ　Ｍ．Ｐａｒｔｉｃｌｅ—ｂａｓｅｄ　ｆｌｕｉｄ　ｓｉｍｕｌａｔｉｏｎ　ｆｏｒ　ｉｎｔｅｒａｃｔｉｖｅ　ａｐｐｌｉｃａｔｉｏｎｓ［ｃ］／７Ｐｒｏｃｅｅｄｉｎｇｓ　ｏｆ　ＡＣＭ　ＳｌＧＯＲＡＰＨ／Ｅｕｒｏｇｒａｐｈｉｃｓ　Ｓｙｍｐｏｓｉｕｍ　ｏｎ　Ｃｏｍｐｕｔｅｒ　Ａｎｉｍａｔｉｏｎ．Ａｉｒｅ　ｌａ　Ｖｉｌｌｅ：Ｅｕｒｏｇｒａｐｈｉｃｓ　Ａｓｓｏｃｉａｔｉｏｎ　Ｐｒｅｓｓ，　２００３：１５４－１５９　ｒ５］Ｆｅ，ｄｋｉｗ　Ｒ，Ｓｔａｍ　Ｊ，Ｊｅｎｓｅｎ　Ｈ　Ｗ．Ｖｉｓｕａｌ　ｓｉｍｕｌａｔｉｏｎ　ｏｆ　ｓｍｏｋｅ　［ｃ］／／Ｃｏｍｐｕｔｅｒ　Ｇｒａｐｈｉｃｓ　Ｐｒｏｃｅｅｄｉｎｇｓ，Ａｎｎｕａｌ　Ｃｏｎｆｅｒｅｎｃｅ　Ｓｅｒｉｅｓ．ＡＣＭ　ＳＩＧ－（ＧＲＡＰＨ．Ｎｅｗ　Ｙｏｒｋ：ＡＣＭ　Ｐｒｅｓｓ，２００１：　１５－２２　ｒ６］Ｉ　，ｏｓａｓｓｏ　Ｆ，Ｇｉｂｏｕ　Ｆ，Ｆｅｄｋｉｗ　Ｒ．Ｓｉｍｕｌａｔｉｎｇ　ｗａｔｅｒ　ａｎｄ　ｓｍｏｋｅ　ｗｉｔｈ　ａｎ　ｏｃｔｒｅｅ　ｄａｔａ　ｓｔｒｕｃｔｕｒｅ［ｃ］／／Ｃｏｍｐｕｔｅｒ　Ｇｒａｐｈｉｃｓ　Ｐｒｏｃｅｅｄｉｎｇｓ，Ａｎｎｕａｌ　Ｃｏｎｆｅｒｅｎｃｅ　Ｓｅｒｉｅｓ，ＡＣＭ　ＳＩＧＧＲＡＰＨ．　Ｎｅｗ　Ｙｏｒｋ：ＡＩ＝Ｍ　Ｐｒｅｓｓ，２００４：４５８～４６１　［７］Ｙａｎ　Ｈ，Ｗａｎｇ　Ｚ　Ｙ，Ｈｅ　Ｊ，ｅｔ　ａ１．Ｒｅａｌ—ｔｉｍｅ　ｆｌｕｉｄ　ｓｉｍｕｌａｔｉｏｎ　ｗｉｔｈ　ａｄａｐｔｉｖｅ　ＳＰＨ［Ｊ］．Ｃｏｍｐｕｔｅｒ　Ａｎｉｍａｔｉｏｎ　ａｎｄ　Ｖｉｒｔｕａｌ　Ｗｏｒｌｄｓ，２００９，２０（３）：４１７—４２６　［８］Ｚｈｏｕ　Ｓｈｉｚｈｅ，Ｍａｎ　Ｊｉａｊｕ．Ｒｅａｌｔｉｍｅ　ｆｌｕｉｄ　ｓｉｍｕｌａｔｉｏｎ　ｂａｓｅｄ　ｏｎ　ｍｕｌｔｉｇｒｉｄ　ｍｅｔｈｏｄ［Ｊ］．Ｊｏｕｒｎａｌ　ｏｆ　Ｃｏｍｐｕｔｅｒ—Ａｉｄｅｄ　Ｄｅｓｉｇｎ＆　Ｃｏｍｐｕｔｅｒ　Ｇｒａｐｈｉｃｓ，２００７，１９（７）：９３５—９４０（ｉｎ　Ｃｈｉｎｅｓｅ）　（周世哲，满家巨．基于多重网格法的实时流体模拟［Ｊ］．计算　机辅助设计与图形学学报．２００７，１９（７）：９３５—９４０）　［９］Ｃａｒｌｓｏｎ　Ｍ，Ｍｕｃｈａ　Ｐ　Ｊ，Ｔｕｒｋ　Ｇ．Ｒｉｇｉｄ　ｆｌｕｉｄ：ａｎｉｍａｔｉｎｇ　ｔｈｅ　ｉｎｔｔｅｒｐｌａｙ　ｂｅｔｗｅｅｎ　ｒｉｇｉｄ　ｂｏｄｉｅｓ　ａｎｄ　ｆｌｕｉｄ口］．　ＡＣＭ　Ｔ＇ｒａｎｓａｃｔｉｏｎｓ　ｏｎ　Ｇｒａｐｈｉｃｓ，２００４，２３（３）：３７７－－＇＇，３８４　ｒｌＯ］Ｐｅｓｋｉｎ　Ｃ　Ｓ．Ｔｈｅ　ｉｍｍｅｒｓｅｄ　ｂｏｕｎｄａｒｙ　ｍｅｔｈｏｄ［Ｊ］．Ａｃｔａ　ＮＩｕｍｅｒｉｃａ．２００２，ｌ１（１）：４７９－５１７　ｒ１１］Ｇｍｅｎｄｅｌｍａｎ　Ｅ，Ｓｅｌｌｅ　Ａ，Ｌｏｓａｓｓｏ　Ｆ，ｅｔ　ａ１．Ｃｏｕｐｌｉｎｇ　ｗａｔｅｒ　ａｎｄ　ｓｍｏｋｅ　ｔｏ　ｔｈｉｎ　ｄｅｆｏｒｍａｂｌｅ　ａｎｄ　ｒｉｇｉｄ　ｓｈｅｌｌｓ　ＬＪ］．ＡＣＭ　Ｔｒａｎｓａｃｔｉｏｎｓ　ｏｎ　Ｇｒａｐｈｉｃｓ．２００５，２４（３）：９７３一－９８１　［１２］Ｋｌｉｎｇｎｅｒ　Ｂ　Ｍ，Ｆｅｌｄｍａｎ　Ｂ　Ｅ，Ｃｈｅｎｔａｎｅｚ　Ｎ，ｅｔ　ａ１．Ｆｌｕｉｄ　附录　本文算法的等价模型为　：　帆　［ｆ。ＩＩ【，　（１）　２Ａｔ　・ｕ］一２Ａｔ　Ｉｌ　Ｕｓ・ｐｎ　Ｊ由ｎＦ　现证明式（１）等价于隐式处理体积守恒条件和　Ｎｅｕｍａｎｎ边界条件．根据变分法，定义二元标量函数　ｆ（ｘ，ｙ）一『ｆ【　［。ｆｌＩＵ＋＿ｒｑｌ－Ｕ　２Ａｔ　（ｐ＋　ｑ２）・（ｕ＋　ｑ１）］一　２Ａｔ　Ｉ１Ｊ　ｎＦ　　Ｕｓ・（户＋Ｙｑ　２）ｎ　其中，ｑ　为任意向量，而ｑ　为任意标量．由于　在（０，０）点达　到极值，因此在（０，０）点，必有ａ＿厂／３Ｘ－－０，可以得到　ｍ　［ｐ（ｕ—ｕ　）＋△ｆＶｐ］・ｑ　一０　（２）　由于ｑ　为任意向量，因此若式（２）恒成立，被积函数中　与ｑ　相乘的部分恒为０，因此有　ａｎｉｍａｔｉｏｎ　ｗｉｔｈ　ｄｙｎａｍｉｃ　ｍｅｓｈｅｓ［Ｊ］．ＡＣＭ　Ｔｒａｎｓａｃｔｉｏｎｓ　ｏｎ　Ｇｒａｐｈｉｃｓ，２００６，２５（３）：８２０一－８２５　［１３］Ｂａｔｔｙ　Ｃ，Ｂｅｒｔａｉｌｓ　Ｆ，Ｂｒｉｄｓｏｎ　Ｒ．Ａ　ｆａｓｔ　ｖａｒｉａｔｉｏｎａｌ　ｆｒａｍｅｗｏｒｋ　ｆｏｒ　ａｃｃｕｒａｔｅ　ｓｏｌｉｄ—ｆｌｕｉｄ　ｃｏｕｐｌｉｎｇ［Ｊ］．ＡＣＭ　Ｔｒａｎｓａｃｔｉｏｎｓ　ｏｎ　Ｇｒａｐｈｉｃｓ．２００７，２６（３）：１００—１０７　ｒｌ４］Ｌｉ　Ｗ，Ｗｅｉ　Ｘ　Ｍ，ａｎｄ　Ｋａｕｆｍａｎ　Ａ．Ｉｍｐｌｅｍｅｎｔｉｎｇ　ｌａｔｔｉｃｅ　Ｂｏｌｔｚｍａｎｎ　ｃｏｍｐｕｔａｔｉｏｎ　ｏｎ　ｇｒａｐｈｉｃｓ　ｈａｒｄｗａｒｅ［Ｊ］．Ｔｈｅ　Ｖｉｓｕａｌ　Ｃｏｍｐｕｔｅｒ，２００３，１９（７／８）：４４４—４５６　［１５］Ｌｉｕ　Ｙｏｕｑｕａｎ，Ｌｉｕ　Ｘｕｅｈｕｉ，Ｚｈｕ　Ｈｏｎｇｂｉｎ，ｅｔ　ａ１．Ｐｈｙｓｉｃａｌｌｙ　ｂａｓｅｄ　ｆｌｕｉｄ　ｓｉｍｕｌａｔｉｏｎ　ｉｎ　Ｃｏｍｐｕｔｅｒ　ａｎｉｍａｔｉｏｎ［Ｊ］．Ｊｏｕｒｎａｌ　ｏｆ　Ｃｏｍｐｕｔｅ￣Ａｉｄｅｄ　Ｄｅｓｉｇｎ＆．Ｃｏｍｐｕｔｅｒ　Ｇｒａｐｈｉｃｓ，２００５，１７　（１２）：２５８１－２５８９（ｉｎ　Ｃｈｉｎｅｓｅ）　（柳有权，刘学慧，朱红斌，等．基于物理的流体模拟动画综　述［Ｊ］．计算机辅助设计与图形学学报，２００５，１７（１２）：２５８卜　２５８９）　［１６］Ｔａｎ　Ｊ，Ｙａｎｇ　Ｘ　Ｂ．Ｐｈｙｓｉｃａｌｌｙ—ｂａｓｅｄ　ｆｌｕｉｄ　ａｎｉｍａｔｉｏｎ：ａ　ｓｕｒｖｅｙ　［Ｊ］．Ｓｃｉｅｎｃｅ　ｉｎ　Ｃｈｉｎａ　Ｓｅｒｉｅｓ　Ｆ：Ｉｎｆｏｒｍａｔｉｏｎ　Ｓｃｉｅｎｃｅｓ，２００８，　５２（１）：７２３—７６９　［１　７］Ｍｏｋｂｅｒｉ　Ｅ，Ｆａｌｏｕｔｓｏｓ　Ｐ．Ａ　ｐａｒｔｉｃｌｅ　ｌｅｖｅｌ　ｓｅｔ　ｌｉｂｒａｒｙ［０Ｌ＿．　［２０１０—０９—２８］．ｈｔｔｐ：／／ｗｗｗ．ｍａｇｉｘ．ｕｃｌａ．ｅｄｕ／ｓｏｆｔｗａｒｅ．ｈｔｍｌ　ｒ１８］Ｂｒｉｄｓｏｎ　Ｒ，Ｇｒｅｉｆ　Ｃ．　Ａ　ｍｕｌｔｉ—ｐｒｅｃｏｎｄｉｔｉｏｎｅｄ　ｃｏｎｊｕｇａｔｅ　ｇｒａｄｉｅｎｔ　ａｌｇｏｒｉｔｈｍ口］．ＳＩＡＭ　Ｊｏｕｒｎａｌ　ｏｎ　Ｍａｔｒｉｘ　Ａｎａｌｙｓｉｓ　ａｎｄ　Ａｐｐｌｉｃａｔｉｏｎｓ，２００６．２７（４）：１０５６—１０６８　［１９］Ｅｎｒｉｇｈｔ　Ｄ，Ｎｇｕｙｅｎ　Ｄ，Ｇｉｂｏｕ　Ｆ，　ａ１．Ｕｓｉｎｇ　ｔｈｅ　ｐａｒｔｉｃｌｅ　ｌｅｖｅｌ　ｓｅｔ　ｍｅｔｈｏｄ　ａｎｄ　ａ　ｓｅｃｏｎｄ　ｏｒｄｅｒ　ａｃｃｕｒａｔｅ　ｐｒｅｓｓｕｒｅ　ｂｏｕｎｄａｒｙ　ｃ：ｏｎｄｉｔｉｏｎ　ｆｏｒ　ｆｒｅｅ　ｓｕｒｆａｃｅ　ｆｌｏｗｓ［ｃ］Ｊ／／Ｐｒｏｃｅｅｄｉｎｇｓ　ｏｆ　ｔｈｅ　４ｔｈ　ＡＳＭＥ－ＪＳＭＥ　Ｊｏｉｎｔ　Ｆｌｕｉｄｓ　Ｃｏｎｆｅｒｅｎｃｅ．．Ｈｏｎｏｌｕｌｕ：　ＡＳＭＥ　Ｐｒｅｓｓ，２００３：１－１５　ｒ２Ｏ］Ｈａｈｎ　ｊ　Ｋ．Ｒｅａｌｉｓｔｉｃ　ａｎｉｍａｔｉｏｎ　ｏｆ　ｒｉｇｉｄ　ｂｏｄｉｅｓ［Ｊ］．ＡＣＭ　ＳＩＧＧＲＡＰＨ　Ｃｏｍｐｕｔｅｒ　Ｇｒａｐｈｉｃｓ，１９８８，２２（４）：２９９—３０８　［２１］Ｐｈａｒｒ　Ｍ，Ｈｕｍｐｈｒｅｙｓ　Ｇ．Ｐｈｙｓｉｃａｌｌｙ　ｂａｓｅｄ　ｒｅｎｄｅｒｉｎｇ：ｆｒｏｍ　ｔｈｅｏｒｙ　ｔｏ　ｉｍｐｌｅｍｅｎｔａｔｉｏｎ［Ｍ］．Ｓａｎ　Ｆｒａｎｃｉｓｃｏ：Ｍｏｒｇａｎ　Ｋａｕｆｍａｎｎ，２００４　ｐ（Ｕ－－Ｕ　）＋Ａｔ　一ｏ　ｕ一【，　一ａｔＶｐ　（３）　同理，在（Ｏ，Ｏ）点，必有３ｆ／３ｙ＝ｏ，可以得到　△　胍　ｕ一Ⅱ。　ｑ２ｎ￣Ｏ　㈩　将式（４）第一个体积分项展开，有　ｍ　ｖｑ。・ｕ＝ｍ　Ｖ．（ｑｚｕ，一ｍ　。ｚ　Ｖ・ｕ＝　Ⅱ　ｕ—ｍ　ｑ２ｖ・ｕ　㈤　将式（５）代人式（４），得到　ｍ　ｚ（ｕ一　ｍ　・Ｕ＝０（６）　类似地，由于口　为任意标量，因此式（６）的面积分和体　积分中以ｑ　为因子的部分恒为０，因此，在固体边界ａｎ，处有　ｕ・　一Ｕｓ・，ｌ　（７）　而在流体域　内有：　Ｖ・ｕ一０　ｆ８）　综上，求解式（１）等价于求解式（３），并且满足体积守恒　条件式（８）与Ｎｅｕｍａｎｎ　界条件式（７）．