!***************************************************************!一维非稳态导热通用程序(不变部分)!Thisisageneralpurposeprogramtosolve1-Ddiffusion!problemintheformof:!ρcdt/dz=1/a(x)d/dx(a(x)Γdt/dx)+s!******************2003.7revised********************************!...................DefineVariables.........................MODULEVARIABLESINTEGER,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!...............................................MainProgram...........................................PROGRAMMAINUSEVARIABLESIMPLICITNONEINTEGERIOPEN(1,FILE="q.dat")OPEN(2,FILE="temp.dat")NF=1!(求解变量指标)KN=1!(非稳态问题输出局部变量,输出一次加1)TU=0!(当前时间)150KT=1!(非线性问题迭代次数)CALLSpeci!FirsttospecifytheproblemCALLGrid!Setupgridpoints200CALLDifsor!Specifythediff-coeffandsourceterm220CALLInterOutput!OutputintermediateresultsCALLCoeff!SetupcoefficientsofdiscretizationequationCALLTDMA!SolvethealgebraicequationbyTDMAIF(LS.EQ.2.OR.LS.EQ.4)THEN!(对非线性问题)IF(DF.GT.EP)THEN!(如果最大偏差大于允许值)DOI=1,M1TA(I)=TA(I)+OM*(T(I)-TA(I))!(采用亚松弛方式将当前解付给上一次迭代值)ENDDODF=0!(最大偏差置零)KT=KT+1!(非线性问题迭代次数加1)GOTO200!(转到DIFSOR模块重新计算扩散系数与源项)ENDIFENDIFCALLGPRINT!(四类问题均要经过一般输出)IF(LS.EQ.3.OR.LS.EQ.4)THEN!(对非稳态问题)IF(TU.LT.TM)THEN!(时间小于设定的最大值)DOI=1,M1TG(I)=T(I)!(当前计算结果付给上一时层)ENDDOKT=1IF(LS.EQ.3)THEN!(非稳态线性问题)GOTO220!(转到中间输出模块)ELSE(非稳态非线性问题)GOTO200!(转到DIFSOR模块重新计算扩散系数与源项)ENDIFENDIFENDIF!.............specialresultsprintout,ifnot,justleaveitopen..........CALLSPRINTIF(NF.NE.KM)THENNF=NF+1GOTO150ENDIFCLOSE(2)CLOSE(1)END!...................Subroutine.................SUBROUTINESETUPUSEVARIABLESREAL,DIMENSION(L1)::P,QENTRYCOEFF!................coefficientsofboundarypoints..........IF(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)=AIENDIFIF(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)=AECENDIF!...................coefficientsofinternalpoints........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)**EXENDIFDOI=2,M2-1AE(I)=RF(I+1)**EX/(XP(I)/GM(I)+XM(I+1)/GM(I+1))AW(I+1)=AE(I)ENDDODOI=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)ENDIFENDDORETURNENTRYTDMA!....................elimination.............P(1)=AE(1)/AP(1)Q(1)=CN(1)/AP(1)DOI=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))ENDDO!..................backsubstitution..........T(M1)=Q(M1)DOI=M2,1,-1T(I)=P(I)*T(I+1)+Q(I)ENDDOIF(LS.EQ.2.OR.LS.EQ.4)THENDOI=1,M1DS=ABS(T(I)-TA(I))IF(T(I).GT.1.E-20)DS=DS/T(I)IF(DF.LT.DS)DF=DSENDDOENDIFRETURNENTRYGPRINTIF(LS.EQ.3.OR.LS.EQ.4)THENM=(TU+0.5*DT)/(K*DT)IF(M.NE.KN)THENTU=TU+DTRETURNENDIFENDIF!...................surfacefluxcalculation..........SELECTCASE(KI)CASE(1)QI=GM(1)*(T(1)-T(2))/XM(2)CASE(2)QI=AICASE(3)QI=AI-BI*T(1)ENDSELECTSELECTCASE(KE)CASE(1)QE=GM(M1)*(T(M1)-T(M2))/XP(M2)CASE(2)QE=AECCASE(3)QE=AEC-BE*T(M1)ENDSELECTIF(LS.EQ.1.OR.LS.EQ.2)THENS=0DOI=2,M2QE=QE*RF(M1)**EXAP(I)=AE(I)+AW(I)-AP(I)S=S+CN(I)+AP(I)*T(I)ENDDOENDIF!................nowitisreadytoprintout.........KN=KN+1IF(KP.NE.2)THENWRITE(*,*)"DependentVariablesDistribution"SELECTCASE(MD)CASE(1)WRITE(*,*)"CartisianCoordinates"CASE(2)WRITE(*,*)"CylindricalCoordinates"CASE(3)WRITE(*,*)"SphericalCoordinates"CASE(4)WRITE(*,*)"Nonuniformcrosssection"ENDSELECTSELECTCASE(LS)CASE(1)WRITE(*,*)"LinearSteadyProblem"CASE(2)WRITE(*,*)"NonlinearSteadyProblem"WRITE(*,*)"Iterative=",KTCASE(3)WRITE(*,*)"LinearUnsteadyProblem"CASE(4)WRITE(*,*)"NonlinearUnsteadyProblem"WRITE(*,*)"Iterative=",KTWRITE(*,*)"AtTime=",TENDSELECTJE=0DOWHILE(JE.LT.M1)JB=JE+1JE=JE+4IF(JE.GT.M1)JE=M1WRITE(*,*)"J"DOJ=JB,JEWRITE(*,*)JENDDOIF(MD.EQ.2.OR.MD.EQ.3)THENWRITE(*,*)"R"ELSEWRITE(*,*)"X"ENDIFDOJ=JB,JEWRITE(*,*)T(J)ENDDOENDDOENDIFQI=QI*RF(2)**EXQE=QE*RF(M1)**EXIF(KF.EQ.1)THENWRITE(*,*)"TotalHeatFlowAtInt.SurfaceQi=",QIWRITE(*,*)"TotalHeatFlowAtExt.SurfaceQe=",QEIF(LS.EQ.1.OR.LS.EQ.2)THENWRITE(*,*)"TotalheatInputFormSourceTermS=",SWRITE(*,*)"HeatBalance:Qi+Qe+ys=",QI+QE+SENDIFENDIFTU=TU+DTRETURNEND