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方法速度稍快，结果基本相同。
通过编程，我更深刻地认识到，只有自己动手，才能加深对问题的理解，才能真正获得自己的知识。
4. 程序语句
程序采用Matlab编写主要分为四部分，即主程序，用于给定主题条件、调用其他函数、循环求解等；网格划分函数Grid.m，用于划分网格；SolutionTDMA.m，用于交替隐式计算；TDMA.m，用于解决三对角矩阵。
14
4.1 Main.m
clear all
clear global
format 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) Sca;
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 )
%ADI-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-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;
%把计算得到的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
展开阅读全文