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

计算方法期末大作业1一、问题分析已知量t0=20℃tf=200℃h1=50W(m2∗K)⁄KmWh*/522ρc=30∗105dx=0.01mλ=60W(m∗K)⁄导出量ca/)4/()(2adxdt2)/(dxdtaFo/dxhBi目标是借助matlab用数值分析的方法求得整块板达到最终稳态的过程中，这块板的温度随时间的变化，以及其中三个代表性点的温度随时间的变化，并通过调整不同的时间步长，空间步长来探索其对结果的影响，从而掌握数值方法的运用，进而验证传热学的理论知识，通过实践及思考获取新知识。由图可得，图形左右对称，即轴线位置属于绝热边界，所以我们可以把板沿轴线左右一分为二，只要计算一半的板即可，这使需要的计算的范围缩小了一半，降低了工作的难度，2也为后续的编程工作节省了计算的时间：要运用数值分析的方法，首先我们要将这块板的这块区域划分成许多互不重叠的单元体。每个单元体的物理量用某一点的物理量来代替，然后利用每个节点的物理量都与它周围的节点的物理量有一定关系，将物理现象描述成一系列代数方程组，再用计算机求解。首先，由于此题属于非稳态导热，选用x-y-t的三维坐标系。然后，将50*100的区域划分为50*100的节点，并列出相应节点的节点方程。最后，用矩阵法直接联解代数方程组。二、节点方程划分为50*100的网格：当(i,j)∈(i=0or100,0j50),(0i25or75i100,j=50)时，ti,j(k+1)=Fo(2ti−1,j(k)+ti,j−1(k)+ti,j+1(k))+(1−4Fo)ti,j(k)3当(i,j)∈(0i100,j=0),(25i75,j=25),(i=25or75,25j75)时，ti,j(k+1)=(1−4Fo−2Bi)ti,j(k)+Fo(ti,j−1(k)+ti,j+1(k)+2ti−1,j(k))+2FoBitf当(i,j)∈(0i25,0j50),(25i75,0j25),(75i100,0j50)时，ti,j(k+1)=Fo(ti,j−1(k)+ti,j+1(k)+ti−1,j(k)+ti+1,j(k))+(1−4Fo)ti,j(k)其中i,j∈Ζ，Fo=a∗dtx2⁄，Bi=hxλ⁄，a=λρc⁄，dt=x24a⁄三、编程思路把50*100的节点的温度当做一个50*100的温度矩阵。因为在非稳态显示格式中，大多数的节点的热平衡方程都是相同的。所以把这个矩阵的全部内容全部往四个方向移动一行或一列，获得四个新的矩阵。再根据热平衡方程，用这四个新的矩阵和原矩阵进行矩阵加减运算，就可以得到一个时间步长后各个节点的温度矩阵。在有边界的地方（M，N，P，Q）依然按照类似方法制造四个矩阵，但是有边界条件的方向使用边界条件给定的特殊值来创造矩阵。然后同样利用显式方程来计算一个时间步长后的温度矩阵。得到一般点和所有点的温度之后，再把各个矩阵中有用的部分拼接起来就得到了完整的一个步长之后的温度矩阵。利用循环语句反复进行上述步骤，直到判定收敛。程序代码见本报告附录。四、结果分析数据说明，h1为外部传热系数，h2为内部传热系数，Bo1，Bo2为相应的Bo数。Fo的特征长度的平方为微元的长与宽的乘积。t1为外部流体的温度，且衡为200摄氏度t2为内部流体和介质的初始温度，且衡为20摄氏度。质量与比热容的乘积为KmJ36/103，除了对时间步长的分析外，时间步长为x^2/(a*4)。1.导热与传热参数对结果的影响为了使数据有一般性，选取1mm*1mm的控制体，则温度场矩阵的阶数为100*50。外部流体温度为200摄氏度，内部流体与对象初始温度为20摄氏度。（1）选取一般性的常见参数，设导热系数为40w/m*K,内部空洞的传热系数为30W/2m*K，外部传热系数为50W/2m*K。计算可得Bi1=0.0125,Bi2=0.0075,Fo=0.000025。获得图像。（这一次运算条件作为标准计算，将会被多次提及）4T1=133.6058T2=80.6383T3=74.5586time=5.3573e+04（2）使导热系数大幅增大，设导热系数为40w/m*K,内部空洞的传热系数为3W/2m*K，外部传热系数为5W/2m*K。计算可得Bi1=0.00125,Bi2=0.00075。为了缩短计算时间，将距离步长改为2.5mm*2.5mm。获得图像。可见在Bi数较小的情况下，物体内的温度梯度较小，几乎可以近似地被认为是集总参数法问题。利用稳态条件下的计算，以由外部热流流入的热量与由内部冷流流出的热量相同为计算依据。)2(22)1(11TTShTTSh。可得物体平均温度应当近似于102摄氏5度。由图像可见与数值解相近。边界层的梯度非常大，可见在这种情况下，总的传热量最大。T1=108.2997T2=100.5442T3=99.3991time=3.8500e+05因为传热系数小，系统吸收热量的速度很慢，所以从初态到稳态花费的时间很多。（3）使用非常小的导热系数，设导热系数为4w/m*K,内部空洞的传热系数为30W/2m*K，外部传热系数为50W/2m*K。计算可得Bi1=0.125,Bi2=0.075。获得图像。因为导热热阻变大了，所以明显物体内所有地方的的整体温度梯度都相对标准计算图变大了。T1=166.2155T2=37.9005T3=31.9094time=1.4946e+05由数据可见因为a变小，热得扩散速度也变慢了。（4）使内部流体的传热系数大幅减小，设导热系数为40w/m*K,内部空洞的传热系数为3W/2m*K，外部传热系数为50W/2m*K。计算可得Bi1=0.0125,Bi2=0.00075,。获得图像。6易于理解，在内部流体的传热系数较小时，整个传热过程被外部流体所主导。内部流体需要更高的边界层温度梯度，来吸收外部流体所带来的热量。T1=172.1985T2=151.5862T3=148.6036time=1.0529e+05因为初始条件与稳态之间的差别较大，导致消耗了很多时间。（5）使外部流体的传热系数大幅减小，设导热系数为40w/m*K,内部空洞的传热系数为30W/2m*K，外部传热系数为5W/2m*K。计算可得Bi1=0.00125,Bi2=0.0075,。获得图像。7可见结果与上一种情况完全相反。T1=31.7660T2=23.6478T3=22.8033time=15540因为从初态到稳态的温度变化很小，所以花费的时间也最小。2．空间步长对结果的影响为了使结果具有一般性，使用在自然空间较为常见的导热和传热系数。设导热系数为40w/m*K,内部空洞的传热系数为30W/2m*K，外部传热系数为50W/2m*K，外部流体温度为200摄氏度，内部流体与对象初始温度为20摄氏度。为了便于阅读，将标准计算的数据列于此。T1=133.6058T2=80.6383T3=74.5586time=5.3573e+04（1）首先取空间步长为2.5*2.5mm，温度场矩阵的阶数为40*20。Fo=1.5625e-04T1=141.7098T2=90.4344T3=84.0126time=1.0257e+05（2）取空间步长为5*5mm，温度场矩阵的阶数为20*10。Fo=6.258T1=144.0982T2=92.0991T3=85.0609time=1.3781e+05由数据可见，往往当使用更加大的微元体时，过程可以进行地更深入。这不难理解,首先，Fo数是代表非稳态进程程度的无量纲时间数，Fo更大，过程进行地越深入。另一方面来讲，因为更长的步长，意味着更加长的时间步数。在相同的收敛判断条件下，由于步长更长的运算单步的温度变化更大，所以进行的程度更深。3.时间步长对结果的影响（1）.首先取时间步长dt=x^2/(a*4)=1.875，Fo=0.25，T1=133.6058T2=80.6383T3=74.5586time=5.3573e+004得到标准数据9(2).取时间步长dt=2，Fo=0.267，发散(3).取时间步长dt=1.7，Fo=0.227，T1=132.6748T2=79.4658T3=73.3294time=5.0849e+004当稳定性条件都满足才能保证整个方程组稳定，本题情况下需要保证1-4Fo≥0，当1-4Fo＜0的时候结果不稳定。满足1-4Fo≥0时，当x一定，dt越小，Fo越小，越能满足稳定性条件。但是相应的计算工作量增大。4.与稳态解的比较稳态过程的计算本质上与非稳态有相似之处，在此不加说明，使用代表稳态运算的m文件可以得到与第一次标准状态相同的条件下的解。T1=110.8181T2=79.4591T3=74.8974显然与标准运算有比较大的出入，原因不明。10附录1．data%设定节点长度，假设整个图形是一个1*0.5米的区域，x为单个微元的长，y为宽。%x与y取0.01m，0.025m，0.05m等，可以使长度方向微元数为4的倍数，宽度方向为2的倍数。t0=20tf=200h1=50h2=30pc=30e5x=0.1y=0.1gama=40%计算得到的常量a=gama/pcdt=x^2/(a*4)Fo=a*dt/x^2Bi1=h1*x/gamaBi2=h2*x/gama%辅助量，X为长度方向的微元数，Y为宽度方向上的微元数X=1/xY=0.5/yA=ones(Y,X)A=t0*A2.steady%只能在x=y=0.01m时使用k=1C4=100whileC40.001%先计算大多数点B=A(1:50,1:99);C=A(1:50,1:1);A1=[CB];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;11A3=[BD];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);12Q2=[C;z];Q3=N1;D=gama*(Q1+Q2+Q4)/(3*gama+h2*x)+(h2*x*Q3)/(3*gama+h2*x);Q=D;%最后把矩阵拆分再拼接X1=G(1:25,1:25);X2=G(1:24,26:75);X3=G(1:25,76:99);X4=G(26:50,1:24);C=ones(25,50);X5=20*C;X6=G(26:50,77:99);