计算传热学程序设计

1中国石油大学(华东)储运与建筑工程学院热能与动力工程系《计算传热学程序设计》设计报告学生姓名:学号：专业班级：指导教师2012年7月7日21、设计题目有一房屋的砖墙厚δ=0.3m，λ=0.85W/(m·℃)，ρc=1.05×106J/(m3·K)，室内温度Tf1保持20℃不变，表面传热系数h1=6W/(m2·℃)。开始时墙的温度处于稳定状态，内墙表面温度Tw1为15℃寒潮入侵后，室外温度Tf2下降为-10℃，外墙的表面传热系数为35W/(m2·℃)。试分析寒潮入侵后多少时间内墙壁面方可感受到外界气温的变化。图1墙壁简化图1.1已知参数壁厚，墙壁导热系数，密度与比热容的乘积，室内和寒潮入侵后室外空气温度，室内空气和外墙的表面传热系数，开始时稳定状态下的内墙表面温度。1.2求解寒潮入侵多少时间后内墙壁面可感受到外界气温的变化？2物理与数学模型2.1物理模型该墙面为常物性，可以假设：（1）其为无限大平面，（2）只有在厚度方向传热，没有纵向传热，则该问题转化为一维常物性无限大平面非稳态导热问题。2.2数学模型以墙外表面为坐标原点，沿厚度方向为坐标正方向，建立坐标系。基于上述模型，取其在x方向上的微元作为研究对象，则该问题的数学模型可描述如下：T()Tcxx（1a）初始条件：𝑤1−ℎ1(𝑓1−𝑤1)(δ−x)/λ（1b）在两侧相应的边界条件是第三类边界条件，分别由傅立叶定律可描述如下：室外寒流入侵室内0x3左边界：0202()xfxThTTX（1c）右边界：11()xfxThTTX（1d）3数值处理与程序设计3.1数值处理采用外点法用均匀网格对求解区域进行离散化，得到的网格系统如图2所示。一共使用了0~N-1共N个节点。节点间距δx为:δx1图2墙壁内的网格划分此例中墙壁导热系数为常值，无源项。则可采用有限体积法对控制方程离散化，得到离散方程为：ppEEWWaTaTaTb（2a）式中：0PWEPaaaa（2b）xaE，xaW，xcaP0（2c）00ppbaT（2d）其中的上标“0”表示此为上一时刻的值,分别为节点所在控制容积左右边界上的导热系数，由于墙壁导热系数不变，故都等于λ，△τ为时间步长。由元体能量平衡法可以得知左右边界节点的离散方程分别为：4左边界节点：(𝑎𝑝02+ℎ2+λ𝑥)0λ𝑥1+ℎ2𝑓2+𝑎𝑝0200（3）右边界节点：(𝑎𝑝02+ℎ1+λ𝑥)1λ𝑥2+ℎ1𝑓1+𝑎𝑝0210（4）离散方程的详细推导过程见附录。3.2程序设计由物理模型可以知道本问题为一维导热问题，一维导热问题的离散方程在取遍所有节点后形成的是三对角的代数方程组，采用追赶法进行求解。程序构成和方法:程序由主程序和一个子程序构成。主程序进行变量定义和各已知参数的输入，以及左右边界节点和内部节点控制方程的输入；子程序tdma实现追赶法用来计算每个节点新的温度。Thomas算法求解过程分为两步:消元和回代。消元是从系数矩阵的第二行起，逐一将每一行的非零元素消去一个，使原来的三元方程化为二元方程。消元进行到最后一行时，二元方程就化为一元方程，直接得到最后一个未知数的值。然后逐一往前回代，由各二元方程求出其它未知解。程序特点：该程序有很强的适应性，一维常物性非稳态平壁导热问题都可以使用此程序，只要适当更改边值条件即可。还可以进行修改解决非常物性问题。程序中对输出节点，最大输出量都进行了控制，对计算结果的分析有很大帮助。而且Thoms算法的优点需要内存小，工作量小，程序设计简单。程序流程图：首先对变量赋值，然后由初始条件建立初始温度场，接着从左边界，内部节点，到右边界进行迭代，直到满足精度要求为止，最后输出结果，程序结束。程序流程如下图3。4、模型与程序验证4.1模型本题简化为厚度为2=0.3m的一维非稳态模型如图4所示，初始温度为15℃，在其中间建立坐标系，左两边为对流换热，且换热系数相同都为h=25W/(m2·℃)，且流体温度Tf=-10℃对于x0，列出其导热微分方程式及定解条件：22(0,0)TTaxx（5）5ac0(,0)(0)TxTx（6）0(,)0xTxx（7）(,)(,)xTxhTTx（8）引入过余温度：(,)TxT（9）开始进行参数赋值用迭代法求解温度场计数并判断是否完成温度场计算输入至屏幕输入至文件程序结束NY图3程序流程图6直接根据公式得到解析解如下：210.(,)exp()cos()nnnnCFo（10）式中，2,axFo，系数nC应该使上述无穷级数在0是满足初始条件，由傅里叶级数理论可得：2sincossinnnnnnC（11）n是超越方程的根，称为特征根。tan,1,2,nnBin…（12）其中hBi。4.2程序验证（1）由模型可以得到相关信息然后进行编程，同等时间下计算出中心处温度的解析解和数值解进行比较，数据记录在表1。然后计算出相对误差，作图5，观察数值解与分析解的比较曲线。由图表中可以发现，平壁中心不同时刻温度值的分析解和数值解相差不是很大，二者吻合的比较好，可以说明所编制的数值解法的程序是正确的。相对误差先增大后减小，增大的原因是此时温度接近零度，相对误差的基数比较小，所以造成相对误差较大，但是此时的绝对误差并不大，在合理范围内，所以除去个别点外，都满足误差小于百分之1。可以验证所编数值解法的程序是正确的。右侧：h=25W/(m2·℃)Tf=-10℃左侧：hTf0δx-δT0=15℃图4一维导热简化模型7（2）空间步长对墙内壁的温度影响如图6及表2。在程序编写过程中用网格节点数对空间步长进行控制，为了观察空间步长对墙内壁温度的影响，表中选择了三个不同的空间步长，分别为选取51，101，201个网格节点，则相应的空间步长为0.006，0.003，0.0015。根据不同步长时温度的变化曲线可以看出，空间步长对内墙壁的影响不大，当空间步长控制在合理范围时可以忽略空间步长的影响。表1分析解与数值解比较时间（h）分析解（℃）数值解（℃）相对误差（%）0151500.55614.85414.811-0.289481.11113.46313.425-0.282261.66711.30611.3-0.053072.2229.0669.0810.1654532.7786.9746.9990.3584743.3335.0825.1120.5903193.8893.3933.4250.9431184.4441.891.9231.7460325．000.5540.5886.1371845.556-0.631-0.598-5.229796.111-1.684-1.651-1.959626.667-2.618-2.586-1.222317.222-3.447-3.417-0.870327.778-4.184-4.154-0.717028.333-4.837-4.809-0.578878.889-5.417-5.391-0.479979.444-5.932-5.907-0.42144-8-6-4-2024681012141618012345678910温度（℃）时间（h）分析解（℃）数值解（℃）图5平壁中心不同时刻的数值解和分析解8图6空间步长对墙内壁温度影响表2空间步长对温度影响数据时间(h)N=51N=101N=2010.00015.00015.00015.0000.50015.00015.00015.0001.00014.99814.99814.9981.50014.98014.98114.9812.00014.92414.92414.9242.50014.82014.82014.8203.00014.67514.67514.6763.50014.50114.50214.5024.00014.31014.31114.3114.50014.11114.11214.1125.00013.91113.91113.9125.50013.71513.71513.7156.00013.52513.52513.5256.50013.34313.34413.3447.00013.17113.17213.1727.50013.00913.01013.0108.00012.85812.85812.8588.50012.71612.71612.7169.00012.58312.58312.5839.50012.46012.46012.4609.91712.36312.36312.364（3）时间步长对温度的影响如图7和表3，根据图中曲线可以看出时间步长选择50s，100s，200s时基本重合，对墙内壁温度影响不大。12.012.513.013.514.014.515.015.501234567891011温度（℃）时间（h）N=51N=101N=2019图7时间步长对温度的影响表3时间步长对温度的影响时间(h)T=50（s）时间(h)T=100（s)时间（h）T=200（s）015015.0000150.694150.69415.0000.667151.38914.9881.38914.9871.33314.9872.08314.9122.08314.910214.9182.77814.7462.77814.7442.66714.773.47214.5133.47214.5113.33314.5584.16714.2454.16714.244414.3084.86113.9664.86113.9664.66714.0445.55613.6925.55613.6935.33313.7816.2513.4316.2513.433613.5276.94413.1886.94413.1906.66713.2887.63912.9647.63912.9667.33313.0668.33312.768.33312.762812.8629.02812.5749.02812.5768.66712.675图8墙内壁温度随时间的变化曲线1212.51313.51414.51515.50246810温度℃时间（h）T=50（s）T=100sT=200s101112131415160246810121416182022温度（℃）时间（h）105计算结果与分析5.1墙内壁温度分析根据题目中要求，计算寒潮入侵多长时间后内墙壁可以感受到外界气温的变化，通过建模，方程离散化，最终通过程序求解方程，得到图8和表4。由图可以看出，开始阶段，内墙壁温不变，随着时间的进一步深入，内壁温度开始降低，当很长时间后，温度变化基本趋于平缓，直到再次平衡。根据图8就可以得到墙内壁温度开始发生变化的时间。表4墙内壁温度随时间变化数据表时间（h)温度（℃）时间（h)温度（℃）时间（h)温度（℃）0.000156.66713.28513.33311.7630.556157.22213.09813.88911.6911.11114.9977.77812.92414.44411.6261.66714.9678.33312.76215.00011.5652.22214.8838.88912.61215.55611.512.77814.7449.44412.47316.11111.4583.33314.56210.00012.34516.66711.4113.88914.35410.55612.22717.22211.3684.44414.13311.11112.11817.77811.3295.00013.91111.66712.01818.33311.2925.55613.69312.22211.92618.88911.2596.11113.48412.77811.84119.44411.2285.2导热系数λ对墙内壁温度的影响墙的导热系数对内