數(shù)值傳熱學(xué)部分習(xí)題答案

1、精選優(yōu)質(zhì)文檔-傾情為你奉上h,Tf321習(xí)題42一維穩(wěn)態(tài)導(dǎo)熱問(wèn)題的控制方程：依據(jù)本題給定條件，對(duì)節(jié)點(diǎn)2采用二階精度的中心差分格式，節(jié)點(diǎn)3采用第三類邊界條件具有二階精度的差分格式，最后得到各節(jié)點(diǎn)的離散方程：節(jié)點(diǎn)1：節(jié)點(diǎn)2：節(jié)點(diǎn)3：求解結(jié)果：，對(duì)整個(gè)控制容積作能量平衡，有：即：計(jì)算區(qū)域總體守恒要求滿足習(xí)題45在42習(xí)題中，如果，則各節(jié)點(diǎn)離散方程如下：節(jié)點(diǎn)1：節(jié)點(diǎn)2：節(jié)點(diǎn)3：對(duì)于節(jié)點(diǎn)3中的相關(guān)項(xiàng)作局部線性化處理，然后迭代計(jì)算；求解結(jié)果：，（迭代精度為10-4）迭代計(jì)算的Matlab程序如下：x=30;x1=20;while abs(x1-x)>0.0001 a=1 0 0;5 -10 5;0

2、-1 1+2*(x-20)(0.25); b=100;-150; 15+40*(x-20)(0.25); t=a(-1)*b; x1=x; x=t(3,1);endtcal=t習(xí)題4-12的Matlab程序%代數(shù)方程形式AiTi=CiTi+1+BiTi-1+Dimdim=10;%計(jì)算的節(jié)點(diǎn)數(shù)x=linspace(1,3,mdim);%生成A、C、B、T數(shù)據(jù)的基數(shù)；A=cos(x);%TDMA的主對(duì)角元素B=sin(x);%TDMA的下對(duì)角線元素C=cos(x)+exp(x); %TDMA的上對(duì)角線元素T=exp(x).*cos(x); %溫度數(shù)據(jù)%由A、B、C構(gòu)成TDMAcoematrix=e

3、ye(mdim,mdim);for n=1:mdim coematrix(n,n)=A(1,n); if n>=2 coematrix(n,n-1)=-1*B(1,n); end if n<mdim coematrix(n,n+1)=-1*C(1,n); endend%計(jì)算D矢量D=(coematrix*T')'%由已知的A、B、C、D用TDMA方法求解T%消元P(1,1)=C(1,1)/A(1,1);Q(1,1)=D(1,1)/A(1,1);for n=2:mdim P(1,n)=C(1,n)/(A(1,n)-B(1,n)*P(1,n-1); Q(1,n)=(D(

4、1,n)+B(1,n)*Q(1,n-1)/(A(1,n)-B(1,n)*P(1,n-1);end%回迭Tcal(1,mdim)=Q(1,mdim);for n=(mdim-1):-1:1 Tcal(1,n)=P(1,n)*Tcal(1,n+1)+Q(1,n);endTcom=T;Tcal;%繪圖比較給定T值和計(jì)算T值plot(Tcal,'r*')hold onplot(T)結(jié)果比較如下，由比較可知兩者值非常切合（在小數(shù)點(diǎn)后8位之后才有區(qū)別）：節(jié)點(diǎn)1節(jié)點(diǎn)2節(jié)點(diǎn)3字段4字段5字段6字段7字段8字段9字段10T的初始值1.1.-.-2.-4.-7.-10.72954-15.03053

5、-19.T的計(jì)算值 1.1.-.-2.-4.-7.-10.72954-15.03053-19.習(xí)題4-14充分發(fā)展區(qū)的溫度控制方程如下：對(duì)于三種無(wú)量綱定義、進(jìn)行分析如下1）由得：由可得：由與無(wú)關(guān)、與無(wú)關(guān)以及、的表達(dá)式可知，除了均勻的情況外，該無(wú)量綱溫度定義在一般情況下是不能用分離變量法的；2）由得：由可得：由與無(wú)關(guān)、與無(wú)關(guān)以及、的表達(dá)式可知，在常見(jiàn)的四種邊界條件中除了軸向及周向均勻熱流的情況外，有，則該無(wú)量綱溫度定義是可以用分離變量法的；3）由得：由可得：同2）分析可知，除了軸向及周向均勻熱流的情況外，有，該無(wú)量綱溫度定義是可以用分離變量法的；習(xí)題4181）采用柱坐標(biāo)分析，寫(xiě)出統(tǒng)一的穩(wěn)態(tài)柱坐標(biāo)

6、形式動(dòng)量方程：RLq=0圖424、和分別是圓柱坐標(biāo)的3個(gè)坐標(biāo)軸，、和分別是其對(duì)應(yīng)的速度分量，其中是管內(nèi)的流動(dòng)方向；對(duì)于管內(nèi)的層流充分發(fā)展有：、，；并且方向的源項(xiàng)：方向的源項(xiàng)：方向的源項(xiàng)：由以上分析可得到圓柱坐標(biāo)下的動(dòng)量方程：方向：方向：方向：邊界條件：，；對(duì)稱線上，不考慮液體的軸向?qū)?，并?jiǎn)化分析可以得到充分發(fā)展的能量方程為：邊界條件：，；，2）定義無(wú)量綱流速：并定義無(wú)量綱半徑：；將無(wú)量綱流速和無(wú)量綱半徑代入方向的動(dòng)量方程得：上式化簡(jiǎn)得：邊界條件：，；對(duì)稱線上，定義無(wú)量綱溫度：其中，是折算到管壁表面上的平均熱流密度，即：；由無(wú)量綱溫度定義可得：將表達(dá)式和無(wú)量綱半徑代入能量方程得：化簡(jiǎn)得：（1）

7、由熱平衡條件關(guān)系可以得：將上式代入式（1）可得：邊界條件：，；，；，單值條件：由定義可知：且：即得單值性條件：3）由阻力系數(shù)及定義有：且：521一維穩(wěn)態(tài)無(wú)源項(xiàng)的對(duì)流擴(kuò)散方程如下所示：（取常物性）邊界條件如下：上述方程的精確解如下：2將分成20等份，所以有：圖示如下： 1 2 3 4 5 6 17 18 19 20 21對(duì)于中心差分、一階迎風(fēng)、混合格式和QUICK格式分別分析如下：1） 中心差分中間節(jié)點(diǎn)： 2） 一階迎風(fēng)中間節(jié)點(diǎn)：3） 混合格式當(dāng)時(shí)，中間節(jié)點(diǎn)：當(dāng)時(shí)，中間節(jié)點(diǎn)： 4） QUICK格式 數(shù)值計(jì)算結(jié)果與精確解的計(jì)算程序如下：%except for HS, any other schem

8、e doesnt take Pe<0 into consideration%expression of exact solution y=dsolve('a*b*Dy=c*D2y','y(0)=y0,y(L)=yL','x')y=subs(y,'L*a*b/c','t')y=simple(subs(y,'a*b/c*x','t*X');ysim=simple(sym(strcat('(',char(y),'-y0)','/(yL-y0)&

9、#39;)y=sym(strcat('(',char(ysim),')*(yL-y0)','+y0')% in the case of Pe=0y1=dsolve('D2y=0','y(0)=y0,y(L)=yL','x')y1=subs(y1,'-(y0-yL)/L*x','(-y0+yL)*X')%grid Pe numbertt=1 5 10;%dimensionless lengthm=20;%mdim is the number of inner nodem

10、dim=m-1;X=linspace(0,1,m+1);%initial value of variable during calculationy0=1;yL=2;%cal exact solutionfor n=1:size(tt,2) t=m*tt(1,n); if t=0 yval1(n,:)=eval(y1); else yval1(n,:)=eval(y); endend%extra treatment because max number in MATLAB is 10308if max(isnan(yval1(:) yval1=yval1' yval1=yval1(:)

11、; indexf=find(isnan(yval1); for n=1:size(indexf,1) if rem(indexf(n,1),size(X,2)=0 yval1(indexf(n),1)=yL; else yval1(indexf(n),1)=y0; end end yval1=reshape(yval1,size(X,2),size(yval1,1)/size(X,2); yval1=yval1'end%CD solutiond=zeros(size(tt,2),mdim);a=repmat(1,size(tt,2),mdim);for n=1:size(tt,2) t

12、=tt(1,n); b(n,:)=repmat(0.5*(1-0.5*t),1,mdim); c(n,:)=repmat(0.5*(1+0.5*t),1,mdim); d(n,1)=0.5*(1+0.5*tt(1,n)*y0; d(n,mdim)=0.5*(1-0.5*tt(1,n)*yL;endc(:,1)=0;b(:,mdim)=0;%numerical cal by using TDMA subfuctionyval2=TDMA(a,b,c,d,mdim);yval2=repmat(1,size(tt,2),1),yval2,repmat(2,size(tt,2),1);Fig(1,X,

13、yval1,yval2,tt);title('CD Vs. Exact Solution')% FUS solutiond=zeros(size(tt,2),mdim);a=repmat(1,size(tt,2),mdim);for n=1:size(tt,2) t=tt(1,n); b(n,:)=repmat(1/(2+t),1,mdim); c(n,:)=repmat(1+t)/(2+t),1,mdim); d(n,1)=(1+tt(1,n)/(2+tt(1,n)*y0; d(n,mdim)=1/(2+tt(1,n)*yL;endc(:,1)=0;b(:,mdim)=0;%

14、numerical cal by using TDMA subfuctionyval3=TDMA(a,b,c,d,mdim);yval3=repmat(1,size(tt,2),1),yval3,repmat(2,size(tt,2),1);Fig(2,X,yval1,yval3,tt);title('FUS Vs. Exact Solution')% HS solutiond=zeros(size(tt,2),mdim);a=repmat(1,size(tt,2),mdim);for n=1:size(tt,2) t=tt(1,n); if t>2 b(n,:)=rep

15、mat(0,1,mdim); c(n,:)=repmat(1,1,mdim); d(n,1)=y0; elseif t<-2 b(n,:)=repmat(1,1,mdim); c(n,:)=repmat(0,1,mdim); d(n,mdim)=yL; else b(n,:)=repmat(0.5*(1-0.5*t),1,mdim); c(n,:)=repmat(0.5*(1+0.5*t),1,mdim); d(n,1)=0.5*(1+0.5*t)*y0; d(n,mdim)=0.5*(1-0.5*t)*yL; endendc(:,1)=0;b(:,mdim)=0;% numerical

16、 cal by using TDMA subfuctionyval4=TDMA(a,b,c,d,mdim);yval4=repmat(1,size(tt,2),1),yval4,repmat(2,size(tt,2),1);Fig(3,X,yval1,yval4,tt);title('HS Vs. Exact Solution')%QUICK Solutiond=zeros(size(tt,2),mdim); a=repmat(1,size(tt,2),mdim);for n=1:size(tt,2) t=tt(1,n); b(n,:)=repmat(1/(2+t),1,mdi

17、m); c(n,:)=repmat(1+t)/(2+t),1,mdim); d(n,1)=(1+tt(1,n)/(2+tt(1,n)*y0; d(n,mdim)=1/(2+tt(1,n)*yL;endc(:,1)=0;b(:,mdim)=0;%numerical cal by using TDMA subfuctionyval5=zeros(size(tt,2),mdim);yval5com=yval5+1;counter=1;%iterativewhile max(max(abs(yval5-yval5com)>10-10 if counter=1 yval5com=TDMA(a,b,

18、c,d,mdim); end for nn=1:size(tt,2) for nnn=1:mdim if nnn=1 d(nn,nnn)=(6*yval5com(nn,nnn)-3*y0-3*yval5com(nn,nnn+1)*tt(1,nn)/(8*(2+tt(1,nn)+(1+tt(1,nn)/(2+tt(1,nn)*y0); elseif nnn=2 d(nn,nnn)=(5*yval5com(nn,nnn)-3*yval5com(nn,nnn+1)-yval5com(nn,nnn-1)-y0)*tt(1,nn)/(8*(2+tt(1,nn); elseif nnn=mdim d(nn

19、,nnn)=(5*yval5com(nn,nnn)-3*yL-yval5com(nn,nnn-1)-yval5com(nn,nnn-2)*tt(1,nn)/(8*(2+tt(1,nn)+(1/(2+tt(1,nn)*yL); else d(nn,nnn)=(5*yval5com(nn,nnn)-3*yval5com(nn,nnn+1)-yval5com(nn,nnn-1)-yval5com(nn,nnn-2)*tt(1,nn)/(8*(2+tt(1,nn); end end end yval5=TDMA(a,b,c,d,mdim); temp=yval5; yval5=yval5com; yv

20、al5com=temp; counter=counter+1;endyval5=yval5com;yval5=repmat(1,size(tt,2),1),yval5,repmat(2,size(tt,2),1);Fig(4,X,yval1,yval5,tt);title('QUICK Vs. Exact Solution')%-TDMA SubFunction-function y=TDMA(a,b,c,d,mdim)%form a b c d resolve yval2 by using TDMA%eliminationp(:,1)=b(:,1)./a(:,1);q(:,1

21、)=d(:,1)./a(:,1);for n=2:mdim p(:,n)=b(:,n)./(a(:,n)-c(:,n).*p(:,n-1); q(:,n)=(d(:,n)+c(:,n).*q(:,n-1)./(a(:,n)-c(:,n).*p(:,n-1);end%iterativey(:,mdim)=q(:,mdim);for n=(mdim-1):-1:1 y(:,n)=p(:,n).*y(:,n+1)+q(:,n);end%-ResultCom SubFunction-function y=ResultCom (a,b,c)for n=1:max(size(c,2) y(2*n-1,:)

22、=a(n,:); y(2*n,:)=b(n,:);end%-Fig SubFunction-function y=Fig(n,a,b,c,d)figure(n);plot(a,b);hold onplot(a,c,'*');str='''legend('for n=1:size(d,2) if n=size(d,2) str=strcat(str,'''''Pe=',num2str(d(1,n),''''')'''); else str

23、=strcat(str,'''''Pe=',num2str(d(1,n),''''','); endendeval(eval(str);精確解與數(shù)值解的對(duì)比圖，其中邊界條件給定，。為了對(duì)比明顯，給出的是的數(shù)值解與精確解的對(duì)比：由圖可以看出，QUICK和CD格式的計(jì)算精度較高，但兩種格式都只是條件穩(wěn)定；HS和FUS格式絕對(duì)穩(wěn)定，但FUS的精度較低；53乘方格式：當(dāng)時(shí)有：因?yàn)椋核裕河上禂?shù)關(guān)系式可得：且： 當(dāng)采用隱式時(shí)，因此可得：同理可得當(dāng)時(shí)有：，55二維穩(wěn)態(tài)無(wú)源項(xiàng)的對(duì)流擴(kuò)散問(wèn)題的控制方程：對(duì)于一階迎風(fēng)

24、、混合、乘方格式的通用離散方程：其中：571）QUICK格式的界面值定義如下：對(duì)（51）式積分可得：對(duì)流項(xiàng)采用QUICK格式的界面插值，擴(kuò)散項(xiàng)采用線性界面插值，對(duì)于及均分網(wǎng)格有：整理得：上式即為QUICK格式離散得到的離散方程；2）要分析QUICK格式的穩(wěn)定性，則應(yīng)考慮非穩(wěn)平流方程：在時(shí)間間隔內(nèi)對(duì)控制容積作積分：得：隨時(shí)間變化采用階梯顯式，隨空間變化采用QUICK格式得：整理得：對(duì)于初始均勻零場(chǎng)，假設(shè)在點(diǎn)有一個(gè)擾動(dòng)；對(duì)點(diǎn)寫(xiě)出QUICK格式的離散方程：可得：對(duì)點(diǎn)分析可得：由于擴(kuò)散對(duì)擾動(dòng)的傳遞恒為正，其值為，所以根據(jù)符號(hào)不變?cè)瓌t有：整理得到QUICK格式的穩(wěn)定性條件為：591）三階迎風(fēng)格式采用上游

25、兩個(gè)節(jié)點(diǎn)和下游一個(gè)節(jié)點(diǎn)的值來(lái)構(gòu)造函數(shù)界面插值形式，所以定義如下：根據(jù)上述定義，在時(shí)對(duì)控制容積內(nèi)的對(duì)流項(xiàng)作積分平均可得：由表21式可知三階迎風(fēng)格式的差分格式：由控制容積積分法得到的對(duì)流項(xiàng)離散格式應(yīng)與Taylor離散展開(kāi)得到的離散格式具有相同的形式和精度，所以比較可得：所以三階迎風(fēng)格式的函數(shù)插值定義為：2）由上述分析可知，得到的三階迎風(fēng)格式的插值定義與給出節(jié)點(diǎn)上導(dǎo)數(shù)表達(dá)式的定義在形式上顯然是一致的；61二維直角坐標(biāo)中不可壓縮流體的連續(xù)方程及動(dòng)量方程如下：假設(shè)常粘性，則；對(duì)公式（2）及（3）分別對(duì)求偏導(dǎo)得：兩式相加得并變換積分順序有：利用連續(xù)方程有：最后即得：64假設(shè)，則有：由連續(xù)性條件有：按SIM

26、PLE算法有：將上兩式代入連續(xù)性方程中有：計(jì)算得：所以： 65假設(shè)，所以各點(diǎn)的流量為：上述流量滿足動(dòng)量方程，但并不滿足連續(xù)性方程，所以對(duì)流量修正：對(duì)節(jié)點(diǎn)3作質(zhì)量守恒有：即得：對(duì)節(jié)點(diǎn)3作質(zhì)量守恒有：即得：聯(lián)立求解上兩式有：，修正后的壓力為：修正后的流量為：由71Matlab計(jì)算程序a=1 2 -2;1 1 1;2 2 1;%the CoeMatrixb=1;3;5;inum=10;%the number of iterationtjacobi=zeros(3,inum+1);tgs=zeros(3,inum+1);%jacobi iterationfor n=2:inum+1 for m=1:s

27、ize(a,1) tjacobi(m,n)=(-1*sum(a(m,:).*tjacobi(:,n-1)')+a(m,m)*tjacobi(m,n-1)+b(m,1)/a(m,m); endend%g-s iterationfor n=2:inum+1 for m=1:size(a,1) if m=1 tgs(m,n)=(-1*sum(a(m,2:end).*tgs(2:end,n-1)')+b(m,1)/a(m,m); elseif m=size(a,1) tgs(m,n)=(-1*sum(a(m,1:m-1).*tgs(1:m-1,n)')+b(m,1)/a(m,m

28、); elsetgs(m,n)=(-1*sum(a(m,1:m-1).*tgs(1:m-1,n)')-sum(a(m,m+1:end).*tgs(m+1:end,n-1)')+b(m,1)/a(m,m); end endend計(jì)算結(jié)果：Jacobi Iteration初值迭代1迭代2迭代3迭代4迭代5迭代6迭代7迭代8迭代9迭代10T101511111111T203-311111111T305-311111111G-S Iteration初值迭代1迭代2迭代3迭代4迭代5迭代6迭代7迭代8迭代9迭代10T101-5-23-71-191-479-1151-2687-6143-13

29、823T2029298120951312172817640114337T30-1-3-7-15-31-63-127-255-511-102374常物性無(wú)內(nèi)熱源的穩(wěn)態(tài)導(dǎo)熱方程如下：對(duì)上式在控制容積內(nèi)積分，界面采用線性插值可得：下邊界采用補(bǔ)充節(jié)點(diǎn)法，可得到二階精度的邊界條件離散格式：由可得：由上述分析可得待求四個(gè)節(jié)點(diǎn)的離散方程：分別用71習(xí)題中的Jacobi、GS迭代程序求解得：Jacobi iteration初值迭代1迭代2迭代3迭代4迭代27迭代28迭代29迭代30T1017.521.87524.26.4019128.28.28.28.73684T2012.517.20.22.24.24.24

30、.24.T30511.15.17.20.20.20.20.T403.9.13.15.18.18.18.18.G-S iteration初值迭代1迭代2迭代3迭代4迭代16迭代17迭代18迭代19T1017.524.27.28.2378228.28.28.28.T2016.87521.23.23.24.24.24.24.T3010.17.19.20.20.20.6842120.6842120.68421T4012.16.17.18.18.18.18.18.由上述計(jì)算結(jié)果可知，Jacobi迭代的速度比GS迭代的速度要慢；76GS點(diǎn)迭代時(shí)，各節(jié)點(diǎn)的離散方程如下示：GS點(diǎn)迭代求解可得：初值迭代1迭代2

31、迭代3迭代4迭代13迭代14迭代15迭代16T1012.524.062530.31.32.32.32.532.5T2025.62538.2812541.42.42.42.542.542.5T3020.62533.2812536.37.37.37.537.537.5T4039.062545.46.47.47.47.547.547.5當(dāng)使用GS線迭代時(shí)，選擇自上而下的迭代，各點(diǎn)離散方程：在求解時(shí)，自上而下同時(shí)求解，即1、2；3、4節(jié)點(diǎn)方程直接求解；GS線迭代求解程序：a=1 -1/4 0 0;-1/4 1 0 0;-1/4 0 1 -1/4;0 -1/4 -1/4 1;b=50/4;90/4;70

32、/4;110/4;inum=30;%the number of iterationtline=zeros(size(a,1),inum+1);temp=tline(:,1);for n=2:inum+1 b1=b+temp; tline(1:2,n)=a(1:2,1:2)b1(1:2,1); b1(3:4,1)=b1(3:4,1)+1/4*tline(1:2,n); tline(3:4,n)=a(3:4,3:4)b1(3:4,1); temp=1/4*tline(3,n);1/4*tline(4,n);0;0;end結(jié)果如下：初值迭代1迭代2迭代3迭代6迭代7迭代8迭代9迭代10T1019.3

33、0.32.32.4997632.32.32.532.5T2027.40.42.42.4997642.42.42.542.5T3032.36.37.37.4999237.37.37.537.5T4042.46.47.47.4999247.47.47.547.5由上述計(jì)算比較可知，線迭代的收斂迭代次數(shù)少于點(diǎn)迭代算法，但線迭代在每個(gè)塊中采用直接求解，所以計(jì)算步驟要多于點(diǎn)迭代，因此兩種算法的計(jì)算速度不能簡(jiǎn)單以收斂次數(shù)衡量，對(duì)于線迭代要綜合考慮塊間直接計(jì)算與收斂迭代次數(shù)；與例1相比，兩者相差體現(xiàn)在邊界條件的給定，但兩者的四角溫度之和相等，最終兩者計(jì)算結(jié)果相同，可以解釋如下：邊界條件的傳入是通過(guò)相關(guān)的內(nèi)節(jié)點(diǎn)實(shí)現(xiàn)的，所以當(dāng)某一內(nèi)節(jié)點(diǎn)相關(guān)的邊界條件溫度值之和相等時(shí)可以視作同一條件，因?yàn)閷?duì)該內(nèi)節(jié)點(diǎn)而言，I類邊界條件的影響效果可以線性疊加；78證明對(duì)于題中給出的