铜芯电缆温度分布MATLAB计算模型.docx
《铜芯电缆温度分布MATLAB计算模型.docx》由会员分享，可在线阅读，更多相关《铜芯电缆温度分布MATLAB计算模型.docx（16页珍藏版）》请在冰豆网上搜索。
铜芯电缆温度分布MATLAB计算模型
1.题目
如图1所示铜芯电缆，电流为5000A，内径为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方法在速度上稍微快一点，结果基本相同。
通过编程，我更加深刻的认识到只有亲自动手才能加深对问题理解，才能真正获得属于自己的知识。
4.
程序语句
程序采用Matlab编写，主要分为4部分，分别是主程序，用于给定题目条件，调用其他函数，循环求解等；网格划分函数Grid.m，用于划分网格；SolutionTDMA.m，用于执行交替隐式计算；TDMA.m，用于求解三对角矩阵。
4.1Main.m
clearall
clearglobal
formatlong
globalX
globalY
globaldX
globaldY
globalDX
globalDXn
globalDXs
globalDY
globalCv
globalCV
globalT
globalT0
globalTf
globalnodX
globalnodY
nodX=200;
nodY=350;
%给出题目参数
X=2*pi;
Y=7E-3;
Grid;%划分网格
Tf=20;
h=10;
%计算导热系数
fori=1:
5/7*nodY
Le（i）=400;%假设铜的导热系数为400W/（m.K）
Lw（i）=400;
end
fori=5/7*nodY+1:
nodY
Le（i）=0.15;
Lw（i）=0.15;
end
fori=1:
5/7*nodY-1
Ln（i）=400;
end
fori=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））;
fori=1:
5/7*nodY
Ls（i）=400;
end
fori=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））;
%设定初始值
fori=1:
nodX
forj=1:
nodY
T（i,j）=23;
end
end
T0=T;
%定义内热源
fori=1:
nodX
forj=1:
5/7*nodY
Sp（i,j）=-28.37;
Sc（i,j）=6525.08+56.74*T0（i,j）;%用T0表示上时刻的值
end
forj=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））;
fori=1:
nodX
aN（i,nodY）=0;
aS（i,1）=0;
end
fori=1:
nodX/4
Sc（i,nodY）=Sc（i,nodY）+Scad;
Sp（i,nodY）=Sp（i,nodY）+Spad;
end
fori=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;
fori=1:
nodX
forj=1:
5/7*nodY
Sc（i,j）=6525.08+56.74*T0（i,j）;%用T0表示上时刻的值
end
forj=5/7*nodY+1:
nodY
Sc（i,j）=0;
end
end
fori=1:
nodX/4
Sc（i,nodY）=Sc（i,nodY）+Scad;
end
fori=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;
fori=1:
nodX+1
forj=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
fori=1:
nodX+1
forj=1:
nodY
T（nodX+1,j）=T（1,j）;
Z（i,j）=T（i,j）;
end
end
surf（X,Y,Z）;
4.2Grid.m
globalX
globalY
globaldX
globaldY
globalDX
globalDXn
globalDXs
globalDY
globalCv
globalCV
globalnodX
globalnodY
%计算半径与尺度系数
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.m
function[W]=TDMA（A,B,C,D,direct）
%TDMAsolver
globalnodX
globalnodY
ifdirect==2
nod=nodY;
elseifdirect==1
nod=nodX;
end
end
P
（1）=B
（1）/A
（1）;
Q
（1）=D
（1）/A
（1）;
fori=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）;
fori=nod-1:
-1:
1
W（i）=P（i）*W（i+1）+Q（i）;
end
end
4.4SolutionTDMA.m
function[]=SolutionTDMA（aE,aW,aN,aS,aP,b）
%ADI-TDMASolver
globalT
globalT0
globalnodX
globalnodY
%首先在Y方向上隐式计算direct==2，X方向为显式，利用T0
direct=2;
i=1;
forj=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）;
forj=1:
nodY
T（i,j）=W（j）;
end
T0=T;
fori=2:
nodX-1
forj=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）;
forj=1:
nodY
T（i,j）=W（j）;
end
T0=T;
end
i=nodX;
forj=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-1,j）+b（i,j）;
end
C
（1）=0;
B（nodY）=0;
W=TDMA（A,B,C,D,direct）;
forj=1:
nodY
T（i,j）=W（j）;
end
T0=T;
%把计算得到的T赋给T0，然后进行X方向隐式计算direct==1
direct=1;
j=1;
fori=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）;
fori=1:
nodX
T（i,j）=W（i）;
end
T0=T;
forj=2:
nodY-1
fori=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）;
fori=1:
nodX
T（i,j）=W（i）;
end
T0=T;
end
j=nodY;
fori=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）;
fori=1:
nodX
T（i,j）=W（i）;
end
end
升级会员 