一种适用于三维流体模拟的三阶高精度对流插值算法的制作方法

本发明涉及计算机动画技术领域，是一种适用于三维流体模拟的三阶高精度对流插值算法，能有效降低对流数值耗散，减少内存开销，增强流体动画的真实感。

背景技术：

在计算机图形学和虚拟现实领域，基于物理的计算机动画技术始终是科研工作者的研究重点，而基于物理的流体动态模拟技术则是其中的一个热点问题。流体模拟技术在多个领域也被广泛应用，工程领域如航天、航空、航海等的发展也离不开流体技术的支撑。

在流体模拟中，数值耗散增加了流体的粘度，使其比预期的更粘稠，并涂抹细节以致视觉质量受损。导致数值耗散的因素中，对流的精度对流体模拟的视觉质量有很大的影响，相关学者为开发精确的对流求解器进行了各种尝试。其中，高阶插值格式广泛应用于计算流体力学(cfd)领域，包括本质非振荡(eno)、加权eno(weno)等，然而由于计算量太大不适用于图形学领域，而且这些插值算法在宽模板上进行，不适用于在非均匀网格上进行仿真。线性方法(linear)、前后向误差补偿矫正方法(bfecc)、maccormack法等对流方法广泛用于图形学中，但计算精度有限，最高只可达二阶精度。

约束插值剖面(cip)方法只在单个网格单元构造插值函数，具备紧模板特性，而且具有三阶精度，但扩展cip求解高维对流方程是一项困难的任务，需要大量的计算时间和内存消耗。另外，现有技术中基于cip的多维对流求解器，如单调cip(mcip)、非分裂型cip(uscip)等，只能部分地解决这些问题，有些甚至以损失数值精度或造成不稳定为代价。如何发展一种高效、高精度的高维cip方案对于提升对流方法的计算效率和精度颇有意义，但仍然是一个具有挑战性的问题。

技术实现要素：

针对现有技术的不足，本发明提供一种适用于二维流体模拟的三阶高精度对流插值算法：三维基于泰勒展开的约束插值剖面法(3dtaylorexpansionbasedconstrainedinterpolationprofile，3dtecip)。现有基于cip的方法应用到高维时，通常存在计算开销高、内存占用大、精度降低或者不稳定等问题。本发明方法只将物理量及其一阶导数作为未知变量，高阶导数并不作为变量保存，可以减少大量内存开销。其次，本发明使用局部泰勒展开按需计算高阶导数，计算速度较快且保持三阶精度，能够更有效保持丰富的流体细节。最后，方法并未破坏原cip方法的紧模板特性，除规则网格外同样适用于非均匀网格。

本发明的技术方案为：一种适用于三维流体模拟的三阶高精度对流插值算法，本发明主要用于提高求解流体运动控制方程中对流项的精度。通过利用半拉格朗日法求解流体运动控制方程的对流项时，针对回退点处物理量的计算采用了高精度插值，降低了数值耗散，同时为减少内存消耗，只存储物理场的值及其一阶导数作为计算变量。

对于三维对流方程，设φ为对流物理量，则在每个网格点i上存储4个变量，即该物理量及其一阶偏导：φi、和回退点上相应值通过高精度插值求得。具体包括以下步骤：

s1)针对回退点p所处的三维网格单元，设网格单元边长为h，八个顶点分别为a、b、c、d、e、f、g、h，网格的所有边均与三维空间的坐标轴平行，网格面abcd和面efgh均平行于xy平面且二者沿z轴正向的间距为h，网格面aehd和面bfgc平行于yz平面且且二者沿x轴正向的间距为h，网格面abfe和面dcgh平行于xz平面且二者沿y轴正向的间距为h。设点m和n分别为回退点p在面abcd和面efgh上的投影点；

s2)、计算八个顶点处φ的二阶偏导的值，具体计算如下：

基于以下三阶精度泰勒展开式，

令y＝z，有

对于网格点a，针对以上公式，令δx＝δz＝h，得到：

对于网格点b，针对以上公式，令δx＝-h，δz＝h，得到：

对于网格点c，针对以上公式，令δx＝-h，δz＝h，得到：

对于网格点d，针对以上公式，令δx＝δz＝h，得到：

对于网格点e，针对以上公式，令δx＝h，δz＝-h，得到：

对于网格点f，针对以上公式，令δx＝-h，δz＝-h，得到：

对于网格点g，针对以上公式，令δx＝-h，δz＝-h，得到：

对于网格点h，针对以上公式，令δx＝h，δz＝-h，得到：

s3)、计算八个顶点处φ的二阶偏导的值，具体计算如下：

基于以下三阶精度泰勒展开式，

令x＝y,y＝z，有

对于网格点a，针对以上公式，令δy＝δz＝h，得到：

对于网格点b，针对以上公式，令δy＝δz＝h，得到：

对于网格点c，针对以上公式，令δy＝-h，δz＝h，得到：

对于网格点d，针对以上公式，令δy＝-h，δz＝h，得到：

对于网格点e，针对以上公式，令δy＝h，δz＝-h，得到：

对于网格点f，针对以上公式，令δy＝h，δz＝-h，得到：

对于网格点g，针对以上公式，令δy＝-h，δz＝-h，得到：

对于网格点h，针对以上公式，令δy＝-h，δz＝-h，得到：

s4)、通过二维基于泰勒展开的约束插值剖面算法计算投影点m和n处φ值及其一阶、二阶偏导：φi、其中i∈{m,n}；

s5)、基于一维约束插值剖面法(cip)插值算法，通过步骤s4)求得的投影点m和n处的φm、φn、计算回退插值点p处的φp、通过步骤s4)求得的投影点m和n处的计算点p处通过步骤s4)求得的投影点m和n处的计算点p处

进一步的，步骤s4)中所述的二维基于泰勒展开的约束插值剖面算法，具体如下：

设φ为插值物理量，则在每个网格点i上存储3个变量，即该物理量及其一阶偏导：φi、插值点上相应值通过高精度插值求得。包括以下步骤：

s401)、针对插值点p′所处的二维网格单元，设网格单元边长为h′，四个顶点分别为a′、b′、c′、d′，网格的所有边均与二维维空间的坐标轴平行，边a′b′和边c′d′平行于x轴且二者沿y轴正向的间距为h′，边a′d′和边b′c′平行于y轴且二者沿x轴正向的间距为h′。设点e′和g′分别为插值目标点p′在边a′b′和a′d′上的投影点；

s402)、计算点a′、b′、d′处φ的二阶偏导的值。具体计算如下：

基于以下三阶精度泰勒展开式，

令δx＝δy＝h′，计算网格点a′处的二阶偏导值，即：

同理，令δx＝-h′，δy＝h′，计算网格点b′处的二阶偏导值，即：

同理，令δx＝h′，δy＝-h′，计算网格点d′处的二阶偏导值，即：

s403)、基于一维cip插值算法，通过点a′和点b′处的φa′、φb′、计算点e′的φe′、通过点a′和点b′处的计算点e′的

通过网格点a′和点d′的φa′、φd′、计算网格点g′的φg′、通过网格点a′和点d′的计算点g′的

s404)、根据以下公式计算得到p′点的φp′、

其中，ξ′和η′分别为a′e′和a′g′的长度。

进一步的，步骤s4)中，针对面abcd，由于点m为回退插值目标点p在平面abcd上的投影，设定m点在ab边上的投影为o，在ad边上的投影为k，ξ和η分别为ao和ak的长度。基于上述二维基于泰勒展开的约束插值剖面算法，以a、b、c、d四点处的φ、值作为输入，将点m作为插值目标点，可得到m处的对应值：φm、具体如下：

同理，以a、b、c、d四点处的值为输入，将点m作为插值目标点，可得到m处的对应值：计算公式如下：

进一步的，步骤s4)中，针对面efgh，由于点n为回退插值目标点p在平面efgh上的投影，设定n点在ef边上的投影为s，在eh边上的投影为j，则ξ和η分别为es和ej的长度。基于前述二维基于泰勒展开的约束插值剖面算法，以e、f、g、h四点处的φ、值作为输入，将点n作为插值目标点，可得到n处的对应值：φn、计算公式如下：

同理，以e、f、g、h四点处的值为输入，将点n作为插值目标点，可得到n处的对应值：计算公式如下：

进一步的，步骤s5)中，基于一维cip插值算法，通过投影点m和n处的φm、φn、计算回退插值点p处的φp、其计算式如下：

φp＝c3(zp-zm)³+c2(zp-zm)²+c1(zp-zm)+c0；

其中，c0＝φm，zp、zm分别为点p和点m的z轴坐标，h为网格边长。

进一步的，步骤s5)中，基于一维cip插值算法，通过投影点m和n处的计算回退插值点p处的其计算式如下：

其中，zp、zm分别为点p和点m的z轴坐标，h为网格边长。

进一步的，步骤s5)中，基于一维cip插值算法，通过投影点m和n处的计算回退插值点p处的其计算公式如下：

其中，zp、zm分别为点p和点m的z轴方向坐标，h为网格边长。

进一步的，步骤s403)中，基于一维cip插值算法，通过点a′和点b′处的φa′、φb′、得到点e′的φe′、其计算公式如下：

φe′＝c3(xe′-xa′)³+c2(xe′-xa′)²+c1(xe′-xa′)+c0；

其中，c0＝φa′，xe′表示e′点横坐标，xa′表示a′点横坐标，h′为网格单元边长。

进一步的，步骤s403)中，基于一维cip插值算法，通过点a′和点b′处的得到点e′的其计算公式如下：

其中，xe′表示e′点横坐标，xa′表示a′点横坐标，h′为网格单元边长。

进一步的，步骤s403)中，基于一维cip插值算法，通过网格点a′和点d′的φa′、φd′、得到网格点g′的φg′、其计算公式如下：

φg′＝c3(yg′-ya′)³+c2(yg′-ya′)²+c1(yg′-ya′)+c0；

其中，c0＝φa′，yg′表示g′点纵坐标，ya′表示a′点纵坐标，h′为网格单元边长。

进一步的，步骤s403)中，基于一维cip插值算法，通过网格点a′和点d′的得到点g′的其计算公式如下：

其中，yg′表示g′点纵坐标，ya′表示a点纵坐标，h′为网格单元边长。

本发明的有益效果为：

本发明能够在时间和内存消耗较少的前提下保持三阶高精度，且具有紧模板特性。对比现有方法，该方法在视觉质量、速度和内存消耗等方面都有明显的改进，能够有效提升流体模拟的对流精度和速度。

附图说明

图1为本发明方法用于三维流体模拟时的插值示意图；

图2为本发明方法中二维基于泰勒展开的约束插值剖面算法的插值示意图；

图3为本发明方法应用于三维流体模拟时的计算流程图；

图4为本发明与相关方法用于三维流体模拟的渲染结果对比，图中(a)、(b)、(c)、(d)、分别为采用linear算法、bfecc算法、mcip和uscip算法的模拟效果图，(e)为采用本发明方法(3dtecip)的模拟效果图。

具体实施方式

下面结合附图对本发明的具体实施方式作进一步说明：

提供一种适用于三维流体模拟的三阶高精度对流插值算法，用于计算流体控制方程的对流项。通过利用半拉格朗日法求解流体运动控制方程的对流项时，针对回退点处物理量的计算采用高精度插值，能够大大降低对流数值耗散，增强流体模拟的真实感。对于三维对流方程，设φ为对流物理量，则在每个网格点i上存储4个变量，即该物理量及其一阶偏导：φi、和回退点上相应值通过高精度插值求得。具体包括以下步骤：

s1)、针对泰勒展开式：

当δx＝0，

当δy＝0，

综合上述5式，并忽略三阶以上的余项，得到三阶精度的泰勒展开式：

整理后得到：

s2)、如图1所示，对于回退插值点p所处的三维网格单元，设网格单元边长为h，点m和n分别为点p在网格面abcd和efgh上的投影点。p点处的值φp、和由网格单元的八个网格点a、b、c、d、e、f、g、h的相应值通过发明的高精度插值求得；

s3)、基于步骤s1)得到的泰勒展开式，令y＝z，有

基于上式计算八个网格点处φ的二阶偏导的值，具体计算如下：

对于网格点a，针对以上公式，令δx＝δz＝h，得到：

对于网格点b，针对以上公式，令δx＝-h，δz＝h，得到：

对于网格点c，针对以上公式，令δx＝-h，δz＝h，得到：

对于网格点d，针对以上公式，令δx＝δz＝h，得到：

对于网格点e，针对以上公式，令δx＝h，δz＝-h，得到：

对于网格点f，针对以上公式，令δx＝-h，δz＝-h，得到：

对于网格点g，针对以上公式，令δx＝-h，δz＝-h，得到：

对于网格点h，针对以上公式，令δx＝h，δz＝-h，得到：

s4)、基于步骤s1)得到的泰勒展开式，令x＝y,y＝z，有

基于上式计算八个网格点处φ的二阶偏导的值，具体计算如下：对于网格点a，针对以上公式，令δy＝δz＝h，得到：

对于网格点b，针对以上公式，令δy＝δz＝h，得到：

对于网格点c，针对以上公式，令δy＝-h，δz＝h，得到：

对于网格点d，针对以上公式，令δy＝-h，δz＝h，得到：

对于网格点e，针对以上公式，令δy＝h，δz＝-h，得到：

对于网格点f，针对以上公式，令δy＝h，δz＝-h，得到：

对于网格点g，针对以上公式，令δy＝-h，δz＝-h，得到：

对于网格点h，针对以上公式，令δy＝-h，δz＝-h，得到：

s5)、通过二维基于泰勒展开的约束插值剖面算法计算投影点m和n处φ值及其一阶、二阶偏导：φi、其中i∈{m,n}；具体计算如下：

s501)、二维基于泰勒展开的约束插值剖面算法的具体实施步骤如下：

s50101)、设φ为插值物理量，则在每个网格点i上存储有3个变量，即该物理量及其一阶偏导：φi、如图2所示，针对插值点p′所处的二维网格单元a′b′c′d′，设网格单元边长为h′，点e′和g′分别为点p′在边a′b′和a′d′上的投影点，a′e′和a′g′的长度分别为ξ′和η′。p′点位置的值φp′、由网格单元的四个顶点a′、b′、c′、d′的相应值通过以下步骤插值求得；

s50102)、计算点a′、b′、d′处φ的二阶偏导的值。具体计算如下：

基于步骤s1)得到的泰勒展开式，

令δx＝δy＝h′，计算网格点a′处的二阶偏导值，即：

同理，令δx＝-h′，δy＝h′，计算网格点b′处的二阶偏导值，即：

同理，令δx＝h′，δy＝-h′，计算网格点d′处的二阶偏导值，即：

s50103)、基于一维cip插值算法，通过点a′和点b′处的φa′、φb′、计算点e′的φe′、设xe′、xa′分别表示e′点、a′点的横坐标，其计算公式如下：

φe′＝c3(xe′-xa′)³+c2(xe′-xa′)²+c1(xe′-xa′)+c0；

其中，c0＝φa′，

通过点a′和点b′处的计算点e′的其计算公式如下：

其中，

s50104)、基于一维cip插值算法，通过网格点a′和点d′的φa′、φd′、计算点g′的φg′、设yg′、ya′分别表示g′点、a′点的纵坐标，其计算公式如下：

φg′＝c3(yg′-ya′)³+c2(yg′-ya′)²+c1(yg′-ya′)+c0；

其中，c0＝φa′，

通过网格点a′和点d′的计算点g′的其计算公式如下：

其中，

s50105)、计算点p′处的φ，具体计算如下：

基于步骤s1)得到的泰勒展开式，

令δx＝ξ′，δy＝η′，可得点a′的值：

令δx＝-ξ′，δy＝η′，可得点e′的值：

令δx＝ξ′，δy＝-η′，可得点g′的值：

重新整理上面三式，得到以下公式以计算点p′处的φ，

s502)、如图1，针对网格面abcd，设m点在ab边上的投影为o，在ad边上的投影为k，ξ和η分别为ao和ak的长度。

基于s501)所述二维基于泰勒展开的约束插值剖面算法，以a、b、c、d四点处的φ、值作为算法输入，将点m作为插值目标点，可得到m处的对应值：φm、具体如下：

基于s501)所述二维基于泰勒展开的约束插值剖面算法，以a、b、c、d四点处的值作为算法输入，将点m作为插值目标点，可得到m处的对应值：具体如下：

s503)、如图1针对网格面efgh，n点在ef边上的投影为s，在eh边上的投影为j，则es和ej的长度分别为ξ和η。

基于s501)所述二维基于泰勒展开的约束插值剖面算法，以e、f、g、h四点处的φ、值作为输入，将点n作为插值目标点，可得到n处的对应值：φn、具体如下：

基于s501)所述二维基于泰勒展开的约束插值剖面算法，以e、f、g、h四点处的值为算法输入，将点n作为插值目标点，可得到n处的对应值：具体如下：

s6)、基于一维cip插值算法，通过投影点m和n处的φm、φn、计算回退插值点p处的φp、其计算式如下：

φp＝c3(zp-zm)³+c2(zp-zm)²+c1(zp-zm)+c0；

其中，c0＝φm，zp、zm分别为点p和点m的z轴坐标，h为网格边长。

s7)、基于一维cip插值算法，通过投影点m和n处的计算回退插值点p处的其计算式如下：

其中，zp、zm分别为点p和点m的z轴坐标，h为网格边长。

s8)、基于一维cip插值算法，通过投影点m和n处的计算回退插值点p处的其计算公式如下：

其中，zp、zm分别为点p和点m的z轴方向坐标，h为网格边长。

图3为三维流体模拟的流程图，其中本发明提供的方法主要用于计算流体控制方程的对流项。图4为二维流体模拟结果渲染图对比，图4(e)为本发明方法用于三维流体模拟的结果渲染图，图4(a)、图4(b)、图4(c)、图4(d)分别为采用低阶精度的linear和bfecc算法、以及mcip和uscip算法的结果渲染图，从图中可以看出，基于cip的三种方法模拟效果优于低阶精度的linear和bfecc算法。由于前述的mcip和uscip方法的缺陷，本发明方法在三维流体模拟中取得的效果优于mcip和uscip，模拟的烟雾细节更丰富。

上述实施例和说明书中描述的只是说明本发明的原理和最佳实施例，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。