水场景目标红外尾迹特性的流体模拟方法

摘要

本发明公开了一种水场景目标红外尾迹特性的流体模拟方法。其步骤为：1)对场景中的模型进行体素化产生固体粒子；2)在存在液体的区域添加液体粒子；3)对每一时间帧，利用SPH对流体的力学特性进行数值模拟；4)对每一时间帧，利用瞬态热传导方程对流体的热力学特性进行数值模拟及红外特征图形绘制。本发明解决了流动介质中计算热力学特性的需求。本方法能够真实模拟水场景中的力学特征和红外特征，并在实验室环境下成功模拟了近年新发现的红外现象：船只在航行过程产生的红外尾迹。

说明书

技术领域

本发明涉及流体模拟的数值计算领域与热力学的数值计算领域。具体涉及一 种水场景目标红外尾迹特性的流体模拟方法。

背景技术

近年来，随着红外热像仪成像分辨率的提高，在遥感及航空拍摄中，发现了 船只航行时的红外新特征：船只在水场景中航行时会留下红外尾迹。然而现有的 热力学计算方法普遍针对于静态物质进行数值模拟计算，并不能很好的模拟这一 特征；而流体模拟的数值计算领域的方法也主要针对于计算模拟流体的力学特征， 无法模拟热力学的热传导特征。因此，难以模拟计算并绘制出该红外特征现象。

下面介绍已有的流体模拟的数值计算方法与热力学的数值计算方法：

1)热力学的数值计算方法

目前在国内外的瞬态热传导问题的数值方法，主要基于网格的有限元法、有 限差分法和有限体积法。近年来也提出了基于粒子的无网格方法。然而研究过程 中，都没有考虑液体运动对其造成的影响，因此也无法模拟水体目标运动造成的 红外现象。

2)流体模拟的数值计算方法

这类方法主要可以细分成为三类：欧拉法、拉格朗日法和混合法。欧拉法是 一种基于网格的方法，它使用2D或者3D的网格，从流体空间中的各个固定点 着手，分析每个固定点上的流体速度、压强、密度等参数随时间和空间的变化。 拉格朗日法是一种基于粒子的方法，用一系列遵循物理规则的粒子来表示流体， 把单个流体粒子作为研究对象，研究其运动要素(位置、速度等)的变化过程，并 通过综合各个流体粒子的运动来获得一定空间内所有流体质点的运动规律。欧拉 法与拉格朗日法混合的主要思想是对三维水体用欧拉网格法表示，再用拉格朗日 法进行小尺度的泡沫和浪花等效果的模拟。这些方法能够很好的模拟水体的力学 性质，并实现良好的可见光水体模拟效果，但无法得到热力学性质，因此无法模 拟红外特征现象。

发明内容

本发明的目的在模拟水场景目标在水中的力学特征和红外特征，解决现有流 体方法无法模拟红外特征和现有热力学方法无法模拟水体流动特征的问题，提供 一种水场景中目标红外尾迹特性的流体模拟方法。

一种水场景目标红外尾迹特性的流体模拟方法包括以下步骤：

1)对场景中的模型进行体素化产生固体粒子；

2)在存在液体的区域添加液体粒子；

3)对每一时间帧，利用SPH对流体的力学特性进行数值模拟；

4)对每一时间帧，利用瞬态热传导方程对流体的热力学特性进行数值模拟 及红外特征图形绘制。

所述步骤1)为：

(2.1)根据模型顶点位置数据计算模型的AABB包围盒C_b；

(2.2)将C_b平均划分为X*Y*Z个子正方体C_i,j,k

其中，X、Y、Z为空间在xyz三个方向上的分辨率，i∈[0,X),j∈[0,Y), k∈[0,Z)；

(2.3)对每个C_i,j,k与模型进行求交，若相交则该区域存在固体体素，否则 该区域固体体素为空，每个求交计算所得的固体体素标记为S_i；

(2.4)根据模型的运动状态，设定S_i每一时刻的速度值。

所述步骤2)为：

(3.1)设定待模拟的水场景液体表面曲面方程为F(p)，其中p∈(R,R,R)， R为有理数集合，F(p)<0表示点p在液体内部，F(p)＝0表示点p在液体边界， F(p)>0表示点p在液体外部，

(3.2)将空间平均划分为X`*Y`*Z`个子正方体C`_i,j,k

其中，X`、Y`、Z`为空间在xyz三个方向上的液体分辨率，i∈[0,X`),j ∈[0,Y`),k∈[0,Z`)；

(3.3)对每个C`_i,j,k，计算F(i,j,k)，若小于等于0则该区域存在液体粒 子，若大于0则该区域不存在液体粒子,每个液体粒子标记为L_i；

(3.4)根据水场景的运动特性与水体温度分布初始化L_i的力学参数与热力 学参数：

(a)根据初始化条件，初始化每个L_i的速度；

(b)根据温度初始化方程T(h)＝T₀–t_a*h，初始化每个L_i的初 始温度；

其中，h为L_i的水深，T₀为液体表面温度，t_a为温度随深度的变化参 数。

所述步骤3)为：

(4.1)对于每一个粒子L_i，根据SPH方程定义，利用所有光滑核半径内的 粒子L_j或S_j计算液体的密度ρ_i：

ρ_i＝ρ(r_i)＝315m/(64πh⁹)∑_j(h²-|r_i-r_j|²)²

其中r_i是L_i在空间中的位置，m是L_i的质量，h为光滑核半径，r_j为粒子L_j或S_j的空间位置；

(4.2)对于每一个粒子L_i，根据理想气体状态方程计算液体的压强：

p_iM＝ρ_iRT_i

其中ρ_i，T_i，p_i分别为L_i的密度、温度、压强，M为平均摩尔质量，R为理 想气体常数；

(4.3)根据SPH公式，利用所有光滑核半径内的粒子L_j或S_j计算每一个粒 子L_i的加速度a_i：

a_i＝a(r_i)＝g+m*45/(πh⁶)∑_j((p_i+p_j)/(2ρ_iρ_j)*(h-r)²*(r_i-r_j)/r)+mμ

*45/(πh⁶)∑_j(u_j-u_i)/(ρ_iρ_j)*(h-r)

其中，r＝|r_i-r_j|，ρ_j是粒子L_j或S_j的密度，a_j为粒子L_j或S_j的加速度，u_j为粒子L_j或S_j的速度，a_i，u_i分别为粒子L_i的加速度和速度；

(4.4)利用加速度与速度，根据牛顿第二定律计算模拟每个粒子L_i或S_i的 运动轨迹，重复步骤(4.1)～(4.4)即可计算得到每一时刻每个粒子的力学特 性。

所述步骤4)为：

(5.1)对于每个粒子L_i，根据笛卡尔坐标系下的热传导方程，利用狄利克 雷边界条件计算以下方程得到温度变化：

${\mathrm{dT}}_i/dt={\Sigma }_j{m}_j/{\rho }_j*(q_j+q_i)(▿W_{ij})+\int (q+q_i)▿Wdx$

$q_i={\Sigma }_j{m}_j/{\rho }_j*(T_j-T_i)(▿W_{ij})+\int (T-T_i)▿Wdx$

其中W为光滑粒子核函数W(x)，W_ij＝W(|r_i-r_j|)，q，T为积分点的热流密度 与温度；

(5.2)根据计算得到的温度场T的分布，利用Planck公式计算以下方程得 到其红外辐射量E：

$E={ϵ}_0·{\int }_{{\lambda }_1}^{{\lambda }_2}\frac{C_1}{{\lambda }^5}{e}^{\frac{-C_2}{\lambda T}}d\lambda$

其中ε₀为材料的表面发射率；C₁为第一辐射常数，其值为3.742×10^-16W·m²； C₂为第二辐射常数，其值为1.4388×10^-2m·K。λ₁和λ₂为红外探测仪的探测波 段的两段波长；

(5.3)若场景中最高和最低辐射温度分别为T_max和T_min，由上式算得的相 应的辐射强度分别为E_max和E_min，则对于辐射强度为E_i的表面，其对应的计算灰 度值为：

$G_i=\frac{E-E_{min}}{E_{max}-E_{\min }}\times 255$

利用计算机图形学中的颜色均衡方法，绘制得到红外特征图形效果。

本发明的有点在于：

传统的流体模拟方法，其基于流体力学分析，从而求得了流体的力学性质， 但无法模拟热力学性质产生的红外现象；传统的热力学模拟方法，基于热传导方 程，可计算物质的瞬态热力学性质，从而模拟红外现象，但无法模拟流体的运动 现象。

本发明的方法提出了水场景目标红外尾迹特性的流体模拟方法。通过体素化 将场景划分为粒子，利用SPH估值近似求解偏微分方程，同时求解流体的力学性 质与热力学性质。这样我们就产生了可以同时具有力学特征和红外特征的水场景。

本方法提出了一种使用体素化统一划分场景的方法，来使得水场景中的液体 和固体能够使用统一的方程进行物理计算，避免了复杂的边界条件计算。

本方法提出了一种使用SPH方法同时计算流体的力学性质与热力学性质的 方法。同时模拟的液体的力学性质与热力学性质，能同时模拟水场景的红外特征。

总之，本发明给出了一种水场景目标红外尾迹特性的流体模拟方法，该方法 对比目前方法，能够同时计算水场景的力学性质与热力学性质，可以模拟绘制出 其尾迹的红外图形特征效果。

附图说明

图1为船只在海面航行时的液体分布图；

图2为船只在海面航行时的红外特征俯视图形。

具体实施方式

水场景目标红外尾迹特性的流体模拟方法包括以下步骤：

1)对场景中的模型进行体素化产生固体粒子；

2)在存在液体的区域添加液体粒子；

3)对每一时间帧，利用SPH对流体的力学特性进行数值模拟；

4)对每一时间帧，利用瞬态热传导方程对流体的热力学特性进行数值模拟 及红外特征图形绘制；

所述步骤1)为：

(1.1)根据模型顶点位置数据计算模型的AABB包围盒C_b；

(1.2)将C_b平均划分为X*Y*Z个子正方体C_i,j,k

其中，X、Y、Z为空间在xyz三个方向上的分辨率，i∈[0,X),j∈[0,Y), k∈[0,Z)；

(1.3)对每个C_i,j,k求与模型进行求交，若相交则该区域存在固体体素，否 则该区域固体体素为空，每个计算所得的固体体素标记为S_i。

(1.4)根据模型的运动状态，设定S_i每一时刻的速度值；

所述步骤2)为：

(2.1)使用方程F(p)表达液体表面方程

其中p∈(R,R,R)，F(p)<0表示点p在液体内部，F(p)＝0表示点p在液体 边界，F(p)>0表示点p在液体外部。F(p)的方程形式可根据实际实验需要选取；

(2.2)将空间平均划分为X`*Y`*Z`个子正方体C`_i,j,k

其中，X`、Y`、Z`为空间在xyz三个方向上的液体分辨率，i∈[0,X`),j ∈[0,Y`),k∈[0,Z`)；

(2.3)对每个C_i,j,k，计算F(i,j,k)，若大于0则该区域存在液体粒子，若 小于0则该区域不存在液体粒子,每个液体粒子标记为L_i；

(2.4)根据水场景的运动特性与水体温度分布初始化L_i的力学参数与热力 学参数：

(a)根据初始化条件，初始化每个L_i的速度。若无特殊要求，L_i的 速度初始化为0(水面静止)；

(b)根据温度初始化方程T(h)＝T₀–t_a*h，初始化每个L_i的初 始温度；

其中，h为L_i的水深，T₀为液体表面温度，t_a为温度随深度的变化参 数，T₀与t_a根据实际实验需要选取。

所述步骤3)为：

(3.1)对于每一个粒子L_i，根据SPH方程定义，利用所有光滑核半径内的 粒子L_j或S_j计算液体的密度：

ρ_i＝ρ(r_i)＝315m/(64πh⁹)∑_j(h²-|r_i-r_j|²)²

其中r_i是L_i在空间中的位置，m是L_i的质量，h为光滑核半径，r_j为粒子L_j或S_j的空间位置；

(3.2)对于每一个粒子L_i，根据理想气体状态方程计算液体的压强：

p_iM＝ρ_iRT_i

其中ρ_i，T_i，p_i分别为L_i的密度、温度、压强，M为平均摩尔质量，R为比 例常数，约为8.31441±0.00026J/(mol·K)；

(3.3)根据SPH公式，利用所有光滑核半径内的粒子L_j或S_j计算每一个粒 子L_i的加速度：

a_i＝a(r_i)＝g+m*45/(πh⁶)∑_j((p_i+p_j)/(2ρ_iρ_j)*(h-r)²*(r_i-r_j)/r)+mμ

*45/(πh⁶)∑_j(u_j-u_i)/(ρ_iρ_j)*(h-r)

其中，r＝|r_i-r_j|，a_i，u_i分别为粒子L_i在的加速度，速度，a_j，u_j分别为粒 子L_j或S_j的加速度或速度；

(3.4)利用加速度与速度，根据牛顿第二定律计算模拟每个粒子L_i或S_i的 运动轨迹，重复步骤(1)～(4)即可计算得到每一时刻每个粒子的力学特性。 海浪形状稳定后的船只在海面航行时的液体分布图如图1所示；

所述步骤4)为：

(4.1)对于每个粒子L_i，根据笛卡尔坐标系下的热传导方程，利用狄利克 雷边界条件计算以下方程得到温度变化：

${\mathrm{dT}}_i/dt={\Sigma }_j{m}_j/{\rho }_j*(q_j+q_i)(▿W_{ij})+\int (q+q_i)▿Wdx$

$q_i={\Sigma }_j{m}_j/{\rho }_j*(T_j-T_i)(▿W_{ij})+\int (T-T_i)▿Wdx$

其中W为光滑粒子核函数W(x)，W_ij＝W(|r_i-r_j|)，可根据实验需要选用，本 利中的计算结果采用lucy核函数进行数值计算。q，T为积分点的热流密度与温 度。

(4.2)根据计算得到的温度场T的分布，利用Planck公式计算以下方程得 到其红外辐射量E：

$E={ϵ}_0·{\int }_{{\lambda }_1}^{{\lambda }_2}\frac{C_1}{{\lambda }^5}{e}^{\frac{-C_2}{\lambda T}}d\lambda$

其中ε₀为材料的表面发射率；C₁为第一辐射常数，其值为3.742×10^-16W·m²； C₂为第二辐射常数，其值为1.4388×10^-2m·K。λ₁和λ₂为红外探测仪的探测波 段的两段波长。

(4.3)若场景中最高和最低辐射温度分别为T_max和T_min，由上式算得的相 应的辐射强度分别为E_max和E_min(ε₀的值为1)。则对于辐射强度为E_i的表面， 其对应的计算灰度值为：

$G_i=\frac{E-E_{min}}{E_{max}-E_{\min }}\times 255$

利用计算机图形学中的颜色均衡方法，本实施例中采用GouraudShading明 暗处理方法，可绘制得到红外特征图形结果，结果如图2所示。