一种基于物理特性的雨滴飞溅模拟算法

，，

(浙江工业大学 计算机科学与技术学院，浙江 杭州 310023)

近年来,有关自然景物的模拟[1-2]已成为计算机图形学领域的一个研究热点.降雨过程及雨滴与固体表面的交互也日益受到该领域研究者的重视[3-4].下落的雨滴跟固体表面接触后,因冲击力而向四周迸射的物理现象叫做飞溅.对雨滴飞溅过程的模拟，有着非常重要的意义.主要表现在：1) 逼真地模拟雨滴飞溅过程，是计算机图形学领域中重要的研究课题，它在数字娱乐、虚拟仿真、影视及广告中具有重要的应用价值和广阔的市场潜力；2) 模拟雨滴在植物冠层内滴落、飞溅与沉积的动态过程，对研究病害在植物冠层中的传播途径，进而确定最佳的防治策略，具有重要的指导意义.

目前，已有不少学者致力于雨滴飞溅现象的研究.100多年前，Worthington第一次通过频闪观测仪拍摄飞溅水滴的照片对飞溅现象进行分析[5-6].Stow C等对飞溅水滴的数量、大小和速度的关系函数进行了数值研究[7].S.Saint-Jean等对雨滴与植物碰撞后飞溅的传输路径进行二维的数值分析[8].但以上的研究并未涉及雨滴的三维可视化模拟.K.Garg等建立一种经验概率模型来模拟基于倾斜角和表面材质的水滴飞溅现象[9]，并采用Kshitiz[10]中的算法来绘制雨景，但未考虑空气阻力的影响.汪继文等[11]实现了单个水滴的溅落过程，但没对雨的溅落效果进行绘制.Thornton等编写了一款可在Maya平台上使用的插件程序，通过手工绘制，允许用户构建卡通风格的液体泼溅场景[12].由于他们主要是针对非真实感动画的绘制和卡通化模拟，所以并没有考虑雨滴行为的物理过程及雨滴的受力情况.针对目前雨水飞溅研究的现状，提出一种基于物理特性的雨滴飞溅跟踪算法.该算法考虑到空气阻力及风力对雨滴的影响，同时基于粒子系统，利用GPU加速雨场景的绘制.这种方法不仅能完善雨场绘制[13]的模拟方法，也对雨滴飞溅路径的研究[14-16]及泼溅效果的模拟提供方法上的借鉴.

1 雨滴的模拟及雨滴受力分析

1.1 雨滴的模拟

粒子系统被认为是模拟不规则模糊物体最为成功的一种图形生成算法[17-19].它由大量不规则的微小粒子图元组成，而每个粒子均具有生命周期和属性值.在生命周期内，经历产生、活动、消亡三个阶段，并按照一定的规则改变形状，不断运动[11].

该算法采用粒子系统对雨场进行建模.在粒子系统中构建2种不同的粒子：落地前的雨滴粒子和落地后的飞溅粒子.其中一个雨滴粒子包含N个飞溅粒子.雨滴粒子在空间的一定高度产生，由随机函数控制，呈均匀分布，粒子数目由全局变量NRain控制.飞溅粒子在雨滴跟地面碰撞的瞬间产生，继承其所属雨滴的当前落地的位置，粒子数目由全局变量N控制.在无风的情况下，雨滴的运动方向默认为竖直向下，而飞溅粒子的方向由其高度角和方位角两个参数决定.两种粒子的数据结构可分别描述为：

Typedef struct{

float pos[3];//雨滴粒子的空间坐标

float vel[3];//雨滴粒子的空间速度

float gra[3];//雨滴粒子的空间加速度

float alpha;//雨滴粒子的透明值

float radius;//雨滴粒子的半径大小

int type;//雨滴粒子的类型

int life;//雨滴粒子的生命周期

PARTICLES pa[Num];//飞溅粒子

}Rain;//雨滴粒子

Typedef struct{

float p[3];//飞溅粒子的空间坐标

float v[3];//飞溅粒子的空间速度

float g[3];//飞溅粒子的空间加速度

float alpha;//飞溅粒子的透明值

float r;//飞溅粒子的半径大小

float elevationAngle;//粒子的高度角

float azimuthalAngle;//粒子的方位角

int life;//飞溅粒子的生命周期

} PARTICLES；//飞溅粒子

为了提高运行效率并节省计算时间，做出如下假设：1) 雨滴在下落过程中只考虑空气阻力、水平风力以及重力的影响；2) 粒子在飞溅过程中只考虑重力对粒子的影响；3) 雨滴粒子和飞溅粒子在运动过程中无能量损失，无水分的蒸发，不考虑形变.

1.2 雨滴受力分析

为简化雨滴在空中的下落过程，假设其初始时刻的速度为零.根据受力情况，分析得到雨滴开始做加速运动，当其受到的空气阻力与重力平衡时，转为匀速下落.图1为雨滴下落过程的受力简化图.

图1 作用在雨滴上的外力

图1中，Fw为雨滴在水平方向上所受的风力，Fa为雨滴在水平方向上所受的空气阻力，Fv为雨滴在竖直方向上遇到的空气阻力，G为雨滴的重力.假设雨滴的形状为圆球形，半径为r，密度为ρ，质量为m，重力系数为g，雨滴速度为v，空气密度ρa=1.21 kg/m3.根据牛顿动力学公式，可以计算出雨滴在下落过程中任意时刻的位移和速度，即

(1)

G=mg

(2)

(3)

(4)

其中：Fd为空气阻力是Fv和Fa的和值；CD为空气阻力系数跟雷诺系数Re成一定比例，计算公式[20]为

(5)

其中：Re=ρavd/μa；d为飞溅粒子的直径；μa为运动粘滞系数.

在雨场模拟中，水平风力的大小由随机函数控制.由于考虑到水平风力和空气阻力的影响，该算法采用二阶龙格库塔法来求解雨滴行为的微分方程，计算公式为

(6)

其中参数h为时间间隔.

二阶龙格库塔法是显示单步式，每前进一步需要计算两个函数值k1和k2.故适当的选取参数h，可使每步计算得到的结果达到所需的精度.这里时间间隔h选取为0.02 s.

2 雨滴飞溅算法

2.1 飞溅粒子属性的初始化和更新

雨滴的飞溅过程与雨滴的物理属性和接触物的表面属性有关.研究表明[21]：雨滴的下落速度越大，飞溅粒子的数目越多；雨滴的半径越大，飞溅粒子的平均半径也越大.而物体表面的粗糙度、硬度、湿润度、倾斜度和亲水性等物理属性也会影响雨滴的飞溅过程.

假设半径为r的雨滴从y轴方向自由落体跟表面S发生碰撞后，破碎过程中产生的飞溅粒子的数量为n.然后n个飞溅粒子将以各自的初速度v朝着不同的方向运动.运动方向由高度角θ和方位角φ组成.图2所示为雨滴粒子跟表面S碰撞后，其中一个飞溅粒子的传播路径.表面S的法向量为N，表面倾斜角为α.

图2 粒子飞溅的示意图

Gary等通过实验得到飞溅现象的数据，发现飞溅粒子的属性值很大程度上符合高斯分布[9].因此，为了描述飞溅粒子的属性，该算法采用高斯分布函数对飞溅粒子进行初始化[22].

1) 飞溅粒子的数目

雨滴碰撞破碎后，产生的飞溅粒子的数目n符合高斯分布函数G(x;u,δ)，u,δ为已知变量，则

n=n0+nvar·rand()

(7)

其中：n0为平均值；nvar为方差.

2) 飞溅粒子的方向

粒子的高度角θ符合高斯分布函数G(x;u,δ)，u,δ为已知变量.方位角φ每隔x°统一分布，其中x=360/n.

θ=θ0+θvar·rand()

(8)

其中：θ0为平均值；θvar为方差.

3) 飞溅粒子的初速度

粒子的初速度大小v符合高斯分布函数G(x;u,δ)，其中u,δ为已知变量，则

v=v0+vvar·rand()

(9)

其中：v0为平均值；vvar为方差.

根据式(9)求得的v值，可分别得到粒子在x,y,z轴上的分速度，即

vx=v·cosθ·cosφ

(10)

vy=v·sinθ

(11)

vz=v·cosθ·sinφ

(12)

为了真实的反映飞溅的物理过程，参数n0，v0，θ0，nvar，vvar和θvar的取值参照文献[9]中的实验数据.从其他学者的物理实验结果[8]可知：雨滴的最小直径为0.5 mm，最大直径不超过6 mm，而直径在0.5 mm和5 mm之间的雨滴能产生飞溅现象，飞溅粒子的半径r不大于1.5 mm，飞溅粒子的速度不大于20 m/s.由此可以限定飞溅粒子属性的取值范围.

4) 飞溅粒子的速度更新

飞溅过程中，粒子保持x和z方向上的速度不变.由于雨滴在飞溅过程中经历的时间较短，所以风力和空气阻力对它的影响很小，可以忽略不计.为增强粒子溅落的效果，使用变量k来控制垂直方向上的速度.故y轴方向上的速度表示为

particle[i]·v[2]=particle[i]·v[2]+

particle[i]·g[2]·t·k

(13)

其中：0≤i≤n-1；k为速度调节因子；t为时间变量.

5) 飞溅粒子的位置更新

飞溅过程中只考虑重力对粒子的影响，根据牛顿第二定律，求得飞溅粒子的运动轨迹为

particle[i]·p[0]=particle[i]·p[0]+

particle[i]·v[0]·t

(14)

particle[i]·p[1]=particle[i]·p[1]+

particle[i]·g[1]·t·t

(15)

particle[i]·p[2]=particle[i]·p[2]+

particle[i]·v[2]·t

(16)

式中：0≤i≤n-1;t为时间变量.

2.2 雨滴路径跟踪算法

跟踪雨滴粒子飞溅路径的算法步骤如下：

Step1初始化雨场的雨量、半径、空间位置、空间加速度、空间速度和生命周期等属性值.

Step2根据牛顿运动定律，计算得到雨滴下落过程的位置和速度.

Step3检测雨滴是否跟表面发生碰撞.如发生碰撞，则将雨滴粒子的生命周期设为0，转至Step 4；若没发生碰撞，则转至Step 2.

Step4雨滴粒子碰撞后破碎成飞溅粒子，修改雨滴的类型type，初始化飞溅粒子各个属性值(坐标、加速度、速度、高度角和方位角、生命周期等).

Step5根据飞溅粒子的初速度和受力情况，更新计算得到飞溅粒子的运动轨迹.

Step6判断飞溅粒子的生命周期.若为0，粒子消亡.

Step7显示所有现存的粒子，即绘制粒子.

3 基于GPU加速技术的场景绘制

3.1 粒子属性的存储方式

在模拟雨场飞溅效果时，该算法将粒子的所有属性都存储到相应的纹理像素中.将速度和位置属性的分量x，y，z分别存储到相应浮点纹理的R，G，B三个分量中[23].为避免GPU对显示存储器访问的读写冲突，该算法采用双缓冲区的办法.即创建一对纹理贴图，其中一块用于输入，以读取信息供cg代码计算，另一块则作为输出，以保存计算结果.粒子的速度和位置信息都存储在一对纹理贴图中.由于粒子的其他属性(大小、颜色、生命周期、透明度等)不需很复杂的计算.因此，可将它们存储在单一纹理内即可.粒子属性的存储方式如图3所示.

图3 粒子属性的存储方式

3.2 粒子的绘制过程

粒子绘制算法分为CPU运行阶段和GPU运行阶段.粒子的产生和消亡需要借助CPU完成.首先，在内存中建立一个空闲堆栈当作索引存储器.然后当新粒子产生时，从索引分配器中取出一个空闲的索引号，并将粒子的属性值写入到索引号指向的粒子堆栈中，这就完成了粒子的产生.当粒子的生命周期结束时，将粒子的索引号压入空闲堆栈中，即完成粒子的消亡.而粒子属性的更新和粒子的绘制需要借助GPU完成.具体实现步骤如下：

Step1以纹理的形式，将粒子属性读入GPU，记为纹理X.

Step2根据前一帧的速度和位置值，运用相应的运动学公式计算出下一帧的速度和位置.

Step3将Step 2中最后得到的结果保存到另一个纹理中，记为纹理Y.

Step4交换纹理X和纹理Y，即将保存结果的纹理Y作为下次计算的输入，而纹理X则作为输出，保存下次计算的结果.

以上粒子属性的更新是以顶点的表示方式在Vertex Shader中完成.最后通过Fragment Shader来完成绘制.

4 算法实现与结果

依据上述建模思想，在Visual Studio 2005的集成开发环境下，该算法使用C++编程语言并采用OpenGL三维图形渲染机制，对雨滴飞溅现象进行可视化模拟.采用的计算机硬件配置：CPU为Inter i5-750；主频为2.67 G；内存为4 GB；显卡为NVIDIA GeForce3.

图4 单个粒子随时间变化的泼溅效果

图5 不同雨量的飞溅效果

图4为单个粒子随时间变化的泼溅效果.图5为绘制不同雨量时的渲染效果.从图5中可以看出：该算法逼真的模拟了不同雨量时下雨的场景和雨滴落地的飞溅效果，具有较好的视觉和实时效果.在雨量和场景几何复杂度不同的情况下，系统运行的性能也有所改变.在图5中，场景的粒子系统共有20 000雨滴粒子，整个模拟过程绘制的速度在10～25帧/s的范围内，能满足实时性的需求.

与以往的工作相比，在绘制效果上，该方法更加逼真地展现了雨滴飞溅水花的大小和形状.在算法上，该方法考虑空气阻力和水平风力对雨滴的影响，构建了两种不同的粒子来描述其运动的物理过程，并对雨滴落地的飞溅路径进行跟踪模拟.在绘制效率上，该算法利用GPU可编程管线的数据并行处理方式以及OPENGL的内存对象，采用双缓存的数据存储结构.跟基于CPU的方法相比，大大提高了绘制效率.

5 结束语

通过粒子系统构建雨场，结合统计学知识和前人的实验数据模拟雨滴的行为轨迹，提出一种基于物理特性的雨滴飞溅算法.实验结果表明：该算法得到了良好的视觉效果，利用GPU的并行性和可编程性，提高了实验的运行速度.在未来的工作中，将进一步完善目前的模型，考虑雨水下落和飞溅过程中雨滴自身形状的变化.此外，表面湿润度、粗糙度、亲水性等属性对飞溅现象的影响也将应用到雨场的绘制中.同时加入更多的随机函数处理方法来提高效果的渲染，这些都有待后续进一步的工作和完善.

参考文献：

[1] 丁维龙，陈敏智，程志君．植物器官可视化建模系统的设计与实现[J]．浙江工业大学学报，2009，37(3)：295-310．

[2] 马培良，丁维龙，古辉．基于OpenGL和双三次贝塞尔曲面的稻叶可视化建模[J]．浙江工业大学学报，2010，38(1)：36-40．

[3] TATARCHUK N, ISIDORO J. Artist-directable real-time rain rendering in city environments[C]//In Proceedings of SIGGRAPH 2006．New York：ACM Press，2006：23-64．

[4] HUBER L, FITT B D L, MCCARTNEY H A. The incorporation of pathogen spores into rain-splash droplets[J]. Plant Pathol，1996，45(3)：506-517．

[5] WORTHINGTON A M. On the forms assumed by drops falling vertically on a horizontal plate[J]．Proc of Royal Society of London，1876，25：261-271．

[6] WORTHINGTON A M. A second paper on the forms assumed by drops of liquids falling vertically on a horizontal plate[J]. Proc of Royal Society of London，1877，25：498-503．

[7] STOW C D, STAINER R D. The physical products of a splashing water drop[J]. Meteorological Society of Japan Journal，1977，55(5)：518-532．

[8] SAINT-JEAN S, CHELLE M, HUBER L. Modelling water transfer by rain-splash in a 3D canopy using Monte Carlo integration[J]．Agricultural and Forest Meteorology，2004，121(1/2)：183-196．

[9] GRAG K, KRISHNAN G, NAYAR S K. Material based splashing of water drops[C]//In Proceedings of Eurographics Symposium on Rendering. Switzerland: Eurographics Association，2007：171-182．

[10] KSHITIZ G, SHREE K N. Photorealistic rendering of rain streaks[J]. ACM Trans on Graphics，2006，25：996-1002．

[11] 汪继文，陆和军．基于粒子系统的水滴溅落模拟[J]．科学技术与工程，2010，10(6)：1547-1550．

[12] THORNTON J. Directable simulation of stylized water splash effects in 3D space[C]//In Proceedings of SIGGRAPH 2006-Material Presented at the ACM SIGGRAPH 2006 Conference. Boston: ACM Press，2006：94-94．

[13] PIERRE R, VINCENT J. Realistic real-time rain rendering[J]. Computer & Graphics，2006，30(4)：507-518．

[14] YANG Meng, HUANG Meng-cheng, YANG Gang, et al. Physically-Based animation for realistic interactions between tree branches and raindrops[C]//In Proc of the 17thACM Symposium on Virtual Reality Software and Technology (ACM VRST 2010)．New York：ACM Press，2010：83-86．

[15] YANG Meng, JIANG Long-sheng, LI Xiao- sheng, et al. Interactive coupling between a tree and raindrops[J]. Computer Animation and Virtual Worlds，2012，23(3/4)：267-277．

[16] WANG H, MUCHA J P, TURK G. Water drops on surfaces[J]. ACM Trans on Graphics，2005，24(3)：921-929．

[17] XU Y Z, LEE K, KIM E, et al. Simulation of smoke to improve unity 3D games engine particle system based on FDS[J]． Soft Technology，2013，19：183-186．

[18] KIM M, CHEN G N. Dynamic particle system for mesh extraction on the GPU[C]//In Proceedings of the 5th Annual Workshop． New York：ACM Press，2012：38-46．

[19] DING Wei, ZHU Zhong-bin. Real-time rain and snow rendering[C]//In 2013 Second International Conference on Agro-Geoinformatics．Fairfax：Agro-Geoinformatics Press，2013：32-35．

[20] MERCER G, SWEATMAN W L, ELVIN A, et al. Process driven models for spray retention by plants[C]//In Proceedings of the 2006 Mathematics in Industry Study Group. New Zealand: Mathematics in Industry Study Group Press，2006：57-85．

[21] STOW C, HADFIELD M. An experimental investigation of fluid flow resulting from the impact of a water drop with an unyielding dry surface[J]. Proc of Royal Society of London 1981,373：419-441．

[22] MUNDO C, SOMMERFELD M, TROPEA C. Droplet-wall collisions: experimental studies of the deformation and breakup process[J]. International Journal Multiphase Flow，1995，21(2)：151-173．

[23] 许楠，郝爱民，王莉莉．一种基于GPU的粒子系统[J]．计算机工程与应用，2006，19(3)：77-79．