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

离步长改为 2.5mm*2.5mm。获得图像。
可见在 Bi 数较小的情况下，物体内的温度梯度较小，几乎可以近似地被认为是集总参数法 问题。 利用稳态条件下的计算， 以由外部热流流入的热量与由内部冷流流出的热量相同为计 算依据。 h1 S1 (T 1 T ) h 2 S 2 (T T 2) 。可得物体平均温度应当近似于 102 摄氏度。
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 更大， 过程进行地越深入。 另一方面来讲， 因为更长的步长，意味着更加长的时间步数。在相同的收敛判断条件下，由于步长更长的运 算单步的温度变化更大，所以进行的程度更深。
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。为了缩短计算时间，将距
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，
其中 i , j ∈ Ζ，Fo = a ∗ dt x 2 ，Bi = hx λ，a = λ ρc，dt = x 2 4a
三、编程思路
把 50*100 的节点的温度当做一个 50*100 的温度矩阵。因为在非稳态显示格式中，大多 数的节点的热平衡方程都是相同的。 所以把这个矩阵的全部内容全部往四个方向移动一行或 一列，获得四个新的矩阵。再根据热平衡方程，用这四个新的矩阵和原矩阵进行矩阵加减运 算，就可以得到一个时间步长后各个节点的温度矩阵。 在有边界的地方（M，N，P，Q）依然按照类似方法制造四个矩阵，但是有边界条件的 方向使用边界条件给定的特殊值来创造矩阵。 然后同样利用显式方程来计算一个时间步长后 的温度矩阵。 得到一般点和所有点的温度之后， 再把各个矩阵中有用的部分拼接起来就得到了完整的一个 步长之后的温度矩阵。 利用循环语句反复进行上述步骤，直到判定收敛。 程序代码见本报告附录。
Bi h dx /
目标是借助 matlab 用数值分析的方法求得整块板达到最终稳态的过程中， 这块板的温度 随时间的变化，以及其中三个代表性点的温度随时间的变化，并通过调整不同的时间步长， 空间步长来探索其对结果的影响，从而掌握数值方法的运用，进而验证传热学的理论知识， 通过实践及思考获取新知识。 由图可得，图形左右对称，即轴线位置属于绝热边界，所以我们可以把板沿轴线左右一 分为二，只要计算一半的板即可，这使需要的计算的范围缩小了一半，降低了工作的难度，
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越小，越能满足稳定性 条件。但是相应的计算工作量增大。
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 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。获得图像。 （这一次运算条件作为标准计算，将会被多次提及）
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,数值解相近。边界层的梯度非常大，可见在这种情况下，总的传热量最大。 T1 =108.2997 T2 =100.5442 T3 =99.3991 time =3.8500e+05 因为传热系数小，系统吸收热量的速度很慢，所以从初态到稳态花费的时间很多。 （3）使用非常小的导热系数，设导热系数为 4w/m*K,内部空洞的传热系数为 30W/ m 2 *K，
计算方法期末大作业
一、问题分析
已知量 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
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);
1
也为后续的编程工作节省了计算的时间：
要运用数值分析的方法， 首先我们要将这块板的这块区域划分成许多互不重叠的单元体。 每个单元体的物理量用某一点的物理量来代替， 然后利用每个节点的物理量都与它周围的节 点的物理量有一定关系，将物理现象描述成一系列代数方程组，再用计算机求解。首先，由 于此题属于非稳态导热，选用 x-y-t 的三维坐标系。然后，将 50*100 的区域划分为 50*100 的节点，并列出相应节点的节点方程。最后，用矩阵法直接联解代数方程组。