!***************************************************************! 一维非稳态导热通用程序(不变部分)! This is a general purpose program to solve 1-D diffusion! problem in the form of:! ρcdt/dz=1/a(x)d/dx(a(x)Γdt/dx)+s!******************2003.7 revised******************************** !...................Define Variables.........................MODULE VARIABLESINTEGER,PARAMETER::L1=130REAL,DIMENSION(L1):: X,XF,XM,XP,R,RF,APREAL,DIMENSION(L1)::AE,AW,CN,T,TAREAL,DIMENSION(L1)::TG,GM,RCINTEGER:: K=1,KM=1,KP=1,OM=1INTEGER:: JB,JE,KE,KI,KF,KN,KR,KT,LS,MD,M1,M2,NFREAL:: AEC,AI,BE,BI,DF,DS,DT,EP,EX,PW,TU,TM,XE,XIREAL:: A1,A2,T2,TC,SC,SP,RO,TE,DN,LMEND!...............................................Main Program...........................................PROGRAM MAINUSE V ARIABLESIMPLICIT NONEINTEGER IOPEN(1,FILE="q.dat")OPEN(2,FILE="temp.dat")NF=1！（求解变量指标）KN=1！（非稳态问题输出局部变量，输出一次加1）TU=0！（当前时间）150 KT=1！（非线性问题迭代次数）CALL Speci !First to specify the problemCALL Grid!Set up grid points200 CALL Difsor!Specify the diff-coeff and source term220 CALL InterOutput!Output intermediate resultsCALL Coeff!Set up coefficients of discretization equationCALL TDMA!Solve the algebraic equation by TDMAIF(LS.EQ.2.OR.LS.EQ.4) THEN ！（对非线性问题）IF(DF.GT.EP) THEN！（如果最大偏差大于允许值）DO I=1,M1TA(I)=TA(I)+OM*(T(I)-TA(I))！（采用亚松弛方式将当前解付给上一次迭代值）END DODF=0！（最大偏差置零）KT=KT+1！（非线性问题迭代次数加1）GOTO 200！（转到DIFSOR模块重新计算扩散系数与源项）END IFEND IFCALL GPRINT !(四类问题均要经过一般输出)IF(LS.EQ.3.OR.LS.EQ.4) THEN！（对非稳态问题）IF(TU.LT.TM) THEN！（时间小于设定的最大值）DO I=1,M1TG(I)=T(I)！（当前计算结果付给上一时层）END DOKT=1IF(LS.EQ.3) THEN！（非稳态线性问题）GOTO 220！（转到中间输出模块）ELSE（非稳态非线性问题）GOTO 200！（转到DIFSOR模块重新计算扩散系数与源项）END IFEND IFEND IF!.............special results print out, if not, just leave it open..........CALL SPRINTIF(NF.NE.KM) THENNF=NF+1GOTO 150END IFCLOSE(2)CLOSE(1)END!...................Subroutine.................SUBROUTINE SETUPUSE V ARIABLESREAL,DIMENSION(L1)::P,QENTRY COEFF!................coefficients of boundary points..........IF(KI.LE.1) THENAP(1)=1AE(1)=0AW(1)=0CN(1)=AIELSEAE(1)=GM(1)/XM(2)AP(1)=AE(1)+BICN(1)=AIEND IFIF(KE.LE.1) THENAP(M1)=1AE(M1)=0AW(M1)=0CN(M1)=AECELSEAW(M1)=GM(M1)/XP(M2)AP(M1)=AW(M1)+BEAE(M1)=0CN(M1)=AECEND IF!...................coefficients of internal points........IF(LS.NE.3.OR.TU.LT.0.5*DT) THENEX=1IF(MD.EQ.3) EX=2AW(2)=GM(2)/XM(2)*RF(2)**EXAE(M2)=GM(M2)/XP(M2)*RF(M1)**EXEND IFDO I=2,M2-1AE(I)=RF(I+1)**EX/(XP(I)/GM(I)+XM(I+1)/GM(I+1))AW(I+1)=AE(I)END DODO I=2,M2AP(I)=AE(I)+AW(I)-AP(I)*(XF(I+1)-XF(I))*R(I)**EXCN(I)=CN(I)*(XF(I+1)-XF(I))*R(I)**EXIF(LS.EQ.3.OR.LS.EQ.4) THENAP(I)=AP(I)+RC(I)*(XF(I+1)-XF(I))*R(I)**EX/DTCN(I)=CN(I)+RC(I)*(XF(I+1)-XF(I))*R(I)**EX/DT*TG(I) END IFEND DORETURNENTRY TDMA !....................elimination.............P(1)=AE(1)/AP(1)Q(1)=CN(1)/AP(1)DO I=2,M1P(I)=AE(I)/(AP(I)-AW(I)*P(I-1))Q(I)=(CN(I)+AW(I)*Q(I-1))/(AP(I)-AW(I)*P(I-1))END DO!..................back substitution..........DO I=M2,1,-1T(I)=P(I)*T(I+1)+Q(I)END DOIF(LS.EQ.2.OR.LS.EQ.4) THENDO I=1,M1DS=ABS(T(I)-TA(I))IF(T(I).GT.1.E-20) DS=DS/T(I)IF(DF.LT.DS) DF=DSEND DOEND IFRETURNENTRY GPRINTIF(LS.EQ.3.OR.LS.EQ.4) THENM=(TU+0.5*DT)/(K*DT)IF(M.NE.KN) THENTU=TU+DTRETURNEND IFEND IF !...................surface flux calculation.......... SELECT CASE(KI)CASE(1)QI=GM(1)*(T(1)-T(2))/XM(2)CASE(2)QI=AICASE(3)QI=AI-BI*T(1)END SELECTSELECT CASE(KE)CASE(1)QE=GM(M1)*(T(M1)-T(M2))/XP(M2) CASE(2)QE=AECCASE(3)QE=AEC-BE*T(M1)END SELECTIF(LS.EQ.1.OR.LS.EQ.2) THENS=0DO I=2,M2QE=QE*RF(M1)**EXAP(I)=AE(I)+AW(I)-AP(I)S=S+CN(I)+AP(I)*T(I)END DOEND IF !................now it is ready to print out.........KN=KN+1IF(KP.NE.2) THENWRITE(*,*) "Dependent Variables Distribution" SELECT CASE(MD)CASE(1)WRITE(*,*) "Cartisian Coordinates"CASE(2)WRITE(*,*) "Cylindrical Coordinates"CASE(3)WRITE(*,*) "Spherical Coordinates"CASE(4)WRITE(*,*) "Nonuniform cross section" END SELECTSELECT CASE(LS)CASE(1)WRITE(*,*) "Linear Steady Problem"CASE(2)WRITE(*,*) "Nonlinear Steady Problem"WRITE(*,*) "Iterative Times=", KTCASE(3)WRITE(*,*) "Linear Unsteady Problem"CASE(4)WRITE(*,*) "Nonlinear Unsteady Problem"WRITE(*,*) "Iterative Times=", KTWRITE(*,*) "At Time=",TEND SELECTJE=0DO WHILE(JE.LT.M1)JB=JE+1JE=JE+4IF(JE.GT.M1) JE=M1WRITE(*,*) "J"DO J=JB,JEWRITE(*,*) JEND DOIF(MD.EQ.2.OR.MD.EQ.3) THENWRITE(*,*) "R"ELSEWRITE(*,*) "X"END IFDO J=JB,JEWRITE(*,*) T(J)END DOEND DOEND IFQI=QI*RF(2)**EXQE=QE*RF(M1)**EXIF(KF.EQ.1) THENWRITE(*,*) "Total Heat Flow At Int.Surface Qi=",QI WRITE(*,*) "Total Heat Flow At Ext.Surface Qe=",QE IF(LS.EQ.1.OR.LS.EQ.2) THENWRITE(*,*) "Total heat Input Form Source Term S=",S WRITE(*,*) "Heat Balance:Qi+Qe+ys=",QI+QE+S END IFEND IFTU=TU+DTRETURNEND。