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

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

A simulation calculation method for frost suppression of precooler

【技术实现步骤摘要】
一种用于抑制预冷器结霜的模拟计算方法
本专利技术提供了预冷器结霜模拟运算
，尤其是一种用于抑制预冷器结霜的模拟计算方法。
技术介绍
结霜现象普遍存在于制冷、低温存储、空气调节、航空航天等领域。因此，掌握具体结霜过程以及合理解释结霜机理对于工程应用来说具有十分重要的研究意义。近年来，随着航空航天事业的不断进步与发展，对于结霜现象的研究又有了新的需求。航空航天设备表面与周围环境之间温差较大(大于100K)，多数情况设备会在大温差工况下运转并发生快速结霜现象。以预冷器为例，主流气体所含水蒸气的结霜现象是设备运行过程中一个较为严重的问题。结霜现象的存在会导致预冷器主流压力损失增大以及热交换量降低等危害，甚至引起传热失效，更有甚者还会在低空加速阶段造成灾难性事故。扩大预冷器管路间距虽然可以控制主流压力损失增大，但仍存在热交换性能下降以及预冷器尺寸增加等问题。过去有学者进行了大量关于热交换器结霜现象的研究，但大多以冷冻和空调用热交换器为研究对象，主要用来解决小温差长时间结霜问题，这与预冷器大温差短时间结霜现象的本质有所不同。在目前的研究成果中，关于预冷器结霜过程的研究十分有限。对于预冷器结霜现象的深入研究有利于其抑霜方法的探索，进而削弱结霜带给预冷器的危害。总结前人研究成果可知，虽然目前已有许多公式与试验数据较为符合，但是其结果精确度不高，应用范围具有一定局限性，无法表征预冷器的局部结霜特性；对于预冷器结霜试验而言，由于试验条件和技术手段的局限性，加之预冷器结霜过程发生的尺度较小，难以进行试验观测，更难以对其中参数进行测量，目前通过试验方法很难观测该结霜过程，因此需要借助数值模拟手段来辅助研究。鉴于结霜机理的复杂性以及结霜方式的多样性，对于预冷器结霜过程的研究而言，实际结霜过程的预测一直是研究中的难题。湿空气在大温差环境下的结霜方式不同于普通的小温差结霜过程，主要由湿空气壁面结霜和湿空气近冷壁面结霜两部分组成。现有的数学模型一般未考虑湿空气近冷壁面结霜，所有关于此类问题的数值模拟研究都是近似将该过程简化为湿空气壁面结霜来进行，这一简化使得预冷器结霜过程的模拟结果有可能产生较大误差。
技术实现思路
为了更好地实现预冷器结霜全过程，揭示湿空气近冷壁面结霜现象的本质及影响因素，完善湿空气结霜过程的研究，本专利技术提供了一种用于抑制预冷器结霜的模拟计算方法。本专利技术给出了一种用于抑制预冷器结霜的模拟计算方法，将改进的多组分焓法格子Boltzmann相变模型与LBM-LGA冰粒沉积模型耦合，使用密度分布函数和温度分布函数分别模拟流体的速度场和温度场，通过求解固相体积分数追踪相变过程中各相的体积分布情况，使用固相体积分数修正焓值以及流场演化方程迁移步，采用改进的多组分伪势多相流格子Boltzmann模型解决双组分两相流大密度比问题，并考虑相变潜热和辐射换热对相变过程的影响，利用LBM-LGA冰粒沉积模型模拟形成的冰粒沉积过程。为解决上述技术问题，本专利技术创造采用的技术方案是：一种用于抑制预冷器结霜的模拟计算方法，主要包括如下步骤：(步骤一)根据物理模型计算所需的计算域和网格尺寸，并进行网格划分，如图1所示；(步骤二)对模型计算域的密度、温度、速度以及固相体积分数进行初始化；(步骤三)输出计算域的流场、温度场、辐射强度以及辐射热流量；流场演化方程为：其中，为组分σ的粒子分布函数，为组分σ的平衡态分布函数，τσ为组分σ的松弛时间；平衡态分布函数为：其中，cs为格子声速，wα为权系数，ρ为密度，eα为离散速度；平衡态速度uσeq(r,t)为：分子间相互作用力Fσ(r,t)为：其中，流体间相互作用力为：式中，ψσ(r,t)＝ρ0[1-exp(-ρσ(r,t))/ρ0]为组分σ的有效密度；G1和G2分别为相邻和次相邻位置的流体间相互作用系数；重力为：双组分混合流体速度u(r,t)为：其中，Mσ为组分σ的分子质量；演化方程迁移步为：其中，为组分σ碰撞后的粒子分布函数；温度场演化方程为：其中，nα(r,t)为温度场的粒子分布函数；为辐射热流量，利用格子Boltzmann方法进行求解，表达式如下：其中，Iα为辐射强度，为权系数；计算辐射强度的演化方程为：其中，τI为求解辐射强度分布函数的松弛时间；相应的平衡态分布函数为：其中，为离散方向α对应的权系数；假设边界是黑体，边界辐射强度为σSB为Stefan-Boltzmann常数，Tbound为边界温度；对应的温度场平衡态分布函数为：(步骤四)输出计算域内每个节点的密度、速度以及温度；密度、速度以及温度的计算式为：(步骤五)输出计算域内每个节点的固相体积分数以及焓值；固相体积分数αs表达式为：焓值表达式为：H＝αscsolidT+(1-αs)cgasT+αsL(步骤六)输出冰粒的位移和速度；当湿空气中水蒸气形成冰粒后，下一时刻，冰粒以一定概率移动到相邻格点或停留在原来位置；假设t时刻，网格节点rp上有一冰粒；t+δt时刻，冰粒的最终格点位置为：r'p＝rp+(v1e1+v2e2+v3e3+v4e4)δt其中，vα为随机布尔变量，取值为1的概率为pα；在D2Q9模型中，t+δt时刻，冰粒移动到东、北、西、南(即p1、p3、p5和p7)的输运概率pα为：其中，δrp为时间步长δt时冰粒的实际位移，考虑曳力以及重力和浮力的冰粒运动方程为：其中，Fd为曳力，Fg为重力和浮力；重力和浮力的表达式为：Fg＝(mp-mg)g其中，为冰粒质量，ρi为冰密度，ri为冰粒半径，为冰粒受到的浮力，ρg为气体密度，g为重力加速度；在粒子跟踪方法中，通过冰粒运动方程求解的冰粒速度和位移为：其中，τp为冰粒的松弛时间，具体数值由Stokes数确定；当冰粒碰撞冷壁面的法向速度小于临界沉积速度Vcr时，冰粒不再运动并开始在冷壁面上沉积，最终成为冷壁面上霜层的一部分；基于Brach和Dunn的工作，临界沉积速度Vcr可以表示为：Vcr＝[2K/(DiR2)]10/7其中，Di为冰粒直径；为有效刚度系数；相关参数Es和Ep分别为冷壁面和冰粒的杨氏模量；νs和νp分别为冷壁面和冰粒的泊松比；R为运动恢复系数，常取0.9；(步骤七)输出霜层厚度；霜层平均厚度为：其中，δf为霜层的平均厚度，At为管束总面积，为霜层总体积，N为计算霜层的总网格数；(步骤八)返回(步骤三)循环计算，直到霜层厚度达到指定值，程序结束。本专利技术具有以下优点和有益效果：本专利技术的本文档来自技高网...

【技术保护点】
1.一种用于抑制预冷器结霜的模拟计算方法，其特征在于，主要包括如下步骤：/n(步骤一)根据物理模型计算所需的计算域和网格尺寸，并进行网格划分；/n(步骤二)对模型计算域的密度、温度、速度以及固相体积分数进行初始化；/n(步骤三)输出计算域的流场、温度场、辐射强度以及辐射热流量；/n流场演化方程为：/n

【技术特征摘要】
1.一种用于抑制预冷器结霜的模拟计算方法，其特征在于，主要包括如下步骤：
(步骤一)根据物理模型计算所需的计算域和网格尺寸，并进行网格划分；
(步骤二)对模型计算域的密度、温度、速度以及固相体积分数进行初始化；
(步骤三)输出计算域的流场、温度场、辐射强度以及辐射热流量；
流场演化方程为：



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



其中，cs为格子声速，wα为权系数，ρ为密度，eα为离散速度；平衡态速度uσeq(r,t)为：



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



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



式中，ψσ(r,t)＝ρ0[1-exp(-ρσ(r,t))/ρ0]为组分σ的有效密度；G1和G2分别为相邻和次相邻位置的流体间相互作用系数；
重力为：



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



其中，Mσ为组分σ的分子质量；
演化方程迁移步为：



其中，为组分σ碰撞后的粒子分布函数；
温度场演化方程为：



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



其中，Iα为辐射强度，为权系数；
计算辐射强度的演化方程为：



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



其中，为离散方向α对应的权系数；假设边界是黑体，边界辐射强度为σSB为Stefan-Boltzmann常数，Tbound为边界温度；
对应的温度场平衡态分布函数为：



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

【专利技术属性】
技术研发人员：赵鑫，陈光，王晓兵，杨沄芃，郭帅帅，郝冬，张妍懿，王仁广，
申请(专利权)人：中汽研汽车检验中心天津有限公司，中国汽车技术研究中心有限公司，
类型：发明
国别省市：天津;12