我国是牛肉生产和消费大国，2019 年我国牛肉产量仅次于美国、巴西和欧盟；牛肉消费量位于世界第二位，仅次于美国[1-2]。随着人们生活水平的提高，牛肉需求量越来越大的同时，大家对牛肉品质要求也随之提高。牛胴体预冷排酸对牛肉品质有重要意义[3-4]。牛被宰杀后体温升高，胴体中心温度可达40 ℃，成为滋生细菌的温床，低温预冷能够很好地抑制细菌生长繁殖，保持牛肉品质。另外，牛在宰杀的时候，因为恐惧等因素，导致肉质坚硬，干躁，缺少弹性，且肌肉中的糖元无氧酵解产生乳酸，不仅影响牛肉品质，还不利于人体健康。低温预冷排酸，可将乳酸分解成二氧化碳、水和酒精然后挥发。由于低温预冷排酸经历了较为充分的解僵、成熟过程，不仅肉质柔软有弹性、味道鲜美，而且安全营养[5-6]。因此，牛宰后胴体预冷排酸是生产高品质牛肉的第一个环节也是必要环节，通过预冷改善牛肉品质成为研究热点。
然而，国内预冷库大多是根据经验来设计，导致预冷阶段牛胴体温度分布不均匀，浪费能源的同时还影响牛肉品质。利用数值模拟技术研究牛肉胴体预冷库温度时空分布，使温度场在预冷库中分布更加均匀，对改善牛肉品质具有重要意义。计算流体动力学(computational fluid dynamics，CFD)是基于计算机技术的一种数值计算工具，用于求解流体的流动和换热问题[7]。CFD技术被国内外研究者广泛应用于研究低温冷却过程中的温度场分布。
文献[8-11]分别模拟了苹果、番茄、蓝莓、茎状蔬菜等冷却时的温度变化情况，为改善果蔬冷却效果提供了参考。刘泽勤等[12]以压差原理的果蔬仓库为研究对象，对稳态下的库内温度分布情况进行模拟仿真，讨论了果蔬间距的改变对库内果蔬温湿度的影响。赵春江等[13]建立了求解短距离运输的冷藏车车厢内温度场分布的计算模型，模拟分析了不同边界条件和货物不同堆栈方式对车厢内温度场分布的影响。赵时等[14]建立了冷藏车厢的仿真模型,利用CFD技术研究梯级送风对空仓时冷藏车厢内温度场的影响。国内虽然针对果蔬包装箱、冷藏车、冷库中温度场的研究比较多，但目前利用CFD技术对牛胴体预冷库研究的很少。而国外对于牛胴体预冷库研究大多是一次预冷一个或几个牛胴体[15-18]，并没有研究实际上的一次预冷几十上百个胴体，无法体现预冷时温度场的均匀性，在提高牛胴体预冷效果方面缺少实用价值。为此，本文以吊挂有76个牛胴体二分体的实际预冷库为研究对象，利用CFD技术进行模拟计算，分析不同预冷时间不同位置的温度分布，对比不同时间不同位置的降温快慢、冷却均匀性，客观评价预冷效果。
1 材料与方法
1.1 物理模型
预冷库物理模型根据实地测量北京某畜牧有限公司牛肉胴体预冷库得出。牛胴体二分体模型参考实际尺寸并做适当简化，应用3 dmax进行三维建模绘制。预冷库空间尺寸(长×宽×高)为1.6 m×5.25 m×3.8 m，4个冷风机空间尺寸(长×宽×高)都为 1.6 m×0.6 m×0.65 m，风机距天花板0.35 m，距墙壁0.6 m。预冷库进风口为圆形，位于冷风机前部，直径0.5 m，风速6 m/s，温度-5 ℃。预冷库出风口为圆形，位于冷风机背面，直径0.5 m。本次模拟实验一次性预冷76个二分体。4根吊挂轨道，每根吊轨吊挂19个二分体。吊轨间距1 m，相邻吊钩间距0.8 m。牛胴体二分体长2.3 m、宽0.75 m、厚0.35 m。预冷库三维图(省略了吊挂牛胴体的吊钩)及牛胴体二分体模型如图1所示。
图1 预冷库结构示意图
Fig.1 Schematic diagram of pre-cooling room
1.2 网格划分
本文采用软件 ICEM CFD 15.0进行网格划分，整个预冷库及内部结构采用六面体非结构化网格进行划分。全局网格最大尺寸不超过0.2 m，壁面及牛胴体最大网格尺寸不超过0.1 m。网格数量为1 621 719个单元，划分结果如图2所示。非结构网格如果质量太差，在Fluent软件计算过程中容易导致计算发散或者收敛困难。网格质量正交性指标分布从0到 1，越接近0网格质量越差。利用mesh quality对网格的正交性进行质量检测，结果表明网格质量均大于0.35，平均网格质量为0.77，表示网格质量良好。
a-预冷库整体网格;b-牛胴体网格
图2 网格划分
Fig.2 Mesh generation
1.3 数学模型
建立数学模型求解牛胴体二分体预冷过程的温度场分布情况时，为了达到预期效果的同时尽量减少计算时间，研究对实际模型进行简化，做出如下假设：空气为不可压缩气体且符合Boussinesq假设；预冷库内管道、铁架等对流场的影响忽略不计；冷风机内部温度场变化对整个预冷库内温度场无影响；忽略墙壁与外界的热交换；全部牛二分体胴体是同时进入预冷库的；牛胴体的物理参数不随温度的变化而改变，忽略牛胴体预冷时的重量损失。
预冷库中空气、牛胴体[19]的具体热物理性能参数如表1 所示。
表1 热物性能参数
Table 1 Parameters of thermal-physical properties
材料密度/[kg·(m3)-1]比热容/[J·(kg·K)-1]导热率/[W·(m·K)-1]牛胴体空气1 0801.293 3201 0050.410.024 3
1.4 初始条件与边界条件
初始条件:当时间 t=0 时，预冷库壁面、预冷库内部空气、冷风机表面的初始温度 T1=0 ℃；牛胴体的初始温度设定为T2=40 ℃，略高于正常体温，这是考虑到牛在被宰杀时肌肉活动增强引起体温升高而设定的[19]。
入口边界:将4个冷风机上靠近牛胴体的8个圆孔设置为速度入口边界条件，风速6 m/s、温度-5 ℃。
出口边界:将4个冷风机上远离牛胴体的8个圆孔设置为出口边界条件，边界上所有物理量梯度为零。
壁面边界:预冷库壁面、牛胴体表面、冷风机表面设置为无滑移壁面条件，壁面上速度为零，且垂直于壁面的速度也为零。
1.5 数值模拟方法
采用基于有限体积法的 CFD 商用软件 Fluent 15.0 进行求解计算和后处理分析。不考虑重力影响，利用非稳态剪切压力传输(shear stress transport，SST)k-ω湍流模型、基于压力的分离式求解器进行计算。在空气流体区采用的控制方程有能量守恒方程、动量守恒方程和质量守恒方程，压力速度耦合方法采用 SIMPLE(semi-implicit method for pressure-linked equations)算法求解。时间步长设置为1 h，步数设为72步，共计72 h。
2 结果与分析
2.1 云图分析
图3计算并监测的是预冷库正中间一平面即x=8 m，不同预冷时间下的温度分布云图。从图3可以看出，牛胴体下部分温度明显高于上部分温度。原因从图4风速场基本稳定后的风速云图可知，冷风从上面吹过后，从下面返回。冷风从上面把牛胴体部分热量带到下面，导致同一时间，冷库下部分冷风温度高于上部分温度。另外胴体下部分相对较宽较大，也是预冷效果不如上部分的原因之一。
a-监测平面x=8 m；b-温度云图
图3 平面x=8 m处不同预冷时间下的温度分布云图
Fig.3 Temperature contours at plane x=8 m
图4 风速分布云图
Fig.4 Velocity distribution
图5显示的是不同时间下胴体表面的温度分布云图，可以看出，正对风机的牛胴体，降温速率明显高于非正对风机的牛胴体。在条件允许范围内，建议定时改变风机吹风方向，或者定时调换正对风机与非正对风机牛胴体的位置，使牛胴体冷却更均匀。
图5 不同预冷时间下的牛胴体表面温度分布云图
Fig.5 Temperature contours of beef carcass surface
2.2 预冷时间
依据《牛羊屠宰与分割车间设计规范》(GB 51225—2017) [20]，牛胴体冷却后中心温度不应高于7 ℃。所以当胴体中心平均温度达到7 ℃时，胴体基本达到预期预冷温度，所用时间为预冷时间。图6是所有牛胴体表面和胴体内部平均温度随预冷时间变化的曲线。由图6可知，所有胴体表面平均温度达到7 ℃ 需要18 h。所有胴体内部平均温度达到7 ℃，需要40 h。通过2.1节图3的温度云图可知，要使胴体中心最高温度低于7 ℃，则需要72 h左右。
图6 牛胴体平均温度变化曲线
Fig.6 Average temperature change curve of beef carcass
图7是预冷库中平均温度低于7 ℃的牛胴体个数随预冷时间增长而增加的曲线。从图7可看出，预冷28 h后开始有1个胴体平均温度低于7 ℃，从 30 h 到44 h，平均温度低于7 ℃胴体个数呈线性增长趋势。38 h后超过一半胴体平均温度低于7 ℃，这是通常所说的二分之一预冷时间。42 h后有超过八分之七的胴体平均温度低于7 ℃，这是通常所说的八分之七预冷时间。44 h后每个胴体平均温度都低于7 ℃。因为28 h后陆续有胴体平均温度低于7 ℃，建议在28 h后，将先预冷好的牛胴体陆续移出预冷库进行下一步加工处理。
图7 平均温度低于7 ℃的牛胴体个数曲线
Fig.7 Curve of beef carcass number with average temperature below 7 ℃
2.3 冷却均匀性
温度离散系数用来评估冷却均匀性，温度离散系数越大，冷库内牛胴体温度数据离散程度越大，温度分布越不均匀;反之，温度离散系数越小，冷库内牛胴体温度数据离散程度越小，温度分布越均匀。采用热力学温度计算温度离散系数[21]，计算如公式(1)所示：
(1)
式中：Cv为温度离散系数，%；n为牛胴体数量；Tave为n个牛胴体的平均温度；Ti为牛胴体i的平均温度。
以离冷风机由近到远分成第1排、第2排、第3排、第4排牛胴体，图8是4排牛胴体和所有牛胴体的温度离散系数曲线图。从图8可以看出，各排和整体温度离散系数都是先升高再降低。当温度离散系数达到峰值后，随着冷却时间的增加，预冷库内胴体间温度差异性越来越小，温度离散系数越来越小。整体温度离散系数峰值为0.5，各排牛胴体温度离散系数变化趋势相近，但后两排比前两排晚18 h达到峰值。前两排在12 h达到峰值，后两排在30 h分别达到峰值。这可能是因为前两排牛胴体离冷风机近，牛胴体降温相对快于后两排，随着温度的降低，前两排胴体间温度差异先减小，温度离散系数先降低。
Cv_1-第1排牛胴体；Cv_2-第2排牛胴体；Cv_3-第3排牛胴体； Cv_4-第4排牛胴体；Cv_all-所有牛胴体
图8 温度离散系数曲线
Fig.8 Curves of temperature variation coefficients
3 结论
风机朝一个方向吹风，导致牛胴体上部分比下部分降温快，非正对冷风机的牛胴体冷却效果不如正对冷风机的牛胴体。在条件允许范围内，建议改进风机，让风机朝不同方向发散吹风或定时改变风机吹风方向，或者定时调换正对风机与非正对风机牛胴体的位置，减少牛胴体温度离散系数，使牛胴体冷却更均匀。
冷却28 h后，开始有牛胴体平均温度低于7 ℃；44 h后，每个胴体平均温度都低于7 ℃；要使胴体中心最高温度低于7 ℃，则需要72 h左右。建议在28 h 后，将先预冷好的牛胴体陆续移出预冷库进行下一步加工处理。
对于不同冷库大小、不同初始温度、不同预冷胴体大小及个数，可参考本研究，适当选择风机(包括风机个数、风速、送风温度等)，适当调整预冷时间。
[1] 曹兵海, 李俊雅, 王之盛, 等.2019年度肉牛牦牛产业技术发展报告[J].中国畜牧杂志, 2020, 56(3):173-178.
CAO B H, LI J Y, WANG Z S, et al.2019 beef yak industry technology development report[J].Chinese Journal of Animal Science, 2020, 56(3):173-178.
[2] 季现秋, 罗欣, 朱立贤, 等.新型牛肉嫩化技术研究进展[J].食品与发酵工业, 2021, 47(1):327-333.
JI X Q, LUO X, ZHU L X, et al.Research progress on new beef tenderization technology[J].Food and Fermentation Industries, 2021, 47(1):327-333.
[3] 王恒鹏, 王引兰, 吴鹏, 等.排酸牛肉的品质变化规律及其烹煮时机的优选[J].食品与发酵工业, 2020, 46(1):222-228.
WANG H P, WANG Y L, WU P, et al.Analysis of quality changes and preferred choice of cooking time of beef during aging[J].Food and Fermentation Industries, 2020, 46(1):222-228.
[4] TRUJILLO F J, PHAM Q T.Modelling the chilling of the leg, loin and shoulder of beef carcasses using an evolutionary method[J].International Journal of Refrigeration, 2003, 26(2):224-231.
[5] 付丽, 杨兆华, 高雪琴, 等.牛肉排酸成熟过程中的品质变化[J].肉类研究, 2018, 32(11):9-14.
FU L, YANG Z H, GAO X Q, et al.Beef quality changes during postmortem aging[J].Meat Research, 2018, 32(11):9-14.
[6] 陈槟颖, 靳烨, 袁倩, 等.不同排酸时间对呼伦贝尔肉羊宰后品质的影响[J].肉类研究, 2014, 28(1):1-3.
CHEN B Y, JIN Y, YUAN Q, et al, Influence of different chilling methods on meat quality of hulunbeir sheep[J].Meat Research, 2014, 28(1):1-3.
[7] 谢晶, 瞿晓华, 徐世琼.冷藏库内气体流场数值模拟与验证[J].农业工程学报, 2005(2):11-16.
XIE J, QU X H, XU S Q.Numerical simulation and verification of airflow in cold-store[J].Transactions of the Chinese Society of Agricultural Engineering, 2005, 21(2):11-16.
[8] 王喜芳, 李保国, 范中阳.果品通风包装预冷效果的数值模拟与实验研究[J].包装工程, 2019, 40(5):19-25.
WANG X F, LI B G, FAN Z Y.Numerical simulation and experiment of precooling efficiency inside fruit ventilated package[J].Packaging Engineering, 2019, 40(5):19-25.
[9] 苏勤, 谌英敏, 柏惠康, 等.基于计算流体力学的开孔均匀性对番茄预冷性能的影响[J].食品与发酵工业, 2020, 46(21):167-172.
SU Q, SHEN Y M, BAI H K, et al.The effect of CFD-based hole uniformity on tomato precooling performance[J].Food and Fermentation Industries, 2020, 46(21):167-172.
[10] 王达, 杨相政, 贾斌广, 等.不同包装结构对蓝莓压差预冷效果的影响[J].浙江大学学报(农业与生命科学版), 2020, 46(1):47-54;63.
WANG D, YANG X Z, JIA B G, et al.Influence of different packaging structures on forced air pre-cooling effect of blueberry[J].Journal of Zhejiang University(Agriculture and Life Sciences), 2020, 46(1):47-54;63.
[11] 阚安康, 黄孜沛, 张婷婷, 等.茎状蔬菜真空预冷热质传递模拟与实验研究[J].上海海事大学学报, 2020, 41(3):122-126.
KAN A K, HUANG Z P, ZHANG T T, et al.Numerical simulation and experimental study on heat and mass transfer of cylindrical vegetables during vacuum precooling[J].Journal of Shanghai Maritime University, 2020, 41(3):122-126.
[12] 刘泽勤, 刘丽芳, 徐静.冷藏库堆码间距对库内温湿度场的影响研究[J].低温与超导, 2020, 48(10):79-84;95.
LIU Z Q, LIU L F, XU J.Study on the influence of stacking space on temperature and humidity fields in cylindrical and conical cold storage[J].Cryogenics & Superconductivity, 2020, 48(10):79-84；95.
[13] 赵春江, 韩佳伟, 杨信廷, 等.基于CFD的冷藏车车厢内部温度场空间分布数值模拟[J].农业机械学报, 2013, 44(11):168-173.
ZHAO C J, HAN J W, YANG X T, et al.Numerical simulation of temperature field distribution in refrigerated truck based on CFD[J].Transactions of the Chinese Society for Agricultural Machinery, 2013, 44(11):168-173.
[14] 赵时, 姜兆亮, 丁兆磊, 等.梯级送风对冷藏车厢内温度场的影响[J].制冷学报, 2020, 41(5):116-123.
ZHAO S, JIANG Z L, DING Z L, et al.Effect of cascade air supply on temperature distribution in refrigerated truck[J].Journal of Refrigeration, 2020, 41(5):116-123.
[15] TRUJILLO F J, PHAM Q T.A computational fluid dynamic model of the heat and moisture transfer during beef chilling[J].International Journal of Refrigeration, 2006, 29(6):998-1 009.
[16] DELELE M A, KUFFI K D, GEERAERD A, et al.Optimizing precooling of large beef carcasses using a comprehensive computational fluid dynamics model[J].Journal of Food Process Engineering, 2019, 42(4):1-13.
[17] PHAM Q T, TRUJILLO F J, MCPHAIL N.Finite element model for beef chilling using CFD-generated heat transfer coefficients[J].International Journal of Refrigeration, 2008, 32(1):102-113.
[18] MIRADE P S, PICGIRARD L.Improvement of ventilation homogeneity in an industrial batch-type carcass chiller by CFD investigation[J].Food Research International,2006, 39(8):871-881.
[19] KUFFI K D, DEFRAEYE T, NICOLAI B M, et al.CFD modeling of industrial cooling of large beef carcasses[J].International Journal of Refrigeration, 2016, 69:324-339.
[20] 中华人民共和国住房和城乡建设部. GB 51225—2017 牛羊屠宰与分割车间设计规范[S].北京:中国计划出版社, 2017.
Ministry of Housing and Urban-Rural Development of the People′s Republic of China(MOHURD).GB 51225—2017 Code for design of cattle and sheep slaughtering and cutting rooms[S]. Beijing:China Planning Press, 2017.
[21] 朱文颖, 史策, 韩帅, 等.基于CFD的苹果隔板包装预冷温度场研究[J].农业机械学报, 2019, 50(1):331-338.
ZHU W Y, SHI C, HAN S, et al.Investigation on pre-cooling temperature field in partition packaging based on CFD[J].Transactions of the Chinese Society for Agricultural Machinery, 2019, 50(1):331-338.