1、!*! 一维非稳态导热通用程序 (不变部分)! 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
2、):TG,GM,RCINTEGER: K=1,KM=1,KP=1,OM=1INTEGER: JB,JE,KE,KI,KF,KN,KR,KT,LS,MD,M1,M2,NF REAL: 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 VARIABLESIMPLICIT NONEINTEGER IOPEN(1,FILE=“q.dat“)OPEN(2,FILE=“temp.dat“)NF=1！（求解变量指标）KN=1！（非稳态问题
3、输出局部变量，输出一次加 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 al
4、gebraic 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！（对非稳态问题）I
5、F(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 openCALL SPRINTIF(NF.NE.KM) THENNF=NF+1GOTO 150END IFCLOSE(2)CLOS
6、E(1)END !.Subroutine.SUBROUTINE SETUPUSE VARIABLESREAL,DIMENSION(L1):P,QENTRY COEFF!coefficients of boundary pointsIF(KI.LE.1) THENAP(1)=1AE(1)=0AW(1)=0CN(1)=AIELSEAE(1)=GM(1)/XM(2)AP(1)=AE(1)+BIAW(1)=0CN(1)=AIEND IFIF(KE.LE.1) THENAP(M1)=1AE(M1)=0AW(M1)=0CN(M1)=AEC ELSEAW(M1)=GM(M1)/XP(M2)AP(M1)=AW
7、(M1)+BEAE(M1)=0CN(M1)=AEC END IF!.coefficients of internal pointsIF(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
8、)*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 IF END 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)/
9、(AP(I)-AW(I)*P(I-1)END DO!back substitutionT(M1)=Q(M1)DO I=M2,1,-1 T(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+DTRETUR
10、NEND IFEND IF!.surface flux calculationSELECT 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(
11、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(*,*) “Nonunifor
12、m 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
13、=“,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 DO END IFQI=QI*RF(2)*EXQE=QE*RF(M1)*EXIF(KF.EQ.1) THENWRITE(*,*) “Total Heat Flow At Int.Surface Qi=“,QIWRITE(*,*) “Total Heat Flow At Ext.Surface Qe=“,QEIF(LS.EQ.1.OR.LS.EQ.2) THENWRITE(*,*) “Total heat Input Form Source Term S=“,SWRITE(*,*) “Heat Balance:Qi+Qe+ys=“,QI+QE+SEND IFEND IFTU=TU+DTRETURNEND
