(1.山东省粮油收储有限公司,济南 250063)
(2.山东建筑大学热能工程学院,济南 250101)
摘 要：多物理场的数值模拟技术是对基于偏微分方程的多物理场模型进行建模和仿真计算一项新技术，多物理场耦合模拟技术不仅可以可以模拟储粮生态系统中的温度、湿度、空气流动等不同物理场的变化规律及相互作用，还可以解释说明其中的生物现象及其演化规律，比如粮仓微环境对储粮发芽率，农药或杀虫剂的衰变率以及昆虫生长率等方面。本文通过多物理场的数值模拟技术在储粮工程中的应用实例,对多物理场的数值模拟技术进行了简要介绍,并对采用多物理场的数值模拟技术在储粮中的通风过程中空气流动、热量水分的迁移和生物现象及其演化规律的研究现状进行了综述,探讨了采用多物理场的数值模拟技术对储粮生态系统研究的优势以及未来发展趋势。(2.山东建筑大学热能工程学院,济南 250101)
关键词：多物理场 数值模拟 粮食储藏 生态系统
传统的储粮生态系统的研究方法主要是实验方法(现场和实验室观测)。实验方法所得到的结果无疑是可靠的，但实验方法需要投入较大的人力物力，成本较高，而且实验结果不具有可重复性。一般来说，储粮生态系统的各种物理现象都不是单独存在的，它包括物理的、化学的和生物的多种变量与粮堆的相互影响和作用，由于生物和非生物因素及其参量之间的相互耦合，储粮生态系统的研究分析要比单独去分析一个物理场要复杂得多。因此，储粮生态系统是一门跨多学科、涉及多物理场应用的领域。因而，储粮生态系统的研究需要一个多物理场分析技术。经过数十年的努力，计算科学的发展以及多学科的交叉为我们提供了更灵巧简洁而又高效的研究手段，使得人们对储粮生态系统的研究和预测成为可能。
多物理场数值模拟技术是国外近年发展起来的一种研究流动、传热、传质及其它物理现象的新方法，它可以形象地再现流动、热湿传递过程及储粮发芽率，农药或杀虫剂的衰变率以及昆虫生长率等方面的情景，因而，在西方发达国家正逐渐被人们所认识和接受。
1 多物理场数值模拟技术原理
多物理场数值模拟技术可以理解为用计算机来做实验，比如某一粮仓内部温度、湿度及粮食水分问题的研究分析，通过建立各种物理现象（温度、湿度和水分）的数学模型(质量守恒方程、动量守恒方程、能量守恒方程)并进行数值模拟计算并显示其计算结果，可以得到粮仓内部各个位置上的基本物理量(如速度、温度、水分等)的分布，以及这些物理量随时间的变化情况 [1]。
采用多物理场数值模拟技术对储粮生态系统进行数值模拟研究，通常包括如下步骤：
(1) 建立反映储粮生态系统演化本质的数学模型。
首先，建立反映储粮生态系统的各个物理量之间关系的微分方程，确定相应的定解条件，这是数值模拟的出发点。储粮生态系统基本控制方程通常包括质量守恒方程、动量守恒方程、能量守恒方程及描述相关生物现象的本构方程，以及这些方程相应的定解条件。
(2) 寻求高效率、高精确度的数值计算方法。
即建立针对控制方程的数值离散化方法，如有限元法。数值计算方法包括微分方程的离散化方法及求解方法，这部分工作包括差分格式、计算网格划分和迭代方法，以及贴体坐标的建立、边界条件的处理等，这些也是多物理场数值模拟技术的核心。
(3)编制程序和进行计算。
这部分工作包括离散后方程组的求解的迭代程序、初始条件和边界条件的输入、控制参数的设定等。
(4)显示计算结果及对结果进行分析。
2 多物理场数值模拟技术在储粮生态系统研究中的应用
多物理场数值模拟技术在粮食储藏中应用的历史并不长，国外从上个世纪八十年代开始应用多物理场数值模拟技术技术对粮食储藏过程中流动、传热、传湿及生物化学等现象进行数值模拟研究，国内则是近几年才开始应用多物理场数值模拟技术对粮食储藏通风过程进行数值模拟分析。
2.1 储粮通风过程中流场、温度和水分场数值模拟
Metzger J F等人[3]1983年提出了针对通风储存小麦的一个模型，并使用基于有限元的方法模拟了强迫对流换热和水分在垂直方向的传递及分布。Chang C S等人[4-5]1993年和1994年提出了一个基于热质平衡理论的严密的数学模型，并采用多物理场数值模拟技术预测了小麦通风储藏过程中温度和水分含量的分布规律，该数值研究是基于有限差分方法，包括代表强迫对流传热传质的源项。Smith E A.[6] 1996年采用有限元的方法对就仓干燥通风时底部具有水平通风道的锥形粮仓内部的速度和压力场进行了模拟分析。SUN DA WEN等人[7-8]1997年提出了一个针对谷物冷却的数学模型，并应用该模型模拟了冬季英国东南部的圆筒仓内谷物冷却过程中温度、水分的变化规律。基于热量和质量守衡定律，2001年Jia Canchun [9]等人模拟了就仓通风过程中粮仓内储藏小麦的温度变化，但是没有模拟小麦水分的变化。Iguaz A.等人[10]2004年等人提出了一个针对周期通风的储藏大米的模型，并采用多物理场数值模拟方法预测了仓内温度的变化。Garg D.和Maier D.E. [11] 2006年使用多物理场数值模拟软件模拟分析了大型圆筒仓中粮粒非均匀分布时通风过程中粮仓内部的流场的分布规律。Daniela de Carvalho Lopes等人[12]等人2006年基于多物理场数值模拟的方法编制一个软件程序（AERO）模拟使用变环境参数做条件的圆筒仓粮食通风过程中的温度的变化规律，而且粮食发热因素也被考虑。Lukasse L.J.S.等人[13]2007年建立了农产品就仓通风时仓内微气候动力学模型，并对土豆储存过程中的温度和水分进行了预测分析。Lukaszuk, J.等人[14]2009年采用数值模拟的方法研究了几种粮粒在不同堆积方式情况下的通风阻力。张忠杰等人[15]2010年采用数值模拟的方法分析了准静态仓储粮堆的温度场。
WANG Yuancheng [16]等人2010年基于多物理场数值模拟技术对就仓降温冷却和干燥通风时粮堆内部的温度和水分的变化过程进行了数值模拟研究并与实验结果进行了比较。研究包括：①根据局部热质平衡原理和吸湿解吸湿理论，建立了深层粮堆内部热湿耦合传递的数学模型而变化的规律进行了数值模拟研究。表1是作者对两种工况下仓储粮堆内部热湿传递过程的数值模拟分析。②采用多物理场数值模拟技术对大型房式仓内粮堆的温度和水分随通风气流温度和湿度改变工况1模拟的是就仓冷却通风过程，数值模拟的粮仓为圆筒仓，如图1所示。粮仓的直径为10m，装粮高度为2.8m。考虑到筒仓的轴对称性，数值计算中采用二维计算模型。粮堆（小麦）初温30℃（303K），干基水分17.65%，（湿基水分15%）。由于模拟的是谷冷机通风降温情况，进风口空气的温度为15℃（288K），进风口空气的相对湿度为56.85%，通风风量为q=14m3/h·ton。图3-6分别为冷却通风过程中粮仓及粮堆中的速度、压力、温度和水分的分布。
工况2模拟的是太阳能/热泵联合就仓干燥通风过程，模拟的实验仓尺寸为1m×2m×0.8m（X﹡Z﹡Y），如图2。小麦的初始温度为20℃，初始水分为17.65%（干基）。小麦的孔隙率为0.4，小麦容重为750 kg/m3，比热容为1871J/（kg·K），导热系数为0.159W/（m·K）。实验仓底部均匀分布的通风孔板，粮堆（小麦）高度为0.3m。模拟的条件为太阳能/热泵联合系统输出的干燥空气，考虑到实际粮仓的围护结构受太阳辐射及仓外空气的影响，在模拟中采用综合温度来考虑太阳辐射及仓外空气的影响。综合温度tz=tw+pl/awo，其中，αw —围护结构外表面与室外空气间的换热系数，W/（m2·K）； tw—室外空气计算温度，℃；τw—围护结构外表面温度，℃；ρ—围护结构外表面对太阳辐射的吸收系数；I—围护结构外表面接受的总的太阳辐射度，W/ m2。综合温度实际上相当于将室外空气温度tw提高了一个由太阳辐射引起的温度的附加值pl/awo[17]。送风口的风速为1.25m/s，送风温度为干燥空气的温度，如图7所示，送风湿度（绝对湿度）为干燥空气的湿度，如图 8所示。图9 为干燥通风12小时粮仓内部温度剖面图，图10 为干燥通风482小时小麦水分剖面图 。
表1 模拟的工况.
图 1. 工况1模拟的物理模型 图 2. 工况2模拟的物理模型
图3 通风时粮仓内部速度分布 图4 通风时粮仓内部压力分布
图5 通风1小时粮堆内温度分布图 图6 通风24小时粮堆内水分分布图
图7 干燥通风时送风温度 图8 干燥通风时送风湿度
图9 通风12小时粮仓内部温度剖面图 图10 通风482小时小麦水分剖面图
2.2 非通风时储粮内部自然对流及温度梯度引起的水分迁移
Abe T.[18]等人基于多物理场数值模拟技术分析了粮仓内谷物温度和水分在储藏过程中变化的规律，研究中未考虑仓内自然对流的作用及影响。PRAKSH Mahesh等人[19-20]1999年对仓储中的吸湿性多孔介质（粮堆）内部的由于自然对流作用引起的热量传递和水分迁移进行了模拟研究，发现了吸湿性多孔介质内部水分由于外界（仓外）大气温度变化而导致局部升高的现象。Jia Can-Chun等人[21-22]于2000年采用数值模拟的方法分别研究了由于粮堆内部发热而引起粮堆温度变化规律以及粮堆温度随仓外气候周期变化而变化的规律。Iguaz A.等人[23] 2004年采用多物理场数值模拟技术模拟分析了仓储粮温和水分随外界气温变化而变化规律。Ali M.S等人[24] 2004年使用二维柱坐标的质扩散方程，借助多物理场数值模拟技术模拟分析了木质圆筒仓内玉米的水分随仓外大气湿度变化而变化的规律，由于忽略了温度梯度对水分扩散的影响，模拟结果低估了仓内玉米的水分。Alabadan B.A.[25]等人2006年、Ruska Laszlo[26]2009年分别采用二维柱坐标的导热微分方程模拟了圆筒仓内谷物温度随仓外气温变化而变化的规律，模拟中忽略了粮粒表面水蒸气的蒸发阻力和吸湿解吸湿的相变潜热的影响。
2.3 储粮生态系统数值模拟（气调、杀虫剂、害虫分布及其他）
THORPE G. R.等人[27-28]于1982年和1997年基于多物理场数值模拟技术分别对通风过程中粮堆内部的谷蠹（Rhyzopertha dominica）分布率以及粮食熏蒸过程中杀虫剂的残留浓度进行了数值模拟研究，并且与实测结果进行了比较，发现数值模拟结果与实验结果基本相符，如图11-12所示。Driscolla R.等人[29]于2000年采用数值模拟的方法预测了粮堆内部害虫发展规律以及发芽率变化情况。S. Mani等人[30]2001年基于多物理场数值模拟技术模拟分析了仓储粮堆中由于害虫诱导发热的温度场。SMITH E. A.等人[31]2001年对CO2气调储粮时粮堆内部的二氧化碳的浓度分布进行了模拟研究，揭示了CO2气调储粮时粮仓内部二氧化碳迁移规律Xanthopoulos G.等人[32]于2003年基于多物理场数值模拟技术对通风储粮时粮堆内部的温度、水分、发芽率和呼吸强度进行了系统的数值分析。
图11 谷蠹的分布率 图12 圆形谷仓中磷化氢的分布
3 多物理场数值模拟技术在粮食储藏中的应用前景
多物理场数值模拟技术是国外近年发展起来的一种研究流动、传热传质及其它物理等现象的新方法，它可以形象地再现流动、热湿传递过程的情景，也为分析和预测粮食储藏、粮仓生态系统的变化规律提供了一个良好的数值分析和优化工具。多物理场数值模拟技术的模拟结果不仅可以提供粮食储藏中空气流动、热湿传递过程中各种热力学参数与流体流动的详细信息，而且还可以为优化储粮通风系统及工艺、探究储粮生态系统的发生发展规律提供了新的方法。随着计算机性能的不断提高和多物理场数值模拟技术 的推广和普及，多物理场数值模拟技术将为更多的仓储技术人员所掌握，并在粮食储藏工程得到更广泛的应用。
参考文献
[1]姚征,陈康民.CFD通用软件综述[J].上海理工大学学报,2002,24(2):137-144
[2]王福军.计算流体动力学分析[M].北京:清华大学出版社,2004
[3]Metzger J F,Muir W E. Computer model of two dimensional conduction and forced convection in stored grain[J].
Canadian Agricultural Engineering,1983,25:119-125
[4]Chang C S,Converse H H,Steele J L. Modeling of temperature of grain during storage with aeration[J].
Transactions of the ASAE,1993,36(2):509-519
[5]Chang C S,Converse H H,Steele J L. Modeling of moisture content of grain during storage with aeration[J].
Transactions of the ASAE,1994,37(6):1891-1898
[6]Smith E A. Pressure and velocity of air during acing and storage of cereal grains[J].
Transport in Porous Media,1996,23(2):197-218
[7]SUN DAWEN, 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.I997,33(4):299-312
[8]Sun DaWen,Woods J.L. SIMULATION OF THE HEAT AND MOISTURE TRANSFER PROCESS DURING DRYING IN DEEP GRAM BEDS[J].
DRYING TECHNOLOGY,1997,15(10):2479-2508
[9]Jia Canchun,Sun Da-Wen,Cao Chongwen.Computer simulation of temperature changes in a wheat storage bin[J].
Journal of Stored Products Research,2001,37(1):165-177
[10]Iguaz A.,Arroqui C.Esnoz A,et al. Modelling and Simulation of Heat Transfer in Stored Rough Rice with Aeration[J].
Biosystems Engineering,2004,89(1):69-77
[11]Garg D.,Maier D.E. Modeling non-uniform airflow distribution in large grain silos using Fluent[C].
Proceedings of the 9th International Working Conference for Stored-Product Protection 754-762.Campinas,SaoPaulo,Brazil,ABRAPOS.2006.
[12]Daniela de Carvalho Lopes, Jose Helvecio Martins, Evandro de Castro Melo,et al. Aeration simulation of stored grain under variable air ambient conditions[J].
Postharvest Biology and Technology,2006,42:115-120
[13]Lukasse L.J.S.,Kramer-Cuppen J.E.,Voor van der A.J. A physical model to predict climate dynamics in ventilated bulk-storage of agricultural produce[J].International Journal of Refrigeration,2007,30:195-204
[14]Lukaszuk,J.,Molenda,M.,Horabik,J.,et al. Variability of pressure drops in grain generated by kernel shape and bedding method[J].
Journal of Stored Products Research, 2009,45:112-118
[15]张忠杰,李琼,杨德勇,等.准静态仓储粮堆温度场的CFD模拟[J].中国粮油学报,2010,25(4):46-50
[16]WANG Yuancheng,DUAN Haifeng,ZHANG Hao,et al. Modeling on Heat and Mass Transfer in Stored Wheat during Forced Cooling Ventilation[J].
Journal of Thermal Science,2010,19(2):167-172
[17]朱颖新.建筑环境学[M].北京:中国建筑工业出版社,2005.
[18]Abe T.,Basunia M.A. Simulation of Temperature and Moisture Changes During Storage of Rough Rice in Cylindrical Bins Owing to Weather Variability[J].
Journal of Agricultural Engineering Research,1996,65:223-233
[19]PRAKSH Mahesh,TURAN Ozden F.,Li Yuguo,et al. A CFD STUDY OF NATIONAL CONVECTION HEAT AND MASS TRANSFER IN RESPIRING HYGROSCOPIC POROUS MEDIA[C].
The Second International Conference CFD in the Minerals and Process Industries,CSIRO,Melbourne, Australia,6-8,December 1999,157-162
[20]PRAKASH M.,TURAN O.F.,THORPE G.R. Natural convection heat and mass transfer in hygroscopic porous media of low permeability[G].
Victoria Univ.Research Report,1999
[21]Jia Canchun,Sun Da-Wen,Cao Chongwen.Mathematical simulation of temperature fields in a stored grain bin due to internal heat generation[J].
Journal of Food Engineering,2000,43(3):227-233
[22]Jia Canchun,Sun Da-Wen,Cao Chongwen.Finite Element Prediction of Transient Temperature Distribution in a Grain Storage Bin[J].
Journal of Agricultural Engineering Research,
2000,76:323-330
[23] Iguaz A.,Arroqui C.Esnoz A,et al. Modelling and Validation of Heat Transfer in Stored Rough Rice without Aeration[J].
Biosystems Engineering,2004,88(4):429-439
[24]Ali M.S.,Al-Amri,Sirelkhatim K. Application of a Mass Transfer Model for Simulation and Prediction of Moisture Distribution in Stored Corn Grains[J].
Scientific Journal of King Faisal University(Basic and Applied Sciences),2004,5:197-213
[25]Alabadan B.A. Temperature Changes in Bulk Stored Maize[J].
Assumption University Journal of Technology,2006,9(3):187-192
[26]Ruska Laszlo,Timar Adrian. Simulation of changes in a wheat storage bin regarding temperature[C].
International Symposia of Risk Factors for Environment and Food Safety & Natural Resources and Sustainable Development,Faculty of Environmental Protection, November 6-7,Oradea,2009,276-285
[27]THORPE G.R.,W.B.ELDER. MODELLING THE EFFECTS OF AERATION ON THE PERSISTENCE OF CHEMICAL PESTICIDES APPLIED TO STORED BULK GRAIN[J].
Journal of Stored Products Research,1982,15:103-114
[28]Thorpe G.R. Modelling ecosystems in ventilated conical bottomed farm grain silos[J].
Ecological Modelling,1997,94:255-286
[29]Driscolla R.,Longsta B.C.,Beckett S.Prediction of insect populations in grain storage[J].
Journal of Stored Products Research,2000,36:131-151
[30]Mani S.,Muir W.E.,Jayas D.S.,et al.Computer modelling of insect-induced hot spots in stored wheat[J].
CANADIAN BIOSYSTEMS ENGINEERING,2001,43:7-14
[31]SMITH E.A.,JAYAS D.S.,ILLE A.DE.Modelling the Flow of Carbon Dioxide Through Beds of Cereal Grains[J].
Transport in Porous Media,2001,44:123-144
[32]Xanthopoulos G.,Woods J.L.Two-dimensional Model of Grain Storage with Dynamic Visualisation: Predictions for Temperature,Moisture Content, Germination and Respiration-a case study for rapeseed[J].AgroThesis,2003,1(1):19-27.
