数值传热学部分习题答案

1、习题 42 习题 42 一维稳态导热问题的控制方程： 0 2 2 =+ S x T 依据本题给定条件，对节点 2 采用二阶精度的中心差分格式， h,Tf 3 2 1 节点 3 采用第三类边界条件具有二阶精度的差分格式，最后得到各节点的离散方程： 节点 1： 100 1 =T 节点 2： 1505105 321 =+TTT 节点 3： 754 32 =+TT 求解结果： ，85 2 =T40 3 =T 对整个控制容积作能量平衡，有： 02150)4020(15)( 3 =+=+xSTThxSq ffB 即：计算区域总体守恒要求满足 习题 45习题 45 在 42 习题中，如果，则各节点离散方程如

2、下： 25. 0 3 )(10 f TTh= 节点 1： 100 1 =T 节点 2： 1505105 321 =+TTT 节点 3： 25. 0 33 25. 0 32 )20(4015)20(21 +=+TTTT 对于节点 3 中的相关项作局部线性化处理，然后迭代计算； 求解结果： ，818.82 2 =T635.35 3 =T（迭代精度为 10 -4） 迭代计算的 Matlab 程序如下： x=30; x1=20; while abs(x1-x)0.0001 a=1 0 0;5 -10 5;0 -1 1+2*(x-20)(0.25); b=100;-150; 15+40*(x-20)(0

3、.25); t=a(-1)*b; x1=x; x=t(3,1); 1 end tcal=t 习题 4-12 的 Matlab 程序 习题 4-12 的 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=ey

4、e(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 nu 对（51）式 dx dx d d dx ud )( )( =积分可得： wewe dx d dx d uu)()()()( = 对流项采用QUICK格式的界面插值，扩散项采用线性界面插值，对于及均 分网格有： 0u )()()(36()(36( 8 1 xx uu WP w PE ewWWPWeWEP =+ 整理得： WWwWwe w Ee e P we we uuu x u xxx uu)( 8 1 )(

5、 4 3 )( 8 1 )( 8 3 )()( 8 3 )( 4 3 + + = + + 上式即为QUICK格式离散得到的离散方程； 2）要分析QUICK格式的稳定性，则应考虑非稳平流方程： x u t = 在时间间隔内对控制容积作积分： t + = tt t e w e w tt t dxdt x udtdx x 得： dtudx tt t we e w ttt + + =)()( 随时间变化采用阶梯显式，随空间变化采用QUICK格式得： tux WWPWWEP t P tt P += + )3636( 8 1 )( 整理得： x u t n i n i n i n i n i n i +

6、+ + 8 733 211 1 对于初始均匀零场，假设在点有一个扰动； ),( ni n i 对点写出QUICK格式的离散方程： 1+i 15 x u t n i n i n i n i n i n i + + + + 8 733 1211 1 1 可得： n i n i x tu = + + 8 7 1 1 对点分析可得： 1i n i n i x tu = + 8 3 1 1 由于扩散对扰动的传递恒为正，其值为 n i x t 2 ，所以根据符号不变原则有： 0)/ ) 8 3 ( 2 + n i n i n i x t x tu 整理得到QUICK格式的稳定性条件为： 3 8 P 59

7、1）三阶迎风格式采用上游两个节点和下游一个节点的值来构造函数界面插值形 式，所以定义如下： += 0 0 ucba ucba EEEPe WPEe 根据上述定义，在时对控制容积内的对流项作积分平均可得： 0u )()( 1 )( 11 WWWPE e w we cbcaba x x dx xx + = = 由表21式可知三阶迎风格式的差分格式： xx n i n i n i n i ni + = + 12 21264 211 , 由控制容积积分法得到的对流项离散格式应与Taylor离散展开得到的离散格式 具有相同的形式和精度，所以比较可得： 6 1 , 6 5 , 3 1 =cba 所以三阶迎

8、风格式的函数插值定义为： += 0 6 1 6 5 3 1 0 6 1 6 5 3 1 u u EEEPe WPEe 16 2）由上述分析可知，得到的三阶迎风格式的插值定义与给出节点上导数表达式 的定义在形式上显然是一致的； 17 61 61 二维直角坐标中不可压缩流体的连续方程及动量方程如下： + + + = + + + + + = + + = + ) 3( )( )( )()()( )2( )( )( )()()( ) 1 (0 v u S y y v x x v y p y vv x vu t v S y y u x x u x p y uv x uu t u y v x u 假设常粘性

9、，则0= vu SS；对公式（2）及（3）分别对yx,求偏导得： + + = + + + + = + + 3 3 2 2 2 2 3 3 )()()( )()()( y v x v yy p yy vv yx vu yt v y y u xx u x p xy uv xx uu xt u x 两式相加得并变换积分顺序有： + + + + + = + + + + + + + y v x u yy v x u xy p x p x v u x u v y v v yy v u y u v x u u xy v x u t 2 2 2 2 2 2 2 2 22 利用连续方程有： + = + + +

10、2 2 2 2 y p x p x v u y v v yy u v x u u x + = + + + 2 2 2 2 2 2 2 2 222 y p x p y v x u y v x u y v x u x v y u 最后即得： = + x v y u y v x u y p x p 2 2 2 2 2 64 64 假设，则有： 5 * = P p 18 5105 * = e u 5 . 3)05(7 . 0 * = n v 由连续性条件有： swne vuvu+=+ 按 SIMPLE 算法有： * 5)( PEPeee pppduu+=+= * 7 . 05 . 3)( PnPnnn

11、 pppdvv+=+= 将上两式代入连续性方程中有： 20507 . 05 . 35 +=+ PP pp 计算得： 06.42 = P p 所以： 06.4706.425 * =+=+= PPP ppp 06.371006.47= EPe ppu 94.32)006.47(7 . 0)(7 . 0= NPn ppv 65 65 假设，所以各点的流量为： 250 * 3 =p150 * 6 =p = = = = = 11)15040(1 . 0 20)150250(2 . 0 24)25010(1 . 0 4)270250(2 . 0 10)250275(4 . 0 * * * * * E D

12、C B A Q Q Q Q Q 上述流量满足动量方程，但并不满足连续性方程，所以对流量修正： += += += += += )(1 . 011 )(2 . 020 )(1 . 024 )(2 . 04 )(4 . 010 6 5 6 3 3 4 2 3 3 1 ppQ ppQ ppQ ppQ ppQ E D C B A 19 对节点 3 作质量守恒有： BDCA QQQQ+=+ 即得： )(2 . 04)(2 . 020)(1 . 024)(4 . 010 2 3 6 3 3 4 3 1 pppppppp+=+ 对节点 3 作质量守恒有： FED QQQ=+ 即得： 20)(1 . 011)(

13、2 . 020 6 5 6 3 =+pppp 联立求解上两式有： 70.48 3 =p， 13.69 6 =p 修正后的压力为： 3 .20170.48250 3 * 33 =+=ppp 87.8013.69150 6 * 66 =+=ppp 修正后的流量为： = = = = = 09. 4)87.8040(1 . 0 09.24)87.803 .201(2 . 0 13.19) 3 .20110(1 . 0 74.13)2703 .201(2 . 0 48.29)3 .201275(4 . 0 E D C B A Q Q Q Q Q 由 )( 76 ppCQ FF = 20 71 71 Ma

14、tlab 计算程序计算程序 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); %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

15、: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 计算结果：计算结果： Jacobi Iter

16、ation 初值初值 迭代迭代1 1 迭代迭代2 2 迭代迭代3 3迭代迭代4 4迭代迭代5 5迭代迭代6 6迭代迭代7 7迭代迭代8 8 迭代迭代9 9 迭代迭代1010 T1 0 15 1111111 1 T2 0 3-3 1111111 1 T3 0 5-3 1111111 1 G-S Iteration 初值初值 迭代迭代1 1 迭代迭代2 2 迭代迭代3 3迭代迭代4 4迭代迭代5 5迭代迭代6 6迭代迭代7 7迭代迭代8 8 迭代迭代9 9 迭代迭代1010 T1 0 1-5 -23-71-191-479-1151-2687 -6143 -13823 T2 0 29 2981209

17、51312172817 6401 14337 T3 0 -1-3 -7-15-31-63-127-255 -511 -1023 74 74 常物性无内热源的稳态导热方程如下： 0 2 2 2 2 = + y T x T 对上式在控制容积内积分，界面采用线性插值可得： 21 () NSEWP TTTTT+= 4 1 下边界采用补充节点法，可得到二阶精度的边界条件离散格式： xqSxx TT B jiji + += +1, 由可得： 0, 0=SqB 1,+ = jiji TT 由上述分析可得待求四个节点的离散方程： )4030( 4 1 321 TTT+= )2030( 4 1 412 TTT+

18、= )15( 4 1 3413 TTTT+= )10( 4 1 4234 TTTT+= 分别用 71 习题中的 Jacobi、GS 迭代程序求解得： Jacobi iteration 初值初值 迭代迭代1 1 迭代迭代2 2 迭代迭代3 3迭代迭代4 4迭代迭代2727迭代迭代2828迭代迭代29 29 迭代迭代3030 T1 0 17.5 21.875 24.91319426.40191 28.736833 28.736836 28.736839 28.73684 T2 0 12.5 17.708333 20.260417 22.032697 24.263148 24.263153 24.2

19、63155 24.263156 T3 0 5 11.944444 15.347222 17.710262 20.684198 20.684203 20.684206 20.684208 T4 0 3.3333333 9.1666667 13.217593 15.202546 18.315777 18.315782 18.315785 18.315787 G-S iteration 初值初值 迭代迭代1 1 迭代迭代2 2 迭代迭代3 3迭代迭代4 4迭代迭代1616迭代迭代1717迭代迭代18 18 迭代迭代1919 T1 0 17.5 24.427083 27.27032728.23782

20、28.736839 28.736841 28.736842 28.736842 T2 0 16.875 21.749132 23.407691 23.972061 24.263156 24.263157 24.263158 24.263158 T3 0 10.833333 17.332176 19.543588 20.296082 20.68420820.6842120.68421 20.68421 T4 0 12.569444 16.360436 17.650426 18.089381 18.315788 18.315789 18.315789 18.315789 由上述计算结果可知，Jac

21、obi 迭代的速度比 GS 迭代的速度要慢； 76 76 GS 点迭代时，各节点的离散方程如下示： )3515( 4 1 321 TTT+= )4050( 4 1 412 TTT+= )4525( 4 1 413 +=TTT )6050( 4 1 234 +=TTT GS 点迭代求解可得： 22 初值初值 迭代迭代1 1 迭代迭代2 2 迭代迭代3 3迭代迭代4 4迭代迭代1313迭代迭代14 14 迭代迭代15 15 迭代迭代1616 T1 0 12.5 24.0625 30.39062531.97265632.49999832.499999 32.5 32.5 T2 0 25.625 38

22、.28125 41.44531342.23632842.49999942.5 42.5 42.5 T3 0 20.625 33.28125 36.44531337.23632837.49999937.5 37.5 37.5 T4 0 39.0625 45.390625 46.97265647.36816447.49999947.5 47.5 47.5 当使用 GS 线迭代时，选择自上而下的迭代，各点离散方程： )3515( 4 1 )1( 3 )( 2 )( 1 += nnn TTT )4050( 4 1 )1( 4 )( 1 )( 2 += nnn TTT )4525( 4 1 )( 4 )

23、( 1 )( 3 += nnn TTT )6050( 4 1 )( 2 )( 3 )( 4 += nnn TTT 在求解时，自上而下同时求解，即 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 iteration tline=zeros(size(a,1),inum+1); temp=tline(:,1); for n=2:inum+1 b1=b+temp; tline(1:2,n)=

24、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 1 迭代迭代2 2 迭代迭代3 3迭代迭代6 6迭代迭代7 7迭代迭代8 8 迭代迭代9 9 迭代迭代1010 T1 0 19.333333 30.965926 32.32670332.4997632.49997332.499997 32.5 32.5 T2 0 27.333333 40.8

25、85926 42.32350342.4997642.49997342.499997 42.5 42.5 T3 0 32.977778 36.983309 37.44202137.4999237.49999137.499999 37.5 37.5 T4 0 42.577778 46.967309 47.44138147.4999247.49999147.499999 47.5 47.5 由上述计算比较可知， 线迭代的收敛迭代次数少于点迭代算法， 但线迭代在每个块中采用直 接求解，所以计算步骤要多于点迭代，因此两种算法的计算速度不能简单以收敛次数衡量， 对于线迭代要综合考虑块间直接计算与收敛迭代次

26、数； 与例 1 相比，两者相差体现在边界条件的给定，但两者的四角温度之和相等，最终两者计算 结果相同，可以解释如下：边界条件的传入是通过相关的内节点实现的，所以当某一内节点 相关的边界条件温度值之和相等时可以视作同一条件，因为对该内节点而言，I 类边界条件 的影响效果可以线性叠加； 78 78 证明证明 23 对于题中给出的线性方程组可以用矩阵记为： bAx = 因为系数矩阵可以分解成： UDLA+= 其中分别主对角阵、严格下三角阵和严格上三角阵： ULD, = = nnnnnn a a a a D aaa aa a L , 33 22 11 1,21 3231 21 0 0 0 0 ? ?

27、? = 0 0 0 0 , 1 ,23 ,2 , 13 , 12, 1 nn n n a aa aaa U ? ? ? ? 采用上述记号，则 Jacobi 迭代与 GS 迭代可以分别记作： Jacobi： )( )()(1)1(kkk UxLxbDx= + 进一步可化为： ， bDxADEx kk1)(1)1( )( + +=E是单位阵； GS： )( )()1(1)1(kkk UxLxbDx= + 进而可得： bLDUxLDx kk1)(1)1( )()( + += 由上述分析可知，Jacobi、GS 迭代的公式均属于形如： BGXX kk += +)()1( 对于某一轮引入的误差矢量，其迭

28、代公式如下（令0=B） ： )()1(kk G= + 对于 Jacobi 迭代，并由严格对角占优条件ADEG 1 = = n iij j ijii aa 1 可得： 1 1 maxmax 1 1 1 1 1 = = = n ij j ij ii ni n ij j ii ij ni a aa a ADE 对于 Jacobi 迭代，误差传递有如下关系： )()()1(kkk G= + 所以 Jacobi 迭代过程中，当严格对角占优满足时，误差传递是衰减的； 24 对于 GS 迭代，采用类似分析可得： ULDG 1 )( += 1)max( 1 max)( 1 1 1 1 1 =+ = = +=

29、n j kk kj k j kk kj n kj kk kj a a a a a a ULD 所以 GS 迭代过程中，如果满足严格对角占优，误差传递也是衰减的； 证毕。 82 82 在原始变量法中，连续方程及动量方程如下： + + + = + + + + + = + + = + )3( )( )( )()()( )2( )( )( )()()( ) 1 (0 v u S y y v x x v y p y vv x vu t v S y y u x x u x p y uv x uu t u y v x u 方程求解的基本变量为，作为源项出现在速度的方程中，通过假设值可以 求解动量方程得到，通过 Poisson 方程： pvu,pvu, 0 p *,v u = + x v y u y v x u y p x p 2 2 2 2 2 可以得到计算的值，但此时的、及并不满足连续方程，不能作进一步计算。实 际上，的关系隐含在连续方程中，所以在得到后，利用连续方程得到的修正 值从而获得新的值，进行下一步计算； * p * u * v * p pvu, *,v up p 在涡量流函数法中只有两个变量w,，的对流扩散