铜芯电缆温度分布MTLAB计算模型

11.题目如图1所示铜芯电缆，电流为5000A，内径为10mm，外包材料聚氯乙烯的厚度为2mm，导热系数为0.15+0.00013{t}K)W/(m。电缆左半边为绝热边界条件，右半边为第三类边界条件，空气温度为20℃，绝缘层表面与环境间的复合表面传热系数为10K)W/(m2。铜的电阻率为2010t-aRRt，m1075.180R，C/004.0a，t的单位为摄氏度。试通过数值方法求解温度分布。图12.编程计算2.1控制方程根据题意，本题为二维稳态导热问题，其控制方程为：011STrrrTrrr边界条件：2~23,2~0007.0r：WfBTThq23~2007.0r：0Bq其中：C20fT。22.2方程离散为建立通用方程，考虑非稳态项的控制方程为：011STrrrTrrrtTc采用全隐格式，在t时间内，对控制容积积分，整理后可得：bTaTaTaTaTaSSNNWWEEPP其中：eeeErra/，/，nnnNrra/，sssSrra/，VSaaaaaaPPSNWEP0tVcaP0，00PPCTaVSb，rrrVsn5.0采用通用表达式，各表达式如下表：表1坐标及系数表达式坐标系极坐标通用表达式东西坐标X南北坐标rY半径rR东西尺度系数rSX东西节点间距rSXX南北节点间距rY东西导热面积rSXYR/南北导热面积rXR3控制体体积rrYXREaWaeerr/eeXSXYR/2NaSannrr/nnYXR/0PatYXRc/bYXRSCPaYXRSaaaaaaPPSNWEP02.3边界条件处理对于北边界，采用附加源项法处理。由于北边界(2~23,2~0007.0r)为第三类边界条件，则最靠近边界的控制容积加入以下附加源项：BnfadCxhTVAS//1,BnadPxhVAS//11,其中：C20fT将附加源项加到相应控制容积后，再令相应的0Na。对于南边界，可认为定温边界条件，由于其导热面积为零，0Sa。对于东西边界，计算时取2计算区域，故东西边界重合，可认为为定温边界条件，温度为上一层相邻控制容积的温度。2.4导热系数与计算取铜导热系数为常数，K)W/(m400。每个控制容积各界面对应导热系数分别为n、s、e、w。对于铜芯或保温层内部控制容积，各导热系数均为常数。两者交接界面的导热系数用调和平均法计算。42.5方程求解方程采用ADI-TDMA方法求解，首先在Y方向进行隐式计算，X方向采用显式计算。各方向对应方程为三对角矩阵，使用TDMA法求解。然后再在X方向进行隐式计算，Y方向采用显式计算。53.结果输出与分析3.1计算结果程序中温度T为二维数组，采用坐标变换方法，将温度表示在极坐标系中。设定温度初场为23℃，循环结束判定条件为9E1/20TTT，网格数为70120条件下，输出结果如图2：图23.2网格独立性考察保持迭代精度8E1/20TTT不变：1.网格数为1424时，计算结果为：6图32.网格数为2848时，计算结果为：图473.网格数为140200时，计算结果为：图5结论：从以上各图可以看出，程序运行结果与网格划分无关，程序具有较好的网格独立性。3.3收获与体会通过这次matlab编程作业，我对二维扩散问题有了更加深刻的理解，对网格划分、通用离散形式、边界条件处理等有了进一步的认识。在编写Matlab程序过程中，我为了直接求解三对角矩阵还曾编写一个Solution.m文件，经过对比后发现此文件相比于TDMA方法在速度上稍微快一点，结果基本相同。通过编程，我更加深刻的认识到只有亲自动手才能加深对问题理解，才能真正获得属于自己的知识。84.程序语句程序采用Matlab编写，主要分为4部分，分别是主程序，用于给定题目条件，调用其他函数，循环求解等；网格划分函数Grid.m，用于划分网格；SolutionTDMA.m，用于执行交替隐式计算；TDMA.m，用于求解三对角矩阵。4.1Main.mclearallclearglobalformatlongglobalXglobalYglobaldXglobaldYglobalDXglobalDXnglobalDXsglobalDYglobalCvglobalCVglobalTglobalT0globalTfglobalnodXglobalnodYnodX=200;nodY=350;%给出题目参数X=2*pi;Y=7E-3;Grid;%划分网格Tf=20;h=10;%计算导热系数fori=1:5/7*nodYLe(i)=400;%假设铜的导热系数为400W/(m.K)Lw(i)=400;endfori=5/7*nodY+1:nodYLe(i)=0.15;Lw(i)=0.15;endfori=1:5/7*nodY-1Ln(i)=400;endfori=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));fori=1:5/7*nodYLs(i)=400;endfori=5/7*nodY+2:nodYLs(i)=0.15;endLs(5/7*nodY+1)=2/(1/Ls(5/7*nodY)+1/Ls(5/7*nodY+2));%设定初始值fori=1:nodXforj=1:nodYT(i,j)=23;endendT0=T;%定义内热源fori=1:nodXforj=1:5/7*nodYSp(i,j)=-28.37;Sc(i,j)=6525.08+56.74*T0(i,j);%用T0表示上时刻的值endforj=5/7*nodY+1:nodYSp(i,j)=0;Sc(i,j)=0;endend%计算系数aE=ones(nodX,1)*(DY.*Le./dX);aW=ones(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));fori=1:nodXaN(i,nodY)=0;aS(i,1)=0;endfori=1:nodX/4Sc(i,nodY)=Sc(i,nodY)+Scad;10Sp(i,nodY)=Sp(i,nodY)+Spad;endfori=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;fori=1:nodXforj=1:5/7*nodYSc(i,j)=6525.08+56.74*T0(i,j);%用T0表示上时刻的值endforj=5/7*nodY+1:nodYSc(i,j)=0;endendfori=1:nodX/4Sc(i,nodY)=Sc(i,nodY)+Scad;endfori=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;fori=1:nodX+1forj=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;endendfori=1:nodX+1forj=1:nodYT(nodX+1,j)=T(1,j);Z(i,j)=T(i,j);endend11surf(X,Y,Z);4.2Grid.mglobalXglobalYglobaldXglobaldYglobalDXglobalDXnglobalDXsglobalDYglobalCvglobalCVglobalnodXglobalnodY%计算半径与尺度系数R=Y/2/nodY:Y/nodY:Y-Y/2/nodY;Rn=R+Y/2/nodY;Rs=R-Y/2/nodY;SX=Y/2/nodY:Y/nodY:Y-Y/2/nodY;%计算节点间距dX=X/nodX.*SX;%东西节点距离数组dY=Y/nodY;%南北节点距离常数%计算导热面积DY=dY;%东西导热面积DX=X/nodX.*R;%南北导热面积DXn=X/nodX*Rn;DXs=X/nodX*Rs;Cv=DY*DX;CV=ones(nodX,1)*Cv;4.3TDMA.mfunction[W]=TDMA(A,B,C,D,direct)%TDMAsolverglobalnodXglobalnodYifdirect==2nod=nodY;elseifdirect==1nod=nodX;endendP(1)=B(1)/A(1);Q(1)=D(1)/A(1);fori=2:nodP(i)=B(i)/(A(i)-C(i)*P(i-1));Q(i)=(D(i)+C(i)*Q(i-1))/(A(i)-C(i)*P(i-1));end12W(nod)=Q(nod);fori=nod-1:-1:1W(i)=P(i)*W(i+1)+Q(i);endend4.4SolutionTDMA.mfunction[]=SolutionTDMA(aE,aW,aN,aS,aP,b)%ADI-TDMASolverglobalTglobalT0globalnodXglobalnodY%首先在Y方向上隐式计算direct==2，X方向为显式，利用T0direct=2;i=1;forj=1:nodYA(j)=aP(i,j);B(j)=aN(i,j);C(j)=aS(i,j);D(j)=aE(i,j)*T0(i+1,j)+aW(i,j)*T0(nodX,j)+b(i,j);endC(1)=0;B(nodY)=0;W=TDMA(A,B,C,D,direct);forj=1:nodYT(i,j)=W(j);endT0=T;fori=2:nodX-1forj=1:nodYA(j)=aP(i,j);B(j)=aN(i,j);C(j)=aS(i,j);D(j)=aE(i,j)*T0(i+1,j)+aW(i,j)*T0(i-1,j)+b(i,j);endC(1)=0;B(nodY)=0;W=TDMA(A,B,C,D,direct);forj=1:nodYT(i,j)=W(j);endT0=T;endi=nodX;13forj=1:nodYA(j)=aP(i,j);B(j)=aN(i,j);C(j)=aS(i,j);D(j)=aE(i,j)*T0(1,j)+aW(i,j)*T0(i-1,j)+b(i,j);endC(1)=0;B(nodY)=0;W=TDMA(A,B,C,D,direct);f