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

1、1.题目如图 1 所示铜芯电缆，电流为5000A，内径为 10mm，外包材料聚氯乙烯的厚度为2mm，导热系数为0.15+0.00013t W/(mK) 。电缆左半边为绝热边界条件，右半边为第三类边界条件，空气温度为20，绝缘层表面与环境间的复合表面传热系数为10 W/(m 2 K) 。铜的电阻率为 RtR0 1a t-20 ， R01.75 10 8m ， a0.004 / C ，t 的单位为摄氏度。试通过数值方法求解温度分布。图 12.编程计算2.1 控制方程根据题意，本题为二维稳态导热问题，其控制方程为：1T1TrrrS 0r rr边界条件：r0.0070 ,3： qBh TfTW 222

2、r0.007 3： qB022其中： Tf20 C。12.2 方程离散为建立通用方程，考虑非稳态项的控制方程为：c T1rT1TS 0tr rrrr采用全隐格式，在t 时间内，对控制容积积分，整理后可得：aP TPaE TEaW TWaN TNaSTSb其中：aEr， aWr，e /rww / wreeaNrn， aSrs，r s /r n / nsaP aEaWaN aSaP0SP VaP0c V ， b SCV aP0 TP0 ， V 0.5 rn rs rt采用通用表达式，各表达式如下表：表 1 坐标及系数表达式坐标系极坐标通用表达式东西坐标X南北坐标rY半径rR东西尺度系数rSX东西节

3、点间距rXSX南北节点间距rY东西导热面积rR Y/SX南北导热面积rRX2控制体体积rrR XYa E aWrRYe r / eSX 2 X e / eaN aSrRXr n /Y n / nnaP0c R X Y /tbSCRXYaPaP aEaW aNaSaP0SPRXY2.3 边界条件处理对于北边界，采用附加源项法处理。由于北边界( r0.0070 ,32 )为22第三类边界条件，则最靠近边界的控制容积加入以下附加源项：SC, adATfV 1/ hx n /BSP ,adA1V 1/ hx n /B其中： Tf 20 C将附加源项加到相应控制容积后，再令相应的aN 0 。对于南边界，

4、可认为定温边界条件，由于其导热面积为零，aS0 。对于东西边界，计算时取 2 计算区域，故东西边界重合，可认为为定温边界条件，温度为上一层相邻控制容积的温度。2.4 导热系数与计算取铜导热系数为常数，400 W/(mK) 。每个控制容积各界面对应导热系数分别为n 、 s 、 e 、 w 。对于铜芯或保温层内部控制容积，各导热系数均为常数。两者交接界面的导热系数用调和平均法计算。32.5 方程求解方程采用 ADI-TDMA 方法求解，首先在 Y 方向进行隐式计算， X 方向采用显式计算。各方向对应方程为三对角矩阵，使用 TDMA 法求解。然后再在 X 方向进行隐式计算， Y 方向采用显式计算。4

5、3. 结果输出与分析3.1 计算结果程序中温度 T 为二维数组，采用坐标变换方法，将温度表示在极坐标系中。设定温度初场为 23，循环结束判定条件为TT0 / T 21E9 ，网格数为 12070条件下，输出结果如图2：图 23.2 网格独立性考察保持迭代精度TT0 /T 21E8不变：1.网格数为 2414 时，计算结果为：5图 32.网格数为 4828时，计算结果为：图 463. 网格数为 200 140 时，计算结果为：图 5结论：从以上各图可以看出，程序运行结果与网格划分无关，程序具有较好的网格独立性。3.3 收获与体会通过这次 matlab 编程作业，我对二维扩散问题有了更加深刻的理解

6、，对网格划分、通用离散形式、边界条件处理等有了进一步的认识。在编写 Matlab 程序过程中，我为了直接求解三对角矩阵还曾编写一个 Solution.m 文件，经过对比后发现此文件相比于 TDMA 方法在速度上稍微快一点，结果基本相同。通过编程，我更加深刻的认识到只有亲自动手才能加深对问题理解，才能真正获得属于自己的知识。74. 程序语句程序采用 Matlab 编写，主要分为4 部分，分别是主程序，用于给定题目条件，调用其他函数，循环求解等；网格划分函数Grid.m，用于划分网格； SolutionTDMA.m ，用于执行交替隐式计算；TDMA.m ，用于求解三对角矩阵。%给出题目参数4.1

7、Main.mX=2*pi;clear allY=7E-3;clear globalGrid; %划分网格format longglobal Xglobal YTf=20;h=10;%计算导热系数global dXfor i=1:5/7*nodYglobal dYglobal DXLe(i)=400;% 假设铜的导热系数为 400W/(m.K)global DXnLw(i)=400;global DXsendglobal DYfor i=5/7*nodY+1:nodYglobal CvLe(i)=0.15;global CVLw(i)=0.15;global Tendglobal T0globa

8、l Tffor i=1:5/7*nodY-1global nodXLn(i)=400;global nodYendfor i=5/7*nodY+1:nodYnodX=200;Ln(i)=0.15;nodY=350;8end上时刻的值Ln(5/7*nodY)=2/(1/Ln(5/7*nodY-1)+1end/Ln(5/7*nodY+1);for j=5/7*nodY+1:nodYSp(i,j)=0;for i=1:5/7*nodYSc(i,j)=0;Ls(i)=400;endendendfor i=5/7*nodY+2:nodYLs(i)=0.15;%计算系数endaE=ones(nodX,1)

9、*(DY.*Le./dX);Ls(5/7*nodY+1)=2/(1/Ls(5/7*nodY)+1aW=ones(nodX,1)*(DY.*Lw./dX);/Ls(5/7*nodY+2);aN=ones(nodX,1)*(DXn.*Ln/dY);aS=ones(nodX,1)*(DXs.*Ls/dY);%设定初始值aP0=0;for i=1:nodXfor j=1:nodY%边界条件的处理，附加源项法T(i,j)=23;Scad=DXn(nodY)/Cv(nodY)*Tf/(1/h+endY/2/nodY/Ln(nodY);endSpad=-DXn(nodY)/Cv(nodY)/(1/h+Y/2

10、/nodY/Ln(nodY);T0=T;for i=1:nodX%定义内热源aN(i,nodY)=0;for i=1:nodXaS(i,1)=0;for j=1:5/7*nodYendSp(i,j)=-28.37;for i=1:nodX/4Sc(i,j)=6525.08+56.74*T0(i,j);% 用 T0 表示Sc(i,nodY)=Sc(i,nodY)+Scad;9Sp(i,nodY)=Sp(i,nodY)+Spad;Sc(i,nodY)=Sc(i,nodY)+Scad;endendfor i=3*nodX/4+1:nodXfor i=3*nodX/4+1:nodXSc(i,nodY)

11、=Sc(i,nodY)+Scad;Sc(i,nodY)=Sc(i,nodY)+Scad;Sp(i,nodY)=Sp(i,nodY)+Spad;endendb=Sc.*CV;b=Sc.*CV;aP=aE+aW+aN+aS-aP0-Sp.*CV;aP=aE+aW+aN+aS-aP0-Sp.*CV;SolutionTDMA(aE,aW,aN,aS,aP,b);endSolutionTDMA(aE,aW,aN,aS,aP,b);cont=1%输出图形norm(T-T0)./T)theta=X/nodX;dR=Y/nodY;while (norm(T-T0)./T)>1E-8)for i=1:no

12、dX+1cont=cont+1for j=1:nodYnorm(T-T0)./T)X(i,j)=j*dR*cos(theta*(i-1);T0=T;Y(i,j)=j*dR*sin(theta*(i-1);for i=1:nodXX(i,1)=0;for j=1:5/7*nodYY(i,1)=0;Sc(i,j)=6525.08+56.74*T0(i,j); % 用 T0end表示上时刻的值endendfor i=1:nodX+1for j=5/7*nodY+1:nodYfor j=1:nodYSc(i,j)=0;T(nodX+1,j)=T(1,j);endZ(i,j)=T(i,j);endend

13、for i=1:nodX/4end10surf(X,Y,Z);DX=X/nodX.*R;% 南北导热面积4.2 Grid.mDXn=X/nodX*Rn;DXs=X/nodX*Rs;global XCv=DY*DX;global YCV=ones(nodX,1)*Cv;global dXglobal dYglobal DX4.3 TDMA.mglobal DXnfunction W = TDMA (A, B, C, D,global DXsdirect )global DY%TDMA solverglobal Cvglobal CVglobal nodXglobal nodXglobal nod

14、Yglobal nodYif direct=2%计算半径与尺度系数nod=nodY;R=Y/2/nodY:Y/nodY:Y-Y/2/nodY;else if direct=1Rn=R+Y/2/nodY;nod=nodX;Rs=R-Y/2/nodY;endSX=Y/2/nodY:Y/nodY:Y-Y/2/nodY;end%计算节点间距dX=X/nodX.*SX;%东西节点距离P(1)=B(1)/A(1);数组Q(1)=D(1)/A(1);dY=Y/nodY;%南北节点for i=2:nod距离 常数P(i)=B(i)/(A(i)-C(i)*P(i-1);%计算导热面积Q(i)=(D(i)+C(i

15、)*Q(i-1)/(A(i)-C(i)*P(i-DY=dY;%东西导热1);面积end11W(nod)=Q(nod);C(1)=0;for i=nod-1:-1:1B(nodY)=0;W(i)=P(i)*W(i+1)+Q(i);W=TDMA(A,B,C,D,direct);endfor j=1:nodYendT(i,j)=W(j);end4.4 SolutionTDMA.mT0=T;function = SolutionTDMA ( aE,aW,aN,aS,aP,b )for i=2:nodX-1%ADI-TDMA Solverfor j=1:nodYA(j)=aP(i,j);global T

16、B(j)=aN(i,j);global T0C(j)=aS(i,j);global nodXglobal nodYD(j)=aE(i,j)*T0(i+1,j)+aW(i,j)*T0(i-1,j)+b(i,j);%首先在 Y 方向上隐式计算enddirect=2， X 方向为显式，利用 T0C(1)=0;direct=2;B(nodY)=0;i=1;W=TDMA(A,B,C,D,direct);for j=1:nodYfor j=1:nodYA(j)=aP(i,j);T(i,j)=W(j);B(j)=aN(i,j);endC(j)=aS(i,j);T0=T;endD(j)=aE(i,j)*T0(

17、i+1,j)+aW(i,j)*T0(nodX,j)+b(i,j);i=nodX;end12for j=1:nodY1,j)+b(i,j);A(j)=aP(i,j);C(1)=0;B(j)=aN(i,j);B(nodX)=0;C(j)=aS(i,j);W=TDMA(A,B,C,D,direct);D(j)=aE(i,j)*T0(1,j)+aW(i,j)*T0(i-1,j)+for i=1:nodXb(i,j);T(i,j)=W(i);endendC(1)=0;T0=T;B(nodY)=0;W=TDMA(A,B,C,D,direct);for j=2:nodY-1for j=1:nodYfor i

18、=1:nodXT(i,j)=W(j);A(i)=aP(i,j);endB(i)=aE(i,j);T0=T;C(i)=aW(i,j);D(i)=aN(i,j)*T0(i,j+1)+aS(i,j)*T0(i,j-1)%把计算得到的 T 赋给 T0，然后进行+b(i,j);X 方向隐式计算direct=1enddirect=1;D(1)=aN(i,j)*T0(i,j+1)+aS(i,j)*T0(i,j-1j=1;)+aW(i,j)*T0(nodX,j)+b(i,j);for i=1:nodXD(nodX)=aN(i,j)*T0(i,j+1)+aS(i,j)*T0(i,j-1)A(i)=aP(i,j);+aE(i,j)*T0(1,j)+b(i,j);B(i)=aE(i,j);C(1)=0;C(i)=aW(i,j);B(nodX)=0;D(i)=aN(i,j)*T0(i,j+1)+b(i,j);W=TDMA(A,B,C,D,direct);endfor i=1:nodXD(1)=aN(i,j)*T0(i,j+1)+aW(i,j)*T0(noT(i,j)=W(i);dX,j)+b(i,j);endD(nodX)=aN(i,