上海交通大学计算方法大作业

2
3
T1 =133.6058 T2 =80.6383 T3 =74.5586 time =5.3573e+04 设导热系数为 40w/m*K,内部空洞的传热系数为 3W/ m 2 *K， （2）使导热系数大幅增大，
2 外部传热系数为 5W/ m *K。计算可得 Bi1=0.00125,Bi2=0.00075。为了缩短计算时间，将距
2
当(i , j)∈(0<i<100 , j=0) , (25<i<75 , j=25) , (i=25 or 75 , 25<j<75)时，
k k k t i k+1 = 1 − 4Fo − 2Bi t i k ,j ,j + Fo t i ,j −1 + t i ,j+1 + 2t i −1 ,j + 2Fo Bi t f
Bi h dx /
目标是借助 matlab 用数值分析的方法求得整块板达到最终稳态的过程中， 这块板的温度 随时间的变化，以及其中三个代表性点的温度随时间的变化，并通过调整不同的时间步长， 空间步长来探索其对结果的影响，从而掌握数值方法的运用，进而验证传热学的理论知识， 通过实践及思考获取新知识。 由图可得，图形左右对称，即轴线位置属于绝热边界，所以我们可以把板沿轴线左右一 分为二，只要计算一半的板即可，这使需要的计算的范围缩小了一半，降低了工作的难度，
3.时间步长对结果的影响
（1）.首先取时间步长dt=x^2/(a*4)=1.875，Fo=0.25，
T1 =133.6058 T2 =80.6383 T3 =74.5586 time =5.3573e+004得到标准数据
8
(2).取时间步长dt=2，Fo=0.267，发散
(3).取时间步长dt=1.7，Fo=0.227，
T1 =141.7098 T2 =90.4344 T3 =84.0126 time =1.0257e+05 （2）取空间步长为 5*5mm，温度场矩阵的阶数为 20*10。Fo=6.25
7
T1 =144.0982 T2 =92.0991 T3 =85.0609 time =1.3781e+05 由数据可见，往往当使用更加大的微元体时，过程可以进行地更深入。这不难理解,首先， Fo 数是代表非稳态进程程度的无量纲时间数， Fo 更大， 过程进行地越深入。 另一方面来讲， 因为更长的步长，意味着更加长的时间步数。在相同的收敛判断条件下，由于步长更长的运 算单步的温度变化更大，所以进行的程度更深。
二、节点方程
划分为 50*100 的网格： 当(i , j)∈(i=0 or 100 , 0<j<50) , (0<i<25 or 75<i<100 , j=50)时，
k k k k t i k+1 = Fo 2t i− ,j 1,j + t i,j −1 + t i,j+1 + 1 − 4Fo t i,j
其中 i , j ∈ Ζ，Fo = a ∗ dt x 2 ，Bi = hx λ，a = λ ρc，dt = x 2 4a
三、编程思路
把 50*100 的节点的温度当做一个 50*100 的温度矩阵。因为在非稳态显示格式中，大多 数的节点的热平衡方程都是相同的。 所以把这个矩阵的全部内容全部往四个方向移动一行或 一列，获得四个新的矩阵。再根据热平衡方程，用这四个新的矩阵和原矩阵进行矩阵加减运 算，就可以得到一个时间步长后各个节点的温度矩阵。 在有边界的地方（M，N，P，Q）依然按照类似方法制造四个矩阵，但是有边界条件的 方向使用边界条件给定的特殊值来创造矩阵。 然后同样利用显式方程来计算一个时间步长后 的温度矩阵。 得到一般点和所有点的温度之后， 再把各个矩阵中有用的部分拼接起来就得到了完整的一个 步长之后的温度矩阵。 利用循环语句反复进行上述步骤，直到判定收敛。 程序代码见本报告附录。
10
A3=[B D]; B=A(1:49,1:100); C=A(1:1,1:100); A4=[C;B]; G=(A1+A2+A3+A4)/4; %然后计算有对流换热的边界的点，首先计算M边 M=A(1:50,100:100); M1=A(1:50,99:99); C=A(2:50,100:100); z=A(50,100); M2=[C;z]; C=A(1:49,100:100); z=A(1,100); M4=[z;C]; C=ones(50,1); M3=tf*C; D=gama*(M1+M2+M4)/(3*gama+h1*x)+(h1*x*M3)/(3*gama+h1*x); M=D; %然后计算N边 N=A(26:50,76:76); N3=A(26:50,77:77); N4=A(25:49,76:76); C=A(27:50,76:76); z=A(50,76); N2=[C;z]; C=ones(25,1); N1=t0*C; D=gama*(N3+N2+N4)/(3*gama+h2*x)+(h2*x*N1)/(3*gama+h2*x); N=D; %计算P边 P=A(25:25,26:75); P4=A(24:24,26:75); P3=A(25:25,27:76); P1=A(25:25,25:74); C=ones(1,50); P2=t0*C; D=gama*(P1+P3+P4)/(3*gama+h2*x)+(h2*x*P2)/(3*gama+h2*x); P=D; %计算Q边 Q=A(26:50,25:25); Q1=A(26:50,24:24); Q4=A(25:49,25:25); C=A(27:50,25:25); z=A(25,50);
4
由图像可见与数值解相近。边界层的梯度非常大，可见在这种情况下，总的传热量最大。 T1 =108.2997 T2 =100.5442 T3 =99.3991 time =3.8500e+05 因为传热系数小，系统吸收热量的速度很慢，所以从初态到稳态花费的时间很多。 （3）使用非常小的导热系数，设导热系数为 4w/m*K,内部空洞的传热系数为 30W/ m 2 *K，
2 2
6
可见结果与上一种情况完全相反。 T1 =31.7660 T2 =23.6478 T3 =22.8033 time =15540 因为从初态到稳态的温度变化很小，所以花费的时间也最小。
2．空间步长对结果的影响
为了使结果具有一般性，使用在自然空间较为常见的导热和传热系数。设导热系数为 40w/m*K,内部空洞的传热系数为 30W/ m 2 *K，外部传热系数为 50W/ m 2 *K，外部流体温度 为 200 摄氏度，内部流体与对象初始温度为 20 摄氏度。 为了便于阅读，将标准计算的数据列于此。 T1 =133.6058 T2 =80.6383 T3 =74.5586 time =5.3573e+04 （1）首先取空间步长为 2.5*2.5mm，温度场矩阵的阶数为 40*20。Fo=1.5625e-04
离步长改为 2.5mm*2.5mm。获得图像。
Biblioteka Baidu
可见在 Bi 数较小的情况下，物体内的温度梯度较小，几乎可以近似地被认为是集总参数法 问题。 利用稳态条件下的计算， 以由外部热流流入的热量与由内部冷流流出的热量相同为计 算依据。 h1 S1 (T 1 T ) h 2 S 2 (T T 2) 。可得物体平均温度应当近似于 102 摄氏度。
2.steady
%只能在x=y=0.01m时使用 k=1 C4=100 while C4>0.001 %先计算大多数点 B=A(1:50,1:99); C=A(1:50,1:1); A1=[C B]; B=A(2:50,1:100); C=A(50:50,1:100); A2=[B;C]; B=A(1:50,2:100); C=ones(50,1); D=tf*C;
2 2
5
易于理解，在内部流体的传热系数较小时，整个传热过程被外部流体所主导。内部流体需要 更高的边界层温度梯度，来吸收外部流体所带来的热量。 T1 =172.1985 T2 =151.5862 T3 =148.6036 time = 1.0529e+05 因为初始条件与稳态之间的差别较大，导致消耗了很多时间。 设导热系数为 40w/m*K,内部空洞的传热系数为 （5）使外部流体的传热系数大幅减小， 30W/ m *K，外部传热系数为 5W/ m *K。计算可得 Bi1=0.00125,Bi2=0.0075,。获得图像。
1.导热与传热参数对结果的影响
为了使数据有一般性，选取 1mm*1mm 的控制体，则温度场矩阵的阶数为 100*50。外 部流体温度为 200 摄氏度，内部流体与对象初始温度为 20 摄氏度。 设导热系数为 40w/m*K,内部空洞的传热系数为 30W/ m 2 （1）选取一般性的常见参数， *K，外部传热系数为 50W/ m *K。计算可得 Bi1=0.0125,Bi2=0.0075,Fo=0.000025。获得图像。 （这一次运算条件作为标准计算，将会被多次提及）
2 外部传热系数为 50W/ m *K。计算可得 Bi1=0.125,Bi2=0.075。获得图像。
因为导热热阻变大了， 所以明显物体内所有地方的的整体温度梯度都相对标准计算图变大了。 T1 =166.2155 T2 =37.9005 T3 =31.9094 time =1.4946e+05 由数据可见因为 a 变小，热得扩散速度也变慢了。 设导热系数为 40w/m*K,内部空洞的传热系数为 （4）使内部流体的传热系数大幅减小， 3W/ m *K，外部传热系数为 50W/ m *K。计算可得 Bi1=0.0125,Bi2=0.00075,。获得图像。
计算方法期末大作业
一、问题分析
已知量 t 0 = 20℃ t f = 200℃ h1 = 50W m2 ∗ K
h2 5W / m2 * K
ρc = 30 ∗ 105 dx = 0.01m λ = 60W m ∗ K 导出量


a / c
dt (dx) 2 /(a 4)
Fo a dt /(dx) 2
1
也为后续的编程工作节省了计算的时间：
要运用数值分析的方法， 首先我们要将这块板的这块区域划分成许多互不重叠的单元体。 每个单元体的物理量用某一点的物理量来代替， 然后利用每个节点的物理量都与它周围的节 点的物理量有一定关系，将物理现象描述成一系列代数方程组，再用计算机求解。首先，由 于此题属于非稳态导热，选用 x-y-t 的三维坐标系。然后，将 50*100 的区域划分为 50*100 的节点，并列出相应节点的节点方程。最后，用矩阵法直接联解代数方程组。
四、结果分析
数据说明，h1为外部传热系数，h2为内部传热系数，Bo1，Bo2为相应的Bo数。Fo的特征长度的平方为 微元的长与宽的乘积。t1为外部流体的温度，且衡为200摄氏度t2为内部流体和介质的初始温度，且衡为20 摄氏度。 质量与比热容的乘积为 3 10
6
除了对时间步长的分析外， 时间步长为x^2/(a*4)。 J / m3 K ，
T1 =132.6748 T2 =79.4658 T3 =73.3294 time =5.0849e+004 当稳定性条件都满足才能保证整个方程组稳定，本题情况下需要保证1-4Fo≥0，当1-4Fo ＜0的时候结果不稳定。满足1-4Fo≥0时，当x一定，dt越小，Fo越小，越能满足稳定性 条件。但是相应的计算工作量增大。
当(i , j)∈(0<i<25 , 0<j<50) , (25<i<75 , 0<j<25) , (75<i<100 , 0<j<50)时，
k k k k t i k+1 = Fo t i k ,j ,j −1 + t i ,j+1 + t i −1 ,j + t i+1 ,j + 1 − 4Fo t i ,j
4.与稳态解的比较
稳态过程的计算本质上与非稳态有相似之处，在此不加说明，使用代表稳态运算的 m 文件可以得到与第一次标准状态相同的条件下的解。 T1 =110.8181 T2 =79.4591 T3 =74.8974 显然与标准运算有比较大的出入，原因不明。
9
附录
1．data
%设定节点长度，假设整个图形是一个 1*0.5 米的区域，x 为单个微元的长，y 为宽。 %x 与 y 取 0.01m，0.025m，0.05m 等，可以使长度方向微元数为 4 的倍数，宽度方向为 2 的倍 数。 t0=20 tf=200 h1 = 50 h2 = 30 pc = 30e5 x = 0.1 y = 0.1 gama = 40 %计算得到的常量 a = gama/pc dt = x^2/(a*4) Fo = a*dt/x^2 Bi1 = h1*x/gama Bi2 = h2*x/gama %辅助量，X为长度方向的微元数，Y为宽度方向上的微元数 X = 1/x Y = 0.5/y A =ones(Y,X) A=t0 * A