数值传热学埋地电缆二维温度场分布的研究

中国石油大学（华东）本科毕业设计（论文）数值传热学大作业题目：埋地电缆二维温度场分布的研究学生姓名：学 号：S专业班级：油气储运工程2016年 6月4日数值传热学大作业数值传热大作业摘要：对于二维电缆的导热问题进行了数值计算，采用通用离散形式，加入初始化条件以提高收敛速度。计算了不同电流下电缆的温度分布，通过网格考核和结果分析认为数值计算结果是可靠的。问题描述：如图所示铜芯电缆，电流为20A，内径为10mm，外包材料为聚氯乙烯的厚度为2mm，导热系数为0.15W/(mK)。铜的电阻率Rt=R01+at-20，R0=1.7510-8m，a=0.0040/，t的单位为。试通过数值方法求解以下几种情况温度分布。（1）电缆完全埋地，土壤导热系数为：0.5 W/(m )，土壤温度为16；（2）电缆一半埋于地下，一半裸露于空气中，空气温度16；（3）电缆完全暴露于空气中，空气温度为16，绝缘层表面与环境间的复合表面传热系数为10。环境条件1环境条件3环境条件2图1 电缆埋地方式示意图解题目标：（1）使用通用离散形式离散控制方程（2）通过ADI TDMA算法求解（3）完成报告（包括控制方程和边界条件、离散过程、边界条件处理、模型验证、网格独立性考核及结果分析）。1. 控制方程和边界条件（以第二种情况为例）1.1 控制方程一般变截面二维稳态导热问题控制方程的通用形式可表示为1A(x)xA(x)Tx+1B(y)yB(y)Ty+S=0该问题宜采用极坐标，形式如下1rrrTr+1rrT+S=01.2 边界条件该计算域关于y轴对称，取其中一半区域计算。由于圆心为奇异点，在中心位置挖去一个半径为m0的小圆面，将该边界设为绝热。边界条件如下r=R，2，Tr=h(tf-tN)32，Tr=0；r=mo，Tr=0=2，T=0；=32，T=02. 方程离散2.1 通用控制方程离散计算采用内节点法，方程离散采用控制容积积分法，假定控制容积界面上热流密度是均匀的，采用全隐格式，在时间内，对控制容积积分，并引入东西尺度系数SX和名义半径R，离散过程如下：扩散项：tt+tsnwe1RYSX1SXXRYSX1SXTX+1RXYRXTYdXdYdt=eReSXe2TE-TPXe-wRwSXw2TP-TWXwYt+nRnTN-TPYe-sRsTP-TSYwXt源项：tt+tsnweSdXdYdt=(SC+SPTP)XYt整理结果：aE=Re(Y)SX2(x)e/e，aW=Rw(Y)SX2(x)w/w，aN=Rn(X)(Y)n/n，aN=Rs(X)(Y)s/saP=aE+aW+aN+aS-SpRXY，b=ScRXY2.2 三种坐标系下系数通用表达式表1三种坐标系中系数的通用表达式 坐标系极坐标通用表达式东西坐标南北坐标半径东西尺度系数东西节点间距南北节点间距东西导热面积南北导热面积控制体体积3. 边界条件处理边界条件处理采用附加源项法。3.1 径向边界r=R，2，Tr=h(tf-tN)32，Tr=0；r=mo，Tr=0（1）r=R，32边界绝热，qB=0，则TN=TPaP-aNTp=aETE+aWTW+aN(TN-TP)+aSTS+baPTp=aETE+aWTW+aSTS+b，aP=aP-aN（2）r=R，2边界为第三类边界条件aP-aNTp=aETE+aWTW+aN(TN-TP)+aSTS+bqB=TN-TP(r)n=hTf-TN=Tf-TP(r)n/+1/h整理得：aP+A1h+(r)nTp=aETE+aWTW+aSTS+b+ATf1h+(r)n，aP=aP-aN（3）r=m边界绝热，则TS=TPaP-aSTP=aETE+aWTW+aNTN+aS(TS-TP)+baPTP=aWTW+aNTN+aETE+b，aP=aP-aS3.2 对称边界=2，T=0；=32，T=0（1）=32边界绝热，则TW=TPaP-aWTP=aETE+aW(TW-TP)+aNTN+aSTS+baPTP=aETE+aNTN+aSTS+b，aP=aP-aW（2）=2边界绝热，则TE=TPaP-aETP=aE(TE-TP)+aWTW+aNTN+aSTS+baPTP=aWTW+aNTN+aSTS+b，aP=aP-aE3.3 内热源电流通过铜芯产生热量，单位长度的热流量Q=I2R=I2R01+a(TP-20)R2电缆横截面上，单位面积内热源表示为S=SC+SPTP=QR2=I2R01+a(TP-20)2R4可得SC=I2R0(1-20a)2R4，SP=I2R0a2R44. 第三类边界条件的推导4.1 分析论证分析电缆导热的物理本质，由能量守恒定律，当导热达到稳态时，内热源产生的热量与边界的散热相同。该物理模型中只有电缆右侧边界向外散热，其它边界均为绝热。单位长度上，内热源产生热量Q=I2R=I2R01+a(TP-20)R2内热源产生热量全部从右侧边界流出，热流密度为q=QRw=I2R01+a(TP-20)2R2Rw第三类边界条件热流密度公式qB=h(tf-tN)内部边界向外散热时温度逐渐降低，而铜的导热系数很大，近似认为铜线内部温度相同。散热边界温度应低于铜芯温度，设温差为0，边界温度为T，则铜芯温度为T+。由能量守恒hT-tf=I2R01+a(T+-20)2R2Rw整理得Th-I2R0a2R2Rw=I2R01+a-202R2Rw+htf4.2 初始化条件的设置计算过程中，需要在网格节点上采用ADI TDMA对离散方程进行迭代。提供一个合理的初场，对迭代的收敛速度非常重要。采用4.1中思想，考虑达到稳态时全场温差不会太大，为简化计算，假设计算域内温度相同。取=0，可以得到全场初始温度T=I2R01-20a2R2Rw+htfh-I2R0a2R2Rw5 . 网格考核数值计算中网格划分采用内节点法，先划分界面，下面mn为x和y方向的界面数。计算仅选择了上半部分，迭代收敛后将计算结果沿y轴镜像后输出。收敛条件为：（1）两次迭代之间全场温度的平均误差小于10-7；（2）若精度不能达到10-7，但两次迭代误差值趋于稳定。收敛条件=1mni,j=0i=mj=n(Ti,j-Ti,j0)210-710-7， =-0括号内数字表示达到收敛所需迭代次数，计算过程中输出迭代次数与温度收敛情况。题目中没有给出铜的导热系数，查阅相关文献，计算中取为400W/(mK)。考虑在铜芯与外包层界面处导热系数有统一的计算方法，取此面为界面，通过挖除奇异点的大小m与Y方向界面数目控制界面的位置。考虑包覆层外侧边界条件不同，不同边界条件采用界面分隔，因此X方向界面数取为奇数。5.1 97网格计算结果（291240）5.2 1914网格计算结果（311089）5.3 4935网格计算结果（132）5.4 9970网格计算结果（401）5.5 199140网格计算结果（1090）可见，网格划分越细，计算结果越精确。当网格划分至4935时，计算结果与更为细密的网格相差不大，认为得到了网格无关解，且4935网格迭代次数较少，因此下面验证计算中采用此网格。6 .模型验证4.1分别取电流I=50A，200A，300A和400A进行数值计算，验证模型的正确性。6.1 I=50A（109）6.2 I=200A（216）6.3 I=300A（772）6.4 I=400A（257116）可见，设置电流越大，达到收敛所需迭代次数越多。在推导范围内取电流值，计算结果均符合物理规律，认为模型是可靠的。7.结果分析上面通过数值计算求解二维电缆导热问题，给出了不同电流下的网格无关解，计算结果与物理规律符合很好，认为模拟结果是可靠的。从上面数值计算结果可得出以下结论：（1）不论电流大小如何，电缆内部温度分布趋势相似，靠近散热边界温度低，铜芯和绝热边界温度高。（2）铜芯热阻小，温度梯度小，整个铜芯温度相近；外包层热阻大，温度梯度大，从与铜芯的接触面到右侧边界温降较快。（3）电流越大，达到稳态后电缆内部温度越高。8.源代码程序语句程序采用Matlab编写，主要分为4部分，分别是主程序，用于给定题目条件，调用其他函数，循环求解等；网格划分函数Grid.m，用于划分网格；SolutionTDMA.m，用于执行交替隐式计算；TDMA.m，用于求解三对角矩阵。8.1 主程序clc;clear %清空工作空间format longlamatP=0.15+0.00013*t;lamatCu=400;%导热系数赋值XN=49;YN=35; %网格划分XL=pi;mo=0.0002;YL=0.007-mo; %挖除奇异点DX=XL/(XN-1);DY=YL/(YN-1);dX=XL/(XN-1);dY=YL/(YN-1);XF(2:XN+1)=dX*(0:XN-1);YF(2:YN+1)=mo+dY*(0:YN-1);X(1)=pi/2+XF(2);X(XN+1)=pi/2+XF(XN+1);X(2:XN)=(XF(2:XN)+XF(3:XN+1)/2+pi/2;Y(1)=YF(2);Y(YN+1)=YF(YN+1);Y(2:YN)=(YF(2:YN)+YF(3:YN+1)/2;R,SX=cordinate(3); %选择坐标系e=10-7; %设置计算精度for i=2:YN+1 %节点坐标Y处导热系数 if YF(i)5.0001/1000 lamatY(i)=lamatP; else lamatY(i)=2*lamatP*lamatCu/(lamatP+lamatCu);%界面处取调和平均 end endendfor i=2:YN %节点坐标X处导热系数 if Y(i)5/1000 lamatX(i)=lamatCu; else lamatX(i)=lamatP; endendr=YF; %aY=subs(R)*DX.*lamatY/dY;r=Y(1:YN); %计算半径aX=subs(R)*DY.*lamatX./(subs(SX).2*dX);aS(2:YN)=aY(2:YN);aN(2:YN)=aY(3:YN+1);aE=aX;aW=aX;I=400;R0=1.75/108;a=0.004;%电流初始化R=subs(R);for i=2:YN %内热源处理 if YF(i)5/1000 Sc(i)=I2*R0*(1-20*a)/(pi2*(0.0054-mo4); Sp(i)=I2*R0*a/(pi2*(0.0054-mo4); else Sc(i)=0;Sp(i)=0; endendb=Sc.*R*DX*DY;aP=aE+aW+aN+aS-Sp.*R*DX*DY;aPE=aP-aE; %附加源项系数处理aPW=aP-aW;aPS=aP(2)-aS(2);aPN1=aP(YN)-aN(YN);Tf=20;h=10; %初始化条件Ts=(I2*R0*(1-20*a)/(pi2*0.0052*0.007)+h*Tf)/(h-I2*R0*a/(pi2*0.0052*0.007);A=YF(end)*DX;DV=R(end)*DX*DY;Scad1=Tf*A/(1/h+DY/(2*lamatP)*DV);%边界条件处理附加源项Spad1=-A/(1/h+DY/(2*lamatP)*DV);aPN2=aPN1-Spad1*DV;bad=Scad1*DV;f=zeros(YN+1,XN+1); %36*50零矩阵a=f;b1=f;c1=f;d=f;b2=f;c2=f;for i=2:YN % for j=2:XN a(i,j)=aP(YN+2-i); b1(i,j)=aE(YN+2-i); c1(i,j)=aW(YN+2-i); b2(i,j)=aS(YN+2-i); c2(i,j)=aN(YN+2-i); d(i,j)=b(YN+2-i); endenda(3:YN-1,XN)=aPE(YN-1:-1:3);%绝热边界a(3:YN-1,2)=aPW(YN-1:-1:3);for i=2:XN if X(i)pi a(2,i)=aPN1; %对流边界 else a(2,i)=aPN2;d(2,i)=bad;%对流边界 endenda(YN,3:XN-1)=aPS;a(2,2)=aPN1-aW(end);a(2,XN)=aPN2-aE(end);a(YN,2)=aPS-aW(2);a(YN,XN)=aPS-aE(2);c1(:,2)=0;b1(:,XN)=0;b2(YN,:)=0;c2(2,:)=0;f(2:YN,2:XN)=Ts;%f=ADI(a,b1,c1,d,b2,c2,f,e); % ADI算法求解f(:,1)=f(:,2);f(:,XN+1)=f(:,XN);for i=1:XN+1 if X(i)e f0=f;f1=f;k=k+1 for i=2:m-1 for j=2:n-1 d1(j)=c2(i,j)*f0(i-1,j)+b2(i,j)*f0(i+1,j)+d(i,j); end f1(i,:)=TDMA(a(i,:),b1(i,:),c1(i,:),d1); end f=f1; for j=2:n-1 for i=2:m-1 d2(i)=b1(i,j)*f1(i,j+1)+c1(i,j)*f1(i,j-1)+d(i,j); end f(:,j)=TDMA(a(:,j),b2(:,j),c2(:,j),d2); end A=norm(f(:)-f0(:)/(m-2)*(n-2)end8.3 坐标系选择函数function R,SX=zuobiaoxi(n)syms rswitch n case 1, R=1;SX=1; case 2, R=r;SX=1; case 3, R=r;SX=r;end8.4 TDMA函数function T=TDMA(a,b,c,d)N=length(a);p(2)=b(2)/a(2);q(2)=d(2)/a(2);for i=3:N-1 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);endT(N-1)=q(N-1);