1、1数值传热学 4-9 章习题答案习题 42一维稳态导热问题的控制方程： 02SxT依据本题给定条件，对节点 2 采用二阶精度的中心差分格式，节点 3 采用第三类边界条件具有二阶精度的差分格式，最后得到各节点的离散方程：节点 1： 10T节点 2： 150532节点 3： 74求解结果： ，82T3对整个控制容积作能量平衡，有： 0215)402(15)(3 xShxSqfB即：计算区域总体守恒要求满足习题 45在 42 习题中，如果 ，则各节点离散方程如下：2.03)(1fTh节点 1：节点 2： 15005321节点 3： 25.0332. )(4)( TTT对于节点 3 中的相关项作局部线
2、性化处理，然后迭代计算；求解结果： ， （迭代精度为 10-4）81.265.3迭代计算的 Matlab 程序如下：x=30;x1=20;while abs(x1-x)0.0001a=1 0 0;5 -10 5;0 -1 1+2*(x-20)(0.25);b=100;-150; 15+40*(x-20)(0.25);t=a(-1)*b;h,Tf3212x1=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 数据的基
3、数；A=cos(x);%TDMA 的主对角元素B=sin(x);%TDMA 的下对角线元素C=cos(x)+exp(x); %TDMA 的上对角线元素T=exp(x).*cos(x); %温度数据%由 A、B、C 构成 TDMAcoematrix=eye(mdim,mdim);for n=1:mdimcoematrix(n,n)=A(1,n);if n=2coematrix(n,n-1)=-1*B(1,n);endif n2b(n,:)=repmat(0,1,mdim);c(n,:)=repmat(1,1,mdim);d(n,1)=y0;elseif t10-10if counter=111y
4、val5com=TDMA(a,b,c,d,mdim);endfor nn=1:size(tt,2)for nnn=1:mdimif nnn=1d(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=2d(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=
5、mdimd(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);elsed(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);endendendyval5=TDMA(a,b,c,d,mdim);temp=yval5;yval5=yval5com;yva
6、l5com=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)=d(:,1)./a(
7、:,1);for n=2:mdimp(:,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:1y(:,n)=p(:,n).*y(:,n+1)+q(:,n);12end%-ResultCom SubFunction-function y=ResultCom (a,b,c)for n=1:max(size(c,2)y(2*n-1,:)=a(n,:);y(2*n,
8、:)=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),);elsestr=strcat(str,Pe=,num2str(d(1,n),);endendeval(eval(str);13精确解与数值解的对比图，其中边界条件给定 ， 。为了对比明显，102L给出的是 的数值解与精确解的对比：10,2P由图可以看出
9、，QUICK 和 CD 格式的计算精度较高，但两种格式都只是条件稳定；HS 和 FUS 格式绝对稳定，但 FUS 的精度较低；53乘方格式： 10,)1.0(,5PPDaeE当 时有：1.0P14951.0).1()1.0(5PDaeE因为： 3./)()(eeee uxux所以： 5297.83091.5.0eEDa由系数关系式 可得：PDewW 53.10)95.10()( weEa且： 2.0txaPp当采用隐式时 ，因此可得：1f 0597.623.15297.80 PWEPafa同理可得当 时有：0， ，0EWPa55二维稳态无源项的对流扩散问题的控制方程： )()()()( yxy
10、vxu对于一阶迎风、混合、乘方格式的通用离散方程： SNWEPaaa其中： 0,)(eeEFPADwwWa,)(nnN150,)(ssSFPADa571）QUICK 格式的界面值定义如下：)36(8WPwEe 0u对（51）式 积分可得：dxud)()(wewe dxu)()()(对流项采用 QUICK 格式的界面插值，扩散项采用线性界面插值，对于 及0u均分网格有： )()()(36()36(81 xxuu WPwPEewWPeWEP 整理得： wweEeePwewe uxxu )(81)(43)(81)(8)()(4 上式即为 QUICK 格式离散得到的离散方程；2）要分析 QUICK 格
11、式的稳定性，则应考虑非稳平流方程： xut在 时间间隔内对控制容积作积分：ttewewt dtdx得： utett )()( 随时间变化采用阶梯显式，随空间变化采用 QUICK 格式得： tux WPWEPtPt )3636(81)(整理得： xt niiniinii 87211 对于初始均匀零场，假设在 点有一个扰动 ；),(ini对 点写出 QUICK 格式的离散方程：1i16xut niiniinii 8731211 可得： ninit1对 点分析可得：1i ninixtu831由于扩散对扰动的传递恒为正，其值为 ，所以根据符号不变原则有：nit20)/83(2niinixttu整理得到
12、 QUICK 格式的稳定性条件为： 38P591）三阶迎风格式采用上游两个节点和下游一个节点的值来构造函数界面插值形式，所以定义如下： 0ucbaEPeWE根据上述定义，在 时对控制容积内的对流项作积分平均可得：0u )()(11 WPEewwe cbabxdx 由表 21 式可知三阶迎风格式的差分格式： xxniiniini 12642,由控制容积积分法得到的对流项离散格式应与 Taylor 离散展开得到的离散格式具有相同的形式和精度，所以比较可得： 6,53cba所以三阶迎风格式的函数插值定义为： 061531uEPeWE172）由上述分析可知，得到的三阶迎风格式的插值定义与给出节点上导数
13、表达式的定义在形式上显然是一致的；1861二维直角坐标中不可压缩流体的连续方程及动量方程如下： )3()()()()() )2()()()()() )1(0 vuSyxvypvxutv yxutyvxu 假设常粘性，则 ；对公式（2）及（3）分别对 求偏导得：0vuS, 3223)()()( yvxypyvyxvytvy uxutx 两式相加得并变换积分顺序有： yvxuyvxuypx uut 2222 利用连续方程有： 22ypxvuyyvxu 222xy最后即得： xvyuypx2264假设 ，则有：5*Pp10eu195.3)0(7.*nv由连续性条件有： swnevu按 SIMPLE
14、算法有： *5)(PEPee ppd 7.03nnnv将上两式代入连续性方程中有： 257.035 PPpp计算得： 06.4Pp所以： 06.4725*31.EPeu9.2)(0)(7.0Nnpv65假设 ， ，所以各点的流量为：250*3p1*61)504(1.224.)7(04.*EDCBAQ上述流量满足动量方程，但并不满足连续性方程，所以对流量修正： )(1.02.4)(.06534231pQpEDCBA对节点 3 作质量守恒有：20BDCAQ即得： )(2.04)(2.0)(1.024)(4.01 23633431 pppp 对节点 3 作质量守恒有： FEDQ即得： 20)(1.0
15、)(2.06563 pp联立求解上两式有：，7.4833.96修正后的压力为： .201.253*3 p87369166修正后的流量为： 09.4)87.40(1. 23213.)(40754.0EDCBAQ由 )(76pCQF2171Matlab 计算程序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+1for m=1:size(a,1)tjac
16、obi(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+1for m=1:size(a,1)if m=1tgs(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); elsetgs(m,n)=(-1*sum(a(m,1:m-1).*t
17、gs(1:m-1,n)-sum(a(m,m+1:end).*tgs(m+1:end,n-1)+b(m,1)/a(m,m);endendend计算结果：Jacobi Iteration初值 迭代1 迭代2 迭代3 迭代4 迭代5 迭代6 迭代7 迭代8 迭代9 迭代10T1 0 1 5 1 1 1 1 1 1 1 1T2 0 3 -3 1 1 1 1 1 1 1 1T3 0 5 -3 1 1 1 1 1 1 1 1G-S Iteration初值 迭代1 迭代2 迭代3 迭代4 迭代5 迭代6 迭代7 迭代8 迭代9 迭代10T1 0 1 -5 -23 -71 -191 -479 -1151 -2
18、687 -6143 -13823T2 0 2 9 29 81 209 513 1217 2817 6401 14337T3 0 -1 -3 -7 -15 -31 -63 -127 -255 -511 -102374常物性无内热源的稳态导热方程如下： 02yTx对上式在控制容积内积分，界面采用线性插值可得：22NSEWPTT41下边界采用补充节点法，可得到二阶精度的边界条件离散格式： xqxBjiji 1,由 可得：0,SqB 1,jijiT由上述分析可得待求四个节点的离散方程： )403(32142TT)5(4313042分别用 71 习题中的 Jacobi、G S 迭代程序求解得：Jacob
19、i iteration初值 迭代1 迭代2 迭代3 迭代4 迭代27 迭代28 迭代29 迭代30T1 0 17.5 21.875 24.913194 26.40191 28.736833 28.736836 28.736839 28.73684T2 0 12.5 17.708333 20.260417 22.032697 24.263148 24.263153 24.263155 24.263156T3 0 5 11.944444 15.347222 17.710262 20.684198 20.684203 20.684206 20.684208T4 0 3.3333333 9.16666
20、67 13.217593 15.202546 18.315777 18.315782 18.315785 18.315787G-S iteration初值 迭代1 迭代2 迭代3 迭代4 迭代16 迭代17 迭代18 迭代19T1 0 17.5 24.427083 27.270327 28.23782 28.736839 28.736841 28.736842 28.736842T2 0 16.875 21.749132 23.407691 23.972061 24.263156 24.263157 24.263158 24.263158T3 0 10.833333 17.332176 19.
21、543588 20.296082 20.684208 20.68421 20.68421 20.68421T4 0 12.569444 16.360436 17.650426 18.089381 18.315788 18.315789 18.315789 18.315789由上述计算结果可知，Jacobi 迭代的速度比 GS 迭代的速度要慢；76GS 点迭代时，各节点的离散方程如下示： )351(432TT042)(136542TGS 点迭代求解可得：23初值 迭代1 迭代2 迭代3 迭代4 迭代13 迭代14 迭代15 迭代16T1 0 12.5 24.0625 30.390625 31.9
22、72656 32.499998 32.499999 32.5 32.5T2 0 25.625 38.28125 41.445313 42.236328 42.499999 42.5 42.5 42.5T3 0 20.625 33.28125 36.445313 37.236328 37.499999 37.5 37.5 37.5T4 0 39.0625 45.390625 46.972656 47.368164 47.499999 47.5 47.5 47.5当使用 GS 线迭代时，选择自上而下的迭代，各点离散方程： )351(41(3)(2)(1 nnTT40)()()(2n5)()(1)(
23、3nn 64)(2)(3)( nTT在求解时，自上而下同时求解，即 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/4;110/4;inum=30;%the number of iterationtline=zeros(size(a,1),inum+1);temp=tline(:,1);for n=2:inum+1b1=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*tlin
24、e(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 迭代10T1 0 19.333333 30.965926 32.326703 32.49976 32.499973 32.499997 32.5 32.5T2 0 27.333333 40.885926 42.323503 42.49976 42.499973 42.499997 42.5 42.5T3 0 32.977778 36.983309 37.44
25、2021 37.49992 37.499991 37.499999 37.5 37.5T4 0 42.577778 46.967309 47.441381 47.49992 47.499991 47.499999 47.5 47.5由上述计算比较可知，线迭代的收敛迭代次数少于点迭代算法，但线迭代在每个块中采用直接求解，所以计算步骤要多于点迭代，因此两种算法的计算速度不能简单以收敛次数衡量，对于线迭代要综合考虑块间直接计算与收敛迭代次数；与例 1 相比，两者相差体现在边界条件的给定，但两者的四角温度之和相等，最终两者计算结果相同，可以解释如下：边界条件的传入是通过相关的内节点实现的，所以当某一内
26、节点相关的边界条件温度值之和相等时可以视作同一条件，因为对该内节点而言，I 类边界条件的影响效果可以线性叠加；78证明24对于题中给出的线性方程组可以用矩阵记为： bAx因为系数矩阵可以分解成： UDL其中 分别主对角阵、严格下三角阵和严格上三角阵：ULD, nnn aaaa ,3211,21321 000 0,1,23,211, nnaU 采用上述记号，则 Jacobi 迭代与 GS 迭代可以分别记作：Jacobi： )()(1)( kkk UxLbDx进一步可化为： ， 是单位阵；bAE1)()( EGS： )()1(1)( kkk xx进而可得： LDUL1)()( 由上述分析可知，Ja
27、cobi、GS 迭代的公式均属于形如： BGXkk)()1(对于某一轮引入的误差矢量 ，其迭代公式如下（令 ）：0)()1(kk对于 Jacobi 迭代， ，并由严格对角占优条件 可得：ADEG1nijjia1maxax1111 nijjinijiji对于 Jacobi 迭代，误差传递有如下关系： )()()1( kkkG所以 Jacobi 迭代过程中，当严格对角占优满足时，误差传递是衰减的；25对于 GS 迭代， ，采用类似分析可得：ULD1)( 1)max(1ax)( 11 njkjkjjnjkj所以 GS 迭代过程中，如果满足严格对角占优，误差传递也是衰减的；证毕。82在原始变量法中，连
28、续方程及动量方程如下： )3()()()()() )2()()()()() )1(0 vuSyxvypvxutv yxutyvxu 方程求解的基本变量为 ， 作为源项出现在速度 的方程中，通过假设 值可p, u, 0p以求解动量方程得到 ，通过 Poisson 方程：*vuxvyyx22可以得到计算的 值，但此时的 、 及 并不满足连续方程，不能作进一步计算。实*p*uv*p际上， 的关系隐含在连续方程中，所以在得到 后，利用连续方程得到 的修vu, *,vup正值从而获得新的 值，进行下一步计算；在涡量流函数法中只有两个变量 ， 的对流扩散方程及 控制方程如下：w,260)()(22wyx
29、ywxvu由假设的流函数 求解涡量控制方程，然后由求解获得的涡量 通过流函数控制方程更0新得到 ，重复直到收敛。因为 有相应控制方程且互相耦合，所以上述步骤可以得*,到合理的 ，最终可以得到速度 。在涡量流函数法中，不存在速度及压力关系隐含vu在连续方程的问题，所以在获得收敛的流函数后，便可以利用 Poisson 方程： 22222 yxxypx 求解 分布；p83坐标右图所示，将 分别在 点 Taylor 展开：3,2,ii),(ni（1）)(!3)(!2321,1,1,1,2, yOyyiiiii （2）)(!)()2( 41,31,21,1,3,yiiiii因为 ， ，所以 得：0,1,
30、iiu,1,1,2iiiwyu)(8)(27841,3,2, yOyiiii 整理即可得： )(8223,1,1, ywiiii84将流函数 在 点展开，有：2,i)1( )(!2)/(231,1,1,2, yOyyiiii 27因为 ， ，所以上式整理为：01,iiuy1,1,2iiiwyu)(8)(31,21,2, yOiii 进而可得 B 节点法中 Thom 公式的表达： )(21,1, yywiii将流函数 在 点展开成具有三阶截差的 Taylor 公式可得：2,i)1( )(!3)2/(!2)/(2 41,31,1,1,2, yOyyyiiiii 且有： )()/(1,2,21,3
31、yywyyiii 将流函数的各阶导数表达式代入 Taylor 公式得： )(48)()2/(8)( 431,1,21,2, yOywiiiii 整理可得 B 节点法中的 Woods 公式： )()(22,21,1, yyiiii 85对直角坐标系中，对流换热的有效压力为： cossingxypcef 对柱坐标系中，对流换热的有效压力为： cef对极坐标系中，对流换热的有效压力为： osgrpcef2891由势流的伯努利方程有： 2upe由工作条件可得空气密度 ，所以可得：205.1 Paupe 75.9843502.12由 ，可得：%5/2u .6)*.()5(222 所以： 53701;up
32、t有效压力： Patef 28.9019842两者差别： %76.5.984237pef92公式（921）中的产生项： zwzvywzuxwyzvyvyuxv xvuxutjiijijt 94由 Prandtl 混合长度理论可知，湍流的切应力公式如下： 22/mwmwlyul 当采用 模型时，有：k /2kct当脉动动量的产生与耗散相平衡时有： 2yut将上述表达式代入有：29（1）)/()/(22mwlkc再考虑 之间关系：k,mlkc2/3将上式代入式（1）中整理可得： 2/1/kcw即证。95脉动动能耗散率 在直角坐标系中展开式为：kixu 222222222 zwyxwzvyxvzyxu由 表达可知其量纲为： 32322 smtltl考查 的控制方程： kcxuxkcxxut ijjijittkk 21 右侧第一项的量纲： 414321 smgtltlm第二项除 外量纲： 1c41423/ ktlltlt第三项除 外量纲： 2 414233/ smgtltlm所以， 和 均为无量纲数；1c2由 可知右侧除 外量纲： ，与 的量纲同，所以 也是一个/ktctll2324/vc无量纲数；即证。
