引言室内外温湿度和气压的差异会使建筑围护结构产生不同的热工性能。水蒸气的存在会导致围护结构隔热性能显著下降[1]。追求围护结构高保温性能的同时，会采取不同的材料组合实现墙体的功能，由此导致霉菌滋生[2]、结露、耐久性下降等与水分相关的变化[3]。气候发生变化后，墙体的热湿传递趋势甚至会增加10%以上[3]。地域与多孔材料不同的情况下，很难找到一个通用的热湿耦合（Heat, Air and Moisture, HAM）传递模型，大量的HAM模型被不断地提出。关于墙体的热湿耦合传递的研究基于Philip[4]等以及Luikov[5]等提出的多孔介质的热质传递耦合理论与唯象的方法，但是其相关控制方程的参数难以获取，且部分参数需要经验系数的确定。Talukdar[6]等的工作为热湿耦合的研究提供了大量的实验数据，常被用以验证模型的正确性。在此基础上，许多学者对多层墙体的HAM耦合模型和实验及工具进行了大量的报道。如Ferroukhi[7]等结合Trnsys与COMSOL软件的联合仿真评估多层墙体对建筑整体热湿度性能的影响。但是大量相关工作主要集中在建筑构件与外部环境的热湿耦合传递或以特隆布墙为主的墙体内对流热湿传递的研究。较少的模型关注建筑围护结构和空气共同作为一个系统的热湿耦合相互作用。给出多孔介质和空气域中的HAM耦合方程的一般形式，采用经典假设得到可供COMSOL软件直接求解的PDE形式。通过与已有软件接口的结果进行对比，量化分析空气层在多层墙体结构中对热湿传递的影响。1控制方程建立由水、空气、干燥多孔材料组成的多相热力学系统，其内部的热质转移可以利用偏微分方程的守恒形式表示。参考常用的建筑领域的经典假设[8-10]：各相间保持局部热平衡与局部质量平衡；内部不考虑碳化等化学反应；将多孔介质近似为连续介质，各物理属性采用体积平均获得。1.1水分输送方程1.1.1多孔介质内的水分输送方程以相对湿度为驱动势，对墙体构建HAM耦合迁移模型。借鉴机械公式[11]，多孔介质的水分输送方程由水蒸气和液态水的质量守恒方程共同组成。∂wφ∂t+ρgug⋅∇ωv+∇⋅gdiff+ul⋅∇ρl+∇⋅glc=0 (1)wφ=εpslρl+εpsgρgωv (2)式中：wφ——总湿度，%；ρg——水蒸气密度，kg/m3；ug——由外部压力梯度引起的水蒸气的对流速度，m/s；ωv——水蒸气质量分数，%；gdiff——水蒸气的扩散通量，kg/(m2·s)；ul——由外部压力梯度引起的液态水的对流速度，m/s；ρl——液态水密度，kg/m3；glc——液态水的毛细管输送，kg/(m2·s)；εp——孔隙率，%；sl——液态水饱和度，%；sg——水蒸气饱和度，%；ρg——湿空气密度，kg/m3。ωv=Mvφcsat(T)ρg=cv·MvPART-cv·Ma+cv·Mv (3)cv=PvRT (4)Pv=φ·Psat(φ,T) (5)式中：Mv——水蒸气摩尔质量，取0.018 kg/mol；csat(T)——饱和蒸汽压力下的蒸汽浓度，mol/m3；cv——蒸汽浓度，mol/m3；PA——绝对压力，取101.325 kPa；Ma——干空气的摩尔质量，取0.028 97 kg/mol；Pv——水蒸气分压力，Pa。根据Fick定律，水蒸气的扩散通量gdiff为：gdiff=-δmat(φ)∇Pv(T,φ) (6)Pv(T,φ)=φPsat(T) (7)δmat=26.1×10-6μRvT1-w/wsat(1-ε)1-w/Wsat+ε (8)式中：Pv(T,φ)——水蒸气分压力，Pa；δmat——材料的水蒸气渗透率[12-13]，多孔介质孔隙率为1时，采用空气-蒸汽渗透率，δ0=1.96e-10s[8]；μ——动力黏度，kg/(m·s)；Rv——水蒸气气体常数，取461.52 J/(kg·K)；w——绝对湿度；wsat——绝对压力下的绝对湿度；ε——多孔介质孔隙率，%；Wsat——饱和压力下的水分含量，kg/m3。为了简化计算，多孔介质的渗透率通过蒸汽阻力因子μfactor(φ)计算。由外部压力梯度引起的水蒸气和液态水的对流速度ug和ul利用达西定律计算。μfactor(φ)=δ0δmat(φ) (9)液态水的毛细管输送glc为：glc=-Dw∂wφ∂φ∇ϕw (10)Dw(θ)=D0enθ (11)式中：Dw——液态水扩散系数；θ——归一化含水量；n、D0——实验拟合常数，参考Lockington[14]等与Hall[15]等给出详细的获取过程。将各项构造为PDE方程形式：ea∂2φ∂t2+da∂φ∂t+∇⋅(-c∇φ-αφ+γ)+β∇φ+aφ=f (12)式中：ea——质量系数，取0；da——阻尼系数；c——扩散系数；α——对流通量系数；γ——守恒通量源，取0；β——对流系数；a——吸收系数；f——源项。da=∂w(φ)∂φ (13)c=δmat(T,PA)Psat(φ,T)+Dw(φ)∂w(φ)∂φ (14)α=δmat(T,PA)dPsat(φ,T)dT∇T (15)β=ρw_airωvuw_air (16)a=Mv∇T⋅uw_airR1TdPsat(φ,T)dT-Psat(φ,T)T2 (17)f=ρlul (18)式中：ρw_air——空气中的水分密度，kg/m3；uw_air——空气中的水分动力黏度，Pa·s。1.1.2空气中的湿控制方程空气中的热湿传递主要由空气热对流和蒸汽扩散组成，多孔介质中存在的与液态水相关的传递现象不需要考虑，其偏微分方程的系数为：ea=0；da=Mvcsat(T)ρg ；c=δmat(T, PA)Psat(φ,T) ；α=δmat(T,PA)dPsat(φ,T)dT∇T；γ=0；β=Mvcsat(T)ρguw_air；a=Mv∇T⋅uw_airR1TdPsat(φ,T)dT-Psat(φ,T)T2；f=0。1.2能量守恒方程1.2.1多孔介质内的能量守恒方程非饱和多孔介质内的单元体能量守恒方程为[16]：∂(ρe)∂t+divρVe+divjq+Q=0 (19)等式左侧依次表示能量变化率、对流项、扩散项、源项。由于墙体为低渗透率材料，内部水分和基体材料的温度满足局部热平衡[17]，因此对流体和固体可以采用单温度方法计算。为了保证模拟计算的准确性，采用混合物定律定义热属性并重写式(18)的各项，得到：dρCpeff∂T∂t+∇⋅qcond+∇⋅qconv+∇⋅qmoist+dQ=0 (20)式中：Cp——材料比热容，J/(K·kg)；qcond、qmoist、qconv——热传导、水分所携带的潜热和显热和干燥空气的热对流，W/m2。将各项替换为以温度为因变量的系数型偏微分方程形式：ea'∂2T∂t2+da'∂T∂t+∇⋅-c'∇T-α'T+γ'+β'∇T+a'T=f' (21)da'=1-εpρmatcmat+εpsl·wφcw+1-slρaircpair1-εpρmat  (22)c'=keff(φ)+(Lv+cpvapT)δmat(T,PA)φdPsat(φ,T)dT (23)α'=0γ'=-(Lv+cpvapT)δmat(T,PA)Psat(T)∇φ (24)β'=ρaircpair +ρw_airωvLv+cpvapTφPsatφ,TdPsatφ,TdT-Lv+cpvapTρw_airωvTuair+ρwcpwuw+cpwDwφ∂wφ∂φ∇φ (25)f'=-LvωvPsat(T)RTu⋅∇φ (26)1.2.2空气中的热量传递在空气域内，不考虑液态水的能量输送，不同的偏微分方程的系数为：da'=ρaircpair；c'=keff(φ)+(Lv+cpvapT)δ0φdPsat(φ,T)dT；γ'=-(Lv+cpvapT)δ0Psat(T)∇φ；β '=ρaircpair+ρw_airωv×Lv+cpvapTφPsatφ,TdPsatφ,TdT-Lv+cpvapTρωvTuair；a'=0；f '=-LvωvPsat(T)RTu⋅∇φ。1.3速度场1.3.1多孔介质内的速度场建筑围护结构中产生的压差是空气传输的驱动势。驱动压力由蒸汽分压力、烟囱效应、器械通风共同作用。多孔介质内的连续性方程为[8]：ε∂ρair∂t=∇·ρairu (27)式中：u——达西速度，由Dupuis-Forcheimer关系得u=εv，其中v为固有速度，m/s。根据不可压缩流的假设，连续性方程简化为：∂∂xjρavj=0 (28)多孔介质内的动量方程通常类比空气内的N-S方程[18]。在建筑领域，速度场的方程描述通常简化为体积力项、压力梯度、达西阻力项构成的动量方程，通常利用达西定律表示。u=-kmatμair∇P-ρairTg (29)式中：kmat——多孔介质渗透率，m2；μair——空气动力黏度，(Pa·s)；g——重力加速度。空气密度可以利用理想气体定律表示：ρair(T)=P0MairRT (30)式中：P0——平衡压力，Pa；Mair——空气的摩尔质量，kg/mol；R——通用气体常数。结合式(26)和式(28)可以得到泊肃叶定律，常用描述多孔内部势流模型的雷诺方程表示：∇δmat∇P*=0 (31)δmat=ρkmatμair (32)P*=P-ρairTgh (33)式中：δmat(s)——多孔介质内的水蒸气渗透率，kg/(m·s·Pa)；h——沿着重力方向的距离，m。因此，多孔内部势流利用对流扩散方程形式的PDE方程表示为：da∂H∂t+∇-c∇H+β·∇H=f (34)方程系数分别为：da=0，c=δmat，β=0，f=0。H为因变量，单位为Pa。1.3.2空气中的速度场空气中的连续性方程利用Navier-Stokes方程表示：ρair∂u∂t+ρair(u·∇)(u)=-∇P+ρairg+μairΔu (35)动量变化率由微元体的所受的体积力、静压力、黏性阻力组成，采用经典偏微分方程中的对流扩散方程形式的PDE方程表示为：da''∂H'∂t+∇·-c''∇H'+β''·∇H'=f'' (36)方程系数分别为：da''=0，c''=δ0，β''=0，f''=0。δ0为空气的扩散系数，取1.96×10-10s。1.4边界条件外表面的外部平均压力pav作为狄利克雷边界条件。通过表面的湿通量g0为：g0=βp×pv_i-pv (37)pv=φ·psatT (38)式中：pv_i——外部环境水蒸气分压力；β——水分迁移系数，s/m；psatT——饱和水蒸气分压力，Pa。通过表面的热通量q0为：q0=h×Text-T+g0×Lv(Text) (39)式中：h——传热系数，W/(m2·K)；Lv(Text)——汽化潜热，kJ/kg。2模型验证2.1物理模型多层墙体物理模型采用附加空气层的多层墙体模型进行验证，COMSOL设置与自编的PDE模型采用相同的材料物性参数。带空气层的多层墙体物理模型如图1所示。各材料的物性性质[19-23]如表1、表2所示。不同材料在不同相对湿度下的含水量、导热系数和扩散系数如表3~表5所示。10.3969/j.issn.1004-7948.2023.07.002.F001图1带空气层的多层墙体物理模型10.3969/j.issn.1004-7948.2023.07.002.T001表1各材料的物理性质材料密度/(kg/m3)热容/[J/(kg·K)]孔隙率渗透率/m2蒸汽阻力因子石膏板5741 1000.6805.40×10-185.2泡沫混凝土6009200.6061.00×10-1792.0玻璃棉板738500.9502.69×10-101.3外饰面1 3608500.2005.80×10-188.110.3969/j.issn.1004-7948.2023.07.002.T002表2不同温度时空气层的物理性质项目公式密度352.716T恒压热容1 010.9+0.043 9T-2.92×10-4T2+6.5×10-7T3扩散系数/(m2/s)2.6×10-5动力黏度/Pa·s1.8×10-5导热系数-8.4×10-4+1.1×10-4T-8.6×10-8T2+6.5×10-7T310.3969/j.issn.1004-7948.2023.07.002.T003表3不同材料在不同相对湿度下的含水量材料相对湿度00.50.70.91.0石膏板05.76.517.0370.0泡沫混凝土01020180290玻璃棉板00.530.582.504.10外饰面01.5811.4013.50220.00kg/m310.3969/j.issn.1004-7948.2023.07.002.T004表4不同材料在不同相对湿度下的导热系数材料相对湿度00.50.70.91.0石膏板0.190.200.210.250.60泡沫混凝土0.14×(1+3×wc(φ)/600玻璃棉板(0.031 8-0.028 3)×T-263.1530+0.028 3外饰面0.9000.9080.9520.9712.060［W/（m·K）］10.3969/j.issn.1004-7948.2023.07.002.T005表5不同材料在不同相对湿度下的扩散系数材料相对湿度00.50.70.91.0石膏板185185185741159 000泡沫混凝土1×10-10+0.99×10-8·φ玻璃棉板0.01外饰面00.3310.2940.46610.000×10-12 m2/s墙体的初始温度为16 ℃，初始相对湿度0.8；室外温度为12 ℃，相对湿度0.7，对流传热系数25 W/(m2·K)，水分迁移系数βout为1.4×10-7 s/m；室内温度为20 ℃，相对湿度为0.55，对流传热系数为7.7 W/(m2·K)，水分迁移系数βout为7.7×10-9 s/m，室内相比室外维持5 Pa的正压。2.2结果分析为了简化模拟过程，将模型转化为二维图形，带空气层的多层墙体简化二维模型如图2所示。中轴线L上的数据代表各层材料的数据的平均值。采用相同的求解自由度数为161 247时，求解器设置求解时间为4 000 h，步长为5，容差设置为5×10-5，采用自有接口的求解时间为348 s。采用相同的设置时，新建模型求解时间仅为163 s，节省了一半的计算时间。墙体内温度、相对湿度随材料层的分布如图3、图4所示。各层达到稳态时水蒸气的平均速度和相对误差如表6所示。10.3969/j.issn.1004-7948.2023.07.002.F002图2带空气层的多层墙体简化二维模型10.3969/j.issn.1004-7948.2023.07.002.F003图3墙体内温度随材料层的分布10.3969/j.issn.1004-7948.2023.07.002.F004图4墙体内相对湿度随材料层的分布10.3969/j.issn.1004-7948.2023.07.002.T006表6各层达到稳态时水蒸气的平均速度和相对误差测点平均速度/(m/s)相对误差/%新建模型COMSOL模型a8.956 0×10-128.937 5×10-120.201b4.808 9×10-94.809 6×10-90.015c1.658 0×10-111.655 0×10-110.199d4.461 2×10-44.452 2×10-40.202e9.618 9×10-129.599 6×10-120.201由图3、图4及表6可知，新建模型的温度和速度的仿真结果与COMSOL软件模拟结果吻合良好。而相对湿度随着材料变化呈现差异。原因为考虑内部的液态水时仅考虑液态水扩散，未考虑在压差作用下的流动传递。对于水蒸气的渗透，δ0为静止空气中的水蒸气渗透率，在受迫流动时实际值应大于δ0。而COMSOL软件在处理空气域的水蒸气扩散时获得的结果明显高于新建模型。COMSOL模型获得的相对湿度斜率大于新建模型。由于COMSOL模型采用含湿多孔介质模型时采用相对渗透率，导致多孔介质区域的湿传递小于构造的PDE模型。数值误差可以通过提高网格质量与减小计算容差等手段不断减小。COMSOL软件在研究含有空气域的多孔介质热湿耦合中，无法采用“建筑材料”模块，仅能采用“吸湿性多孔介质”模块计算。由于使用孔隙率条件，无法将模型简化为一维的模型。而自行构造的PDE形式不受此限制。由于软件自有接口的二维模型考虑水蒸气的横向扩散传递与自然对流作用，导致水蒸气的传输大于空气域中不考虑热效应引起的水分扩散。3结语构建含空气层的多层墙体共轭的模型，并通过测试案例验证模型的可行性。与COMSOL软件自有接口相比，新建模型相同条件下的计算速度提升50%。新建的PDE方程具有更好的扩展性，可以适应从一维到三维的不同场景，补充了对含有空气层的建筑多物理场缺失的不足。此外，COMSOL软件自有模型接口无法将空气域与建筑热物理模型耦合，需要额外构建对流扩散方程的PDE形式以施加对流所需的速度场，新建的PDE模型不受空气域热扩散等物理现象对维度的制约，文中甚至可以采用1维的线段使问题进一步简化。