并按所給邊界條件進(jìn)行導(dǎo)熱計算其中介質(zhì)導(dǎo)熱系數(shù)λ取20

1、一、 對圖 1 所示偏心半圓環(huán)二維區(qū)域生成適體計算網(wǎng)格，并按所給邊界條件進(jìn)行導(dǎo)熱計算，其中介質(zhì)導(dǎo)熱系數(shù)取 20。邊界條件： T(AD)= 5； T(BC)=5 0；TAB0；TCD 500 wm 2 ；n ()網(wǎng)格劃分：采用 Thomas&Middlecoeff 法對網(wǎng)格線密度和邊界正交性的控制，得到：偏心半圓導(dǎo)熱問題由于偏心半圓環(huán)形狀較為復(fù)雜， 無法直接在物理平面下進(jìn)行計算， 而需要在與物理平面相對應(yīng)的計算平面下進(jìn)行計算， 因而其控制方程也需相應(yīng)的轉(zhuǎn)換到計算平面中。假定 、（ ）均為常數(shù)，首先將控制方程2T2T0轉(zhuǎn)換到 x2y2 平面。得到如下控制方程：12T 2( T )12T0J

2、 22J 2J222T(T2T簡化后得：2)20其中： Jxyxy ，x2y2 ,x xy y ,x2y2 。邊界條件通用方程 :T q ，轉(zhuǎn)換后有：T( )qnJ當(dāng)邊界為絕熱邊界時，則q 0；上邊界有熱源輸入時，q 500。將上述方程離散后即可編程進(jìn)行跌代計算。方程的離散控制方程離散后得到：j2j(T jT j) 2j (T j 1T j 1)ij (T j 1T j 1T j 1T j 1)Tiii 1i 1iiii 1i 1i 1i 14(i ji j )j( xij 1xij 1 )2( yij 1yi j 1) 2i4其中：j( xij1xij1)( xij 1xij 1 ) ( y

3、i j1yi j1 )( yi j 1yij 1)i4j(xij1xi j 1 )2( yij1yij1) 2i4邊界條件離散后得到AB 、CD 邊界點的通用方程：j(T jT j1)ij ( 2qJ2T j 1 )ii 1iiTij2(i ji j )其中絕熱邊界則簡化為：Tiji j (Tij1 Tij1 )i j(0 2Tij 1 )2( ijij )得到溫度場：程序：#include <iostream.h>#include <math.h>#include <stdio.h>void circleboundry(); /偏心半圓環(huán)邊界void TM

4、methodmesh ();/ Thomas&Middlcoeff法網(wǎng)格（改善正交性，并傳遞邊界的疏密度）void heattransfer();/ 穩(wěn)態(tài)導(dǎo)熱計算（上下表面有熱源（或絕熱q 0)）void exportdata();/輸出網(wǎng)格點計算數(shù)據(jù)至.dat 文本#define pi 3.1415926#define Tlb 5.0#define Trb 50.0double const qtop=500;/ 左邊界條件（溫度）/ 右邊界條件（溫度）/頂部熱源double const qbottom=0;double const heatcoe=5;/底部熱源/介質(zhì)導(dǎo)熱系數(shù)int

5、i,j;int const m=41,n=41;double xmn,ymn;/計算平面節(jié)點代號,表示第 i 列第 j 行點分別為 方向 (或 xy 方向 )總節(jié)點數(shù) (取奇數(shù)計算平面i,j 點對應(yīng)物理平面點的x,y 坐標(biāo))double Tmn;int const pre=5000;int step;/ 預(yù)設(shè)總迭代次數(shù)/ 迭代次數(shù)double totalpre;/ 殘差值double const min=1e-3;/ 約定殘差精度void main()circleboundry();/加載二維偏心半圓環(huán)邊界TMmethodmesh ();/ Thomas&Middlcoeff法網(wǎng)格he

6、attransfer();exportdata();/導(dǎo)熱計算（上下表面絕熱）/ 輸出網(wǎng)格點坐標(biāo)數(shù)據(jù)void circleboundry ()/偏心半圓環(huán)邊界/#define xscale 1.02#define yscale 1.08/ 弧長方向步長增長比例/徑向步長增長比例double anglemn;double const r1=2.0;double const r2=4*r1;/內(nèi)環(huán)半徑/外環(huán)半徑for (i=0;i<m;i+)/內(nèi)外圓弧邊界if (xscale=1.0)xi0=r1*(1-cos(pi*i/(m-1);yi0=r1*sin(pi*i/(m-1);xin-1=r

7、2*(0.5-cos(pi*i/(m-1);yin-1=r2*sin(pi*i/(m-1);elseif(i<=(m-1)/2)anglei0=0.5*pi*(1-pow(xscale,i)/(1-pow(xscale,(m-1)/2);xi0=r1*(1-cos(anglei0);yi0=r1*sin(anglei0);xin-1=r2*(0.5-cos(anglei0);yin-1=r2*sin(anglei0);elseanglei0=pi*(1-pow(1/xscale,(i-(m-1)/2)/(2*(1-pow(1/xscale,(m-1)/2)+pi/2; xi0=r1*(1

8、-cos(anglei0);yi0=r1*sin(anglei0);xin-1=r2*(0.5-cos(anglei0);yin-1=r2*sin(anglei0);for(j=0;j<n;j+)/左右邊界if (yscale=1.0)x0j=-0.5*r2*j/(n-1);y0j=0;xm-1j=0.5*r2+r2*j/(n-1);ym-1j=0;elseif(j<=(n-1)/2)x0j=(1-pow(yscale,j)*(-0.5*r2)/(2*(1-pow(yscale,(n-1)/2);y0j=0;xm-1j=0.5*r2+(1-pow(yscale,j)*r2/(2*(

9、1-pow(yscale,(n-1)/2); ym-1j=0;elsex0j=(1-pow(1/yscale,(j-(n-1)/2)*(-0.5*r2)/(2*(1-pow(1/yscale,(n-1)/2)+x0(n-1)/2; y0j=0;xm-1j=(1-pow(1/yscale,(j-(n-1)/2)*r2/(2*(1-pow(1/yscale,(n-1)/2)+xm-1(n-1)/2; ym-1j=0;for (j=1;j<n-1;j+)/內(nèi)點賦初值for(i=1;i<m-1;i+)xij=i;yij=j;void TMmethodmesh ()/ Thomas&

10、Middlcoeff法 /#define w1.8double xxmn,yymn;/ 松弛因子，取值范圍(0,2)/ 計算平面i,j 點對應(yīng)物理平面點的x,y坐標(biāo)中間值double xxxmn,yyymn;double A,B,C;double Wmn,Vmn;/ 超松弛迭代法x,y 坐標(biāo)中間值/ 網(wǎng)格結(jié)構(gòu)參數(shù)/ 源項分別代表 和 double deltax,deltay,totalx,totaly; for (step=1;step<pre;step+)/T&M/ 殘差變量法 +SOR 迭代totalx=0;totaly=0;for (j=1;j<n-1;j+)for(

11、i=1;i<m-1;i+)xxij=xij;yyij=yij;Wi0=(-2)*(yi+10-yi-10)*(yi+10-2*yi0+yi-10)+(xi+10-xi-10)*(xi+10-2*xi0+xi-10)/(pow(xi+10-xi-10),2)+pow(yi+10-yi-10),2); / 0 邊界上的源項Win-1=(-2)*(yi+1n-1-yi-1n-1)*(yi+1n-1-2*yin-1+yi-1n-1)+(xi+1n-1-xi-1 n-1)*(xi+1n-1-2*xin-1+xi-1n-1)/(pow(xi+1n-1-xi-1n-1),2)+pow(yi+1n-1-

12、yi-1n-1),2); / n-1 邊界上的源項 V0j=(-2)*(y0j+1-y0j-1)*(y0j+1-2*y0j+y0j-1)+(x0j+1-x0j-1)*(x0j+1-2*x0j+x0j-1)/(pow(x0j+1-x0j-1),2)+pow(y0j+1-y0j-1),2); / 0Vm-1j=(-2)*(ym-1j+1-ym-1j-1)*(ym-1j+1-2*ym-1j+ym-1j-1)+(xm-1j+1-xm-1j-1)*(xm-1j+1-2*xm-1j+xm-1j-1)/(pow(xm-1j+1-xm-1j-1),2)+pow(y m-1j+1-ym-1j-1),2); /

13、m-1Wij=(Wi0*(n-1-j)+Win-1*j)/(n-1);/內(nèi)點上源項Vij=(V0j*(m-1-i)+Vm-1j*i)/(m-1);A=(pow(xij+1-xij-1),2)+pow(yij+1-yij-1),2)/4.0;B=(xi+1j-xi-1j)*(xij+1-xij-1)+(yi+1j-yi-1j)*(yij+1-yij-1)/4.0; C=(pow(xi+1j-xi-1j),2)+pow(yi+1j-yi-1j),2)/4.0;xxxij=(2*A*(xi+1j+xi-1j)-B*(xi+1j+1-xi+1j-1-xi-1j+1+xi-1j-1)+2*C*(xi j

14、+1+xij-1)+A*Wij*(xi+1j-xi-1j)+C*Vij*(xij+1-xij-1)/(4*(A+C);yyyij=(2*A*(yi+1j+yi-1j)-B*(yi+1j+1-yi+1j-1-yi-1j+1+yi-1j-1)+2*C*(yi j+1+yij-1)+A*Wij*(yi+1j-yi-1j)+C*Vij*(yij+1-yij-1)/(4*(A+C);xij=w*xxxij+(1-w)*xij;yij=w*yyyij+(1-w)*yij;deltax=fabs(xij-xxij);totalx=totalx+deltax;deltay=fabs(yij-yyij);tot

15、aly=totaly+deltay;totalstep=totalx+totaly;printf("step=%d residual=%fn",step,totalstep);if(totalstep<min) break;printf(" 迭代次數(shù) =%dn",step);printf("residual=%fn",total);voidheattransfer()/ 穩(wěn)態(tài)導(dǎo)熱溫度分布計算（上下表面有熱源（或絕熱q 0) ）/#define w1.8double TTmn,TTTmn;double Amn,Bmn,Cmn;/

16、超松弛迭代因子，取值范圍(0,2)/ 溫度計算的中間值/ 分別代表轉(zhuǎn)換至計算平面后方程中的系數(shù), ,double Jmn;/Jacobi因子double heatmn;double delta;double totalpre;/中間值/ 殘差中間值/ 每個迭代步驟的殘差總值step=0;for(j=0;j<n;j+)T0j=Tlb;Tm-1j=Trb;for(j=0;j<n;j+)for(i=1;i<m-1;i+)Tij=Tlb;for(step=1;step<pre;step+)totalstep=0.0;for(j=0;j<n;j+)for(i=1;i<

17、m-1;i+)TTij=Tij;if(j=0)Aij=pow(xij+1-xij,2)+pow(yij+1-yij,2);Cij=(pow(xi+1j-xi-1j,2)+pow(yi+1j-yi-1j,2)/4.0;Jij=0.5*(xi+1j-xi-1j)*(yij+1-yij)-0.5*(xij+1-xij)*(yi+1j-yi-1j); heatij=2*qbottom*Jij/(heatcoe*sqrt(Cij);TTTij=(Aij*(Ti+1j+Ti-1j)+Cij*(2.0*Tij+1+heatij)/(2.0*(Aij+Cij);else if(j=n-1)Aij=pow(xi

18、j-1-xij,2)+pow(yij-1-yij,2);Cij=(pow(xi+1j-xi-1j,2)+pow(yi+1j-yi-1j,2)/4.0;Jij=0.5*(xi+1j-xi-1j)*(yij-yij-1)-0.5*(xij-xij-1)*(yi+1j-yi-1j); heatij=2*qtop*Jij/(heatcoe*sqrt(Cij);TTTij=(Aij*(Ti+1j+Ti-1j)+Cij*(2.0*Tij-1+heatij)/(2.0*(Aij+Cij);elseAij=(pow(xij+1-xij-1,2)+pow(yij+1-yij-1,2)/4.0;Bij=(xi+1

19、j-xi-1j)*(xij+1-xij-1)+(yi+1j-yi-1j)*(yij+1-yij-1)/4.0;Cij=(pow(xi+1j-xi-1j,2)+pow(yi+1j-yi-1j,2)/4.0;TTTij=(2.0*Aij*(Ti+1j+Ti-1j)+2.0*Cij*(Tij+1+Tij-1)-Bij*(Ti+1 j+1-Ti+1j-1-Ti-1j+1+Ti-1j-1)/(4.0*(Aij+Cij);Tij=w*TTTij+(1-w)*Tij;delta=fabs(Tij-TTij);totalstep=totalstep+delta;printf("step=%d tem

20、presidual=%fn",step,totalstep);if(totalstep<min) break;printf("Iteration Complete!n");printf("Iteration=%dn",step);printf("Residual=%fn",totalstep);void exportdata()FILE*fp;if(fp=fopen("temperature.dat","w")=NULL)printf("cannot open file

21、n");return;fprintf(fp,"title=temperaturen");fprintf(fp,"V ARIABLES = X,Y,Tn");fprintf(fp,"ZONE I=%d,J=%d,F=POINTn",m,n);for(j=0;j<n;j+)for(i=0;i<m;i+)fprintf(fp,"%f %f %fn",xij,yij,Tij);fclose(fp);二、 2，采用 SIMPLE 算法求解圖 2 所示區(qū)域（二維平板流動）中的層流流動。圖 2 平板流動區(qū)域

22、由于該邊界形狀較為簡單， 為簡化計算，將直接在其物理平面上采用均分網(wǎng)格進(jìn)行計算。平板層流計算網(wǎng)格示意圖對于壓力驅(qū)動的平板層流問題， 控制方程為動量方程和連續(xù)性方程 （密度為常數(shù)）：(uu)(uv)2u2u1P動量方程：xyx2y2x(vu)(vv)2v2v1Pxyx2y2yu v連續(xù)性方程：0xy編程得出： u 的分布圖：速度 v 的分布圖：壓力分布圖：程序：#include <iostream.h>#include <math.h>#include <fstream.h>#include <stdio.h>#define I 62#define

23、 J 42void TDMA(double*F,double*A,double*B,double*C,double*D,int M);doubleXJI,YJI,xJI,yJI,uJ-1I-1,vJ-1I-1,PJI,UI-1,VJ-1,PPI,xxJI,yyJI;double AI,BI,CI,DI,TJI,TTI;double aPJI,aEJI,aSJI,aNJI,aWJI;double aP1JI,aE1JI,aS1JI,aN1JI,aW1JI;double apJI,aeJI,awJI,anJI,asJI,bJI,pJI;double dx,dy,dxw,dxe,dyn,dys,es

24、p,esp1,esp2;int i,j,k,k1;double m,n,ca;double uuJI,vvJI;double rou=1000,miu=10,alfa1=0.95,alfa2=0.01;int tt=0;/double Max(double a,double b)double z;if(a>b)z=a;elsez=b;return z;void main()fstream outfile1;outfile1.open("file11.dat",ios:out);for(i=0;i<I;i+)X0i=(double)20/(I-2)*i;/ =0/

25、Y0i=0;for(i=0;i<I;i+)XJ-2i=(double)20/(I-2)*i;/ =1/YJ-2i=2;for(j=0;j<J;j+)Xj0=0;Yj0=(double)2/(J-2)*j;for(j=0;j<J;j+)XjI-2=20;YjI-2=(double)2/(J-2)*j;/超限插值 /for(i=1;i<I-2;i+)for(j=1;j<J-2;j+)m=(double)i/(I-2);n=(double)j/(J-2);Xji=(1-m)*Xj0+m*XjI-2+(1-n)*X0i+n*XJ-2i-(1-m)*(1-n)*X00-(1

26、-m)*n*XJ-20-m*(1-n)*X0I-2-m*n*XJ-2I-2;Yji=(1-m)*Yj0+m*YjI-2+(1-n)*Y0i+n*YJ-2i-(1-m)*(1-n)*Y00-(1-m)*n*YJ-20-m*(1-n)*Y0I-2-m*n*YJ-2I-2;dx=(double)20/(I-2);dy=(double)2/(J-2);dxe=dx;dxw=dx;dyn=dy;dys=dy;/計算節(jié)點坐標(biāo)/x00=X00;y00=Y00;yJ-10=YJ-20;xJ-10=XJ-20;xJ-1I-1=XJ-2I-2;yJ-1I-1=YJ-2I-2;x0I-1=X0I-2;y0I-1=Y

27、0I-2;for(i=1;i<I-1;i+)x0i=(X0i-1+X0i)/2;y0i=(Y0i-1+Y0i)/2;xJ-1i=(XJ-2i-1+XJ-2i)/2;yJ-1i=(YJ-2i-1+YJ-2i)/2;for(j=1;j<J-1;j+)yj0=(Yj-10+Yj0)/2;xj0=(Xj-10+Xj0)/2;yjI-1=(Yj-1I-2+YjI-2)/2;xjI-1=(Xj-1I-2+XjI-2)/2;for(j=1;j<J-1;j+)for(i=1;i<I-1;i+)xji=(Xj-1i-1+Xji-1+Xj-1i+Xji)/4;yji=(Yj-1i-1+Yj

28、i-1+Yj-1i+Yji)/4;/速度場 /for(j=1;j<J-1;j+)for(i=1;i<I-2;i+)uji=10;for(i=1;i<I-1;i+)for(j=1;j<J-2;j+)vji=1;/ 邊界條件 /for(j=1;j<J-1;j+)uj0=20;for(j=1;j<J-2;j+)vj0=0;vjI-1=0;for(i=1;i<I-2;i+)u0i=0;uJ-1i=0;for(i=1;i<I-1;i+)v0i=0;vJ-2i=0;/壓力場 /for(j=0;j<J;j+)for(i=0;i<I;i+)Pji=0

29、;pji=0;/doesp1=0;esp2=0;k1+;for(j=1;j<J-1;j+)ujI-2=ujI-3;for(j=1;j<J-1;j+)for(i=1;i<I-2;i+)aWji=miu*dy/dx+Max(rou*(uji-1+uji)/2*dy,0);aEji=miu*dy/dx+Max(-rou*(uji+1+uji)/2*dy,0);aNji=miu*dxe/dyn+Max(-rou*(vji+1+vji)/2*dx,0);aSji=miu*dxe/dys+Max(rou*(vj-1i+1+vj-1i)/2*dx,0);if(i=1)aWji=miu*dy

30、/dx+Max(rou*uji-1*dy,0); aNji=1.5*miu*dxe/dyn+Max(-rou*(vji*0.75+vji+1*0.25)*dx,0); aSji=1.5*miu*dxe/dys+Max(rou*(vj-1i*0.75+vj-1i+1*0.25)*dx,0); /*aNji=1.5*miu*dxe/dyn+Max(-rou*(vj0+vji+1)/2*dx,0); aSji=1.5*miu*dxe/dys+Max(rou*(vj-10+vj-1i+1)/2*dx,0);*/if(i=I-3)aEji=0;/*aEji=miu*dy/dx+Max(-rou*(uji

31、+1+uji)/2*dy,0);*/ aNji=1.5*miu*dxe/dyn+Max(-rou*(vji+1*0.75+vji*0.25)*dx,0); aSji=1.5*miu*dxe/dys+Max(rou*(vj-1i+1*0.75+vj-1i*0.25)*dx,0); /*aNji=1.5*miu*dxe/dyn+Max(-rou*(vji+vji+2)/2*dx,0); aSji=1.5*miu*dxe/dys+Max(rou*(vj-1i+vj-1i+2)/2*dx,0);*/if(j=J-2)aNji=2*miu*dxe/dyn+Max(-rou*(vji+1+vji)/2*d

32、x,0);if(j=1)aSji=2*miu*dxe/dys+Max(rou*(vj-1i+1+vj-1i)/2*dx,0);aPji=aWji+aEji+aNji+aSji;for(i=1;i<I-1;i+)for(j=1;j<J-2;j+)aW1ji=miu*dyn/dxw+Max(rou*(uj+1i-1+uji-1)/2*dy,0);aE1ji=miu*dyn/dxe+Max(-rou*(uj+1i+uji)/2*dy,0);aN1ji=miu*dx/dy+Max(-rou*(vj+1i+vji)/2*dx,0);aS1ji=miu*dx/dy+Max(rou*(vji+v

33、j-1i)/2*dx,0);if(j=1)aW1ji=1.5*miu*dyn/dxw+Max(rou*(uj+1i-1*0.25+uji-1*0.75)*dy,0); aE1ji=1.5*miu*dyn/dxe+Max(-rou*(uj+1i*0.25+uji*0.75)*dy,0); /*aW1ji=1.5*miu*dyn/dxw+Max(rou*(uj+1i-1+u0i-1)/2*dy,0); aE1ji=1.5*miu*dyn/dxe+Max(-rou*(uj+1i+u0i)/2*dy,0);*/ aS1ji=miu*dx/dy+Max(rou*vj-1i*dx,0);if(j=J-3)

34、aW1ji=1.5*miu*dyn/dxw+Max(rou*(uj+1i-1*0.75+uji-1*0.25)*dy,0); aE1ji=1.5*miu*dyn/dxe+Max(-rou*(uj+1i*0.75+uji*0.25)*dy,0);/*aW1ji=1.5*miu*dyn/dxw+Max(rou*(uj+1i-1+uj-1i-1)/2*dy,0); aE1ji=1.5*miu*dyn/dxe+Max(-rou*(uj+1i+uj-1i)/2*dy,0);*/ aN1ji=miu*dx/dy+Max(-rou*vj+1i*dx,0);if(i=1)aW1ji=2*miu*dyn/dxw

35、+Max(rou*(uj+1i-1+uji-1)/2*dy,0);if(i=I-2)/*aE1ji=2*miu*dyn/dxe+Max(-rou*(uj+1i+uji)/2*dy,0);*/ aE1ji=0;aP1ji=aW1ji+aE1ji+aN1ji+aS1ji;for(j=1;j<J-1;j+)for(i=1;i<I-2;i+)Ai=aPji/alfa1;Bi=aEji;Ci=aWji;Di=aNji*uj+1i+aSji*uj-1i+(Pji-Pji+1)*dy+(1-alfa1)*aPji/alfa1*uji;if(i=1)Ci=0;Di=aWji*uj0+aNji*uj

36、+1i+aSji*uj-1i+(Pji-Pji+1)*1.5*dy+(1-alfa1)*aPji/alfa1*uji;if(i=I-3)Bi=0;Di=aNji*uj+1i+aSji*uj-1i+(Pji-Pji+1)*1.5*dy+(1-alfa1)*aPji/alfa1*uji;TDMA(U,A,B,C,D,I-2);for(i=1;i<I-2;i+)uji=Ui;for(i=1;i<I-1;i+)for(j=1;j<J-2;j+)Aj=aP1ji/alfa1;Bj=aN1ji;Cj=aS1ji;Dj=aW1ji*vji-1+aE1ji*vji+1+(Pji-Pj+1i)

37、*dx+(1-alfa1)*aP1ji/alfa1*vji;if(j=1)Cj=0;Dj=aW1ji*vji-1+aE1ji*vji+1+aS1ji*vj-1i+(Pji-Pj+1i)*1.5*dx+(1-alfa1)*aP1ji/alfa1*vji;if(j=J-3)Bj=0;Dj=aW1ji*vji-1+aE1ji*vji+1+aN1ji*vj+1i+(Pji-Pj+1i)*1.5*dx+(1-alfa1)*aP1ji/alfa1*vji;TDMA(V,A,B,C,D,J-2);for(j=1;j<J-2;j+)vji=Vj;/壓力修正值方程/for(j=1;j<J-1;j+)

38、aPj0=1;aPjI-2=1;for(i=1;i<I-1;i+)aP10i=1;aP1J-2i=1;for(j=1;j<J-1;j+)for(i=1;i<I-1;i+)aeji=rou*dy*dy/aPji;awji=rou*dy*dy/aPji-1;anji=rou*dx*dx/aP1ji;asji=rou*dx*dx/aP1j-1i;bji=dy*rou*(uji-1-uji)+dx*rou*(vj-1i-vji);if(i=1)awji=0;if(i=I-2)aeji=0;if(j=1)asji=0;if(j=J-2)anji=0;apji=aeji+awji+anji+asji;for(k=0;k<1000;k+)for(j=1;j<J-1;j+)for(i=1;i<I-1;i+)Ai=apji;Bi=ae