本发明涉及计算机图形学，尤其涉及一种基于叠层结构的三维球坐标系下不可压缩流体仿真方法。
背景技术：
1、在有限厚度的球壳内的流动运动能呈现出十分有趣的视觉效果。例如，地球上的流动以及与复杂地形的相互作用，构成了各种有趣而壮观的气象现象，这些气象现象经常出现在科幻电影中，并引起人们对它们的模拟的关注。在气象科学界，人们对地球上的大气层在某些确定的初始条件下如何演化的问题进行了深入研究。气象科学领域存在一些广泛应用的通用大气环流模型(或气候模型)，如echam、ncar-ccsm、gfs模型等。这些模型除了复杂性之外，计算成本也非常高，这对计算机图形学界来说是一个巨大的挑战。另一方面，在计算机图形学中，人们可能希望在动画中流体域的厚度在某种程度上与地面半径相当，这将比地球大气层厚度与其半径的比率大得多。对于此类模拟，需要一个用于三维域中球坐标流动运动的流体动力学系统的实用物理模型。与基于笛卡尔坐标系的模拟相比，使用球坐标系的模拟更适合径向切向流运动的精确恢复和控制。
2、在计算机图形学领域中，现有的在球坐标系下流体仿真的相关工作主要集中在二维的球面上流体仿真，即忽略球坐标系的径向方向，只在一个二维球表面上进行流体仿真，没有径向的运动。而仅有的三维球坐标系仿真工作主要使用谱方法将流体用大量正交基函数进行投影，然后在该基函数组成的空间中进行求解，过程十分复杂，对于美术人员来说难以实现流体运动控制，且不能很好地与固体几何进行耦合。因此，计算机图形学领域需要一种使用三维球坐标系在三维球壳空间中仿真流体的新方法。
技术实现思路
1、本发明旨在解决使用三维球坐标系在球壳空间中进行不可压缩流体仿真的问题，因此提出了一种基于叠层结构的三维球坐标系流体仿真方法。
2、本发明通过如下技术方案实现
3、一种基于叠层结构的三维球坐标系下不可压缩流体仿真方法，其特征在于，使用叠层结构将三维球壳空间离散为多个二维球面层，作为一种欧拉视角下的流体仿真算法，根据所设置场景对各层包括速度、密度在内的物理场进行初始化，包括如下步骤：
4、步骤1、进行对流计算，更新旧的速度场为新的中间速度场；
5、步骤2、根据前向欧拉积分方法，更新步骤1计算得到的中间速度场为新的速度场；
6、步骤3、根据截断的球谐函数展开数值求解由流体不可压缩条件推导的泊松方程，计算得到压力场，根据该压力场更新由步骤2计算得到的速度场的径向分量；
7、步骤4、根据重平衡步骤，调整步骤3计算得到的速度场；
8、步骤5、根据快速傅里叶展开数值求解各离散层由三维不可压缩条件推导的二维泊松方程，计算得到压力场，根据该压力长更新由步骤4计算得到的速度场的经纬方向分量；
9、回到循环的步骤1，进行下一迭代步的计算过程，直到完成模拟。
10、根据本发明提供的模拟循环内流程，所述对流计算过程为将各速度分量分开求解，其三个分量的对流方程分别为：
11、
12、
13、
14、其中r，θ，分别为三维球坐标系下的径向方向和纬度经度方向的坐标，u为速度场，角标为该物理量在对应方向上的分量，即ur、uθ和分别为速度u在方向、方向和方向上的分量。求解时使用半拉格朗日方法，先线性插值得到该位置的经纬方向速度分量，根据该速度沿着该层的大圆方向前向回溯一个时间步到达中间位置点，再在该中间位置点线性插值得到径向速度ur，根据ur沿着径向方向向前回溯一个时间步到达新的位置点，再该位置点上线性插值得到新的速度，再将该插值得到的新速度投影到起始点的r，θ，方向上作为更新的中间速度场三个分量。
15、所述的使用前向欧拉法更新由外力导致的速度变化为：
16、u**＝u*+δtf,
17、其中u**为更新后的速度场，u*为上一步骤对流求解后得到的中间速度场，δt为仿真时间步长，f为外力。
18、所述的根据流体不可压缩条件推导得到的泊松方程为：
19、
20、其中p为压力场，m为对应离散单元中流体体积分数，ρ为密度，为球坐标系下的nabla算符。使用截断球谐展开对上述方程进行求解，先将等式两边进行球谐展开，球谐展开公式为：
21、
22、其中为在球坐标系下的标量物理场，为连带勒让德多项式，n为阶数，m为度数，n为指定的最高的连带勒让德多项式阶数，和分别为的实系数和虚系数。n越大，球谐展开所需时间越长，泊松方程求解的精度越高。将泊松方程等式左边的p和等式右边的q进行球谐展开后，通过球谐函数的特性将原式简化为一系列关于多项式系数的三对角线性方程组：
23、
24、其中和为对应连带勒让德多项式在第i层的系数，ri＝rmin+(i+0.5)δr为球坐标径向方向的坐标。求解后得到的各系数通过代入球谐展开公式可得到求解后的压力场，使用该压力场对步骤2更新后的速度场的径向方向分量进行更新：
25、
26、所述的基于叠层结构的三维球坐标系下不可压缩流体仿真方法，其特征在于，所述的重平衡策略为通过强制执行层间零通量边界条件来确保步骤3更新后的速度径向分量在流体层内的不可压缩性。该策略首先计算各层在径向方向上的通量为fk，然后计算该层球面总面积sk，其中k为第k层，然后对各速度径向方向分量进行更新：
27、
28、所述的基于叠层结构的三维球坐标系下不可压缩流体仿真方法，其特征在于，所述的根据快速傅里叶展开数值求解各离散层由三维不可压缩条件推导的二维泊松方程，该方程为：
29、
30、其中代表在球面方向上进行的nabla算符，u为经过步骤3后得到的速度场，为经过步骤4后得到的速度径向分量，分别为对应速度分量的流体体积分数，为周边各两个的平均值。
31、使用快速傅里叶变换将等式两边从时域变换到频域上，将原二维泊松方程分解为各不同频率下的三对角线性方程组，与步骤3类似，求解后通过傅里叶逆变换得到压力场解，根据该压力场对速度的经纬方向分量进行更新：
32、
33、
34、本发明的优点和有益效果：
35、本发明提供了一种计算机图形学中在三维球壳空间内不可压缩无粘流体流动的求解器，可对任意固体地形进行单向的流固耦合。本发明提供一种叠层结构算法，能够高效执行三维球形不可压缩无粘流体流动模拟，同时由于使用三维-二维结合的投影方法，使得仿真器能够以不高的计算开销达到不错的视觉效果。本发明提供一种灵活的流体控制方法，能够方便美术人员进行流体效果的调试。
1.一种基于叠层结构的三维球坐标系下不可压缩流体仿真方法，其特征在于，使用叠层结构将三维球壳空间离散为多个二维球面层，作为一种欧拉视角下的流体仿真算法，根据所设置场景对各层包括速度、密度在内的物理场进行初始化，包括如下步骤：
2.根据权利要求1所述的基于叠层结构的三维球坐标系下不可压缩流体仿真方法，其特征在于，步骤1中在求解对流方程时，对各分量分开求解，其三个分量的对流方程分别为：
3.根据权利要求1所述的基于叠层结构的三维球坐标系下不可压缩流体仿真方法，其特征在于，步骤2中使用前向欧拉法更新由外力导致的速度变化为：
4.根据权利要求1所述的基于叠层结构的三维球坐标系下不可压缩流体仿真方法，其特征在于，步骤3中根据流体不可压缩条件推导得到的泊松方程为：
5.根据权利要求1所述的基于叠层结构的三维球坐标系下不可压缩流体仿真方法，其特征在于，步骤4中的重平衡策略为通过强制执行层间零通量边界条件来确保步骤3更新后的速度径向分量在流体层内的不可压缩性，该策略首先计算各层在径向方向上的通量为fk，然后计算该层球面总面积sk，其中k为第k层，然后对各速度径向方向分量进行更新：
6.根据权利要求1所述的基于叠层结构的三维球坐标系下不可压缩流体仿真方法，其特征在于，步骤5所述的根据快速傅里叶展开数值求解各离散层由三维不可压缩条件推导的二维泊松方程，该方程为：
