2018, Vol. 41
相变材料通过相态变化实现潜热存储与释放，解决能源供应时间与空间矛盾，提高能源的利用效率[1-3]。在建筑及暖通领域，相变材料多以平板的形式嵌入墙体，以实现建筑围护结构的蓄能隔热[4-5]。文献[6]提出了一种新型双层相变墙系统，通过模拟发现该相变墙系统能够降低建筑35.4%的峰值负荷与12.8%的全年能耗。文献[7]将相变材料置于房间墙体，对波兰、法国马赛等不同气候区下房间的夏季得热量模拟进行模拟，发现该相变墙体能够减少房间20%的得热量。文献[8]提出了一种结合太阳能集热器的相变蓄能地板采暖系统，实验结果表明，该相变蓄能地板能有效的维持房间温度，提高房间热舒适性。
近年来，不少学者对相变材料的传热模型进行了研究。文献[6]基于有效焓的方法建立了相变平板的传热模型，采用有限差分法对模型进行求解计算。文献[9]基于有效热容假设对相变墙房间进行了数值模拟，进一步分析了相变墙表面温度及热流响应。文献[10]分别利用有效热容法及有效焓法建立了相变材料的传热模型，通过数值模拟对相变温度区间的2种模型的适用性及准确性进行了验证。传统数值求解方法求解复杂，计算效率低，文献[11]探讨了将相变板与墙体结合在一起采用动态简化模型应用于轻型建筑结构、中型建筑结构、重型建筑结构的传热计算的准确性。文中指出该模型在相变板与轻型墙体、中型墙体结合在一起时可获得更好的准确度。但相变板与重型墙体结合在一起时，采用该简化模型计算结果准确性较差，该文献并未就相变平板本身简化模型的准确性及适用性进行分析。
文中分别建立了相变平板可变热容热阻的一维一阶、二阶与三阶动态简化模型，利用数值模型的计算结果辨识模型的相关参数，并对相变平板的传热特性进行模拟。利用数值结果对不同的简化模型进行对比分析，进一步验证简化模型的准确性与适用性，并指出二阶模型计算量适中、计算结果准确性也很高，可为相变板或相变板与普通墙体集成的传热计算提供参考。
1 动态简化传热模型相变平板的结构如图 1(a)所示，图 1(b)~图 1(d)分别为相变平板一阶(2R1C)、二阶(4R2C)与三阶(6R3C)一维简化RC动态传热模型原理图，该简化模型能够模拟相变材料在不同边界条件下的动态传热过程，计算相变材料表面的温度与通过相变材料的热流，进而模拟相变平板材料的传热特性。以4R2C简化传热模型为例，采用集总参数[12]的思想将相变材料的总热容分为2个集总等效热容(C1，C2)，总热阻分为2个集总等效热阻(R1、R2、R3，R4)，进而根据利用热力学第三定律建立相变平板的一维动态简化RC传热模型。
|
图 1 相变平板结构与动态简化模型 Figure 1 Structure and dynamic simplified models of planar PCM |
在对流边界下，4R2C简化传热模型如式(1)与式(2)所示，通过左、右表面的热流q1与q2，可根据式(3)进行计算，相变平板左、右表面的温度T1、T2可根据式(4)~式(5)进行计算，其中Tout1与Tout2分别为相变平板左、右侧环境温度，℃；Tp1与Tp2分别表示材料内节点的温度，℃；h1与h2分别为左右两边的对流换热系数，W/(m2·℃)；t为时间变量，s。同理，可依次建立对流边界下，相变平板的2R1C与6R3C简化传热模型，如图 1所示。在简化模型中，假设相变材料是各向同性的均匀连续介质，其R、C参数通过辨识获取。
| $ {C_1}\frac{{{\rm{d}}{T_{p1}}}}{{{\rm{d}}t}} = \frac{{{T_{{\rm{out1}}}} - {T_{{\rm{p1}}}}}}{{1/{h_1} + {R_1}}} - \frac{{{T_{{\rm{p1}}}} - {T_{{\rm{p2}}}}}}{{{R_1} + {R_3}}}。$ | (1) |
| $ {C_2}\frac{{{\rm{d}}{T_{{\rm{p2}}}}}}{{{\rm{d}}t}} = \frac{{{T_{{\rm{p1}}}} - {T_{{\rm{p2}}}}}}{{{R_1} + {R_3}}} - \frac{{{T_{{\rm{p2}}}} - {T_{{\rm{out2}}}}}}{{{R_4} + 1/{h_2}}}。$ | (2) |
| $ {q_1} = \frac{{{T_{{\rm{out1}}}} - {T_{{\rm{p1}}}}}}{{1/{h_1} + {R_1}}},{q_2} = \frac{{{T_{{\rm{p2}}}} - {T_{{\rm{out2}}}}}}{{{R_4} + 1/{h_2}}}。$ | (3) |
| $ {T_1} = \left( {{h_1}{R_1}{T_{{\rm{out1}}}} + {T_{{\rm{p1}}}}} \right)/\left( {{h_1}{R_1} + 1} \right)。$ | (4) |
| $ {T_2} = \left( {{h_2}{R_4}{T_{{\rm{out2}}}} + {T_{{\rm{p2}}}}} \right)/\left( {{h_2}{R_4} + 1} \right)。$ | (5) |
在非稳态边界条件下，获取相变材料热特性的理论解十分困难，多采用数值方法[6, 9-10]求解，其数值结果能够较为准确地反映相变材料的热特性。相变平板的总热阻与总热容Rto、Cto，可根据平板节点温度利用式(6)与式(7)进行计算[13]，简化模型R、C参数辨识可以简化为寻找1组最优的R、C分配比例(βi)，进而确定各节点之间的热阻Ri、热容Ci参数，使得相变平板的温度与热流响应与数值解的误差最小。以一阶与二阶模型为例，(β)2R1C为2R1C简化模型R、C参数的分配比例，其Ri、Ci参数可利用式(8)计算；在4R2C简化模型中，将相变平板分为节点p1与节点p2的2部分，(β2)4R2C表示其节点分配比，(β2)4R2C与(β3)4R2C分别表示每个节点部分前后的热阻参数配比，模型Ri、Ci参数如式(9)~式(11)所示。同理，可确定6R3C简化模型的Ri、Ci参数，进而依次对各个模型的R、C参数进行辨识。
| $ {R_{{\rm{to}}}}\left( {{T_{{\rm{p}}i}}} \right) = \delta /{K_{{\rm{pcm}}}}\left( {{T_{{\rm{p}}i}}} \right),{K_{{\rm{pcm}}}}\left( {{T_{{\rm{p}}i}}} \right) = \left\{ \begin{array}{l} {K_{\rm{s}}},{T_{{\rm{p}}i}} < {T_{\rm{s}}};\\ \left( {T - {T_{\rm{s}}}} \right)\frac{{\left( {{K_{\rm{l}}} - {K_{\rm{s}}}} \right)}}{{\left( {{T_{\rm{l}}} - {T_{\rm{s}}}} \right)}},{T_{\rm{s}}} \le {T_{{\rm{p}}i}} \le {T_{\rm{l}}};\\ {K_{\rm{l}}},{T_{{\rm{p}}i}} > {T_{\rm{s}}}。\end{array} \right. $ | (6) |
| $ {C_{{\rm{to}}}}\left( {{T_{{\rm{p}}i}}} \right) = \delta {C_{{\rm{pcm}}}}\left( {{T_{{\rm{p}}i}}} \right),{C_{{\rm{pcm}}}}\left( {{T_{{\rm{p}}i}}} \right) = \left\{ \begin{array}{l} {\rho _{\rm{s}}}{c_{{\rm{ps}}}},{T_{{\rm{p}}i}} < {T_{\rm{s}}};\\ \frac{{{\rho _{\rm{s}}}{c_{{\rm{ps}}}} + {\rho _{\rm{l}}}{c_{{\rm{pl}}}}}}{2} + \frac{{{\rho _{\rm{s}}} + {\rho _{\rm{l}}}}}{2}\left( {\frac{\lambda }{{{T_{\rm{l}}} - {T_{\rm{s}}}}}} \right),{T_{\rm{s}}} \le {T_{{\rm{p}}i}} \le {T_{\rm{l}}};\\ {\rho _{\rm{l}}}{c_{{\rm{pl}}}},{T_{{\rm{p}}i}} > {T_{\rm{s}}}。\end{array} \right. $ | (7) |
| $ \begin{array}{l} {C_1} = {C_{to}}\left( {{T_{{\rm{pl}}}}} \right),{R_1} = {\left( {{\beta _1}} \right)_{{\rm{2R1C}}}} * {R_{to}}\left( {{T_{{\rm{pl}}}}} \right),\\ {R_2} = \left( {1 - {{\left( {{\beta _1}} \right)}_{{\rm{2R1C}}}} * {R_{to}}\left( {{T_{{\rm{pl}}}}} \right),} \right. \end{array} $ | (8) |
| $ \begin{array}{l} {R_1} = {\left( {{\beta _1}} \right)_{{\rm{4R2C}}}} * {\left( {{\beta _2}} \right)_{{\rm{4R2C}}}} * {R_{{\rm{to}}}}\left( {{T_{{\rm{pl}}}}} \right),\\ {R_2} = {\left( {{\beta _1}} \right)_{{\rm{4R2C}}}} * \left( {1 - {{\left( {{\beta _2}} \right)}_{{\rm{4R2C}}}} * {R_{{\rm{to}}}}\left( {{T_{{\rm{pl}}}}} \right)} \right., \end{array} $ | (9) |
| $ \begin{array}{l} {R_3} = \left( {1 - {{\left( {{\beta _1}} \right)}_{{\rm{4R2C}}}} * \left( {_{\beta 3} * {R_{{\rm{to}}}}\left( {{T_{{\rm{p2}}}}} \right)} \right.} \right.,\\ {R_4} = \left( {1 - {{\left( {{\beta _1}} \right)}_{{\rm{4R2C}}}} * \left( {1 - {{\left( {{\beta _2}} \right)}_{{\rm{4R2C}}}} * {R_{{\rm{to}}}}\left( {{T_{{\rm{p2}}}}} \right)} \right.,} \right. \end{array} $ | (10) |
| $ \begin{array}{l} {C_1} = {\left( {{\beta _1}} \right)_{{\rm{4R2C}}}} * {C_{{\rm{to}}}}\left( {{T_{{\rm{pl}}}}} \right),\\ {C_2} = \left( {1 - {{\left( {{\beta _1}} \right)}_{{\rm{4R2C}}}} * {C_{{\rm{to}}}}\left( {{T_{{\rm{p2}}}}} \right),} \right. \end{array} $ | (11) |
式中：δ为相变平板的厚度，m；Ks、Kl分别为相变材料固相与液相时的导热系数，W/(m·℃)；cps、cpl分别为相变材料固相与液相时的比热，J/(kg·℃)；ρs、ρl分别为相变材料固相与液相时的密度，kg/m3；λ为相变潜热，J/kg；Tpi表示平板节点i(1~3)的温度，℃，Ts、Tl分别为相变材料的凝固温度与熔化温度，℃。
文中以时域内通过简化模型计算相变平板的表面温度、热流与数值模型求解结果的误差最小值作为优化目标，利用遗传算法(GA)[14-15]分别对不同RC简化模型参数进行辨识。图 2表示RC简化模型参数的辨识过程，其中，J为优化目标函数，如式(12)所示，Tsz2、qsz2表示相变平板右壁面温度以及通过右壁面热流的数值解，N-1为误差样本总数，W为权重，可根据材料热阻确定。βi优化变量(R、C分配比例)，R、C参数越多，优化变量的数量也越多，其中，2R1C模型优化变量为1个、4R2C模型优化变量为3个、6R3C模型优化变量为5个。
| $ J\left( {{\beta _i}} \right) = \sqrt {\frac{{\sum\limits_{k = 1}^N {{{\left( {{q_2} - {q_{sz2}}} \right)}^2}} }}{{N - 1}}} + W\sqrt {\frac{{\sum\limits_{k = 1}^N {{{\left( {{T_2} - {T_{sz2}}} \right)}^2}} }}{{N - 1}}} 。$ | (12) |
|
图 2 模型参数辨识原理示意图 Figure 2 Schematic of parameters identification |
文中以某相变平板为例，利用ANSYS建立该相变平板的数值模型[16]，模拟不同边界条件下，相变平板材料的动态传热过程。相变平板材料选用SP29[17]，厚度为0.01 m，热物性参数如表 1所示。对数值模型进行网格无关性验证，确定网格数量为2万个，模拟时间步长为1.0 s，相变模型的迭代收敛精度为10-6，数值模型计算结果作为RC简化模型的参数辨识与验证依据。
|
表 1 相变平板材料物性参数 Table 1 Physical properties of Planar PCM |
根据式(1)~式(5)依次建立上述相变平板的简化RC模型，设计2种工况对简化模型的参数辨识及模型准确性进行研究，其周期性边界条件如式(13)~式(14)所示，各工况材料初始温度均为25 ℃，左、右表面对流换热系数分别为60 W/(m2·℃)与10 W/(m2·℃)。利用数值模型模拟工况1边界条件下相变平板的温度与热流响应，建立式(12)所示目标函数J，利用遗传算法依次辨识2R1C、4R1C及6R3C简化模型的R、C参数，并对不同简化模型的计算结果进行对比分析。此外，进一步利用数值模型与简化模型模拟了工况2(工况1与工况2，周期与极值不同)边界条件下的相变平板的温度与热流，如式(3)~式(5)，进而验证简化模型在不同边界条件下的准确性及适用性。
工况1(参数辨识)：
| $ \begin{array}{*{20}{c}} {{T_{{\rm{out1}}}} = 29 + 8 \times \sin \left( {\frac{{\rm{ \mathsf{ π} }}}{{12 \times 3\;600}}t + \frac{{\rm{ \mathsf{ π} }}}{3}} \right),}\\ {{T_{{\rm{out2}}}} = 27。} \end{array} $ | (13) |
工况2(模型验证)：
| $ \begin{array}{*{20}{c}} {{T_{{\rm{out1}}}} = 27 + 8 \times \sin \left( {\frac{{\rm{ \mathsf{ π} }}}{{6 \times 3\;600}}t + \frac{{\rm{ \mathsf{ π} }}}{3}} \right),}\\ {{T_{{\rm{out2}}}} = 27。} \end{array} $ | (14) |
利用ANSYS数值模型模拟工况1边界条件下相变平板材料的温度与热流响应，图 3为ANSYS数值模型计算收敛图，数值模型在进行模拟计算时，由于边界的变化，计算残差呈现周期性变化(存在非相变与相变过程)，但都小于设定的迭代收敛精度(10-6)，表明数值模拟结果收敛(计算时间为1.800 0 e+05 s，约2.1个周期)。
|
图 3 数值模型计算收敛图 Figure 3 Convergence graph of numerical model |
根据工况1边界条件下数值模型的计算结果，对简化模型的R、C参数进行辨识，结果如表 2所示，总热容Cto、热阻Rto可由式(6)~式(7)计算得到。
|
|
表 2 动态简化模型RC参数辨识结果 Table 2 RC parameters of dynamic simplified models |
根据表 2中的R、C参数，分别利用2R1C、4R2C及6R2C动态简化模型模拟工况1边界条件下，相变平板的表面温度及热流，并与数值模型的计算结果进行对比，如表 3所示。表面热流平均相对误差与温度平均误差可按式(15)与式(16)进行计算，qsz、Tsz表示采用数值模型计算的表面热流与温度，qsm、Tsm表示采用简化模型计算的表面热流与温度。结果表明，采用不同的简化模型模拟时，通过左表面的热流平均相对误差基本一致；4R2C简化模型计算得到的表面温度与右表面热流误差明显小于2R1C简化模型的计算结果，准确性更好；6R2C简化模型与4R2C简化模型计算结果误差相差不大，但采用6R2C简化模型时，需要辨识5个参数(4R2C辨识3个)，模型更复杂，计算成本更大。因此，在保证准确性的同时，4R2C简化模型计算成本更小，适用性更好。
| $ \Delta q = \sum\limits_{k = 1}^N {\left| {{q_{{\rm{sz}}}} - {q_{{\rm{sm}}}}} \right|} /\sum\limits_{k = 1}^N {\left| {{q_{{\rm{sz}}}}} \right|} \times 100\% 。$ | (15) |
| $ \Delta T = \sum\limits_{k = 1}^N {\left| {{T_{{\rm{sz}}}} - {T_{{\rm{sm}}}}} \right|} /\sum\limits_{k = 1}^N {\left| {{T_{{\rm{sz}}}}} \right|} \times 100\% 。$ | (16) |
|
|
表 3 工况1边界条件下的简化模型与数值模型误差对比 Table 3 Comparisons between simplified model and numerical model in case 1 |
为进一步验证动态简化模型在不同边界条件下的适用性，利用4R2C简化模型对工况2边界条件下相变平板材料的传热过程进行模拟。图 4、图 5分别为工况2边界条件下，相变平板材料左、右表面的温度与热流(纵坐标)随时间(横坐标)的变化曲线。结果可知，4R2C简化模型计算结果与数值结果十分吻合，相变平板左、右表面的温度平均误差分别为0.26 ℃(0.7%)与0.12 ℃(0.4%)；热流平均相对误差分别为10.4%与5.1%，进一步表明在不同的边界条件下，4R2C简化模型仍具有很好的准确性与适用性。
|
图 4 工况2下相变平板表面温度 Figure 4 Surface temperature of Planar PCM in case 2 |
|
图 5 工况2下通过相变平板左表面的热流(q1)与右表面热流(q2) Figure 5 Heat flux of left surface (q1) and Heat flux of right surface (q2) in case 2 |
文中所建立的动态简化模型能够有效提高计算效率。以工况2边界为例，在同一台计算机(i7处理器，8核16线程并行)上，利用ANSYS数值模型模拟5个周期(单个周期为12 h)内相变平板的传热特性所需耗时约2 d，采用简化模型在进行模拟求解时仅需1 min左右，小于数值模型计算所需时间的0.1%，计算效率明显提高。
5 结论针对相变平板分别建立了一阶(2R1C)、二阶(4R2C)及三阶(6R3C)热容热阻动态简化传热模型，根据具体工况(工况1)下数值模型的求解结果对简化模型的RC参数进行辨识，进而模拟相变平板材料的动态传热过程，并根据模拟结果对不同RC简化模型的准确性进行了对比分析。此外，进一步研究了4R2C简化模型在不同边界条件(工况2)下的适用性。
1) 在进行参数辨识后，4R2C与6R3C简化模型的准确性优于2R1C简化模型，在工况1边界条件下，4R2C与6R3C简化模型模拟的表面温度误差及热流相对误差分别为0.2 ℃(0.7%)与10%左右。在保证准确性的前提下，4R2C简化模型所需辨识参数较少，计算成本更小，适用性更好。
2) 在工况2边界条件下，采用二阶(4R2C)简化模型计算得到的相变平板左、右表面的温度及通过左、右表面的热流与数值模型的结果基本一致，其左、右表面温度温差均小于0.3 ℃，通过左、右表面热流平均相对误差分别为10.4%与5.1%，表明在不同边界条件下，4R2C简化模型仍具有很好的准确性与适用性。
3) 采用简化模型计算时，计算时间远小于数值模型的计算时间，计算效率明显提高。
| [1] | Pielichowska K, Pielichowski K. Phase change materials for thermal energy storage[J]. Progress in Materials Science, 2014, 65(10): 67–123. |
| [2] | Soares N, Costa J J, Gaspar A R, et al. Review of passive PCM latent heat thermal energy storage systems towards buildings' energy efficiency[J]. Energy and Buildings, 2013, 59: 82–103. DOI:10.1016/j.enbuild.2012.12.042 |
| [3] | Rodriguez U E, RuizV L, Vega S, et al. Applications of phase change material in highly energy-efficient houses[J]. Energy and Buildings, 2012, 50(7): 49–62. |
| [4] | Kuznik F, David D, Johannes K, et al. A review on phase change materials integrated in building walls[J]. Renewable and Sustainable Energy Reviews, 2011, 15(1): 379–391. DOI:10.1016/j.rser.2010.08.019 |
| [5] | Zhou D, Zhao C Y, Tian Y. Review on thermal energy storage with phase change materials (PCMs) in building applications[J]. Applied Energy, 2012, 92(4): 593–605. |
| [6] | Diaconu B M, Cruceru M. Novel concept of composite phase change material wall system for year-round thermal energy savings[J]. Energy and Buildings, 2010, 42(10): 1759–1772. DOI:10.1016/j.enbuild.2010.05.012 |
| [7] | Kosny J, Kossecka E, Brizezinski A, et al. Dynamic thermal performance analysis of fiber insulations containing bio-based phase change materials (PCMs)[J]. Energy and Building, 2012, 52(3): 122–131. |
| [8] |
肖伟, 王馨, 张群力, 等.
结合太阳能空气集热器的定形相变蓄能地板采暖系统实验研究[J]. 太阳能学报, 2008, 29(11): 1319–1323.
XIAO Wei, WANG Xin, ZHANG Qunli, et al. Experimental study on underfloor air supply system with air solar collector and shape-stabilized PCM[J]. Acta Energiae solaris sinica, 2008, 29(11): 1319–1323. DOI:10.3321/j.issn:0254-0096.2008.11.001 (in Chinese) |
| [9] | Kant K, Shukla A, Sharma A. Heat transfer studies of building brick containing phase change materials[J]. Solar Energy, 2017, 155(C): 1233–1242. |
| [10] | Jin X, Hu H, Shi X, et al. Comparison of two numerical heat transfer models for phase change material board[J]. Applied Thermal Engineering, 2018, 128: 1331–1339. DOI:10.1016/j.applthermaleng.2017.09.015 |
| [11] | Zhu N, Wang S W, Xu X H, et al. A simplified dynamic model of building structures integrated with shaped-stabilized phase change materials[J]. International Journal of Thermal Sciences, 2010, 49(9): 1722–1731. DOI:10.1016/j.ijthermalsci.2010.03.020 |
| [12] |
章熙民, 任泽霈. 传热学[M]. 第四版. 北京: 中国建筑工业出版社, 2001.
ZHANG Ximin, REN Zepei. Heat transfer[M]. fourth edition. Beijing: China Architecture & Building Press, 2001. (in Chinese) |
| [13] | Alawadhi E M. Thermal analysis of a building brick containing phase change material[J]. Energy and Buildings, 2008, 40(3): 351–357. DOI:10.1016/j.enbuild.2007.03.001 |
| [14] |
周明, 孙树栋. 遗传算法原理及应用[M]. 北京: 国防工业出版社, 1999.
ZHOU Ming, SUN Shudong. Principle and application of genetic algorithm[M]. Beijing: National Defence Industry Press, 1999. (in Chinese) |
| [15] | Wang S W, Xu X H. Simplified building model for transient thermal performance estimation using GA-based parameter identification[J]. International Journal of Thermal Sciences, 2006, 45(4): 419–432. DOI:10.1016/j.ijthermalsci.2005.06.009 |
| [16] |
柴国荣.
基于ANSYS的相变墙体传热特性计算分析[J]. 新型建筑材料, 2011, 38(7): 79–81.
CHAI Guorong. Numerical analysis of the heat transfer characteristics of phase change wall through ansys[J]. New Building Materials, 2011, 38(7): 79–81. (in Chinese) |
| [17] | Meng E, Yu Hang, Zhan G Y, et al. Experimental and numerical study of the thermal performance of a new type of phase change material room[J]. Energy Conversion and Management, 2013, 74: 386–394. DOI:10.1016/j.enconman.2013.06.004 |