张晓静1，王远成1，高 帅1，赵会义2，魏 雷2
(1.山东建筑大学热能工程学院，山东 济南 250101；2.国家粮食局科学研究院，北京 100037)
摘 要：小麦储藏过程中，温度和水分是安全储粮的关键因素。迄今为止，相关研究人员设计了多种数学模型来预测通风过程中小麦堆的温度和水分，然而这些模型一般需要占用大量的内存且计算周期很长。基于局部热质平衡原理采用Fortran程序设计了一套新的数学模型，有效地解决了以上问题，并且首次提出R值，使水分传递更加接近实际情况，通过与实验结果的比较，验证了该模型的准确性与可行性，对以后指导安全储粮具有非常重要的意义。
关键词：Fortran程序；储粮通风；传热传质
经过长期的实践和研究，人们发现当小麦堆温度低于 15℃(小麦安全储藏的温度)，小麦堆水分在 12%～12.5%(安全水分)时，可以有效地避免虫害的发生，抑制粮堆中生物体的生命活动，延缓储粮品质的劣变[1]。就仓机械通风在小麦储藏中具有降温效果显著、费用较低等特点，在确保储粮安全方面，发挥着至关重要的作用。为了掌握通风过程中储粮生态系统变化规律，储藏工作者进行了大量的实地测试工作，以研究特定的天气条件下谷物通风系统的性能和空气流速以及粮堆内部温湿度变化规律。但是通风实验需要投入较大的人力物力，成本较高，实验结果还不具有可重复性。为了减少工作量，Bloome和Shove[2]、Thorpson[3]、Alam 和Shove[4]和Sharma[5]开发了仿真模型来描述在可变参数通风下热与质量传递现象。由于模型的限制，这些模型不能模拟谷物在整个通风段的温度。Sun Dawen 等[6-7]提出了一个谷物冷却的数学模型，并应用该模型模拟了冬季英国东南部的圆筒仓内谷物冷却过程中温度、水分的变化规律。基于热量和质量守衡定律，Jia Canchun等[8]模拟了就仓通风过程中粮仓内储藏小麦的温度变化，但是没有模拟小麦堆水分的变化。王远成等[9]对大型房式仓地上笼粮堆温度和水分变化规律进行了数值模拟研究。
以上发展的模型虽然能比较准确地模拟通风过程中温度和水分的变化，但是一般都存在计算周期长，占用计算机内存大的缺陷。本实验通过集总参数的方法结合Fortran程序设计了一套新的数学模型。该模型不针对某一特种粮食，而是普遍适用于各种储粮品种。为了与前人模拟结果进行对比，对储粮过程中小麦堆的温度和水分变化进行模拟。获得在不同通风条件下小麦堆的温度和水分变化，通过与试验结果的比较也证明这个模型的准确性与可行性。
将粮堆划分为相等厚度的谷物薄层，ΔH为空间增量(cm)，Ti为进风温度(℃)，Wi为进风空气湿度比率(kg水/kg干空气)，To为出风温度(℃)，Tw为小麦堆温度(℃)，Wo为出风空气湿度比率(kg水/kg干空气)。这些粮层间的空气流动的方向如图1所示，建立仿真模型的目的是预测在规定的时间内每一层粮层最终的温度和水分含量。从上一层排出的空气状态被用作输入下一层粮层的条件，此过程被重复，直到一段时间的增量上整个谷物柱上的模拟完成。然后该模型以这段时间增量末期的谷物状态为初始条件，开始模拟在进风状态下下一个时间增量的第一层的小麦状态，然后使用和第一次时间增量相同的程序，模拟了在第二个时间增量的整个小麦堆。不断循环直到完成小麦的整个通风过程。
图1 粮层间空气流动方向示意图
2.1 焓湿关系
Brooke[10]给出以下焓湿关系来计算空气的属性：
当干球温度Tdb大于或等于0℃时，饱和蒸气压ps(g/cm2)由下式算得：
ps=exp[58.885 8-12 301.69/(491.69+1.8Tdb)-5.169 23ln(491.19+1.8Tdb)]
(1)
当干球温度Tdb低于0℃时，饱和蒸气压ps(g/cm2)由下式算得：
ps=exp[27.645 2-11 286.148 9/(491.69+1.8Tdb)-0.460 57ln(491.19+1.8Tdb)]
(2)
在给定湿度下的气体压力pv(g/cm2)定义为：
pv=(RH)PS
(3)
其中RH是进风的相对湿度，用小数表示。
空气的湿度比率W(每千克空气所含的水)由下式给出：
W=0.621 9pv/(patm-pv)
(4)
其中是patm是标准大气压，单位是g/cm2。
空气的比体积vsa(m3/kg)(干空气)定义为：
vsa=0.007 573(491.69+1.8Tdb)(1+1.605 5W)
(5)
空气的质量流速ma(kg 干空气/h)由下式给出：
ma=Qa/vsa
(6)
其中Qa是空气的体积流量，单位是m3/h。
空气的比热cpa(kcal/kg℃)由下式给出：
cpa=0.240 5+0.44W
(7)
自由水分的蒸发潜热ls(kcal/kg水)由下式给出：
ls=597.768-0.569 83T
(8)
T是粮堆的温度，单位是℃
空气焓h(kcal/kg干空气)由ASHRAE基本手册(1972)[11]给出为：
h=4.267+0.240Tdb+W(597.386+0.444Tdb)
(9)
2.2 小麦基本属性
Othmer Huang[12]推导出了小麦达到平衡时空气的相对湿度RHe
RHe=exp[(l/ls)ln(2.048 16ps)+c]/2.048 16ps
(10)
l/ls=1+23exp(-0.40M)
(11)
c=-3.34×104M-4.0
(12)
其中l是小麦中水分的蒸发潜热，单位是kcal/kg水；c是常数，M是小麦的水分，单位是干基百分比。
Viravanichai[13]得到了热量的以下关系，硬红春小麦在温度范围为
-33.5℃至21.8℃的小麦的比热cpw(kcal/kg℃)值：
cpw=0.250+0.007 43M,(-33.5℃--21.4℃)
(13)
cpw=0.279+0.005 80M,(-21.4℃--10.8℃)
(14)
cpw=0.297+0.007 08M,(-10.8℃-3.6℃)
(15)
cpw=0.245+0.010 58M,(0.6℃-8.9℃)
(16)
cpw=0.273+0.009 33M,(8.9℃-21.8℃)
(17)
2.3 热湿传递公式
2.3.1 热平衡公式
根据局部热平衡，可以得到以下公式：
空气损失的热+粮堆损失的热+蒸发潜热=0
maΔtcpa(Ti-To)+mgcpw(Tw-To)+(hi-ho)maΔtl/ls=0
(18)
其中ma是空气的质量流速，单位是kg干空气/h；mg是粮层中的小麦质量，单位是kg；hi是进风的焓，单位是kcal/kg干空气；ho是出风的焓，单位是kcal/kg干空气。
对于干燥过程，可变形为：
空气损失的热=粮堆获得的热+蒸发潜热
2.3.2 水分平衡公式
根据质量守恒原理，可以得到以下公式：
空气损失的水分+粮堆吸收的水分=0
(Wi-Wo)100maΔt+(Mo-Mi)md=0
(19)
其中md是粮层中的固体质量，单位是kg；Mo是在时间增量Δt内模拟的小麦堆水分，单位是%干基；Mi是小麦堆的初始水分，单位是%干基；Δt是时间增量，单位是h。
对于干燥过程，可变形为：
空气获得的水分=粮堆失去的水分
本模型以空气逐时的温度和相对湿度、粮堆的深度、粮层的厚度、空气流速、谷物的初始水分和初始温度为初始输入数据，各种基础属性由上述公式计算得到。
进风的初始蒸汽压pvi(g/cm2)，可使用公式1或2和3利用温度Ti算得。类似地，谷物的水分蒸气压pg(g/cm2)，可以利用温度Tw和水分含量Mi，使用公式1或2，3，10，11和12计算得到。
比较pvi，pg，Tw，Ti的关系，存在下面四种可能：
加热干燥小麦 Ti>Tw,Pvi<Pg
冷却加湿小麦 Ti<Tw,Pvi>Pg
冷却干燥小麦 Ti<Tw,Pvi<Pg
加热加湿小麦 Ti>Tw,Pvi>Pg
在以上四个过程中，都达到了热质平衡。当水分传递未达到平衡时，出风空气湿度比率可以通过下面的等式计算得到：
(20)
其中
是空气和小麦在湿度为mi传质温度为Tw或Ti下达到平衡时的空气湿度比率，单位是(kg水/kg干空气)。
R(R以百分比表示)定义为：粮层内实际传递的水分与在平衡状态下传递水分之比。
干燥和调质实验是在温度为6.8℃至25.0℃的范围内进行的，通过反复的实验，R值应为80。在-19.0℃至0℃的温度范围润湿测试使用的R为经验数值。
本实验模拟了120 cm深的粮堆冷却调质和冷却干燥时的温度和湿度。初始条件如表1所示：
表1 用于通风模型的变量
通风实验的实验设备是Sharma[5]在实验室构造出来的如图2。这个设备可以被分成三部分：储存粮食的垂直通风柱，调节空气的空气室，调节穿过粮食的空气的空气分配室。下面是这个设备的简介，更多的细节Sharma[5]已经给出。
通风测试使用的小麦放在由丙烯腈-丁二烯-苯乙烯(ABS)制造的直径为10.16 cm的观察管内。观察管的高度是24.13 cm，粮层的厚度是15.2 cm。沿着观察管设置了6个直径为1.27 cm的孔洞，通过孔洞里抽样检测小麦，沿着观察管轴线分布有热电偶。Brooks玻璃管流量计(型号1110)安装在通风柱和排风支管之间的流动管线上，用于测量空气流速。实验开始前后各从取样室抽取小麦样品，后者在烘干箱中130℃下干燥19 h，水分含量通过下面公式以干重为基础计算：
粮食水分
图2 实验原理图
将利用程序计算出来的结果(程序值)与实验值和Sharma(1973)的模拟结果(模拟值)作比较，得到以下图形。
5.1 工况1模拟结果比较与分析
工况1是在进风参数恒定的状态下进行模拟的，模拟结果如图3、图4。
图3 小麦堆各层温度随时间变化图
由模拟结果可以发现，随着通风时间的增加，小麦堆温度逐渐降低，但是小麦没有冷却到空气的干球温度(图3a)，这是因为在吸附过程中潜热降低。同时由于局部热平衡的作用，随着粮层高度的增加，温度变化越来越缓和(图3b，图3c，图3d)。在通风的前几个小时里，预测温度比测量温度要高，然而在通风的后几个小时里，预测温度比测量温度要低。预测温度和测量温度差距的最大值是0.5℃。通风15 h时，粮层内水分出现一个峰值(图4a)，随着通风时间的增加，峰值升高且峰值后水分含量普遍增加(图4b)，说明当温度达到平衡时，水分尚未达到平衡。
图4 不同时刻下小麦堆各层水分变化图
5.2 工况2模拟结果比较与分析
工况2是在进风参数随时间不断变化的状态下进行模拟的，模拟结果如图5、图6、图7。
图5 不同R值对小麦堆水分含量的影响(80h)
低R值表征水分传递较少，即从空气和小麦中带走少量显热。而高R值表征水分传递较多，即从空气和小麦中带走大量显热。由模拟结果可知最终模拟R值等于80最切合实际(图5)。通风前期，第一粮层的预测温度比测量温度要高。这个差值最大到2℃(图6a)。这是因为当空气流过粮层时，它的湿度趋向于平衡湿度。因此，在较上面的粮层里，空气的湿度与平衡时的湿度差值很小。这就会使得达到平衡时水分传递的实际速率比最大传递速率的80%大，即上层粮层的R值应该比80大，这也是通风前期预测水分低于测量水分(图7)的原因。模型预测的小麦堆温度取决于在吸附过程中空气潜热的增加。在通风的前几个小时内，最高可传送水分的百分比是最低的。因此，R值应该比80小。可能是在通风初期的几个小时里预测温度比测量温度低(图6b，图6c)的原因。
图6 小麦堆各层温度随时间变化图
图7 80 h后小麦堆各层水分变化图
在这些实验中测量和预测水分的最大差值为1.19%。产生差值的原因有可能和以下因素有关：
由于在每个实验中小麦和空气的状态不同，取相同的R值是不准确的。R值的变化预示着水分含量的变化。
粮层的边缘效应。
进风相对湿度和小麦堆水分含量的测量误差。
采用局部热平衡的原理，建立了一套模拟小麦储藏过程中温度和水分的数学模型，借助Fortran程序可以方便地在不同的进风状态下进行模拟，得到不同时刻不同粮层厚度的小麦的温度和水分。由模拟结果得知，储粮通风过程中温度与湿度是相互耦合的。与前人的模型相比，本模型在局部平衡的基础上首次提出了R值，使水分传递更接近于实际情况。通过与测量结果的比较，验证了该模型的准确性与可行性。该模型的突出优点就是占用内存小，运行速度快，操作简单，普适性强，与之前的数学模型相比，有更好的应用价值。
参考文献：
[1]王远成，段海峰，张来林.就仓通风时粮堆内部热湿耦合传递过程的数值预测[J].河南工业大学学报，2009，30(6)：76-77.
[2]Bloome P D，Shove G C.Near equilibrium simulation of shelled corn drying.ASAE，St.Joseph，Michigan，1970:70-325.
[3]Thompson T L.Temporary storage of high moisture shelled corn Using continuous aeration.ASAE，St.Joseph，Michigan，1970:70-825.
[4]Alam A，Shove G C.Simulated drying of soybeans.ASAE，St.Joseph，l，Michigan，1971：71-816.
[5]Sharma S C.Simulation of heat and mass transfer during aeration of wheat and rapeseed bu1ks.Unpublished M.Sc.Thesis，University of Manitoba，Winnipeg，1973.
[6]Sun D W，Woods J L.Deep-bed simulation of the cooling of stored grain with ambient air：a test bed for ventilation control strategies[J].Journal of Stored Products Research.1997，33(4)：299-312
[7]Sun D W，Woods J L.Simulation of heat and moisture transfer process during drying in deep grain beds[J].Drying Technology，1997，15(10)：2479-2508
[8]Jia C C，Sun D W，Cao C W.Computer simulation of temperature changes in a wheat storage bin[J].Journal of Stored Products Research，2001，37(1)：165-177
[9]王远成，张忠杰，吴子丹，等.计算流体力学技术在小麦储藏中的应用[J].中国粮油学报，2008，27(5)：86-90.
[10]Brooker D B.Mathematical Model of the Psychrometric Chart.Trans.ASAE，1967：558-560.
[11]ASHRAE.Handbook of fundamentals.american society of heating，Refrigerating and air-conditioning engineers[M]，New York，N.Y.1972.
[12]Othmer D F，Huang H.Correlating vapour pressure and latent heat data，Ind.Eng.Chem.1940：32，841-846.
[13]Viravanichai S.Effect of moisture content and temperature on specific heat of wheat.Unpublished M.Sc.Thesis，University of manitoba，Winnipeg.1971:52.
ZHANG Xiao-jing1,WANG Yuan-cheng1,GAO Shuai1, ZHAO Hui-yi2,WEI Lei2
(1.School of Thermal Energy Engineering , Shandong Jianzhu University, Jinan Shandong 250101; 2.Academy of State Administration of Grain，Beijing 100037)
Abstract：The temperature and moisture content are critical to achieve safe wheat storage.The researchers have designed a variety of mathematical models to predict the temperature and moisture content of wheat during ventilation, which generally take a lot of memory and time.The problem can be solved effectively by a new mathematical model based on local heat and mass balance principle and Fortran program.And the proposing of R makes moister transmission closer to the actual condition.The accuracy and feasibility of the model was verified by comparison with the experimental results and it is very significant for the guidance of wheat storage.
Key words：Fortran program; ventilation for stored grain; heat and mass transfer
收稿日期：2015-09-10
基金项目：国家自然基金(51276102)、国家粮食公益专项(201313001，2015449-001-03)
作者简介：张晓静，1990年出生，硕士研究生.
通信作者：王远成，1963年出生，教授，博士.
中图分类号：S 379.5
文献标识码：A
文章编号：1007-7561(2016)02-0085-05