资源描述：
《一维非稳态导热通用程序(好)》由会员上传分享，免费在线阅读，更多相关内容在工程资料-天天文库。
1、!一维非稳态导热通用程序(不变部分)!Thisisageneralpurposeprogramtosolve1・Ddiffusion!problemintheformof:pcdt/dz=l/a(x)d/dx(a(x)rdt/dx)+s!DefineVariablesMODULEVARIABLESINTEGERTARAMETER::L1=130REALQIMENSION(Ll)::X,XF,XM,XP,R,RF,APREAL,DIMENSION(L1)::AE,AW,CN,T,TAREAL,DIM
2、ENS1ON(L1)::TG,GM,RCINTEGER::K=1,KM=1,KP=1,OM=1INTEGER::JB,JE,KE,KI,KF,KNKR,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!MainProgramPROGRAMMAINUSEVARIABLESIMPLICITNONEINTEGERIOPEN(l,FILE=”q.da
3、t”)OPEN(2,FILE=”temp.dat”)NF=1!(求解变量指标)KN=1!(非稳态问题输出局部变量,输出一次加1)TU=0!(当前时间)150KT=1!(非线性问题迭代次数)CALLSpeciXFirsttospecifytheproblemCALLGrid!Setupgridpoints200CALLDifsor!Specifythediff-coeffandsourceterm220CALLInterOutput!OutputintermediateresultsCALLCocf
4、T!SetupcoefficientsofdiscretizationequationCALLTDMA!SolvethealgebraicequationbyTDMAIF(LS.EQ.2.OR.LS.EQ.4)THEN!(对非线性问题)IF(DF.GT.EP)THEN!(如果最大偏差大于允许值)DO1=1,MlTA(1)=TA(1)+OM*(T(1)-TA(I))!(采用亚松弛方式将当前解付给上一次迭代值)ENDDODF=0!(最大偏差置零)KT=KT+1!（非线性问迭代次数加1）GOTO200!
5、（转到DIFSOR模块重新计算扩散系数与源项）ENDIFENDIFCALLGPRINT!（四类问题均要经过一般输出）IF（LS.EQ.3.OR.LS.EQ.4）THEN!（对非稳态问题）IF（TU.LT.TM）THEN!（时间小于设定的最大值）DO1=1,MlTG（I）=T（I）!（当前计算结果付给上一时层）ENDDOKT=1IF（LSJEQ.3）THEN!（非稳态线性问题）GOTO220!（转到中间输出模块）ELSE（非稳态非线性问题）GOTO200!（转到DIFSOR模块重新计算扩散系数与源项
6、）ENDIFENDIFENDIF!specialresultsprintout,ifnot,justleaveitopenCALLSPRINTIF(NF.NE.KM)THENNF=NF+1GOTO150ENDIFCLOSE⑵CLOSE(l)END!SubroutineSUBROUTINESETUPUSEVARIABLESREAL,D1MENS1ON(L1)::P,QENTRYCOEFF!coefficientsofboundarypointsIF(KI.LE.1)THENAP(1)=1AE(l)=
7、0AW(l)=0CN(1)=AIELSEAE(1)=GM(1)/XM(2)AP(1)=AE⑴+BIAW(1)=OCN(1)=AIENDIFIF(KE.LE.1)THENAP(M1)=1AE(M1)=OAW(M1)=OCN(M1)=AECELSEAW(M1)=GM(M1)/XP(M2)AP(M1)=AW(M1)+BEAE(M1)=OCN(M1)=AECENDIF!coefficientsofinternalpointsIF(LS.NE.3.OR.TU.LT.0.5*DT)THENEX=1IF(MD.
8、EQ.3)EX=2AW(2)=GM(2)/XM(2)*RF(2)**EXAE(M2)=GM(M2)/XP(M2)*RF(M1)**EXENDIFDO1=2,M2-1AE(I)=RF(I+1)**EX/(XP(I)/GM(I)+XM(I+1)/GM(I+1))AW(I+1)=A玖I)ENDDODO1=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
