数值传热学部分习题答案

1、精选优质文档-倾情为你奉上h,Tf321习题42一维稳态导热问题的控制方程：依据本题给定条件，对节点2采用二阶精度的中心差分格式，节点3采用第三类边界条件具有二阶精度的差分格式，最后得到各节点的离散方程：节点1：节点2：节点3：求解结果：，对整个控制容积作能量平衡，有：即：计算区域总体守恒要求满足习题45在42习题中，如果，则各节点离散方程如下：节点1：节点2：节点3：对于节点3中的相关项作局部线性化处理，然后迭代计算；求解结果：，（迭代精度为10-4）迭代计算的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习题4-12的Matlab程序%代数方程形式AiTi=CiTi+1+BiTi-1+Dimdim=10;%计算的节点数x=linspace(1,3,mdim);%生成A、C、B、T数据的基数；A=cos(x);%TDMA的主对角元素B=sin(x);%TDMA的下对角线元素C=cos(x)+exp(x); %TDMA的上对角线元素T=exp(x).*cos(x); %温度数据%由A、B、C构成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%计算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值和计算T值plot(Tcal,'r*')hold onplot(T)结果比较如下，由比较可知两者值非常切合（在小数点后8位之后才有区别）：节点1节点2节点3字段4字段5字段6字段7字段8字段9字段10T的初始值1.1.-.-2.-4.-7.-10.72954-15.03053

5、-19.T的计算值 1.1.-.-2.-4.-7.-10.72954-15.03053-19.习题4-14充分发展区的温度控制方程如下：对于三种无量纲定义、进行分析如下1）由得：由可得：由与无关、与无关以及、的表达式可知，除了均匀的情况外，该无量纲温度定义在一般情况下是不能用分离变量法的；2）由得：由可得：由与无关、与无关以及、的表达式可知，在常见的四种边界条件中除了轴向及周向均匀热流的情况外，有，则该无量纲温度定义是可以用分离变量法的；3）由得：由可得：同2）分析可知，除了轴向及周向均匀热流的情况外，有，该无量纲温度定义是可以用分离变量法的；习题4181）采用柱坐标分析，写出统一的稳态柱坐标

6、形式动量方程：RLq=0图424、和分别是圆柱坐标的3个坐标轴，、和分别是其对应的速度分量，其中是管内的流动方向；对于管内的层流充分发展有：、，；并且方向的源项：方向的源项：方向的源项：由以上分析可得到圆柱坐标下的动量方程：方向：方向：方向：边界条件：，；对称线上，不考虑液体的轴向导热，并简化分析可以得到充分发展的能量方程为：边界条件：，；，2）定义无量纲流速：并定义无量纲半径：；将无量纲流速和无量纲半径代入方向的动量方程得：上式化简得：边界条件：，；对称线上，定义无量纲温度：其中，是折算到管壁表面上的平均热流密度，即：；由无量纲温度定义可得：将表达式和无量纲半径代入能量方程得：化简得：（1）

7、由热平衡条件关系可以得：将上式代入式（1）可得：边界条件：，；，；，单值条件：由定义可知：且：即得单值性条件：3）由阻力系数及定义有：且：521一维稳态无源项的对流扩散方程如下所示：（取常物性）边界条件如下：上述方程的精确解如下：2将分成20等份，所以有：图示如下： 1 2 3 4 5 6 17 18 19 20 21对于中心差分、一阶迎风、混合格式和QUICK格式分别分析如下：1） 中心差分中间节点： 2） 一阶迎风中间节点：3） 混合格式当时，中间节点：当时，中间节点： 4） QUICK格式 数值计算结果与精确解的计算程序如下：%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);精确解与数值解的对比图，其中边界条件给定，。为了对比明显，给出的是的数值解与精确解的对比：由图可以看出，QUICK和CD格式的计算精度较高，但两种格式都只是条件稳定；HS和FUS格式绝对稳定，但FUS的精度较低；53乘方格式：当时有：因为：所以：由系数关系式可得：且： 当采用隐式时，因此可得：同理可得当时有：，55二维稳态无源项的对流扩散问题的控制方程：对于一阶迎风

24、、混合、乘方格式的通用离散方程：其中：571）QUICK格式的界面值定义如下：对（51）式积分可得：对流项采用QUICK格式的界面插值，扩散项采用线性界面插值，对于及均分网格有：整理得：上式即为QUICK格式离散得到的离散方程；2）要分析QUICK格式的稳定性，则应考虑非稳平流方程：在时间间隔内对控制容积作积分：得：随时间变化采用阶梯显式，随空间变化采用QUICK格式得：整理得：对于初始均匀零场，假设在点有一个扰动；对点写出QUICK格式的离散方程：可得：对点分析可得：由于扩散对扰动的传递恒为正，其值为，所以根据符号不变原则有：整理得到QUICK格式的稳定性条件为：591）三阶迎风格式采用上游

25、两个节点和下游一个节点的值来构造函数界面插值形式，所以定义如下：根据上述定义，在时对控制容积内的对流项作积分平均可得：由表21式可知三阶迎风格式的差分格式：由控制容积积分法得到的对流项离散格式应与Taylor离散展开得到的离散格式具有相同的形式和精度，所以比较可得：所以三阶迎风格式的函数插值定义为：2）由上述分析可知，得到的三阶迎风格式的插值定义与给出节点上导数表达式的定义在形式上显然是一致的；61二维直角坐标中不可压缩流体的连续方程及动量方程如下：假设常粘性，则；对公式（2）及（3）分别对求偏导得：两式相加得并变换积分顺序有：利用连续方程有：最后即得：64假设，则有：由连续性条件有：按SIM

26、PLE算法有：将上两式代入连续性方程中有：计算得：所以： 65假设，所以各点的流量为：上述流量满足动量方程，但并不满足连续性方程，所以对流量修正：对节点3作质量守恒有：即得：对节点3作质量守恒有：即得：联立求解上两式有：，修正后的压力为：修正后的流量为：由71Matlab计算程序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计算结果：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常物性无内热源的稳态导热方程如下：对上式在控制容积内积分，界面采用线性插值可得：下边界采用补充节点法，可得到二阶精度的边界条件离散格式：由可得：由上述分析可得待求四个节点的离散方程：分别用71习题中的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.由上述计算结果可知，Jacobi迭代的速度比GS迭代的速度要慢；76GS点迭代时，各节点的离散方程如下示：GS点迭代求解可得：初值迭代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当使用GS线迭代时，选择自上而下的迭代，各点离散方程：在求解时，自上而下同时求解，即1、2；3、4节点方程直接求解；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结果如下：初值迭代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由上述计算比较可知，线迭代的收敛迭代次数少于点迭代算法，但线迭代在每个块中采用直接求解，所以计算步骤要多于点迭代，因此两种算法的计算速度不能简单以收敛次数衡量，对于线迭代要综合考虑块间直接计算与收敛迭代次数；与例1相比，两者相差体现在边界条件的给定，但两者的四角温度之和相等，最终两者计算结果相同，可以解释如下：边界条件的传入是通过相关的内节点实现的，所以当某一内节点相关的边界条件温度值之和相等时可以视作同一条件，因为对该内节点而言，I类边界条件的影响效果可以线性叠加；78证明对于题中给出的