[发明专利]一种用于抑制预冷器结霜的模拟计算方法

权利要求书

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为计算霜层的总网格数；

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