内科大冶金传输原理课程程序设计上机指导及说明书

冶金传输原理课程程序设计上机指导及说明书冶金传输原理数值计算部分共计8学时，包含有：1.速度边界层计算，2学时2.二维稳态导热计算，2学时3.一维非稳态导热计算（隐式差分方程），2学时4.二维非稳态导热计算（显式差分方程），2学时边界层微分方程数值解法1、目的和意义通过该程序设计，让学生掌握常微分方程初值问题数值解的求解方法及程序设计。2、方程的推导根据普朗特微分方程和布拉休斯解，可以得到如下方程： 运用常微分方程的数值解法进行求解。为了方便，在推导时用函数y代替，自变量x代替，将上述方程改写为：令，则上列方程可以变为：利用四阶标准Longe-Kutta法，可以写成一下的计算式。其中： 3、程序设计由于该问题不完全满足常微分方程初值问题的要求。缺少的值，但是多一个的值。需要通过间接求。程序设计上采用选定区间，然后用二分法来寻找。具体程序见教科书119页（Basic语言）、120页（Fortran语言）。4、结果分析计算结果见教科书122页。结果与豪沃斯计算表一致，而且精度更高一些。

二、二维稳态导热问题程序设计1、目的和意义通过该上机实践，让学生初步掌握偏微分方程数值解的基本方法，以及迭代法求解线性方程组的基本方法。2、内部节点差分方程的建立 3、边界节点的差分方程的建立1）第三类边界条件（左边界）2）第二类边界条件（上边界）3）第一类边界条件（右边界）当Bi时，第三类边界条件转化为第一类边界条件。即4、C语言程序#include<stdio.h>#include<math.h>main(){intI,J;intNI=52,NJ=52;intNIM1=NI-1,NJM1=NJ-1;intTF=1000,T0=100,NITER=0,NUM;floatLAMD,A0=20,ALFA=200,B=0.01,EPS=0.01,T2;floatX[52],Y[52];floatSEW[52],SNS[52],DXEP[52],DXPW[52],DYNP[52],DYPS[52];floatAE[52][52],AW[52][52],AS[52][52],AN[52][52],AP[52][52];floatT[52][52],T1[52][52],SU[52][52],SP[52][52],LAM[52][52];floatDX1=0.05,DX2=0.03,DX3=0.01,DY1=0.01,DY2=0.05;FILE*fp1,*fp2;fp1=fopen("d:\y1","w");fp2=fopen("d:\y2","w");X[1]=-DX1/2;for(I=2;I<=10;I++)X[I]=X[I-1]+DX1;for(I=11;I<=18;I++)X[I]=X[I-1]+DX2;for(I=19;I<=NI;I++)X[I]=X[I-1]+DX3;Y[1]=-DY1/2;for(J=2;J<=5;J++)Y[J]=Y[J-1]+DY1;for(J=6;J<=NJ;J++)Y[J]=Y[J-1]+DY2;fprintf(fp1,"%d,%d",NI,NJ);for(I=1;I<=NI;I++)fprintf(fp1,"%7.4f",X[I]);for(J=1;J<=NJ;J++)fprintf(fp1,"%7.4f",Y[J]);DXPW[1]=0.0;DXEP[NI]=0.0;for(I=1;I<=NIM1;I++){DXEP[I]=X[I+1]-X[I];DXPW[I+1]=DXEP[I];}DYPS[1]=0.0;DYNP[NJ]=0.0;for(J=1;J<=NJM1;J++){DYNP[J]=Y[J+1]-Y[J];DYPS[J+1]=DYNP[J];}SEW[1]=0.0;SEW[NI]=0.0;for(I=2;I<=NIM1;I++)SEW[I]=0.5*(DXEP[I]+DXPW[I]);SNS[1]=0.0;SNS[NJ]=0.0;for(J=2;J<=NJM1;J++)SNS[J]=0.5*(DYNP[J]+DYPS[J]);for(I=2;I<=NIM1;I++){for(J=2;J<=NJM1;J++)T[I][J]=T0;}lable_1:NITER=NITER+1;for(I=2;I<=NIM1;I++){for(J=2;J<=NJM1;J++)LAM[I][J]=LAMDA0;}for(I=2;I<=NIM1;I++){for(J=2;J<=NJM1;J++){AE[I][J]=2*LAM[I][J]*LAM[I+1][J]/(LAM[I][J]+LAM[I+1][J])*SNS[J]/DXEP[I];AW[I][J]=2*LAM[I][J]*LAM[I-1][J]/(LAM[I][J]+LAM[I-1][J])*SNS[J]/DXPW[I];AS[I][J]=2*LAM[I][J]*LAM[I][J-1]/(LAM[I][J]+LAM[I][J-1])*SEW[I]/DYPS[J];AN[I][J]=2*LAM[I][J]*LAM[I][J+1]/(LAM[I][J]+LAM[I][J+1])*SEW[I]/DYNP[J];SU[I][J]=0;SP[I][J]=0;}}I=2;for(J=2;J<=NJM1;J++){AW[I][J]=ALFA*SNS[J];T[1][J]=TF;SU[I][J]=1.0E30*TF;SP[I][J]=-1.0E30;}I=NIM1;for(J=2;J<=NJM1;J++){SU[I][J]=1.0E30*T0;SP[I][J]=-1.0E30;}J=2;for(I=2;I<=NIM1;I++){AS[I][J]=0;T[I][1]=T[I][J];}J=NJM1;for(I=2;I<=NIM1;I++){AN[I][J]=ALFA*SEW[I];T[I][NJ]=TF;}for(I=2;I<=NIM1;I++){for(J=2;J<=NJM1;J++){AP[I][J]=AE[I][J]+AW[I][J]+AS[I][J]+AN[I][J]-SP[I][J]*SEW[I]*SNS[J];T2=AE[I][J]*T[I+1][J]+AW[I][J]*T[I-1][J]+AS[I][J]*T[I][J-1];T2=T2+AN[I][J]*T[I][J+1]+SU[I][J]*SEW[I]*SNS[J];T1[I][J]=T2/AP[I][J];}}NUM=0;for(I=2;I<=NIM1;I++){for(J=2;J<=NJM1;J++){IF(fabs(T1[I][J]-T[I][J])>EPS){NUM=NUM+1;T[I][J]=T1[I][J];}}}printf("%d,%d",NITER,NUM);if(NUM>0)gotolable_1;for(I=2;I<=NIM1;I++){for(J=2;J<=NJM1;J++){fprintf(fp1,"%6.1f",T1[I][J]);fprintf(fp2,"%9.4f,%9.4f,%9.4f",X[I],Y[J],T1[I][J]);}fprintf(fp1,"\n");}}

三、一维非稳态导热问题隐式格式程序设计1、目的和意义通过该上机练习，让学生更进一步掌握偏微分方程数值解的知识，特别是非稳定态导热问题隐式差分方法的解法，为将来进一步掌握二维和三维非稳定态导热问题隐式格式计算奠定基础。2、内部节点差分方程的建立；；；；；多步法该方法与稳态导热的形式完全相同，可以采用同样的方法求解。3、边界节点的差分方程的建立1）第三类边界条件（左边界）整理得：2）第二类边界条件（上边界），附加源项法3）第一类边界条件（右边界）3、C语言程序#include<stdio.h>#include<math.h>main(){intI;intNI=22;intNIM1=NI-1;floatX[22];floatSEW[22],DXEP[22],DXPW[22];floatAE[22],AW[22],AP[22],AP0[22];floatT0[22],T1[22],SU[22],SP[22],AP1[22],T10[22];floatLAM[22],DEN[22],CP[22];floatDX1=0.05,DX2=0.03,DX3=0.01;floatTF=1000,TI=100,LAMDA0=20,ALFA=200,B=0.0,EPS=.01,DENSIT=7000,CP0=15,DT=5;intNITER=0,TIME=0,TIMEEND=3600,NUM;FILE*fp1,*fp2;fp1=fopen("d:\y1","w");fp2=fopen("d:\y2","w");X[1]=-DX1/2;for(I=2;I<=10;I++)X[I]=X[I-1]+DX1;for(I=11;I<=18;I++)X[I]=X[I-1]+DX2;for(I=19;I<=NI;I++)X[I]=X[I-1]+DX3;fprintf(fp1,"%d,%d",NI);for(I=1;I<=NI;I++)fprintf(fp1,"%7.4f",X[I]);DXPW[1]=0.0;DXEP[NI]=0.0;for(I=1;I<=NIM1;I++){DXEP[I]=X[I+1]-X[I];DXPW[I+1]=DXEP[I];}SEW[1]=0.0;SEW[NI]=0.0;for(I=2;I<=NIM1;I++)SEW[I]=0.5*(DXEP[I]+DXPW[I]);for(I=2;I<=NIM1;I++)T0[I]=TI;}lable_1:TIME=TIME+DT;for(I=2;I<=NIM1;I++)LAM[I]=LAMDA0*(1+B*T0[I]);DEN[I]=DENSIT;CP[I]=CP0;}for(I=2;I<=NIM1;I++){AE[I]=2*LAM[I]*LAM[I+1]/(LAM[I]+LAM[I+1])/DXEP[I];AW[I]=2*LAM[I]*LAM[I-1]/(LAM[I]+LAM[I-1])/DXPW[I];AP0[I]=DEN[I]*CP[I]/DT*SEW[I];AP1[I]=DEN[I]*CP[I]/DT*SEW[I];SU[I]=0;SP[I]=0;}}I=2;{AW[I]=ALFA*SNS;T0[1]=TF;}I=NIM1;{AE[I]=ALFA*SNS;T0[NI]=TI;}for(I=2;I<=NIM1;I++){AP[I]=AP1[I]+AE[I]+AW[I]-SP[I]*SEW[I];SU[I]=SU[I]*SEW[I]+AP0[I]*T0[I];T10[I]=T0[I];}}lable_2:NITER=NITER+1;NUM=0;for(I=2;I<=NIM1;I++){T1[I]=(AE[I]*T10[I+1]+AW[I]*T10[I-1]+SU[I])/AP[I];if(fabs(T1[I]-T0[I])>EPS){NUM=NUM+1;T10[I]=T1[I];}}printf("%d,%d",NITER,NUM);if(NUM>0)gotolable_2;fprintf(fp1,"%f,%f,%f,%f,%f",TIME,T1[2],T1[NIM1]);if(TIME<TIMEEND)gotolable_1;for(I=2;I<=NIM1;I++)fprintf(fp2,"%6.1f",T1[I]);fprintf(fp2,"\n");}}

四、二维非稳态导热问题显式格式程序设计1、目的和意义使得学生初步掌握显式格式的求解过程。同时了解显示格式的稳定性问题。2、内部节点差分方程的建立 ；；；；；； 稳定性判断条件：3、边界节点的差分方程的建立1）第三类边界条件（左边界）第Ⅲ类边界条件；整理得：2）第二类边界条件（上边界），附加源项法3）第一类边界条件（右边界）4、C语言程序#include<stdio.h>#include<math.h>main(){intI,J;intNI=22,NJ=12;intNIM1=NI-1,NJM1=NJ-1;floatX[22],Y[12];floatSEW[22],SNS[22],DXEP[22],DXPW[22],DYNP[22],DYPS[22];floatAE[22][12],AW[22][12],AS[22][12],AN[22][12],AP[22][12],AP0[22][12];floatT0[22][12],T1[22][12],SU[22][12],SP[22][12],AP1[22][12];floatLAM[22][12],DEN[22][12],CP[22][12];floatDX1=0.05,DX2=0.03,DX3=0.01,DY1=0.01,DY2=0.05;floatTF=1000,TI=100,LAMDA0=20,ALFA=200,B=0.0,EPS=.01,DENSIT=7000,CP0=15,DT=5;intNITER=0,TIME=0,TIMEEND=3600,NUM;FILE*fp1,*fp2;fp1=fopen("d:\y1","w");fp2=fopen("d:\y2","w");X[1]=-DX1/2;for(I=2;I<=10;I++)X[I]=X[I-1]+DX1;for(I=11;I<=18;I++)X[I]=X[I-1]+DX2;for(I=19;I<=NI;I++)X[I]=X[I-1]+DX3;Y[1]=-DY1/2;for(J=2;J<=5;J++)Y[J]=Y[J-1]+DY1;for(J=6;J<=NJ;J++)Y[J]=Y[J-1]+DY2;fprintf(fp1,"%d,%d",NI,NJ);for(I=1;I<=NI;I++)fprintf(fp1,"%7.4f",X[I]);for(J=1;J<=NJ;J++)fprintf(fp1,"%7.4f",Y[J]);DXPW[1]=0.0;DXEP[NI]=0.0;for(I=1;I<=NIM1;I++){DXEP[I]=X[I+1]-X[I];DXPW[I+1]=DXEP[I];}DYPS[1]=0.0;DYNP[NJ]=0.0;for(J=1;J<=NJM1;J++){DYNP[J]=Y[J+1]-Y[J];DYPS[J+1]=DYNP[J];}SEW[1]=0.0;SEW[NI]=0.0;for(I=2;I<=NIM1;I++)SEW[I]=0.5*(DXEP[I]+DXPW[I]);SNS[1]=0.0;SNS[NJ]=0.0;for(J=2;J<=NJM1;J++)SNS[J]=0.5*(DYNP[J]+DYPS[J]);for(I=2;I<=NIM1;I++){for(J=2;J<=NJM1;J++)T0[I][J]=TI;}lable_1:TIME=TIME+DT;for(I=2;I<=NIM1;I++){for(J=2;J<=NJM1;J++)LAM[I][J]=LAMDA0*(1+B*T0[I][J]);DEN[I][J]=DENSIT;CP[I][J]=CP0; T0[I][J]=T1[I][J]}for(I=2;I<=NIM1;I++){for(J=2;J<=NJM1;J++){AE[I][J]=2*LAM[I][J]*LAM[I+1][J]/(LAM[I][J]+LAM[I+1][J])*SNS[J]/DXEP[I];AW[I][J]=2*LAM[I][J]*LAM[I-1][J]/(LAM[I][J]+LAM[I-1][J])*SNS[J]/DXPW[I];AS[I][J]=2*LAM[I][J]*LAM[I][J-1]/(LAM[I][J]+LAM[I][J-1])*SEW[I]/DYPS[J];AN[I][J]=2*LAM[I][J]*LAM[I][J+1]/(LAM[I][J]+LAM[I][J+1])*SEW[I]/DYNP[J];AP0[I][J]=DEN[I][J]*CP[I][J]/DT*SEW[I]*SNS[J];AP1[I][J]=DEN[I][J]*CP[I][J]/DT*SEW[I]*SNS[J];SU[I][J]=0;SP[I][J]=0;}}I=2;for(J=2;J<=NJM1;J++){AW[I][J]=ALFA*SNS[J];T0[1][J]=TF;}I=NIM1;for(J=2;J<=NJM1;J++){AE[I][J]=ALFA*SNS[J];T0[NI][J]=TI;}J=2;for(I=2;I<=NIM1;I++){AS[I][J]=0;T0[I][1]=T0[I][J];}J=NJM1;for(I=2;I<=NIM1;I++){AN[I][J]=ALFA*SEW[I];T0[I][NJ]=TF;}for(I=2;I<=NIM1;I++){for(J=2;J<=NJM1;J++){AP[I][J]=AP0[I][J]-(AE[I][J]+AW[I][J]+AS[I][J]+AN[I][J])+SP[I][J]*SEW[I]*SNS[J];SU[I][J]=SU[I][J]*SEW[I]*SNS[J];T1[I][J]=(AE[I][J]*T0[I+1][J]+AW[I][J]*T0[I-1][J]+AS[I][J]*T0[I][J-1]+AN[I][J]*T0[I][J+1]+AP[I][J]*T0[I][J]+SU[I][J])/AP1[I][J];}}printf("%d,%d",NITER,NUM);if(NUM>0)gotolable_2;fprintf(fp1,"%f,%f,%f,%f,%f",TIME,T1[2][2],T1[2][NJM1],T1[NIM1][2],T1[NIM1][NJM1]);if(TIME<TIMEEND)gotolable_1;for(I=2;I<=NIM1;I++){for(J=2;J<=NJM1;J++)fprintf(fp2,"%6.1f",T1[I][J]);fprintf(fp2,"\n");}}二、一维非稳态导热问题隐式格式程序设计1、内部节点差分方程的建立；；；；；多步法该方法与稳态导热的形式完全相同，可以采用同样的方法求解。2、边界节点的差分方程的建立1）第三类边界条件（左边界）整理得：2）第二类边界条件（上边界），附加源项法3）第一类边界条件（右边界）3、C语言程序#include<stdio.h>#include<math.h>main(){intI;intNI=22;intNIM1=NI-1;floatX[22];floatSEW[22],DXEP[22],DXPW[22];floatAE[22],AW[22],AP[22],AP0[22];floatT0[22],T1[22],SU[22],SP[22],AP1[22],T10[22];floatLAM[22],DEN[22],CP[22];floatDX1=0.05,DX2=0.03,DX3=0.01;floatTF=1000,TI=100,LAMDA0=20,ALFA=200,B=0.0,EPS=.01,DENSIT=7000,CP0=15,DT=5;intNITER=0,TIME=0,TIMEEND=3600,NUM;FILE*fp1,*fp2;fp1=fopen("d:\y1","w");fp2=fopen("d:\y2","w");X[1]=-DX1/2;for(I=2;I<=10;I++)X[I]=X[I-1]+DX1;for(I=11;I<=18;I++)X[I]=X[I-1]+DX2;for(I=19;I<=NI;I++)X[I]=X[I-1]+DX3;fprintf(fp1,"%d,%d",NI);for(I=1;I<=NI;I++)fprintf(fp1,"%7.4f",X[I]);DXPW[1]=0.0;DXEP[NI]=0.0;for(I=1;I<=NIM1;I++){DXEP[I]=X[I+1]-X[I];DXPW[I+1]=DXEP[I];}SEW[1]=0.0;SEW[NI]=0.0;for(I=2;I<=NIM1;I++)SEW[I]=0.5*(DXEP[I]+DXPW[I]);for(I=2;I<=NIM1;I++)T0[I]=TI;}lable_1:TIME=TIME+DT;for(I=2;I<=NIM1;I++)LAM[I]=LAMDA0*(1+B*T0[I]);DEN[I]=DENSIT;CP[I]=CP0;}for(I=2;I<=NIM1;I++){AE[I]=2*LAM[I]*LAM[I+1]/(LAM[I]+LAM[I+1])/DXEP[I];AW[I]=2*LAM[I]*LAM[I-1]/(LAM[I]+LAM[I-1])/DXPW[I];AP0[I]=DEN[I]*CP[I]/DT*SEW[I];AP1[I]=DEN[I]*CP[I]/DT*SEW[I];SU[I]=0;SP[I]=0;}}I=2;{AW[I]=ALFA*SNS;T0[1]=TF;}I=NIM1;{AE[I]=ALFA*SNS;T0[NI]=TI;}for(I=2;I<=NIM1;I++){AP[I]=AP1[I]+AE[I]+AW[I]-SP[I]*SEW[I];SU[I]=SU[I]*SEW[I]+AP0[I]*T0[I];T10[I]=T0[I];}}lable_2:NITER=NITER+1;NUM=0;for(I=2;I<=NIM1;I++){T1[I]=(AE[I]*T10[I+1]+AW[I]*T10[I-1]+SU[I])/AP[I];if(fabs(T1[I]-T0[I])>EPS){NUM=NUM+1;T10[I]=T1[I];}}printf("%d,%d",NITER,NUM);if(NUM>0)gotolable_2;fprintf(fp1,"%f,%f,%f,%f,%f",TIME,T1[2],T1[NIM1]);if(TIME<TIMEEND)gotolable_1;for(I=2;I<=NIM1;I++)fprintf(fp2,"%6.1f",T1[I]);fprintf(fp2,"\n");}}三、二维稳态导热问题程序设计1、内部节点差分方程的建立2、边界节点的差分方程的建立1）第三类边界条件（左边界）2）第二类边界条件（上边界）3）第一类边界条件（右边界）当Bi时，第三类边界条件转化为第一类边界条件。即3、C语言程序#include<stdio.h>#include<math.h>main(){intI,J;intNI=52,NJ=52;intNIM1=NI-1,NJM1=NJ-1;intTF=1000,T0=100,NITER=0,NUM;floatLAMD,A0=20,ALFA=200,B=0.01,EPS=0.01,T2;floatX[52],Y[52];floatSEW[52],SNS[52],DXEP[52],DXPW[52],DYNP[52],DYPS[52];floatAE[52][52],AW[52][52],AS[52][52],AN[52][52],AP[52][52];floatT[52][52],T1[52][52],SU[52][52],SP[52][52],LAM[52][52];floatDX1=0.05,DX2=0.03,DX3=0.01,DY1=0.01,DY2=0.05;FILE*fp1,*fp2;fp1=fopen("d:\y1","w");fp2=fopen("d:\y2","w");X[1]=-DX1/2;for(I=2;I<=10;I++)X[I]=X[I-1]+DX1;for(I=11;I<=18;I++)X[I]=X[I-1]+DX2;for(I=19;I<=NI;I++)X[I]=X[I-1]+DX3;Y[1]=-DY1/2;for(J=2;J<=5;J++)Y[J]=Y[J-1]+DY1;for(J=6;J<=NJ;J++)Y[J]=Y[J-1]+DY2;fprintf(fp1,"%d,%d",NI,NJ);for(I=1;I<=NI;I++)fprintf(fp1,"%7.4f",X[I]);for(J=1;J<=NJ;J++)fprintf(fp1,"%7.4f",Y[J]);DXPW[1]=0.0;DXEP[NI]=0.0;for(I=1;I<=NIM1;I++){DXEP[I]=X[I+1]-X[I];DXPW[I+1]=DXEP[I];}DYPS[1]=0.0;DYNP[NJ]=0.0;for(J=1;J<=NJM1;J++){DYNP[J]=Y[J+1]-Y[J];DYPS[J+1]=DYNP[J];}SEW[1]=0.0;SEW[NI]=0.0;for(I=2;I<=NIM1;I++)SEW[I]=0.5*(DXEP[I]+DXPW[I]);SNS[1]=0.0;SNS[NJ]=0.0;for(J=2;J<=NJM1;J++)SNS[J]=0.5*(DYNP[J]+DYPS[J]);for(I=2;I<=NIM1;I++){for(J=2;J<=NJM1;J++)T[I][J]=T0;}lable_1:NITER=NITER+1;for(I=2;I<=NIM1;I++){for(J=2;J<=NJM1;J++)LAM[I][J]=LAMDA0;}for(I=2;I<=NIM1;I++){for(J=2;J<=NJM1;J++){AE[I][J]=2*LAM[I][J]*LAM[I+1][J]/(LAM[I][J]+LAM[I+1][J])*SNS[J]/DXEP[I];AW[I][J]=2*LAM[I][J]*LAM[I-1][J]/(LAM[I][J]+LAM[I-1][J])*SNS[J]/DXPW[I];AS[I][J]=2*LAM[I][J]*LAM[I][J-1]/(LAM[I][J]+LAM[I][J-1])*SEW[I]/DYPS[J];AN[I][J]=2*LAM[I][J]*LAM[I][J+1]/(LAM[I][J]+LAM[I][J+1])*SEW[I]/DYNP[J];SU[I][J]=0;SP[I][J]=0;}}I=2;for(J=2;J<=NJM1;J++){AW[I][J]=ALFA*SNS[J];T[1][J]=TF;SU[I][J]=1.0E30*TF;SP[I][J]=-1.0E30;}I=NIM1;for(J=2;J<=NJM1;J++){SU[I][J]=1.0E30*T0;SP[I][J]=-1.0E30;}J=2;for(I=2;I<=NIM1;I++){AS[I][J]=0;T[I][1]=T[I][J];}J=NJM1;for(I=2;I<=NIM1;I++){AN[I][J]=ALFA*SEW[I];T[I][NJ]=TF;}for(I=2;I<=NIM1;I++){for(J=2;J<=NJM1;J++){AP[I][J]=AE[I][J]+AW[I][J]+AS[I][J]+AN[I][J]-SP[I][J]*SEW[I]*SNS[J];T2=AE[I][J]*T[I+1][J]+AW[I][J]*T[I-1][J]+AS[I][J]*T[I][J-1];T2=T2+AN[I][J]*T[I][J+1]+SU[I][J]*SEW[I]*SNS[J];T1[I][J]=T2/AP[I][J];}}NUM=0;for(I=2;I<=NIM1;I++){for(J=2;J<=NJM1;J++){IF(fabs(T1[I][J]-T[I][J])>EPS){NUM=NUM+1;T[I][J]=T1[I][J];}}}printf("%d,%d",NITER,NUM);if(NUM>0)gotolable_1;for(I=2;I<=NIM1;I++){for(J=2;J<=NJM1;J++){fprintf(fp1,"%6.1f",T1[I][J]);fprintf(fp2,"%9.4f,%9.4f,%9.4f",X[I],Y[J],T1[I][J]);}fprintf(fp1,"\n");}}四、二维非稳态导热问题显式格式程序设计1、内部节点差分方程的建立；；；；；；稳定性判断条件：2、边界节点的差分方程的建立1）第三类边界条件（左边界）第Ⅲ类边界条件；整理得：2）第二类边界条件（上边界），附加源项法3）第一类边界条件（右边界）3、C语言程序#include<stdio.h>#include<math.h>main(){intI,J;intNI=22,NJ=12;intNIM1=NI-1,NJM1=NJ-1;floatX[22],Y[12];floatSEW[22],SNS[22],DXEP[22],DXPW[22],DYNP[22],DYPS[22];floatAE[22][12],AW[22][12],AS[22][12],AN[22][12],AP[22][12],AP0[22][12];floatT0[22][12],T1[22][12],SU[22][12],SP[22][12],AP1[22][12];floatLAM[22][12],DEN[22][12],CP[22][12];floatDX1=0.05,DX2=0.03,DX3=0.01,DY1=0.01,DY2=0.05;floatTF=1000,TI=100,LAMDA0=20,ALFA=200,B=0.0,EPS=.01,DENSIT=7000,CP0=15,DT=5;intNITER=0,TIME=0,TIMEEND=3600,NUM;FILE*fp1,*fp2;fp1=fopen("d:\y1","w");fp2=fopen("d:\y2","w");X[1]=-DX1/2;for(I=2;I<=10;I++)X[I]=X[I-1]+DX1;for(I=11;I<=18;I++)X[I]=X[I-1]+DX2;for(I=19;I<=NI;I++)X[I]=X[I-1]+DX3;Y[1]=-DY1/2;for(J=2;J<=5;J++)Y[J]=Y[J-1]+DY1;for(J=6;J<=NJ;J++)Y[J]=Y[J-1]+DY2;fprintf(fp1,"%d,%d",NI,NJ);for(I=1;I<=NI;I++)fprintf(fp1,"%7.4f",X[I]);for(J=1;J<=NJ;J++)fprintf(fp1,"%7.4f",Y[J]);DXPW[1]=0.0;DXEP[NI]=0.0;for(I=1;I<=NIM1;I++){DXEP[I]=X[I+1]-X[I];DXPW[I+1]=DXEP[I];}DYPS[1]=0.0;DYNP[NJ]=0.0;for(J=1;J<=NJM1;J++){DYNP[J]=Y[J+1]-Y[J];DYPS[J+1]=DYNP[J];}SEW[1]=0.0;SEW[NI]=0.0;for(I=2;I<=NIM1;I++)SEW[I]=0.5*(DXEP[I]+DXPW[I]);SNS[1]=0.0;SNS[NJ]=0.0;for(J=2;J<=NJM1;J++)SNS[