基于COMSOL软件的深埋地下工程围护结构热湿耦合传递模拟

2020-06-23刘顺波

科学技术与工程 2020年14期

韦明，汪波，刘顺波

(火箭军工程大学作战保障学院，西安 710025)

围护结构是建筑物的重要组成部分，同时也是建筑物热湿负荷的主要来源，对围护结构进行传热传湿计算，是进行地下工程通风空调系统设计的基础和重要依据。相比地面建筑而言，由于岩土的热惰性和较高的含湿量使得围护结构对地下工程内空气温湿度的影响更为显著。针对围护结构的热湿负荷计算，大多集中在地面建筑的单层或复合墙体，研究墙体热、湿和空气传递对室内环境的影响。郭兴国等[1]以含湿率和温度为驱动势建立了多层墙体的一维瞬态热湿耦合传递模型，通过与单纯的导热方程进行对比，表明湿传递引起的潜热会显著改变墙体热负荷，但是在模型中忽略了湿传递对材料本身导热和传湿性能的改变。对于空气渗透性多孔建筑材料，刘向伟等[2]在热湿传递方程中引入空气对流项，建立了墙体内非稳态的热、湿及空气耦合传递模型。Su等[3]、张燎原[4]研究了在不同气候环境下多层墙体的热湿行为，总结了墙体内部的凝结特点和复合墙体在不同气候条件下的受潮程度。为解决复合墙体中每层材料物性参数的不连续性，苏向辉等[5]提出了有限差分逼近方法，避免了边界条件设定为常数带来的计算误差。谭斯鹏[6]在Luikov热湿耦合传递方程的基础上，提出了以水蒸气含量为驱动势的传递模型，解决了实际中在界面处含湿量不连续的问题。

对于地下工程壁面的热湿传递过程，在研究围护结构的热湿负荷时大都将材料的导热率、水蒸气扩散系数等看作常数。张华玲等[7]在热湿传递模型中同时考虑了水蒸气和液态水的传递，指出墙体传热达到稳定状态的时间远小于传湿过程；肖光华[8]、王琴等[9]忽略地下工程传湿过程对传热过程的影响，单独计算地下工程的传热量。王莹莹等[10]、李玮等[11]通过对多种建筑材料进行研究表明湿迁移对传热的影响不可忽略，在混凝土中虽然是潜热热流占比较少，但仍达到8.3%；相比于传湿引起的潜热传递量，Abahri等[12]指出传湿引起的显热传递量相对较少，在模拟计算时可以不加考虑。Chu等[13]虽然全面考虑了地下工程中热湿耦合作用，研究地下工程中两个相互连通区域内的热湿传递过程，没有体现不同含湿量对建筑材料的传热传湿性能的影响。大多数建筑材料均为多孔介质，热传递不仅依靠温度梯度驱动的导热，更多依赖于水分输送引起的显热和潜热变化[14]，在寒冷地区的水蒸气相变传热量能达到总传热量的26.1%，这对于含湿量较高的地下岩土更为明显。同时含湿量的变化会改变建筑材料的有效孔隙率，进而影响材料的水蒸气扩散系数和液态水渗透率。

地下工程围护结构引起的能耗不仅来源于热负荷，湿负荷也占有很大比重。以深埋地下工程中的贴壁式衬砌结构和离壁式衬砌结构为研究对象，在空气-岩石、空气-混凝土和岩石-混凝土3个内部界面上施加第四类边界条件，将建筑材料的传热和传湿特性看作是其含湿量的函数，以温度和相对湿度为驱动势，构建空气-墙体瞬态热湿耦合传递模型。在地下工程内部空气温湿度呈年周期性变化的情况下，利用COMSOL Multiphysics对围护结构和空气之间在两年内的热湿传递过程进行逐时模拟，同时求解温度和相对湿度分布，得到墙体热湿负荷变化规律，为地下工程的通风空调系统设计和优化其运行策略，降低工程能耗提供参考。

1 热湿耦合传递模型

1.1 热湿传递理论

各向同性多孔建筑材料内部存在大量的毛细孔径。假设在代表性单元体积(REV)内材料各项参数相等，在没有空气流动的情况下，多孔材料内部的传湿过程包括分子扩散和液态水流动，其中由浓度差引起的分子扩散符合Fick定律，液态水流动由不平衡的毛细吸力引起，可用Dacy定律描述。传热过程包括温度差引起的热传导和传湿引起的热量传递[15]。水分的存在会显著增强多孔材料的导热能力，同时水蒸气的扩散以焓变的方式引起潜热变化。

在建立传湿控制方程时，常见的驱动势有含湿量、水蒸气分压、毛细压力，对数毛细压力和相对湿度。但是以含湿量作为驱动势时，在不同材料的界面处存在含湿量的跳变；以水蒸气分压作为驱动势时不能反映水蒸气的相变；而以毛细压力作为驱动势，当多孔介质含湿量从零变化完全饱和时，毛细压力的变化范围为0～109Pa，这些都会引起控制方程数值解的不稳定甚至很难达到收敛；以对数毛细压力作为驱动势时在求解控制方程前需要先借助MATLAB求解方程的各项系数[16]。材料的各项性质均可以表示为相对湿度的函数，在热湿传递计算中得到广泛应用。现以温度和相对湿度为驱动势，建立热湿耦合传递模型。

1.2 几何模型

图1所示为离壁式和贴壁式衬砌结构示意图，左侧为厚度d1=0.3 m的混凝土结构，中间空气夹层厚度为d2=0.3 m，对于图1(b)中所示的贴壁式衬砌结构没有空气夹层，左侧的混凝土结构直接被覆于岩石表面，右侧为地下岩土，对于深埋地下工程d3可取为10 m[17]。在所模拟计算的时间内，图中的半无限大几何结构可按一维问题处理。

图1 两种围护结构Fig.1 Two envelopes

1.3 控制方程

1.3.1 多孔介质中的湿传递方程

根据质量守恒定律，多孔介质内任一控制体积内的含湿量变化可表示为

(1)

式(1)中：w为材料的含湿量，kg/m3；φ为相对湿度；t为时间, s;gdiff、gliq分别为水蒸气扩散通量和液态水通量，kg/(m2/s)。根据Fick定律，水蒸气扩散通量可写为

gdiff=δmat(φ)Pv(φ,T)

(2)

式(2)中：δmat为多孔材料的水蒸气渗透系数，s；Pv为水蒸气分压，Pa；T为温度，K；在各个材料的数据库中经常给出的是蒸气阻力系数μ(φ)，其表达式为

(3)

式(3)中：δ0为空气中的水蒸气渗透系数，s；相对湿度定义为水蒸气分压与该温度下饱和水蒸气分压之比，则有

Pv(φ,T)=φPsat(T)

(4)

式(4)中：Psat为只与温度相关的饱和水蒸气分压，Pa，对于不同的温度范围，其近似表达式略不相同，对于本文涉及的温度范围在零度以上，选取文献[18]中推荐的表达式：

(5)

多孔材料中存在着气相压力和液相压力，两者之间的压力差称为毛细压力，毛细压力梯度的作用下，液态水向梯度相反的方向运动，根据Darcy定律，液态水通量可表示为

gliq=-Kl(φ)Pc

(6)

式(6)中：Kl为液态水渗透系数，s；Pc为毛细压力，Pa，反映了多孔介质中气液两项之间的压力关系，根据Kelvin定律，将Pc表示为相对湿度的函数：

Pc=ρwRvTln(φ)

(7)

式(7)中：ρw为水的密度，kg/m3；Rv为水蒸气气体常数，J/(kg·K)。当以液态水中的含湿量作为驱动势时，类比于Fick定律，液态水通量可表示为

gliq=-Dl(φ)w(φ)

(8)

式(8)中：Dl为液态水扩散系数，m2/s。

1.3.2 多孔介质中的传热控制方程

多孔介质中的焓主要存在于固体基质和水分中[19]，根据能量守恒定律和多孔介质中热量的传递形式，控制体积的能量守恒方程表示如下：

(9)

式(9)中：H为控制体积的总焓，J/m3；qcond、qmoist分别为传导热通量和水分传递引起的热通量，W/m2;根据焓与温度的关系，式(9)左边一项可写为

(10)

式(10)中：ρmat为多孔基质密度，kg/m3；Cpmat、Cpw分别为固体基质和水的热容，J/(kg·K)。由Fourier定律可知，传导热通量的方向和温度梯度相反，大小成正比例关系

qcond=-λ(φ)T

(11)

式(11)中：λ为与含湿量相关的导热系数，W/(m·K)。当水分在控制体积间发生传递时，热量也随之转移，水分传递引起的热通量可写为

qmoist=gliqCpwT+gdiff(CpvT+Heva)

(12)

式(12)中：Cpv为水蒸气热容，J/(kg·K)；Heva为蒸发潜热，J/kg。在建筑物理中，湿分传递引起的显热传递远小于潜热的传递量，尤其当多孔材料内相对湿度大于50%的时候，可将上式简化为

qmoist=gdiffHeva

(13)

1.3.3 空气中湿传递方程

水分在空气中只以水蒸气的形式存在，水蒸气的扩散是唯一的传递途径，传湿控制方程可写为

(14)

式(14)中：ρvap为水蒸气密度，kg/m3；gadiff为空气中的水蒸气扩散通量，kg/(m2·s)；假设空气和水蒸气均为理想气体，则由理想气体公式可知：

(15)

空气中的水蒸气扩散通量表达式和式(2)相同，水蒸气渗透系数变为δ0。

1.3.4 空气中热传递方程

与热量在多孔介质中的传递相似，空气中的热量传递方程可写为如下形式：

(16)

式(16)中：ρair为空气密度，kg/m3；Cpair为空气热容，J/(kg·K)；qacond、qmoist分别为空气传导热通量和水蒸气扩散引起的热通量，W/m2，其表达式分别和式(11)、式(13)相似，即

qacond=λ0T

(17)

qamoist=gadiffHeva

(18)

式中：λ0为空气热导率，W/(m·K)；

1.4 定解条件

热湿传递方程在本质上是一组相互耦合的偏微分方程，为获得符合条件的特定解，需要对模型定义初始条件和边界条件。深埋地下工程位于地下恒温层，地表温度波动已不能影响该处温度，温度恒定为当地地表的年平均温度。以西安地区为研究背景，模型的初始温度为15 ℃，衬砌结构初始相对湿度可设为70%。根据文献[20]的描述，地下工程内部温湿度以年为周期波动，夏季最高气温不得高于30 ℃，冬季不得低于10 ℃，相对湿度控制在50%～80%。

(19)

(20)

式中：Tin为工程内部的空气的逐时温度，℃；φin为工程内部逐时相对湿度。混凝土壁面和空气之间的传热量和传湿速率根据第三类边界条件可知：

qs=h(Tin-Ts)

(21)

gs=β(Pvin-Pvs)

(22)

式中：qs为壁面热通量，W/m2；h为表面传热系数，W/(m2·K)；Ts为壁面温度，K；gs为壁面传湿速率；β为壁面传湿系数，s/m；Pvin、Pvs分别为内部和壁面水蒸气分压，Pa，可由式(4)求得。根据ASHRAE 2017手册，表面传热系数和传湿系数分别可设为h=7.7 W/m2、β=2×10-8s/m。由于采用单域方法，将空气夹层和混凝土以及岩石看作是一个整体，因此在衬砌结构内部的空气夹层和岩土以及混凝土的接触界面上不需要额外施加边界条件，默认此处温湿度连续，相关温度梯度和相对湿度连续。

2 模型验证

2.1 偏微分方程在COMSOL中的应用形式

COMSOL主要以求解偏微分方程(partial differential equation,PDE)为基本的仿真形式，在进行仿真求解的过程实质是根据不同的问题特点利用不同求解器逐步对PDE进行离散迭代求解，对于一些高度耦合的非线性问题可以进行特定设置来加强仿真过程中的收敛性。COMSOL针对不同的研究领域，预定义了一些常用模块，模块中的方程形式固定，一个方程只能求解一个变量，可操作性不强。对于本文中推导的热湿传递方程，由于温度和湿度的相互耦合，在一个偏微分方程中需要同时求解温度和相对湿度两个变量，因此需要借助COMSOL中的自定义系数形式对偏微分方程进行求解。自定义系数形式PDE的基本形式为

au=f

(23)

式(23)中，u为要求解的变量，其余参数含义可以参阅文献[21]。为使热湿传递方程符合以上形式，并采用一个求解器同时求解温度和相对湿度，将各个通量表达式分别代入式(1)、式(9)并转化为以温度和相对湿度为因变量，则多孔介质中的热湿传递控制方程可写为

(24)

式(24)中：

dp4=ρmatCpmat+w(φ)Cpw

cp3=Hevaδmat(φ)Psat(T)

空气中控制方程同样可写为

(25)

式(25)中：

da4=ρairCpair+ρvapCpv；

ca1=δ0(φ)Psat(T)；

ca3=Hevaδ0Psat(T)；

2.2 模型验证

为了检验模型在处理实际问题的正确性，需要利用实验数据对模型进行检验。欧洲标准EN15026[22]详细规范了模拟建筑热湿传递的各种细节，并提供了相关实验数据，为验证模型的正确性，利用标准中的实验数据对比模型在相同情况下的模拟结果。标准中实验对象为一可视为半无限大结构的混凝土墙壁，其几何结构如图2所示。

图2 实验中的几何结构Fig.2 The gecmetry of the experiment

墙体宽度为20 m，初始时刻墙体温度为20 ℃，相对湿度为50%。在t=0 s时在内边界上施加温度为30 ℃，相对湿度为95%的第一类边界条件，在外边界上保持初始温湿度不变，上下界面绝热隔湿，热通量和壁面水分通量均为0，因此模拟时将几何结构简化为一维。相关材料参数请参阅此标准。图3所示为在7、30和365 d模型仿真结果和试验数据的对比。

图3 模拟结果和试验结果对比Fig.3 Comparison of simulation andexperiment results

从图3中可以看出，模型计算处的温湿度场和试验结果基本吻合，两者间的最大误差不超过1%，说明构建的模型能正确地模拟建筑热湿传递过程。

3 模型计算及结果分析

将控制方程组式(24)、式(25)输入COMSOL中的PDE接口。由于材料的热物性参数以及与湿分传递有关的参数均为相对湿度和温度的函数，因此需要将方程组中的各项系数以函数的形式提前定义。模型中用到的材料参数以及与含湿量的关系如表1和图4所示，表1中的导热系数为干物质导热系数，表1中数据来源于WUFI材料数据库。按上文所述设置定解条件后模拟衬砌结构在边界条件作用下，两年内温湿度变化情况。

表1 材料热物理性能参数

图4 材料性质随含湿量变化Fig.4 Material properties vary with moisture content

3.1 传热分析

图5所示为壁面处的热流密度随时间的变化关系。在模拟的初始阶段，由于壁面和空气之间较大的温差，壁面处热流变化比较剧烈。在经过20 d的换热后，两种衬砌结构的壁面传热量开始以近似正弦波型变化。在贴壁式衬砌结构中，传热量变化的幅度约为13 W/m2，平均传热量只有3 W/m2，在传热量最大时为9 W/m2，最小时为-4 W/m2。在离壁式衬砌结构中，壁面平均传热量只有1 W/m2，壁面热流密度的变化幅度只有4 W/m2，其中最大时壁面热流密度为3 W/m2，最小时为-1 W/m2。从热流密度的变化幅度中可以反映出，离壁式衬砌结构的壁面温度可以更好地跟随外界温度的变化，使两者之间的温度差保持在一个较小的范围内，进而使壁面的热负荷处于一个较低的水平。

图5 壁面处热流密度Fig.5 Heat flux at the wall surface

随着壁面换热量的周期性变化，围护结构内部温度也大致呈现出周期性变化规律。图6所示为在20、200、380和560 d，两种衬砌结构内部的温度变化情况。从两幅图中可以看出，在一定时间后围护结构内部温度同样具有一定的周期性，并且在20 d时，离壁式衬砌结构的温度已经由初始阶段的不断升高进入这种周期性变化；图6(b)中20 d和380 d所对应的温度有较大的差别，说明，贴壁式衬砌结构温度进入周期性变化的时间要大于20 d。在图6(a)所示的空气夹层中存在较大的温度梯度，岩石温度变化范围在2 ℃之内，而在图6(b)所示的贴壁式衬砌结构中，岩石温度升高较大，最高时上升了7 ℃。这是由于混凝土和岩石的导热系数远大于空气导热系数。由表1可知，混凝土和岩石的导热系数远大于空气的导热系数0.026 W/(m·K)，空气夹层在围护结构传热中能够很好地大幅度地增加热量传递的阻力。因此离壁式衬砌结构有利于减小地下工程的热负荷。

图6 衬砌结构内部温度Fig.6 Temperature inlining structure

3.2 传湿分析

地下工程中室内空气和壁面之间的传热、传湿的变化规律具有相似的特点。图7显示墙体和壁面之间的传湿速率均成周期性变化，其中贴壁式的壁面传湿速率比离壁式略大。

图7 壁面传湿速率Fig.7 Moisture flux on the wall surface

贴壁式衬砌壁面和离壁式衬砌壁面两者最大传湿速率分别为1.5×10-7kg/(m2·s)和1×10-7kg/(m2·s)，相差约50%，两种衬砌结构壁面最小传湿速率分别为-0.9×10-7kg/(m2·s)和-0.8×10-7kg/(m2·s)，相差12.5%，其中贴壁式衬砌结构壁面传湿速率的变化幅度约比离壁式衬砌高出30%。这说明离壁式衬砌结构更能降低混凝土墙体和室内空气之间的传湿速率，减小围护结构的湿负荷。贴壁式衬砌结构中，混凝土墙体内的含湿量受到室内空气相对湿度和岩土含湿量两个因素控制，在离壁式衬砌结构中，混凝土墙体内的含湿量主要受室内空气和夹层内空气影响。图8所示为两种围护结构中，混凝土墙体和岩石中的含湿量变化情况。图8(a)显示，贴壁式衬砌中混凝土的含湿量不断减小，从刚开始的72 kg/m3逐渐降低到69 kg/m3，而离壁式衬砌中混凝土含湿量以一定的频率波动，升高的幅度很小。结合图8(b)可以得出，由于空气夹层的存在，使得混凝土中的水分传递到岩石中。

图8 混凝土和岩石中含湿量Fig.8 Moisture content in concrete and stone

在建筑材料传热传湿中，热量的传递速度远大于水分的传递速度，夹层中空气的相对湿度相比混凝土而言更容易受到温度的影响。图9中展示了两种衬砌结构中相对湿度的变化情况。图9(a)为第380 d，两种衬砌结构中左侧0.6 m范围内的相对湿度分布。结合图5可知，受到温度的影响，夹层内空气相对湿度发生较大改变，并在混凝土-夹层空气界面和夹层空气-岩石界面处和两者进行水分交换，混凝土一侧界面的较高温度使得此处的相对湿度大于岩石一侧界面，因此水分不断从混凝土传递到岩石中，而对于贴壁式衬砌结构，温度对相对湿度的影响只能从混凝土墙体逐步延伸到岩石中。由式(3)可知，水蒸气在空气中传递速率远大于在多孔介质中的传递速率，图9(b) 为两种衬砌结构中混凝土墙体表面处相对湿度和室内空气相对湿度之间的差异。