CN115114815A - 一种利用霜层表面性质预测冷表面结霜的模拟方法 - Google Patents
一种利用霜层表面性质预测冷表面结霜的模拟方法 Download PDFInfo
- Publication number
- CN115114815A CN115114815A CN202210569667.8A CN202210569667A CN115114815A CN 115114815 A CN115114815 A CN 115114815A CN 202210569667 A CN202210569667 A CN 202210569667A CN 115114815 A CN115114815 A CN 115114815A
- Authority
- CN
- China
- Prior art keywords
- frost layer
- frost
- grid
- density
- water vapor
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 34
- 238000004088 simulation Methods 0.000 title claims abstract description 31
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Chemical compound O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims abstract description 67
- 238000012546 transfer Methods 0.000 claims abstract description 24
- 230000008569 process Effects 0.000 claims abstract description 10
- 238000004458 analytical method Methods 0.000 claims abstract description 5
- 238000004364 calculation method Methods 0.000 claims description 31
- 230000008859 change Effects 0.000 claims description 25
- 229920006395 saturated elastomer Polymers 0.000 claims description 16
- 239000012530 fluid Substances 0.000 claims description 9
- 230000006870 function Effects 0.000 claims description 7
- 230000008093 supporting effect Effects 0.000 claims description 6
- 238000012795 verification Methods 0.000 claims description 6
- 238000010586 diagram Methods 0.000 claims description 4
- 238000005516 engineering process Methods 0.000 claims description 4
- 230000008878 coupling Effects 0.000 claims description 3
- 238000010168 coupling process Methods 0.000 claims description 3
- 238000005859 coupling reaction Methods 0.000 claims description 3
- 239000006071 cream Substances 0.000 claims description 3
- 238000009792 diffusion process Methods 0.000 claims description 3
- 230000000704 physical effect Effects 0.000 claims description 3
- 238000007711 solidification Methods 0.000 claims description 3
- 230000008023 solidification Effects 0.000 claims description 3
- 230000007704 transition Effects 0.000 claims description 3
- 238000013523 data management Methods 0.000 claims description 2
- 239000013078 crystal Substances 0.000 description 6
- 238000005057 refrigeration Methods 0.000 description 4
- 230000015572 biosynthetic process Effects 0.000 description 2
- 238000013461 design Methods 0.000 description 2
- 238000005265 energy consumption Methods 0.000 description 2
- 230000007246 mechanism Effects 0.000 description 2
- 230000000149 penetrating effect Effects 0.000 description 2
- 239000011148 porous material Substances 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000010257 thawing Methods 0.000 description 2
- YCKRFDGAMUMZLT-UHFFFAOYSA-N Fluorine atom Chemical compound [F] YCKRFDGAMUMZLT-UHFFFAOYSA-N 0.000 description 1
- 238000003556 assay Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000001413 cellular effect Effects 0.000 description 1
- 229910052731 fluorine Inorganic materials 0.000 description 1
- 239000011737 fluorine Substances 0.000 description 1
- 238000009413 insulation Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000006911 nucleation Effects 0.000 description 1
- 238000010899 nucleation Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 239000000758 substrate Substances 0.000 description 1
- 239000013589 supplement Substances 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/08—Thermal analysis or thermal optimisation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Computing Systems (AREA)
- Fluid Mechanics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明属于结霜模拟运算技术领域，尤其涉及一种利用霜层表面性质预测冷表面结霜的模拟方法，包括以下步骤：S1、建立冷表面物理模型；S2、获得满足模拟精度的网格数量；S3、计算流场参数；S4、标记并获得水蒸气与冰的传质速率；S5、获得霜层的厚度及密度；S6、构建用于冷平板的霜层生长数值分析方法；S7、同工况下与相关文献实验结果进行对比，验证模型的准确性；S8、将构建霜层生长模型用于实际，分析冷表面上霜层生长的分布情况。本发明通过对前人忽视了的霜层表面性质，尤其是冰项体积分数的研究，具体对霜层厚度的增大过程进行描述，提出了更细致考虑更全面的霜层生长数值模型。计算的霜层生长情况与实际的霜层厚度、密度更为接近。
Description
技术领域
本发明属于结霜模拟运算技术领域，尤其涉及一种利用霜层表面性质预测冷表面结霜的模拟方法。
背景技术
结霜是制冷、低温和航空航天领域中非常常见的物理现象。当温度较高的湿空气遇到温度较低的冷表面时，如果冷表面的温度低于水蒸气对应分压力下的露点温度且低于零度，那么就会发生凝华现象，即结霜。研究发现，少量的霜能提高换热器的换热性能，但是随着霜层逐渐增厚，霜层内部会充斥着孔隙，孔隙中为热导率极低的空气，从而使整个霜层具有显著的隔热性。由于霜的多孔性，其有效热导率仅能达到冰的十分之一。过厚的霜堵塞流通通道，使换热器的总传热系数大大降低，严重时甚至会使换热器停机。
据报道，结霜工况下空气源热泵换热能力下降可达35％以上，全世界50％～ 75％的热泵换热器性能降低由结霜引起，结霜带来的能耗损失约占热泵总能耗的10.2％。此外，除、化霜也同样需要消耗大量能源，且会给正在使用制冷设备的人们体感上的不适。所以，为了抑制制冷设备的结霜，首先有必要对冷表面上的霜层生长进行合理预测。而现有的预测方法多集中于制冷***研究，且主要集中于调整物理模型的结构，对于结霜现象的机理描述不清。
发明内容
本发明提供一种利用霜层表面性质预测冷表面结霜的模拟方法，目的在于对霜层厚度的增大过程进行描述，提出了更细致考虑更全面的霜层生长数值模型，计算的霜层生长情况与实际的霜层厚度、密度更为接近。
为实现上述目的，本发明采用的技术方案为：
一种利用霜层表面性质预测冷表面结霜的模拟方法，其特征在于：包括以下步骤：
S1、建立冷表面物理模型；
利用ICEM软件建立冷表面物理模型，采用四边形结构化网格对计算域网格进行划分，并对近壁面网格在满足Y+的要求下进行加密处理；
S2、获得满足模拟精度的网格数量；
对划分的计算域网格进行网格无关性验证，通过自适应网格技术调整网格，同时分别划分几种不同数量的网格进行对比验证，随着网格数量的增加，霜层的总质量不在变化且差值逐渐趋于稳定时，即可得到满足模拟精度的网格数量；
S3、计算流场参数；
利用连续性、动量、能量和组分方程对整个流场区域进行模拟计算，各控制方程通用形式如下：
其中，ρ为流体密度；Φ为广义变量；Γ为广义扩散系数；S为广义源相；
S4、标记并获得水蒸气与冰的传质速率；
通过用户自定义UDM把能否相变的区域分别储存标记为1和0；由于靠近壁面的第一层网格始终存在着冷表面的支撑作用，在整个计算时间内始终被标记为1；而那些处于同一列，且正下方接触着冰项体积分数大于临界体积分数的网格，也被标记为1；在t＝0的初始状态下，仅靠近壁面的第一层网格被标记，而每次迭代完成都更新UDM，标记满足条件的网格；整个计算时间内，流场域内是一个随时迭代、标记更新的过程；
水蒸气与冰的传质速率由两部分组成：
一部分为主要增大霜层厚度的体积传质速率
具体可以表示为：
其中，τv表示与相变速率有关的项，hsub＝2.833×106kJ/kg表示相变的固化潜热值，αa及ρa分别表示湿空气的体积分数和湿空气的密度，ω及ωsat分别表示水蒸气的实际质量分数和水蒸气处于饱和分压力时的饱和质量分数，
表示控制体的表面体积之比，
表示控制体到冷表面的导热传热量；
另一部分为主要提高霜层密度的密度传质速率
具体表示为：
其中，a表示与相变速率有关的项，ωa表示水空气的质量分数，ρa表示湿空气的密度，
表示控制体的表面体积之比，ρf,avg表示霜层的平均密度；
两者加和即为总传质速率：
S5、获得霜层的厚度及密度；
霜层的平均厚度及平均密度，可以用用户自定义函数计算，并用UDM储存记录；
霜层的平均厚度计算公式如下：
其中，δj中表示第j列网格中，冰项体积分数大于10-6部分的网格在y方向上的最大坐标，Δx＝1mm表示网格的长度，L＝100mm表示冷表面的总长度；对霜层厚度乘以网格x方向的增量加集合，再除以霜层总长L，即可求得霜层的平均厚度
霜层的平均密度计算公式如下：
其中，mcv表示控制体内霜的质量，Vcv表示控制体的体积，Vfl表示霜层的总体积；
S6、构建用于冷平板的霜层生长数值分析方法；
首先将物理模型离散化处理为数量众多的网格单元，然后在模拟过程中提取某一网格单元Cellx内流体的物性参数；求解动量、连续性、能量、组分、湍流、欧拉方程，得到体积分数、温度等物性参数，计算总相变速率，更新动量、质量、能量、组分源相；在霜层厚度以及霜层密度与流体速度场、温度场之间建立耦合关系；对流场内所有网格单元进行循环，并迭代至下一时间步长；反复循环直至所有网格单元均达到所需计算的时间要求；
S7、同工况下与相关文献实验结果进行对比，验证模型的准确性；
S8、将构建霜层生长模型用于实际，分析冷表面上霜层生长的分布情况。
进一步地，步骤S3中，所述广义源相公式如下：
冰项的质量源相
湿空气项的质量源相
冰项的能量源相
湿空气项的能量源相
湿空气项的x方向动量源相
湿空气项的y方向动量源相
湿空气项的组分源相
由欧拉多相流及组分输运方程可得，干空气的质量变化Sa＝Sm-a-Sn＝0。
进一步地，步骤S8中，由于霜层生长时，以水蒸气过饱和量为主体的相变驱动力和相变热阻分布不均匀，在冷平板的局部会产生较厚较密的霜层，而在其余位置霜层较薄较稀疏；霜层厚度随着流通距离的增长，呈现先增大再减小的趋势，霜层的密度在入口处远大于其余位置；霜层内部流线矢量图表明，在霜层的锋区只有很少的潮湿气流，在霜层的其他区域几乎没有潮湿气流，这也验证了相变驱动力霜层生长模型准确性。
进一步地，步骤S4中，当湿空气温度处于173.15K～273.15K的范围时，水蒸气饱和分压力计算公式如下：
psat(T)＝exp(a0/T+a1+a2T+a3T2+a4T3+a5T4+a6lnT)
而当湿空气温度处于273.15K～473.15K的范围时，水蒸气饱和分压力计算公式为：
psat(T)＝exp(b0/T+b1+b2T+b3T2+b4T3+b5lnT)
当控制体温度为T时，水蒸气对应的饱和分压力下饱和质量分数可以通过湿空气的热物理性质来确定，式中a0＝-5.674×103；a1＝6.392；a2＝-9.678× 10-3；a3＝6.211×10-7；a4＝2.075×10-9；a5＝-9.484×10-13；a6＝4.164；b0＝-5.800 ×10+3；b1＝1.391；b2＝-4.864×10-2；b3＝4.176×10-5；b4＝-1.445×10-8；b5＝ 6.546；
与T对应的饱和潮湿空气的湿度比为：
式中，p0为大气压力；
水蒸气的饱和质量分数由下式给出：
导热传热量
近似为：
式中，kf.cv和kf.p分别为控制体中霜的热导率和冷板附近霜的热导率；Tf.cv和Tp分别为控制体中的霜温度和冷板温度；dp为控制体与冷板之间的距离；
霜的热导率由以下一系列公式给出：
kpar＝(1-ψ)ki+ψka
其中，ki和ka分别为冰和空气的热导率，ψ为霜层的孔隙率；
其余的物理量ω、αa、ρa，由UDF直接提取。
与现有技术相比，本发明具有以下有益效果：
本发明的核心点在于，突出体现霜层表面体积分数的作用。在过去的研究方法中，普遍忽视了霜层表面对结霜的影响，但考虑到结霜实质上是水蒸气穿透霜层表面、渗透进霜层内部的传热传质行为，霜层表面的物理状况对结霜意义重大。由于结霜的发生实质是从表面晶核形成开始的，随即相邻的晶核之间相互联系构建起冰桥。这需要两方面的协同作用：首先，底部霜层需要能够支撑起上层冰晶体生长，这要求霜层的下方具有一定致密性，体现的是下层霜的支撑依托作用；其次，冰晶体继续生长时不同的晶体之间距离不能太远，否则冰桥就无法生成，冰晶体之间不能连接，这同样要求处于下方的霜层具有一定的致密性。综上所述，本方法补充了结霜的充要条件：同一列的网格中，对于即将相变的上层网格，其接触的正下方含霜网格的冰项体积分数必须大于临界阈值。
由于本发明这种需要上下层霜晶体相互依托、共同生长的特点，能够有效解决前期霜层在厚度上急剧增大的问题。与现有技术相比，本发明的前期误差相对较小，且误差较大的时间范围可由30分钟缩减至10分钟，大大提高了计算精准度。在结霜后期，霜层平均厚度及平均密度误差均不足5，对比同类技术也具有相当的优势。
本发明利用计算机辅助技术方法来实现，提出了一种冷表面结霜过程的模拟分析方法。使用欧拉多相流建立平板结霜工况的计算模型，并且用C编程语言编制仿真计算程序UDF对计算区域的控制方程进行反复迭代计算求解，模拟霜层增长情况。本发明可以准确的反应霜层在平板不同区域的结霜程度，可以有效地指导换热器的结构设计，从而有效地避免局部结霜过多而导致的制冷量下降以及换热器利用率不高的问题。
附图说明
图1为霜层生长机制示意图；
图2为霜层厚度方向生长示意图；
图3为利用霜层表面性质预测冷表面结霜的模拟方法的流程图；
图4为霜层平均厚度计算方法示意图；
图5为利用霜层表面性质预测冷表面结霜的模拟方法的实验数据厚度对比图；
图6为利用霜层表面性质预测冷表面结霜的模拟方法的实验数据密度对比图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图3所示，一种利用霜层表面性质预测冷表面结霜的模拟方法，包括以下步骤：
S1、建立冷表面物理模型；
由于平板的二维物理模型并不复杂，可利用ICEM软件直接建立。采用四边形结构化网格对计算域网格进行划分，并对近壁面网格在满足Y+的要求下进行加密处理。
S2、获得满足模拟精度的网格数量；
对划分的计算域网格进行网格无关性验证，通过自适应网格技术调整网格，同时分别划分几种不同数量的网格进行对比验证，随着网格数量的增加，霜层的总质量不在变化且差值逐渐趋于稳定时，即可得到满足模拟精度的网格数量。
S3、计算流场参数；
利用连续性、动量、能量和组分方程对整个流场区域进行模拟计算，各控制方程通用形式如下：
其中，ρ为流体密度；Φ为广义变量；Γ为广义扩散系数；S为广义源相。
其中，涉及计算所采用的7个源相公式如下：
冰项的质量源相
湿空气项的质量源相
冰项的能量源相
湿空气项的能量源相
湿空气项的x方向动量源相
湿空气项的y方向动量源相
湿空气项的组分源相
由欧拉多相流及组分输运方程可得，干空气的质量变化Sa＝Sm-a-Sn＝0。
S4、标记并获得水蒸气与冰的传质速率；
本发明的设计思路与元胞自动机的运行原理类似，通过用户自定义储存单元(UDM)把能否相变的区域分别储存标记为1和0。由于靠近壁面的第一层网格始终存在着冷表面的支撑作用，在整个计算时间内始终被标记为1；而那些处于同一列，且正下方接触着冰项体积分数大于临界体积分数的网格，也被标记为1。在t＝0的初始状态下，仅靠近壁面的第一层网格被标记，而每次迭代完成都更新UDM，标记满足条件的网格。整个计算时间内，流场域内是一个随时迭代、标记更新的过程。如图2所示，在流动时间为n时，紧挨壁面的网格1～9以及域中冰项体积分数大于临界体积分数的网格10～21被标记，这预示着经历一个时间步长后，处于这些网格的同一列且接触着的上方网格可以相变成霜。经过一个时间步长后，可以看到已有网格22～30结霜，其中网格22～27 中冰项的体积分数大于临界体积分数，而其余小于临界体积分数；再经过一个时间步长，即流动时间处于n+2时，网格28～30因为霜层不够致密，不能支撑上霜晶体的生成，所以只有网格31～36新生成了霜。计算过程中这三步依次循环进行，直到达到最大计算时间才会结束。
水蒸气与冰的传质速率由两部分组成：
一部分为主要增大霜层厚度的体积传质速率
具体可以表示为：
其中，τv表示与相变速率有关的项，hsub＝2.833×106kJ/kg表示相变的固化潜热值，αa及ρa分别表示湿空气的体积分数和湿空气的密度，ω及ωsat分别表示水蒸气的实际质量分数和水蒸气对应分压力下的饱和质量分数，
表示控制体的表面体积之比，
表示控制体到冷表面的导热传热量。
另一部分为主要提高霜层密度的密度传质速率
具体表示为：
其中，a表示与相变速率有关的项，ωa表示水空气的质量分数，ρa表示湿空气的密度，
表示控制体的表面体积之比，ρf,avg表示霜层的平均密度。
两者加和即为总传质速率：
上述所涉及的水蒸气饱和分压力和水蒸气处于饱和分压力时的饱和质量分数计算方式如下：
当湿空气温度处于173.15K～273.15K的范围时，水蒸气饱和分压力计算公式为：
psat(T)＝exp(a0/T+a1+a2T+a3T2+a4T3+a5T4+a6lnT)
而当湿空气温度处于273.15K～473.15K的范围时，水蒸气饱和分压力计算公式为：
psat(T)＝exp(b0/T+b1+b2T+b3T2+b4T3+b5lnT)
当控制体温度为T时，水蒸气对应的饱和分压力下饱和质量分数可以通过湿空气的热物理性质来确定，式中a0＝-5.674×103；a1＝6.392；a2＝-9.678× 10-3；a3＝6.211×10-7；a4＝2.075×10-9；a5＝-9.484×10-13；a6＝4.164；b0＝-5.800 ×10+3；b1＝1.391；b2＝-4.864×10-2；b3＝4.176×10-5；b4＝-1.445×10-8；b5＝ 6.546。
与T对应的饱和潮湿空气的湿度比为：
式中，p0为大气压力。
水蒸气的饱和质量分数由下式给出：
导热传热量
近似为：
式中，kf.cv和kf.p分别为控制体中霜的热导率和冷板附近霜的热导率；Tf.cv和Tp分别为控制体中的霜温度和冷板温度；dp为控制体与冷板之间的距离。
霜的热导率由以下一系列公式给出：
kpar＝(1-ψ)ki+ψka
其中，ki和ka分别为冰和空气的热导率，ψ为霜层的孔隙率。
其余的物理量，例如ω、αa、ρa，可由UDF直接提取。上述所涉及的全部公式，均用C语言编写成代码，经UDF加载进FLUENT中。
S5、获得霜层的厚度及密度；
霜层的平均厚度及平均密度，可以用用户自定义函数(UDF)计算，并用 UDM储存记录。
霜层的平均厚度计算公式如下：
其中，δj中表示第j列网格中，冰项体积分数大于10-6部分的网格在y方向上的最大坐标，Δx＝1mm表示网格的长度，L＝100mm表示冷表面的总长度。如图4所示，对霜层厚度乘以网格x方向的增量加集合，再除以霜层总长L，即可求得霜层的平均厚度
霜层的平均密度计算公式如下：
其中，mcv表示控制体内霜的质量，Vcv表示控制体的体积，Vfl表示霜层的总体积。
S6、构建用于冷平板的霜层生长数值分析方法；
首先将物理模型离散化处理为数量众多的网格单元，然后在模拟过程中提取某一网格单元Cellx内流体的物性参数。求解动量、连续性、能量、组分、湍流、欧拉方程，得到体积分数、温度等物性参数，计算总相变速率，更新动量、质量、能量、组分源相。在霜层厚度以及霜层密度与流体速度场、温度场之间建立耦合关系。对流场内所有网格单元进行循环，并迭代至下一时间步长。反复循环直至所有网格单元均达到所需计算的时间要求。
S7、同工况下与相关文献实验结果进行对比，验证模型的准确性；
主要是与相关实验文献结果的平均厚度、平均密度随时间的变化相对比，并与同类模拟文献的云图趋势对比，以验证模型的准确性。
与对应的参考文献(HermesC,PiuccoRO,BarbosaJR,etal.Erratum to:"Astudyoffrostgrowthanddensificationonflatsurfaces", ExperimentalThermalandFluidScience33(2009)371–379[J].2009, 33(6):1035-1035.)数据进行对比。在湿空气入口温度Tin＝289.15K，入口流速Uin＝0.7m/s，水蒸气相对湿度RH＝80％，冷表面Tp＝257.15K、261.15K及265.15K 的边界条件下，以1s为时间步长，计算了7200s。为了方便收敛，所有项的松弛因子取0.3，残差因子取10-6。
在7200s的计算时间内，计算所得的平均厚度和平均密度均与实验数据相近。如图5-6所示，在计算的前600s为大误差区间，厚度误差小于5％而密度最大误差不超过10％，在计算的后期密度及厚度误差均不超过5％。模拟结果与实验数据基本趋于一致，验证霜层生长方法的准确性。
S8、将构建霜层生长模型用于实际，分析冷表面上霜层生长的分布情况。
由于霜层生长时，以水蒸气过饱和量为主体的相变驱动力和相变热阻分布不均匀，在冷平板的局部会产生较厚较密的霜层，而在其余位置霜层较薄较稀疏。霜层厚度随着流通距离的增长，呈现先增大再减小的趋势，霜层的密度在入口处远大于其余位置。霜层内部流线矢量图表明，在霜层的锋区只有很少的潮湿气流，在霜层的其他区域几乎没有潮湿气流，这也验证了相变驱动力霜层生长模型准确性。
以上公开的本发明优选实施例只是用于帮助阐述本发明。优选实施例并没有详尽叙述所有的细节，也不限制该发明仅为所述的具体实施方式。显然，根据本说明书的内容，可作很多的修改和变化。本说明书选取并具体描述这些实施例，是为了更好地解释本发明的原理和实际应用，从而使所属技术领域人员能很好的理解和利用本发明。本发明仅受权利要求书及其全部范围和等效物的限制。
Claims (4)
1.一种利用霜层表面性质预测冷表面结霜的模拟方法，其特征在于：包括以下步骤：
S1、建立冷表面物理模型；
利用ICEM软件建立冷表面物理模型，采用四边形结构化网格对计算域网格进行划分，并对近壁面网格在满足Y+的要求下进行加密处理；
S2、获得满足模拟精度的网格数量；
对划分的计算域网格进行网格无关性验证，通过自适应网格技术调整网格，同时分别划分几种不同数量的网格进行对比验证，随着网格数量的增加，霜层的总质量不在变化且差值逐渐趋于稳定时，即可得到满足模拟精度的网格数量；
S3、计算流场参数；
利用连续性、动量、能量和组分方程对整个流场区域进行模拟计算，各控制方程通用形式如下：
其中，ρ为流体密度；Φ为广义变量；Γ为广义扩散系数；S为广义源相；
S4、标记并获得水蒸气与冰的传质速率；
通过用户自定义UDM把能否相变的区域分别储存标记为1和0；由于靠近壁面的第一层网格始终存在着冷表面的支撑作用，在整个计算时间内始终被标记为1；而那些处于同一列，且正下方接触着冰项体积分数大于临界体积分数的网格，也被标记为1；在t＝0的初始状态下，仅靠近壁面的第一层网格被标记，而每次迭代完成都更新UDM，标记满足条件的网格；整个计算时间内，流场域内是一个随时迭代、标记更新的过程；
水蒸气与冰的传质速率由两部分组成：
一部分为主要增大霜层厚度的体积传质速率
具体可以表示为：
其中，τv表示与相变速率有关的项，hsub＝2.833×106kJ/kg表示相变的固化潜热值，αa及ρa分别表示湿空气的体积分数和湿空气的密度，ω及ωsat分别表示水蒸气的实际质量分数和水蒸气处于饱和分压力时的饱和质量分数，
表示控制体的表面体积之比，
表示控制体到冷表面的导热传热量；
另一部分为主要提高霜层密度的密度传质速率
具体表示为：
其中，a表示与相变速率有关的项，ωa表示水空气的质量分数，ρa表示湿空气的密度，
表示控制体的表面体积之比，ρf,avg表示霜层的平均密度；
两者加和即为总传质速率：
S5、获得霜层的厚度及密度；
霜层的平均厚度及平均密度，可以用用户自定义函数计算，并用UDM储存记录；
霜层的平均厚度计算公式如下：
其中，δj中表示第j列网格中，冰项体积分数大于10-6部分的网格在y方向上的最大坐标，Δx＝1mm表示网格的长度，L＝100mm表示冷表面的总长度；对霜层厚度乘以网格x方向的增量加集合，再除以霜层总长L，即可求得霜层的平均厚度
霜层的平均密度计算公式如下：
其中，mcv表示控制体内霜的质量，Vcv表示控制体的体积，Vfl表示霜层的总体积；
S6、构建用于冷平板的霜层生长数值分析方法；
首先将物理模型离散化处理为数量众多的网格单元，然后在模拟过程中提取某一网格单元Cellx内流体的物性参数；求解动量、连续性、能量、组分、湍流、欧拉方程，得到体积分数、温度等物性参数，计算总相变速率，更新动量、质量、能量、组分源相；在霜层厚度以及霜层密度与流体速度场、温度场之间建立耦合关系；对流场内所有网格单元进行循环，并迭代至下一时间步长；反复循环直至所有网格单元均达到所需计算的时间要求；
S7、同工况下与相关文献实验结果进行对比，验证模型的准确性；
S8、将构建霜层生长模型用于实际，分析冷表面上霜层生长的分布情况。
2.根据权利要求1所述的利用霜层表面性质预测冷表面结霜的模拟方法，其特征在于：所述步骤S3中，所述广义源相公式如下：
冰项的质量源相
湿空气项的质量源相
冰项的能量源相
湿空气项的能量源相
湿空气项的x方向动量源相
湿空气项的y方向动量源相
湿空气项的组分源相
由欧拉多相流及组分输运方程可得，干空气的质量变化Sa＝Sm-a-Sn＝0。
3.根据权利要求2所述的利用霜层表面性质预测冷表面结霜的模拟方法，其特征在于：所述步骤S4中，水蒸气的饱和分压力和水蒸气对应分压力下的饱和质量分数计算方式如下：
当湿空气温度处于173.15K～273.15K的范围时，水蒸气饱和分压力计算公式为：
psat(T)＝exp(a0/T+a1+a2T+a3T2+a4T3+a5T4+a6lnT)
而当湿空气温度处于273.15K～473.15K的范围时，水蒸气饱和分压力计算公式为：
psat(T)＝exp(b0/T+b1+b2T+b3T2+b4T3+b5lnT)
当控制体温度为T时，水蒸气对应的饱和分压力下饱和质量分数可以通过湿空气的热物理性质来确定，式中a0＝-5.674×103；a1＝6.392；a2＝-9.678×10-3；a3＝6.211×10-7；a4＝2.075×10-9；a5＝-9.484×10-13；a6＝4.164；b0＝-5.800×10+3；b1＝1.391；b2＝-4.864×10-2；b3＝4.176×10-5；b4＝-1.445×10-8；b5＝6.546；
与T对应的饱和潮湿空气的湿度比为：
式中，p0为大气压力；
水蒸气的饱和质量分数由下式给出：
导热传热量
近似为：
式中，kf.cv和kf.p分别为控制体中霜的热导率和冷板附近霜的热导率；Tf.cv和Tp分别为控制体中的霜温度和冷板温度；dp为控制体与冷板之间的距离；
霜的热导率由以下一系列公式给出：
kpar＝(1-ψ)ki+ψka
其中，ki和ka分别为冰和空气的热导率，ψ为霜层的孔隙率；
其余的物理量ω、αa、ρa，由UDF直接提取。
4.根据权利要求3所述的利用霜层表面性质预测冷表面结霜的模拟方法，其特征在于：所述步骤S8中，由于霜层生长时，以水蒸气过饱和量为主体的相变驱动力和相变热阻分布不均匀，在冷平板的局部会产生较厚较密的霜层，而在其余位置霜层较薄较稀疏；霜层厚度随着流通距离的增长，呈现先增大再减小的趋势，霜层的密度在入口处远大于其余位置；霜层内部流线矢量图表明，在霜层的锋区只有很少的潮湿气流，在霜层的其他区域几乎没有潮湿气流，这也验证了相变驱动力霜层生长模型准确性。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210569667.8A CN115114815B (zh) | 2022-05-24 | 2022-05-24 | 一种利用霜层表面性质预测冷表面结霜的模拟方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210569667.8A CN115114815B (zh) | 2022-05-24 | 2022-05-24 | 一种利用霜层表面性质预测冷表面结霜的模拟方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115114815A true CN115114815A (zh) | 2022-09-27 |
| CN115114815B CN115114815B (zh) | 2024-06-07 |
Family
ID=83326581
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210569667.8A Active CN115114815B (zh) | 2022-05-24 | 2022-05-24 | 一种利用霜层表面性质预测冷表面结霜的模拟方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115114815B (zh) |
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116502426A (zh) * | 2023-04-19 | 2023-07-28 | 东北电力大学 | 基于凝固融化过程结霜模拟方法、装置、终端和介质 |
| CN117172076A (zh) * | 2023-10-18 | 2023-12-05 | 中国市政工程华北设计研究总院有限公司 | 一种lng空温式气化器在结霜条件下的性能模拟方法 |
| CN117272876A (zh) * | 2023-11-23 | 2023-12-22 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种强对流条件下平板结霜霜层物性关联式的建立方法 |
| CN117408054A (zh) * | 2023-10-18 | 2024-01-16 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种基于结霜相似律的圆管结霜平均厚度预测方法 |
| CN117407634A (zh) * | 2023-10-18 | 2024-01-16 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种基于结霜特征曲线的平板结霜厚度快速预测方法 |
| CN117408053A (zh) * | 2023-10-18 | 2024-01-16 | 中国空气动力研究与发展中心计算空气动力研究所 | 强对流条件下低温平板干模态结霜特征曲线的建立方法 |
| CN117407635A (zh) * | 2023-10-18 | 2024-01-16 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种基于结霜相似律的平板结霜厚度预测方法 |
| CN117494400A (zh) * | 2023-10-18 | 2024-02-02 | 中国空气动力研究与发展中心计算空气动力研究所 | 强对流条件下低温平板干模态结霜相似因素获得方法 |
Citations (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20160292327A1 (en) * | 2015-04-03 | 2016-10-06 | Airbus Group India Private Limited | Estimating frost mass formed in a design component of a multi-component structure |
| US20160292328A1 (en) * | 2015-04-03 | 2016-10-06 | Airbus Group India Private Limited | Estimating frost mass formed in a design component of a multi-component structure |
| CN106682376A (zh) * | 2017-04-01 | 2017-05-17 | 国网河南省电力公司电力科学研究院 | 参数随工况变化实际特性的全过程汽轮机建模及辨识方法 |
| US20180312763A1 (en) * | 2016-04-07 | 2018-11-01 | Worthington Industries, Inc. | System for reclaiming liquefied petroleum gas |
| CN109709138A (zh) * | 2018-12-29 | 2019-05-03 | 北京卫星环境工程研究所 | 真空低温凝华结霜的测试***和测试方法 |
| CN111159890A (zh) * | 2019-12-28 | 2020-05-15 | 中汽研汽车检验中心（天津）有限公司 | 一种用于抑制预冷器结霜的模拟计算方法 |
| CN111695242A (zh) * | 2020-05-19 | 2020-09-22 | 东南大学 | 一种湿饱和烟气蒸汽凝结的数值模拟方法 |
| CN112082146A (zh) * | 2020-08-14 | 2020-12-15 | 东北电力大学 | 火电机组旁路供热用减温减压器液滴蒸发段长度的确定方法 |
| CN112800700A (zh) * | 2021-04-13 | 2021-05-14 | 中国空气动力研究与发展中心计算空气动力研究所 | 低温表面干模态结霜模拟方法、装置、电子设备和介质 |
| CN113379925A (zh) * | 2021-07-30 | 2021-09-10 | 广东工业大学 | 一种冷表面结霜模拟图像生成方法及装置 |
| CN114322422A (zh) * | 2021-12-09 | 2022-04-12 | 西安交通大学 | 一种冷表面结霜量测量方法及应用 |
-
2022
- 2022-05-24 CN CN202210569667.8A patent/CN115114815B/zh active Active
Patent Citations (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20160292327A1 (en) * | 2015-04-03 | 2016-10-06 | Airbus Group India Private Limited | Estimating frost mass formed in a design component of a multi-component structure |
| US20160292328A1 (en) * | 2015-04-03 | 2016-10-06 | Airbus Group India Private Limited | Estimating frost mass formed in a design component of a multi-component structure |
| US20180312763A1 (en) * | 2016-04-07 | 2018-11-01 | Worthington Industries, Inc. | System for reclaiming liquefied petroleum gas |
| CN106682376A (zh) * | 2017-04-01 | 2017-05-17 | 国网河南省电力公司电力科学研究院 | 参数随工况变化实际特性的全过程汽轮机建模及辨识方法 |
| CN109709138A (zh) * | 2018-12-29 | 2019-05-03 | 北京卫星环境工程研究所 | 真空低温凝华结霜的测试***和测试方法 |
| CN111159890A (zh) * | 2019-12-28 | 2020-05-15 | 中汽研汽车检验中心（天津）有限公司 | 一种用于抑制预冷器结霜的模拟计算方法 |
| CN111695242A (zh) * | 2020-05-19 | 2020-09-22 | 东南大学 | 一种湿饱和烟气蒸汽凝结的数值模拟方法 |
| CN112082146A (zh) * | 2020-08-14 | 2020-12-15 | 东北电力大学 | 火电机组旁路供热用减温减压器液滴蒸发段长度的确定方法 |
| CN112800700A (zh) * | 2021-04-13 | 2021-05-14 | 中国空气动力研究与发展中心计算空气动力研究所 | 低温表面干模态结霜模拟方法、装置、电子设备和介质 |
| CN113379925A (zh) * | 2021-07-30 | 2021-09-10 | 广东工业大学 | 一种冷表面结霜模拟图像生成方法及装置 |
| CN114322422A (zh) * | 2021-12-09 | 2022-04-12 | 西安交通大学 | 一种冷表面结霜量测量方法及应用 |
Non-Patent Citations (4)
| Title |
|---|
| ZHIMIN HAN 等: "A frost model based on the frost layer\'s supporting function", 《JOURNAL OF HEAT AND MASS TRANSFER》, 31 March 2023 (2023-03-31), pages 1 - 12 * |
| 崔静;杨霆浩;杨帆;李虎林;杨广峰;: "基于改进格子Boltzmann焓法模型的霜层生长数值研究", 计算物理, no. 03, 8 June 2018 (2018-06-08), pages 75 - 86 * |
| 文力;金旭;刘忠彦;张阔;吴哲;: "喷射特征的双级压缩***压缩过程及喷射参数耦合特性", 流体机械, no. 11, 30 November 2019 (2019-11-30), pages 80 - 86 * |
| 李翠;庄钰涵;厉彦忠;赵志翔;: "低温管路管外结霜过程的数值模拟研究", 低温工程, no. 04, 15 August 2018 (2018-08-15), pages 5 - 11 * |
Cited By (15)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116502426A (zh) * | 2023-04-19 | 2023-07-28 | 东北电力大学 | 基于凝固融化过程结霜模拟方法、装置、终端和介质 |
| CN117407635A (zh) * | 2023-10-18 | 2024-01-16 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种基于结霜相似律的平板结霜厚度预测方法 |
| CN117494400A (zh) * | 2023-10-18 | 2024-02-02 | 中国空气动力研究与发展中心计算空气动力研究所 | 强对流条件下低温平板干模态结霜相似因素获得方法 |
| CN117408054A (zh) * | 2023-10-18 | 2024-01-16 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种基于结霜相似律的圆管结霜平均厚度预测方法 |
| CN117407634A (zh) * | 2023-10-18 | 2024-01-16 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种基于结霜特征曲线的平板结霜厚度快速预测方法 |
| CN117408053A (zh) * | 2023-10-18 | 2024-01-16 | 中国空气动力研究与发展中心计算空气动力研究所 | 强对流条件下低温平板干模态结霜特征曲线的建立方法 |
| CN117172076A (zh) * | 2023-10-18 | 2023-12-05 | 中国市政工程华北设计研究总院有限公司 | 一种lng空温式气化器在结霜条件下的性能模拟方法 |
| CN117407635B (zh) * | 2023-10-18 | 2024-05-14 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种基于结霜相似律的平板结霜厚度预测方法 |
| CN117172076B (zh) * | 2023-10-18 | 2024-05-10 | 中国市政工程华北设计研究总院有限公司 | 一种lng空温式气化器在结霜条件下的性能模拟方法 |
| CN117408054B (zh) * | 2023-10-18 | 2024-04-12 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种基于结霜相似律的圆管结霜平均厚度预测方法 |
| CN117494400B (zh) * | 2023-10-18 | 2024-04-19 | 中国空气动力研究与发展中心计算空气动力研究所 | 强对流条件下低温平板干模态结霜相似因素获得方法 |
| CN117407634B (zh) * | 2023-10-18 | 2024-05-03 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种基于结霜特征曲线的平板结霜厚度快速预测方法 |
| CN117408053B (zh) * | 2023-10-18 | 2024-05-07 | 中国空气动力研究与发展中心计算空气动力研究所 | 强对流条件下低温平板干模态结霜特征曲线的建立方法 |
| CN117272876A (zh) * | 2023-11-23 | 2023-12-22 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种强对流条件下平板结霜霜层物性关联式的建立方法 |
| CN117272876B (zh) * | 2023-11-23 | 2024-01-26 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种强对流条件下平板结霜霜层物性关联式的建立方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115114815B (zh) | 2024-06-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN115114815A (zh) | 一种利用霜层表面性质预测冷表面结霜的模拟方法 | |
| Cui et al. | A new time-and space-dependent model for predicting frost formation | |
| Aschwanden et al. | An enthalpy formulation for glaciers and ice sheets | |
| US20220027539A1 (en) | Method and system for manufacturing a heat exchanger for supercritical pressure fluid | |
| Gadgil | On convective heat transfer in building energy analysis | |
| Bartrons et al. | A finite volume method to solve the frost growth using dynamic meshes | |
| Chen et al. | Simulation of frost growth and densification on horizontal plates with supersaturated interface condition | |
| Wu et al. | Numerical simulation of subcooled nucleate boiling by coupling level-set method with moving-mesh method | |
| Tahavvor et al. | Experimental and numerical study of frost formation by natural convection over a cold horizontal circular cylinder | |
| Fu et al. | Transient laminar natural convection in an enclosure partitioned by an adiabatic baffle | |
| Wang et al. | Research on falling film dehumidification performance of microencapsulated phase change materials slurry | |
| Yuan et al. | Multiscale and multilayer structural modeling and simulation on drying of grain packing porous media | |
| Xu et al. | A frost model based on the frost layer's supporting function | |
| Mobli et al. | Estimating bubble interfacial heat transfer coefficient in pool boiling | |
| Zhao et al. | Exploring the limits of condensation heat transfer: A numerical study of microscale-confined condensation between parallel surfaces having wetting contrast | |
| CN117172076B (zh) | 一种lng空温式气化器在结霜条件下的性能模拟方法 | |
| Chen et al. | Simulation of heat and mass transfer in a grain pile on the basis of a 2D irregular pore network | |
| Navid et al. | A new approach to delay or prevent frost formation in membranes | |
| Betchen et al. | An investigation of the effects of a linear porosity distribution on non-equilibrium heat transfer in high-conductivity graphitic foam | |
| CN116882253A (zh) | 基于Modelica的热构件建模方法、设备及介质 | |
| Green et al. | Marginal stability in inhomogeneous porous media | |
| Harges et al. | Modeling of frost growth on surfaces with varying contact angle | |
| Yang et al. | A numerical simulation of pool boiling using CAS model | |
| Hervatte | CFD Simulation of a Fin-Tube evaporator under icing | |
| CN113051849B (zh) | 考虑气体成分变化的低温推进剂贮箱增压性能预测方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant |