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

数值传热学4-9章习题答案
习题4－2
一维稳态导热问题的控制方程：
022=+∂∂S x
T
λ依据本题给定条件，对节点2
节点3采用第三类边界条件具有二阶精度的差分格式，最后得到各节点的离散方程：
节点1：100
1=T 节点2：1505105321-=+-T T T 节点3：75
432=+-T T 求解结果：，852=T 40
3=T 对整个控制容积作能量平衡，有：
2150)4020(15)(3=⨯+-⨯=∆+-=∆+x S T T h x S q f f B 即：计算区域总体守恒要求满足
习题4－5
在4－2习题中，如果，则各节点离散方程如下：
25
.03)
(10f T T h -⨯=节点1：100
1=T 节点2：150
5105321-=+-T T T 节点3：25
.03325.032)20(4015])20(21[-⨯+=-⨯++-T T T T 对于节点3中的相关项作局部线性化处理，然后迭代计算；
求解结果：
，（迭代精度为10-4）
818.822=T 635.353=T 迭代计算的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.25)]; t=a^(-1)*b;
x1=x;
x=t(3,1);
end
tcal=t
习题4-12的Matlab程序
%代数方程形式A i T i=C i T i+1+B i T i-1+D i
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);
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值
plot(Tcal,'r*')
hold on
plot(T)
n g
i
n t
h a r e 结果比较如下，由比较可知两者值非常切合（在小数点后8位之后才有区别）：
习题4-14
充分发展区的温度控制方程如下：
)(1r
T
r r r x T u
c p ∂∂∂∂=∂∂λρ对于三种无量纲定义、、进行分析如下
w b w T T T T --=
Θ∞∞
--=ΘT T T T w w
w T T T T --=Θ∞1）由得：
w
b w
T T T T --=
Θw
w b T T T T +Θ-=)(由可得：
T x T x T x T T T x T w b w w b ∂∂Θ-+∂∂Θ=∂+Θ-∂=∂∂)1(]
)[(r
T r T T r T T T r T w w b w w b ∂∂Θ-+∂Θ∂-=∂+Θ-∂=∂∂)1()(])[(由与无关、与无关以及
、的表达式可知，除了均匀的情况外，该无量b T r Θx x T ∂∂r
T
∂∂w T 纲温度定义在一般情况下是不能用分离变量法的；2）由得：
∞
∞
--=
ΘT T T T w ∞
∞+Θ-=T T T T w )(由可得：
T x
T x T T T x T w w ∂∂Θ=∂+Θ-∂=∂∂∞∞])[(r
T r T T r T T T r T w w w ∂∂Θ+∂Θ∂-=∂+Θ-∂=∂∂∞∞∞)(])[(由与无关、与无关以及
、的表达式可知，在常见的四种边界条件中除了b T r Θx x T ∂∂r
T ∂∂轴向及周向均匀热流的情况外，有，则该无量纲温度定义是可以用分const q w =0=∂∂r
T w
离变量法的；3）由得：
w
w
T T T T --=
Θ∞w
w T T T T +Θ-=∞)(
由可得：
T x
T x T T T x T w w w ∂∂Θ-=∂+Θ-∂=∂∂∞)1(]
)[(
r T T r T T T r T w w w -+∂Θ∂-=∂+Θ-∂=∂∂∞∞1()(])[(同2）分析可知，除了轴向及周向均匀热流const q w =温度定义是可以用分离变量法的；
习题4－18
1）采用柱坐标分析，写出统一的稳态柱坐标形式动量方程：
S r r r r r r x x w r v r r r u x +∂∂∂∂+∂∂∂∂+∂∂∂∂=∂∂+∂∂+∂∂(1)(1)()(1)(1)(θφλθφλφλφρθφρφρ、和分别是圆柱坐标的3个坐标轴，、和分别是其对应的速度分量，其中x r θu v w 是管内的流动方向；
x 对于管内的层流充分发展有：
、，
；0=v 0=w 0=∂∂x
u
并且方向的源项：x x p
S ∂∂-
=方向的源项：r r p
S ∂∂-
=方向的源项：θθ∂∂-
=p
r S 1由以上分析可得到圆柱坐标下的动量方程：
方向：
x 0)(1)(1=∂∂-∂∂∂∂+∂∂∂∂x p
u r r r u r r r θλθλ方向：
r 0=∂∂r p 方向：
θ0=∂∂θ
p 边界条件：
，R r =0
=u ，；对称线上，0=r 0=∂∂r u 0=∂∂θu 不考虑液体的轴向导热，并简化分析可以得到充分发展的能量方程为：
)(
1(1θ
λθλρ∂∂∂∂+∂∂∂∂=∂∂T
r r r T r r r x T u
c p 边界条件：
，；，R r =w q r T =∂∂λ0=r 0
=∂∂r
T
，πθ/0=0=∂∂-θ
λT
2）定义无量纲流速：
dx
dp R u
U 2
-=
λ并定义无量纲半径：；将无量纲流速和无量纲半径代入方向的动量方程得：
R r /=ηx 0))1
((1)1((122
=∂∂-∂-∂∂∂+∂-∂∂∂
x
p U dx dp R R R R U dx dp R R
R R θληλθηηλληη
η上式化简得：
011(1(1=+∂∂∂∂+∂∂∂∂θ
ηθηηηηηU U 边界条件：
，1=η0
=U ，
；对称线上，0=η0=∂∂ηU 0=∂∂θ
U
定义无量纲温度：
λ
/0R q T T b
-=Θ其中，是折算到管壁表面上的平均热流密度，即：；0q R
q q w
π=
0由无量纲温度定义可得：
b
T R
q T +Θ=
λ
0将表达式和无量纲半径代入能量方程得：
T η(1)(100θ
ληλθηηλληηηρ∂Θ
∂∂∂+∂Θ∂∂∂=∂∂R q R R R R q R R R x T u
c b p 化简得：
（1）
)1(1)(10θ
ηθηηηηηρ∂Θ∂∂∂+∂Θ∂∂∂=∂∂x T u c q R b p 由热平衡条件关系可以得：
m
m m b m p b p p RU U q R u u R q A u u dx dT A u c x T u c x T u
c 02022
1221)(===∂∂=∂∂ππρρρ将上式代入式（1）可得：
)1(1)(12θ
ηθηηηηη∂Θ
∂∂∂+∂Θ∂∂∂=m U U 边界条件：，；，
0=η0=∂Θ
∂η
1=ηR q q w πη10==∂Θ∂，
；，0=θ0=∂Θ∂θπθ=0=∂Θ
∂θ
单值条件：由定义可知： 且： 0/0=-=ΘλR q T T b b b ⎰⎰Θ=
ΘA
A
b UdA
UdA 即得单值性条件：
=Θ⎰⎰A
A UdA
UdA 3）由阻力系数及定义有：
f Re 228)(21/Re ⎪⎭⎫ ⎝⎛=⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡-=D D U D u u dx dp D f e m e m m
e νρ且：m W b m W b m W R q T T D T T q Nu ,0,,0~2
)/(2Θ=-=-=
λ
λ
5－2
1．一维稳态无源项的对流－扩散方程如下所示：
（取常物性）
x
x u 2
2∂∂Γ=∂∂φφρ边界条件如下：
L
L x x φφφφ====,;,00上述方程的精确解如下：
1
1
)/(00--=
--⋅Pe L x Pe L e e φφφφΓ=/uL Pe ρ2．将分成20等份，所以有：
L ∆
=P Pe 20 1 2 3 4 5 6
…………
…………… 17 18 19 20 21
对于中心差分、一阶迎风、混合格式和QUICK 格式分别分析如下：1）中心差分
中间节点：
2
)5.01()5.01(1
1-∆+∆++-=
i i i P P φφφ20,2 =i 2）一阶迎风
中间节点： ∆
-∆++++=P P i i i 2)1(1
1φφφ20
,2 =i 3）混合格式
当时，中间节点：
1=∆P 2
)5.01()5.01(1
1-∆+∆++-=
i i i P P φφφ 20
,2 =i 当时，中间节点： 10,5=∆P 1-=i i φφ20,2 =i 4）QUICK 格式
*
12111)35(8122121
⎥⎦
⎤⎢⎣⎡---++++++=+--∆
∆
-∆∆+∆i i i i i i i P P P P P φφφφφφφ2
≠i
*
1111)336(8122121
⎥⎦
⎤
⎢⎣⎡--++++++=+-∆∆
-∆∆+∆i i i i i i P P P P P φφφφφφ2
=i 数值计算结果与精确解的计算程序如下：
%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')
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 yval1(n,:)=eval(y); end end
%extra treatment because max number in MATLAB is 10^308if 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';
end
%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);
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
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)]; 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)))>10^-10
if counter==1
yval5com=TDMA(a,b,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,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;
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)
%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
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,:);
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));
a n d A l l t h i n g s i n t h e
i r b e i n g a r e g 13
精确解与数值解的对比图，其中边界条件给定，。

为了对比明显，10=φ2=L φ给出的是的数值解与精确解的对比：
10,2,1=∆P 由图可以看出，QUICK 和CD 格式的计算精度较高，但两种格式都只是条件稳定；HS 和FUS 格式绝对稳定，但FUS 的精度较低；
5－3
乘方格式：⎪⎪⎩⎪⎪⎨
⎧<-≤≤--+≤≤->=∆∆∆∆∆∆∆∆10
,010,
)1.01(100,
)1.01(10,055
P P P P P P P P D a e E 当时有：
1.0=∆P
14
951.0)1.01.01()1.01(55=⨯-=-=∆P D a e
E
因为：301.0/3)()()()()()(===Γ=Γ=
∆e
e
e e e e e e e P u x u u x D ρδρρδ所以：5297
.2830951.0951.0=⨯==e E D a 由系数关系式可得：∆=-P D a D a e
E
w W 53.3130)951.01.0((=⨯+=⨯+
=∆w e
E
W D D a P a 且：
205
.01
.010=⨯=
∆∆=
t
x
a P p ρ当采用隐式时，因此可得：
1=f 0597
.62253.315297.280
=++=++=P W E P a fa fa a 同理可得当时有：
10=∆P ，，0
=E a 3=W a 5
=P a 5－5
二维稳态无源项的对流－扩散问题的控制方程：
)()()()(y
y x x y v x u ∂∂Γ∂∂+∂∂Γ∂∂=∂∂+∂∂φ
φφρφρφφ对于一阶迎风、混合、乘方格式的通用离散方程：
S
S N N W W E E P P a a a a a φφφφφ+++=其中：
[]0,)(e e e E F P A D a -+=∆[]0,)(w w w W F P A D a +=∆[]
0,)(n n n N F P A D a -+=∆
[]
0,)(s s s S F P A D a +=∆5－7
1）QUICK 格式的界面值定义如下：
⎪⎪⎩
⎪⎪⎨
⎧-+=-+=)36(81)36(81WW P W w W E P e φφφφφφφφ0
>u 对（5－1）式
积分可得：dx
dx d d dx u d )
()(φ
φρΓ
=w
e w e dx
d dx d u u )()()()(φ
φφρφρΓ-Γ=-对流项采用QUICK 格式的界面插值，扩散项采用线性界面插值，对于及0>u 均分网格有：
)](([]))(36())(36[(81
x
x u u W P w P E e w WW P W e W E P ∆-Γ-∆-Γ=-+--+φφφφρφφφρφφφ整理得：
WW w W w e w E e e P w e w e u u u x u x x x u u φρφρρφρφρρ)(8
1
])(43)(81[])(83[)]()(83)(43[-++∆Γ+-∆Γ=∆Γ+∆Γ+-上式即为QUICK 格式离散得到的离散方程；
2）要分析QUICK 格式的稳定性，则应考虑非稳平流方程：
x
u t ∂∂-=∂∂φφ在时间间隔内对控制容积作积分：
t ∆⎰⎰⎰⎰∆+∆+∂∂-=∂∂t t t e w e w t t t dxdt
x u dtdx x φ
φ得： dt
u dx t
t t
w e e
w
t t t ⎰
⎰∆+∆+--=-)()(φφφφ随时间变化采用阶梯显式，随空间变化采用QUICK 格式得：
φt
u x WW P W W E P t
P t t P ∆+---+-=∆-∆+)]3636(8
1[)(φφφφφφφφ整理得：
x
u
t
n
i n i n i n i n
i n i ∆+-+-∆---++87332111φφφφφφ对于初始均匀零场，假设在点有一个扰动；),(n i n i ε对点写出QUICK 格式的离散方程：
1+i
x
u
t
n
i n i n i n i n i n i ∆+-+-∆--+++++87331211
11φφφφφφ可得：
n
i
n i x t u εφ∆∆=
++8711对点分析可得：
1-i n
i
n i x
t u εφ∆∆-
=+-8311由于扩散对扰动的传递恒为正，其值为
，所以根据符号不变原则有：n
i x t ερ2
∆Γ∆0)/)83(2
≥∆Γ∆+∆∆-
n
i n i n i x
t x t u εερε整理得到QUICK 格式的稳定性条件为：
3
8≤
∆P 5－9
1）三阶迎风格式采用上游两个节点和下游一个节点的值来构造函数界面插值形式，所以定义如下：
⎩⎨
⎧<++=>++=0
0u c b a u c b a EE
E P e W P E e φφφφφφφφ根据上述定义，在时对控制容积内的对流项作积分平均可得：
0>u ])()([1)
(11WW W P E e w w e c b c a b a x
x
dx x x φφφφφφφ--+-+∆=-∆=∂∂∆⎰由表2－1式可知三阶迎风格式的差分格式：
x
x
n
i n i n i n i n
i ∆+-+=
∂∂--+1221264211,φφφφφ由控制容积积分法得到的对流项离散格式应与Taylor 离散展开得到的离散格式具有相同的形式和精度，所以比较可得：
6
1,65,31-
===c b a 所以三阶迎风格式的函数插值定义为：
⎪⎪⎩
⎪⎪⎨
⎧<-+=>-+=0
6165310616531u u EE E P e W P E e φφφφφφφφ
2）由上述分析可知，得到的三阶迎风格式的插值定义与给出节点上导数表达式的定义在形式上显然是一致的；
6－1
二维直角坐标中不可压缩流体的连续方程及动量方程如下：
⎪⎪⎪⎪⎪⎩
⎪⎪⎪⎪⎪⎨
⎧+∂∂∂∂+∂∂∂∂+∂∂-=∂∂+∂∂+∂∂+∂∂∂∂+∂∂∂∂+∂∂-=∂∂+∂∂+∂∂=∂∂+∂∂)
3()
()()()()()
2()()()()()()
1(0v
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 ηηρρρηηρρρ假设常粘性，则；对公式（2）及（3）分别对求偏导得：
0==v u S S y x ,⎪⎪⎩⎪⎪⎨
⎧∂∂+⎪⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂∂∂-=⎪⎪⎭
⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂⎪⎪⎭
⎫
⎝⎛∂∂∂∂+∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂-=⎪⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂33222233)()()()()()(y v x v y y p y y vv y x vu y t v y y u x x u x p x y uv x x uu x t u x ηηρρρηηρρρ两式相加得并变换积分顺序有：
⎥⎦⎤⎢⎣
⎡⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂+⎥⎦⎤⎢⎣⎡⎪⎪⎭⎫
⎝⎛∂∂+∂∂∂∂+⎪⎪⎭⎫
⎝⎛∂∂+∂∂-=⎥⎦⎤
⎢⎣
⎡⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂y v x u y y v x u x y p x p x v u x u v y v v y y v u y u v x u u x y v x u t 2222222222ηηρρ利用连续方程有：
⎪⎪⎭⎫ ⎝⎛∂∂+∂∂-=⎥⎦⎤
⎢⎣⎡⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂2222y p x p x v u y v v y y u v x
u u x ρ⎪⎪⎭
⎫ ⎝⎛∂∂+∂∂-=⎦⎤
⎢⎣⎡∂∂∂∂-∂∂∂∂+∂∂+∂∂+∂∂∂∂22222222222y 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 ρ222226－4
假设，则有：
5*
=P p 5
105*
-=-=e u
5
.3)05(7.0*=-⨯=n v 由连续性条件有：
s
w n e v u v u +=+按SIMPLE 算法有：
'
''*5)(P E P e e e p p p d u u +-=-+='''*7.05.3)(P
n P n n n p p p d v v +=-+=将上两式代入连续性方程中有：
20
507.05.35''+=+++-P P p p 计算得：
06.42'=P p 所以：
06
.4706.425'
*
=+=+=P P P p p p 06
.371006.47=-=-=E P e p p u 94
.32)006.47(7.0)(7.0=-⨯=-=N P n p p v 6－5
假设，，所以各点的流量为：
250*
3=p 150*
6=p ⎪⎪⎪⎩
⎪⎪⎪⎨⎧-=-⨯==-⨯=-=-⨯=-=-⨯==-⨯=11
)15040(1.020)150250(2.024)25010(1.04)270250(2.010)250275(4.0*****E D
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'1p p Q p p Q p p Q p p Q p p Q E
D C B A
对节点3
作质量守恒有：
B
D C A Q Q Q Q +=+即得：
)
(2.04)(2.020)(1.024)(4.010'
2'3'6'3'3'4'3'1p p p p p p p p -⨯+--⨯+=-⨯+--⨯+对节点3作质量守恒有：
F
E D Q Q Q =+即得：
20
)(1.011)(2.020'
6'5'6'3=-⨯+--⨯+p p p p 联立求解上两式有：
，70.48'3-=p 13
.69'
6-=p 修正后的压力为：
3.20170.48250'3*33=-=+=p p p 87
.8013.69150'6*66=-=+=p p p 修正后的流量为：
⎪⎪⎪⎩⎪⎪⎪⎨⎧-=-⨯==-⨯=-=-⨯=-=-⨯==-⨯=09
.4)87.8040(1.009.24)87.803.201(2.013.19)3.20110(1.074.13)2703.201(2.048.29)3.201275(4.0E
D C B A Q Q Q Q Q 由)
(76p p C Q F F -=
21
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);%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
计算结果：Jacobi Iteration
G-S Iteration
7－4
常物性无内热源的稳态导热方程如下：
02
222=∂∂+∂∂y T
x T 对上式在控制容积内积分，界面采用线性插值可得：
22
()N S E W P T T T T T +++=
4
1
下边界采用补充节点法，可得到二阶精度的边界条件离散格式：
λ
δλδx
q S x x T T B j i j i ⋅+⋅∆⋅+
=+1,,由可得：
0,0==S q B 1
,,+=j i j i T T 由上述分析可得待求四个节点的离散方程：
)
4030(41
321T T T +++=)
2030(41
412T T T +++=)
15(41
3413T T T T +++=)
10(4
1
4234T T T T +++=分别用7－1习题中的Jacobi 、G －S 迭代程序求解得：
Jacobi iteration
G-S iteration
由上述计算结果可知，Jacobi 迭代的速度比G －S 迭代的速度要慢；
7－6
G －S 点迭代时，各节点的离散方程如下示：
)3515(41
321T T T +++=
)
4050(41
412T T T +++=)
4525(41
413+++=T T T )
6050(4
1
234+++=T T T G －S 点迭代求解可得：
23
当使用G －S 线迭代时，选择自上而下的迭代，各点离散方程：
)3515(41
)1(3)(2)(1-+++=
n n n T T T )
4050(41
)1(4)(1)(2-+++=n n n T T T )
4525(41
)(4)(1)(3+++=n n n T T T )
6050(41
)(2)(3)(4+++=n n n T T T 在求解时，自上而下同时求解，即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
结果如下：
由上述计算比较可知，线迭代的收敛迭代次数少于点迭代算法，但线迭代在每个块中采用直接求解，所以计算步骤要多于点迭代，因此两种算法的计算速度不能简单以收敛次数衡量，对于线迭代要综合考虑块间直接计算与收敛迭代次数；
与例1相比，两者相差体现在边界条件的给定，但两者的四角温度之和相等，最终两者计算结果相同，可以解释如下：边界条件的传入是通过相关的内节点实现的，所以当某一内节点相关的边界条件温度值之和相等时可以视作同一条件，因为对该内节点而言，I 类边界条件的影响效果可以线性叠加；
7－8
证明
对于题中给出的线性方程组可以用矩阵记为：
b
Ax =因为系数矩阵可以分解成：
U
D L A ++=其中分别主对角阵、严格下三角阵和严格上三角阵：
U L D ,,⎥⎥
⎥
⎥⎥⎥⎦
⎤⎢⎢⎢⎢
⎢⎢⎣
⎡=⎥⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢
⎢⎢⎣
⎡=-n n n n n n a a a a D a a a a a a L ,33
22
11
1
,2
1
3231210000
⎥⎥
⎥
⎥
⎥⎥⎦
⎤
⎢⎢⎢⎢
⎢⎢⎣⎡=-0000,1,23,2,13,12
,1n n n n a a a a a a U 采用上述记号，则Jacobi 迭代与G －S 迭代可以分别记作：Jacobi ： )
()()(1)1(k k k Ux Lx b D x
--=-+进一步可化为： ，　是单位阵；
b D x A D E x k k 1)(1)
1()(--++-=E G －S ： )()()1(1)1(k k k Ux Lx b D x --=--+进而可得： b
L D Ux L D x
k k 1)(1)1()()(--++++-=由上述分析可知，Jacobi 、G －S 迭代的公式均属于形如：
B
GX X k k +=+)()1(对于某一轮引入的误差矢量，其迭代公式如下（令）：
ε0=B )
()1(k k G εε=+对于Jacobi 迭代，，并由严格对角占优条件可得：
A D E G 1
--=∑≠=>
n
ii j j ij
ii a
a 1
11max max 11111<⎪⎪⎪⎭
⎫ ⎝
⎛==-∑∑
≠=≤≤≠=≤≤-n
i j j ij ii
n i n
i
j j ii ij
n
i a a a a A D E 对于Jacobi 迭代，误差传递有如下关系：
)
()()1(k k k G εεε<=+所以Jacobi 迭代过程中，当严格对角占优满足时，误差传递是衰减的；
对于G －S 迭代，，采用类似分析可得：
U L D G 1
)(-+-=1)max(1max )(11111
<<⎪⎪⎪⎪⎪
⎭
⎫
⎝⎛-=+-∑∑∑=-=+=-n j kk
kj
k j kk kj n k j kk kj a a a a a a U L D 所以G －S 迭代过程中，如果满足严格对角占优，误差传递也是衰减的；证毕。

8－2
在原始变量法中，连续方程及动量方程如下：
⎪⎪⎪⎪⎪⎩
⎪⎪⎪⎪⎪⎨⎧+∂∂∂∂+∂∂∂∂+∂∂-=∂∂+∂∂+∂∂+∂∂∂∂+∂∂∂∂+∂∂-=∂∂+∂∂+∂∂=∂∂+∂∂)
3()()()()()()
2(()()()()()
1(0v 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 ηηρρρηηρρρ方程求解的基本变量为，作为源项出现在速度的方程中，通过假设值可
p v u ,,p v u ,0
p 以求解动量方程得到，通过Poisson 方程：
*
*
,v u ⎦⎤⎢⎣⎡∂∂∂∂-∂∂∂∂=⎪
⎪⎭
⎫ ⎝⎛∂∂+∂∂x v y u y v x u y p x p ρ22222可以得到计算的值，但此时的、及并不满足连续方程，不能作进一步计算。

实*
p *
u *
v *
p 际上，的关系隐含在连续方程中，所以在得到后，利用连续方程得到的修p v u ,,*
*
,v u p 正值从而获得新的值，进行下一步计算；
p 在涡量－流函数法中只有两个变量，的对流－扩散方程及控制方程如下：
w ,ψw ψ
26
⎪⎪⎩
⎪⎪⎨⎧=-∂∂+∂∂⎪⎪⎭
⎫
⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂=∂∂+∂∂0)()(22
22
w y x y w y x w x y vw x uw ψψηηρρ由假设的流函数求解涡量控制方程，然后由求解获得的涡量通过流函数控制方程更0
ψw 新得到，重复直到收敛。

因为有相应控制方程且互相耦合，所以上述步骤可以得*ψw ,ψ到合理的，最终可以得到速度。

在涡量－流函数法中，不存在速度及压力关系隐含ψv u ,在连续方程的问题，所以在获得收敛的流函数后，便可以利用Poisson 方程：
⎥⎥⎦
⎤
⎢⎢⎣⎡⎪⎪⎭
⎫ ⎝⎛∂∂∂-⎪⎪⎭⎫ ⎝⎛∂∂⎪⎪⎭⎫ ⎝
⎛∂∂=∂∂+∂∂2
2222
222222y x y x y p
x p ψψ
ψρ求解分布；
p 8－3
坐标右图所示，将分别在点Taylor 展开：
3,2,,i i ψψ),(n i （1）)(!3)(!2)(32
1
,3321,2
21
,1,2
,y O y y y y y y
i i i i i δδψ
δψδψ
ψψ+∂∂+∂∂+∂∂+
=　（2）)(!3)2(!2)2()2(43
1
,3
321
,2
21
,1,3
,y O y y y y y y
i i i i i δδψ
δψδψ
ψψ+∂∂+∂∂+∂∂+
=因为
，，所以得：
01
,1
,==∂∂i i u y
ψ1,1
,1
,2
2i i i w y
u
y
=∂∂=
∂∂ψ
)2()1(8-⨯)
()(27841,22,3,2,y O w y i i i i δδψψψ++=-整理即可得：
)
(28722
3
,2,1,1,y O y
w i i i i δδψψψ+-+-=
8－4
将流函数在点展开，有：
2,i ψ)1,(i )(!2)2/(232
1
,221
,1,2
,y O y y
y
y
i i i i δδψδψ
ψψ+∂∂+∂∂+
=
27
因为
，，所以上式整理为：
01
,1
,==∂∂i i u y
ψ1,1
,1
,2
2i i i w y
u y
=∂∂=
∂∂ψ)
(8
)(31,2
1,2
,y O w y i i i δδψψ++=进而可得B 节点法中Thom 公式的表达：
)
()
(82
1,2,1,y O y
w i i i δδψψ+-=
将流函数在点展开成具有三阶截差的Taylor 公式可得：
2,i ψ)1,(i )(!3)2/(!2)2/(243
1
,3
321,221
,1,2
,y O y y y y
y
y
i i i i i δδψ
δψδψ
ψψ+∂∂+∂∂+∂∂+
=且有：
)()2/(1
,2,2
21
,3
3y O y w w y w y y y i i i δδψ
ψ+-=∂∂=⎪⎪⎭
⎫ ⎝⎛∂∂∂∂=∂∂将流函数的各阶导数表达式代入Taylor 公式得：
)(48
)()()2/(8)(4
31,2,1,21,2
,y O y y O y w w w y i i i i i δδδδδψψ+⎪⎪⎭⎫ ⎝⎛+-++=整理可得B 节点法中的Woods 公式：
)(2
1
)()
(1222,2
1,2,1,y O w y w i i i i δδψψ+--=
8－5
对直角坐标系中，对流换热的有效压力为：
θ
ρθρcos sin gx gy p p c c eff ++=对柱坐标系中，对流换热的有效压力为：
gx
p p c eff ρ+=对极坐标系中，对流换热的有效压力为：
θ
ρcos gr p p c eff +=
9－1
由势流的伯努利方程有：
22
1u p p e ρ+
=由工作条件可得空气密度，所以可得：
205.1=ρPa u p p e 75.9849350205.12
1
1000002122=⨯⨯-=-
=ρ由，可得：
%5/2
'=u u
25
.6)50*05.0()%5(222'==⨯=u u 所以： 53
.725.6*205.12
;===u
p t ρ有效压力：
Pa
p p p t eff 28.9850153.775.98492=+=+=两者差别：
%
0076.075
.9849253
.7==
-=
p
p
p eff δ9－2
公式（9－21）中的产生项：
⎥⎦⎤
⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂+⎪⎪⎭
⎫ ⎝⎛∂∂+∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂+⎢⎣⎡⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂+∂∂∂∂=⎪⎪⎭
⎫ ⎝⎛∂∂+∂∂∂∂z w z w z w z v y w y w z u x w x w y w z v z v y v y v y v y u x v x v x w z u z u x v y u y u x u x u x u x u x u x u t j
i
i j i j t
ηη9－4
由Prandtl 混合长度理论可知，湍流的切应力公式如下：
2
2
2/m w m w l y u y u l ρτρτ=∂∂⇒⎪⎪⎭
⎫ ⎝⎛∂∂=当采用模型时，有：
ε-k ε
νμ/2k c t =当脉动动量的产生与耗散相平衡时有：
2
⎪
⎪⎭
⎫ ⎝⎛∂∂=y u t νε将上述表达式代入有：
（1）
)/()/(2
2m w l k c ρτεεμ⋅=再考虑之间关系：k ,εm
l k c 2
/3μ
ε=将上式代入式（1）中整理可得：
2
/12/1/k
c w μρτ=即证。

9－5
脉动动能耗散率在直角坐标系中展开式为：⎪⎪⎭
⎫ ⎝⎛∂∂=k
i x
u '
νε⎥⎥⎦
⎤⎢⎢⎣
⎡⎪⎪⎭⎫ ⎝⎛∂∂+⎪⎪⎭⎫ ⎝
⎛∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+⎪⎪⎭⎫ ⎝⎛∂∂+⎪⎪⎭⎫ ⎝⎛∂∂=2'2'
2
'2
'2
'2
'2
'2
'2
'
z w y w x w z v y v x v z u y u x
u
νε由表达可知其量纲为：
ε3
23
22
2
-⋅=⎪⎭⎫
⎝⎛⋅⋅s m t l l t l t
l 考查的控制方程：
εk c x u x
u x
u k c x x x u t i j j i j i t k t
k
k k 221ερ
εεσηηερερε
-⎪⎪⎭
⎫ ⎝⎛∂∂+∂∂∂∂+⎥⎦⎤⎢⎣⎡∂∂⎪⎪⎭⎫
⎝
⎛+∂∂=∂∂+∂∂右侧第一项的量纲： 4
14321--⋅⋅⋅=⎪⎪⎭⎫ ⎝⎛⋅⋅⋅s m kg t
l m
l t l l t m l 第二项除外量纲：
1c ()
()()4
14
2
32///--⋅⋅⋅=⋅⋅⋅⋅
⋅s m kg t l m
l t l l t l l t m t l t l 第三项除外量纲： 2c (
)
()4
14
22
323//--⋅⋅⋅=
⋅s m kg t l m
t l t l l m 所以，和均为无量纲数；
1c 2c 由可知右侧除外量纲：，与的量纲同，所以也是一个ενμ/2
k c t =μc ()()
t
l t l
t l 2
3
2
4
//=v μc 无量纲数；
即证。