关于 COMSOLMultiphysics垃圾填埋场温压耦合数值模拟探索

作者；张毅

  目前对垃圾填埋场内温度场的时空变化特征研究主要采用现场监测的方法，国外在这方面的研究较多，而国内的研究却有所欠缺。对填埋气体数值模拟研究，前人亦做了很多工作。但是，前人对垃圾填埋场的研究大都以单个物理场为主，较少将填埋场内温度场与气体压力场综合考虑，尤其是关于温度场及热量的研究相当欠缺，在热量的产生释放和运移传输等领域仍然存在着一定的争议。

  本文基于传热理论和流体力学原理，尝试对填埋场温度场和气体压力场进行耦合数值模拟研究，既可以极大地节省生产成本与时间，又能预测填埋场内温度和气压的变化情况，从而为填埋场的设计、运营、改进及场区内各种资源能源的回收利用提供充足可靠的热学与流体力学性能指标，旨在有效地提高填埋场的生产、回收和利用效率，预防与治理其产生的安全隐患问题。

1温度一气体压力耦合数学模型

1.1模型假设条件

  垃圾填埋场中热量伴随填埋气体（CH4、CO2等温室气体为主）而生，并随着气体和渗滤液的传输而发生传递，故温度的变化过程实质上是一个温度场与气体运移场的耦合问题，建立耦合方程时引进以下假设：

  (1)考虑到整个填埋场尺寸远远大于填埋气体的分子单元，模拟时将垃圾土等效成各向同性的均质介质材料，便于确定相关参数。

  (2)实际上，填埋气体运移过程是气体和液体等多相流动过程的一个部分，不同相态的物质流动时必定存在一定影响。考虑到气体的运移速率比其他相态的要快很多，并且渗滤液的运移状态相对填埋气体要稳定得多，因此在模拟时可把多相流动的问题简化为填埋气体单相流的问题。

  (3)填埋气体为理想气体，其运移规律遵循达西定律。

  (4)填埋气体视为不可压缩气体。

1.2温度场控制方程

  从垃圾土中选取任意一个平行六面体微元，热流量矢量分解为沿x、y和z3个坐标轴方向的热流分量。如图l所示，Φx、Φy、Φz分别代表这3个方向的分量。

根据傅里叶定律和能量守恒定律，可得：

  式(1)是笛卡尔坐标下填埋场的三维瞬态热传导微分方程，p、c、Φ、λ均为变量。

  填埋场内热量分布受气体流动的影响，将此项考虑后，方程变为：

  式(2)中，T为温度(K)，t为时间(s)，c为填埋体的热容(kJ/(m3．K))A为垃圾体的热传导系数(J/(ms·K))，cg为气体热容(kj/(m3.K)),Vg为气体流速(m/s)，Φ是热源项。

  垃圾填埋场内，温度的上升主要由于垃圾降解放热引起。

热源项，即垃圾降解放热速率：

  式(3)中，Φ代表放热速率( W/m3)；ΔHb为有机物分解产生的热量( J/mol)；Rb为产气速率(mol/( m3·s))。

  垃圾被填埋直至产气稳定的过程，气体呈非稳态，气体流速满足达西定律：

  式(4)中，k为垃圾土的导气渗透系数(m/s)；μs为填埋气体的动力粘度(Pa-S)。

  式(5)即为考虑气体迁移后的温度场数值方程。

1.3填埋气体迁移方程

  根据多孔介质的流体动力学原理，在假设填埋场内介质变形很小及液态水的运移比气体稳定多时，填埋气体的连续性方程可写作如下形式：

  式(6)中，SG为垃圾土体的饱气率；ne为垃圾体的有效孔隙度；p为填埋气体的密度（kg/m3）；f为时间(d)；V为相对介质骨架的达西流速（md）；g为内部热源项，即气体产生速率( d-l)。

将填埋气体看作理想气体，有：

  式(7)中，Wm为摩尔质量(kg/mol)；p为压力(Pa)；R为理想气体常数(m3.Pa/(mol-K));T为温度(K);达西流速V为：

  式(8)中，K为导气系数张量(m/d)，含气率OG的函数(OG-SGne)；φ是流体的势能(m)；k是导气率张量( m2)。

导气系数写成如下经验公式：

  式(9)中，Ks是饱气状态下的导气系数；m是与介质特性相关的参数。

  综合上式，有：

2广州兴丰垃圾场温度一气体压力耦合数值模拟

2.1广州兴丰垃圾场概况

  广州兴丰垃圾场位于白云区太和镇兴丰村南侧，填埋区面积47.5 hm2，可填埋处理垃圾2 560万t。

   场地地形条件复杂，地貌单元属低山丘陵，场地山地一般坡度为(35±5)0，沿山脊平缓处5°～l0°，陡

坎、隘口处700~800。场地及邻近山地高程在60～200 m，东侧的菠箩山顶为场地最高点，高程为192.88m，场地最南侧水塘为最低点，高程为73.74 m，高差为119.14 m。

  场地地层主要为第四系覆盖层和下卧混合花岗岩My3。第四系覆盖层分布整个场地，主要为坡积土和残积土，厚度为0.50～31.70 m，平均渗透系数为

2.67xl0-5～1.69xl0-3 cm／s，透水性较好。

  场地地下水来源主要为大气降水及地表水，按赋存方式分为第四系孔隙水和基岩裂隙水。孔隙水主要分布在沟谷的冲积层中，受大气降水影响变化明显。基岩裂隙水主要分布在岩石强风化带和岩石中等风化带中。地表水主要为大气降水，降雨后，以山脊山坳为分水岭。地表水大部分汇集到沟谷区，部分地表水渗入地下形成地下水，由高处向低处，以渗流或泉水形式向谷底排泄。

2.2概念模型

  物理现象可通过偏微分方程进行描述，上文得到的温度一气压耦合模型表明填埋场内部热量与气体之间存在非线性、强耦合性及时变性的关系。热量伴随气体而生，气体在扩散与迁移过程会携带大量热量，从而改变填埋场内部的温度分布；温度的重新分布又会影响温度梯度，在不同的梯度下气体的运动方向与轨迹会发生变化。

  然而，非线性强耦合偏微分方程组解析解的求解十分困难，实际工程应用中多采用数值解。本文采用的COMSOL Multiphysics大型有限元数值仿真软件内嵌了大多数物理场偏微分方程，并且支持用户自定义物理场的偏微分方程或方程组，从而实现多物理场的直接强耦合分析。

  本研究中，填埋气体被看成不可压缩的的流体，运移规律遵循达西定律，使用软件内嵌的达西模块来描述(如式(14)所示）；将垃圾体视为各向同性的多孔介质，热量伴随气体而生，热量的传播采用内嵌的带热源项的多孔介质热传导模块(如式(13)所示）来描述。COMSOL在进行数值求解时会根据

这3项建立起2个偏微分方程之间的强耦合关系，只要把模型参数、边界条件和初始条件根据实际状况设置好，便能求解出该方程组的数值解。

2.2.1几何模型

  根据相关地形地质资料，结合实地调研考察，采用COMSOL Multiphysics有限元软件建立一个二维平面几何模型，一般而言，渗滤液及地下水的水头位置约在距底部1/3H处，由于水的比热较大，水体所覆盖的区域温度变化较小，可以看作是恒温区，且区域内部水相处于饱和状态，填埋气体很难渗入，因此认为场区内地下水位以下的区域是一个非渗气的恒温

区，计算时可以忽略不计，故本文把渗滤液及地下水的水位线看作是填埋场的底边界。

  模型具体尺寸见图2，填埋场底部长L=150 m，填埋总体高度H=120 m，有效高度H’=2/3H=80 m，周边斜坡自然坡度比约1:2，每填埋20 m高的垃圾加铺一层Im厚的覆盖层，以防雨水下渗造成渗滤液污染地下水以及垃圾填埋气体逃逸到大气环境中。

2.2.2边界条件和初始条件

  分别对压力边界和温度边界进行设置。填埋区左右两侧边界为零通量边界，对温度、压力皆适用。压力边界：填埋场顶部设置为恒压边界，压力值为一个大气压( 101 325 Pa)；底部是渗滤液及地下水水头位置，填埋气体无法通过，故气体的通量为零，为非渗透边界。温度边界：顶部为与大气接触，主要受地表气温的影响，计算时可取广州兴丰区的年平均值22.5℃，底部边界为水头位置，垃圾体与渗滤液长期接触，考虑到长期受渗滤液及地下水滞热效应的影响，并结合常年实测资料，温度取值为37℃。垃圾体内部初始压力设置为一个大气压值(101 325 Pa)，内部初始温度设为年平均值为22.5℃。

2.2.3有限元模型

  图3所示为填埋场有限元模型，采用自由网格方法划分创建网络，单元类型为拉格朗日二次单元，网格形状为三角形，网格自由度数目为49 766个，网格点数6 304个，网格单元数12 276个，边界单元数1 536个，端点单元数18个，最小单元质量0.726，单元面积比0.013。

2.2.4数值计算模型参数的取值与验证

  各种材料的热力学及流体动力学参数的取值参照岩土工程勘察报告、国内学者关于兴丰垃圾场的研究报告、垃圾产气产热资料，具体的数值计算参数如表1所示。

  为了验证数学模型的适用性与正确性，对兴丰场表层垃圾体进行了为期4个月的监测。为了更加贴切地反映真实场景，验算时必须首先根据垃圾土表层实测温度数据确定模型上边界的温度条件，如图4所示为2012年10月30日至2013年3月2日的表层温度监测数据，数据在监测期间波动性较大，采取高斯( GAUSS)函数进行拟合计算（如图蓝线所示），取其平均值18.41℃作为验算模型的顶部温度边界条件。

  图5为2012年10月30日至2013年3月2日期间垃圾土表层以下Im深度处温度监测值与模拟值的对比图，图中实测拟合曲线与模拟曲线之间的最大相对误差率均小于15%，由此可见参数和模型的选取具有一定的合理性。

3温度一气体耦合场模拟结果分析

3.1  温度场时空分布特征

  图6所示为兴丰场封场后50、100以及150 a内部温度分布状况。由图可知，150 a后的温度值最低，其次依次为100、50 a，反映出垃圾降解能力逐渐减弱以及热量不断往外界散失的过程。水平方向上，温度分布均匀，温度值大致相等，仅仅在左右两侧不透气不透热边界附近出现了较小范围的波动，这主要与垃圾土的均布性、左右两侧边界的对称性以及地形坡度等原因有关。垂直方向上，垃圾体顶部温度相对较低，基本维持在295～300 K(22～

27℃)之间，主要由于与大气接触，受环境温度的影响较大；垃圾体中部区域(10～70 m)的温度较高，随着时间的推移，“高温区”逐渐往下向场区底部转移，直至150 a整个垃圾体的温度分布自上而下呈现递增的规律，此时底部附近区域温度约为310~315 K(37～42℃)之间，表明场区内温度场已渐渐接近稳定状态。

  箭头方向指示垃圾体内总热通量（对流热通量和扩散热通量，辐射热通量占少部分，在垃圾土内可以忽略不计）的运动方向，由于两侧为绝缘边界，热流通量主要集中在垂直方向上进行传播。

  图7(a)~(c)是垃圾体内一定深度处温度随时间变化的趋势图。由图可知，封场后数年内，温度迅猛上升，在5～7 a的时候达至峰值温度，随后的3～4 a内维持在峰值附近，之后缓慢下降，直至最后热量收支达至平衡，温度稳定，温度下降速率远小于上升速率。下降速率随着时间的推移逐渐变小。其后揭示了这样一个过程：垃圾被填埋后便会发生降解，一开始，以好氧反应为主，垃圾土产生大量的热量，温度显著上升，随着氧气逐渐消耗，好氧反应减弱，转为厌氧反应，厌氧反应相对于好氧反应是一个低放热的过程，会持续数十甚至数百年，温度在此期间变化甚微。

  图7中2条曲线分别代表考虑填埋气体耦合的温度变化状况及没有考虑耦合的情况，可见考虑耦合后的温度比未考虑耦合的情况要高，差值在4～5 K之间。垃圾降解阶段产生的热量与产气量是相互影响、相互制约的。气体的产生与迁移是一个伴随着热量变化的过程，而热量的传递反过来会影响气体的释放传输。

  为了进一步了解垃圾体内温度在垂向上的分布特征及规律，特作封场不同时间内温度随深度变化的分布图以供分析。如图8所示，在温度场趋向稳定化的过程中，温度在垂向上的分布特征呈现“中间高，两头低”的规律，温度自顶部沿深度的增加会有显著上升，尔后在中部区域(30～60 m)相对较稳定，往下逐渐减小，在底部渗滤液水头位置附近达至恒定值310

K(37℃)左右。这是由于大量热量在顶底部散失，中部区域垃圾稳定放热而且填埋气体较平稳运移造成的。此外，图中所示同一深度处的温度随着填埋年龄的增加而降低，再次反映了垃圾中可降解有机物的生物化学反应正逐渐衰弱，而热量却不断往外散失的过程。

3.2气体压力场时空分布特征

  垃圾填埋后，微生物有氧呼吸开始，产生各种高温气体，填埋场垃圾体的温度上升，由于覆盖层或新鲜垃圾隔绝了氧气的来源，氧气消耗完毕后，微生物的厌氧呼吸便开始，此过程生成大量甲烷和二氧化碳（温室气体），并伴有热量的产生。此后，气体会携带着部分热量通过扩散和对流的方式在垃圾土内运动，重新影响温度场的变化，而温度的重分布又影响微生物降解的速率与活性，反过来作用于产气量。因此，气体与热量是密不可分的，下面具体研究填埋气体的分布规律。

  填埋场内的气压分布主要与产气量及气体的运动有关，而温度是影响产气量的重要因素，通过影响微生物的降解活性从而控制产气量；填埋气体的运动方式以对流和扩散为主，对流是以压力梯度为动力，而扩散则是以浓度梯度为动力。图9(a)~(c)分别反映了封场50 a、100 a及150 a后填埋场气压分布状况，如图所示，水平向上的分布特征与温度场的类似，气压值大致相等，只在两侧边坡附近出现了轻微的波动；垂直方向上，以覆盖层为界，形成4个颜色深浅不同的“气压带”，每个“气压带”内的气压值大致相等，说明覆盖层对气体流动的阻挡作用较明显，气体产生后基本被局限在2个覆盖层所封闭的区域中，形成如图所示的色带；而不同“气压带”之间，气压随着深度的增加而增加，最大值出现在底部附近，约为104 kPa左右，最小值在顶部附近，约为一个大气压，主要是因为底部为不渗透边界，气体产生后聚集在底部扩散不开，只能向顶部往大气中逃逸。

  图9中箭头方向指示填埋气体速度流动的方向，由图可知，由于填埋场底部和边坡两侧均为绝缘边界，气体自产生后，便沿着垃圾中的孔隙向填埋场顶部流动，并从顶部处逃逸至大气中，此过程中会携带部分热量，从而影响温度场，致使其重新分布。

  图10反映了-15、-30、-50和-70 m处填埋气体在封场后压力随时间的变化情况。由图可知，垃圾填埋后1—2 a内，气压显著上升至峰值压力104 kPa左右，之后维持在峰值压力附近较短时间（约1～2 a左右），最后缓慢下降，下降速率远小于上升速率，反映出如下一个过程：垃圾一经填埋，微生物便会利用被封闭的有限氧气进行好氧呼吸，产生大量高温气体（C02为主），形成骤增的高压，这个时期一般不超过1a，随后氧气消耗殆尽，微生物进入厌氧呼吸阶段，产生大量甲烷气体，随着时间的推移，厌氧呼吸逐渐减弱，产气量明显减弱，加之覆盖层的阻挡作用，气压在随后的数十甚至上百年内便会缓慢衰减。

  图11中的曲线分别代表封场后不同时间气压沿深度的变化规律。由图可知，气压总体上会随着深度的增加而上升，但其过程并不是线性连续的，在覆盖土层的位置会骤增，增幅最大者可达3 500 Pa，最小则为200 Pa，增幅自顶部往底部逐渐减小，表明覆盖层对填埋气体的流动具有明显的阻挡作用，填埋气体到达覆盖层后不能有效地运动，部分气体便会积聚起来，形成覆盖层上下部垃圾土之间的“气压断层”，越往上接近顶部，气体量越少，其逃逸速度越大，这种差异性会越明显。对比不同填埋年龄段的垃圾体，可以发现气压值随着填埋时间的增加而减小，差值在2 000—4 000 Pa之间，表明垃圾体内部有机质成分的降解反应正在逐渐衰减，填埋气体的产生量越来越少，并不断往外界散失，造成气压值有所降低。

4结论

  本文以多孔介质热力学原理及气体渗流理论为指导，给出了垃圾场温度一气压耦合的数学模型，并以广州兴丰垃圾填埋场为研究背景，结合野外监测资料，采用COMSOL Multiphysics有限元软件对场区内部温度场与气体压力场进行了数值模拟分析，得到以下主要结论。

  (1)场区内一定深度处的温度随时间变化规律为：在封场后5～7 a内快速上升至峰值并保持3~4 a左右，随后缓慢减少直至稳定。封场后某一时刻场区内的温度在空间上变化规律为：水平方向上温度基本一致，只在左右边界附近发生较小的波动；热量和温度主要沿垂直方向传播扩散，在稳定化的过程中中部区域的温度一般高于顶底部，直至最后达到稳定状态，此时温度自顶部往下部逐渐增大。

(2)场区内一定深度处的压力随时间变化规律为：在封场后1～2 a内快速上升至峰值并保持1～2 a左右，最后缓慢下降，下降速率远小于增加速率。封场后某一时刻场区内的压力随时间变化规律为：水平方向上的分布规律与温度场相当；气体主要沿垂直方向传播，整体上气压自顶部往底部减小，但由于覆盖层的阻挡作用，气压会在覆盖层的位置发生突变，形成“阶坎”，每两个覆盖层所包围的区域内气压变化不明显，可认为基本不变。

4摘要：温度与气体压力是衡量垃圾填埋场稳定性的2个重要指标，进行温度一气体压力耦合研究有利于揭示二者在填埋场内部的时空演化规律。以广州兴丰垃圾填埋场为研究背景，推导了多孔介质温度一气体压力耦合数学模型。采用COMSOL Multiphysics有限元模拟软件，结合监测数据，检验了数学模型的适用性和合理性，预测了封场后至稳定化过程中兴丰垃圾填埋场温度与气体压力的时空动态演化特征。研究结果表明：场区内一定深度处的温度在封场后5~7 a内快速上升至峰值并保持3~4 a左右，随后缓慢减少直至稳定，呈现“中间低，两头高”的规律；封场后某一时刻场区内的温度水平方向上基本一致，只在左右边界附近发生较小的波动；场区内一定深度处的压力在封场后1~2 a内快速上升至峰值并保持1~2 a左右，最后缓慢下降；某一时刻场区内的压力在水平方向上的分布规律与温度场相当；气体主要沿垂直方向传播，整体上气压自顶部往底部减小。分析成果可为填埋场设计、施工、运营及管理提供指导性意见。