CN111125922B - 一种冷空间内下落液滴冻结过程的计算方法 - Google Patents

Abstract

本发明公开了一种冷空间内下落液滴冻结过程的计算方法，在焓法格子Boltzmann相变模型的基础上进行修正，构建了冷空间内下落液滴冻结过程计算方法。该方法使用密度分布函数和温度分布函数分别模拟流体的速度场和温度场，通过求解固相体积分数追踪相变过程中各相的体积分布情况，使用固相体积分数修正焓值以及流场演化方程迁移步，并考虑相变潜热和温差相对较大(大40K)时辐射换热对相变过程的影响。

Description

一种冷空间内下落液滴冻结过程的计算方法

技术领域

本发明属于液滴冻结过程计算领域，尤其是涉及一种冷空间内下落液滴冻结过程的计算方法。

背景技术

下落液滴冻结现象普遍存在于自然界、制冷、空气调节、航空航天等领域。因此，掌握具体冻结过程对于工程应用来说具有十分重要的研究意义。以喷雾结晶和人工降雪为研究背景，从传热传质以及相变理论出发分析液滴冻结过程。当液滴温度被冷却至一定温度以下时，液滴会发生相变，即由液态转变为固态，该过程称为液滴冻结过程。液滴由液态转变为固态时的转变温度称为液滴冻结温度。在湿空气近冷壁面结霜机理中，液滴冻结过程是整个结霜过程的开始，具体可以概括为：液滴在冷环境中被冷却到冻结温度以下，液滴表面先形成冰壳，液滴由外向内逐渐发生冻结现象，最终整个液滴被冻结成冰粒。当液滴在冷空间内下落时，液滴被周围环境冷却发生冻结现象，液滴表面先形成冰壳，假定该冰壳为刚性，冰壳的存在将阻碍流体流动，但是换热过程依然继续。当液滴与周围环境温差相对较大(大于40K)时，辐射换热对相变过程有一定影响。因此，运动液滴冻结过程将同时受热传导、热对流以及热辐射影响。此外，液滴内部固液相界面处也将释放相变潜热。然而，由于液滴内部固液相界面不规则变化、固液两相区结构复杂以及温度场影响，液滴冻结过程的试验研究是一项极具挑战性和高成本的任务。虽然试验方法可以得到液滴冻结过程的温度分布，但是试验设备造成的误差或干扰对实验结果有很大影响，例如，当利用热电偶测量体积较小的液滴表面或内部温度时，热电偶的热阻会影响液滴的换热过程，因此热电偶对实验结果的准确性将产生一定影响。此外，观测用来判断液滴是否完全冻结的固液相界面位置也十分困难。因此，需要借助其他手段进行辅助研究。目前用来计算非等温冻结过程的方法主要是Jiaung提出的焓法格子Boltzmann模型。对于运动液滴冻结过程而言，该过程是一个涉及流体流动的非等温相变过程，目前用来解决非等温相变过程的数学模型大多未将冻结过程对流体流动的影响考虑在内。此外，当液滴与周围环境间的温差相对较大(大于40K)时，辐射换热对冻结过程也会产生一定影响，但大部分计算相变过程的格子Boltzmann模型都未考虑辐射换热。在目前的研究成果中，关于冷空间内下落液滴冻结过程中温度场、速度场以及液滴冻结状态的研究也十分有限。因此需要建立新的数学模型来解决这一问题。

发明内容

为了更好地解决冷空间内悬浮液滴以及下落液滴冻结问题，得到不同时刻液滴内温度场、液滴冻结状态、冻结参数及影响因素，本发明给出了一种冷空间内下落液滴冻结过程的计算方法。

针对上述问题，本发明给出了一种冷空间内下落液滴冻结过程的计算方法，在焓法格子Boltzmann相变模型的基础上进行修正，构建了冷空间内下落液滴冻结过程计算方法。该方法使用密度分布函数和温度分布函数分别模拟流体的速度场和温度场，通过求解固相体积分数追踪相变过程中各相的体积分布情况，使用固相体积分数修正焓值以及流场演化方程迁移步，并考虑相变潜热和温差相对较大(大于40K)时辐射换热对相变过程的影响。

为解决上述技术问题，本发明创造采用的技术方案是：

一种冷空间内下落液滴冻结过程的计算方法，主要包括如下步骤：

(步骤一)根据物理模型计算所需的计算域和网格尺寸，进行网格划分，如图1所示；

(步骤二)对模型计算域的密度、温度、速度以及固相体积分数进行初始化；

(步骤三)输出计算域的流场、温度场、辐射强度以及辐射热流量；

模型的流场演化方程为：

其中，f为粒子分布函数，τ为松弛时间。

平衡态分布函数(x,t)为：

其中，c_s为格子声速，ω_i为权系数，ρ为密度，e_i为离散速度。平衡态速度u^eq为：

分子间相互作用力为：

F(x)＝F_f(x)+F_g(x)

其中，流体间相互作用力为：

式中，ψ(ρ)＝ρ^*[1-exp(-ρ/ρ^*)]为有效密度。G₁和G₂分别为相邻和次相邻位置的流体间相互作用系数。

重力为：

F_g(x)＝ρ(x)g

演化方程迁移步为：

f_i(x,t+δ_t)＝α_sf_i(x,t)+(1-α_s)f_i(x-e_iδ_t,t)

其中，f_i为碰撞后的粒子分布函数。

温度场演化方程为：

其中，g_i为温度场的粒子分布函数。为辐射热流量，采用格子Boltzmann方法进行求解，表达式如下：

其中，I_i为辐射强度，为权系数。

计算辐射强度的演化方程为：

其中，τ_I为求解辐射强度分布函数的松弛时间。相应的平衡态分布函数为：

其中，为离散方向i对应的权系数。假设边界是黑体，边界辐射强度为σ_SB为Stefan-Boltzmann常数，T_bound为边界温度。

对应的温度场平衡态分布函数为：

(步骤四)输出计算域内每个节点的密度、速度以及温度；

密度、速度以及温度的计算式为：

(步骤五)输出计算域内每个节点的固相体积分数以及焓值；

固相体积分数α_s表达式为：

焓值表达式为：

H＝α_sc_solidT+(1-α_s)c_gasT+α_sL

(步骤六)返回(步骤三)循环计算，直到液滴内部固相体积分数全为1，程序结束。

本发明具有以下优点和有益效果：

本发明所述的一种冷空间内下落液滴冻结过程的计算方法，可以解决冷空间内悬浮液滴以及下落液滴冻结问题，得到不同时刻液滴内温度场、液滴冻结状态、冻结参数及影响因素，使得液滴冻结过程的研究更加完善，并得到如下结论：液滴下部空气冷却效果较强，液滴下表面优先发生冻结现象；垂直初速度对液滴形变有一定影响，随着下落速度的增大，形成的冰粒形状从圆形变为椭圆形；液滴冻结时间随液滴初始温度降低而缩短，随下落速度增大而缩短；上述结论可以用于指导实际生产，即提高冷空间内下落液滴冻结效率的最有效方法是改变其初始温度。

附图说明

图1体现了D2Q9模型的网格划分和速度矢量；

其中，图1(a)为D2Q9模型的网格划分图；图1(b)为D2Q9模型的速度矢量图。

图2为焓法格子Boltzmann相变模型与下落液滴冻结过程计算方法模拟得到的液滴冻结状态的对比图；

其中，图2(a)为焓法格子Boltzmann相变模型；图2(b)为下落液滴冻结过程计算方法。

图3为冷空间内下落液滴冻结过程的物理模型。

图4给出了不同网格体系中下落液滴内A点的固相体积分数。

图5体现了不同时刻冷空间的温度场和速度场；

其中，图5(a)t＝0s；图5(b)t＝0.6s；图5(c)t＝0.75s；图5(d)t＝0.9s；图5(e)t＝1.2s；

图5(f)t＝1.5s；图5(g)t＝1.8s；图5(h)t＝2.22s；图5(i)为温度云图图例。

图6体现了不同时刻冷空间内下落液滴的固相体积分数；

图6(a)t＝0s；图6(b)t＝0.6s；图6(c)t＝0.75s；图6(d)t＝0.9s；图6(e)t＝1.2s；图6(f)t＝1.5s；图6(g)t＝1.8s；图6(h)t＝2.22s；图6(i)为固相体积分数云图图例。

图7为t＝0.75s和t＝1.2s时下落液滴内O点到B点的固相体积分数对比图。

图8体现了相同时间间隔下落液滴内O点和B点的固相体积分数增量对比图。

图9为t＝0.9s时不同垂直初速度的下落液滴冻结形状对比图；

其中，图9(a)v₀＝0.001m/s；图9(b)v₀＝0.005m/s；图9(c)v₀＝0.01m/s；图9(d)v₀＝0.015m/s；图9(e)t＝0.02m/s。

图10为不同垂直初速度的下落液滴长短轴比对比图。

图11为不同液滴直径时液滴冻结时间与液滴初始温度的关系对比图。

具体实施方式

下面通过实施例对本发明作进一步详细说明。

实施例1

本实施例分别采用焓法格子Boltzmann模型和下落液滴冻结过程计算方法模拟冷空间内下落液滴冻结过程。选取计算域为0.03m×0.07m，计算网格数为600×1400，四周采用周期性边界。冷空间内有一初始温度为T_l,in＝273K的圆形液滴，液滴初始位置为(0.015m,0.06m)，液滴初始直径为D＝0.003m，液滴内O点为中心点，B点为边界点，A点为O点和B点的中点，液滴垂直初速度为v₀＝0.001m/s，水密度为ρ_l＝1000kg/m³，液滴初始固相体积分数为α_s＝0。冷空间内空气初始温度为T_g,in＝223K，空气密度为ρ_g＝1.584kg/m³。选取水的冻结温度为T_f＝268K，水比热容为c_liquid＝4.2×10³J/(kg·K)，冰比热容为c_solid＝2.01×10³J/(kg·K)，水导热率为k＝0.55W/(m·K)，空气Prandtl数为Pr＝0.71，Stefan数为Ste＝5.625。如图2所示，由对比结果可知，使用焓法格子Boltzmann相变模型模拟得到的下落液滴在表面形成冰壳后依然发生形变，这与实际下落液滴冻结过程不符，进一步说明了下落液滴冻结过程计算方法的正确性以及优势。

实施例2

利用本发明的冷空间内下落液滴冻结过程计算方法模拟冷空间内下落液滴冻结过程。图3给出了相应的物理模型，选取计算域为0.03m×0.07m，计算网格数为600×1400，四周采用周期性边界。冷空间内有一初始温度为T_l,in＝273K的圆形液滴，液滴初始位置为(0.015m,0.06m)，液滴初始直径为D＝0.003m，液滴内O点为中心点，B点为边界点，A点为O点和B点的中点，液滴垂直初速度为v₀＝0.001m/s，水密度为ρ_l＝1000kg/m³，液滴初始固相体积分数为α_s＝0。冷空间内空气初始温度为T_g,in＝223K，空气密度为ρ_g＝1.584kg/m³。选取水的冻结温度为T_f＝268K，水比热容为c_liquid＝4.2×10³J/(kg·K)，冰比热容为c_solid＝2.01×10³J/(kg·K)，水导热率为k＝0.55W/(m·K)，空气Prandtl数为Pr＝0.71，Stefan数为Ste＝5.625。

为了确保网格密度对计算精度没有影响，使用4套均匀网格分别进行数值模拟计算(网格1：150×350，网格2：300×700，网格3：600×1400，网格4：1200×2800)。图4给出了不同网格体系中下落液滴内A点的固相体积分数。模拟结果显示，网格3和网格4得到的模拟结果之间相对误差小于0.1％。考虑计算成本，采用网格3进行数值模拟计算。

图5和图6分别给出了不同时刻冷空间的温度场和速度场以及下落液滴的固相体积分数。冷空间内下落液滴冻结过程可以概括如下：液滴在下落过程中被周围环境冷却，液滴温度随时间不断降低，当液滴温度低于水的冻结温度时，液滴发生冻结现象。对于整个液滴来说，冻结现象首先发生在液滴表面，随着时间的推移，液滴内固液相界面位置逐渐向液滴中心移动，直到液滴完全冻结成冰粒。由模拟结果可知，液滴内部固相体积分数在下落过程后期呈不对称分布，液滴下部先于上部发生冻结现象，主要因为液滴下部空气冷却效果强于上部。

图7给出了t＝0.75s和t＝1.2s时下落液滴内O点(液滴中心)到B点(液滴边缘)的固相体积分数。模拟结果显示，同一时刻液滴表面的固相体积分数高于液滴中心，主要因为热量的传导由外向内进行，液滴表面先发生冻结现象。图8给出了相同时间间隔下落液滴内O点和B点的固相体积分数增量。初始时刻，液滴未开始冻结，液滴内O点和B点的固相体积分数增量为零。一段时间后，液滴开始冻结，此时液滴内O点和B点的固相体积分数增量先增大后减小。在液滴冻结过程后期，液滴冻结趋势逐渐减缓，主要因为相变驱动力在液滴冻结过程后期逐渐减小。

当下落液滴发生冻结现象时，液滴表面形成冰壳，冰壳的存在将阻碍下落液滴发生形变，随后液滴继续下落但不再发生形变。图9显示了t＝0.9s时不同垂直初速度的下落液滴冻结形状。模拟结果显示，随着液滴垂直初速度的增大，液滴形状从圆形变为椭圆形。图10给出了不同垂直初速度(从v₀＝0.001m/s到v₀＝0.02m/s)的下落液滴长短轴比。由模拟结果可知，液滴在垂直初速度较大时形变比较明显，因为垂直初速度越大，液滴下落速度越大，液滴下落过程中受到的阻力越大，液滴形变越明显。

冻结时间是决定冻结效率的重要因素之一，通过求解固相体积分数来预测液滴冻结时间。图11给出了不同液滴直径时液滴冻结时间与液滴初始温度的关系。液滴冻结时间随液滴初始温度降低而缩短，主要原因为相变驱动力随液滴初始温度降低而增大，较高的相变驱动力将促进相变现象发生，因此缩短了液滴冻结时间，提高了液滴冻结效率。

通过模拟冷空间内不同下落方式液滴冻结过程，获得了不同时刻冷空间的温度场、速度场以及下落液滴冻结状态及其相应的冻结参数，考察了下落方式和液滴初始温度对液滴冻结现象的影响，得到如下结论：液滴下部空气冷却效果较强，液滴下表面优先发生冻结现象；垂直初速度对液滴形变有一定影响，随着下落速度的增大(从v₀＝0.001m/s到v₀＝0.02m/s)，形成的冰粒形状从圆形变为椭圆形；液滴冻结时间随液滴初始温度降低而缩短，随下落速度增大而缩短。

上述结论可以用于指导实际生产，即可以通过延长冷却时间来降低液滴初始温度(从269K到278K)。此外，同一液滴初始温度时，减小液滴直径也可以缩短液滴冻结时间(以275K为例，将直径由3mm减小到1mm时，可以缩短2秒左右的冻结时间)，主要因为减小液滴直径将导致液滴单位重量表面积增大以及液滴终态速度降低，进而缩短液滴冻结时间。然而，在实际应用中(如人工造雪)，如果液滴直径太小，液滴会因终态速度过小而被风吹走。因此，提高冷空间内下落液滴冻结效率的最有效方法是改变其初始温度。

以上对本发明创造的一个实施例进行了详细说明，但所述内容仅为本发明创造的较佳实施例，不能被认为用于限定本发明创造的实施范围。凡依本发明创造申请范围所作的均等变化与改进等，均应仍归属于本发明创造的专利涵盖范围之内。

Claims (1)

1.一种冷空间内下落液滴冻结过程的计算方法，其特征在于，主要包括如下步骤：

步骤一、根据物理模型计算所需的计算域和网格尺寸，并进行网格划分；

步骤二、对模型计算域的密度、温度、速度以及固相体积分数进行初始化；

步骤三、输出计算域的流场、温度场、辐射强度以及辐射热流量；

模型的流场演化方程为：

其中，f为粒子分布函数，τ为松弛时间；

平衡态分布函数f_i ^(eq)(x,t)为：

其中，c_s为格子声速，ω_i为权系数，ρ为密度，e_i为离散速度；平衡态速度u^eq为：

分子间相互作用力为：

F(x)＝F_f(x)+F_g(x)

其中，流体间相互作用力为：

式中，ψ(ρ)＝ρ^*[1-exp(-ρ/ρ^*)]为有效密度；G₁和G₂分别为相邻和次相邻位置的流体间相互作用系数；

重力为：

F_g(x)＝ρ(x)g

演化方程迁移步为：

f_i(x,t+δ_t)＝α_sf_i(x,t)+(1-α_s)f_i(x-e_iδ_t,t)

其中，f_i为碰撞后的粒子分布函数；

温度场演化方程为：

其中，g_i为温度场的粒子分布函数；为辐射热流量，采用格子Boltzmann方法进行求解，表达式如下：

其中，I_i为辐射强度，为权系数；

计算辐射强度的演化方程为：

其中，τ_I为求解辐射强度分布函数的松弛时间；相应的平衡态分布函数为：

其中，为离散方向i对应的权系数；假设边界是黑体，边界辐射强度为σ_SB为Stefan-Boltzmann常数，T_bound为边界温度；

对应的温度场平衡态分布函数为：

步骤四、输出计算域内每个节点的密度、速度以及温度；

密度、速度以及温度的计算式为：

步骤五、输出计算域内每个节点的固相体积分数以及焓值；

固相体积分数α_s表达式为：

焓值表达式为：

H＝α_sc_solidT+(1-α_s)c_gasT+α_sL

步骤六、返回步骤三循环计算，直到液滴内部固相体积分数全为1，程序结束。