1、11.题目如图 1 所示铜芯电缆，电流为 5000A，内径为 10mm，外包材料聚氯乙烯的厚度为2mm，导热系数为 0.15+0.00013t 。电缆左半边为绝热边界条件，右半边为K)W/(m第三类边界条件，空气温度为 20，绝缘层表面与环境间的复合表面传热系数为 10。铜的电阻率为 ，K)W/(m2 2010t-aRt， ，t 的单位为摄氏度。试通过数值方法求解温度分1075.80RC/4.a布。图 12.编程计算2.1 控制方程根据题意，本题为二维稳态导热问题，其控制方程为： 011STrrTr边界条件：23,07.r WfBThq：. 02其中： 。C20fT2.2 方程离散为建立通用方
2、程，考虑非稳态项的控制方程为： 011STrrTrtTc采用全隐格式，在 时间内，对控制容积积分，整理后可得：bTaTaTaSNWEP其中：， ，eeEra/wWra/， ，nN/sS/VaaPSNWEP 0， ，tVc0CTbrsn5.0采用通用表达式，各表达式如下表：表 1 坐标及系数表达式坐标系 极坐标 通用表达式东西坐标 X南北坐标 rY半径 R东西尺度系数 rS东西节点间距 X南北节点间距 rY3东西导热面积 rSXYR/南北导热面积 控制体体积 rEaWe/eXSYR/2NSnr/n/0PatYXRcbSCPa YXRSaaPNWEP02.3 边界条件处理对于北边界，采用附加源项法
3、处理。由于北边界( )为23,07.r第三类边界条件，则最靠近边界的控制容积加入以下附加源项： BnfadCxhTVAS/1, BnadP/,其中： C20fT将附加源项加到相应控制容积后，再令相应的 。0Na对于南边界，可认为定温边界条件，由于其导热面积为零， 。0Sa对于东西边界，计算时取 计算区域，故东西边界重合，可认为为定温边界条件，2温度为上一层相邻控制容积的温度。42.4 导热系数与计算取铜导热系数为常数， 。K)W/(m40每个控制容积各界面对应导热系数分别为 、 、 、 。对于铜芯或保温层内nsew部控制容积，各导热系数均为常数。两者交接界面的导热系数用调和平均法计算。2.5
4、方程求解方程采用 ADI-TDMA 方法求解，首先在 Y 方向进行隐式计算， X 方向采用显式计算。各方向对应方程为三对角矩阵，使用 TDMA 法求解。然后再在 X 方向进行隐式计算，Y 方向采用显式计算。53. 结果输出与分析3.1 计算结果程序中温度 T 为二维数组，采用坐标变换方法，将温度表示在极坐标系中。设定温度初场为 23，循环结束判定条件为 ，网格数为9E1/20T条件下，输出结果如图 2：7012图 23.2 网格独立性考察保持迭代精度 不变：8E1/20T1. 网格数为 时，计算结果为：426图 32. 网格数为 时，计算结果为：284图 473. 网格数为 时，计算结果为：1
5、402图 5结论：从以上各图可以看出，程序运行结果与网格划分无关，程序具有较好的网格独立性。3.3 收获与体会通过这次 matlab 编程作业，我对二维扩散问题有了更加深刻的理解，对网格划分、通用离散形式、边界条件处理等有了进一步的认识。在编写 Matlab 程序过程中，我为了直接求解三对角矩阵还曾编写一个 Solution.m 文件，经过对比后发现此文件相比于TDMA 方法在速度上稍微快一点，结果基本相同。通过编程，我更加深刻的认识到只有亲自动手才能加深对问题理解，才能真正获得属于自己的知识。84. 程序语句程序采用 Matlab 编写，主要分为 4 部分，分别是主程序，用于给定题目条件，调
6、用其他函数，循环求解等；网格划分函数 Grid.m，用于划分网格； SolutionTDMA.m，用于执行交替隐式计算；TDMA.m，用于求解三对角矩阵。4.1 Main.mclear allclear globalformat longglobal Xglobal Yglobal dXglobal dYglobal DXglobal DXnglobal DXsglobal DYglobal Cvglobal CVglobal Tglobal T0global Tfglobal nodXglobal nodYnodX=200;nodY=350;%给出题目参数X=2*pi;Y=7E-3;Grid
7、; %划分网格Tf=20;h=10;%计算导热系数for i=1:5/7*nodYLe(i)=400;%假设铜的导热系数为 400W/(m.K) Lw(i)=400;endfor i=5/7*nodY+1:nodYLe(i)=0.15;Lw(i)=0.15;endfor i=1:5/7*nodY-1Ln(i)=400;endfor i=5/7*nodY+1:nodYLn(i)=0.15;9endLn(5/7*nodY)=2/(1/Ln(5/7*nodY-1)+1/Ln(5/7*nodY+1);for i=1:5/7*nodYLs(i)=400;endfor i=5/7*nodY+2:nodYL
8、s(i)=0.15;endLs(5/7*nodY+1)=2/(1/Ls(5/7*nodY)+1/Ls(5/7*nodY+2);%设定初始值for i=1:nodXfor j=1:nodYT(i,j)=23;endendT0=T;%定义内热源for i=1:nodXfor j=1:5/7*nodYSp(i,j)=-28.37; Sc(i,j)=6525.08+56.74*T0(i,j);%用 T0 表示上时刻的值 endfor j=5/7*nodY+1:nodYSp(i,j)=0;Sc(i,j)=0;endend%计算系数aE=ones(nodX,1)*(DY.*Le./dX);aW=ones(
9、nodX,1)*(DY.*Lw./dX);aN=ones(nodX,1)*(DXn.*Ln/dY);aS=ones(nodX,1)*(DXs.*Ls/dY);aP0=0;%边界条件的处理，附加源项法Scad=DXn(nodY)/Cv(nodY)*Tf/(1/h+Y/2/nodY/Ln(nodY);Spad=-DXn(nodY)/Cv(nodY)/(1/h+Y/2/nodY/Ln(nodY);for i=1:nodXaN(i,nodY)=0;aS(i,1)=0;endfor i=1:nodX/410Sc(i,nodY)=Sc(i,nodY)+Scad;Sp(i,nodY)=Sp(i,nodY)+
10、Spad;endfor i=3*nodX/4+1:nodXSc(i,nodY)=Sc(i,nodY)+Scad;Sp(i,nodY)=Sp(i,nodY)+Spad;endb=Sc.*CV;aP=aE+aW+aN+aS-aP0-Sp.*CV;SolutionTDMA(aE,aW,aN,aS,aP,b);cont=1norm(T-T0)./T)while (norm(T-T0)./T)1E-8)cont=cont+1norm(T-T0)./T)T0=T;for i=1:nodXfor j=1:5/7*nodYSc(i,j)=6525.08+56.74*T0(i,j); %用 T0表示上时刻的值
11、endfor j=5/7*nodY+1:nodYSc(i,j)=0;endendfor i=1:nodX/4Sc(i,nodY)=Sc(i,nodY)+Scad;endfor i=3*nodX/4+1:nodXSc(i,nodY)=Sc(i,nodY)+Scad;endb=Sc.*CV;aP=aE+aW+aN+aS-aP0-Sp.*CV;SolutionTDMA(aE,aW,aN,aS,aP,b);end%输出图形theta=X/nodX;dR=Y/nodY;for i=1:nodX+1for j=1:nodYX(i,j)=j*dR*cos(theta*(i-1);Y(i,j)=j*dR*sin(theta*(i-1);X(i,1)=0;Y(i,1)=0;endendfor i=1:nodX+1for j=1:nodYT(nodX+1,j)=T(1,j);Z(i,j)=T(i,j);end
