一维非稳态导热通用程序
传热学传热学--第三章 第三节 一维非稳态导热问题
传热学--第三章第三节一维非稳态导热问题§3 — 3 一维非稳态导热的分析解本节介绍第三类边界条件下：无限大平板、无限长圆柱、球的分析解及应用。
如何理解无限大物体，如：当一块平板的长度、宽度>> 厚度时，平板的长度和宽度的边缘向四周的散热对平板内的温度分布影响很少，以至于可以把平板内各点的温度看作仅是厚度的函数时，该平板就是一块“无限大”平板。
若平板的长度、宽度、厚度相差较小，但平板四周绝热良好，则热量交换仅发生在平板两侧面，从传热的角度分析，可简化成一维导热问题。
一、无限大平板的分析解已知：厚度的无限大平板，初温t0，初始瞬间将其放于温度为的流体中，而且>t0，流体与板面间的表面传热系数为一常数。
试确定在非稳态过程中板内的温度分布。
解：如图3-5 所示，平板两面对称受热，所以其内温度分布以其中心截面为对称面。
对于x 0 的半块平板，其导热微分方程：(0<x< , )定解条件：t(x,0)= t0(0 x )（边界条件）（边界条件）引入过余温度：则（0<x< , ）(3-9)(x,0)= (0 x ) （初始条件）（边界条件）（边界条件）对偏微分方程分离变量求解得：（3-10 ）其中离散值是下列超越方程的根，称为特征值。
其中Bi 是以特征长度为的毕渥数。
由此可见：平板中的无量纲过余温度与三个无量纲数有关：以平板厚度一半为特征长度的傅立叶数、毕渥数及即：（3-12)二、非稳态导热的正规状况阶段1 、平板中任一点的过余温度与平板中心的过余温度的关系前述得到的分析解是一个无穷级数，计算工作量大，但对比计算表明，当Fo>0.2 时，采用该级数的第一项与采用完整的级数计算平板中心温度的误差小于1% ，因此，当Fo>0.2时，采用以下简化结果：（3-13 ）其中特征值之值与Bi 有关。
由上式（3-13 ）可知：Fo>0.2 以后平板中任一点的过余温度(x ，τ) 与平板中心的过余温度(0 ，τ)=（τ ）之比为：（3-14 ）此式反映了非稳态导热过程中一种很重要的物理现象：即当Fo>0.2 以后，虽然(x ，τ) 与（τ ）各自均与τ 有关，但其比值则与τ 无关，而仅取决于几何位置（x/ ）及边界条件（Bi ）。
传热学一维稳态和非稳态导热
11.1 通过平壁的一维稳态导热
一、第一类边界条件：表面温度
为常数
单层平壁
a. 几何条件：单层平板；s；
b. 物理条件：、cp、 已知；无内热源；
c. 时间条件：稳态导热， ∂T/∂t=0；
d. 边界条件：第一类。
微分方程式可简化：
2T x2
0
直接积分得： TC1xC2
带入边界 条件：
CC12
Tw1 (Tw2
流量，并记为qL
qL
Q Tw1 Tw2 L 1 d2 ln 传热学一维稳态和非稳态导热
2 d1
单位长度导热热阻
11.2 通过圆筒壁的一维稳态导热
多层圆筒壁
不同材料构成的多层圆筒壁，其导热 热流量可按总温差和总热阻计算
热流量
Q
Tw1 Twn+1 n 1 ln di1
i1 2i L di
单位长度的热流量
2
C1xC2
式中积分常数C1和C2可由边界条件确定，它们分别为：
C2Tw2qv s2; C10
所以，平壁内温度分布为： TTw2qv s2x2
• 可见，该条件下平壁内温度是按抛物线规律分布。令 温度分布关系式中的x=0，则得平壁中心温度为：
qv 2 T T s w
2 传热学一维稳态和非稳态导热
• 设在一管道外面包上一层绝热层(如图所示)。
• 此时单位管长的总热阻γΣ为：
d 1 11 2 1 传热1 学ln 一维d d 稳1 2 态 和非2 稳态1 导x 热lnd d 2 xd 1 x2
11.2 通过圆筒壁的一维稳态导热
• 式中：α1为管内流体与管内壁之间的给热系数，W/m2℃；α2为绝
qL
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)::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。
一维非稳态导热 圆柱体 matlab
一维非稳态导热圆柱体 matlab导热是物体中传热作用的一种。
热传导是指物质内部热量的传递及传递机制，描述的是能量（热量）在空间和时间上的传输。
而非稳态导热是指物体内部温度场和热流密度随时间和空间的发展过程。
一维非稳态导热问题是指导热物理学中，只考虑热量沿一个方向传导的问题。
圆柱体是一种常见的几何体，因其在工程领域的广泛应用，研究圆柱体的导热问题具有重要意义。
在研究一维非稳态导热圆柱体问题时，matlab是一种常用的数学软件，其强大的数学运算和可视化功能使得它成为了工程热传导问题求解的重要工具。
通过使用matlab，可以方便地求解一维非稳态导热圆柱体问题，并进行可视化展示。
下面我们将通过matlab来求解一维非稳态导热圆柱体问题，并对结果进行分析。
1. 问题建模假设圆柱体材料均匀，热传导系数为k，密度为ρ，比热容为c。
设圆柱体半径为R，长度为L。
假设圆柱体表面维持恒定的温度T0，初始时刻整个圆柱体的温度场分布为T(x,0) = f(x)，其中f(x)为已知函数。
根据热传导方程，我们可以得到一维非稳态导热圆柱体的数学模型。
2. 热传导方程根据一维热传导方程，我们可以得到圆柱体内部温度场满足的偏微分方程：ρc∂T/∂t = k∇²T3. 离散化为了利用计算机进行求解，我们需要将偏微分方程进行离散化处理。
这里我们可以使用有限差分法（finite difference method）对空间和时间进行离散化。
将圆柱体划分为若干个网格点，并采用显式差分法进行时间推进，就可以得到圆柱体温度场随时间的演化过程。
4. matlab求解在matlab中，我们可以编写程序来实现离散化求解。
首先可以定义圆柱体以及热传导材料的参数，然后通过循环计算每个时间步长内圆柱体温度场的演化，最终得到温度在空间和时间上的分布情况。
借助matlab强大的可视化功能，我们可以直观地展示圆柱体温度场的变化过程。
5. 结果分析得到圆柱体温度场的数值解之后，我们可以对结果进行分析。
数值传热学一维非稳态导热
数值传热学一维非稳态导热
数值传热学一维非稳态导热是一个拟表达热量输运多方面考虑下的相关分析技术，例如光斑热传递，带有间断层热传导，恒定物质热传导等等。
本文将重点简要介绍一维非稳态导热模型中的理论方法，为解决该问题提供重要基础。
首先，我们讨论的一维非稳态导热模型是一维的，在这种模型中，温度的变化是由上下相邻的单元格热传导加权平均值决定的，从一个单元格到另一个单元格的变化必须满足偏微分方程的通用表达式。
其次，根据以上的假设，一维非稳态导热的数值解将以定义的步长迭代，用于求解温度在不同单元中的变化。
在数值模拟中，需要对边界条件、热导率和温度输入进行有效描述，以确定最终的解答模式。
同时，本次分析中，利用有限差分和蒙特卡罗方法来求解温度场。
这种有趣且可行的做法，不但实现了所需求解的模式，而且能够精确地给出结果。
此外，在电脑指令中，采取该方法对数值运算很有效，从而提高了计算机解的精度和实现的质量。
最后，一维非稳态导热模型是在一定物理场中进行计算的，通用性很强，其能够很好地模拟简单模型中物理场的变化。
因此，它经常被用于诸如热管道传热、滑动轴热传导、负载温度场仿真等多种领域的研究。
总而言之，一维非稳态导热的数值模拟具有良好的数学基础、使用简单的算法以及电脑指令，从而实现快速求解热传导问题的目的，是今后研究的重要课题。
一维圆柱非稳态导热方程求解
一维圆柱非稳态导热方程求解引言在工程和科学领域中，导热方程是研究物体内部温度分布变化的重要方程。
本文将介绍一维圆柱非稳态导热方程的求解方法。
问题描述考虑一个半径为 R 的圆柱体，其内部由一个材料填充。
我们想要求解该圆柱体内部的温度分布，其中的热传导过程遵循非稳态导热方程。
在一维情况下，非稳态导热方程可以表示为：equation1equation1其中，T 是温度关于时间和半径的函数，α 是热扩散系数。
数值方法为了求解该方程，我们将应用有限差分法，离散化时间和半径的变量。
首先，我们将时间区域 [0, T] 离散化为 N 个等间距的时间步长，得到：equation2equation2其中，Δt 是时间步长，N 是总的时间步数。
接下来，我们将半径区域 [0, R] 离散化为 M 个等间距的半径步长，得到：equation3equation3其中，Δr 是半径步长，M 是总的半径步数。
将温度T 在时间和半径的离散点上进行近似，我们可以写出离散的导热方程为：equation4equation4其中，上标 n 表示时间步数，下标 i 表示半径步数。
由于我们已经将时间和半径离散化，上述方程可以重写为一个递推关系：equation5equation5数值实验我们将通过一个简单的数值实验来验证我们的数值方法。
设定圆柱体半径 R = 1，热扩散系数α = 1。
令边界条件为半径为 R 处的温度始终为 0，初始条件为半径为 R/2 处的温度为 1，其他位置的温度为 0。
在使用递推关系计算之前，需要选择合适的时间步长Δt 和半径步长Δr。
一般来说，选择小的步长可以提高数值解的准确性，但同时也会增加计算时间。
在这里，我们选择Δt = 0.001，Δr = 0.1。
我们使用一个循环来逐步计算每个时间步骤的温度分布，并在每个时间步骤结束后，使用画图软件绘制温度分布图。
结果与讨论经过数值计算和绘图，我们得到了随时间演化的温度分布图。
一维圆柱非稳态导热方程求解
一维圆柱非稳态导热方程求解（最新版）目录一、引言二、一维圆柱非稳态导热方程的建立三、求解方法四、数值计算与结果分析五、结论正文一、引言在工程热物理领域，非稳态导热问题是一个重要研究内容。
在实际应用中，许多热传导过程都是非稳态的，例如，物体在温度场作用下的热膨胀、热应变等过程。
对于这类问题，建立正确的非稳态导热方程并求解，具有重要的理论和实际意义。
本文针对一维圆柱非稳态导热问题，进行了深入研究，以期为解决实际工程问题提供理论依据。
二、一维圆柱非稳态导热方程的建立考虑一个一维圆柱体，其半径为 r，高为 h。
在非稳态热传导过程中，圆柱体内部各点的温度随时间发生变化。
根据热力学第一定律和热传导基本定律，可以建立一维圆柱非稳态导热方程。
假设圆柱体内部材质均匀，热导率恒定，则其非稳态导热方程可表示为：T(r,z,t)/t = α(T(r,z,t)/z)其中，T(r,z,t) 表示圆柱体内部某点的温度，α表示热导率，t 表示时间。
三、求解方法为了求解上述非稳态导热方程，采用有限差分法。
首先将圆柱体内部划分为若干个网格点，每个网格点的坐标为 (r_i, z_j)，其中 i 和 j 为网格点在圆柱体内的位置。
然后，对每个网格点，根据热传导方程，建立有限差分方程。
通过求解有限差分方程，可以得到圆柱体内部各点的温度场随时间的变化。
四、数值计算与结果分析为了验证所建立的非稳态导热方程的正确性，以及有限差分法的有效性，进行了数值计算。
首先，设定圆柱体的半径 r、高度 h、热导率α等参数。
然后，采用不同的时间步长和网格划分，进行数值计算。
计算结果表明，随着时间步长和网格划分的减小，计算结果趋于一致，说明所建立的非稳态导热方程和求解方法的正确性。
通过分析计算结果，可以得到以下结论：在非稳态导热过程中，圆柱体内部温度场呈现出明显的径向和轴向分布特征。
在径向方向，温度分布呈现出从内到外的逐渐降低趋势；在轴向方向，温度分布呈现出从上到下的逐渐降低趋势。
一维圆柱非稳态导热方程求解
一维圆柱非稳态导热方程求解摘要：1.问题背景及意义2.一维圆柱非稳态导热方程3.求解方法及步骤4.数值例子与分析5.应用及拓展正文：一、问题背景及意义在工程热传导领域，一维圆柱非稳态导热问题具有广泛的应用背景。
例如在金属加工、半导体制造、建筑节能等许多实际问题中，都需要研究热传导过程。
对于一维圆柱非稳态导热问题，了解其求解方法及过程有助于分析和优化实际工程中的热传导现象。
二、一维圆柱非稳态导热方程考虑一维圆柱非稳态导热问题，我们可以得到以下偏微分方程：$$frac{partial u}{partial t} = kfrac{partial^2 u}{partial r^2}$$其中，$u(r, t)$ 表示圆柱的热量分布，$k$ 为热传导系数，$r$ 为圆柱半径，$t$ 为时间。
三、求解方法及步骤1.边界条件处理：根据实际问题，给出圆柱的边界条件，如温度边界条件、热流密度边界条件等。
2.离散化：将空间坐标和时间坐标进行离散化，得到离散化的偏微分方程。
3.数值求解：采用有限差分法、有限元法等数值方法求解离散化方程。
4.分析结果：根据求解结果，分析圆柱热量分布随时间和空间的变化规律。
四、数值例子与分析以一个简单的一维圆柱非稳态导热问题为例，给出数值求解过程及结果分析。
1.设定参数：热传导系数$k = 10$，圆柱半径$r = 0.1$，时间步长$Delta t = 0.01$，空间步长$Delta r = 0.001$。
2.边界条件：圆柱表面受到均匀热流密度作用，即$q = 100$。
3.数值求解：采用有限差分法求解离散化方程。
4.结果分析：通过观察不同时刻的热量分布图，可以发现热量在圆柱内部传播的速度和规律。
根据分析，可以提出优化热传导过程的方法和措施。
五、应用及拓展一维圆柱非稳态导热方程求解方法不仅适用于金属加工、半导体制造等热传导领域，还可以拓展到其他领域，如地球物理学、生物医学等。
通过改变边界条件和相关参数，可以研究不同条件下一维圆柱非稳态导热现象的特点。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
!***************************************************************! 一维非稳态导热通用程序(不变部分)! 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。