銅芯電纜溫度分布MATLAB計(jì)算模型

1、1.題目如圖 1 所示銅芯電纜，電流為5000A，內(nèi)徑為 10mm，外包材料聚氯乙烯的厚度為2mm，導(dǎo)熱系數(shù)為0.15+0.00013t W/(mK) 。電纜左半邊為絕熱邊界條件，右半邊為第三類邊界條件，空氣溫度為20，絕緣層表面與環(huán)境間的復(fù)合表面?zhèn)鳠嵯禂?shù)為10 W/(m 2 K) 。銅的電阻率為 RtR0 1a t-20 ， R01.75 10 8m ， a0.004 / C ，t 的單位為攝氏度。試通過(guò)數(shù)值方法求解溫度分布。圖 12.編程計(jì)算2.1 控制方程根據(jù)題意，本題為二維穩(wěn)態(tài)導(dǎo)熱問(wèn)題，其控制方程為：1T1TrrrS 0r rr邊界條件：r0.0070 ,3： qBh TfTW 222

2、r0.007 3： qB022其中： Tf20 C。12.2 方程離散為建立通用方程，考慮非穩(wěn)態(tài)項(xiàng)的控制方程為：c T1rT1TS 0tr rrrr采用全隱格式，在t 時(shí)間內(nèi)，對(duì)控制容積積分，整理后可得：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采用通用表達(dá)式，各表達(dá)式如下表：表 1 坐標(biāo)及系數(shù)表達(dá)式坐標(biāo)系極坐標(biāo)通用表達(dá)式東西坐標(biāo)X南北坐標(biāo)rY半徑rR東西尺度系數(shù)rSX東西節(jié)

3、點(diǎn)間距rXSX南北節(jié)點(diǎn)間距rY東西導(dǎo)熱面積rR Y/SX南北導(dǎo)熱面積rRX2控制體體積rrR XYa E aWrRYe r / eSX 2 X e / eaN aSrRXr n /Y n / nnaP0c R X Y /tbSCRXYaPaP aEaW aNaSaP0SPRXY2.3 邊界條件處理對(duì)于北邊界，采用附加源項(xiàng)法處理。由于北邊界( r0.0070 ,32 )為22第三類邊界條件，則最靠近邊界的控制容積加入以下附加源項(xiàng)：SC, adATfV 1/ hx n /BSP ,adA1V 1/ hx n /B其中： Tf 20 C將附加源項(xiàng)加到相應(yīng)控制容積后，再令相應(yīng)的aN 0 。對(duì)于南邊界，

4、可認(rèn)為定溫邊界條件，由于其導(dǎo)熱面積為零，aS0 。對(duì)于東西邊界，計(jì)算時(shí)取 2 計(jì)算區(qū)域，故東西邊界重合，可認(rèn)為為定溫邊界條件，溫度為上一層相鄰控制容積的溫度。2.4 導(dǎo)熱系數(shù)與計(jì)算取銅導(dǎo)熱系數(shù)為常數(shù)，400 W/(mK) 。每個(gè)控制容積各界面對(duì)應(yīng)導(dǎo)熱系數(shù)分別為n 、 s 、 e 、 w 。對(duì)于銅芯或保溫層內(nèi)部控制容積，各導(dǎo)熱系數(shù)均為常數(shù)。兩者交接界面的導(dǎo)熱系數(shù)用調(diào)和平均法計(jì)算。32.5 方程求解方程采用 ADI-TDMA 方法求解，首先在 Y 方向進(jìn)行隱式計(jì)算， X 方向采用顯式計(jì)算。各方向?qū)?yīng)方程為三對(duì)角矩陣，使用 TDMA 法求解。然后再在 X 方向進(jìn)行隱式計(jì)算， Y 方向采用顯式計(jì)算。4

5、3. 結(jié)果輸出與分析3.1 計(jì)算結(jié)果程序中溫度 T 為二維數(shù)組，采用坐標(biāo)變換方法，將溫度表示在極坐標(biāo)系中。設(shè)定溫度初場(chǎng)為 23，循環(huán)結(jié)束判定條件為T(mén)T0 / T 21E9 ，網(wǎng)格數(shù)為 12070條件下，輸出結(jié)果如圖2：圖 23.2 網(wǎng)格獨(dú)立性考察保持迭代精度TT0 /T 21E8不變：1.網(wǎng)格數(shù)為 2414 時(shí)，計(jì)算結(jié)果為：5圖 32.網(wǎng)格數(shù)為 4828時(shí)，計(jì)算結(jié)果為：圖 463. 網(wǎng)格數(shù)為 200 140 時(shí)，計(jì)算結(jié)果為：圖 5結(jié)論：從以上各圖可以看出，程序運(yùn)行結(jié)果與網(wǎng)格劃分無(wú)關(guān)，程序具有較好的網(wǎng)格獨(dú)立性。3.3 收獲與體會(huì)通過(guò)這次 matlab 編程作業(yè)，我對(duì)二維擴(kuò)散問(wèn)題有了更加深刻的理解

6、，對(duì)網(wǎng)格劃分、通用離散形式、邊界條件處理等有了進(jìn)一步的認(rèn)識(shí)。在編寫(xiě) Matlab 程序過(guò)程中，我為了直接求解三對(duì)角矩陣還曾編寫(xiě)一個(gè) Solution.m 文件，經(jīng)過(guò)對(duì)比后發(fā)現(xiàn)此文件相比于 TDMA 方法在速度上稍微快一點(diǎn)，結(jié)果基本相同。通過(guò)編程，我更加深刻的認(rèn)識(shí)到只有親自動(dòng)手才能加深對(duì)問(wèn)題理解，才能真正獲得屬于自己的知識(shí)。74. 程序語(yǔ)句程序采用 Matlab 編寫(xiě)，主要分為4 部分，分別是主程序，用于給定題目條件，調(diào)用其他函數(shù)，循環(huán)求解等；網(wǎng)格劃分函數(shù)Grid.m，用于劃分網(wǎng)格； SolutionTDMA.m ，用于執(zhí)行交替隱式計(jì)算；TDMA.m ，用于求解三對(duì)角矩陣。%給出題目參數(shù)4.1

7、Main.mX=2*pi;clear allY=7E-3;clear globalGrid; %劃分網(wǎng)格format longglobal Xglobal YTf=20;h=10;%計(jì)算導(dǎo)熱系數(shù)global dXfor i=1:5/7*nodYglobal dYglobal DXLe(i)=400;% 假設(shè)銅的導(dǎo)熱系數(shù)為 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上時(shí)刻的值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;%計(jì)算系數(shù)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);%設(shè)定初始值aP0=0;for i=1:nodXfor j=1:nodY%邊界條件的處理，附加源項(xiàng)法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%定義內(nèi)熱源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表示上時(shí)刻的值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;% 南北導(dǎo)熱面積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%計(jì)算半徑與尺度系數(shù)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%計(jì)算節(jié)點(diǎn)間距dX=X/nodX.*SX;%東西節(jié)點(diǎn)距離P(1)=B(1)/A(1);數(shù)組Q(1)=D(1)/A(1);dY=Y/nodY;%南北節(jié)點(diǎn)for i=2:nod距離 常數(shù)P(i)=B(i)/(A(i)-C(i)*P(i-1);%計(jì)算導(dǎo)熱面積Q(i)=(D(i)+C(i

15、)*Q(i-1)/(A(i)-C(i)*P(i-DY=dY;%東西導(dǎo)熱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 方向上隱式計(jì)算enddirect=2， X 方向?yàn)轱@式，利用 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)%把計(jì)算得到的 T 賦給 T0，然后進(jìn)行+b(i,j);X 方向隱式計(jì)算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,