CN107025332A - 一种基于sph的织物表面微观水扩散过程可视化方法 - Google Patents

Abstract

本发明公开一种基于SPH的织物表面微观水扩散过程可视化方法，该方法包括如下步骤：构建SPH微观水流体；查找SPH微观水流体的空间网格；获取基于SPH微观水流体的表面张力和接触角；设计织物结构；设计基于织物结构的扩散通道；设计微观水流体与织物的碰撞与吸收。本发明从小尺度水体出发，主要针对三维空间内不同小尺度流体和不同材质的织物发生接触时的反应，和流体分子在织物内部的扩散过程做仿真，并将结果可视化显示出来。

Description

一种基于SPH的织物表面微观水扩散过程可视化方法

技术领域

本发明涉及已知流体与织物材质的具体物理特性下的仿真过程计算和可视化显示技术领域，具体涉及一种基于SPH(smoothed particle hydrodynamics,光滑粒子动力学)的织物表面微观水扩散过程可视化方法。

背景技术

目前的流体仿真技术主要有两种，一种是无网格流体仿真，也被称为基于粒子的流体仿真，另外一种是基于网格的流体仿真。

基于粒子的流体仿真技术，使用的多是拉格朗日方法，考虑的流体内微团，或者说质点随时间的位置变化规律，用质点的运动规律来表示整个流体的运动情况。

光滑粒子法(“Smoothed Particle Hydrodynamics Method”)，就是近20年来发展起来的一种仿真计算方法。最初由Lucy等人在1977年提出，并在后来不断完善和改进，也是目前流体仿真常用的力学方法。

如前文所述，光滑粒子法，也称光滑粒子流体动力学，用粒子代表流体内的物理质点，承载流体内具体位置的密度，粘度，速度，加速度等物理属性，通过粒子的各自运动来得到整个粒子系统的运动情况。在粒子数量足够多的情况下，能够更准确地描述流体中各个位置在各个时刻的属性变化，可以更准确地描述整个力学过程。而且由于粒子间不依赖网格，可以避免流体变形时发生的网格扭曲问题。

但是拉格朗日方法也有其缺陷，虽然说粒子数足够多的情况下，可以完全表示流体的整体运动，但事实上很难做到对成千上万的流体质点进行实时的跟踪计算，计算成本太高。

基于网格的流体仿真主要使用的是欧拉方法，不考虑具体流体内质点，以流场空间为对象，研究室的各个时刻流场内各空间点上的属性变化，记录的是流体质点流经这些空间点的运动，把带属性的空间点综合考虑，从而得到流体的整体运动情况。在某一时刻某个空间节点的属性即经由该空间节点的流体质点的属性，它们在属性上是一致的。

与拉格朗日方法记录粒子的位置变化不同，欧拉记录的是经由空间点的流体粒子的速度向量。欧拉方法因为这种固定的网格节点观测的便利性，不仅仅应用于流体仿真计算，同时还广泛应用在气象观测中。

欧拉方法中，不同于牛顿第二定律的加速度表示，而是使用随体导数描述流体质点的 速度矢量随时间的变化情况。原因有两方面，一是流体内场的时间的非定常性会引起速度变化，二是流体内部的空间是不均匀的，当流体质点随时间变化而运动时，质点周围的空间也在发生变化，这种不均匀的变化会引起流体质点的速度变化。

根据时间的非定常性，通过求导可以得到局部导数，描述时间的非定常性引起的速度变化。而空间的不均匀性的影响同样可以用导数表示，称为对流导数。局部导数和对流导数之和即欧拉方法中的随体导数，表示流体质点的加速度，描述流体质点的速度变化过程。

欧拉方法更关注整个流体场的变化，因此更多的用于较大规模的流体仿真计算，在划分出网格后可以根据流体的边界条件开始从边界网格开始计算流体内部的流动状态。当所有网格点的数据变化收敛时，所有网格点一起反映的就是流体的运动状态。目前比较有名的工程分析软件CFX，Fluent用的就是欧拉法。

拉格朗日法更注重细节的表达，通过更多的粒子的追踪来得到更具体的流体信息。例如Realflow用的就是拉格朗日方法中的光滑粒子流体动力学方法。

在小尺度水体仿真方面，由于水体规模较小，更要求水体的细致表达，而且，在水体与织物的碰撞接触时会有比较大的形态变化，所以选用基于粒子的拉格朗日方法会更合适一些。

对于拉格朗日方法中具有代表性的sph光滑粒子流体动力学方法，在近20年来也被广泛使用和改进。有改进光滑粒子流体动力学中光滑核函数的，针对传统光滑核函数面对巨大应力变化时不稳定提出改进，也有水面近场爆炸方面的应用的，还有模拟绘画颜料逐渐干涸凝固的等等。

这些成果，多是应用在较大尺度的、形状变化较剧烈的水体的上的，很少有应用到小尺度的、从较微观角度上应用和模拟计算水体变化的，也很少有三维空间内不同小尺度水体和不同织物间发生吸收和扩散的可视化过程实现的。

随着在计算机虚拟现实中对可视化仿真的需求，越来越多的的流体仿真计算方法被提出和改进，其中，比较有名的软件就有CFX，Fluent，RealFlow等。

但是，一方面，部分此类软件，或者方法更关注水体本身的可视化表达，而没有更多的考虑流体与外界接触物的一些反应和变化，尤其是被吸收后的后续扩散过程。

另一方面，就是部分软件或方法考虑的是较大尺度水体的反应，很少考虑到小尺度水体的仿真计算。例如RealFlow,主要就是用于水面波浪的模拟，包括物体落入水面后的浮动和水花、水面高速运动物体的尾迹仿真等等。

发明内容

本发明的目的在于克服现有技术存在的上述不足，提供一种基于SPH的织物表面微观 水扩散过程可视化方法，具体技术方案如下。

一种基于SPH的织物表面微观水扩散过程可视化方法，其包括如下步骤：

(7)构建SPH微观水流体；

(8)查找SPH微观水流体的空间网格；

(9)获取基于SPH微观水流体的表面张力和接触角；

(10)设计织物结构；

(11)设计基于织物结构的扩散通道；

(12)设计微观水流体与织物的碰撞与吸收。

进一步地，步骤(1)具体包括：

(a)流体质点区域密度计算

流体质点粒子的区域密度通过搜索质点粒子附近的粒子数目累计计算得出；流体中的每一个质点粒子都有其作用力的影响距离极限标记为R，即影响力的作用半径，也是SPH中的光滑核半径；通过对以R为半径的球形范围内的其它流体质点粒子进行统计，确定区域流体的质量和大致体积，从而得到以目标质点粒子为中心的区域流体平均密度，作为该流体质点的密度属性；

(b)基于密度的流体质点区域压力计算

流体内的压力由密度计算得到，流体质点粒子所受到的压力加速度，由流体质点粒子的区域密度决定，由该区域内所有粒子的密度属性影响累计计算得到；

(c)基于密度和流体质点速度的流体质点区域粘度计算

流体内部存在阻碍流体发生相对运动的粘力，粘力方向总是与流体质点粒子的运动方向相反；。而流体粘力也与流体的密度相关，流体越稀，粘度越小，流体越稠，粘度越大；一个流体质点粒子的粘力加速度，由其影响半径内的所有粒子与其的相对运动速度，及流体粒子的区域密度累计计算得到；

(d)流体质点的加速度计算与流体粒子速度、位置的更新

通过计算得到流体内部质点粒子的压力加速度和粘力加速度，加上重力加速度，得到流体质点粒子的合力加速度；通过合力加速度更新质点粒子在当前帧的速度，再通过速度更新流体质点粒子当前帧的位置；新的位置将更新下一帧将要进行的区域密度计算，而新的速度则是下一帧新的粘力加速度的计算基础。

进一步地，步骤(2)具体包括：

(a)用空间网格分割质点粒子

根据流体的边界动态建立空间网格，网格步长为粒子的影响半径，即光滑核半径；根据粒子在流体中的相对位置，将粒子分配到空间网格内；对于单个质点粒子，其影响半径内的其它粒子称为邻居粒子，位于该单个质点粒子所处的网格，及该网格的26邻域网格内，不需要检索整个流体的全部质点粒子；

(b)动态分配网格空间

如果将流体所在最小矩形包围空间全部划分网格，并分配空间，在时间效率上能够提高质点粒子的检索速度，但在空间效率上会有大量空闲空间，特别是在流体体积较大，较不规则，粒子影响半径较小时，大量的网格数目可能会导致溢出；因此，定义一个映射表，根据粒子在流体中的相对位置计算出对应空间网格的xyz值，计算得到网格id，如果该网格已存在，则将粒子安置在该网格中，如果该网格不存在，在现有网格队列后插入该新网格，并得到队列序号，将网格id和队列序号绑定映射，从而使该网格可检索。

进一步地，步骤(3)具体包括：

(a)流体边界粒子判定与表面张力计算

将区域内有质点粒子的位置标记为1，没有质点粒子的位置标记为0，即可得到一个0-1高度图，对这个三维曲面求导，即可得到曲面的变化率；在流体内部，由于流体质点粒子密度较大，变化不明显；在流体的边界，流体质点粒子稀疏，梯度变化明显；位于明显梯度变化位置的质点粒子，即为该流体的边界粒子，受流体表面张力的约束；表面张力的作用方向，即该质点所在0-1场空间曲面梯度所指方向；

(b)接触角受力分析

将步骤(a)中得到的表面张力与已知的接触角属性结合，得到气液张力分量，作用于步骤(a)得到的边界粒子上，约束流体粒子的无限向外扩张，表现接近实际的接触角现象。

进一步地，步骤(4)具体包括：

(a)确定织物的分层结构

织物从整体到局部分为三部分，第一部分为织物整体，第二部分为纤维纱线，分为水平纱线和垂直纱线，第三部分为分段纤维，是对纱线进行的细分，作为织物仿真计算的最小对象；织物整体由水平和垂直纤维纱线组成，织物的面积大小由纱线数目决定；水平纱线与垂直纱线之间交错编织，形成织物导水 的纤维通道；纱线的长度由分段纤维数目决定，分段纤维具有固定长度，分段纤维数量越多，纱线越长；分段纤维作为最小的计算对象，由三角面片组成其表面，交互计算通过三角面完成；

(b)织物的储水设计

织物的储水主要在分段纤维内完成，分段纤维是固定长度、也是固定体积的最小计算对象，每个分段纤维的最大含水量由织物材质决定，即单位体积的织物材质能够储存的最大液体体积；

织物吸收流体后的饱和度变化，为分段纤维定义不同饱和度时的颜色表示，并根据每一帧分段纤维的饱和度改变该分段纤维的颜色，从而达到表达流体被吸收及在织物内部扩散的可视化效果。

进一步地，步骤(5)具体包括：

(a)同纤维扩散过程

在流体被织物吸收后，根据菲克定律从浓度高的区域向浓度低的区域在织物内扩散，首先发生的，也是速度最快的是同纤维扩散，即在同一条纱线纤维内扩散；

织物内的纱线由多个分段纤维组成，分段纤维是吸收流体的主体对象，也是储水的主体对象；流体分子在纱线内从高浓度区域向低浓度区域扩散的过程，即在饱和度高的分段纤维向饱和度低的分段纤维扩散的过程；

分段纤维之间两两首尾相连，相连的同纱线分段纤维能够直接发生流体分子扩散传递过程；

(b)相邻纤维扩散过程

相邻纤维扩散发生在不同纱线的分段纤维之间，从高浓度的分段纤维，向低浓度的分段纤维发生流体分子扩散；因为相邻分段纤维之间存在空气介质，扩散速率低于同纤维内扩散；

(c)高浓度差水分析出扩散过程

在织物远离流体与织物接触中心的位置，织物的纱线间隙较为干燥，每一个间隙由8条分段纤维围成，高浓度的分段纤维在纱线间隙会有少量流体分子析出，并通过间隙进入其它低浓度的分段纤维之中。

进一步地，步骤(6)具体包括：

(a)流体质点粒子与织物的碰撞设计

流体质点粒子与织物的碰撞，主要是与织物中分段纤维发生碰撞，而跟分 段纤维进行碰撞，即与构成分段纤维表面的三角面进行碰撞检测；

对质点粒子到三角面所在平面的垂直距离进行计算，判断距离是否满足碰撞误差条件；质点粒子到三角面所在平面的投影点，是否位于三角面的范围内；

如果质点粒子距离三角面所在平面不在误差允许范围内，说明质点粒子距离三角面还较为遥远，不足以发生碰撞；如果质点粒子到三角面所在平面的距离满足要求，然而投影不在三角面内，则粒子不会与三角面发生碰撞；进一步地，还包括以下步骤：

首先，要排除从纤维内部向外运动、从而与三角面发生碰撞的质点粒子；在检测碰撞时，需要检测质点粒子的速度方向，为了统一计算三角面法向量，三角面顶点全部按逆时针顺序存储；当质点粒子的速度方向与三角面所在平面的法向量夹角小于90度时，说明质点粒子是朝远离分段纤维的方向运动，碰撞不会发生；当质点粒子的速度方向与三角面的法向量夹角大于90度，则质点粒子是朝纤维内部运动的，可能发生碰撞，可以进行距离和投影判断；当质点粒子的速度方向与三角面的法向量垂直时，说明质点粒子平行于三角面进行运动，根据粒子到三角面的距离判断粒子是在三角面上滚动，或者只是在空中运动；

其次，流体内的质点粒子有时候运动速度较快，可能在单一帧内存在直接穿过纤维可能，从而造成漏判；预碰撞的处理方法在处理碰撞时，会先预计在当前帧速度下，粒子在下一帧有可能出现的空间位置，看看是否满足穿透的情况，如果满足穿透，对该流体粒子做上标记，在下一帧的时候，对该流体粒子进行吸收，从而在保持原有帧密度的情况下，减少因为小尺度仿真而可能导致的碰撞误差；

(b)流体质点粒子与织物的吸收判定

流体质点粒子与织物满足碰撞判定后，会进行吸收判定。吸收判定的主体仍是织物结构内的最小对象分段纤维；每个分段纤维都有自己的饱和吸收量和帧最大吸收量；

在满足流体质点粒子与分段纤维三角面的碰撞判定后，检查该碰撞三角面所在的分段纤维吸收是否达到饱和，如果达到饱和状态，吸收判定不成立，碰撞判定成立。质点粒子根据碰撞的速度被弹飞或者沿三角面在重力、摩擦力和张力的共同作用下滚动；

如果该分段纤维尚未达到饱和，但帧吸收量已经达到上限，质点粒子依旧不满足吸收判定，但是满足碰撞判定；

如果该分段纤维尚未达到饱和，帧吸收量也未达到上限，则质点粒子被该分段纤维吸收，该分段纤维的当前吸收量增加。

与现有技术相比，本发明具有如下优点和技术效果：从小尺度水体出发，主要针对三维空间内不同小尺度流体和不同材质的织物发生接触时的反应，和流体分子在织物内部的扩散过程做仿真，并将处理结果可视化显示出来。流体部分的仿真方法主要是基于拉格朗日方法中的光滑粒子流体动力学，并考虑了不同材质的织物和流体间的各种物理因素，如接触角，导水系数等，根据输入的计算帧数，计算流体从空中自由落体，到与织物发生碰撞吸收，吸收部分在织物纤维中扩散的全过程。最后根据每秒设定的播放帧数对计算结果进行显示，从而可以很直观地看到不同小尺度流体和不同织物的交互过程。

附图说明

图1为基于SPH的织物表面微观水扩散过程可视化方法的流程图。

图2为流体质点粒子的生命周期示意图。

图3为流体粒子与织物碰撞过程示意图。

具体实施方式

以下结合附图对本发明的具体实施作进一步说明，但本发明的实施和保护不限于此。

一种基于SPH的织物表面微观水扩散过程可视化方法，其包括如下步骤：

(13)构建SPH微观水流体；

(14)查找SPH微观水流体的空间网格；

(15)获取基于SPH微观水流体的表面张力和接触角；

(16)设计织物结构；

(17)设计基于织物结构的扩散通道；

(18)设计微观水流体与织物的碰撞与吸收。

步骤(13)具体包括：

(a)流体质点区域密度计算

流体质点粒子的区域密度通过搜索质点粒子附近的粒子数目累计计算得出；流体中的每一个质点粒子都有其作用力影响范围，随着距离的增大而影响力衰减，这个距离有一个极限，流体中的每一个质点粒子都有其作用力的影响距离极限标记为R，即影响力的作用半径，也是SPH中的光滑核半径；通过对以R为半径的球形范围内的其它流体质点粒子进行统计，确定区域流体的质量和大致体积，从而得到以目标质点粒子为中心的区域流体平均密度，作为该流体质点的密度属性；

(b)基于密度的流体质点区域压力计算

流体内的压力由密度计算得到，流体内部的压力总是从密度大的位置指向密度小的位置。流体质点粒子所在区域密度较大，会形成向外的正压，流体质点粒子所在区域密度较小，会形成向内的正压。所以，流体质点粒子所受到的压力加速度，由流体质点粒子的区域密度决定，由该区域内所有粒子的密度属性影响累计计算得到；

(c)基于密度和流体质点速度的流体质点区域粘度计算

流体内部存在阻碍流体发生相对运动的粘力，粘力方向总是与流体质点粒子的运动方向相反。流体质点粒子的相对运动速度越快，粒子间的粘力越大。而流体粘力也与流体的密度相关，流体越稀，粘度越小，流体越稠，粘度越大；一个流体质点粒子的粘力加速度，由其影响半径内的所有粒子与其的相对运动速度，及流体粒子的区域密度累计计算得到；

(d)流体质点的加速度计算与流体粒子速度、位置的更新

通过计算得到流体内部质点粒子的压力加速度和粘力加速度，加上重力加速度，得到流体质点粒子的合力加速度；通过合力加速度更新质点粒子在当前帧的速度，再通过速度更新流体质点粒子当前帧的位置；新的位置将更新下一帧将要进行的区域密度计算，而新的速度则是下一帧新的粘力加速度的计算基础。

步骤(14)具体包括：

(a)用空间网格分割质点粒子

为了提高影响半径内的粒子的搜索效率，根据流体的边界动态建立空间网格，网格步长为粒子的影响半径，即光滑核半径；根据粒子在流体中的相对位置，将粒子分配到空间网格内；对于单个质点粒子，其影响半径内的其它粒子 称为邻居粒子，位于该单个质点粒子所处的网格，及该网格的26邻域网格内，不需要检索整个流体的全部质点粒子；

(b)动态分配网格空间

如果将流体所在最小矩形包围空间全部划分网格，并分配空间，在时间效率上能够提高质点粒子的检索速度，但在空间效率上会有大量空闲空间，特别是在流体体积较大，较不规则，粒子影响半径较小时，大量的网格数目可能会导致溢出；因此，定义一个映射表，根据粒子在流体中的相对位置计算出对应空间网格的xyz值，计算得到网格id，如果该网格已存在，则将粒子安置在该网格中，如果该网格不存在，在现有网格队列后插入该新网格，并得到队列序号，将网格id和队列序号绑定映射，从而使该网格可检索。

步骤(15)具体包括：

(a)流体边界粒子判定与表面张力计算

将区域内有质点粒子的位置标记为1，没有质点粒子的位置标记为0，即可得到一个0-1高度图，对这个三维曲面求导，即可得到曲面的变化率；在流体内部，由于流体质点粒子密度较大，变化不明显；在流体的边界，流体质点粒子稀疏，梯度变化明显；位于明显梯度变化位置的质点粒子，即为该流体的边界粒子，受流体表面张力的约束；表面张力的作用方向，即该质点所在0-1场空间曲面梯度所指方向；

(b)接触角受力分析

液体与固体接触时，过固气液交界点作气液界面的切线，该切线与固体的夹角为该液体与固体的接触角。固液接触角是固气液交界处三力平衡的结果，包括固液张力，气液张力和固气张力力，气液张力的接触角余弦值与固液张力的和通固气张力达到平衡；在缺少气液张力分量的时候，或者说接触角为90度的时候，气液张力分量大小为0，流体会向外无限扩张；将步骤(a)中得到的表面张力与已知的接触角属性结合，得到气液张力分量，作用于步骤(a)得到的边界粒子上，约束流体粒子的无限向外扩张，表现接近实际的接触角现象。

步骤(16)具体包括：

(a)确定织物的分层结构

织物从整体到局部分为三部分，第一部分为织物整体，第二部分为纤维纱线，分为水平纱线和垂直纱线，第三部分为分段纤维，是对纱线进行的细分，作为织物仿真计算的最小对象；织物整体由水平和垂直纤维纱线组成，织物的 面积大小由纱线数目决定；水平纱线与垂直纱线之间交错编织，形成织物导水的纤维通道；纱线的长度由分段纤维数目决定，分段纤维具有固定长度，分段纤维数量越多，纱线越长；分段纤维作为最小的计算对象，由三角面片组成其表面，交互计算通过三角面完成；

(b)织物的储水设计

织物的储水主要在分段纤维内完成，分段纤维是固定长度、也是固定体积的最小计算对象，每个分段纤维的最大含水量由织物材质决定，即单位体积的织物材质能够储存的最大液体体积；

织物吸收流体后的饱和度变化，为分段纤维定义不同饱和度时的颜色表示，并根据每一帧分段纤维的饱和度改变该分段纤维的颜色，从而达到表达流体被吸收及在织物内部扩散的可视化效果。

步骤(17)具体包括：

(a)同纤维扩散过程

在流体被织物吸收后，根据菲克定律从浓度高的区域向浓度低的区域在织物内扩散，首先发生的，也是速度最快的是同纤维扩散，即在同一条纱线纤维内扩散；

织物内的纱线由多个分段纤维组成，分段纤维是吸收流体的主体对象，也是储水的主体对象；流体分子在纱线内从高浓度区域向低浓度区域扩散的过程，即在饱和度高的分段纤维向饱和度低的分段纤维扩散的过程；

分段纤维之间两两首尾相连，相连的同纱线分段纤维能够直接发生流体分子扩散传递过程；

(b)相邻纤维扩散过程

织物整体由水平纱线和垂直纱线交错编织而成，而纱线又由分段纤维构成，所以，处于不同纱线的分段纤维也有交织的可能。相邻纤维扩散发生在不同纱线的分段纤维之间，从高浓度的分段纤维，向低浓度的分段纤维发生流体分子扩散；因为相邻分段纤维之间存在空气介质，扩散速率低于同纤维内扩散；

(c)高浓度差水分析出扩散过程

在织物远离流体与织物接触中心的位置，织物的纱线间隙较为干燥，每一个间隙由8条分段纤维围成，高浓度的分段纤维在纱线间隙会有少量流体分子析出，并通过间隙进入其它低浓度的分段纤维之中。

步骤(18)具体包括：

(a)流体质点粒子与织物的碰撞设计

流体质点粒子与织物的碰撞，主要是与织物中分段纤维发生碰撞，而跟分段纤维进行碰撞，即与构成分段纤维表面的三角面进行碰撞检测；

对质点粒子到三角面所在平面的垂直距离进行计算，判断距离是否满足碰撞误差条件；质点粒子到三角面所在平面的投影点，是否位于三角面的范围内；

如果质点粒子距离三角面所在平面不在误差允许范围内，说明质点粒子距离三角面还较为遥远，不足以发生碰撞；如果质点粒子到三角面所在平面的距离满足要求，然而投影不在三角面内，则粒子不会与三角面发生碰撞；进一步地，还包括以下步骤：

首先，要排除从纤维内部向外运动、从而与三角面发生碰撞的质点粒子；在检测碰撞时，需要检测质点粒子的速度方向，为了统一计算三角面法向量，三角面顶点全部按逆时针顺序存储；当质点粒子的速度方向与三角面所在平面的法向量夹角小于90度时，说明质点粒子是朝远离分段纤维的方向运动，碰撞不会发生；当质点粒子的速度方向与三角面的法向量夹角大于90度，则质点粒子是朝纤维内部运动的，可能发生碰撞，可以进行距离和投影判断；当质点粒子的速度方向与三角面的法向量垂直时，说明质点粒子平行于三角面进行运动，根据粒子到三角面的距离判断粒子是在三角面上滚动，或者只是在空中运动；

其次，流体内的质点粒子有时候运动速度较快，可能在单一帧内存在直接穿过纤维可能，从而造成漏判；预碰撞的处理方法在处理碰撞时，会先预计在当前帧速度下，粒子在下一帧有可能出现的空间位置，看看是否满足穿透的情况，如果满足穿透，对该流体粒子做上标记，在下一帧的时候，对该流体粒子进行吸收，从而在保持原有帧密度的情况下，减少因为小尺度仿真而可能导致的碰撞误差；

(b)流体质点粒子与织物的吸收判定

流体质点粒子与织物满足碰撞判定后，会进行吸收判定。吸收判定的主体仍是织物结构内的最小对象分段纤维；每个分段纤维都有自己的饱和吸收量和帧最大吸收量；

在满足流体质点粒子与分段纤维三角面的碰撞判定后，检查该碰撞三角面所在的分段纤维吸收是否达到饱和，如果达到饱和状态，吸收判定不成立，碰撞判定成立。质点粒子根据碰撞的速度被弹飞或者沿三角面在重力、摩擦力和张力的共同作用下滚动；

如果该分段纤维尚未达到饱和，但帧吸收量已经达到上限，质点粒子依旧不满足吸收判定，但是满足碰撞判定；

如果该分段纤维尚未达到饱和，帧吸收量也未达到上限，则质点粒子被该分段纤维吸收，该分段纤维的当前吸收量增加。

本实例所用到的部分算法和模型举例具体如下。

1.1光滑粒子流体动力学模型：

1.1.1应用基本场域函数：

及哈密顿形式与拉普拉斯形式：

其中As表示目标计算属性值，m_j表示流体质点粒子质量，A_j是流体质点粒子与目标计算属性值同属性的数值，ρ_j表示流体质点粒子j周围的区域密度，W为光滑核函数，r为目标位置的空间坐标的向量形式，r_j为邻居粒子的空间坐标的向量形式，h为光滑核半径。

其中光滑核函数满足W(r,h)＝W(-r,h)和∫W(r)dr＝1的条件。

1.1.2步骤(1)中的区域密度计算方法：

其中ρ表示流体区域密度，m_j表示流体质点粒子质量，r为目标位置的空间坐标的向量形式，r_j为邻居粒子的空间坐标的向量形式，h为光滑核半径。

1.1.3步骤(1)中的压力加速度计算方法：

其中a_i ^pressure表示流体质点粒子i的压力加速度，p为目标区域压力大小，m_j表示流体质点粒子质量，ρ_i和ρ_j表示流体质点粒子i和j周围的区域密度，r为目标位置的空间坐标的向量形式，r_j为邻居粒子的空间坐标的向量形式，h为光滑核半径。

1.1.4步骤(1)中的粘度加速度计算应用：

其中a_i ^viscosity表示流体质点粒子i的黏力加速度，F_i为目标区域黏力大小，μ是黏力常量，m_j表示流体质点粒子质量，u_i和u_j表示流体质点粒子i和j相对于目标位置的运动速度。ρ_i和ρ_j表示流体质点粒子i和j周围的区域密度，r为目标位置的空间坐标的向量形式，r_j为邻居粒子的空间坐标的向量形式，h为光滑核半径。

1.1.5步骤(2)中的空间网格方法：

对于一个流体在空间中的位置可用一个长方体包围，根据光滑粒子流体动力学模型中的光滑核半径对长方体进行切分，将流体所在空间切分成x*y*z的复数个小长方体。根据流体质点在这个长方体空间内的相对位置计算得到流体质点应在的网格序号id，通过哈希表寻找网格。如果网格存在，在该网格内添加该流体质点id，如果网格不存在，创建网格并添加该流体质点id。

查找流体质点的邻接粒子，首先根据流体在空间内的相对位置得到网格id(x,y,z)，首先验算该网格内的其余粒子的位置关系，然后，依次验算26邻域的网格中的粒子：(x-1,y,z),(x+1,y,z),(x,y-1,z)....,(x+1,y+1,z+1)。

1.1.6步骤(3)中的表面张力计算方法：

通过建立标量场计算，将有流体质点粒子的位置标志为1，不存在质点粒子的位置标志为0，由此建立表示空间位置是否存在流体粒子的标量场：

其中Cs表示0-1标量场数值，m_j表示流体质点粒子质量，ρ_j表示流体质点粒子j周围的区域密度，r为目标位置的空间坐标的向量形式，r_j为邻居粒子的空间坐标的向量形式，h为光滑核半径。

对标量场函数求导得到标量场梯度值：

梯度变化剧烈的位置即标量场的0-1边界，即流体边界，由此得到流体边界粒子，并由梯度指向得到表面张力方向。

最终表面张力的计算结果可得：

其中a_i ^{surfaceTension}表示流体质点粒子i的表面张力，σ是张力常量，κ为标量场Cs的拉普拉斯形式，表示曲率大小用以衡量张力大小，ρ_i表示流体质点粒子i周围的区域密度。

步骤(4)的扩散方法：

根据菲克第二定律可以得到流体粒子在织物内部的扩散速率：

其中，D是扩散系数，主要取决于织物对流体的亲和力大小，J是扩散通量，，C是浓度，x是距离，负号表示扩散方向为浓度梯度的反方向，浓度梯度越大，扩散通量越大。

由于单位距离扩散通量越大，浓度下降越快，有代入上式可得

步骤(5)中的碰撞处理：

1.2.包围盒：

包围盒为AABB包围盒，通过统计纤维分段内三角面片的x,y,z方向上的最大值和最小值，得到纤维分段的包围盒，通过统计纤维内纤维分段的包围盒顶点数值得到纤维的包围盒，通过统计织物内纤维包围盒的顶点属性得到织物的整体包围盒。

1.3.四层碰撞检测筛选：

首先与织物包围盒进行检测，然后与纤维包围盒进行检测，然后与纤维分段包围盒进行检测，最后与三角面片进行碰撞检测。

反向碰撞排除：

如果粒子速度矢量V与三角面片所在平面的法向量N夹角小于90度，即粒子从三角面片背面与三角面发生碰撞，反之，从正面发生碰撞。系统内否认从背面发生碰撞的合理性。

1.4.粒子投影是否在三角面内的判断方法：

三角面顶点A,B,C，粒子在三角面上投影点P，将P与A,B,C分别连线，得到三角形PAB，PBC，PCA，同时可以很直接获得线段AB，BC，AC，PA，PB，PC的长。根据海伦秦九韶公式可以计算得到三角形面积S_ABC，S_PAB，S_PBC，S_PAC。如果S_PAB+S_PBC+S_PAC>S_ABC,则点P位于三角面ABC之外，反正，则在三角面ABC内。

2.方法的前期准备：

需要准备实验对象流体和织物材质的相关物理性质，包括流体密度，粘力，流体与织物的摩擦力，流体与织物的接触角，流体在织物内的扩散系数，接触角的取值范围一般为0-180度，并计算根据流体质点密度计算得流体质点质量。

3.试验结果：

根据织物与流体的接触角主要有以下几种情况：

a.θ＝0，完全润湿

b.θ<π/2，部分润湿或者完全润湿

c.θ＝π/2，是否润湿临界线

d.θ>π/2，不润湿

e.θ＝π，完全不润湿

不浸润的情况下，织物的纤维分段不会出现颜色改变。当织物与流体的接触角满足浸润条件时，根据纤维分段的含水量多少会有从浅蓝到深蓝色的梯度变化，颜色越深，代表含水量越多。

根据织物与流体的接触角以及织物与流体的扩散系数，应该有不同的扩散速率，在扩散通量较大时，流体在织物内的扩散速率会更快，最终都会在织物内达到平衡。

Claims (7)

1.一种基于SPH的织物表面微观水扩散过程可视化方法，其特征在于包括如下步骤：

(1)构建SPH微观水流体；

(2)查找SPH微观水流体的空间网格；

(3)获取基于SPH微观水流体的表面张力和接触角；

(4)设计织物结构；

(5)设计基于织物结构的扩散通道；

(6)设计微观水流体与织物的碰撞与吸收。

2.根据权利要求1所述的基于SPH的织物表面微观水扩散过程可视化方法，其特征在于步骤(1)具体包括：

(a)流体质点区域密度计算

流体质点粒子的区域密度通过搜索质点粒子附近的粒子数目累计计算得出；流体中的每一个质点粒子都有其作用力的影响距离极限标记为R，即影响力的作用半径，也是SPH中的光滑核半径；通过对以R为半径的球形范围内的其它流体质点粒子进行统计，确定区域流体的质量和大致体积，从而得到以目标质点粒子为中心的区域流体平均密度，作为该流体质点的密度属性；

(b)基于密度的流体质点区域压力计算

流体内的压力由密度计算得到，流体质点粒子所受到的压力加速度，由流体质点粒子的区域密度决定，由该区域内所有粒子的密度属性影响累计计算得到；

(c)基于密度和流体质点速度的流体质点区域粘度计算

流体内部存在阻碍流体发生相对运动的粘力，粘力方向总是与流体质点粒子的运动方向相反；而流体粘力也与流体的密度相关，流体越稀，粘度越小，流体越稠，粘度越大；一个流体质点粒子的粘力加速度，由其影响半径内的所有粒子与其的相对运动速度，及流体粒子的区域密度累计计算得到；

(d)流体质点的加速度计算与流体粒子速度、位置的更新

通过计算得到流体内部质点粒子的压力加速度和粘力加速度，加上重力加速度，得到流体质点粒子的合力加速度；通过合力加速度更新质点粒子在当前帧的速度，再通过速度更新流体质点粒子当前帧的位置；新的位置将更新下一帧将要进行的区域密度计算，而新的速度则是下一帧新的粘力加速度的计算基础。

3.根据权利要求1所述的基于SPH的织物表面微观水扩散过程可视化方法，其特征在于步骤(2)具体包括：

(a)用空间网格分割质点粒子

根据流体的边界动态建立空间网格，网格步长为粒子的影响半径，即光滑核半径；根据粒子在流体中的相对位置，将粒子分配到空间网格内；对于单个质点粒子，其影响半径内的其它粒子称为邻居粒子，位于该单个质点粒子所处的网格，及该网格的26邻域网格内，不需要检索整个流体的全部质点粒子；

(b)动态分配网格空间

如果将流体所在最小矩形包围空间全部划分网格，并分配空间，在时间效率上能够提高质点粒子的检索速度，但在空间效率上会有大量空闲空间，特别是在流体体积较大，较不规则，粒子影响半径较小时，大量的网格数目可能会导致溢出；因此，定义一个映射表，根据粒子在流体中的相对位置计算出对应空间网格的xyz值，计算得到网格id，如果该网格已存在，则将粒子安置在该网格中，如果该网格不存在，在现有网格队列后插入该新网格，并得到队列序号，将网格id和队列序号绑定映射，从而使该网格可检索。

4.根据权利要求1所述的基于SPH的织物表面微观水扩散过程可视化方法，其特征在于步骤(3)具体包括：

(a)流体边界粒子判定与表面张力计算

将区域内有质点粒子的位置标记为1，没有质点粒子的位置标记为0，即可得到一个0-1高度图，对这个三维曲面求导，即可得到曲面的变化率；在流体内部，由于流体质点粒子密度较大，变化不明显；在流体的边界，流体质点粒子稀疏，梯度变化明显；位于明显梯度变化位置的质点粒子，即为该流体的边界粒子，受流体表面张力的约束；表面张力的作用方向，即该质点所在0-1场空间曲面梯度所指方向；

(b)接触角受力分析

将步骤(a)中得到的表面张力与已知的接触角属性结合，得到气液张力分量，作用于步骤(a)得到的边界粒子上，约束流体粒子的无限向外扩张，表现接近实际的接触角现象。

5.根据权利要求1所述的基于SPH的织物表面微观水扩散过程可视化方法，其特征在于步骤(4)具体包括：

(a)确定织物的分层结构

织物从整体到局部分为三部分，第一部分为织物整体，第二部分为纤维纱线，分为水平纱线和垂直纱线，第三部分为分段纤维，是对纱线进行的细分，作为织物仿真计算的最小对象；织物整体由水平和垂直纤维纱线组成，织物的面积大小由纱线数目决定；水平纱线与垂直纱线之间交错编织，形成织物导水的纤维通道；纱线的长度由分段纤维数目决定，分段纤维具有固定长度，分段纤维数量越多，纱线越长；分段纤维作为最小的计算对象，由三角面片组成其表面，交互计算通过三角面完成；

(b)织物的储水设计

织物的储水主要在分段纤维内完成，分段纤维是固定长度、也是固定体积的最小计算对象，每个分段纤维的最大含水量由织物材质决定，即单位体积的织物材质能够储存的最大液体体积；

织物吸收流体后的饱和度变化，为分段纤维定义不同饱和度时的颜色表示，并根据每一帧分段纤维的饱和度改变该分段纤维的颜色，从而达到表达流体被吸收及在织物内部扩散的可视化效果。

6.根据权利要求1所述的基于SPH的织物表面微观水扩散过程可视化方法，其特征在于步骤(5)具体包括：

(a)同纤维扩散过程

在流体被织物吸收后，根据菲克定律从浓度高的区域向浓度低的区域在织物内扩散，首先发生的，也是速度最快的是同纤维扩散，即在同一条纱线纤维内扩散；

织物内的纱线由多个分段纤维组成，分段纤维是吸收流体的主体对象，也是储水的主体对象；流体分子在纱线内从高浓度区域向低浓度区域扩散的过程，即在饱和度高的分段纤维向饱和度低的分段纤维扩散的过程；

分段纤维之间两两首尾相连，相连的同纱线分段纤维能够直接发生流体分子扩散传递过程；

(b)相邻纤维扩散过程

相邻纤维扩散发生在不同纱线的分段纤维之间，从高浓度的分段纤维，向低浓度的分段纤维发生流体分子扩散；因为相邻分段纤维之间存在空气介质，扩散速率低于同纤维内扩散；

(c)高浓度差水分析出扩散过程

在织物远离流体与织物接触中心的位置，织物的纱线间隙较为干燥，每一个间隙由8条分段纤维围成，高浓度的分段纤维在纱线间隙会有少量流体分子析出，并通过间隙进入其它低浓度的分段纤维之中。

7.根据权利要求1所述的基于SPH的织物表面微观水扩散过程可视化方法，其特征在于步骤(6)具体包括：

(a)流体质点粒子与织物的碰撞设计

流体质点粒子与织物的碰撞，主要是与织物中分段纤维发生碰撞，而跟分段纤维进行碰撞，即与构成分段纤维表面的三角面进行碰撞检测；

对质点粒子到三角面所在平面的垂直距离进行计算，判断距离是否满足碰撞误差条件；质点粒子到三角面所在平面的投影点，是否位于三角面的范围内；

如果质点粒子距离三角面所在平面不在误差允许范围内，说明质点粒子距离三角面还较为遥远，不足以发生碰撞；如果质点粒子到三角面所在平面的距离满足要求，然而投影不在三角面内，则粒子不会与三角面发生碰撞；进一步地，还包括以下步骤：

首先，要排除从纤维内部向外运动、从而与三角面发生碰撞的质点粒子；在检测碰撞时，需要检测质点粒子的速度方向，为了统一计算三角面法向量，三角面顶点全部按逆时针顺序存储；当质点粒子的速度方向与三角面所在平面的法向量夹角小于90度时，说明质点粒子是朝远离分段纤维的方向运动，碰撞不会发生；当质点粒子的速度方向与三角面的法向量夹角大于90度，则质点粒子是朝纤维内部运动的，可能发生碰撞，可以进行距离和投影判断；当质点粒子的速度方向与三角面的法向量垂直时，说明质点粒子平行于三角面进行运动，根据粒子到三角面的距离判断粒子是在三角面上滚动，或者只是在空中运动；

其次，流体内的质点粒子有时候运动速度较快，可能在单一帧内存在直接穿过纤维可能，从而造成漏判；预碰撞的处理方法在处理碰撞时，会先预计在当前帧速度下，粒子在下一帧有可能出现的空间位置，看看是否满足穿透的情况，如果满足穿透，对该流体粒子做上标记，在下一帧的时候，对该流体粒子进行吸收，从而在保持原有帧密度的情况下，减少因为小尺度仿真而可能导致的碰撞误差；

(b)流体质点粒子与织物的吸收判定

流体质点粒子与织物满足碰撞判定后，会进行吸收判定。吸收判定的主体仍是织物结构内的最小对象分段纤维；每个分段纤维都有自己的饱和吸收量和帧最大吸收量；

在满足流体质点粒子与分段纤维三角面的碰撞判定后，检查该碰撞三角面所在的分段纤维吸收是否达到饱和，如果达到饱和状态，吸收判定不成立，碰撞判定成立。质点粒子根据碰撞的速度被弹飞或者沿三角面在重力、摩擦力和张力的共同作用下滚动；

如果该分段纤维尚未达到饱和，但帧吸收量已经达到上限，质点粒子依旧不满足吸收判定，但是满足碰撞判定；

如果该分段纤维尚未达到饱和，帧吸收量也未达到上限，则质点粒子被该分段纤维吸收，该分段纤维的当前吸收量增加。