铜芯电缆温度分布MATLAB计算模型
1. 题目 如图1所示，铜芯电缆电流为5万A，内径为10mm，聚氯乙烯外包材料厚度为2mm，导热系数为0.15 0.00013{t}。电缆左侧为绝缘边界条件，右侧为第三种边界条件，空气温度为20℃，复合表面在绝缘层表面与环境之间的传热系数为10。铜的电阻率为，，，t单位为摄氏度。温度分布试图通过数值法解决。 图 1 2. 编程计算 2.1 控制方程 根据问题的含义，本问题为二维稳态导热问题，其控制方程为： 边界条件： ： ： 其中：。 2.2 方程离散 考虑到非稳态项的控制方程是： 采用全隐格式，及时控制容积分，整理后可得： 其中： ，， ，， ，， 通用表达式，如下表所示： 表1 坐标和系数表达式 坐标系 极坐标 通用表达式 东西坐标 南北坐标 半径 东西尺度系数 东西节点间距 南北节点间距 东西导热面积 南北导热面积 控制体体积 2.3 边界条件处理 对北边界，采用附加源法处理。由于北边界()是第三种边界条件，最接近边界的控制体积增加了以下附加源： 其中： 将附加源项添加到相应的控制容积中，然后使相应的。由于北边界()是第三种边界条件，最接近边界的控制体积增加了以下附加源： 其中： 将附加源项添加到相应的控制容积中，然后使相应的。 由于其导热面积为零，可以将定温边界条件视为南边界。 对于东西边界，计算时取计算区域，因此东西边界重叠，可以认为是恒温边界条件，温度体积的温度。 2.4 导热系数及计算 铜导热系数为常数。 各控制体积各界面对应的导热系数分别为、、、。铜芯或保温层的内部控制体积为常数。通过调和平均法计算两个交接界面的导热系数。 2.5 方程求解 方程采用ADI-TDMA方法求解，首先在Y方向进行隐式计算，X方向采用显式计算。三对角矩阵用于各方向的对应方程TDMA法求解。然后在X方向进行隐式计算，Y方向采用显式计算。 3. 结果输出与分析 3.1 计算结果 在程序中，温度T二维数组，采用坐标变换法在极坐标系中表示温度。 初始设定温度为23℃，在网格数的条件下，输出结果如图2所示： 图 2 3.2 网格独立性考察 迭代精度保持不变： 1. 计算结果为： 图 3 2. 计算结果为： 图 4 3. 计算结果为： 图 5 结论：从上图可以看出，程序运行结果与网格划分无关，网格独立性好。 3.3 收获与体会 通过这次matlab在编程操作中，我对二维扩散有了更深入的了解，对网格划分、一般离散形式、边界条件处理等有了进一步的了解。在编写Matlab在程序过程中，为了直接解决三对角矩阵，我还写了一个Solution.m经过比较，发现该文件与文件相比TDMA方法在速度上稍微快一点，结果基本相同。在编写Matlab在程序过程中，为了直接解决三对角矩阵，我还写了一个Solution.m经过比较，发现该文件与文件相比TDMA方法速度稍快，结果基本相同。 通过编程，我更深刻地认识到，只有自己动手，才能加深对问题的理解，才能真正获得自己的知识。 4. 程序语句 程序采用Matlab编写主要分为四部分，即主程序，用于给定主题条件、调用其他函数、循环求解等；网格划分函数Grid.m，用于划分网格；SolutionTDMA.m，用于交替隐式计算；TDMA.m，用于解决三对角矩阵。 14 4.1 Main.m clear all clear global at long global X global Y global dX global dY global DX global DXn global DXs global DY global Cv global CV global T global T0 global Tf global nodX global nodY nodX=200; nodY=350; 给出主题参数% X=2*pi; Y=7E-3; Grid; %划分网格 Tf=20; h=10; 导热系数计算% for i=1:5/7*nodY Le(i)=400；%假设铜的导热系数为400W/(m.K) Lw(i)=400; end for i=5/7*nodY 1:nodY Le(i)=0.15; Lw(i)=0.15; end for i=1:5/7*nodY-1 Ln(i)=400; end for i=5/7*nodY 1:nodY Ln(i)=0.15; end Ln(5/7*nodY)=2/(1/Ln(5/7*nodY-1) 1/Ln(5/7*nodY 1)); for i=1:5/7*nodY Ls(i)=400; end for i=5/7*nodY 2:nodY Ls(i)=0.15; end Ls(5/7*nodY 1)=2/(1/Ls(5/7*nodY) 1/Ls(5/7*nodY 2)); 初始值设置% for i=1:nodX for j=1:nodY T(i,j)=23; end end T0=T; %定义内热源 for i=1:nodX for j=1:5/7*nodY Sp(i,j)=-28.37; Sc(i,j)=6525.08 56.74*T0(i,j);%用T0表示上一刻的值 end for j=5/7*nodY 1:nodY Sp(i,j)=0; Sc(i,j)=0; end end %计算系数 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)); for i=1:nodX aN(i,nodY)=0; aS(i,1)=0; end for i=1:nodX/4 Sc(i,nodY)=Sc(i,nodY) Scad; Sp(i,nodY)=Sp(i,nodY) Spad; end for i=3*nodX/4 1:nodX Sc(i,nodY)=Sc(i,nodY) Scad; Sp(i,nodY)=Sp(i,nodY) Spad; end b=Sc.*CV; aP=aE aW aN aS-aP0-Sp.*CV; SolutionTDMA(aE,aW,aN,aS,aP,b); cont=1 norm((T-T0)./T) while (norm((T-T0)./T)>1E-8) cont=cont 1 norm((T-T0)./T) T0=T; for i=1:nodX for j=1:5/7*nodY Sc(i,j)=6525.08 56.74*T0(i,j); %用T0表示上一刻的值 end for j=5/7*nodY 1:nodY Sc(i,j)=0; end end for i=1:nodX/4 Sc(i,nodY)=Sc(i,nodY) Scad; end for i=3*nodX/4 1:nodX Sc(i,nodY)=Sc(i,nodY) Scad; end b=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 1 for j=1:nodY X(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; end end for i=1:nodX 1 for j=1:nodY T(nodX 1,j)=T(1,j); Z(i,j)=T(i,j); end end surf(X,Y,Z); 4.2 Grid.m global X global Y global dX global dY global DX global DXn global DXs global DY global Cv global CV global nodX global nodY %计算半径和尺度系数 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.3 TDMA.m function [W] = TDMA (A, B, C, D, direct ) %TDMA solver global nodX global nodY if direct==2 nod=nodY; else if direct==1 nod=nodX; end end P(1)=B(1)/A(1); Q(1)=D(1)/A(1); for i=2:nod P(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)); end W(nod)=Q(nod); for i=nod-1:-1:1 W(i)=P(i)*W(i 1) Q(i); end end 4.4 SolutionTDMA.m function [ ] = SolutionTDMA ( aE,aW, aN,aS,aP,b ) ?I-TDMA Solver global T global T0 global nodX global nodY 首先，Y方向上的隐藏计算direct==2，X方向显式，利用T0 direct=2; i=1; for j=1:nodY A(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); end C(1)=0; B(nodY)=0; W=TDMA(A,B,C,D,direct); for j=1:nodY T(i,j)=W(j); end T0=T; for i=2:nodX-1 for j=1:nodY A(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); end C(1)=0; B(nodY)=0; W=TDMA(A,B,C,D,direct); for j=1:nodY T(i,j)=W(j); end T0=T; end i=nodX; for j=1:nodY A(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-,j)+b(i,j); end C(1)=0; B(nodY)=0; W=TDMA(A,B,C,D,direct); for j=1:nodY T(i,j)=W(j); end T0=T; %把计算得到的T赋给T0，然后进行X方向隐式计算direct==1 direct=1; j=1; for i=1:nodX A(i)=aP(i,j); B(i)=aE(i,j); C(i)=aW(i,j); D(i)=aN(i,j)*T0(i,j+1)+b(i,j); end D(1)=aN(i,j)*T0(i,j+1)+aW(i,j)*T0(nodX,j)+b(i,j); D(nodX)=aN(i,j)*T0(i,j+1)+aE(i,j)*T0(1,j)+b(i,j); C(1)=0; B(nodX)=0; W=TDMA(A,B,C,D,direct); for i=1:nodX T(i,j)=W(i); end T0=T; for j=2:nodY-1 for i=1:nodX A(i)=aP(i,j); B(i)=aE(i,j); C(i)=aW(i,j); D(i)=aN(i,j)*T0(i,j+1)+aS(i,j)*T0(i,j-1)+b(i,j); end D(1)=aN(i,j)*T0(i,j+1)+aS(i,j)*T0(i,j-1)+aW(i,j)*T0(nodX,j)+b(i,j); D(nodX)=aN(i,j)*T0(i,j+1)+aS(i,j)*T0(i,j-1)+aE(i,j)*T0(1,j)+b(i,j); C(1)=0; B(nodX)=0; W=TDMA(A,B,C,D,direct); for i=1:nodX T(i,j)=W(i); end T0=T; end j=nodY; for i=1:nodX A(i)=aP(i,j); B(i)=aE(i,j); C(i)=aW(i,j); D(i)=aS(i,j)*T0(i,j-1)+b(i,j); end D(1)=aS(i,j)*T0(i,j-1)+aW(i,j)*T0(nodX,j)+b(i,j); D(nodX)=aS(i,j)*T0(i,j-1)+aE(i,j)*T0(1,j)+b(i,j); C(1)=0; B(nodX)=0; W=TDMA(A,B,C,D,direct); for i=1:nodX T(i,j)=W(i); end end