【优秀硕士博士论文】基于navier-stokes方程的烟雾湍流效果模拟

硕士学位论文题目基于NAVIERSTOKES方程的烟雾湍流效果模拟专业计算机应用技术研究方向计算机图形学授予学位日期2014年6月30日基于NAVIERSTOKES方程的烟雾湍流效果模拟计算机应用技术专业湍流是自然界最常见的一种流体现象，它使流体产生丰富的流动细节。流体湍流的模拟被广泛地应用于系统仿真、游戏特效、电影特效等领域中。基于物理的湍流模拟十分困难，虽然它能生成逼真的效果，但是在实际模拟中，其实时性与真实性之间的矛盾还远远没有得到解决。本文即尝试着在保证效率的同时，从提高烟雾真实性入手，来模拟烟雾的湍流效果。本文以湍流为研究对象，详细分析了湍流的生成条件和原因，在充分研究相关噪声、湍流模型、湍流能量谱、湍流流动工作的基础上，建立了以网格系统和粒子系统相结合的方式来模拟湍流流动的算法，并且减少了计算湍流时粒子系统与网格系统之间的耦合，在满足湍流效果的真实性的同时保证了湍流效率。本文算法的核心思想是通过一个真实的基于统计学上的K湍流模型结合柯尔莫戈洛夫定律、高效的粒子系统和小波噪声技术来精确模拟烟雾流动的湍流细节。相比已有的算法，本文工作的优点在于减少湍流粒子系统与网格系统之间的耦合，在保证湍流流动真实性的同时也兼顾了湍流的效率，使流体湍流的实时应用成为可能。本文工作的主要内容包括（1）运用小波噪声技术替代柏林噪声。克服柏林噪声渲染场景存在的细节损失和锯齿等问题。（2）在湍流模拟中引入了K两方程湍流模型，并且使用柯尔莫戈洛夫定律定义的湍流能量谱来合成湍流。（3）最后把网格系统与粒子系统结合起来。应用K湍流模型的粒子只需要采样所在位置网格系统的速度。这样降低了网格系统与粒子系统之间的耦合，使湍流在低分辨率的网格下仍然能够顺利地产生，进而通过控制粒子数目的大小来人工控制湍流产生的规模。在保证烟雾效率的同时，通过加入湍流提高了烟雾湍流的真实性。实验结果证明，在低分辨率三维网格下，通过使用基于NAVIERSTOKES方程结合K湍流模型的粒子方法模拟出流体湍流的效果能够满足应用需求。当网格在低分辨率下发生变化时，湍流效果仍然十分稳定。在湍流效果相同时，效率也比高分辨率下网格方法有了明显的提高。关键词湍流，粒子，K模型，网格，真实性，效率SMOKETURBULENCEEFFECTSIMULATIONBASEDONNAVIERSTOKESEQUATIONMAJORCOMPUTERAPPLICATIONSTECHNOLOGYPOSTGRADUATEXIANGYUTUTORZHANGYANCITURBULENCEISTHEMOSTCOMMONPHENOMENONINTHENATUREITPRESENTSABUNDANTFLUIDDETAILSTHESIMULATIONOFTURBULENCEFLUIDISWIDELYUSEDINSYSTEMSIMULATION,COMPUTERGAMEEFFECTS,MOVIEEFFECTSANDSOMEOTHERFIELDSITISVERYDIFFICULTTOSIMULATETURBULENCEFLUIDBASEDONPHYSICALMETHODSALTHOUGHITCANPRODUCEREALEFFECTS,CONTRADICTIONBETWEENREALTIMEANDREALITYISFARBEYONDSOLVEDTHISPAPERTRIESTOIMPROVESMOKEREALITYINCONDITIONOFKEEPINGITSEFFECTSWHENSIMULATINGSMOKETURBULENCETHISPAPERMAKESASTUDYOFTURBULENCEALGORITHMS,ANALYZESITSGENERATINGCONDITIONSANDFACTORS,ESTABLISHESASIMULATIONALGORITHMCOMBINEDGRIDSYSTEMWITHPARTICLESYSTEMAFTERSUFFICIENTRESEARCHESONRELATEDNOISE,TURBULENCEMODELS,TURBULENCEENERGYSPECTRUMANDTURBULENCEFLUID,ANDREDUCESCOUPLINGBETWEENPARTICLESYSTEMANDGRIDSYSTEMWHENCOMPUTINGTURBULENCETHEREFORE,ITMEETSBOTHREALITYANDEFFICIENCYOFTURBULENCEATTHESAMETIMETHECOREALGORITHMOFTHISPAPERISTHATTOSIMULATETURBULENCEDETAILSOFSMOKEFLUIDACCURATELYBYKTURBULENCEMODELCOMBINEDWITHKOLMOGOROVLAW，EFFICIENTPARTICLESYSTEMANDWAVELETNOISETECHNOLOGYCOMPAREDWITHEXISTINGALGORITHM,THEADVANTAGESOFTHISPAPERAREASFOLLOWSITREDUCEDCOUPLINGBETWEENTURBULENCEPARTICLESYSTEMANDGRIDSYSTEM,ENSUREDREALITYANDREALTIMEOFTURBULENCEFLOWATTHESAMETIMEANDMADEITSREALTIMEAPPLICATIONOFFLUIDTURBULENCEPOSSIBLETHEMAINCONTENTOFTHISPAPERINCLUDES1DISCARDINGPERLINNOISE,ITUSESWAVLETNOISEINSTEADSOITCOULDOVERCOMEDETAILLOSSANDALIASINGOFPERLINNOISEWHENRENDERINGSCENE2LEADINGKTWOEQUATIONSTURBULENCEMODELINTOTURBULENCESIMULATION,ITSYNTHESISTURBULENCEWITHENERGYSPECTRUMDEFINEDBYKOLMOGOROVSLAW3FINALLY,ITCOMBINESGRIDSYSTEMANDPARTICLESYSTEMTOGETHERPARTICLESINKTURBULENCEMODELONLYNEEDTOSAMPLEVELOCITIESOFGRIDSYSTEMINTHEIRCORRESPONDINGPOSITIONSINTHISWAY,ITREDUCESCOUPLINGBETWEENGRIDSYSTEMANDPARTICLESYSTEMANDGENERATESTURBULENCEINLOWERGRIDRESOLUTIONSUCCESSFULLYANDTHENITCANCONTROLTURBULENCESCALEARTIFICIALLYBYDECIDINGNUMBEROFPARTICLESWHILEMAINTAININGSMOKEEFFICIENCYITIMPROVESSMOKEREALITYBYADDINGTURBULENCEEFFECTTHEEXPERIMENTALRESULTSSHOWTHATITCANMEETAPPLICATIONREQUIREMENTBASEDONUSINGNAVIERSTOKEEQUATIONCOMBINEDWITHKMODELPARTICLESWHENGRIDRESOLUTIONCHANGING,TURBULENCEEFFECTISSTILLSTABLEFRAMERATEISALSOIMPROVEDCOMPAREDWITHHIGHRESOLUTIONGRIDMETHODATTHESAMEEFFECTKEYWORDSTURBULENCE,PARTICLE,MODEL,GRID,REALITY,EFFICIENCY目录目录11绪论511课题研究背景512国内外研究现状6121历史发展6122湍流方法7123最新进展10124湍流模拟的难点1113本文主要工作1114本文组织结构122流体计算算法概述1421背景简介1422NAVIERSTOKES方程的介绍1523NAVIERSTOKES方程的数值求解16231HELMHOLTZHODGE分解17232计算平流项18233计算粘度扩散项19234初始化和边界条件20235优缺点2024湍流模型介绍21241湍流粘度恒定模型21242混合尺寸模型22243单方程模型2225小结233湍流小波噪声的算法与实现2431噪声函数在湍流模拟中的作用2432小波噪声生成步骤2633小波算法理论实现26331正交性27332上采样28333小波和下采样30334构建噪声带32335小数部分的缩放与平移36336在多维空间中的噪声36337投影噪声3734特别注意的问题38341噪声方差38342噪声求值3935小结404实时湍流算法的实现4141湍流K模型4142能量谱42421能量频谱43422湍流大小尺寸43423柯尔莫戈洛夫（KOLMOGOROV）定律4443湍流的合成44431细节湍流合成的重要性质44432湍流速度场的合成4544流场的组合4745烟雾的绘制5046实现细节51461稳定性51462初始状态52463湍流相干性5247小结545实验结果及其分析5551实验环境5552实验结果与分析55521小波噪声55522后处理方法流动效果57523对比实验一62524对比实验二64525对比实验三6653本文局限性讨论676总结与展望69参考文献71作者在读期间的科研成果74声明75致谢76图11现实天空中的湍流5图21JOSSTAM的网格半拉格朗日平流法14图22粒子的平流15图23求解NS方程示意图21图31小波噪声四个步骤25图32可改善的均匀二次B样条曲线29图33一维上的小波噪声34图342J3时一维影响函数的衰减过程36图41能量的流动43图42频谱、噪声与速度场的结合46图43湍流的合成示意图47图44算法流程示意图48图45烟雾后处理方法50图46纹理坐标T1或T2变形示意图53图51两种噪声的比较56图52后处理方法下的烟雾效果图57图53后处理方法下粒子对应时刻的速度值58图54湍流运动60图55整个烟雾场景图59图56后处理方法第100帧烟雾不同视图61图57323232普通网格方法求解NS方程三维烟雾62图58646464普通网格方法求解NS方程三维烟雾63图59不同分辨率下的烟雾64图510不同分辨率下的烟雾湍流效率走势图65图511三种烟雾不同条件下帧率对比661绪论11课题研究背景无论在自然界或者是电影中，烟雾、水和火的效果随处可见。它们不可思议的流动细节使得这些物理现象看上去十分生动有趣。如图11所示，天空中的云产生的湍流和飞机降落时产生的尾流，它们的旋涡状效果十分明显。日常生活中，物体燃烧生成的火焰就是由不均匀的热空气生成的小漩涡产生的。在火山或者是大型爆炸的烟气羽流中，可以看到不同规模大小从数米到数百米的气浪或者浓烟。在流体中，比较经典的例子就是激流了，从细小的浪花到巨大的漩涡，它们的运动极其复杂。这些复杂和炫目的运动使得它们成为了电影和电脑游戏中最具挑战性的部分。大多数流动的紊乱效果都可以被描述为它们的湍流性质。在高速运动和低粘度的流动中，流动会迅速地变得非常不稳定，一连串的漩涡运动会把稳流撕裂成湍流流动。实际上，即使流体处于很小的不稳定状态也会导致在大规模流动时产生明显的变化，所以湍流流动是难以控制和预测的，这是流体模拟面临的一个普遍问题，即需要非常高的模拟分辨率和大量的计算时间。虽然当今的电影画面可以展示出相当真实的动画场景，但是流体模拟中实时性和真实性之间的矛盾还远远没有解决。为了获得高质量的视觉效果，动画设计师们必须对实现进行严格的控制，为此他们不得不进行许多特别的操作。图11现实天空中的湍流流体每一帧的流动模拟效果通常时间非常短暂而且极大地依赖于需要人工调整的粒子系统和额外的层次渲染。在计算机图形学中，为了解决上述问题，科学工作者使用了许多方法。其中，大量的前期工作集中于减少数值粘度，因为数值粘度的存在会减弱小型湍流的效果。这可以通过流体中精确的平流来解决1。虽然这种方法已经有成功的运用，但是流体实时性和真实性之间的矛盾这个基本问题却依然没有解决。12国内外研究现状流体模拟在工程和物理相关领域已经有了相当长的历史。它们在计算机图形学中的应用也变得越来越重要。但是流体模拟在计算机图形学和计算流体力学的应用中存在着巨大的差异。因此，虽然以上两者在理论方法上十分类似，但是，它们在遇到的问题时所站的角度和解决问题的方法却完全不同。本章将总结流体模拟在计算机图形学和计算流体力学中的发展历史和最新进展，并且主要介绍流体湍流的模拟。WILCOX2和POPE3的著作对湍流理论有详细的讲解，而BRIDON4也为流体模拟在计算机图形学中的应用写了一本详细的专著。121历史发展NAVIERSTOKES（纳维斯托克斯）方程对流体力学做了详细的解释，是现代流体力学的基础。为了在数值计算中解出这些方程，必须对这些方程进行离散化处理。许多工程上的应用都是基于有限元法（FEM）5或者有限差分法（FDM）6，以及有限体积法（FVM）。其中，有限差分法可以运用于结构固定的网格（GRIDS）上，而有限元法更适合在不规则的网格（MESHES）上进行计算。如果能够提前预知流体的运动，后者可以把计算资源集中于对流体流动细节要求很高的最关键部分，但是它的缺点在于需要前期的大量计算和分离网格这个多余的步骤。MAC（MARKERANDCELL）方法7改进了有限差分方法，他把流体的速度信息存储在格子（CELL）的表面而不是中心来进行计算，这提高了计算导数的精度。尽管大多数方法都是通过直接离散化NAVIERSTOKES方程来对流体模拟进行数值计算，但是不同的方法之间存在不同的规则。基于全拉格朗日的光滑粒子流体动力学（SMOOTHEDPARTICLEHYDRODYNAMICS，SPH）方法通过粒子之间的交互并结合可压缩的内核来模拟了流体力学8。MLLER等人用这个方法解决了自由表面问题，这使得SPH在计算机图形学领域得以流行9。但是为了保证粒子在复杂场景中的稳定性，就不得不选择更小的时间步长。另一方面，离散格子玻尔兹曼方法（LATTICEBOLTZMANNMETHODS，LBM）基于的是玻尔兹曼气体动力学10，它直接在规则的网格上使用模型模拟流量和流动的碰撞。虽然这种方法在计算机图形学中并不常用，但是它们已经被成功地运用于液体的模拟中11。122湍流方法在流体动力学中，湍流是流体流动的常见状态之一，它描述的是流动区域中一种紊乱和随机性的变化过程。当流体流速很小的时候，流体的流动稳定，流线清晰可见且没有混合，这被称为层流；当流速增大时，流体开始慢慢混合，流体流线变得模糊，直到流线彻底消失，流动中形成了许多小漩涡，此时的流体就被称为湍流或者乱流。这其中包括了低动能扩散、高动能对流、压力的快速变化以及速度在时间和空间上的变化。小尺度下的湍流模拟在数值计算上付出的高额代价使得研究者不得不寻找合适的湍流模拟的模型。在这些模型中，他们把湍流流动从流体整体的流动中分离出来，在基于流体整体流动的基础上，再单独从统计学的角度来分析湍流的特点。湍流模型理论在工程的应用上已经有了很长的历史，它们通常是用来估计湍流对流体整体流动的影响。比如风力作用于机翼上或者是湍流在管道流动衰减过程中的影响。而在计算机图形学中，湍流模拟主要的目的是为了增强流体本身的细节，使得流体在视觉上更加真实。下面将简单介绍湍流的建模和湍流的合成模型。1221湍流模拟在计算流体力学（COMPUTERFLUIDDYNAMICS，CFD）中经典的流体模型是模拟湍流的粘度，即在平均流量中引入虚拟扩散的湍流12。最简单经典的湍流模型可以通过混合不同尺寸大小的模型来描述非平凡的流动1314。混合尺寸模型已经成功地应用于描述边界层的流动中，比如航空工业。但是它们依赖手工输入混合尺寸，这使得模型并不能够独立于场景之外而且它只能够模拟几种已知类型的流动。为了解决这个问题，有人提出了对流扩散模型，在不考虑具体场景的情况下，此模型用一系列微分方程来模拟湍流的动态状态。在航空领域中广泛应用的SPALARTALMARAS模型15直接用一个偏微分方程（PDE）来模拟湍流粘度。两方程模型（TWOEQUATIONMODEL），模拟了代表湍流性质的参数，比如湍流势能和湍流耗散系数。K模型16和K模型17通用性强而且很简单，所以直到现在仍然是在工程上和计算流体力学中使用最为广泛的湍流模型。但是，这些模型在预测湍流和稳定性方面严重地依赖于流体的类型。因此，JONES和CHEN于1994年对此模型进行了改进，这使得它们的稳定性有了较大幅度的提高18。经典湍流模型的类型受到如下两个基本要素的限制。1）它依赖于雷诺兹平均数（REYNOLDSAVERAGING）来分离平均流动和湍流，但是这样可能并不完善，特别是在描述诸如湍流变化和相干漩涡（COHERENTEDDIES）的时候。鉴于此，SMAGORINSKY在湍流模型19中引入了大涡模拟（LARGEEDDYSIMULATIONS，LES）。LES模型用更精确的频谱滤波器来分离平均流动和湍流，它可以直接模拟相干漩涡。LES模型最初用于气象学中描述空气的内部流动，但是现在经过改进可以用于计算流体力学的诸多应用中了20。LES模型最新的一个改进是分离了涡模拟（DETACHEDEDDYSIMULATIONSDES）21，这个模型更加适合于模拟相干漩涡。2）经典模型所提供的信息。当湍流粘度足以描述平均流动的虚拟扩散时，却并没有提供各向异性流动和能量交换的信息。雷诺兹应力传输模型（REYNOLDSSTRESSTRANSPORTMODELS）解决了上面的这个问题。它不仅模拟湍流粘度，而且还模拟了整个湍流的应力张量（TURBULENTSTRESSTENSOR）2223。在此基础上，如同通用郎之万函数（GENERALLANGEVINEQUATION）24这样的密度方程直接使用粒子方法描述了湍流的高阶统计学性质。而椭圆松弛方法（ELLIPTICRELAXATION）25却在湍流模型中用全局优化法代替了流动局部对流扩散的处理方法。虽然，所有的这些方法比经典的湍流模型在理论上拥有更高阶的预测精度，但是它们的计算复杂度也随之增加，这也导致了程序开销的增大。因此，它们主要是用于计算流体力学的研究领域和高品质的模拟仿真中。随着近些年来计算机计算能力的不断增长，LES和雷诺兹应力传输模型的应用逐渐变得广泛起来。在上世纪70年代到80年代，计算流体力学中的湍流模型的基本原理陆续被提出并且得到了较好的发展，接下来，将会结合具体的话题进行讨论，特别是在湍流比较成熟的领域。近些年来，比较流行的湍流话题有可压缩流体中的湍流模拟26，湍流转变模拟2728。其中，非常吸引人一种的趋势是用LES方法实现多尺寸湍流模拟技术29，以及用自适应尺寸来模拟湍流（SCALEADAPTIVESIMULATIONSAS）30。1222湍流合成在计算机图形学中，早期的湍流方法是用统一的场来合成湍流，再对合成的湍流放大或者对基础流体模拟中合成的湍流进行替换。JOSSTAM在1993年把柯尔莫戈洛夫频谱（KOLMOGOROVSPECTRUM）产生的无散度湍流引入了湍流合成中31。这个方法曾经被用来模拟核弹的爆炸32和火焰的燃烧33。BRIDSON提出用旋度矢量噪声场来产生无散度的速度场34。他们通过潜在的场来显式地计算物体周围的流动。但是，这些方法却没有考虑湍流空间和时间的分布以及湍流本身的动力学问题。因此，这可能会导致合成的湍流不太真实。近期，为了获得更好的湍流特性，研究者引入了更加复杂的流体统计方法。KIM用小波分解来确定局部湍流的强度35，并且使用与小波频谱匹配的旋度噪声来合成湍流。这个方法假设流体基础解法可以解出湍流动力学的方程，但是实际情况却非完全如此，比如在复杂场景或者是非常不精确的粗网格下求解湍流会存在问题。SCHECHTER36和NARAIN37用传输模型来导出湍流参数。这虽然提高了湍流动力学性质，但是他们需要极大简化能量传输模型或者是做更强的条件假设，比如用不同尺度混合模型来逼近它们的单方程（ONEEQUATION）能量模型。THEODOREKIM提出了一种把网格和粒子以及两方程模型结合的算法38。在此算法中，首先通过模拟低分辨率下的小波分解来确定湍流细节的位置和能量特性；然后，再把湍流的细节合并到大规模的流体中，形成带有湍流流动的流体。但是，此方法却无法实时模拟湍流，这需要后期一步一步演示渲染。因此，能够模拟一种实时流动的湍流是十分重要的。123最新进展流体模拟中的湍流方法和漩涡表示法是计算机图形学中很活跃的一个领域。下面主要讨论湍流合成的最新进展。当一些文章还在关注和讨论基于旋度噪声纹理和漩涡的合成时，其它合成湍流的方法也在最近几年被提了出来。CHEN提出了一种混合湍流的方法，这种方法在基于粒子系统的基础上合成湍流39。但是，该方法并不是基于罗纳德均匀NAVIERSTOKES方程（REYNOLDSAVERAGEDNAVIERSTOKESEQUATION，RANS）湍流模型，他用PDF方法在粒子系统上解出了郎之万方程40。最后，更新的速度直接作为粒子细节运动的速度，这避免了旋度噪声合成中存在的相干性以及两方程湍流模型中存在的稳定性问题。但是，通过插值随机郎之万方程速度求得的速度作为实体粒子的运动速度时，合成的运动既不是无散度的又不是空间相干的。这无法满足产生真实漩涡或者是湍流的必备条件。另外，ZHAO提出了作用于上采样网格上的随机外力来合成湍流的方法41。他可以通过预计算一个给定的能量频谱来得到无散度的外力场，这种方法不需要一系列分离的标记物。例如，合成旋度纹理需要的漩涡基元或者是纹理坐标。但是，合成的外力不能够直接作用于流场，相反，它是用于预计算的场。如果这种方法用于尺寸稍大一些的湍流中，这种暂时的不一致将会导致湍流中出现人工痕迹和人工的扰动。虽然如此，如果外力场可以在程序运行时更新，这个方法也为湍流合成提供了一种比较有吸引力的思路。THEODOREKIM提出了一种利用湍流提高存在于流体模拟中的显式空间分辨率的方法。先前的方法集中于提高湍流的潜在速度场。KIM发现当流动被自由表面湍流驱动，并且超过某个频率后，这会直接在液体表面形成波浪的效果，这样只需液体模拟的表面，进而减少了流体模拟的维度。他通过利用CLOSESTPOINTMETHOD（CPM）42技术成功地跨过了这个生成液体表面时所面临的巨大障碍。CPM方法直接对三维扩张场域进行操作，这个操作需要三维运算符，由于表面运算符不是天然的三维生成形态，因此，可以利用逆阿贝尔变换（INVERSEABELTRANSFORM）来构建一个这种可行的运算符。然后，再利用扩张场域的平流方法来减少巨大的模糊性。最后，利用两个湍流耦合方法来生成高分辨率的液体表面波浪。124湍流模拟的难点从上文湍流方法介绍和实际中所遇到的困难可以看出，湍流模拟的难点主要有以下三点1湍流噪声的选取。柏林噪声（PERLINNOISE）是计算机图形学领域中最为常用的噪声，但是它频率带宽较宽，倘若应用于湍流或者其它图像图形处理中会产生细节失真和锯齿等不好的效果。因此需要找到一种新型的噪声替代柏林噪声；2湍流流动模型和能量谱的选择。众所周知，湍流流动模型有许多种类，选择模型对湍流进行模拟直接关系到湍流流动的最终效果，而能量谱的正确选择能够准确反映湍流的能量变化。因此，需要选择正确的湍流流动模型和能量谱；3湍流模拟的具体方法。流体湍流效果模拟中实时性和真实性之间存在着的巨大矛盾。目前，在网格方法中，无法直接求解湍流流动方程，因此，改方法并不能够实时和逼真地模拟出湍流的效果。为此，后来又引入了两方程模型，但是它在实际应用中也因为过于复杂而缺少进展。网格和粒子方法是模拟流体流动的主流方法，而模拟湍流又不能够缺少两方程模型的参与。所以，需要找到一种能够把流体流动模拟与两方程模型结合起来的办法。13本文主要工作传统的湍流流动模拟方法虽然能够模拟出湍流效果，但是它在实时性和真实性之间的矛盾却没有有效地得到解决。绝大多数时候，效果的实现者不得不在这二者之间选择一边或者更有侧重点。本文尝试初步找到这两者之间的平衡点。本文主要工作包括（1）本文模拟流体使用了经典的NAVIERSTOKES方程，然后以JOSSTAM求解NAVIERSTOKES方程的经典欧拉网格法方法为基础，并且把流体的不可压缩性作为NAVIERSTOKES方程求解的假设条件，用计算机在GPU上求解出加入了限制条件NAVIERSTOKES方程，这是本文的基础。（2）运用小波噪声（WAVELETNOISE）技术替代了传统流体模拟中所使用的柏林噪声。本文所运用的小波噪声技术可以很好地解决柏林噪声存在的模拟流动的细节损失和锯齿较大等问题，克服了柏林噪声存在的缺陷。（3）在流体模拟中引入了K湍流模型这个两方程模型，它使用了与湍流粘度类似的定义，即混合长度模型，但是它却用了湍流势能项K和湍流扩散项来定义。K湍流模型可以有效地模拟小规模的湍流，这些湍流无法直接用高分辨率的速度网格直接实现，因为它与流体的计算效率和网格分辨率有直接关系。相反，用基于统计学性质的湍流来合成湍流可以克服湍流在分辨率上的问题。（4）在结合前面（1）、（2）、（3）的基础上，运用GPU强大的并行处理能力，运用拉格朗日方法把粒子系统作为载体，模拟出真实可见的湍流效果。本文使用朗格朗日粒子系统的本质是直接把拉格朗日标记方法用于湍流的模拟，这样可以直接绕开网格高分辨率的效率问题。在粒子系统中，它可以精确地通过粒子数量来控制湍流细节的模拟，从而可以直接精确地在需要的地方合成湍流。（5）本文最后使用常用的后处理方法把以上计算湍流所得到的数据进行模拟渲染和可视化。本文使用拉格朗日方法来表示湍流的细节，但是忽略了相邻湍流之间的相互影响。这大大地提高了程序的效率，这样在现有计算条件下，能够通过改变粒子的相对数量来更好地有效计算流体运动的细节。本文把湍流细节和大规模地基础流动完全分开大规模地基础流动是在低分辨率网格下进行计算的，而湍流的细节则是在粒子系统中进行计算的。而能量传输模型对大规模流动只有非常低的耦合度，这保证了低分辨率的基础流动不会对细节产生不利影响。14本文组织结构全文一共分为六个章节。第一章为绪论，从计算机图形学基础和湍流基础以及它们的发展现状和存在的问题等方面入手，主要阐述了本文研究背景、国内外研究现状、本文主要工作以及当前存在的一些问题；第二章对流体模拟以及湍流模拟的相关算法和基本原理做了比较详细的解释，以及对NAVIERSTOKES方程的一般解法以及湍流的基本模型进行了介绍，它们是本文后续工作的基础和关键；第三章首先详细阐述噪声函数和小波噪声在湍流模拟中的作用，然后解释生成湍流所需小波噪声的原理和生成的算法和过程以及需要注意的问题；第四章在第三章小波噪声的基础上，详细阐述了湍流模型、湍流能量谱、湍流的合成算法的详细流程，并且给出了湍流实现的几个关键细节；第五章是对实验结果的分析和对比以及对本文工作不足之处的讨论；第六章将对本文工作做一个全面的总结，并指出本文所用方法存在的不足之处，并对湍流模拟的进一步发展做出展望。2流体计算算法概述21背景简介在数学上模拟流体，必须要有一个正确的表达式来表示流体的瞬时状态。其中，流体的速度是最为重要的一个变量，因为它决定了流体本身以及它所携带的物质如何流动。流体的速度是关于空间和时间的变量，因此，它是一个矢量，构建一个速度场就是构建一个矢量场。欧拉方法用网格表示速度场，如图21所示。因此，流体模拟最关键的部分是及时，准确地确定当前瞬时的速度场。22NAVIERSTOKES方程的介绍DXDY0,0N,0N,N0,NUX,TTUX,TX图21JOSSTAM的网格半拉格朗日平流法1PVTFUU（21）0（22）NAVIERSTOKES方程包括了流体动力学的三个基本控制方程连续性方程、动量方程、能量方程。这三个方程都遵循基本的物理学定律1质量守恒定律；2牛顿第二定律；3能量守恒定律。它们满足粘度不可压缩流体动量守恒的运动方程。它们是由CLAUDELOUISNAVIER和GEORGEGABRIELSTOKES分别于1821年和1845年推导得出的。由于三种方程在形式上是类似的，故本文仅给出了NAVIERSTOKES的连续性方程，如公式（21）所示，它是描述流动物质运动的方程。其中，U保存的是流体的速度信息，这里的速度可以是一维、二维或者三维矢量信息；表示的是流体的密度信息，流体的密度会随着流体的流动时刻变化；V是流体的粘度，一般流体，比如水的粘度非常小，如果流体没有粘度，那么它的值为0；P表示流体的压力；F表示流体受到的任何外力，这对流体的持续流动也十分重要；表示的是梯度算子；表示的是散度算子；2表示的是拉普拉斯算子。公式（22）是流体的不可压缩的约束条件。从公式（21）可以看出等号左边是流体运动速度关于时间的加速度。而等号右边第一项是平流项，其中流体的速度可以使流体沿着流动输送物体、流体本身的密度以及其它的物理量。比如流动的水里面携带有红色的液体颜料，液体本身携带了另外的液体，这时流体本身输送了其它的流体，所以它被称为平流项。如图22所示。图22粒子的平流等号右边第二项为压力项，流体分子在运动时，倾向于彼此相互挤压和碰撞。当对流体施加外力时，外力并不会被立即地传导到整个流体中，而是依靠分子间的相互碰撞来传导力，进而产生了压力。压力是流体受力的一种，而等号右边第二项代表的是压力的加速度，称为压力项。等号右边第三项被称为扩散项，比如根据日常生活经验可以得知，蜂蜜比其它很多液体更加黏稠。因此，蜂蜜流动得十分缓慢，而自来水或者酒精却流动得很快。粘度是液体本身流动性快慢的衡量标准，粘度越大，在相同受力时，液体运动越缓慢。当没有粘度时，粘度设为0。粘度对液体运动产生阻碍，造成了动量的扩散，从而引起速度的扩散。所以，这第三项被称为扩散项。等号右边第四项被称为外力项，这是由外力作用于流体上而产生的加速度。外力有很多种，它可以是局部的，也可以是整体的。局部所施加的外力影响流体的某个特定区域，例如，空调风对指定方向空气的吹动，而比如像重力，是属于整体力，它均匀地作用于所有物体。23NAVIERSTOKES方程的数值求解JOSSTAM于1999年用MAC离散方法，引入了带有一阶压力投影的半拉格朗日平流法（SEMILAGRANGIANADVECTION）43。这是无条件稳定的，虽然这个解法在计算机图形学模拟技术中数值耗散十分严重，但是它仍然是本文的基本方法和基础引用。JOSSTAM为了把NAVIERSTOKES方程应用于计算机图形学流体模拟中，首先需要把它转换成用数值积分法更容易求解的形式。从公式（21）以及（22）中可以看出NAVIERSTOKES方程是关于速度U、粘度V、压力P和外力F的方程式。数值求解这个方程需要用到以下一些方法231HELMHOLTZHODGE分解根据矢量微积分学的基础知识，在笛卡尔坐标系上的任意矢量，都可以被分解为基本矢量分量的和。同样的，对于一个在D中的矢量场W，能够唯一地分解成为公式（23）的形式PWU（23）其中，U的散度为0，而且与D平行，即在上UN0。（关于它的证明细节，详见CHORIN和MARSDEN1993的第13节）HELMHOLTZHODGE分解已经证明了任何一个矢量场可以分解为一个无散度的矢量场和一个标量场的梯度之和，而且无散度的矢量场在边界上为0。因此，可以利用以上两个结果来求解NAVIERSTOKES方程，在每个时间点上计算平流、扩散和外力来更新速度的值。得到一个散度为非0的速度场W。但是，从NAVIERSTOKES方程中不可压缩性的假设出发，结果必须是一个散度为0的速度。因此，运用HELMHOLTZHODGE分解减去压力场的梯度就可以得到需要的散度为0的速度PU（24）根据公式（23），把散度算子应用于公式（23）两边，可以得到2W（25）因为公式（22）中0U，所以公式（25）可以简化为2P（26）公式（26）是流体压力的POISSON（泊松）方程，有被称为POISSON压力方程。这说明已知W，可以通过公式（26）得到P。所以，如公式（26）所示，还需要找到一个计算W的方法。众所周知，有两个矢量M和N，可以通过计算M和N的点积，求得矢量M在N的一个单位矢量上的投射。所以，点积就是把矢量M投射到N方向上所得到的分量。因此，使用HELMHOLTZHODGE分解定理可以定义一个投影算符P，它的意思是把一个矢量场W投影到它无散度的分量U上。把P带入方程（23）中，就可以得到PP根据P的定义，。所以，0。这样就可以简化NAVIERSTOKES方程了。把投影算符带入公式（21）得12VTFUU（27）由于0，公式（27）左边的导数也是无散度的，所以，有/PTTU。同时，P0P，所以去掉压力项得到简化的方程（28）。2VTFUU（28）从公式（28）可以看出，等号右边的公式从左到右依次为平流项、扩散项和外力项。通过这三项得到一个散度速度场W，再对它进行投影运算得到一个新的无散度的速度场U。只需要从公式（26）中解出P，然后再带入公式（24）中从W减去梯度P就行了。但是这并不是像公式（28）那样，计算它们各个分量，然后再把它们加起来，而是经过状态上的变换来合成并求得最终的解。或者说，以上每个分量都是求解最终结果的一个步骤，都是取上一个结果的场作为输入，而产生的新场作为输出。可以为每个步骤定义各自的算子，它们分别是平流项A、扩散项D、外力项F以及投影项P。SPAU（29）从公式（29）可以看出从右到左分别依次计算首先，计算平流项，然后是扩散项，接着是外力项，最后是投影项，如图23所示。232计算平流项流体本身的速度使得流体沿着流动传送流体本身以及在流体中的其它量，比如密度、压力等。如果要计算平流量，就必须更新每个坐标节点上的值。在计算平流时，JOSSTAM把每个网格单元看成是一个粒子。所以，计算平流就变成了更新一个粒子系统了，这样就会使计算变得更加容易。假设每个粒子在一个时间段T中，沿着速度场的方向移动了R距离，那么TTTRU（210）这就是著名的欧拉方法，对于显式积分的微分方程来说，它是一个简单可行的方法。但是，欧拉方法有两个问题第一是在流体模拟中，使用显式的方法计算平流会造成流体在T比较大的时候不能够保持其本身的稳定性，如果TU比单个网格单元的尺寸还要大，稳定性问题还会被“放大”。第二个是与GPU本身有关，在模拟片元程序时，它不能够“移动”正在写的片元的位置。因此，JOSSTAM使用了逆向和隐式积分法。他不是通过粒子“移动”来计算平流量，而是通过每个网格单元，逆时追踪粒子轨迹，得到它先前的位置，再复制先前位置的量（这个量可以是任何需要求解的物理量）。假设Q是需要求得的物理量，那么,QXTXTU（211）这不仅对于任何时间段和速度都是稳定的，而且，在GPU上都能够很容易地实现。由于速度传送流体本身或是在流体中的其他量。因此，这个过程就是平流。JOSSTAM把这个方法命名为半拉格朗日平流法，如图21所示。233计算粘度扩散项正如前面介绍的，液体的粘度会对它自身的流动造成阻碍，这样最终会导致速度的扩散。粘度扩散的公式如下2VTU（212）同平流一样，要找到求解这个方程的方法。最容易的方法就是显式的方法2,XTTVXTUU（213）在方程（213）中，2是方程（23）拉普拉斯算子的不连续形式。如同计算平流的显式欧拉方法一样，对于T和的值比较大时就很不稳定。根据JOSSTAM的方法，对公式（213）使用隐式解法得2,VTXTIU（214）在公式（214）中，I是单位矩阵。此公式对任何时间和粘度都是稳定的。它是速度的泊松方程。因为泊松方程式的形式是矩阵方程式AXB，这里的X是要解的数值矢量，在泊松压力方程和粘度扩散方程中，X是压力项P或速度项U，而A是一个矩阵，B是一个常数矢量。求解压力和粘度扩散时都会用到泊松方程，所以对泊松方程的求解十分关键。在这里最好的方法还是迭代法。在数值计算中，雅克比（JACOBI）迭代是最常用的方法，因此把方程（23）离散后，重写为1,1,KKIJIJIJIJIJIJXXAB（215）其中，A和是常数。这个方程需要在GPU上反复循环运行很多次，但是开销并不是很大。234初始化和边界条件任何有限域的微积分都需要考虑到边界条件。边界条件确定了积分或者微分在起点和终点上数值计算所得到的值。与此同时，任何流体都必须有一个初始条件来推动整个流体模拟的运转。对于本文来说，已经假定好了初始流体的初始速度和各处的压强。在不停地循环求解NAVIERSTOKES方程每个时间步长中，都需要求解速度项和压力项。本文使用了最简单的边界条件，把速度项和压力项在边界上都设置为0，这样保证了当烟雾湍流流动到边界后不会对边界里的烟雾湍流流动造成任何影响。235优缺点以上就是对NAVIERSTOKES方程基础解法的简单介绍。图23是求解NS方程过程的简单示意图。NAVIERSTOKES方程的基础解法，简化和省去了NAVIERSTOKES方程中某些比较复杂而又不易解出的各项，使得NAVIERSTOKES方程第一次在计算机中真外力平流扩散散度压力梯度流体模拟图23求解NS方程示意图正有了解。但是，这样求解出来的流体，由于简化和省去了一些项使得流体的流动效果看上去不真实，分层比较严重。因此，为了丰富流动的细节，使流体流动看起来更加真实，需要使用湍流来增加流体流动的细节。24湍流模型介绍湍流模型致力于预测湍流的波动性质，它以平均速度场作为基础。对于诸如罗纳德均匀NAVIERSTOKES方程等平均的模拟方法来说，雷诺兹应力张量（REYNOLDSSTRESSTENSOR）是其中一个非常重要的变量。因此当今研究湍流的很大一部分工作便会集中于雷诺兹压力张量的精确计算和其通用的算法。在本小节中，经典的湍流模型基于的是湍流粘度。它在计算流体力学工程应用上使用十分地频繁。对湍流粘度进行建模有许多不同的方法，其中绝大多数方法都是基于经验和主观上的假设。这些方法可以按照建模变量的数量和它们的完整度来分类。241湍流粘度恒定模型本文模型都是基于对湍流流体中粘度的计算而建立。最简单的建模方法是假设流体中的粘度VT为恒定，但该模型只能够模拟简单的流动。一种改进的方法是将粘度VT表示为多个不同尺寸大小的漩涡模型系数的混合。但是，此方法却仍然受到漩涡尺寸以及比例的限制并无法解决湍流中湍流尺寸不一的问题。242混合尺寸模型SMAGORINSKY、BALDWIN和LOMAX共同提出的混合尺寸模型是模拟湍流中很常见的一种模型。它是通过对边界层流中显式湍流粘度的通用化而得来的。它被定义为2|TMVL（216）其中，LM是混合尺寸长度，05/IJIJJIUX是平均流动场的旋转张量（ROTATIONTENSOR）。混合尺寸长度对前面提到的问题的几何性进行了编码，它的精确性取决于准确的尺寸值。LM的解析式只能适用于几种简单类型的湍流，比如适用于具有线性性质的边界条件中。而对于复杂的边界条件，用此方法生成湍流的运动将是不可预测的。因此，这个模型并不完整，还需要对它进行进一步地完善。243单方程模型SPALARTALLMARAS单方程模型是湍流模型中比较完整使用非常广泛的一种模型。它是可以直接解出湍流的粘度的一种传输方程，可以解决模拟一些比较复杂的流动，它在相当范围内的各种流动中使用都能够得到较好的效果。它是在分子水平上的粘度使用经验主义、维度分析理论、伽利略不变式和选择性依赖建模而成的。此模型并不依赖于较高的网格分辨率，它只要求能够精确捕捉到速度场梯度的袋鼠模型，因此，在粗网格下，以及湍流精度要求不高的情况下，计算湍流能够取得较好的效果。它的微分方程形式如公式15TJIJDJKJVKKUTCTXLXX（217）的流动。其中，粘度系数被定义为05TV（218）SPALARTALLMARAS对于大量的分离流动，无剪切流以及衰变湍流效果不好。由于是单方程模型，相对于两方程模型来说，使得使用者对模型变量的控制不能做到随心所欲。25小结本章首先对NAVIERSTOKES方程的解法和特点进行了较为详细地阐述，它是本文网格流体流动的基础。然后再对与本文相关的湍流模拟中历史上使用的几种模型做了简单的介绍，并且简单地分析了各个湍流模型的特点。为本文后面章节中引入的两方程K湍流模型做了必要的铺垫。3湍流小波噪声的算法与实现31噪声函数在湍流模拟中的作用噪声的使用在流体模拟乃至湍流模拟中都十分重要。它的作用是使在计算机图形学中模拟的运动过程显得更加接近自然界的无序过程和随机过程。由于本文需要生成湍流效果，而这本身就属于流体流动细节的范畴，这使得对噪声的要求就变得更加重要，本文使用的噪音需要能够体现湍流细节和尽量减少人工锯齿痕迹。现今计算机图形学中使用最广泛的噪声是柏林噪声44。它是由KENPERLIN发明的，自从被引入计算机图形学以来，它就成为了计算机图形学中不可或缺的重要噪声。它由在不同频谱范围内的有限噪声带构成。这些噪声带被当作模块用来构建复杂的噪声类型，这可以通过在频谱中控制不同噪声带的权重来实现。例如，粉红噪声中它的每个频带与它的频率成反比。因为这个方法中每个频带使用的是单一的权重，所以它并不支持细节频谱的成形。相反，LEWIS45引入了WIENER插值法（WIENERINTERPOLATIONMETHOD）来构建一个包含随机频谱的纹理。这使控制噪声的特点变得更加容易，但是也显著地增加了很多开销。尽管研究者已经掌握构建和定义噪声的许多方法46，柏林噪声最原始的版本依然被广泛应用。因为，它效率很高，生成简单，频带覆盖了足够宽的频谱可控制的应用领域。每个柏林带宽是带宽受限的，所以它只包含在2的N次方周围的频带。但是，每个频带又包含了较宽的频率。由于柏林噪声的带宽限制不是非常严格，所以它会导致存在一些问题。柏林带宽在接近那奎斯特极限的地方包含了低频（包含了图像中的细节）和高频（造成锯齿）两种频率。如果把接近那奎斯特极限的这部分频率剔除会导致细节失真，但是包含这部分又会造成人工锯齿的存在。平衡这两者之间的关系使很多程序员进退两难。由于一般情况下锯齿比细节的丢失更加明显，很多时候大家都会选择减小频带的宽度来保证没有锯齿。这造成的后果就是，当拉近场景时，纹理的细节和几何细节并不同步。LEWIS曾经引入过稀疏卷积法来提高噪声的带宽受限性质，但是，这仍然不能根本解决细节丢失和产生锯齿之间的矛盾。更值得注意的一个问题是，当渲染场景的时候，即使三维的噪声满足带宽受限的要求，但是倘若在二维表面对三维噪声进行采样时，采样的结果在二维表面上并不满足带宽受限的性质。由于柏林噪声所存在的锯齿和细节丢失问题在湍流模拟中会