一种粮食床通风过程中温度和水分变化的模拟分析方法

摘要

本发明公开了一种粮食床通风过程中温度和水分变化的模拟分析方法，所述方法包括：基于粮堆焓平衡方程和水分平衡方程，采用离散化方法将粮食床分成N个节段，采用粮食平衡水分方程计算第1个节段的粮层含水率，然后采用迎风差分格式将温度和水分分配到第2、3、4、······N‑1、N节点的过程。能够合理估计变化的粮食温度和含水率，通风空气的温度和湿度通常在通风过程连续变化，本发明的方法模拟通风粮堆的热量和水分迁移，给出热量和质量平衡方程的逐步推导，就一维系统假定均一流速的气流通过粮食床，热量也仅在气流方向流动，既可阐明解决问题的原理，又可保证软件程序中简单而整体的代数式。

说明书

技术领域

本发明涉及食品科学技术领域，具体涉及一种粮食床通风过程中温度和水分变化的模拟分析方法。

背景技术

目前，国家加大了对粮食储运领域应用基础研究，一些学者开始探索粮堆的温度和水分迁移过程，当前国内公开的试验结果主要是采用商业计算流体动力学(CFD)软件进行模拟。该软件包含了连续性方程、动量守恒方程、对流传热方程和水分迁移方程，其运行计算机语言环境是Fluent软件中的Fortran语言,该软件非常昂贵，一般科技工作者和粮库无法利用。

探索粮堆温度和水分迁移过程的算法，意义重大。粮堆生态系统中生物与非生物成分之间相互作用，是动态的，每个成分持续地影响其它成分。通风的作用就是运动合适质量的空气穿过粮堆以调节储粮存在的条件。粮食储藏作业的目标就是保持粮食的品质，减少昆虫、螨类、霉菌危害粮食的可能性。为了成功实现这个目标，就要调控储粮生态系统微生境，造成粮食籽粒内、粮堆(空间或籽粒之间)空隙内的物理和化学条件对有害生物是敌对的，而有益于粮食品质保持。在实践应用方面，国家粮食局2002年发布了《储粮机械通风技术规程》(LS/T 1202-2002)，提出了允许降温通风的温度条件和湿度条件。在这个储粮技术规程及我国粮食储藏学科的专著中，对通风过程缺乏描述。粮堆降温通风作业往往伴随降水，温度迁移前沿区带早于水分迁移前沿区带，国内对粮堆降温通风期间温度和水分迁移的前沿区带和尾巴区带缺乏研究。与国外粮食品种单一、储存时间短的储粮比较，考虑我国的粮仓仓型、储粮区域、粮食种类品种、储存时间、粮食含水率及杂质含量等因素到通风过程模拟中，有利于节约能耗，降低粮堆水分损失，达到优粮优储的目标。

粮堆物理特性的知识提供了粮堆通风作业的原理基础。粮堆的物理特性是内在联系的，并被其它的物理、生物因素所影响。粮堆颗粒间隙的条件包括空气的湿度和温度、有毒气体如熏蒸剂、接触式杀虫剂。颗粒间隙的条件强烈影响粮堆籽粒的条件。粮食颗粒之间的空气温度和湿度是最关键的非生物因子变量，通过通风作业过程调控。粮堆内湿度和干球温度对应于籽粒间隙空气的湿球温度，湿球温度与昆虫群体生长的速率直接相关。在粮堆中霉菌倾向于如何发展，也是粮堆温度和湿度的函数。

发明内容

为此，本发明提供一种粮食床通风过程中温度和水分变化的模拟分析方法。

为了实现上述目的，本发明提供如下技术方案：

本发明提供一种粮食床通风过程中温度和水分变化的模拟分析方法，所述方法包括：基于粮堆焓平衡方程和水分平衡方程，采用离散化方法将粮食床分成N个节段，采用粮食平衡水分方程计算第1个节段的粮层含水率，然后采用迎风差分格式将温度和水分分配到第2、3、4、······N-1、N节点的过程。

本发明的一个实施例中，所述离散化方法为：在通风气流方向将粮食床分为N-1层，粮食床的高度是L，每个层长度是Δx，则方程为

Δx＝L/(N-1)

所述第1个节段的粮层含水率的计算方程为：

其中，M是平衡含水率，小数干基表示，ERH是平衡相对湿度，小数表示；系数A、B、C与之对应的国产粮食种类的如表1所示。

本发明的一个实施例中，所述粮堆焓平衡方程为：

式中，M，粮食含水率，干基小数表示；T，粮堆温度，℃；w，粮粒间隙空气的湿度比率，kg水/kg空气；x,一维方向的距离，m；t，通风时间，hr；H

C

本发明的一个实施例中，所述粮堆水分平衡方程为:

式中，μ

将f

式中，假定水蒸汽和空气具有相同的速率。

本发明的一个实施例中，所述粮食床水分平衡方程通过下述质量方程解释水分通过粮食床任何附加的运输可能来自水蒸气的扩散：

式中，D

本发明的一个实施例中，所述粮层水分平衡方程进一步处理为：

式中，

本发明的一个实施例中，所述迎风差分格式为：

式中，w

本发明的一个实施例中，所述方法包括采用中心差分格式计算粮食床湿度比率、温度和水分速率变化的二阶导数。

本发明的一个实施例中，所述方法还包括通过算法中对高含水率玉米粮堆涉及储藏真菌引起的呼吸热。

本发明具有如下优点：

本发明的粮食床通风过程中温度和水分变化的模拟分析方法，能够合理估计变化的粮食温度和含水率，通风空气的温度和湿度通常在通风过程连续变化，本发明的方法模拟通风粮堆的热量和水分迁移，给出热量和质量平衡方程的逐步推导，就一维系统假定均一流速的气流通过粮食床，热量也仅在气流方向流动，既可阐明解决问题的原理，又可保证软件程序中简单而整体的代数式。

本发明的模拟分析方法试验证明：采用湿度比率0.013kg/kg和20℃空气降温通风35℃和含水率13.5％的稻谷粮堆，稻谷粮堆平均温度和含水率在通风150小时分别是23.57℃和12.42％，在通风200小时分别是21.09℃和12.25％，在通风300小时分别是20.03℃和12.49％。

采用湿度比率0.013kg/kg和20℃空气通风35℃和含水率12.0％的大豆粮堆，粮堆平均温度和含水率在通风150小时分别是23.22℃和11.30％，在通风200小时分别是20.89℃和11.20％，在通风300小时分别是20.02℃和11.48％。

附图说明

为了更清楚地说明本发明的实施方式或现有技术中的技术方案，下面将对实施方式或现有技术描述中所需要使用的附图作简单地介绍。显而易见地，下面描述中的附图仅仅是示例性的，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图引伸获得其它的实施附图。

本说明书所绘示的结构、比例、大小等，均仅用以配合说明书所揭示的内容，以供熟悉此技术的人士了解与阅读，并非用以限定本发明可实施的限定条件，故不具技术上的实质意义，任何结构的修饰、比例关系的改变或大小的调整，在不影响本发明所能产生的功效及所能达成的目的下，均应仍落在本发明所揭示的技术内容得能涵盖的范围内。

图1为本发明粮食床的离散化方法的结构示意图；

图2为本发明的稻谷和大豆粮食床降温通风期间各节段温度变化；

图3为本发明的稻谷和大豆粮食床降温通风期间各节段含水率变化。

具体实施方式

以下由特定的具体实施例说明本发明的实施方式，熟悉此技术的人士可由本说明书所揭露的内容轻易地了解本发明的其他优点及功效，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例1、通风粮堆的焓平衡方程和水分平衡方程的建立

粮堆通风过程包含了水蒸气速率(μ

式中ρ

通过一维粮堆具有横切面积1m

式中，μ

就焓平衡方程，包含粮食温度(T)、粮食含水率(M)和粮粒间隙空气湿度比率(w)，表达为：

式中，M，粮食含水率，干基小数表示；T，粮堆温度，℃；w，粮粒间隙空气的湿度比率，kg水/kg空气；x,一维方向的距离，m；t，通风时间，h；H

C

粮堆水分平衡方程为:

代入f

推导方程5，假定水蒸气和空气具有相同的速率。但是Thorpe在1995年指出，水分通过粮食床任何附加的运输可能来自水蒸汽的扩散，能够被加入一个额外的术语质量守恒方程来解释：

式中，D

数学方法描述的难度是，通风空气的温度和湿度随着时间发生任意变化。获得通风粮堆温度和粮食含水率响应的方法是利用数值解分析，解主导的偏微分方程。我们采用Kreyszig在1998年提出的简化数值解方法，用于解方程3和6。

本实施例中，粮堆和粮食床为相同客体。本发明的方法在通风粮堆控制热和水分迁移的方程数值解，如图1所示，粮食床的离散化方法，在通风气流方向粮食床分为N-1粮层，即具有N个节段，床的高度是L(m)，因此每个粮层的长度是Δx(m)，则方程为：

节点被步长距离Δx分开，总计N个节段，在第1个节段，通风所用的气流的温度和湿度比率是已知的，它们对应进入粮堆的空气温度和湿度比率。假定空气和粮食之间质量和热平衡，方程8能够计算第1个节段的粮层含水率。

式中，A、B和C为系数。发明人发现，现有技术中系数A、B和C的值不能够覆盖国内的同一种类粮食品种多的问题，且所述A、B和C的值与所述粮堆粮食的种类、吸着类型(吸附、解吸或吸着平均)有关。为此，在本公开的一个优选示例中，发明人采用静态称重的平衡水分测定方法，选定国产的17个品种的稻谷、14个品种的小麦、16个品种的玉米和10个品种的大豆，在10～35℃的温度范围内和11％～96％平衡相对湿度范围内，经过大量的测试，测定其相应的平衡水分等温线；进一步采用方程8，通过非线性回归拟合获得系数A、B和C的值。并给出了国外部分粮食种类的A、B和C的值，以便用于进口粮食的粮堆水分迁移和温度迁移的分析。

表1粮食种类的方程(8)的系数

在通风开始点，粮食床内粮食含水率和温度的分布(它们在整个粮食床内是不均匀的)是已知的。因此，温度和含水率被分配到节点2、3、4、······N-1、N。下一步是计算每个节点粮食温度和粮食含水率随时间变化。为了完成这个目标，对方程3和6出现的一阶和二阶导数作简单近似处理。例如，

式中w

在

式中，

把它们代入方程10，结果是：

通过依次类推，温度随迁移距离x变化的导数可以近似化为：

正在通风的粮食床内，粮食含水率和温度从点到点是连续变化的，可从方程3和方程6中可以捕捉到。按照上述处理，含水率和温度这些变量仅在节点计算。换句话说，正在通风粮食床的条件仅在不连续点计算，这个问题已经被离散化。当Δx变得越来越小，数值解通常成为连续方程3和6的更好近似值。

粮食的含水率和温度也随时间连续变化，在这个方法中，它们仅在一定的不连续的时间间隔被评价。当通风系统打开之后，如果粮食的含水率和温度在时间间隔Δt之后被计算，这些值在时刻Δt、2Δt、3Δt、···、(p-1)Δt、pΔt、(p+1)Δt···被更新。这里p表示第p时刻的更新。仅在一定时刻计算变量的这个思路，用于离散化方程3和6的时间导数，这样，

式中

解方程3和6，通过把方程9、13、14、15代入它们，给出离散形式。

水分守恒方程6接着变换为：

解方程18最简单的方法，是明确地更新粮食含水率，这样，

在解程序的开始，值

在内部节点的粮食温度在下一个时刻步态之后相似地计算，同时粮食含水率被更新，即，

式中，C

这个解的方法编成Visual Basic程序给出。需要注意的是，编写的计算机程序和上述方程之间的主要差异是，程序中不包含

实施例2、计算粮粒间隙空气的湿度比率H(i)及粮食含水率M(i)

(1)粮食水分汽化的潜热(h

h

(2)粮食籽粒间隙空气的相对平衡湿度ERH、湿度比率H(i)及湿化积分热H

P＝P

式中M

(3)计算粮食含水率M(i)

实施例3、呼吸作用对通风粮食床的热和质量传递的影响

焓守恒方程为：

式中，S

按照化学成分组成，S

式中，S

水分守恒方程变成：

水分守恒方程的一维形式是：

式中，S

呼吸热通过对试验数据拟合数学方程获得真菌呼吸热的计算方法。呼吸作用被认为是碳水化合物完全燃烧形成CO

dm＝8.83×10

方程41用于15.5℃、含水率25％(湿基)、破碎率30％的脱粒玉米。由于粮食温度和含水率随粮仓内位点发生变化，需要将真实时刻t值作图变为生理时刻t

式中，M

(1)温度修正因子M

当T≤15.5℃和M≤19％时，M

当T＞15.5℃和19＜M＜28％时，

M

当T＞15.5℃和M＞28％时，

M

(2)水分修正因子M

式中M

方程41到46用于计算生理时刻t

如图2所示，采用本发明的算法分析国产稻谷与大豆粮堆降温通风期间的温度和水分迁移变化，分析结果如图2和图3所示。

表2给出了国产稻谷和大豆粮堆降温通风期间，粮堆各节段降低到20℃的时数及对应含水率。如表3所示，在粮堆含水率保持不变的情况下降温通风，稻谷和大豆粮堆由35℃降低到20℃，分别需要123和126小时。

如图2和图3所示，在连续通风条件下，在降温通风前100小时，稻谷粮堆前4层已经降低到20℃，粮堆平均温度和含水率分别是27.08℃和12.84％。降温通风到150小时，稻谷粮堆前7层已经降低到20℃，粮堆平均温度和含水率分别是23.57℃和12.42％。降温通风到200小时，稻谷粮堆前11层已经降低到20℃，粮堆平均温度和含水率分别是21.09℃和12.25％。降温通风到250小时，稻谷粮堆前15层已经降低到20℃，粮堆平均温度和含水率分别是20.20℃和12.33％。降温通风到300小时，稻谷粮堆第21层已经降低到20℃，粮堆平均温度和含水率分别是20.03℃和12.49％。

表2国产粮食降温通风各层(节段)降低到20℃的时数及对应含水率

表3国产粮食降温通风期间风机运转时数预测

在连续通风条件下，在降温通风前100小时，大豆粮堆前4层已经降低到20℃，粮堆平均温度和含水率分别是26.81℃和11.64％。降温通风到150小时，大豆粮堆前8层已经降低到20℃，粮堆平均温度和含水率分别是23.22℃和11.30％。降温通风到200小时，大豆粮堆前11层已经降低到20℃，粮堆平均温度和含水率分别是20.89℃和11.20％。降温通风到250小时，大豆粮堆前15层已经降低到20℃，粮堆平均温度和含水率分别是20.14℃和11.31％。降温通风到300小时，大豆粮堆第21层已经降低到20℃，粮堆平均温度和含水率分别是20.02℃和11.48％。

本发明是基于CFD的粮食床通风过程中温度和水分变化的模拟分析方法，能够合理估计变化粮食温度和含水率，通风空气的温度和湿度通常在通风过程连续变化，本发明的方法模拟通风粮堆的热量和水分迁移，给出热量和质量平衡方程的逐步推导，就一维系统假定均一流速的气流通过粮食床，热量也仅在气流方向流动，既可阐明解决问题的原理，又可保证软件程序中简单而整体的代数式。

虽然，上文中已经用一般性说明及具体实施例对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。