CN111159890B - 一种用于抑制预冷器结霜的模拟计算方法 - Google Patents

Abstract

本发明公开了一种用于抑制预冷器结霜的模拟计算方法，可以实现对预冷器结霜全过程的运算，揭示湿空气近冷壁面结霜现象的本质及影响因素，使得湿空气结霜过程的研究更加完善；通过计算可得，在一定时间范围内，增大管间距，减小管径，降低湿空气进口流速，提高管壁壁温或者降低湿空气进口相对湿度可以有效抑制预冷器极低温管束表面霜层的形成，起到抑霜作用。

Description

一种用于抑制预冷器结霜的模拟计算方法

技术领域

本发明提供了预冷器结霜模拟运算技术领域，尤其是一种用于抑制预冷器结霜的模拟计算方法。

背景技术

结霜现象普遍存在于制冷、低温存储、空气调节、航空航天等领域。因此，掌握具体结霜过程以及合理解释结霜机理对于工程应用来说具有十分重要的研究意义。近年来，随着航空航天事业的不断进步与发展，对于结霜现象的研究又有了新的需求。航空航天设备表面与周围环境之间温差较大(大于100K)，多数情况设备会在大温差工况下运转并发生快速结霜现象。以预冷器为例，主流气体所含水蒸气的结霜现象是设备运行过程中一个较为严重的问题。结霜现象的存在会导致预冷器主流压力损失增大以及热交换量降低等危害，甚至引起传热失效，更有甚者还会在低空加速阶段造成灾难性事故。扩大预冷器管路间距虽然可以控制主流压力损失增大，但仍存在热交换性能下降以及预冷器尺寸增加等问题。过去有学者进行了大量关于热交换器结霜现象的研究，但大多以冷冻和空调用热交换器为研究对象，主要用来解决小温差长时间结霜问题，这与预冷器大温差短时间结霜现象的本质有所不同。在目前的研究成果中，关于预冷器结霜过程的研究十分有限。对于预冷器结霜现象的深入研究有利于其抑霜方法的探索，进而削弱结霜带给预冷器的危害。总结前人研究成果可知，虽然目前已有许多公式与试验数据较为符合，但是其结果精确度不高，应用范围具有一定局限性，无法表征预冷器的局部结霜特性；对于预冷器结霜试验而言，由于试验条件和技术手段的局限性，加之预冷器结霜过程发生的尺度较小，难以进行试验观测，更难以对其中参数进行测量，目前通过试验方法很难观测该结霜过程，因此需要借助数值模拟手段来辅助研究。鉴于结霜机理的复杂性以及结霜方式的多样性，对于预冷器结霜过程的研究而言，实际结霜过程的预测一直是研究中的难题。湿空气在大温差环境下的结霜方式不同于普通的小温差结霜过程，主要由湿空气壁面结霜和湿空气近冷壁面结霜两部分组成。现有的数学模型一般未考虑湿空气近冷壁面结霜，所有关于此类问题的数值模拟研究都是近似将该过程简化为湿空气壁面结霜来进行，这一简化使得预冷器结霜过程的模拟结果有可能产生较大误差。

发明内容

为了更好地实现预冷器结霜全过程，揭示湿空气近冷壁面结霜现象的本质及影响因素，完善湿空气结霜过程的研究，本发明提供了一种用于抑制预冷器结霜的模拟计算方法。

本发明给出了一种用于抑制预冷器结霜的模拟计算方法，将改进的多组分焓法格子Boltzmann相变模型与LBM-LGA冰粒沉积模型耦合，使用密度分布函数和温度分布函数分别模拟流体的速度场和温度场，通过求解固相体积分数追踪相变过程中各相的体积分布情况，使用固相体积分数修正焓值以及流场演化方程迁移步，采用改进的多组分伪势多相流格子Boltzmann模型解决双组分两相流大密度比问题，并考虑相变潜热和辐射换热对相变过程的影响，利用LBM-LGA冰粒沉积模型模拟形成的冰粒沉积过程。

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

一种用于抑制预冷器结霜的模拟计算方法，主要包括如下步骤：

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

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

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

流场演化方程为：

其中，为组分σ的粒子分布函数，/>为组分σ的平衡态分布函数，τ_σ为组分σ的松弛时间；

平衡态分布函数为：

其中，c_s为格子声速，w_α为权系数，ρ为密度，e_α为离散速度；平衡态速度u_σ ^eq(r,t)为：

分子间相互作用力F^σ(r,t)为：

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

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

重力为：

双组分混合流体速度u(r,t)为：

其中，M^σ为组分σ的分子质量；

演化方程迁移步为：

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

温度场演化方程为：

其中，n_α(r,t)为温度场的粒子分布函数；为辐射热流量，利用格子Boltzmann方法进行求解，表达式如下：

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

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

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

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

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

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

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

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

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

焓值表达式为：

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

(步骤六)输出冰粒的位移和速度；

当湿空气中水蒸气形成冰粒后，下一时刻，冰粒以一定概率移动到相邻格点或停留在原来位置；假设t时刻，网格节点r_p上有一冰粒；t+δ_t时刻，冰粒的最终格点位置为：

r'_p＝r_p+(v₁e₁+v₂e₂+v₃e₃+v₄e₄)δ_t

其中，v_α为随机布尔变量，取值为1的概率为p_α；

在D2Q9模型中，t+δ_t时刻，冰粒移动到东、北、西、南(即p₁、p₃、p₅和p₇)的输运概率p_α为：

其中，δ_rp为时间步长δ_t时冰粒的实际位移，

考虑曳力以及重力和浮力的冰粒运动方程为：

其中，F_d为曳力，F_g为重力和浮力；

重力和浮力的表达式为：

F_g＝(m_p-m_g)g

其中，为冰粒质量，ρ_i为冰密度，r_i为冰粒半径，/>为冰粒受到的浮力，ρ_g为气体密度，g为重力加速度；

在粒子跟踪方法中，通过冰粒运动方程求解的冰粒速度和位移为：

其中，τ_p为冰粒的松弛时间，具体数值由Stokes数确定；

当冰粒碰撞冷壁面的法向速度小于临界沉积速度V_cr时，冰粒不再运动并开始在冷壁面上沉积，最终成为冷壁面上霜层的一部分；基于Brach和Dunn的工作，临界沉积速度V_cr可以表示为：

V_cr＝[2K/(D_iR²)]^10/7

其中，D_i为冰粒直径；为有效刚度系数；相关参数E_s和E_p分别为冷壁面和冰粒的杨氏模量；ν_s和ν_p分别为冷壁面和冰粒的泊松比；R为运动恢复系数，常取0.9；

(步骤七)输出霜层厚度；

霜层平均厚度为：

其中，δ_f为霜层的平均厚度，A_t为管束总面积，为霜层总体积，N为计算霜层的总网格数；

(步骤八)返回(步骤三)循环计算，直到霜层厚度达到指定值，程序结束。

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

本发明的计算方法，可以实现对预冷器结霜全过程的运算，揭示湿空气近冷壁面结霜现象的本质及影响因素，使得湿空气结霜过程的研究更加完善；通过计算可得，在一定时间范围内，增大管间距，减小管径，降低湿空气进口流速，提高管壁壁温或者降低湿空气进口相对湿度可以有效抑制预冷器极低温管束表面霜层的形成，起到抑霜作用。

附图说明

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

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

图2为预冷器极低温管束通道结霜过程的物理模型。

图3为t＝1s和t＝2.5s时不同网格体系的管束平均霜层厚度图。

图4体现了不同时刻预冷器极低温管束通道的结霜工况；

其中，图4(a)为t＝0s时管束结霜工况图，图4(b)为t＝1s时管束结霜工况图，图4(c)为t＝2s时管束结霜工况图，图4(d)为t＝2.5s时管束结霜工况图，图4(e)为t＝4.5s时管束结霜工况图，图4(f)为t＝2.5s时圆管局部放大图，图4(g)为固相体积分数云图图例。

图5体现了t＝0s和t＝2.5s时预冷器极低温管束通道的温度场和压力场；

其中，图5(a)为t＝0s时的温度场，图5(b)为t＝2.5s时的温度场，图5(c)为t＝0s时的压力场，图5(d)为t＝2.5s时的压力场，图5(e)为温度云图图例，图5(f)为压力云图图例。

图6体现了管间距对管束平均霜层厚度、湿空气侧平均换热系数和湿空气侧压力损失系数的影响(模拟工况：湿空气进口速度v_in＝10m/s，管壁温度T_t＝150K，管径D_t＝0.002m，湿空气进口相对湿度)；

其中，图6(a)为管间距对管束平均霜层厚度影响图；图6(b)为管间距对湿空气侧平均换热系数影响图；图6(c)为管间距对湿空气侧压力损失系数影响图。

图7体现了管径对管束平均霜层厚度、湿空气侧平均换热系数和湿空气侧压力损失系数的影响(模拟工况：管间距s_t＝0.003m，湿空气进口速度v_in＝10m/s，管壁温度T_t＝150K，湿空气进口相对湿度)；

其中，图7(a)为管径对管束平均霜层厚度的影响图；图7(b)为管径对湿空气侧平均换热系数的影响图；图7(c)为管径对湿空气侧压力损失系数的影响图。

图8体现了湿空气进口流速对管束平均霜层厚度、湿空气侧平均换热系数和湿空气侧压力损失系数的影响(模拟工况：管间距s_t＝0.003m，管壁温度T_t＝150K，管径D_t＝0.002m，湿空气进口相对湿度)；

其中，图8(a)为湿空气进口流速对管束平均霜层厚度影响图；图8(b)为湿空气进口流速对湿空气侧平均换热系数影响图；图8(c)为湿空气进口流速对湿空气侧压力损失系数影响图。

图9体现了管壁温度对管束平均霜层厚度、湿空气侧平均换热系数和湿空气侧压力损失系数的影响(模拟工况：管间距s_t＝0.003m，湿空气进口速度v_in＝10m/s，管径D_t＝0.002m，湿空气进口相对湿度)；

其中，图9(a)为管壁温度对管束平均霜层厚度的影响图；图9(b)为管壁温度对湿空气侧平均换热系数的影响图；图9(c)为管壁温度对湿空气侧压力损失系数的影响图。

图10体现了湿空气进口相对湿度对管束平均霜层厚度、湿空气侧平均换热系数和湿空气侧压力损失系数的影响(模拟工况：管间距s_t＝0.003m，湿空气进口速度v_in＝10m/s，管壁温度T_t＝150K，管径D_t＝0.002m)；

其中，图10(a)为湿空气进口相对湿度对管束平均霜层厚度的影响图；图10(b)为湿空气进口相对湿度对湿空气侧平均换热系数的影响图；图10(c)为湿空气进口相对湿度对湿空气侧压力损失系数的影响图。

具体实施方式

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

本实施例使用预冷器结霜计算方法模拟预冷器极低温管束通道的结霜过程。相应的物理模型如图2所示，计算域长和高分别为0.04m和0.006m，计算网格数为1600×240。计算域内有十列二十五根直径为D_t＝0.002m的呈等边三角形交叉排列的圆管，管间距为s_t＝0.003m，左右边界分别为速度边界和流出边界，上下边界为对称边界，圆管管壁为无滑移壁面边界。湿空气初始温度为T_v,i＝295K，湿空气进口相对湿度为75％，湿空气中空气与水蒸气密度比为ρ_a/ρ_w＝80，湿空气进口流速为v_in＝10m/s。管壁温度为T_t＝150K且保持不变。选取水蒸气冻结温度为T_f＝268K，水蒸气比热容为c_gas＝1.87×10³J/(kg·K)，冰比热容为c_solid＝2.01×10³J/(kg·K)，湿空气中水蒸气与空气运动粘度比为ν_w/ν_a＝33，湿空气Prandtl数为Pr＝0.705，Reynolds数为Re＝66，Stefan数为Ste＝21，Stokes数为0.01。初始固相体积分数α_s＝0。

湿空气侧平均换热系数为：

其中，Q为计算域入口到出口的热流量，ΔT为湿空气流与管束壁面的温差。

湿空气侧压力损失系数为：

其中，Δp为计算域入口到出口的压降，p_d为入口压力。

为了确保网格密度对计算精度没有影响，需要对模型进行网格无关性验证。本节分别模拟了4套网格体系下预冷器极低温管束通道的结霜过程(网格1：60×400，网格2：120×800，网格3：240×1600，网格4：480×3200)。图3给出了t＝1s和t＝2.5s时不同网格体系的管束平均霜层厚度。可以看出，网格3与网格4的模拟结果之间相对误差不超过0.1％。考虑计算成本，本节采用网格3进行数值模拟计算。

图4给出了不同时刻预冷器极低温管束通道的结霜工况。湿空气进入管束通道后被低温管束冷却，当湿空气温度被冷却到水蒸气冻结温度以下时，湿空气中的水蒸气会发生结霜现象。由于湿空气与周围环境之间温差较大，湿空气中的水蒸气不仅在低温管束壁面上发生结霜现象，在距离低温管束壁面或霜层表面0.15mm处也有可能形成细小冰粒，随后这些细小冰粒在低温管束壁面或霜层表面发生沉积。模拟结果显示，圆管管壁上的霜层厚度随时间不断增加(如图4(a)-(e)所示)。由图4(f)可知，单根圆管管壁上的霜层分布不均，圆管迎风侧结霜量略小于背风侧，主要因为大温差环境湿空气结霜过程中湿空气壁面结霜占主导，而在湿空气壁面结霜过程中圆管背风侧结霜量受强制对流换热效应影响大于迎风侧结霜量。此外，如图4(e)所示，后列圆管管壁上的霜层厚度大于前列圆管，因为后列圆管处的湿空气温度相对较低，湿空气中的水蒸气更容易发生结霜现象，因此结霜量相对较大。同时，部分湿空气中水蒸气形成的小冰粒也会随湿空气流运动到后列圆管，撞击圆管管壁发生沉积。如图4(e)所示，当霜层增加到一定厚度时，管路通道被形成的霜层堵塞，进而发生传热失效，由于管束结构主要用作换热部件，因此仅分析有效传热时间内的模拟结果。图5给出了t＝0s和t＝2.5s时预冷器极低温管束通道的温度场和压力场。可见，通道出口湿空气温度随时间逐渐降低，入口与出口压差随时间逐渐增大。

为了更好地分析结构参数(管间距和管径)以及热力学参数(湿空气进口流速、管壁温度和湿空气进口相对湿度)对预冷器极低温管束结霜工况和换热特性的影响，对预冷器极低温管束结霜过程进行变参数研究，如表1所示，具体工况总结如下：

表1预冷器极低温管束结霜过程的模拟工况

(1)管间距

定义管束密实度为湿空气流经管束总面积与换热器总体积的比值，表征换热器的结构特性。如表2所示，相同尺寸通道内管间距越小，管束换热面积越大。图6给出了管间距对管束平均霜层厚度、湿空气侧平均换热系数和湿空气侧压力损失系数的影响。相同工况下，0～0.75s内，减小管束管间距将增加管束霜层厚度。因为当管束管间距减小时，管束换热面积增加，相应的换热量增大，进而导致霜层厚度增加。当霜层厚度增加到一定厚度(管壁间距的1/4左右)时，流道被形成的霜层堵塞，此时传热失效，霜层堵塞流道的时间随管间距减小而缩短。如图6(b)所示，在0～0.75s范围内，对于管间距较小的工况，湿空气侧平均换热系数有所降低。当管束管间距减小时，管束霜层厚度增加，霜层的存在增加了热阻，因此湿空气侧平均换热系数降低。同时，在0～0.75s内，管束管间距越小，湿空气侧压力损失系数越大。因为减小管束管间距会导致管束霜层厚度增加，进而引起湿空气侧压力损失增大。因此，结霜工况下不建议使用过小的管间距。

表2相同尺寸通道内不同管间距的管束密实度

(2)管径

表3给出了相同尺寸通道内不同管径的管束密实度。相同尺寸通道内管径越大，管束换热面积越大。图7给出了管径对管束平均霜层厚度、湿空气侧平均换热系数和湿空气侧压力损失系数的影响。由模拟结果可知，对于相同尺寸通道而言，0～0.7s内，管径越大，管束霜层厚度越厚。因为管径越大，管束换热面积越大，相应的换热量越大，形成的霜层厚度越厚。当霜层厚度增加到一定厚度(管壁间距的1/4左右)时，流道被形成的霜层堵塞，霜层堵塞流道的时间随管径增大而缩短。此外，0～0.7s范围内，管径越大，湿空气侧平均换热系数越小。当管径增大时，霜层厚度的增加将导致热阻增大，进而引起换热性能下降。同时，湿空气侧压力损失系数也随管径增大而增大。因此，在选取换热器管束时，应选择较细管束。

表3相同尺寸通道内不同管径的管束密实度

(3)湿空气进口流速

图8给出了湿空气进口流速对管束平均霜层厚度、湿空气侧平均换热系数和湿空气侧压力损失系数的影响。0～1.7s内，当湿空气进口流速增大时，同一时刻管束霜层厚度明显增加。这是因为湿空气进口流速越大，其强制对流换热效应越强，进而导致换热量增加，霜层厚度增大。当霜层厚度增加到一定厚度(管壁间距的1/4左右)时，流道被形成的霜层堵塞，此时传热失效，霜层堵塞流道的时间随湿空气进口流速增大而缩短。此外，0～1.7s内，湿空气进口流速越大，湿空气侧平均换热系数越大。因为湿空气进口流速的增大将导致强制对流换热效应的增强。如图8(c)所示，湿空气侧压力损失系数随湿空气进口流速增大而增加，这是由湿空气侧压力损失系数随霜层厚度增加而增大所导致的。

(4)管壁温度

图9给出了管壁温度对管束平均霜层厚度、湿空气侧平均换热系数和湿空气侧压力损失系数的影响。1s内，管壁温度越低，管束霜层厚度越厚。管壁温度的降低有利于霜层的形成，当管壁温度降低时，湿空气流与管壁温差增大，相变驱动力增强，进而导致结霜量增加，霜层厚度增大。当霜层厚度增加到一定厚度(管壁间距的1/4左右)时，流道被形成的霜层堵塞，霜层堵塞流道的时间随管壁温度降低而缩短。1s范围内，管壁温度越低，湿空气侧平均换热系数越低。结霜工况下，一方面，降低管壁温度可以强化换热；另一方面，较低的管壁温度将导致更大的结霜量，相应的热阻增加，湿空气侧平均换热系数降低。模拟结果表明，在管壁温度对湿空气流经极低温管束换热性能的影响方面，后者占主导地位。此外，随着管壁温度的降低，湿空气侧压力损失系数也相应地有所增大。因为管壁温度越低，霜层厚度越厚，湿空气侧压力损失越大。

(5)湿空气进口相对湿度

图10给出了湿空气进口相对湿度对管束平均霜层厚度、湿空气侧平均换热系数和湿空气侧压力损失系数的影响。0～1.6s范围内，湿空气进口相对湿度越大，管束霜层厚度越厚。当湿空气进口相对湿度增大时，湿空气中的水蒸气含量升高，形成的霜层厚度增加。当霜层厚度增加到一定厚度(管壁间距的1/4左右)时，流道被形成的霜层堵塞，霜层堵塞流道的时间随湿空气进口相对湿度增大而缩短。湿空气进口相对湿度越大，湿空气侧平均换热系数越低。1.6s内，当湿空气进口相对湿度增大时，霜层厚度增加，相应的热阻增大，进而导致湿空气侧平均换热系数降低。此外，湿空气侧压力损失系数随湿空气进口相对湿度增加而增大。

根据上述模拟结果分析可以得出：相同工况下，0～0.75s内增大管束管间距(从0.025m增大到0.004m)，0～0.7s内减小管径(从0.0025m减小到0.001m)，0～1.7s内降低湿空气进口流速(从12m/s降低到8m/s)，0～1s内提高管壁壁温(从130K提高到150K)或者0～1.6s内降低湿空气进口相对湿度(从85％降低到65％)可以有效抑制预冷器极低温管束表面霜层的形成，起到抑霜的作用。

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

Claims (1)

1.一种用于抑制预冷器结霜的模拟计算方法，其特征在于，主要包括如下步骤：

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

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

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

流场演化方程为：

其中，为组分σ的粒子分布函数，/>为组分σ的平衡态分布函数，τ_σ为组分σ的松弛时间；

平衡态分布函数为：

其中，c_s为格子声速，w_α为权系数，ρ为密度，e_α为离散速度；平衡态速度u_σ ^eq(r,t)为：

分子间相互作用力F^σ(r,t)为：

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

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

重力为：

双组分混合流体速度u(r,t)为：

其中，M^σ为组分σ的分子质量；

演化方程迁移步为：

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

温度场演化方程为：

其中，n_α(r,t)为温度场的粒子分布函数；为辐射热流量，利用格子Boltzmann方法进行求解，表达式如下：

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

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

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

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

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

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

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

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

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

焓值表达式为：

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

步骤六、输出冰粒的位移和速度；

当湿空气中水蒸气形成冰粒后，下一时刻，冰粒以一定概率移动到相邻格点或停留在原来位置；假设t时刻，网格节点r_p上有一冰粒；t+δ_t时刻，冰粒的最终格点位置为：

r'_p＝r_p+(v₁e₁+v₂e₂+v₃e₃+v₄e₄)δ_t

其中，v_α为随机布尔变量，取值为1的概率为p_α；

在D2Q9模型中，t+δ_t时刻，冰粒移动到东、北、西、南即p₁、p₃、p₅和p₇的输运概率p_α为：

α＝1,3,5,7

其中，为时间步长δ_t时冰粒的实际位移，/>

考虑曳力以及重力和浮力的冰粒运动方程为：

其中，F_d为曳力，F_g为重力和浮力；

重力和浮力的表达式为：

F_g＝(m_p-m_g)g

其中，为冰粒质量，ρ_i为冰密度，r_i为冰粒半径，/>为冰粒受到的浮力，ρ_g为气体密度，g为重力加速度；

在粒子跟踪方法中，通过冰粒运动方程求解的冰粒速度和位移为：

其中，τ_p为冰粒的松弛时间，具体数值由Stokes数确定；

当冰粒碰撞冷壁面的法向速度小于临界沉积速度V_cr时，冰粒不再运动并开始在冷壁面上沉积，最终成为冷壁面上霜层的一部分；基于Brach和Dunn的工作，临界沉积速度V_cr可以表示为：

V_cr＝[2K/(D_iR²)]^10/7

其中，D_i为冰粒直径；为有效刚度系数；相关参数E_s和E_p分别为冷壁面和冰粒的杨氏模量；ν_s和ν_p分别为冷壁面和冰粒的泊松比；R为运动恢复系数，常取0.9；

步骤七、输出霜层厚度；

霜层平均厚度为：

其中，δ_f为霜层的平均厚度，A_t为管束总面积，为霜层总体积，N为计算霜层的总网格数；

步骤八、返回步骤三循环计算，直到霜层厚度达到指定值，程序结束。