基于有限元模拟的咖啡粉充填密度对咖啡胶囊内冲泡流场的影响

李园锦 潘嘹 钱静 杨硕

摘 要：目的：研究咖啡粉孔隙率对咖啡胶囊内冲泡流场及其达到稳态时间的影响，寻找咖啡胶囊内最佳咖啡粉充填密度。方法：基于自由流动和多孔介质流动理论，应用Comsol Multiphysics软件耦合Navier-Stokes方程和Brinkman方程，对咖啡胶囊不同孔隙率下胶囊内冲泡流场进行有限元模拟，研究流场分布规律。利用示踪法计算冲泡流场达到稳态的平衡时间，并构建咖啡粉不同孔隙率对平衡时间的经验公式。采用相同方法模拟手冲咖啡机中的冲泡流场，通过比较意式咖啡机出口流量和模拟出口流量验证方法准确性。结果：意式咖啡机出口流量和模拟出口流量基本一致，误差小于6.10%，表明本研究使用的方法能够用于表征咖啡粉冲泡流场分析。对咖啡胶囊内冲泡流场的模拟结果显示，咖啡胶囊孔隙率越小，咖啡胶囊稳态冲泡流场整体流速越小，且流场达到稳态时间越长。在注水孔、通孔膜小孔和出口附近流速和速度梯度显著增大，靠近杯壁方向流速减小。构建的咖啡粉孔隙率对平衡时间的经验公式能够准确预测不同孔隙率下咖啡粉中流场的平衡时间。结论：通过Comsol软件耦合Navier-Stokes方程与Brinkman方程能够用于模拟咖啡冲泡过程中的流场分布，通过咖啡孔隙率与平衡时间的经验公式有助于确定咖啡胶囊冲泡风味最佳的孔隙率，为咖啡胶囊风味改进中咖啡粉充填密度的选择提供科学依据。

关键词：咖啡胶囊;孔隙率;充填密度;有限元模拟

随着经济高速增长和生活质量提升，咖啡为越来越多消费者所喜爱[1]。在消费升级背景下，普通的速溶咖啡已经不能满足消费者对咖啡品质的追求，而味道纯正的手冲咖啡需要专业人员花费大量时间操作才能得到，无法适应快节奏生活下消费者快速冲泡的需求。咖啡胶囊是一种结合了速溶咖啡与手冲咖啡两者优点的新型咖啡产品，不仅能得到味道纯正的咖啡，而且咖啡胶囊机的操作过程方便快捷[2]。咖啡胶囊的冲泡品质受到咖啡豆处理工艺、咖啡粉充填密度和冲泡工艺等因素的多重影响，因此需要研究影响咖啡胶囊冲泡品质的重要参数并进行优化。

近十几年来，国内外很多学者针对影响手冲咖啡风味的因素开展了大量研究，比如胡荣锁等[3]、何余勤等[4]通过电子鼻对感官风味定性表达，研究了研磨时间与烘焙程度对咖啡风味物质的影响。Natnicha Bhumiratana等[5]通过量化15种香味感官属性，发现相较于咖啡豆种类，烘焙程度和研磨工艺对咖啡风味影响更大。虞健等[6]研究表明，不同烘焙程度主要对咖啡酸味、辛辣味和刺激味感官属性有影响。总体上来说，这些研究工作主要集中在咖啡豆处理工艺对咖啡风味的影响上[7-10]。手冲咖啡的风味依赖于咖啡师的经验与操作，而咖啡胶囊只能通过提前设定机器中各种参数来控制咖啡的品质，目前有关咖啡胶囊重要参数对咖啡品质的影响缺乏研究。咖啡粉的充填密度决定了其孔隙率的大小，而孔隙率决定了冲泡介质（即热水）对咖啡风味物质的萃取程度，进而影响咖啡的冲泡品质。咖啡粉过实孔隙率过小，则热水通过咖啡粉难度变大，导致萃取时间长而使咖啡味酸苦，咖啡粉过松孔隙率过大，热水将快速流过咖啡，导致萃取不足而使咖啡味寡淡。

流体在多孔介质内流动情况十分复杂，难以通过试验发现其内部规律，而基于流体力学理论建立合适方程组采用数值计算方法，可以有效研究多孔介质流动规律。王冬平[11]基于达西定律耦合N-S方程和Brinkman方程，建立非线性渗流耦合模型描述了地下水突水全过程，并分析了突水过程中水压和流速的分布规律，找到诱发突水的关键部位与实际情况一致。杨贝贝等[12]通过构建导数形式的波动方程，求解了饱和多孔介质动力问题，并发现饱和多孔地基中空沟隔振效果随孔隙率和泊松比增大而变差。本研究针对咖啡胶囊不同孔隙率下冲泡流場进行数值模拟，通过将咖啡粉均质化为多孔介质，分析自由流动和多孔介质流动理论，应用Comsol软件耦合Navier-Stokes方程与Brinkman方程模拟出咖啡胶囊稳态冲泡流场，分析稳态时水在咖啡胶囊内的速度与分布规律，在稳态结果基础上，利用示踪法对冲泡流场达到稳态的平衡时间进行模拟计算，并拟合孔隙率与平衡时间函数，为提升咖啡胶囊风味中孔隙率的选择提供科学依据。

1 理论依据

雷诺数计算公式如式（1）：

式（1）中，ρ表示流体密度、υ表示特征速度、d表示特征长度尺寸、μ表示流体动力黏度。计算本文研究的咖啡胶囊内流体流动的雷诺数为171，所以咖啡胶囊内流体流动状态为层流。

针对咖啡胶囊内冲泡流场的数值模拟，有如下基本假设：（1）流体为连续介质[13];（2）在咖啡粉中流体渗流为饱和多孔介质渗流。咖啡胶囊内流体流动包含自由流动和多孔介质流动两种，在咖啡粉中为多孔介质流动，其余部分均为自由流动，所以采用耦合Navier-Stokes方程与Brinkman方程的模拟方法[14]。整个冲泡过程流场温度几乎不变，模型为等温流动模型。

咖啡胶囊内水的自由流动规律选用Navier-Stokes方程描述，方程如式（2）、（3）：

式（2）、（3）中，ρ为流体密度、p为压强、u为速度矢量、μ为动力黏度、F为体积力源项、g为重力加速度、I为单位矩阵。

咖啡胶囊内多孔介质流动规律选用Brinkman方程描述，方程如式（4）、（5）：

式（4）、（5）中，p为压强、u为速度矢量、μ为动力黏度、K为多孔介质渗透系数、F为体积力项。

2 数值模拟

Comsol Multiphysics软件（Comsol）具有强大的多物理场仿真计算能力，以其独特优势被广泛应用于流体力学的计算当中。Comsol数值模拟计算的流程为选择物理场模块、建立几何模型、设置边界条件与参数、划分网格、利用求解器求解、结果后处理。本研究采用Comsol进行数值模拟，采用Solidworks建立三维模型。

2.1 三维模型建立

Comsol内置LiveLink接口，可以直接导入Solidworks中建立的三维模型。本研究以鼎加弘思公司咖啡胶囊为研究对象建立三维模型。模拟方法常用模拟实验来验证其可行性与准确性，用意式咖啡机进行模拟实验，以意式咖啡机粉碗为对象建立模拟实验模型。

2.1.1 咖啡胶囊三维模型 咖啡胶囊由高阻隔杯体、封口膜、通孔膜、咖啡粉、滤纸、泄压导流结构件组成。图1咖啡胶囊整体高度为49 mm、封口膜直径为48 mm、注水孔直径为3 mm、通孔膜小孔直径为1 mm、相邻小孔间距10 mm、泄压导流结构件出口通孔直径为1 mm。忽略冲泡过程中咖啡粉底部滤纸对流场的影响。

2.1.2 模拟实验模型 图2所示意式咖啡机粉碗整体高度为31 mm、碗口直径为52 mm，滤网上筛孔直径为0.5 mm、出口直径为5 mm。

2.2 边界条件设置

边界条件设置的准确性对咖啡胶囊内冲泡流场的数值模拟十分重要。以实际咖啡胶囊机冲泡时注水压力作为入口边界条件，将实际情况与三维模型相结合，设置边界条件见表1。

2.3 参数设置

选择9组不同咖啡粉孔隙率进行数值模拟，分别编号为C1～C9，对应渗透率见表2。咖啡胶囊机冲泡温度为92 ℃，查表得此温度下水的密度为963.96 kg/m3，水的动力黏度为0.000 31 Pa·s。

2.4 网格划分

网格划分对于模拟计算的收敛性和准确性有较大影响。由于咖啡胶囊通孔膜和泄压导流件结构复杂，所以采用非结构化网格，并在自由流动和多孔介质流动边界面进行网格加密。整个咖啡胶囊三维模型最终网格数量为73万。

3 模拟方法验证

为了证明本研究所应用的耦合Navier-Stokes方程与Brinkman方程的模拟方法可以用于咖啡胶囊内冲泡流场的数值模拟及其准确性，选用一个工况来进行模拟实验，在实验室进行咖啡冲泡实验，以出口流量为判断依据，比较实验与数值模拟两者所得稳态出口流量结果。

3.1 材料与仪器

材料：咖啡粉（埃塞俄比亚风味，鼎加弘思公司，真密度为1.001 kg/m3）、怡宝矿泉水。仪器：意式半自动咖啡机（DL-KF6001，广东东菱电器有限公司）、电子天平（ML6001E，梅特勒仪器上海有限公司）。

3.2 实验方案

利用意式半自动咖啡机对咖啡粉进行加压冲泡，实验咖啡粉孔隙率为68%，冲泡入口压力为251 kPa恒压。实验步骤如下：

（1）取下咖啡机把手及粉碗。将咖啡机开关调至萃取档，预热咖啡机。（2）称量14.9 g咖啡粉装入粉碗，用填压器压实咖啡粉饼至指定体积。（3）安好粉碗，锁上把手。确认工作指示灯亮起，机器已可以使用，按下萃取按钮，开始冲泡。（4）待咖啡稳定流出后，用烧杯接住并计时20 s。（5）选择合适量筒测量烧杯中咖啡体积。孔隙率计算公式为式（6）：

式（6）中，ρ*为多孔介质的表观密度、ρs为多孔介质的真密度。

4 结果与讨论

4.1 模拟方法验证结果

实验测得当稳定流出咖啡时，平均出口流量为3.524 mL/s、数值模拟计算出口流量为3.309 mL/s、误差范围为6.10%。该实验证明，在Comsol软件中Navier-Stokes方程与Brinkman方程耦合的模拟方法可以用于咖啡胶囊冲泡流场的数值模拟，并且在适当参数值与边界条件设置下模拟精度较高。在数值模擬计算中，假设咖啡粉分布是均匀的，而实际中咖啡粉的分布略有差别，所以计算会产生误差。

4.2 数值模拟稳态流场结果分析

为了方便观察数据及分析稳态流场速度的分布规律，选取模型中心截面图作为研究对象。由图3可以看出，孔隙率对咖啡胶囊模型中流场的影响十分显著，入口压力一定时，随着孔隙率的减小，流场整体速度随之减小，但流速分布规律趋于一致，在多孔介质区域流体速度较低且变化较小，在注水孔、通孔膜小孔和出口附近流速和速度梯度显著增大，这是由于流体在流经通孔膜小孔和流向出口时，流通横截面积突然减小，所以流速会极速增大[16]。当水经过通孔膜后，越靠近杯壁，整体流速越小。

4.3 数值模拟平衡时间结果分析

选用示踪法计算咖啡胶囊冲泡流场平衡时间。在初始阶段，出口浓度保持为0，此时还没有水从出口流出，随后出口浓度开始迅速增加，到了一定时间后出口浓度缓慢增长，此时咖啡胶囊内流场近似达到稳态，因此以出口浓度/入口浓度达到99%作为流场达到稳态的判断标准。图4展示了不同孔

图4 不同孔隙率下流场达到稳态的平衡时间隙率下出口浓度与入口浓度之比随时间变化的曲线，反映了咖啡胶囊内流场随时间逐渐达到稳态的时间长短，可以得到孔隙率φ=70%、68%、65%、62%、59%、55%、50%、43%、38%时，咖啡胶囊内流场达到稳态的平衡时间分别为18.8、36、126、164、297、382、892、9 720、98 520 s。当入口压力一定时，随着咖啡粉孔隙率的减小，渗透系数会变小，水渗流通过咖啡粉的难度加大，流场达到稳态所需的平衡时间就越长。

根据图4曲线趋势，构建平衡时间关于孔隙率的经验公式。利用Matlab软件拟合模型参数，结果如式（7）所示：

式（7）中，y为咖啡胶囊冲泡流场达到稳态的平衡时间（s）、x为咖啡粉孔隙率（%）。

由图5可知，当入口压力一定时，随着咖啡粉孔隙率的减小，咖啡胶囊内冲泡流场达到稳态的平衡时间逐渐增大，而且在孔隙率小于60%时，平衡时间加速增长。平衡时间过长，会造成萃取过度，导致咖啡苦味过重。

5 结论

本研究耦合Navier-Stokes方程与Brinkman方程对咖啡胶囊不同孔隙率下的冲泡流场进行了模拟，理论模型和实验结果吻合良好。模拟结果表明，孔隙率对咖啡胶囊内冲泡流场影响显著，当入口压力一定时，随着孔隙率的减小，流场整体速度随之减小;在注水孔、通孔膜小孔和出口附近流速和速度梯度显著增大，靠近杯壁方向流速减小;咖啡粉孔隙率与咖啡胶囊冲泡流场达到稳态的平衡时间两者之间存在分式型函数关系，当入口压力一定时，平衡时间随孔隙率的减小而增加。本文研究结果为咖啡胶囊机的设计和工艺参数的选择提供了科学依据。

参考文献

[1]饶建平.即饮咖啡的国内外市场现状和发展趋势[J].食品工业科技，2017，38（23）：346-351.

[2]唐晓双，刘飞，汪才华，等.胶囊式咖啡概述及其安全性分析[J].食品安全质量检测学报，2013，4（3）：961-965.

[3]胡荣锁，谷风林，宗迎，等.不同研磨时间对咖啡感官风味的影响[J].中国农学通报，2012，28（9）：259-263.

[4]何余勤，胡榮锁，张海德，等.基于电子鼻技术检测不同焙烤程度咖啡的特征性香气[J].农业工程学报，2015，31（18）：247-255.

[5]Natnicha Bhumiratana，Koushik Adhikari.Evolution of sensory aroma attributes from coffee beans to brewed coffee[J]. Food Science and Technology，2011，44（10）：2185-2192.

[6]虞健.不同烘焙程度咖啡感官特征及主要化学成分分析[D].江苏无锡：江南大学，2014.

[7]翟晓娜，杨剀舟，柴智，等.电子鼻电子舌在咖啡风味研究中的应用[J].食品工业科技，2016，37（5）：365-370.

[8]吕文佳，刘云，杨剀舟，等.咖啡主要烘焙风味物质的形成及变化规律[J].食品工业科技，2015，36（3）：394-400.

[9]Czerny M，Grosch W.Potent odorants of raw arabica coffee.their changes during roasting[J]. Journal of Agricultural and Food Chemistry，2000，48（3）：868-872.

[10]Prakash M，Ravi R，Sarvamangala G K，et al.Sensory profiling and product positioning of roasted and ground （brew）coffee and soluble （instant）coffee with and without added flavor[J]. Journal of Sensory Studies，2000，15（1）：101-117.

[11]王冬平.采动诱发隐伏陷落柱突水机制与渗流耦合模型研究[J].煤炭工程，2016，48（7）：93-96.

[12]杨贝贝，李伟华，赵成刚，等.基于Comsol的饱和多孔介质动力方程的数值模拟及应用[J].地震工程学报，2017，39（2）：321-328.

[13]李亨，张锡文，何枫.论多孔介质中流体流动问题的数值模拟方法[J].石油大学学报（自然科学版），2000（5）：111-116.

[14]杨天鸿，陈仕阔，朱万成，等.矿井岩体破坏突水机制及非线性渗流模型初探[J].岩石力学与工程学报，2008（7）：1411-1416.

[15]刘正.低渗透煤层微观孔隙结构研究及应用[D].西安：西安科技大学，2013.

[16]黄腾瑶，周晓华，胡延伟，等.基于Comsol的煤层注水压力对湿润半径的影响研究[J].矿业安全与环保，2018，45（4）：49-53.