下载文档原格式
东南大学能源与环境学院课程作业报告课程名称：传热学作业名称：传热学大作业——利用matlab程序解决热传导问题院（系）：能源与环境学院专业：热能与动力工程姓名：姜学号：完成时间：2012 年11 月8日评定成绩：审阅教师：目录一．题目及要求 (3)二．各节点离散化的代数方程..............................3&13 三．源程序......................................................5&16 四．不同初值时的温度分布情况...........................7&18 五．冷量损失的计算.......................................12&24 六．计算小结 (27)传热大作业——利用matlab 程序解决复杂热传导问题姓名:姜 学号: 班级:成绩：____________________一、题目及要求计算要求：一个长方形截面的冷空气通道的尺寸如附图所示。
假设在垂直于纸面的方向上冷空气及通道墙壁的温度变化很小，可以忽略。
试用数值方法计算下列两种情况下通道壁面中的温度分布及每米长度上通过壁面的冷量损失：（1） 内、外壁面分别维持在10℃及30℃；（2） 内、外壁面与流体发生对流传热，且有110f t C =︒、2120/()h W m K =⋅，230f t C =︒、224/()h W m K =⋅。
(取管道导热系数为0.025/()W m K λ=⋅)二、各节点的离散化的代数方程1、基本思想：将导热问题的温度场，用有限个离散点上的值的集合来代替，通过求解按一定方法建立起来的关于这些值的代数方程，来获得离散点上被求物理量的值。
2、基本步骤：(1)建立控制方程以及定解条件：对于(1)问有：2.2m3m 2m1.2m h 1、t f1h 1、t f2导热微分方程22220t t x y ∂∂+=∂∂定解条件为第一类边界条件对(2)问有： 导热微分方程22220t t x y ∂∂+=∂∂定解条件为第三类边界条件(2)区域离散化：如下图所示，用一系列与坐标轴平行的网格线把求解区域划分成许多子区域，以网格线的交点作为需要确定温度值的空间位置，称为节点。
通过在室内的某些位置布置适当的节点，采集回来室内的温湿度以及空气质量等实际参数。
首先对室内空间建模，用一个无限细化的三维矩阵来模拟出室内的温度分布情况，针对采集回来的数据，采用插值法和适当次数的拟合函数的拟合，得出三维矩阵的实际值的分布，最后结合matlab软件绘制出计算出的温度场的三维图像。
一．数据的采集与处理因为影响人的舒适感的温度层只是室内的某一高度范围内的温度，而温度传感器虽然是布置在一个平面内，但是采用插值法和拟合函数法是可以大致再现出影响人的舒适感的温度层的温度变化的。
同时，在构建出的三维模型中，用第三维表示传感器层面的温度。
在传感器层面，传感器分布矩阵如下：X=【7.5 36.5 65.5】（模型内单位为cm）Y=【5.5 32.5 59.5】Z=【z1 z2 z3；z4 z5 z6；z7 z8 z9；】（传感器采集到的实时参数）采用meshgrid（xi，yi，zi，…）产生网格矩阵；首先按照人的最小温度分辨值，将室内的分布矩阵按照同样的比例细化，均分，使取值点在坐标一定程度上也是接近于连续变化的，从而才能最大程度上使处理数据得来的分布值按最小分辨值连续变化！根据人体散热量计算公式：C=hc（tb-Ta）其中hc为对流交换系数；结合Gagge教授提出的TSENS热感觉指标可以计算出不同环境下人的对环境温度变化时人体温度感知分辨率，作为插值法的一个参考量，能使绘制出的温度场更加的符合人体的温度变化模式。
例如按照10cm的均差产生网格矩阵（实际上人对温度的分辨率是远远10cm大于这个值的，但是那样产生的网格矩阵也是异常庞大的，例如以0.5cm为例，那么就可以获得116*108=12528个元素，为方便说明现已10cm为例）：[xi yi]=meshgrid(7.5:10:65.5,5.5:10:59.5)xi =7.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.5000yi =5.5000 5.5000 5.5000 5.5000 5.5000 5.500015.5000 15.5000 15.5000 15.5000 15.5000 15.500025.5000 25.5000 25.5000 25.5000 25.5000 25.500035.5000 35.5000 35.5000 35.5000 35.5000 35.500045.5000 45.5000 45.5000 45.5000 45.5000 45.500055.5000 55.5000 55.5000 55.5000 55.5000 55.5000产生网格矩阵之后，就可以在测得的实时数据的基础上，通过相关的温度场的专业的估算函数，以及相关的数值处理函数来估计整个分布面（有最小的分辨率）上的温度了。
大作业报告【9-16】题目略 分析过程：● 第一步，建立节点温度差分方程。
有题目可以看出，本题有7类不同的节点，需要7个不同的方程。
由于y x ∆=∆在原始方程中将它们消去后有： 1、底边的中间5个点满足的方程：x q t t t t w j i j i j i j i ∆+++=-+-1,,1,1,222λλλλ2、右边对流换热的中间3个点满足的方程：∞∆+--∆+++=+t B t t t t B i j i j i j i j i i 22)42(1,1,,1,3、左边绝热部分的中间3个点满足的方程：1,1,,1,222+-+++=j i j i j i j i t t t t λλλλ4、上边恒温点满足的方程：w j i t t =.5、中间5×3的网格里的点满足的方程：1,1,,1,1,4+-+-+++=j i j i j i j i j i t t t t t6、左下角的既属于绝热层又接受加热的点：0)(2)(22,1,,,1=-+-+∆⋅++j i j i j i j i w t t t t x q λλ 7、右下角的既被加热又与流体有对流换热的点：0)(2)(2)(22,,1,,,1=-⋅∆⋅+-+-+∆⋅+-j i f j i j i j i j i w t t y h t t t t x q λλ● 第二步，采用高斯-塞德尔迭代法进行迭代：建立矩阵，左下点为原点(1,1)，i 值向右增加，j 值向上增加。
按照上面7种类型的分组依次进行计算，每次计算直接用计算结果代替原值，完成一轮迭代。
程序设定200次迭代。
根据最终结果判断，100次迭代就已经结果收敛。
● 结果分析： Matlab 运算结果：表1 节点的温度分布（i为横行，j为竖列）5 400 400 400 400 400 400 4004 397.1890 396.8375 395.7156 393.5919 389.9486 383.6080 371.52583 394.2443 394.4454 392.4332 388.7034 382.5945 372.9575 357.94582 394.2443 393.4299 390.8685 386.1939 378.7685 367.6818 351.89201 395.0363 394.1616 391.4168 386.4354 378.6037 367.1090 351.1356 节点/˚C 123456 7做出三维分布图像：结果分析：由数值计算结果，可以看出：1、最顶部的点温度恒定，满足题设条件。
数学物理方法论文（圆柱体齐次边界条件以及matlab的可视化）学院: 信息院年级: 2011级班级: 通信一班姓名: ***一、应用背景:现在由于公寓楼采用集体供暖的措施,需要在地下铺设圆柱形长管道,而对于管道的半径选取需要进行一定的计算,才能使得热水在传送过程中不会因温度过低而达不到用户的要求。
对于以上的问题，我们可以简化到应用数学中来解决。
二、简化例题：若一供暖公司采用的运输管道为标准圆柱体，其半径为ρ。
，长为L ，设管入口有均匀分布的强度为q 。
的热流流入，出口有相同的热流流出，管道侧面保持温度为0℃。
求解管道内的稳恒温度。
三、例题解答：解： 因为上下底非齐次边界的非齐次项是常数，故可较易化成齐次边界。
这样本征值问题就变成傅里叶级数本征问题，而不是贝塞尔函数本征问题，同时系数的求解也较为简单。
令 z 01kq =νv kq u +=0则v 的定解问题为：⎪⎪⎪⎩⎪⎪⎪⎨⎧=-===∆===z k q v L z z z z 00200,00ρρννν 分离变数得到的本征值问题为：⎪⎩⎪⎨⎧===+==0',0'0"02L z z Z Z Z h Z 解得 ()......2,1,0,==n Ln h n πz Ln A Z n n πcos =，问题在柱内的有限解为：()⎪⎪⎭⎫⎝⎛∞=⋅=∑ρππρνL n n n I z L n A z 00cos ， 由初始条件，得z kq IA n z L n L n 00ncos -=∑∞=⋅⎪⎪⎭⎫⎝⎛ππρ将上式右端展为傅里叶余弦级数，则有()dz L zn z k q L n LI n LA ππρcos 2000⋅-=⎰L L z n z n L L z n n L Ln L k I q 02000si n c o s 2-⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛⋅+=πππππρ()()[]112-2000--⎪⎭⎫ ⎝⎛=nn L L n L k Iq ππρ ()()()⎪⎩⎪⎨⎧+===12,402000220m n L n I n k Lq m n πρπ但不等于，()()()00,2-01000000=-==⎰I kL q zdz k q LI A L ()∑∞=+⋅⎪⎭⎫ ⎝⎛++⎪⎭⎫⎝⎛++=∴0002020012cos 12121242-m z L m L m I m L m I k L q k L q ππρπρπνν+=∴kqu四、以上问题的Matlab 的可视化：建立圆柱体：t=linspace(-pi,pi,200); y=linspace(0,4,200); [T,Y]=meshgrid(t,y); X=sin(T); Z=cos(T); mesh(X,Y,Z); axis equal由于本题中需要四维的数学建模（三维立体图形，再加一维的温度坐标），对于专业知识要求较高，以现有的知识水平还无法对四维进行分析求解，因此采用三维坐标（二维管道的横截面，再加一维的平面温度分布），进行求解。
数值计算大作业一、用数值方法求解尺度为100mm×100mm 的二维矩形物体的稳态导热问题。
物体的导热系数λ为1.0w/m·K。
边界条件分别为： 1、上壁恒热流q=1000w/m2; 2、下壁温度t1=100℃； 3、右侧壁温度t2=0℃; 4、左侧壁与流体对流换热，流体温度tf=0℃，表面传热系数 h 分别为1w/m2·K、10 w/m2·K、100w/m2·K 和1000 w/m2·K;要求：1、写出问题的数学描述；2、写出内部节点和边界节点的差分方程；3、给出求解方法；4、编写计算程序(自选程序语言)；5、画出4个工况下的温度分布图及左、右、下三个边界的热流密度分布图；6、就一个工况下（自选）对不同网格数下的计算结果进行讨论；7、就一个工况下（自选）分别采用高斯迭代、高斯——赛德尔迭代及松弛法（亚松弛和超松弛）求解的收敛性（cpu 时间，迭代次数）进行讨论；8、对4个不同表面传热系数的计算结果进行分析和讨论。
9、自选一种商业软件（fluent 、ansys 等）对问题进行分析，并与自己编程计算结果进行比较验证（一个工况）。
（自选项）1、写出问题的数学描述 设H=0.1m微分方程 22220t tx y∂∂+=∂∂x=0，0<y<H ：()f th t t xλ∂-=-∂ 定解条件 x=H ，0<y<H ：t=t 2 y=0，0<x<H ：t=t1t 1t 2h ；t fq=1000 w/m 2y=H ，0<x<H ：tq yλ∂-=∂ 2、写出内部节点和边界节点的差分方程 内部节点：()()1,,1,,1,,122220m n m n m nm n m n m n t t t t t t x y -+-+-+-++=∆∆左边界： (),1,,1,1,,,022m n m n m n m nm n m n f m n t t t t t t x x h y t t y y y xλλλ-++---∆∆∆-+++∆=∆∆∆右边界： t m,n =t 2上边界： 1,,1,,,1,022m n m n m n m nm n m n t t t t t t y y q x x x x yλλλ-+----∆∆∆+++∆=∆∆∆ 下边界： t m,n =t 13、求解过程利用matlab 编写程序进行求解，先在matlab 中列出各物理量，然后列出内部节点和边界节点的差分方程，用高斯-赛德尔迭代法计算之后用matlab 画图。
传热学数值计算大作业传热学数值计算大作业一选题《传热学》第四版P179页例题 4-3二相关数据及计算方法1．厚2δ=0.06m的无限大平板受对称冷却，故按一半厚度作为模型进行计算2. δ=0.03m，初始温度t0=100℃,流体温度t∞=0℃；λ=40W/(m.K),h=1000W/(m2.K),Bi=h*△x/λ=0.25；3.设定Fo=0.25和Fo=1两种情况通过C语言编程（源程序文件见附件）进行数值分析计算；当Fo=0.25时，Fo<1/(2*(1+Bi)),理论上出现正确的计算结果；当Fo=1时，Fo>1/(2*(1+Bi))，Fo>0.5,理论上温度分布出现振荡，与实际情况不符。
三网格划分将无限大平面的一半划分为6个控制体，共7个节点。
△x=0.03/N=0.03/6=0.005，即空间步长为0.005m四节点离散方程绝热边界节点即i=1时，tij+1=2Fo△ti+1j+(1-2Fo△)tij 内部节点即0tij+1=tij（1-2Fo△Bo△-2Fo△）+2Fo△ti-1j+2Fo△Bo△tf五温度分布线图（origin）六结果分析1 空间步长，时间步长对温度分布的影响空间步长和时间步长决定了Bo和Fo，两者越小计算结果越精确，但同时计算所需的时间就越长。
2 Fo数的大小对计算结果的影响编程时对Fo=1及0.25的情况分别进行了计算，发现当Fo=1时，各点温度随时间发生振荡，某点的温度高反而会使下一时刻的温度变低，违反了热力学第二定律，因此在计算中对Fo的选取有限制。
为了保证各项前的系数均为正值，对于内节点，Fo>0.5；对于对流边界节点，Fo<1/(2*(1+Bi))。
3 备注在Fo=0.25时，为了反映较长时间后温度的分布，取T=600,并选取了其中部分时刻的温度输出进行画图。
图像显示，随着时间的增长，各点温度趋向一致。
而当Fo=1时由于结果会出现振荡，只取T=6观察即可。
传热学二维稳态导热问题的数值解法杨达文2011151419赵树明2011151427杨文晓2011151421吴鸿毅2011151416第一题：a=linspace(0,0.6,121);t1=[60+20*sin(pi*a/0.6)];t2=repmat(60,[80 121]);s=[t1;t2]; %构造矩阵for k=1:10000000 %理论最大迭代次数，想多大就设置多大S=s;for j=2:120for i=2:80S(i,j)=0.25*(S(i-1,j)+S(i+1,j)+S(i,j-1)+S(i,j+1));endendif norm(S-s)<0.0001break; %如果符合精度要求，提前结束迭代elses=S;endendS %输出数值解数值解数据量太大，这里就不打印出来，只画出温度分布。
画出温度分布：figure(1)xx=linspace(0,0.6,121);yy=linspace(0.4,0,81);[x,y]=meshgrid(xx,yy);surf(x,y,S)axis([0 0.6 0 0.4 60 80])grid onxlabel('L1')ylabel('L2')zlabel('t（温度）').60.66666777778L 1L 2t （温度）A0=[S(:,61)];for k=1:81B1(k)=A0(81-k+1);endB1 %x=L1/2时y方向的温度A1=[S(41,:)] %y=L2/2时x方向的温度x=0:0.005:0.6;y=0:0.005:0.4;A2=60+20*sin(pi*x/0.6)*((exp(pi*0.2/0.6)-exp(-pi*0.2/0.6))/2)/((exp(pi*0.4/0.6)-exp(-pi*0.4/0.6) )/2) %计算y=L2/2时x方向的解析温度B2=60+20*sin(pi*0.3/0.6)*((exp(pi*y/0.6)-exp(-pi*y/0.6))/2)/((exp(pi*0.4/0.6)-exp(-pi*0.4/0.6))/ 2) %计算x=L1/2时y方向的解析温度figure(2)subplot(2,2,1);plot(x,A1,'g-.',x,A2,'k:x'); %画出x=L1/2时y方向的温度场、画出x=L1/2时y方向的解析温度场曲线xlabel('L1');ylabel('t温度');title('y=L2/2');legend('数值解','解析解');subplot(2,2,2);plot(x,A1-A2); %画出具体温度场与解析温度场的差值曲线xlabel('L1');ylabel('差值');title('y=L2/2时，比较=数值解-解析解');subplot(2,2,3);plot(y,B1,'g-.',y,B2,'k:x'); %画出y=L2/2时x方向的温度场、画出y=L2/2时x方向的解析温度场曲线xlabel('L2');ylabel('t温度');title('x=L1/2');legend('数值解','解析解');subplot(2,2,4);plot(y,B1-B2); %画出具体温度场与解析温度场的差值曲线xlabel('L2');ylabel('差值');title('x=L1/2时，比较=数值解-解析解');y=L2/2时x方向的温度：60 60.1635347276130 60.3269574318083 60.4901561107239 60.653018915996160.8154342294146 60.9772907394204 61.1384775173935 61.298884093677961.4584005332920 61.6169175112734 61.7743263876045 61.930519281669662.0853891461909 62.2388298405943 62.3907362037523 62.541004126057762.6895306207746 62.8362138946214 62.9809534175351 63.123649991570263.2642058188844 63.4025245687647 63.5385114436490 63.672073244095163.8031184326565 63.9315571966177 64.0573015095482 64.180265191631864.3003639687311 64.4175155301449 64.5316395850212 64.642657917384664.7504944397430 64.8550752452343 64.9563286582797 65.054185283707565.1485780543131 65.2394422768254 65.3267156762441 65.410338438521565.4902532515567 65.5664053444751 65.6387425251668 65.707215216057165.7717764880854 65.8323820928694 65.8889904930310 65.941562890665265.9900632539310 66.0344583417471 66.0747177265744 66.110813815270166.1427218680003 66.1704200151959 66.1938892725421 66.213113553990066.2280796827826 66.2387774004857 66.2451993740203 66.247341200688866.2452014111934 66.2387814706441 66.2280857775556 66.213121660833566.1938993747528 66.1704320919304 66.1427358942990 66.110829762085766.0747355608048 66.0344780262737 65.9900847476605 65.941586148577365.8890154662295 65.8324087286383 65.7718047299493 65.707245003846265.6387737950858 65.5664380291767 65.4902872802189 65.410373736929465.3267521668755 65.2394798789402 65.1486166840471 65.054224854168964.9563690796505 64.8551164248743 64.7505362822981 64.642700324897664.5316824570463 64.4175587638655 64.3004074590802 64.180308831415964.0573451895733 63.9316008058186 63.8031618582281 63.672116371626463.5385541572596 63.4025667512431 63.2642473518283 63.123690755529062.9809932921539 62.8362527587866 62.6895683527611 62.541040603677462.3907713045038 62.2388634418130 62.0854211252013 61.930549515936761.7743547548873 61.6169438897778 61.4584248018242 61.298906131798361.1384972055701 60.9773079591820 60.8154488635041 60.653030848523060.4901652273162 60.3269636197632 60.1635378760476 60x=L1/2时y方向的温度：60 60.1308958471008 60.2618814819943 60.3930468323419 60.524481948785060.6562770664196 60.7885226663977 60.9213095376979 61.054728839108661.1888721614654 61.3238315901874 61.4596997681540 61.596569958966661.7345361106384 61.8736929197574 62.0141358961654 62.155961428198162.2992668485325 62.4441505006859 62.5907118062120 62.739051332642462.8892708622179 63.0414734614594 63.1957635516239 63.352246980097063.5110310927684 63.6722248074423 63.8359386883315 64.002285021688564.1713778926236 64.3433332631650 64.5182690516120 64.696305213238964.8775638224022 65.0621691561100 65.2502477791090 65.441928630549065.6373431122839 65.8366251788694 66.0399114293203 66.247341200688866.4590566635297 66.6752029193167 66.8959280998773 67.121383468913967.3517235256817 67.5871061108928 67.8276925149213 68.073647588380968.3251398551535 68.5823416279436 68.8454291264398 69.114582598162569.3899864420822 69.6718293350911 69.9603043614169 70.255609145064670.5579459853794 70.8675219958221 71.1845492460516 71.509244907413471.8418314019312 72.1825365549057 72.5315937512233 72.889242095483173.2557265760494 73.6312982331452 74.0162143310978 74.410738534857774.8151410909089 75.2296990126956 75.6546962706925 76.090423987246276.5371806363247 76.9952722483076 77.4650126199600 77.946723529732178.4407349585321 78.9473853161230 79.4670216732992 8066666666L 1t 温度y =L 2/2--1.--0.-3L 1差值y =L 2/2时，比较=数值解-解析解66778L 2t 温度x =L 1/200.050.10.150.20.250.30.350.4--1.--0.-3L 2差值x =L 1/2时，比较=数值解-解析解。
传热学数值计算大作业航14 艾迪2011011537 如图所示，有一个正方形截面的无限长的水泥柱，热导率为，密度为，比热容为。
水泥柱的边长为。
水泥柱的左侧靠墙，可以认为保持温度为。
水泥柱被包围在温度为°的热空气中。
三个面上均只考虑对流换热，并且对流换热系数分别为，，。
请编写程序数值求解该稳态导热问题（可使用Fortran 或C 或Matlab 语言）。
作业要求提交源代码和报告，报告内容包括：（1） 给出该导热问题的数学描述； （2） 描述所采用的差分格式和求解过程；（3） 验证求解结果的准确性，给出网格无关性验证； （4） 给出求解结果（温度云图、边界热流、平均温度等）； （5） （选做）讨论对流换热系数、热导率等参数对求解结果的影响。
解：(1)、因为无内热源，温度分布：222201230(0,0)(x,0)t(0,y)t ,((x,0))(,y)(x,)((,y)),((x,H))f f f t tx H y H x ydt h t t dx dt H dt H h t H t h t t dx dxλλλ∂∂+=<<<<∂∂⎧=-=-⎪⎪⎨⎪-=--=-⎪⎩（2）、采用热平衡法建立内节点和边界节点的离散方程，x 、y 方向各取n 个节点，即()()11n n -⨯- 个网格，且x y ∆=∆ 。
对于任意内节点（i ，j ），有:,1,1,,1,1t (t t t t )/4i j i j i j i j i j -+-+=+++D边界三边界一边界节点：边界1、 1,0(1j )j t t n =≤≤边界2、11,1,21,11,1h 2h 2(2)t 2t t t (1k n)k k k k f xxt λλ-+∆∆+=+++<<边界3、22,k n 1,k n,k 1,k 1h 2h 2(2)t 2t t t (1k n)n n f xxt λλ--+∆∆+=+++<<边界4、33k,n k,n 11,n k 1,n h 2h 2(2)t 2t t t (1k n)k fxxt λλ--+∆∆+=+++<<C 点、2121n,1n 1,1n,2(h h )(h )(2)t t t f xh xt λλ-+∆+∆+=++D 点、2323n,n ,n 11,n (h h )(h )(2)t t t n n f xh xt λλ--+∆+∆+=++(3)、由于各个节点都写成了差分显示表达，可用高斯—赛德尔迭代法求解。
东南大学能源与环境学院课程作业报告作业名称：传热学大作业——利用matlab程序解决热传导问题院系：能源与环境学院专业：建筑环境与设备工程学号：姓名：2014年11月9日一、题目及要求1.原始题目及要求2.各节点的离散化的代数方程3.源程序4.不同初值时的收敛快慢5.上下边界的热流量（λ=1W/(m℃））6.计算结果的等温线图7.计算小结题目：已知条件如下图所示：二、各节点的离散化的代数方程各温度节点的代数方程ta=(300+b+e)/4 ; tb=(200+a+c+f)/4; tc=(200+b+d+g)/4; td=(2*c+200+h)/4 te=(100+a+f+i)/4; tf=(b+e+g+j)/4; tg=(c+f+h+k)/4 ; th=(2*g+d+l)/4ti=(100+e+m+j)/4; tj=(f+i+k+n)/4; tk=(g+j+l+o)/4; tl=(2*k+h+q)/4tm=(2*i+300+n)/24; tn=(2*j+m+p+200)/24; to=(2*k+p+n+200)/24; tp=(l+o+100)/12 三、源程序【G-S迭代程序】【方法一】函数文件为：function [y,n]=gauseidel(A,b,x0,eps)D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);G=(D-L)\U;f=(D-L)\b;y=G*x0+f;n=1;while norm(y-x0)>=epsx0=y;y=G*x0+f;n=n+1;end命令文件为：A=[4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0,0;-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0;0,-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0;0,0,-2,4,0,0,0,-1,0,0,0,0,0,0,0,0;-1,0,0,0,4,-1,0,0,-1,0,0,0,0,0,0,0;0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0,0;0,0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0;0,0,0,-1,0,0,-2,4,0,0,0,-1,0,0,0,0;0,0,0,0,-1,0,-1,0,4,0,0,0,-1,0,0,0;0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0,0;0,0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0;0,0,0,0,0,0,0,-1,0,0,-2,4,0,0,0,-1;0,0,0,0,0,0,0,0,-2,0,0,0,24,-1,0,0;0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1,0;0,0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1;0,0,0,0,0,0,0,0,0,0,0,-1,0,0,-1,12];b=[300,200,200,200,100,0,0,0,100,0,0,0,300,200,200,100]';[x,n]=gauseidel(A,b,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]',1.0e-6) xx=1:1:4;yy=xx;[X,Y]=meshgrid(xx,yy);Z=reshape(x,4,4);Z=Z'contour(X,Y,Z,30)Z =139.6088 150.3312 153.0517 153.5639108.1040 108.6641 108.3119 108.1523 84.1429 67.9096 63.3793 62.4214 20.1557 15.4521 14.8744 14.7746 【方法2】>> t=zeros(5,5);t(1,1)=100;t(1,2)=100;t(1,3)=100;t(1,4)=100;t(1,5)=100;t(2,1)=200;t(3,1)=200;t(4,1)=200;t(5,1)=200;for i=1:10t(2,2)=(300+t(3,2)+t(2,3))/4 ;t(3,2)=(200+t(2,2)+t(4,2)+t(3,3))/4;t(4,2)=(200+t(3,2)+t(5,2)+t(4,3))/4;t(5,2)=(2*t(4,2)+200+t(5,3))/4;t(2,3)=(100+t(2,2)+t(3,3)+t(2,4))/4;t(3,3)=(t(3,2)+t(2,3)+t(4,3)+t(3,4))/4; t(4,3)=(t(4,2)+t(3,3)+t(5,3)+t(4,4))/4; t(5,3)=(2*t(4,3)+t(5,2)+t(5,4))/4;t(2,4)=(100+t(2,3)+t(2,5)+t(3,4))/4;t(3,4)=(t(3,3)+t(2,4)+t(4,4)+t(3,5))/4;t(4,4)=(t(4,3)+t(4,5)+t(3,4)+t(5,4))/4;t(5,4)=(2*t(4,4)+t(5,3)+t(5,5))/4;t(2,5)=(2*t(2,4)+300+t(3,5))/24;t(3,5)=(2*t(3,4)+t(2,5)+t(4,5)+200)/24;t(4,5)=(2*t(4,4)+t(3,5)+t(5,5)+200)/24;t(5,5)=(t(5,4)+t(4,5)+100)/12;t'endcontour(t',50);ans =100.0000 200.0000 200.0000 200.0000 200.0000 100.0000 136.8905 146.9674 149.8587 150.7444 100.0000 102.3012 103.2880 103.8632 104.3496 100.0000 70.6264 61.9465 59.8018 59.6008 100.0000 19.0033 14.8903 14.5393 14.5117【Jacobi迭代程序】函数文件为：function [y,n]=jacobi(A,b,x0,eps)D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);B=D\(L+U);f=D\b;y=B*x0+f;n=1;while norm(y-x0)>=epsx0=y;y=B*x0+f;n=n+1;end命令文件为：A=[4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0,0;-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0; 0,-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0; 0,0,-2,4,0,0,0,-1,0,0,0,0,0,0,0,0;-1,0,0,0,4,-1,0,0,-1,0,0,0,0,0,0,0; 0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0,0; 0,0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0;0,0,0,-1,0,0,-2,4,0,0,0,-1,0,0,0,0;0,0,0,0,-1,0,-1,0,4,0,0,0,-1,0,0,0;0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0,0;0,0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0;0,0,0,0,0,0,0,-1,0,0,-2,4,0,0,0,-1;0,0,0,0,0,0,0,0,-2,0,0,0,24,-1,0,0;0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1,0;0,0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1;0,0,0,0,0,0,0,0,0,0,0,-1,0,0,-1,12];b=[300,200,200,200,100,0,0,0,100,0,0,0,300,200,200,100]'; [x,n]=jacobi(A,b,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]',1.0e-6); xx=1:1:4;yy=xx;[X,Y]=meshgrid(xx,yy);Z=reshape(x,4,4);Z=Z'contour(X,Y,Z,30)n =97Z =139.6088 150.3312 153.0517 153.5639108.1040 108.6641 108.3119 108.152384.1429 67.9096 63.3793 62.421420.1557 15.4521 14.8744 14.7746四、不同初值时的收敛快慢1、[方法1]在Gauss 迭代和Jacobi 迭代中，本程序应用的收敛条件均为norm(y-x0)>=eps ，即使前后所求误差达到e 的-6次方时，跳出循环得出结果。
将误差改为0.01时，只需迭代25次，如下[x,n]=gauseidel(A,b,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]',0.01)运行结果为 将误差改为0.1时，需迭代20次，可见随着迭代次数增加，误差减小，变化速度也在减小。
[方法2]通过 i=1:10判断收敛，为迭代10次，若改为1:20，则迭代20次。
2、在同样的误差要求下，误差控制在e 的-6次方内，Gauss 迭代用了49次达到要求，而Jacobi 迭代用了97次，可见，在迭代中尽量采用最新值，可以大幅度的减少迭代次数，迭代过程收敛快一些。
在Gauss 中，初值为100，迭代46次达到精确度1.0e-6,初值为50时，迭代47次，初值为0时，迭代49次，初值为200时迭代50次，可见存在一个最佳初始值，是迭代最快。
这一点在jacobi 迭代中表现的尤为明显。
五、上下边界的热流量：上边界t=200℃,∞t =10℃，所以，热流量Φ1=λ*[2*100-200x y ∆∆+x y a∆∆t -200+x y ∆∆b t -200+x y ∆∆c t -200+2*t -200d x y ∆∆] =1*（100/2+(200-139.6088)+(200-150.3312)+(200-153.0517)+(200-153.5639)/2) =230.2264W 下边界热流量Φ2=|λ*[x y ∆∆mi t -t +x y∆∆o j t -t +x y ∆∆p k t -t +2*t -t q l x y ∆∆]- h*(2*10-100x y ∆∆+x *t -t n ∆∆∞y +x *t -t o ∆∆∞y +x *t-t m ∆∆∞y +2*t -t p x y ∆∆∞)|=|1*((84.1429-20.1557)+(67.9096-15.4521)+(63.3793-14.8744)+(62.4214- 14.7746)/2)-10*(90/2+(20.1557-10)+(15.4521-10)+(14.8744-10)+(14.7746-10)/2)| = |-489.925|W =489.25W六、温度等值线Gauss:Yacobi:七、计算小结导热问题进行有限差分数值计算的基本思想是把在时间、空间上连续的温度场用有限个离散点温度的集合来代替，即有限点代替无限点，通过求解根据傅里叶定律和能量守恒两大法则建立关于控制面内这些节点温度值的代数方程，获得各个离散点上的温度值。
