计算传热学程序设计
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
图1墙壁简化图
1.1已知参数
壁厚，墙壁导热系数，密度与比热容的乘积，室内和寒潮入侵后室外空气温度，室内空气和外墙的表面传热系数，开始时稳定状态下的内墙表面温度。
1.2求解
寒潮入侵多少时间后内墙壁面可感受到外界气温的变化？
2物理与数学模型
2.1物理模型
该墙面为常物性，可以假设：（1）其为无限大平面，（2）只有在厚度方向传热，没有纵向传热，则该问题转化为一维常物性无限大平面非稳态导热问题。
采用外点法用均匀网格对求解区域进行离散化，得到的网格系统如图2所示。一共使用了0~N-1共N个节点。
节点间距δx为:
图2墙壁内的网格划分
此例中墙壁导热系数为常值，无源项。则可采用有限体积法对控制方程离散化，得到离散方程为：
（2a）
式中：
（2b）
， ， （2c）
（2d）
其中的上标“0”表示此为上一时刻的值,分别为节点所在控制容积左右边界上的导热系数，由于墙壁导热系数不变，故都等于λ，△τ为时间步长。由元体能量平衡法可以得知左右边界节点的离散方程分别为：
直接根据公式得到解析解如下：
（10）
式中， ，系数 应该使上述无穷级数在 是满足初始条件，由傅里叶级数理论可得：
（11）
是超越方程的根，称为特征根。
…（12）
其中 。
4.2程序验证
（1）由模型可以得到相关信息然后进行编程，同等时间下计算出中心处温度的解析解和数值解进行比较，数据记录在表1。然后计算出相对误差，作图5，观察数值解与分析解的比较曲线。
11.368
4.444
14.133
11.111
12.118
17.778
11.329
5.000
13.911
11.667
12.018
18.333
11.292
5.556
13.693
12.222
11.926
18.889
11.259
6.111
13.484
12.778
11.841
19.444
11.228
5.2导热系数λ对墙内壁温度的影响
表2空间步长对温度影响数据
时间(h)
N=51
N=101
N=201
0.000
15.000
15.000
15.000
0.500
15.000
15.000
15.000
1.000
14.998
14.998
14.998
1.500
14.980
14.981
14.981
2.000
14.924
14.924
14.924
2.500
墙的导热系数对内表面的影响，在图9和表5中发现，导热系数对内壁温度影响比较大，λ=1.2时，温度下降的趋势会更快，要比λ=0.85时下降快的多，下降速度更快，更短时间内达到稳态，因此导热系数越大温度扩散越快，导热系数越小则温度变化越慢，需要更长时间到达稳态，但是这时对于要求恒温的空间有好处，受波动影响更小。
13.715
6.000
13.525
13.525
13.525
6.500
13.343
13.344
13.344
7.000
13.171
13.172
13.172
7.500
13.009
13.010
13.010
8.000
12.858
12.858
12.858
8.500
12.716
12.716
12.716
9.000
12.583
（2）空间步长对墙内壁的温度影响如图6及表2。在程序编写过程中用网格节点数对空间步长进行控制，为了观察空间步长对墙内壁温度的影响，表中选择了三个不同的空间步长，分别为选取51，101，201个网格节点，则相应的空间步长为0.006，0.003，0.0015。根据不同步长时温度的变化曲线可以看出，空间步长对内墙壁的影响不大，当空间步长控制在合理范围时可以忽略空间步长的影响。
14.266
13.808
3.889
14.354
13.968
13.398
4.444
14.133
13.666
12.999
5
13.911
13.371
12.62
5.556
13.693
13.088
12.267
6.111
13.484
12.823
11.942
6.667
13.285
12.575
11.644
7.222
8.333
12.762
15.000
11.565
2.222
14.883
8.889
12.612
15.556
11.51
2.778
14.744
9.444
12.473
16.111
11.458
3.333
14.562
10.000
12.345
16.667
11.411
3.889
14.354
10.556
12.227
17.222
中国石油大学(华东)
储运与建筑工程学院热能与动力工程系
《计算传热学程序设计》
设计报告
学生姓名:
学号：
专业班级：
指导教师
2012年7月7日
1、设计题目
有一房屋的砖墙厚δ=0.3m，λ=0.85W/(m·℃)，ρc=1.05×106J/(m3·K)，室内温度Tf1保持20℃不变，表面传热系数h1=6W/(m2·℃)。开始时墙的温度处于稳定状态，内墙表面温度Tw1为15℃寒潮入侵后，室外温度Tf2下降为-10℃，外墙的表面传热系数为35W/(m2·℃)。试分析寒潮入侵后多少时间内墙壁面方可感受到外界气温的变化。
2.2数学模型
以墙外表面为坐标原点，沿厚度方向为坐标正方向，建立坐标系。基于上述模型，取其在x方向上的微元作为研究对象，则该问题的数学模型可描述如下：
（1a）
初始条件：
（1b）
在两侧相应的边界条件是第三类边界条件，分别由傅立叶定律可描述如下：
左边界：
（1c）
右边界：
（1d）
3数值处理与程序设计
3.1数值处理
13.433
6
13.527
6.944
13.188
6.944
13.190
6.667
13.288
7.639
12.964
7.639
12.966
7.333
13.066
8.333
12.76
8.333
12.762
8
12.862
9.028
12.574
9.028
12.576
8.667
12.675
图8墙内壁温度随时间的变化曲线
5计算结果与分析
5.1墙内壁温度分析
根据题目中要求，计算寒潮入侵多长时间后内墙壁可以感受到外界气温的变化，通过建模，方程离散化，最终通过程序求解方程，得到图8和表4。由图可以看出，开始阶段，内墙壁温不变，随着时间的进一步深入，内壁温度开始降低，当很长时间后，温度变化基本趋于平缓，直到再次平衡。根据图8就可以得到墙内壁温度开始发生变化的时间。
-2.618
-2.586
-1.22231
7.222
-3.447
-3.417
-0.87032
7.778
-4.184
-4.154
-0.71702
8.333
-4.837
-4.809
-0.57887
8.889
-5.417
-5.391
-0.47997
9.444
-5.932
-5.907
-0.42144
图6空间步长对墙内壁温度影响
程序特点：该程序有很强的适应性，一维常物性非稳态平壁导热问题都可以使用此程序，只要适当更改边值条件即可。还可以进行修改解决非常物性问题。
程序中对输出节点，最大输出量都进行了控制，对计算结果的分析有很大帮助。而且Thoms算法的优点需要内存小，工作量小，程序设计简单。
程序流程图：首先对变量赋值，然后由初始条件建立初始温度场，接着从左边界，内部节点，到右边界进行迭代，直到满足精度要求为止，最后输出结果，程序结束。程序流程如下图3。
0.358474
3.333
5.082
5.112
0.590319
3.889
3.393
3.425
0.943118
4.444
Leabharlann Baidu1.89
1.923
1.746032
5．00
0.554
0.588
6.137184
5.556
-0.631
-0.598
-5.22979
6.111
-1.684
-1.651
-1.95962
6.667
T=50（s）
时间(h)
T=100（s)
时间（h）
T=200（s）
0
15
0
15.000
0
15
0.694
15
0.694
15.000
0.667
15
1.389
14.988
1.389
14.987
1.333
14.987
2.083
14.912
2.083
14.910
2
14.918
2.778
14.746
2.778
14.744
左边界节点：
（3）
右边界节点：
（4）
离散方程的详细推导过程见附录。
3.2程序设计
由物理模型可以知道本问题为一维导热问题，一维导热问题的离散方程在取遍所有节点后形成的是三对角的代数方程组，采用追赶法进行求解。
程序构成和方法:程序由主程序和一个子程序构成。主程序进行变量定义和各已知参数的输入，以及左右边界节点和内部节点控制方程的输入；子程序tdma实现追赶法用来计算每个节点新的温度。Thomas算法求解过程分为两步:消元和回代。消元是从系数矩阵的第二行起，逐一将每一行的非零元素消去一个，使原来的三元方程化为二元方程。消元进行到最后一行时，二元方程就化为一元方程，直接得到最后一个未知数的值。然后逐一往前回代，由各二元方程求出其它未知解。
表4墙内壁温度随时间变化数据表
时间（h)
温度（℃）
时间（h)
温度（℃）
时间（h)
温度（℃）
0.000
15
6.667
13.285
13.333
11.763
0.556
15
7.222
13.098
13.889
11.691
1.111
14.997
7.778
12.924
14.444
11.626
1.667
14.967
2.667
14.77
3.472
14.513
3.472
14.511
3.333
14.558
4.167
14.245
4.167
14.244
4
14.308
4.861
13.966
4.861
13.966
4.667
14.044
5.556
13.692
5.556
13.693
5.333
13.781
6.25
13.431
6.25
13.098
12.346
11.372
7.778
12.924
12.136
11.125
8.333
12.762
11.942
10.901
8.889
12.612
11.764
10.698
9.444
12.473
11.602
10.514
图9导热系数对墙内壁温度的影响
5.3墙外换热系数h的影响
墙外表面传热系数对温度分布的影响，如图10和表6影响不大。
14.820
14.820
14.820
3.000
14.675
14.675
14.676
3.500
14.501
14.502
14.502
4.000
14.310
14.311
14.311
4.500
14.111
14.112
14.112
5.000
13.911
13.911
13.912
5.500
13.715
13.715
表5导热系数对墙内壁温度的影响
时刻（h）
λ=0.85
λ=1.0
λ=1.2
0
15
15
15
0.556
15
15
15
1.111
14.997
14.99
14.971
1.667
14.967
14.924
14.833
2.222
14.883
14.77
14.564
2.778
14.744
14.542
14.207
3.333
14.562
4、模型与程序验证
4.1模型
本题简化为厚度为2 =0.3m的一维非稳态模型如图4所示，初始温度为15℃，在其中间建立坐标系，左两边为对流换热，且换热系数相同都为h=25 W/(m2·℃)，且流体温度Tf=-10℃对于x 0，列出其导热微分方程式及定解条件：
（5）
（6）
（7）
（8）
引入过余温度：
（9）
由图表中可以发现，平壁中心不同时刻温度值的分析解和数值解相差不是很大，二者吻合的比较好，可以说明所编制的数值解法的程序是正确的。相对误差先增大后减小，增大的原因是此时温度接近零度，相对误差的基数比较小，所以造成相对误差较大，但是此时的绝对误差并不大，在合理范围内，所以除去个别点外，都满足误差小于百分之1。可以验证所编数值解法的程序是正确的。
12.583
12.583
9.500
12.460
12.460
12.460
9.917
12.363
12.363
12.364
（3）时间步长对温度的影响如图7和表3，根据图中曲线可以看出时间步长选择50s，100s，200s时基本重合，对墙内壁温度影响不大。
图7时间步长对温度的影响
表3时间步长对温度的影响
时间(h)
图10墙外换热系数对温度影响
表6墙外换热系数对温度影响
时刻（h）
h=20W/(m2·℃)
h=35
W/(m2·℃)
h=50
W/(m2·℃)
0
15
15
15
1.111
14.998
14.997
14.996
2.222
14.914
14.883
14.865
3.333
14.662
14.562
14.507
4.444
14.311
表1分析解与数值解比较
时间（h）
分析解（℃）
数值解（℃）
相对误差（%）
0
15
15
0
0.556
14.854
14.811
-0.28948
1.111
13.463
13.425
-0.28226
1.667
11.306
11.3
-0.05307
2.222
9.066
9.081
0.165453
2.778
6.974
6.999
1.1已知参数
壁厚，墙壁导热系数，密度与比热容的乘积，室内和寒潮入侵后室外空气温度，室内空气和外墙的表面传热系数，开始时稳定状态下的内墙表面温度。
1.2求解
寒潮入侵多少时间后内墙壁面可感受到外界气温的变化？
2物理与数学模型
2.1物理模型
该墙面为常物性，可以假设：（1）其为无限大平面，（2）只有在厚度方向传热，没有纵向传热，则该问题转化为一维常物性无限大平面非稳态导热问题。
采用外点法用均匀网格对求解区域进行离散化，得到的网格系统如图2所示。一共使用了0~N-1共N个节点。
节点间距δx为:
图2墙壁内的网格划分
此例中墙壁导热系数为常值，无源项。则可采用有限体积法对控制方程离散化，得到离散方程为：
（2a）
式中：
（2b）
， ， （2c）
（2d）
其中的上标“0”表示此为上一时刻的值,分别为节点所在控制容积左右边界上的导热系数，由于墙壁导热系数不变，故都等于λ，△τ为时间步长。由元体能量平衡法可以得知左右边界节点的离散方程分别为：
直接根据公式得到解析解如下：
（10）
式中， ，系数 应该使上述无穷级数在 是满足初始条件，由傅里叶级数理论可得：
（11）
是超越方程的根，称为特征根。
…（12）
其中 。
4.2程序验证
（1）由模型可以得到相关信息然后进行编程，同等时间下计算出中心处温度的解析解和数值解进行比较，数据记录在表1。然后计算出相对误差，作图5，观察数值解与分析解的比较曲线。
11.368
4.444
14.133
11.111
12.118
17.778
11.329
5.000
13.911
11.667
12.018
18.333
11.292
5.556
13.693
12.222
11.926
18.889
11.259
6.111
13.484
12.778
11.841
19.444
11.228
5.2导热系数λ对墙内壁温度的影响
表2空间步长对温度影响数据
时间(h)
N=51
N=101
N=201
0.000
15.000
15.000
15.000
0.500
15.000
15.000
15.000
1.000
14.998
14.998
14.998
1.500
14.980
14.981
14.981
2.000
14.924
14.924
14.924
2.500
墙的导热系数对内表面的影响，在图9和表5中发现，导热系数对内壁温度影响比较大，λ=1.2时，温度下降的趋势会更快，要比λ=0.85时下降快的多，下降速度更快，更短时间内达到稳态，因此导热系数越大温度扩散越快，导热系数越小则温度变化越慢，需要更长时间到达稳态，但是这时对于要求恒温的空间有好处，受波动影响更小。
13.715
6.000
13.525
13.525
13.525
6.500
13.343
13.344
13.344
7.000
13.171
13.172
13.172
7.500
13.009
13.010
13.010
8.000
12.858
12.858
12.858
8.500
12.716
12.716
12.716
9.000
12.583
（2）空间步长对墙内壁的温度影响如图6及表2。在程序编写过程中用网格节点数对空间步长进行控制，为了观察空间步长对墙内壁温度的影响，表中选择了三个不同的空间步长，分别为选取51，101，201个网格节点，则相应的空间步长为0.006，0.003，0.0015。根据不同步长时温度的变化曲线可以看出，空间步长对内墙壁的影响不大，当空间步长控制在合理范围时可以忽略空间步长的影响。
14.266
13.808
3.889
14.354
13.968
13.398
4.444
14.133
13.666
12.999
5
13.911
13.371
12.62
5.556
13.693
13.088
12.267
6.111
13.484
12.823
11.942
6.667
13.285
12.575
11.644
7.222
8.333
12.762
15.000
11.565
2.222
14.883
8.889
12.612
15.556
11.51
2.778
14.744
9.444
12.473
16.111
11.458
3.333
14.562
10.000
12.345
16.667
11.411
3.889
14.354
10.556
12.227
17.222
中国石油大学(华东)
储运与建筑工程学院热能与动力工程系
《计算传热学程序设计》
设计报告
学生姓名:
学号：
专业班级：
指导教师
2012年7月7日
1、设计题目
有一房屋的砖墙厚δ=0.3m，λ=0.85W/(m·℃)，ρc=1.05×106J/(m3·K)，室内温度Tf1保持20℃不变，表面传热系数h1=6W/(m2·℃)。开始时墙的温度处于稳定状态，内墙表面温度Tw1为15℃寒潮入侵后，室外温度Tf2下降为-10℃，外墙的表面传热系数为35W/(m2·℃)。试分析寒潮入侵后多少时间内墙壁面方可感受到外界气温的变化。
2.2数学模型
以墙外表面为坐标原点，沿厚度方向为坐标正方向，建立坐标系。基于上述模型，取其在x方向上的微元作为研究对象，则该问题的数学模型可描述如下：
（1a）
初始条件：
（1b）
在两侧相应的边界条件是第三类边界条件，分别由傅立叶定律可描述如下：
左边界：
（1c）
右边界：
（1d）
3数值处理与程序设计
3.1数值处理
13.433
6
13.527
6.944
13.188
6.944
13.190
6.667
13.288
7.639
12.964
7.639
12.966
7.333
13.066
8.333
12.76
8.333
12.762
8
12.862
9.028
12.574
9.028
12.576
8.667
12.675
图8墙内壁温度随时间的变化曲线
5计算结果与分析
5.1墙内壁温度分析
根据题目中要求，计算寒潮入侵多长时间后内墙壁可以感受到外界气温的变化，通过建模，方程离散化，最终通过程序求解方程，得到图8和表4。由图可以看出，开始阶段，内墙壁温不变，随着时间的进一步深入，内壁温度开始降低，当很长时间后，温度变化基本趋于平缓，直到再次平衡。根据图8就可以得到墙内壁温度开始发生变化的时间。
-2.618
-2.586
-1.22231
7.222
-3.447
-3.417
-0.87032
7.778
-4.184
-4.154
-0.71702
8.333
-4.837
-4.809
-0.57887
8.889
-5.417
-5.391
-0.47997
9.444
-5.932
-5.907
-0.42144
图6空间步长对墙内壁温度影响
程序特点：该程序有很强的适应性，一维常物性非稳态平壁导热问题都可以使用此程序，只要适当更改边值条件即可。还可以进行修改解决非常物性问题。
程序中对输出节点，最大输出量都进行了控制，对计算结果的分析有很大帮助。而且Thoms算法的优点需要内存小，工作量小，程序设计简单。
程序流程图：首先对变量赋值，然后由初始条件建立初始温度场，接着从左边界，内部节点，到右边界进行迭代，直到满足精度要求为止，最后输出结果，程序结束。程序流程如下图3。
0.358474
3.333
5.082
5.112
0.590319
3.889
3.393
3.425
0.943118
4.444
Leabharlann Baidu1.89
1.923
1.746032
5．00
0.554
0.588
6.137184
5.556
-0.631
-0.598
-5.22979
6.111
-1.684
-1.651
-1.95962
6.667
T=50（s）
时间(h)
T=100（s)
时间（h）
T=200（s）
0
15
0
15.000
0
15
0.694
15
0.694
15.000
0.667
15
1.389
14.988
1.389
14.987
1.333
14.987
2.083
14.912
2.083
14.910
2
14.918
2.778
14.746
2.778
14.744
左边界节点：
（3）
右边界节点：
（4）
离散方程的详细推导过程见附录。
3.2程序设计
由物理模型可以知道本问题为一维导热问题，一维导热问题的离散方程在取遍所有节点后形成的是三对角的代数方程组，采用追赶法进行求解。
程序构成和方法:程序由主程序和一个子程序构成。主程序进行变量定义和各已知参数的输入，以及左右边界节点和内部节点控制方程的输入；子程序tdma实现追赶法用来计算每个节点新的温度。Thomas算法求解过程分为两步:消元和回代。消元是从系数矩阵的第二行起，逐一将每一行的非零元素消去一个，使原来的三元方程化为二元方程。消元进行到最后一行时，二元方程就化为一元方程，直接得到最后一个未知数的值。然后逐一往前回代，由各二元方程求出其它未知解。
表4墙内壁温度随时间变化数据表
时间（h)
温度（℃）
时间（h)
温度（℃）
时间（h)
温度（℃）
0.000
15
6.667
13.285
13.333
11.763
0.556
15
7.222
13.098
13.889
11.691
1.111
14.997
7.778
12.924
14.444
11.626
1.667
14.967
2.667
14.77
3.472
14.513
3.472
14.511
3.333
14.558
4.167
14.245
4.167
14.244
4
14.308
4.861
13.966
4.861
13.966
4.667
14.044
5.556
13.692
5.556
13.693
5.333
13.781
6.25
13.431
6.25
13.098
12.346
11.372
7.778
12.924
12.136
11.125
8.333
12.762
11.942
10.901
8.889
12.612
11.764
10.698
9.444
12.473
11.602
10.514
图9导热系数对墙内壁温度的影响
5.3墙外换热系数h的影响
墙外表面传热系数对温度分布的影响，如图10和表6影响不大。
14.820
14.820
14.820
3.000
14.675
14.675
14.676
3.500
14.501
14.502
14.502
4.000
14.310
14.311
14.311
4.500
14.111
14.112
14.112
5.000
13.911
13.911
13.912
5.500
13.715
13.715
表5导热系数对墙内壁温度的影响
时刻（h）
λ=0.85
λ=1.0
λ=1.2
0
15
15
15
0.556
15
15
15
1.111
14.997
14.99
14.971
1.667
14.967
14.924
14.833
2.222
14.883
14.77
14.564
2.778
14.744
14.542
14.207
3.333
14.562
4、模型与程序验证
4.1模型
本题简化为厚度为2 =0.3m的一维非稳态模型如图4所示，初始温度为15℃，在其中间建立坐标系，左两边为对流换热，且换热系数相同都为h=25 W/(m2·℃)，且流体温度Tf=-10℃对于x 0，列出其导热微分方程式及定解条件：
（5）
（6）
（7）
（8）
引入过余温度：
（9）
由图表中可以发现，平壁中心不同时刻温度值的分析解和数值解相差不是很大，二者吻合的比较好，可以说明所编制的数值解法的程序是正确的。相对误差先增大后减小，增大的原因是此时温度接近零度，相对误差的基数比较小，所以造成相对误差较大，但是此时的绝对误差并不大，在合理范围内，所以除去个别点外，都满足误差小于百分之1。可以验证所编数值解法的程序是正确的。
12.583
12.583
9.500
12.460
12.460
12.460
9.917
12.363
12.363
12.364
（3）时间步长对温度的影响如图7和表3，根据图中曲线可以看出时间步长选择50s，100s，200s时基本重合，对墙内壁温度影响不大。
图7时间步长对温度的影响
表3时间步长对温度的影响
时间(h)
图10墙外换热系数对温度影响
表6墙外换热系数对温度影响
时刻（h）
h=20W/(m2·℃)
h=35
W/(m2·℃)
h=50
W/(m2·℃)
0
15
15
15
1.111
14.998
14.997
14.996
2.222
14.914
14.883
14.865
3.333
14.662
14.562
14.507
4.444
14.311
表1分析解与数值解比较
时间（h）
分析解（℃）
数值解（℃）
相对误差（%）
0
15
15
0
0.556
14.854
14.811
-0.28948
1.111
13.463
13.425
-0.28226
1.667
11.306
11.3
-0.05307
2.222
9.066
9.081
0.165453
2.778
6.974
6.999