CN109214113A - 一种建筑围护结构动态传热的新模型降阶求解方法 - Google Patents

Abstract

本发明公开了一种建筑围护结构动态传热的新模型降阶求解方法，针对建筑动态传热问题，基于热平衡原理，建立围护结构动态传热的控制方程；针对室内外气象参数和围护结构热工性能参数，分析偏微分方程的特征，建立边界条件和初始条件；利用稳健离散方法，获得偏微分方程离散后的状态空间形式的原始系统；针对原始系统，基于Laguerre正交多项式，构建一种新型模型降阶求解方法；最后，将降阶系统的解还原至原始系统的解。相对于Krylov子空间模型降阶方法在频率域上进行构造降阶矩阵，新方法在时间域上处理更加直观简单；而相对于平衡截断模型降阶方法，计算资源消耗减少了。避免了本征正交分解(POD)模型降阶方法要求事先得到系统状态的数据集。

Description

一种建筑围护结构动态传热的新模型降阶求解方法

技术领域

本发明属于建筑性能模拟技术领域，涉及建筑动态热性能数值模拟的一种快速精准求解方法，具体涉及一种建筑围护结构动态传热的新模型降阶求解方法。

背景技术

在建筑热环境调控方案的制定及其系统的设计、建筑全年能耗准确评估以及动态建筑信息模型(BIM)的构建中，迫切需要准确快速计算建筑在一年中任意时刻围护结构动态传热量及其温度分布。尽管目前许多软件可进行建筑动态能耗计算，如BLAST、DOE-2、EnergyPlus、ES、TRNSYS、DeST等，然而由于计算方法的不同，针对同一建筑使用不同建筑能耗模拟软件的计算结果会不同。基于热平衡的偏微分方程数值求解是其中精准求解方法，不仅能求解简单的一维线性问题，而且能求解任何复杂围护结构的多维瞬态非线性传热问题。然而，针对实际室外气象参数非周期性变化，围护结构非线性传热规律，进行全年逐时或逐分的精准计算，直接数值求解方法的计算量巨大。

模型降阶(MOR)方法是将一个较大复杂系统转化为一个近似的较小系统,从而降低大型复杂系统的理论分析难度，减少数据存储量和运算量,加速系统的模拟计算，因此在大规模集成电路、自动化控制和机械等工程领域得到许多应用。目前建筑性能模拟领域开始涉及的模型降阶方法主要有：krylov子空间、平衡截断、本征正交分解(POD)三类模型降阶方法。krylov子空间模型降阶方法是一种比较成熟的算法，它可以匹配原始系统传递函数一定数量的矩，实现起来相对简单，计算资源消耗较小。平衡截断模型降阶方法能保持原始系统的稳定性，可以直接得到原始系统与降阶系统之间的误差关系，这是Krylov子空间所无法比拟的，但它最致命的缺点是计算耗费较大。POD模型降阶方法是先通过数值模拟或是数据采样等方式获得系统状态的数据集,然后根据数据集构造一个标准正交投影矩阵,最后利用正交矩阵通过投影变换得到降阶系统，其降阶精度完全依赖于数据集本身。基于POD的模型降阶方法由于过程简单且其与系统本身结构无关,使其在线性系统和非线性系统的模型降阶中都有应用。

尽管常规模型降阶方法在建筑性能模拟领域得到初步探讨，但目前很不成熟，因此，针对建筑动态传热问题，寻求性能优越而计算量较少的新型模型降阶方法，必将有助于建筑技术的进步。

发明内容

本发明的目的是提供一种建筑围护结构动态传热模拟的新型模型降阶求解方法，以时间域上直观处理方式实现建筑动态热性能的准确模拟及计算资源消耗的显著减少。该新型模型降阶方法是基于Laguerre正交多项式，在时间域上对系统的状态变量在以正交多项式为基底的空间中展开,然后利用正交多项式相应系数向量或矩阵构造所需降阶矩阵,并且通过该方法获得的降阶系统可以保持原始系统输出变量展开式中的系数,从而保证该方法的计算结果可靠而且计算量大为减少。该正交多项式模型降阶方法，不像Krylov子空间模型降阶方法需在频率域上进行构造降阶矩阵及只能求解线性系统，比平衡截断模型降阶方法的计算量少，也避免了POD模型降阶方法要求事先知道系统状态的数据集。

本发明所采用的技术方案是，

针对建筑动态传热问题，基于热平衡原理，建立围护结构动态传热的控制方程(关于时间和空间的偏微分方程)；针对室内外气象参数和围护结构热工性能参数，分析偏微分方程的特征，建立边界条件和初始条件；利用稳健离散方法，获得偏微分方程离散后的状态空间形式的原始系统；针对原始系统，基于Laguerre正交多项式，构建一种新型模型降阶求解方法；最后，将降阶系统的解还原至原始系统的解。

本发明主要包括以下步骤：

步骤一：分析围护结构构造，建立其非稳态传热控制方程。一般围护结构主体，因厚度相对于长宽尺寸很小，故可看作一维问题，其传热控制方程如下：

ρ—围护结构密度，kg·m^-3；

c—围护结构比热容，kJ·(kg·℃)^-1；

T—温度，℃；

λ—围护结构导热系数，W·(m·℃)^-1；

S—围护结构内热源,W·m^-2；

x—空间坐标，m；

t—时间,s；

步骤二：边界条件和初始条件的建立与处理。

为了数值求解控制方程(1)，还需补充其边界条件和初始条件。一般围护结构内、外侧边界条件如下:

T_i,a、T_o,a—围护结构内、外侧空气温度/℃，考虑到太阳辐射，室外空气温度取室外综合温度；n—壁面外法线方向向量；h—壁面处的换热系数/W·(m²·℃)^-1；下角标i,o—分别代表围护结构的内外侧；

一般围护结构内部的初始温度分布即初始条件是未知的，这里依据室外参数类似日周期性变化的特点，以第一天气象参数迭代计算来获得初始值T₀(x)，即初始条件。

步骤三：采用控制体积法，对步骤一建立的围护结构动态传热偏微分方程(1)进行空间离散，结合步骤二建立的边界条件(2)和(3)，得到状态空间形式描述的围护结构动态传热问题数学模型的原始系统，如下(4)式。

式中，为系统的状态变量，即为围护结构各空间离散节点的温度,其中n为空间离散节点数；为系统输入变量，即室内外温度(T_i,a、T_o,a)；为系统输出变量，即所需求解的问题，可以是围护结构任意位置的温度、热流密度或它们的统计值等等，m则随着输出变量的不同而不同。系数矩阵其中A主要与围护结构热工性能参数和空间离散尺度有关，B与围护结构内外侧传热特性参数有关，C和D则随着输出变量的不同而不同，它们的具体计算详见后面的“具体实施方式”。

步骤四：对步骤三得到的n阶原始系统(4)式,采用Laguerre多项式降阶方法进行降阶，得到r阶降阶系统如下，使得r《n。

式中，为降阶系统相应的系数矩阵；—降阶系统的状态变量，—降阶系统的输出变量；

步骤五：求解降阶系统(5a)式,得到降阶系统状态变量的解。

步骤六：通过步骤四求得的降阶系统系数矩阵和步骤五求得的降阶系统状态变量由式(5b)计算得到输出变量即为原始系统(4)式输出变量y(t)的近似解，即为控制方程式(1)在边界条件(2)、(3)下输出变量的数值解。

与现有技术相比，本发明的正交多项式模型降阶方法是在时间域下对状态变量以Laguerre正交多项式为基底下展开，降阶过程实际上是在求解状态变量T(t)在正交多项式下的展开系数，找出展开系数所在空间，并求出相应空间下的标准正交基，该方法较直观自然，且因降阶系统可以保持原始系统输出变量展开式中的系数从而保证计算精度高。如果直接计算原始系统(n阶)，则计算规模十分庞大，而通过此降阶方法，我们只需求解一个规模较小的系统(r阶，r《n)，从而节约了计算资源，大大缩短计算程序的运行时间。

与建筑性能现有讨论的模型降阶技术相比，本发明具有以下有益效果。相对于Krylov子空间模型降阶方法在频率域上进行构造降阶矩阵，新方法在时间域上处理更加直观简单；而相对于平衡截断模型降阶方法，计算资源消耗减少了。避免了本征正交分解(POD)模型降阶方法要求事先得到系统状态的数据集。

附图说明

图1是建筑围护结构正交多项式模型降阶方法的技术路线图；

图2是围护结构节点划分示意图；

图3是室外综合温度的逐时变化图；

图4是墙体外表面计算温度的逐时变化图；

图5是墙体内表面计算温度的逐时变化图；

图6是多项式模型降阶方法关于墙体内、外表面计算温度的相对误差；

图7是多项式模型降阶方法与直接算法的计算总耗时间对比图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

一种建筑围护结构动态传热的新模型降阶求解方法，是基于Laguerre正交多项式，直接在时间域上对系统的状态变量在以正交多项式为基底的空间中展开,然后利用正交多项式相应系数向量或矩阵构造所需降阶矩阵,其技术路线如图1所示，详细步骤如下。

步骤一：分析围护结构构造，建立围护结构动态传热物理数学模型。

对于一般围护结构主体，因厚度相对于长宽尺寸很小，可看做一维问题；室外气候参数具有24小时周期性昼夜变化的特点，同时一年四季逐日也在变化，故围护结构主体的传热近似为一维瞬态过程，其动态控制方程为：

ρ—围护结构密度，kg·m^-3

c—围护结构比热容，kJ·(kg·℃)^-1

T—围护结构温度，℃

λ—围护结构导热系数，W·(m·℃)^-1

S—围护结构内热源，W·m^-2

x—空间坐标，m

t—时间，s

步骤二：边界条件和初始条件的建立与处理。

围护结构内、外侧边界条件如下:

T_i,a—围护结构室内空气温度/℃；T_o,a—围护结构室外综合温度/℃；n—壁面外法线方向向量；h—壁面处的换热系数/W·(m²·℃)^-1；下角标i,o—分别代表围护结构的内外侧表面。

一般围护结构内部的初始温度分布是未知的，本发明依据室外参数日周期性变化的特点，以第一天气象参数迭代计算来获得初始值T₀(x)，即初始条件。

步骤三：围护结构动态传热控制方程的空间离散，得到求解问题的原始系统。

如上所示，围护结构动态传热过程的控制方程(1)为关于时间和空间的偏微分方程，无法直接求解，首先需要进行空间离散，转化为关于时间的常微分方程组。围护结构的节点划分如图2所示,共分为n个子域，即n个内节点。其中(Δx)_i表示第i个子域的大小，(δx)_i表示第i个节点到第(i+1)个节点的距离。对于均分网格，则有(Δx)_i＝L/n，i＝1,…,n,其中L为围护结构总厚度；(δx)_i＝L/n，i＝1,…,n-1,但(δx)₀＝(δx)_n＝L/2n。这里采用控制体积法，对围护结构动态传热偏微分方程(1)进行空间离散，结合边界条件(2)和(3)，得到状态空间形式描述的围护结构动态传热问题数学模型的原始系统，如下式(4)。

式中,为系统的状态变量，即为围护结构各空间离散节点的温度,其中n为空间离散节点数；为系统输入变量，可以是室内外温度(T_i,a、T_o,a)或围护结构中内热源S引起的热量扰动，其中p为输入变量数，对于一般外围护结构，其内无内热源，故这里只讨论双输入变量(T_i,a、T_o,a)的原始系统(p＝2)；为系统输出变量，可以是围护结构任意位置的温度、热流密度或它们的统计值等等，m是输出变量数，这里以输出围护结构内、外侧壁面温度(即m＝2)为例。

为系数矩阵。当围护结构内热源s为0、以围护结构内外侧壁面温度(T_i、T_o)为输出变量时，变量T(t)、u(t)、y(t)的表达式和系数矩阵A、B、C、D计算式如下：

其中，a₁'_,W,a'_n,E,a_i,E,a_i,W,a_i,P(i＝1,n)分别为：

步骤四：对原始系统采用Laguerre多项式降阶方法进行降阶，得到降阶系统。如下式(5)，使得降阶系统的阶数r远小于原始系统的阶数n。一般r减小，计算误差增大而计算量减少，因此对于每一具体问题，先需实验找到合理的降阶系统的阶数。

式(5)中各系数矩阵的求解过程如下：

首先对矩阵B,D进行分块处理，B＝[B₁B₂],B₁,B₂∈Rⁿ，D＝[D₁D₂],D₁,D₂∈R²，对原始系统的输入函数分别取

其中u₁(t)＝T_o,a∈R,u₂(t)＝T_i,a∈R将之分别代入到原始系统(4)中，得到：

和

假设原始系统(4)的解为T(t)∈Rⁿ，系统(6)和系统(7)的解分别为T₁(t)∈Rⁿ，T₂(t)∈Rⁿ，易知T(t)＝T₁(t)+T₂(t)，y(t)＝y₁(t)+y₂(t),即系统(6)和系统(7)的解的和即为原始系统的解。系统(6)和系统(7)为线性单输入系统，由此求解原始系统的多输入问题可以转化为通过求解系统(6)和系统(7)的问题即求解单输入系统问题。

现仅给出系统(6)利用Laguerre多项式的处理过程，系统(7)的处理方式同系统(6).具体实施如下，将系统(6)式在时间区域[0,t]上积分，得到：

假设T₁(t)是系统的状态方程(8)的解，这里将目标函数T₁(t)利用Laguerre多项式作为基底对它进行近似计算，即：

其中：a_1k∈Rⁿ,k＝0,1,2,…,M。M的选取决定了降阶系统的阶数r，对于双输入系统，它们之间关系为注意，对Laguerre多项式，有下列积分公式成立：

同样对输入变量u(t)在区间[0,t]的积分值，利用Laguerre多项式作为基底来对它进行近似计算，如下：

其中：u_1k∈R，k＝0,1,2,…,M。将状态变量展开式(9)和输入变量积分展开式(11)代入式(8a)的第一个公式,得到：

因为L₀(t),L₁(t),…,L_N(t)线性无关，所以通过利用上式两端对应系数相等，得到如下方程组：

通过(13)式求解出a_1k，得到如下等式：

其中：

易知当u₁₀不为零时，矩阵U是可逆矩阵，由此可知

并且有：

同样，对于系统(7)做系统(6)相同处理可得到：

其中：

因此只需对的基底向量进行处理，即分别针对上述的两个Krylov子空间中的每一个向量实施Arnoldi过程，由此可得到基于Laguerre降阶算法：过程归纳如下。

(1)：输入系数矩阵A,B,C,D和降阶阶数r；

(2)：对矩阵B，D进行分块处理：B＝[B₁B₂],B₁,B₂∈Rⁿ，D＝[D₁D₂],D₁,D₂∈R²；

(3)：分别计算得到空间：其中

(4)：对上述得到的两个Krylov子空间分别进行Arnoldi处理，即：

(a)初始化向量

(b)计算

(c)依次下去，计算第M+1项标准正交向量

(d)得到系统(6)，(7)的降阶矩阵：V₁＝[v₁₀,v₁₁,…,v_1M]，V₂＝[v₂₀,v₂₁,…,v_2M]

(5)：对矩阵V₁，V₂进行标准正交化处理得到原始系统(4)的降阶矩阵V。

(6)分别通过关系式求解降阶系统的系数矩阵，输出降阶系统的系数矩阵：

步骤五：求解降阶系统式(5a),得到降阶系统的解

步骤六：通过步骤四求得的降阶系统系数矩阵和步骤五求得的降阶系统状态变量由式(5b)计算得到输出变量即为原始系统(4)式输出变量y(t)的近似解；可通过步骤四求得的降阶矩阵V以及关系式求得原始系统式(4)的状态变量T(t)，即为控制方程式(1)在边界条件(2)、(3)下的数值解。

下面结合具体案例，依据上述技术路线的详细步骤，编写并运行计算程序，通过计算结果分析本发明的建筑围护结构正交多项式模型降阶方法的可靠性及有益效果。

不失一般性，以某多层外墙为例。该墙由四层结构组成，从外到内分别是：0.01m的建筑保温浆料粉刷层，0.1m的岩棉板保温层，0.3m的现浇钢筋混凝土结构和0.01m的建筑保温浆料粉刷层，其热物性参数见表1。

表1某外墙各层材料热物性参数

该外墙内外壁面的综合换热系数分别取为h_o＝23W/(m²·K)和h_i＝8.7W/(m²·K)。假设空气室内为空调，温度维持在26℃。依据某地典型年中6、7月的室外温度和南侧太阳辐射强度，计算的逐时室外综合温度分布如图3所示。

这里为了便于计算结果的形象对比以及提供高精度的数据，在控制方程离散时，空间网格间距均取0.1mm，因墙体总厚度为4200mm,则内节点数n＝4200，即原始系统的阶数为4200。时间步长取3.6秒。一般工程计算中，为了提高速度，时间和空间步长都可取得较高些，依据研究问题或工程应用的计算精度要求而确定。在利用Laguerre多项式进行降阶处理时，M取为20，则降阶系统的阶数r＝42。

针对该案例的上述已知条件(墙体热工参数和墙体内外侧环境参数)，依据图1技术路线的详细步骤，编写并运行计算程序，得到分别采用Laguerre正交多项式模型降阶方法和原始系统直接求解方法模拟该墙在夏季6、7月的动态传热过程的计算结果。所谓原始系统直接求解方法，即没有上述步骤四(降阶矩阵及降阶系统系数矩阵的计算)和步骤五(降阶系统的求解)，直接求解步骤三得到的原始系统。计算程序采用Matlab2016a编写，用于运行程序的工作站主要配置为：Core^TMI7-4790，内存16G(DDR3 1600MHz)，显卡NVIDIA QuadroK2000(2G)。

基于正交多项式模型降阶方法与直接求解方法，该外墙的计算结果对比如下。

图4和图5分别是墙体外、内表面计算温度的逐时变化图，两种算法的相对误差如图6；图中有两种线型，但因为二者非常接近，图例中有标出；同时更说明此方法效果良好。另外，由于数据点非常多，每个线性加标记符号将影响显示效果，请用修订模式打开文档，可见添加标记符号后图像效果。。其中，相对误差定义如下：

式中y(i)、y′(i)分别是直接求解方法、正交多项式模型降阶方法关于输出量的计算值。

可以看出该Laguerre正交多项式模型降阶方法与直接求解方法的计算结果吻合程度高，两者之间误差很小，均小于4×10^-6。由此可见，该Laguerre正交多项式模型降阶方法的计算准确度高。

该正交多项式模型降阶方法与直接求解方法的程序运行耗时(即计算耗时)对比如图7所示，可见正交多项式模型降阶方法的计算耗时远小于直接求解方法的，并且随着计算的持续时间(天数)增加，两者差别更大。这是因为直接求解方法的计算耗时基本与计算的持续时间(天数)成正比；而关于模型降阶方法，当计算的持续时间(天数)从7天增加到56天时，其计算耗时仅增加0.8秒。进一步比较，计算8周(56天)时，正交多项式模型降阶方法的计算耗时仅为直接求解方法的2.8％。由此可见，该正交多项式模型降阶方法不仅计算准确度高，而且能有效缩短程序运行时间及减少内存量，计算效率明显提升。依此推算，当计算的持续时间为一年时，计算耗时的减少将是十分巨大。

Claims (5)

1.一种建筑围护结构动态传热的模型降阶求解方法，其特征在于，针对建筑动态传热问题，首先基于热平衡原理建立围护结构动态传热的控制方程，即关于时间和空间的偏微分方程；然后针对室内外气象参数和围护结构热工性能参数，分析偏微分方程的特征，建立边界条件和初始条件；利用稳健离散方法，获得偏微分方程离散后的状态空间形式的原始系统；针对原始系统，基于Laguerre正交多项式构建新型模型降阶求解方法；最后，将降阶系统的解还原至原始系统的解。

2.根据权利要求1所述的一种建筑围护结构动态传热的模型降阶求解方法，其特征在于，包括以下步骤：

步骤一：分析围护结构构造，建立其非稳态传热控制方程：

ρ—围护结构密度，kg·m^-3；

c—围护结构比热容，kJ·(kg·℃)^-1；

T—温度，℃；

λ—围护结构导热系数，W·(m·℃)^-1；

S—围护结构内热源,W·m^-2；

x—空间坐标，m；

t—时间,s；

步骤二：边界条件的建立和初始条件的处理；

补充控制方程(1)边界条件和初始条件，围护结构内、外侧边界条件如下：

T_i,a、T_o,a—围护结构内、外侧空气温度/℃，考虑到太阳辐射，室外空气温度取室外综合温度；n—壁面外法线方向向量；h—壁面处的换热系数/W·(m²·℃)^-1；下角标i—代表围护结构的内侧；下角标o—代表围护结构的外侧；依据室外参数类似日周期性变化的特点，以第一天气象参数迭代计算来获得初始值T₀(x)，即为初始条件；

步骤三：采用控制体积法，对步骤一中建立的围护结构动态传热偏微分方程(1)进行空间离散，结合步骤二建立的边界条件(2)和(3)，得到状态空间形式描述的围护结构动态传热问题数学模型的原始系统，即公式(4)：

式中，为系统的状态变量，即为围护结构各空间离散节点的温度,其中n为空间离散节点数；为系统输入变量，即室内外温度(T_i,a、T_o,a)；为系统输出变量，即所需求解的问题，可以是围护结构任意位置的温度、热流密度或它们的统计值，m则随着输出变量的不同而不同；系数矩阵其中A与围护结构热工性能参数和空间离散尺度有关，B与围护结构内外侧传热特性参数有关，C和D则随着输出变量的不同而不同；

步骤四：对步骤三得到的n阶原始系统公式(4),采用Laguerre多项式降阶方法进行降阶，得到r阶降阶系统如下，使得r＜＜n；

式中，为降阶系统相应的系数矩阵；—降阶系统的状态变量，—降阶系统的输出变量；

步骤五：求解降阶系统(5a)式,得到降阶系统状态变量的解；

步骤六：通过步骤四求得的降阶系统系数矩阵和步骤五求得的降阶系统状态变量由式(5b)计算得到输出变量即为原始系统(4)式输出变量y(t)的近似解，即为控制方程式(1)在边界条件(2)、(3)下输出变量的数值解。

3.根据权利要求2所述的一种建筑围护结构动态传热的模型降阶求解方法，其特征在于，所述步骤二中，围护结构内部的初始温度分布是未知的，这里依据室外参数类似日周期性变化的特点，以第一天气象参数迭代计算来获得初始值T₀(x)。

4.根据权利要求2所述的一种建筑围护结构动态传热的模型降阶求解方法，其特征在于，所述步骤四中，系数矩阵A,B,C,D的求解过程包括如下步骤：

(1)：输入系数矩阵A,B,C,D和降阶阶数r；

(2)：对矩阵B，D进行分块处理：B＝[B₁ B₂],B₁,B₂∈Rⁿ，D＝[D₁ D₂],D₁,D₂∈R²；

(3)：分别计算得到空间：其中

(4)：对上述得到的两个Krylov子空间分别进行Arnoldi处理得到降阶矩阵；

(5)：对步骤(4)得到的降阶矩阵进行标准正交化处理得到原始系统(4)的降阶矩阵V；

(6)：分别通过关系式求解降阶系统的系数矩阵，输出降阶系统的系数矩阵：

5.根据权利要求4所述的一种建筑围护结构动态传热的模型降阶求解方法，其特征在于，所述步骤(4)包括如下步骤：

(a)初始化向量

(b)计算

(c)依次下去，计算第M+1项标准正交向量

(d)得到系统(6)，(7)的降阶矩阵：V₁＝[v₁₀,v₁₁,…,v_1M]，V₂＝[v₂₀,v₂₁,…,v_2M]。