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

1、51題目如圖1所示銅芯電纜，電流為5000A，內(nèi)徑為10mm,外包材料聚氯乙烯的厚度為2mm,導(dǎo)熱系數(shù)為0.15+0.00013t W/(m K)。電纜左半邊為絕熱邊界條件，右半邊為第三類(lèi)邊界條件，空氣溫度為20E，絕緣層表面與環(huán)境間的復(fù)合表面?zhèn)鳠嵯禂?shù)為10W/(m2 K)。銅的電阻率為 Rt 二 R°1 at-20 1, R)=1.75 10* 門(mén) m，a = 0.004/C， t的單位為攝氏度。試通過(guò)數(shù)值方法求解溫度分布。維執(zhí)L M " T. *、圖12編程計(jì)算2.1控制方程根據(jù)題意,本題為二維穩(wěn)態(tài)導(dǎo)熱問(wèn)題，其控制方程為:1 “侗'1 M"t+s=0邊

2、界條件：二 3 二r =0.007v _ 0, 2-:2 2a313:小r =0.007:qB - 022其中：Tf=20 C。qB = h Tf -Tw2.2方程離散為建立通用方程，考慮非穩(wěn)態(tài)項(xiàng)的控制方程為:-c汀.:t1譏町+s = or扎I + 曲丿r胡7胡丿采用全隱格式，在t時(shí)間內(nèi)，對(duì)控制容積積分，整理后可得:其中:apTp = QeTe OwTw QnTn a$Ts b3er：U e / 'eAr5 宀 w/ 'waNrn- r n/'nap = aE aw aN as aP -S VaPV，b =Sc ：V a；TP0,V = 0.5 rn rs r心二 t

3、米用通用表達(dá)式，各表達(dá)式如下表：表1坐標(biāo)及系數(shù)表達(dá)式坐標(biāo)系極坐標(biāo)通用表達(dá)式東西坐標(biāo)eX南北坐標(biāo)rY半徑rR東西尺度系數(shù)rSX東西節(jié)點(diǎn)間距r56(SXjSX)南北節(jié)點(diǎn)間距5r6Y東西導(dǎo)熱面積ArR 也Yysx南北導(dǎo)熱面積rAQr9x )控制體體積R(AX QY )a E (aW )irR3)(30 er"e(SX張X/丸eaN (aS )rA0RQX )屮n(6Y / 打0aP(Pc RXAY MtbSC R也X心丫apaP = aE + aW * a【N +as +a； -SPR占XAY2.3邊界條件處理對(duì)于北邊界，采用附加源項(xiàng)法處理。由于北邊界(r二0.007 八0，3 2二)為2

4、 2第三類(lèi)邊界條件，則最靠近邊界的控制容積加入以下附加源項(xiàng)：Sc,adATf二 V 1/ h 亠 EX n / ' BA1V 1/ h 亠 2x n / ' B其中：Tf =20 C將附加源項(xiàng)加到相應(yīng)控制容積后，再令相應(yīng)的aN =0對(duì)于南邊界，可認(rèn)為定溫邊界條件，由于其導(dǎo)熱面積為零，as =0對(duì)于東西邊界，計(jì)算時(shí)取2二計(jì)算區(qū)域，故東西邊界重合，可認(rèn)為為定溫邊界條件, 溫度為上一層相鄰控制容積的溫度。2.4導(dǎo)熱系數(shù)與計(jì)算取銅導(dǎo)熱系數(shù)為常數(shù)，一 400W/(m K) o每個(gè)控制容積各界面對(duì)應(yīng)導(dǎo)熱系數(shù)分別為n、 s > e、'w o對(duì)于銅芯或保溫層內(nèi)部控制容積，各導(dǎo)熱系

5、數(shù)均為常數(shù)。兩者交接界面的導(dǎo)熱系數(shù)用調(diào)和平均法計(jì)算。2.5 方程求解X 方向采用顯式計(jì)X 方向進(jìn)行隱式計(jì)方程采用 ADI-TDMA 方法求解，首先在 Y 方向進(jìn)行隱式計(jì)算， 算。各方向?qū)?yīng)方程為三對(duì)角矩陣，使用 TDMA 法求解。然后再在 算，丫方向采用顯式計(jì)算。3.結(jié)果輸出與分析3.1計(jì)算結(jié)果將溫度表示在極坐標(biāo)系中。-To /T2 ：佢-9 ，網(wǎng)格數(shù)為120 70程序中溫度T為二維數(shù)組，采用坐標(biāo)變換方法，設(shè)定溫度初場(chǎng)為23 C,循環(huán)結(jié)束判定條件為 T條件下，輸出結(jié)果如圖2:93.2網(wǎng)格獨(dú)立性考察保持迭代精度 T -To /T 2 : 1E -8不變:1.網(wǎng)格數(shù)為24 14時(shí)，計(jì)算結(jié)果為：22

6、 9522 922 8522 622 7622 722.65圖32.網(wǎng)格數(shù)為48 28時(shí)，計(jì)算結(jié)果為：x 10S W|III”呂IIIIIIII36-202468x 10'2322 9522 922 8522 8227522 722 653.網(wǎng)格數(shù)為200 140時(shí)，計(jì)算結(jié)果為:-2Z8圖5-二?-21'?結(jié)論：從以上各圖可以看出，程序運(yùn)行結(jié)果與網(wǎng)格劃分無(wú)關(guān)，程序具有較好的網(wǎng)格 獨(dú)立性。3.3收獲與體會(huì)通過(guò)這次matlab編程作業(yè)，我對(duì)二維擴(kuò)散問(wèn)題有了更加深刻的理解，對(duì)網(wǎng)格劃分、通用離散形式、邊界條件處理等有了進(jìn)一步的認(rèn)識(shí)。在編寫(xiě)Matlab程序過(guò)程中，我為了直接求解三對(duì)角矩陣

7、還曾編寫(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í)。4. 程序語(yǔ)句程序采用 Matlab 編寫(xiě)，主要分為 4 部分，分別是主程序，用于給定題目條件，調(diào) 用其他函數(shù)，循環(huán)求解等；網(wǎng)格劃分函數(shù)Grid.m，用于劃分網(wǎng)格；SolutionTDMA.m ,用于執(zhí)行交替隱式計(jì)算； TDMA.m ，用于求解三對(duì)角矩陣。4.1 Main.m%給出題目參數(shù)X=2*pi;clear allY=7E-3;clear globalGrid; %劃分網(wǎng)格format lo

8、ngTf=20; h=10;global Xglobal Y%計(jì)算導(dǎo)熱系數(shù)global dXfor i=1:5/7*nodYglobal dYLe(i)=400;% 假設(shè)銅的導(dǎo)熱系數(shù)global DX為 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 T0global Tffor i=1:5/7*nodY-1global nodXLn(i)=400;global nodYendfor

9、 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)+1/Ln(5/7*nodY+1);for i=1:5/7*nodYLs(i)=400;endfor i=5/7*nodY+2:nodYLs(i)=0.15;endLs(5/7*nodY+1)=2/(1/Ls(5/7*nodY)+1/Ls(5/7*nodY+2);%設(shè)定初始值for i=1:nodXfor j=1:nodYT(i,j)=23;endendT0=T;%定義內(nèi)熱源for i=1:nodXfor j=1:5/7

10、*nodYSp(i,j)=-28.37;endfor j=5/7*nodY+1:nodYSp(i,j)=0;Sc(i,j)=0;endend%計(jì)算系數(shù)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;%邊界條件的處理，附加源項(xiàng)法Scad=DXn(nodY)/Cv(nodY)*Tf/(1/h+Y/2/nodY/Ln(nodY);Spad=-DXn(nodY)/Cv(nodY)/(1/h+Y/2/n

11、odY/Ln(nodY);for i=1:nodXaN(i,nodY)=0;aS(i,1)=0;endfor i=1:nodX/4#Sc(i,j)=6525.08+56.74*T0(i,j);% 用 T0 表示Sc(i,nodY)=Sc(i,nodY)+Scad;Sp(i,nodY)=Sp(i,nodY)+Spad;Sc(i,nodY)=Sc(i,nodY)+Scad;15endfor i=3*nodX/4+1:nodXSc(i,nodY)=Sc(i,nodY)+Scad;Sp(i,nodY)=Sp(i,nodY)+Spad; endb=Sc.*CV; aP=aE+aW+aN+aS-aP0-S

12、p.*CV;SolutionTDMA(aE,aW,aN,aS,aP,b); cont=1norm(T-T0)./T)while (norm(T-T0)./T)>1E-8) cont=cont+1 norm(T-T0)./T)T0=T;for i=1:nodXfor j=1:5/7*nodYSc(i,j)=6525.08+56.74*T0(i,j); % 用 T0 表示上時(shí)刻的值endfor j=5/7*nodY+1:nodYSc(i,j)=0;endendfor i=1:nodX/4endfor i=3*nodX/4+1:nodX Sc(i,nodY)=Sc(i,nodY)+Scad;e

13、ndb=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+1for j=1:nodYX(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;endendfor i=1:nodX+1for j=1:nodYT(nodX+1,j)=T(1,j);Z(i,j)=T(i,j);end endsurf(X,Y,Z);4.2 Grid.mgloba

14、l X global Y global dX global dY global DXglobal DXn global DXs global DY global Cv global CVglobal nodX global nodY%計(jì)算半徑與尺度系數(shù)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.*R;% 南北導(dǎo)熱面積DXn=X/nodX*Rn;DXs=X/nodX*Rs;Cv=DY*DX;CV=ones(nodX,1)*Cv;4.3 TDMA.m

15、function W = TDMA (A, B, C, D, direct )%TDMA solverglobal nodXglobal nodYif direct=2nod=nodY;else if direct=1nod=nodX;endend%東西節(jié)點(diǎn)距離%南北節(jié)點(diǎn)%東西導(dǎo)熱面積end%計(jì)算節(jié)點(diǎn)間距 dX=X/nodX.*SX;數(shù)組dY=Y/nodY;距離 常數(shù)%計(jì)算導(dǎo)熱面積DY=dY;P(1)=B(1)/A(1);Q(1)=D(1)/A(1);for i=2:nodP(i)=B(i)/(A(i)-C(i)*P(i-1);Q(i)=(D(i)+C(i)*Q(i-1)/(A(i)-C(i)

16、*P(i-1);W(nod)=Q(nod);for i=nod-1:-1:1W(i)=P(i)*W(i+1)+Q(i);endend4.4 SolutionTDMA.mfunction = SolutionTDMA ( aE,aW, aN,aS,aP,b )%ADI-TDMA Solverglobal Tglobal T0global nodXglobal nodY%首 先 在 Y 方 向上隱式 計(jì)算 direct=2， X 方向?yàn)轱@式，利用 T0direct=2;i=1;for j=1:nodYA(j)=aP(i,j);B(j)=aN(i,j);C(j)=aS(i,j);D(j)=aE(i,

17、j)*T0(i+1,j)+aW(i,j)*T0(nodX,j)+ b(i,j);C(1)=0;B(nodY)=0;W=TDMA(A,B,C,D,direct);for j=1:nodYT(i,j)=W(j);endT0=T;for i=2:nodX-1for j=1:nodYA(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);endC(1)=0;B(nodY)=0;W=TDMA(A,B,C,D,direct);for j=1:nodYT(i,j)=W(j);endT0=T

18、;endendi=nodX;for j=1:nodY1,j)+b(i,j);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);endC(1)=0;B(nodY)=0;W=TDMA(A,B,C,D,direct);for j=1:nodYT(i,j)=W(j);endT0=T;%把計(jì)算得到的T賦給TO,然后進(jìn)行X 方向隱式計(jì)算 direct=1direct=1;j=1;for i=1:nodXA(i)=aP(i,j);B(i)=aE(i,j);C(i)=aW(i,j);D(i)=aN(i,j)*TO(i,j+1)+b(i,j);endD(1)=aN(i,j)*TO(i,j+1)+aW(i,j)*TO(no dX,j)+b(i,j);D(nodX)=aN(i,j)*TO(i,j+1)+aE(i,j)*TO(C(1)=O;B(nodX)=O;W=TDMA(A,B,C,D,direct);for