数值传热学陶文铨主编第二版习题答案

数值传热学陶文铨主编第二版习题答案 数值传热学4-9章习题答案 习题4,2 一维稳态导热问题的控制方程: 2T x2,S 0 依据本题给定条件，对节点2 节点3采用第三类边界条件具有二阶精度的差分格式，最后得到各节点的离散方程:节点1: T1 100 节点2: 5T1～10T2,5T3 ～150 节点3: ～T2,4T3 75 求解结果: T2 85，T3 40 对整个控制容积作能量平衡，有: qB,S x hf(Tf～T3),S x 15 (20～40),150 2 0 即:计算区域总体守恒 要求 对教师党员的评价套管和固井爆破片与爆破装置仓库管理基本要求三甲医院都需要复审吗 满足 1 习题4,5 在4,2习题中，如果h 10 (T0.25 3～Tf)，则各节点离散方程如下: 节点1: T1 100 节点2: 5T1～10T2,5T3 ～150 节点3: ～T2,[1,2 (T3～20)0.25]T3 15,40 (T3～20)0.25 对于节点3中的相关项作局部线性化处理，然后迭代计算; 求解结果: T2 82.818，T3 35.635(迭代精度为10-4) Matlab程序 %代数方程形式AiTi=CiTi+1+BiTi-1+Di mdim=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构成TDMA coematrix=eye(mdim,mdim); for n=1:mdim coematrix(n,n)=A(1,n); if n>=2 coematrix(n,n-1)=-1*B(1,n); 2 end if n<mdim coematrix(n,n+1)=-1*C(1,n); end end %计算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(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); end Tcom=[T;Tcal]; %绘图比较给定T值和计算T值 3 plot(Tcal,?r*?) hold on plot(T) 2 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网92to.com,您的在线图书馆 结果比较如下，由比较可知两者值非常切合(在小数点后8位之后才有区别): 习题4-14 充分发展区的温度控制方程如下: cpu 对于三种无量纲定义 T1 T ( r) xr r rT～TwT～TwT～T 、 、 进行分析如下 Tb～TwTw～T T ～Tw T (Tb～Tw) ,Tw 1)由 T～Tw得: Tb～Tw 由T可得: T T T [(Tb～Tw) ,Tw] b,(1～ )w x x x x T T [(Tb～Tw) ,Tw] (Tb～Tw),(1～ )w r r r r 由Tb与r无关、 与x无关以及 T T、的 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 达式可知， 4 除了Tw均匀的情况外，该无量 x r 纲温度定义在一般情况下是不能用分离变量法的; 2)由 T～T 得: T (Tw～T ) ,T Tw～T 由T可得: T T [(Tw～T ) ,T ] w x x x T T [(Tw～T ) ,T ] (Tw～T ), w r r r r 由Tb与r无关、 与x无关以及 T T、的表达式可知， 在常见的四种边界条件中除了 x r Tw 0，则该无量纲温度定义是可以用分 r轴向及周向 均匀热流qw const的情况外，有 离变量法的; 3)由 T～Tw得: T (T ～Tw) ,Tw T ～Tw 3 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网 92to.com,您的在线图书馆 由T可得: T T [(T ～Tw) ,Tw] (1～ )w x x x T T [(T ～Tw) ,Tw] (T ～Tw),(1～w r r r 5 同2)分析可知，除了轴向及周向均匀热流qw const温度定义是可以用分离变量法的; 习题4,18 1)采用柱坐标分析，写出统一的稳态柱坐标形式动量方程: 1 1 1 1 ( u ),(r v ),( w ) ( ),( r),(),S xr rr x xr r rr r x、r和 分别是圆柱坐标的3个坐标轴，u、v和w分别是其对应的速度分量，其中x是管内的流动方向; 对于管内的层流充分发展有: u 0; x p并且x方向的源项:S ～ x p r方向的源项:S ～ r 1 p 方向的源项:S ～r v 0、w 0， 由以上分析可得到圆柱坐标下的动量方程: x方向: r方向: 方向: 边界条件: 1 u1 u p( r),()～ 0 r r rr r x p 0 r p 0 r R，u 0 u u 0;对称线上， 0 r 0， r T1 T1 ( r), xr r rr 6 T T qw;r 0，r R， r r T 0 0/ ，～ 不考虑液体的轴向导热，并简化分析可以得到充分发展的能量方程为: cpu( T) r 边界条件: 0 4 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网92to.com,您的在线图书馆 2)定义无量纲流速: U u dp～R2 dx 并定义无量纲半径: r/R;将无量纲流速和无量纲半径代入x方向的动量方程得: 1 (～R2dp1dp1 U) (～R2U) p RR ( RR ),1 R ( R )～ x 0 上式化简得: 1 ( U ),1 7 (1 U ),1 0 边界条件: 1，U 0 0， U 0;对称线上，U 0 定义无量纲温度: T～Tb q 0R/ 其中，qqw 0是折算到管壁表面上的平均热流密度，即:q0 R; 由无量纲温度定义可得: T q0R ,Tb 将T表达式和无量纲半径 代入能量方程得: c Tb1 qR pu x RR ( R0 R ),1 R ( q0R R ) 化简得: R 8 q c Tb1 1 1 pu ( ),() 0 x (1) 由热平衡条件关系可以得: c T TdTu12 Ru12q0U pu x cpub x cpumAb dx(u) q0 mA2um R2RUm 2 将上式代入式(1)可得: 5 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网 92to.com,您的在线图书馆 2U1 1 1 ( ),() Um q1 边界条件: 0， 9 0; 1， wq 0 R 0， 0; ， 0 单值条件: 由定义可知: Tb～Tbq 0 且: A UdAbb 0R/ A UdA 即得单值性条件: A UdA A UdA 0 3)由阻力系数f及Re定义有: fRe ～Ddp/dx u2mDe8 De e12 ( ) 2 um UD m 且: Nu q0DT 22 W,m～Tb (T W,m～Tb W,mqR/ 10 )0 6 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网92to.com,您的在线图书馆 5,2 1(一维稳态无源项的对流,扩散方程如下所示: 2 u 2 x x (取常物性) 边界条件如下: x 0, 0;x L, L 上述方程的精确解如下: ～ 0e(Pe x/L)～1 Pe uL/ L～ 0ePe～1 2(将L分成20等份，所以有: Pe 20P 1 2 3 4 5 6 ………… …………… 17 18 19 20 21 对于中心差分、一阶迎风、混合格式和QUICK格式分别分析如下: 1) 中心差分 中间节点: i 2) 一阶迎风 中间节点: i 3) 混合格式 11 当P 1时，中间节点: i 2, 20 i (1～0.5P ) i,1,(1,0.5P ) i～12(1～0.5P ) i,1,(1,0.5P ) i～1 i 2, 20 2 i,1,(1,P ) i～12,P i 2, 20 当P 5,10时，中间节点: i i～1 i 2, 20 4) QUICK格式 1,P P1 i i,1, i～1, 2,P 2,P 2,P 1 (5 ～ ～ ～3 )ii～1i～2i,1 8 i 2 1 (6 ～3 ～3 )ii～1i,1 i 2 8 **1,P P1 i i,1, i～1, 2,P 2,P 2,P 7 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网 92to.com,您的在线图书馆 数值计算结果与精确解的计算程序如下: %except for HS, any other scheme 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?) 12 y=simple(subs(y,?a*b/c*x?,?t*X?)); ysim=simple(sym(strcat(„(„,char(y),?-y0)?,?/(yL-y0)?))) y=sym(strcat(„(„,char(ysim),?)*(yL-y0)?,?+y0?)) % in the case of Pe=0 y1=dsolve(„D2y=0?,?y(0)=y0,y(L)=yL?,?x?) y1=subs(y1,?-(y0-yL)/L*x?,?(-y0+yL)*X?) %grid Pe number tt=[1 5 10]; %dimensionless length m=20; %mdim is the number of inner node mdim=m-1; X=linspace(0,1,m+1); %initial value of variable during calculation y0=1; yL=2; %cal exact solution for n=1:size(tt,2) t=m*tt(1,n); if t==0 yval1(n,:)=eval(y1); else 13 yval1(n,:)=eval(y); end end %extra treatment because max number in MATLAB is 104 if max(isnan(yval1(:))) yval1=yval1?; yval1=yval1(:); 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?; 8 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网92to.com,您的在线图书馆 end 14 %CD solution d=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([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; end c(:,1)=0; b(:,mdim)=0; %numerical cal by using TDMA subfuction yval2=TDMA(a,b,c,d,mdim); yval2=[repmat([1],size(tt,2),1),yval2,repmat([2],size(tt,2),1) ]; Fig(1,X,yval1,yval2,tt); title(„CD Vs. Exact Solution?) % FUS solution d=zeros(size(tt,2),mdim); a=repmat([1],size(tt,2),mdim); 15 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; end c(:,1)=0; b(:,mdim)=0; %numerical cal by using TDMA subfuction yval3=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 solution d=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 9 16 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网 92to.com,您的在线图书馆 b(n,:)=repmat([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; end end c(:,1)=0; b(:,mdim)=0; % numerical cal by using TDMA subfuction yval4=TDMA(a,b,c,d,mdim); yval4=[repmat([1],size(tt,2),1),yval4,repmat([2],size(tt,2),1) 17 ]; Fig(3,X,yval1,yval4,tt); title(„HS Vs. Exact Solution?) %QUICK Solution d=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; end c(:,1)=0; b(:,mdim)=0; %numerical cal by using TDMA subfuction yval5=zeros(size(tt,2),mdim); yval5com=yval5+1; counter=1; %iterative while max(max(abs(yval5-yval5com)))>1010 if counter==1 18 yval5com=TDMA(a,b,c,d,mdim); end for nn=1:size(tt,2) 10 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网 92to.com,您的在线图书馆 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,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); 19 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; yval5com=temp; counter=counter+1; end yval5=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) 20 %form a b c d resolve yval2 by using TDMA %elimination p(:,1)=b(:,1)./a(:,1); q(:,1)=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 %iterative y(:,mdim)=q(:,mdim); for n=(mdim-1):-1:1 11 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网 92to.com,您的在线图书馆 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,:)=a(n,:); 21 y(2*n,:)=b(n,:); end %-------------Fig SubFunction------------------ function y=Fig(n,a,b,c,d) figure(n); plot(a,b); hold on plot(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=strcat(str,?????Pe=„,num2str(d(1,n)),?????,?); end end eval(eval(str)); 12 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网92to.com,您的在线图书馆 22 精确解与数值解的对比图，其中边界条件给定 0 1， L 2。为了对比明显，给出的是P 1,2,10的数值解与精确解的对比: 由图可以看出，QUICK和CD格式的计算精度较高，但两种格式都只是条件稳定;HS和FUS格式绝对稳定，但FUS的精度较低; 5 ,3 0, 5aE (1～0.1P ),乘方格式: De (1,0.1P )5～P , ～P, P 100 P 10～10 P 0P 10 当P 0.1时有: 13 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网92to.com,您的在线图书馆 aE (1～0.1P )5 (1～0.1 0.1)5 0.951 De 因为: De e( u)e e( u)e 3/0.1 30 ( x)e( u)e( x)e(P )e 所以: 23 aE 0.951De 0.951 30 28.5297 由系数关系式aWaE～ P 可得: DwDe aW (P ,aE) Dw (0.1,0.951) 30 31.53 De 且: a0 p 当采用隐式时f 1，因此可得: P x t 1 0.1 2 0.05 0aP faE,faW,aP 28.5297,31.53,2 62.0597 同理可得当P 10时有: aE 0 ，aW 3，aP 5 5,5 二维稳态无源项的对流,扩散问题的控制方程: ( u ) ( v ) , ( ),( ) x y x x y y 对于一阶迎风、混合、乘方格式的通用离散方程: aP P aE E,aW W,aN N,aS S 其中: aE DeA(P e),～Fe,0 aW DwA(P w),Fw,0 aN DnA(P n),～Fn,0 24 14 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网 92to.com,您的在线图书馆 aS DsA(P s), Fs,0 5,7 1)QUICK格式的界面值定义如下: 1 (6 P,3 E～ W) e8 1 (6 ,3 ～ )wWPWW 8 u 0 d d( ) d( u )积分可得: 对(5,1)式 dxdx d d ( u )e～( u )w ( )e～( )w dxdx 对流项采用QUICK格式的界面插值，扩散项采用线性界 面插值，对于u 0及均分网格有: ～ W ～ P1 [(6 P,3 E～ W)( u)e～(6 W,3 P～ WW)( u)w] [ e( E)～ w(P)] 8 x x整理得: 25 333131 [( u)e～( u)w,(e,w)] P [e～( u)e] E,[w,( u)e,( u)w] W～ ( u)w WW48 x x x8 x848 上式即为QUICK格式离散得到的离散方程; 2)要分析QUICK格式的稳定性，则应考虑非稳平流方程: ～u t x 在 t时间间隔内对控制容积作积分: et, t t, te dtdx ～u w t x t w xdxdt 得: e w ( t, t～ t)dx ～u t, t t ( e～ w)dt 随时间变化采用阶梯显式，随空间变化采用QUICK格式得: 1t, tt ( P～ P) x ～u[(6 P,3 E～ W～6 W～3 P, WW)] t 26 8 整理得: in,1～ in tnn 3 in,3 in,1～7 i～1, i～2～u 8 x 对于初始均匀零场，假设在(i,n)点有一个扰动 in; 对i,1 点写出QUICK格式的离散方程: 15 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网 92to.com,您的在线图书馆 ,1n in,1～ i,1 tnnn3 in,1,3 i,2～7 i, i～1 ～u8 x 可得: ,1 in,1 7u tn i 8 x 3u tn i 8 x对i～1点分析可得: ,1 in～1 ～ 由于扩散对扰动的传递恒为正，其值为 tn i，所以根 据符号不变原则有: x2 (～3u tn tnn i, )/ ) 0 ii28 x x 整理得到QUICK格式的稳定性条件为: 27 P 8 3 5,9 1)三阶迎风格式采用上游两个节点和下游一个节点的值来构造函数界面插值形式，所以定义如下: e a E,b P,c W e a P,b E,c EEu 0u 0 根据上述定义，在u 0时对控制容积内的对流项作积分平均可得: 1e 1dx ( e～ w) x w x x 1 [a E,(b～a) P,(c～b) W～c WW] x 由表2,1式可知三阶迎风格式的差分格式: xnnn4 in,6 ～12 ,2 ,1ii～1i～2 12 xi,n 由控制容积积分法得到的对流项离散格式应与Taylor离散展开得到的离散格式具有相同的形式和精度，所以比较可得: 151a ,b ,c ～ 366 所以三阶迎风格式的函数插值定义为: 151 , ～ WeEP 366 1 ,5 ～1 ePEE E 366 u 0 u 016 28 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网92to.com,您的在线图书馆 2)由上述分析可知，得到的三阶迎风格式的插值定义与给出节点上导数表达式的定义在形式上显然是一致的; 17 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网92to.com,您的在线图书馆 6,1 二维直角坐标中不可压缩流体的连续方程及动量方程如下: u v x, y 0(1) u ( u) ( uu) p ( ) ( u y)(2) t, x,( uv) y ～ x x, x, y,Su v ( v) ( ) ( v) t , ( vu) x, ( vv) y ～ p y, x, y y,Sv(3) 假设常粘性，则Su Sv 0;对公式(2)及(3)分别对x,y求偏导得: ( u) ( uu) ( uv) p 3u 2 , , u 29 x t x x x y ～ x x , x3, x y2 ( vv) p 2v 3 v ( v) y t , y ( vu) x , y y ～ y y , y x2 , y3两式相加得并变换积分顺序有: u v u t x, y , x 2u x,vu y,u v y , y 2v v y,v u x,u v x 2 2 2 ～ p p u v 2 x2, y2 , x 2 u v x, y , y2 x, y 利用连续方程有: u u v v 30 2p 2p x u x ,v y , y v y,u x ～ x2, y2 2 u v y x, 2u 2v u v u v 2p 2p x2, y2,2 x y～2 x y ～ x2, y2 最后即得: 2p 2p u v u x2, y2 2 v x y～ y x 6,4 假设p* P 5，则有: 18 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网92to.com,您的在线图书馆 *ue 5～10 ～5 31 *vn 0.7 (5～0) 3.5 由连续性条件有: ue,vn uw,vs 按SIMPLE算法有: *???ue ue,de(pP～pE) ～5,pP *???vn vn,dn(pP～pn) 3.5,0.7pP 将上两式代入连续性方程中有: „„～5,pP,3.5,0.7pP 50,20 „pP 42.06 *?计算得: 所以: pP pP,pP 5,42.06 47.06 ue pP～pE 47.06～10 37.06 vn 0.7(pP～pN) 0.7 (47.06～0) 32.94 6,5 假设p3 250，p6 150，所以各点的流量为: * QA * QB * QC * QD Q* E** 0.4 (275～250) 10 0.2 (250～270) ～4 0.1 (10～250) ～24 0.2 (250～150) 20 0.1 (40～150) ～11 上述流量满足动量方程，但并不满足连续性方程，所以对流量修正: QA QB QC QD 32 Q E „ 10,0.4 (p1?～p3)?? ～4,0.2 (p3～p2)?? ～24,0.1 (p4～p3) „„ 20,0.2 (p3～p6)?? ～11,0.1 (p5～p6)19 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网 92to.com,您的在线图书馆 对节点3作质量守恒有: QA,QC QD,QB 即得: „„„„„„„10,0.4 (p1?～p3)～24,0.1 (p4～p3) 20,0.2 (p3～p6)～4,0.2 (p 3～p2) 对节点3作质量守恒有: QD,QE QF 即得: 20,0.2 (p?p??? 3～6)～11,0.1 (p5～p6) 20 联立求解上两式有: p?.70，p? 3 ～486 ～69.13 修正后的压力为: p p*p? 33,3 250～48.70 201.3 33 p*? 6 p6,p6 150～69.13 80.87 修正后的流量为: QA 0.4 (275～201.3) 29.48 QB 0.2 (201.3～270) ～13.74 Q 0.1 (10～201.3) ～19.13 C QD 0.2 (201.3～80.87) 24.09 QE 0.1 (40～80.87) ～4.09 由QF CF(p6～p7) 20 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网92to.com,您的在线图书馆 7,1 Matlab计算程序 a=[1 2 -2;1 1 1;2 2 1];%the CoeMatrix b=[1;3;5]; inum=10;%the number of iteration tjacobi=zeros(3,inum+1); tgs=zeros(3,inum+1); 34 %jacobi iteration for n=2:inum+1 for m=1:size(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); end end %g-s iteration for 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); else tgs(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 end end 计算结果: 35 Jacobi Iteration G-S Iteration 7,4 常物性无内热源的稳态导热方程如下: 2T 2T, 0 x2 y2 对上式在控制容积内积分，界面采用线性插值可得: 21 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网92to.com,您的在线图书馆 TP 1,TW,TE,TS,TN, 4 下边界采用补充节点法，可得到二阶精度的边界条件离散格式: Ti,j Ti,j,1, 由qB 0,S 0可得: x x SqB x , Ti,j Ti,j,1 由上述分析可得待求四个节点的离散方程: 1T1 (30,40,T2,T3) 4 1T2 (T1,30,20,T4) 4 1T3 (15,T1,T4,T3) 4 36 1T4 (T3,T2,10,T4) 4 分别用7,1习题中的Jacobi、G,S迭代程序求解得: Jacobi iteration G-S iteration 由上述计算结果可知，Jacobi迭代的速度比G,S迭代的速度要慢; 7,6 G,S点迭代时，各节点的离散方程如下示: 1(15,35,T2,T3) 4 1T2 (T1,50,40,T4) 4 1T3 (25,T1,T4,45) 4 1T4 (T3,T2,50,60) 4T1 G,S点迭代求解可得: 22 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网92to.com,您的在线图书馆 当使用G,S线迭代时，选择自上而下的迭代，各点离散方程: 1 37 (15,35,T2(n),T3(n～1)) 41 T2(n) (T1(n),50,40,T4(n～1)) 41 T3(n) (25,T1(n),T4(n),45) 41 T4(n) (T3(n),T2(n),50,60) 4T1(n) 在求解时，自上而下同时求解，即1、2;3、4节点方程直 接求解; G,S线迭代求解程序: 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/4;110/4]; inum=30;%the number of iteration tline=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 结果如下: 由上述计算比较可知，线迭代的收敛迭代次数少于点迭代 算法，但线迭代在每个块中采用直接求解，所以计算步骤要 多于点迭代，因此两种算法的计算速度不能简单以收敛次数 38 衡量，对于线迭代要综合考虑块间直接计算与收敛迭代次数; 与例1相比，两者相差体现在边界条件的给定，但两者的四角温度之和相等，最终两者计算结果相同，可以解释如下:边界条件的传入是通过相关的内节点实现的，所以当某一内节点相关的边界条件温度值之和相等时可以视作同一条件，因为对该内节点而言，I类边界条件的影响效果可以线性叠加; 7,8 证明 23 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网92to.com,您的在线图书馆 对于题中给出的线性方程组可以用矩阵记为: Ax b 因为系数矩阵可以分解成: A L,D,U 其中D,L,U分别主对角阵、严格下三角阵和严格上三角阵: 0 a a11 210 a22 39 L a31a 320 D a 33 an1an2 an,n～10 an,n 0a1,2a1,3 a1,n a 2,n U 0a2,3 0 a n～1,n 0 采用上述记号，则Jacobi迭代与G,S迭代可以分别记作: Jacobi: x(k,1) D～1(b～Lx(k)～Ux(k)) 进一步可化为: x(k,1) (E～D～1A)x(k),D～1b， E是单位阵;G,S: x(k,1) D～1(b～Lx(k～1)～Ux(k)) 进而可得: x(k,1) ～(D,L)～1Ux(k),(D,L)～1b 由上述分析可知，Jacobi、G,S迭代的公式均属于形如: X(k,1) GX(k),B 40 对于某一轮引入的误差矢量 ，其迭代公式如下(令B 0): (k,1) G (k) 对于Jacobi迭代，G E～D～1A，并由严格对角占优条件aii n aij可得: jj 1ii nn E～D～1A maxa ij max 1 1 i n ja a 1 i n j 1iii ij aiij 1 j 1i 对于Jacobi迭代，误差传递有如下关系: (k,1) (k) (k) 所以Jacobi迭代过程中，当严格对角占优满足时，误差传递是衰减的; 24 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网92to.com,您的在线图书馆 41 对于G,S迭代，G ～(D,L)U，采用类似分析可得: ～1 nakj j k,1akk～1～(D,L)U max k～1a kj1～ j 1 akk na kj) 1 max( aj 1 kk 所以G,S迭代过程中，如果满足严格对角占优，误差传 递也是衰减的; 证毕。 8,2 在原始变量法中，连续方程及动量方程如下: u v x, y 0 u u ( ) ( ) p y ( u) ( uu) ( uv),,, ～,,Su t x y x x y v v ( ) ( ) ( v) ( vu) ( vv) p y, ,, ～,,Sv t x y y x y (1)(2) (3) 0方程求解的基本变量为u,v,p，p作为源项出现在速度u,v的方程中，通过假设p值可以 求解动量方程得到u,v，通过Poisson方程: ** 2p 2p u v u v , 2 x y～ y x x2 y2 可以得到计算的p值，但此时的u、v及p并不满足连续 方程，不能作进一步计算。实际上，u,v,p的关系隐含在连续 方程中，所以在得到u,v后，利用连续方程得到p的修正值 42 从而获得新的p值，进行下一步计算; 在涡量,流函数法中只有两个变量 ,w，w的对流,扩散方程及 控制方程如下: ****** 25 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网92to.com,您的在线图书馆 ( uw) ( vw) w w , , y x x y y x 22 , ～w 0 x2 y2 由假设的流函数 求解涡量控制方程，然后由求解获得的涡量w通过流函数控制方程更新得到 ，重复直到收敛。因为 ,w有相应控制方程且互相耦合，所以上述步骤可以得到合理的 ，最终可以得到速度u,v。在涡量,流函数法中，不存在速度及压力关系隐含在连续方程的问题，所以在获得收敛的流函数后，便可以利用Poisson方程: *0 2 2 2 2 2p 2p, 2 x2 y2 ～ x y x2 y2 求解p分布; 8,3 43 坐标右图所示，将 i,2, i,3分别在(i,n)点Taylor展开: i,2 i,1, y 2 y,2 yi,1i,1( y)2 3 ,32! yi, 1( y)2,O( y3) (1) 3! (2 y)3 ,O( y4) (2) 3!i,1 i,3 i,1, y 2 (2 y),2 yi,1 2 0，2 yi,1(2 y)2 3 ,32! y 因为 y ui,1 i,1 i,1 u y wi,1，所以8 (1)～(2)得: i,1 8 i,2～ i,3 7 i,2,2( y)2wi,1,O( y4) 整理即可得: wi,1 ～7 i,1,8 i,2～ i,32 y2,O( y2) 8,4 将流函数 i,2在(i,1)点展开，有: i,2 i,1, y yi,1 2 ,2 y2( y/2)2,O( y3) 2!i,1 26 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网92to.com,您的在线图书馆 因为 u 44 y ui,10， 2 wi,1，所以上式整理为: i,1 y2 i,1 yi,1 ( y)2 i,2i,1,8wi,1,O( y3) 进而可得B节点法中Thom公式的表达: w8( i,2～ i,1) i,1 y2,O( y) 将流函数 i,2在(i,1)点展开成具有三阶截差的Taylor公式可得: y2)2 3 ( y/2)3i,2i,1, y, 2 ( y/ i,2 y212!, y3,O(i,1i,13! y4) 且有: 3 2 y3 i,1 wwi,2～wi,1 y y2 y ( y/2),O( y) 将流函数的各阶导数表达式代入Taylor公式得: ,( y)2 wi,2～wi,1 ( y)3i,2i,18wi,1, y/2) 4 (,O( y) 48,O( y) 整理可得B节点法中的Woods公式: w12( i,2～ i,1) i,1 ( y)2～12wi,2,O( y2) 45 8, 5 对直角坐标系中，对流换热的有效压力为: peff p, cgysin , cgxcos 对柱坐标系中，对流换热的有效压力为: peff p, cgx 对极坐标系中，对流换热的有效压力为: peff p, cgrcos 27 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网 92to.com,您的在线图书馆 9,1 由势流的伯努利方程有: pe p,12 u 2 由工作条件可得空气密度 1.205，所以可得: p pe～ 由u?2121 u 100000～ 1.205 502 98493.75Pa 22/u 5%，可得: u?2 (5% u)2 (0.05*50)2 6.25 46 所以: pt u 有效压力: ;2 1.205*6.25 7.53 peff p,pt 98492.75,7.53 98501.28Pa 两者差别: peff～p p 7.53 0.0076% 98492.75 9,2 公式(9,21)中的产生项: t uj uj ui , xi xi xj u u u u v u u w t u , , , , , x x x y y x z z x v v u v v v v v w w w u w w v w w w , , , y y, y , z z, y , x x, z , y y, z , z z, z x x y 9,4 由Prandtl混合长度理论可知，湍流的切应力公式如下: u2 u 2 w lm / l wm y y 当采用k～ 模型时，有: 2 47 t c k2/ 当脉动动量的产生与耗散相平衡时有: 2 u t y 将上述表达式代入有: 28 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网 92to.com,您的在线图书馆 2 (c k2/ ) ( w/ lm) (1) k3/2 再考虑 ,k之间关系: c lm 将上式代入式(1)中整理可得: /21/2 w/ c1 k 即证。 9,5 ui? 脉动动能耗散率 x k 48 ? u x 在直角坐标系中展开式为: 222222222 u? u? v? v? v? w? w? w? , y , z , x , y , z , x , y , z l2 由 表达可知其量纲为: t 考查 的控制方程: l2 l 3 t t l t , 2 49 m2 s～3 , uk t xk xk c1 ui x ,k t x j k ui uj , x j xi2 ～c2 k 1 ml2 m 右侧第一项的量纲: 34 l t lt l l t 第二项除c1外量纲: kg m～1 s～4 l2t3 l/t2 m,l/t,,l/t,m t llll t4 kg m～1 s～4 ml2/t3m 第三项除c2外量纲: 3 ll/t2l t4 所以，c1和c2均为无量纲数; 由 t c k/ 可知右侧 除c 量纲数; 即证。 2 50 ,, 2 kg m～1 s～4 4 ,l/t,外量纲: l 2 /t3 l2 ，与v的量纲同，所以c 也是一个无t 29 百度搜索“就爱阅读”,专业资料,生活学习,尽在就爱阅读网 92to.com,您的在线图书馆 51