T178.8采用区域离散方法A时;内点采用中心差分T277
T3269.9d2TTi+12Ti+Ti1T=0 有 Ti0 22dxx将2点,3点带入 T2T2+T11 3 即T0T2T20 232x9T2T3+T2T42T3+T211 4 即T0T0T2T3T20 334x2x29边界点4
dT1 (1)一阶截差 由x=1 1,得 T4T3
dx3qxxVx (2)二阶截差 TM1TM1SB
111.1. 所以 T4T336T43
1112 即 2T42T3
93采用区域离散方法B
d2TT=0 由控制容积法 dx2dTdTTx0 dTdTwe 所以代入2点4点有
TTTT19 3221T20 即 T2T30
1128336TTTT3199 544T3T4 T50 T40 即
112828363 对3点采用中心差分有
T42T3+T2132T30 即
99T2T3T40 1919dT11,得 T5T4 dx6(1)精确解求左端点的热流密度
e由 T2exex
e1 对于点5 由x=1
所以有 qdTdxx0e2exxee0.648069 22e1e1x0(2)由A的一阶截差公式
q dTdxx0T2T10.247730.7431 13(3)由B的一阶截差公式
q dTdxx00.216400.6492 13(4)由区域离散方法B中的一阶截差公式:
T2T1dT0.108460.6504 dxB(x)B通过对上述计算结果进行比较可得:区域离散B有控制容积平衡法建立的离散方程与区域离散方程A中具有二阶精度的格式精确度相当!
4-3
解: 对平板最如下处理:
1 2 3 4
由左向右点分别表述为1、2、3、4点,x的正方向为由左向右; 控制方程为
λ𝑑𝑥2+𝑆=0 (1)
边界条件为
𝑑2𝑡
X=0,T=75℃;X=0.1,λ𝑑𝑥+ℎ(T−Tf)=0;
则2、3点采用二阶截差格式,有 则有以下两式:
𝑑𝑇
λλ𝑇3−2𝑇2+𝑇1∆𝑥2𝑇4−2𝑇3+𝑇2
∆𝑥2T4−T3∆x
+𝑆=0 (2) +𝑆=0 (3)
𝑑𝑇𝑑𝑥
一阶截差公式可由λλ(
+ℎ(T−Tf)=0变形得到
)=h(T4−Tf)再变形得到
h×∆xλ
T4=[T3+Tf]/(1+
∆𝑥2h×∆xλ
) (4)
T5−T32∆x
二阶截差公式可以联立λTf),可得以下公式
𝑇5−2𝑇4+𝑇3
+𝑆=0和λ(
)=h(T4−
T4=[T3+
∆𝑥2𝑆2λ
+
h×∆xλ
]/(1+
h×∆xλ
) (5)
分别联立2、3、4式与2、3、5式,把S=50×103W/m3,λ=10W/m∙℃,h=50 W/m∙℃,Tf=25℃,𝑇1=75℃,∆x=1/30带入到式子中,则有 联立2、3、4式的解为:
𝑇2=78.58℃,𝑇3=76.59℃,𝑇4=69.03℃ 联立3、4、5式的解为:
𝑇2=80.42℃,𝑇3=80.28℃,𝑇4=74.58℃
对控制方程进行积分,并将边界条件带入,则有关于T的方程 𝑇=−2500𝑥2+250𝑥+75 (6) 把x2=
130
,x3=
230
,x3=0.1代入上述6式则有:
𝑇2=80.56℃,𝑇3=80.56℃,𝑇4=75.1℃
相比之下,对右端点采用二阶截差的离散更接近真实值
4-4
解:对平板作如下分析:
1 2 3 4 5 由左向右分别对点编号为1、2、3、4、5 控制方程与4-3相同,为
λ𝑑𝑥2+𝑆=0 (1)
边界条件为
𝑑2𝑡
X=0,T=75℃;X=0.1,λ𝑑𝑥+ℎ(T−Tf)=0;
设1点和2点的距离为∆x,另1点对2点进行泰勒展开,有
𝑑2𝑡𝑑𝑥2𝑑𝑇
=(𝑇1−𝑇2+𝑑𝑥∆x)∆x2
𝑑𝑇𝑑𝑥
𝑑𝑇2
其中
=
𝑇3−𝑇22∆x
,则有
λλ2𝑇1−3𝑇2+𝑇3
∆x2+𝑆=0 (2)
对3点进行离散有
𝑇4−2𝑇3+𝑇2
∆𝑥2+𝑆=0 (3)
𝐴𝑇𝑓+(𝛿𝑥)5对右端点有: [𝑎𝑝+1𝐴+
(𝛿𝑥)5ℎλ]𝑇4=𝑎𝑤𝑇3+[𝑆/∆𝑥+1]
ℎλ代入数据有
𝑇3−3𝑇2+155.56=0 𝑇4−2𝑇3+𝑇2=−5.56
342.85𝑇4-300𝑇3=1681
解得:𝑇2=78.1℃,𝑇3=78.7℃,𝑇4=73.8℃ 由导热定律有
𝑇4−𝑇3𝑇5−𝑇4
=2 ∆𝑥∆𝑥则有𝑇5=71.35℃
4—12 编写程序:
M=rand(10,3)
A=M(:,1);B=M(:,2);C=M(:,3); B(10)=0;C(1)=0;T=12:21; D(1)=A(1)*T(1)-B(1)*T(2) for i=2:9;
D(i)= A(i)*T(i)-B(i)*T(i+1)-C(i)*T(i-1) end
D(10)= A(10)*T(10)-C(10)*T(9); P(1)=B(1)/A(1);Q(1)= D(1)/A(1); for i=2:10;
P(i)=B(i)/(A(i)-C(i)*P(i-1));
Q(i)=(D(i)+C(i)*Q(i-1))/(A(i)-C(i)*P(i-1)); end
for i=10:-1:2; t(10)=Q(10);
t(i-1)=P(i-1)*t(i)+Q(i-1); end
disp(D(1:10)) disp(T(1:10)) disp(t(1:10))
运行结果:
由运行结果可知:无论系数怎样变化,T与t都是一致的。
编程题:
具体编程如下:
n=input('请输入网格数:n=');
dt=input('请输入时间步长: dt ='); % 设定区域的长宽
L=input('请输入边界长度:L ='); a=input('请输入材质热扩散率:a ='); S=input('请输时间层数:S ='); N=n+1 ; % 节点数 dx=L/n; % 空间步长 F=a*dt/(dx^2); % 傅里叶数
% %%%%%%%%%%%%%%%%%%%%%%初始条件%%%%%%%%%%%%%%%%%%%%%%%%% for t=1:1:S for i=1:1:N if t==1
T(i,t)=100; else if i==1||i==N T(1,t)=0; T(N,t)=10; end end end end
% %%%%%%%%%%%%%%%%%%%%%初始条件%%%%%%%%%%%%%%%%%%%%%%%%% A=2*F+1;B=F;C=F; for t=2:1:S
% %%%%%%%%%%%%%%%%%%%计算各个时层的P、Q%%%%%%%%%%%%%%%%%%% for i=1:1:N-1
D(i,2)=T(i,1); if i==1
P(1,t)=0;Q(1,t)=T(1,t); else
P(i,t)=B/(A-C*P((i-1),t));
Q(i,t)=(D(i,t)+C*Q((i-1),t))/(A-C*P((i-1),t)); end end
% %%%%%%%%%%%%%%%%%%%%%计算各个时层的P、Q%%%%%%%%%%%%%%%%% % %%%%%%%%%%%%%%%%%%%计算各个时层的温度%%%%%%%%%%%%%%%%%%% for i=N-1:-1:2
T(i,t)=P(i,t)*T((i+1),t)+Q(i,t); D(i,(t+1))=T(i,t); end end
x=0:0.01:1;
% 生成其中几个时层图像
plot(x,T(1:101,1),':bo',x,T(1:101,2),':ro',x,T(1:101,3),':bo',x,T(1:101,4),':ro',x,T(1:101,5),':bo',x,T(1:101,6),':ro') xlabel('x'); ylabel('T');
title('一维非稳态热传导方程的数值解')
输入节点数为100,时间步长为0.1,边界长度为1,热扩散率为1,时间层数为50
输入节点数为100,时间步长为1,边界长度为1,热扩散率为1,时间层数为50
与显式比较
输入节点数为100,时间步长为0.1
例题:
温度为Tf =90 80 ℃的水进入一管径为0.1 0.12m的长圆管作层流流动,
he=8w/(m2 ℃)。 平均速度为0.6m/s。管外受到温度为T=30 ℃的流体的冷却,
设换热已进入充分发展阶段:1、试确定管内对流换热系数he;2、给出x=1m截面上水温沿半径方向的分布;3、平均水温沿x方向的分布。要求:确定网格独
立解的节点数。
由题意可列出方程:
T1Tcpu(r)xrrr
Tr0,0r边界条件:TrR,he(TT)r
无量纲化:TT2Rumurx22(1) ,X,,P,eTbTumRaRPedd2()(1)0ddd0,0dd)1Biwd11/(4(12)d)0
对控制方程积分,写成追赶法的形式:
2r1rr1Ti()Ti1()Ti1*(12)r r2rr2T1T2边界条件:
(1Br)TTiL1L117查得90 80 ℃水的导热系数λ=0.68,导温系数1.6210,um=0.6,
R=0.05,
*设初始的Φ*=1,代入求解得到Φ,直到103
dTb(TbT)Tb30CeX dX
x0.37105由边界条件X=0时,Tb=90 C=60
Tb3060e
编程为:
取x=1时,解得Tb T(Tb30)30
%phi为φ,eta为η,phi1为φ*
L1=input('输入节点个数:'); R=0.05;
%外节点法划分网格 n=R/(L1-1); i=2:L1; X(i)=(i-1)*n;
X(1)=0; X(L1)=R; eta=X/R; eta(1)=0;
%迭代部分
phi1=ones(1,L1); for i=2:L1-1 A=2*X/n; B=0.5+X/n; C=X/n-0.5;
Lambda=1/(4*sum(phi1.*eta.*(1-eta.^2)*(n/R))); S=Lambda.*eta.*(1-eta.^2).*phi1; D=S.*n; end A(1)=1; B(1)=1; C(1)=0; D(1)=0;
P(1)=B(1)/A(1); Q(1)=D(1)/A(1); A(L1)=1+8*n/0.68; B(L1)=0; C(L1)=1; D(L1)=0; %追赶法解方程组 for i=2:L1
P(i)=B(i)/(A(i)-C(i)*P(i-1));
Q(i)=(D(i)+C(i)*Q(i-1))/(A(i)-C(i)*P(i-1)); end
phi(L1)=Q(L1); for i=L1:-1:2
phi(i-1)=P(i-1)*phi(i)+Q(i-1); end
while abs((phi-phi1)./phi)>=10e-3 phi1=phi;
Lambda=1/(4*sum(phi1.*eta.*(1-eta.^2).*(n/R))); S=Lambda.*eta.*(1-eta.^2).*phi1; D=S.*n; for i=2:L1
P(i)=B(i)/(A(i)-C(i)*P(i-1));
Q(i)=(D(i)+C(i)*Q(i-1))/(A(i)-C(i)*P(i-1)); end
phi(L1)=Q(L1); for i=L1:-1:2
phi(i-1)=P(i-1)*phi(i)+Q(i-1); end end
Tb=30+60*exp(-Lambda/3.7037e+5); T=phi*Lambda*(Tb-30)+30;
h=-0.68*(T(L1-1)-T(L1))/n/(T(L1)-Tb) figure(1); plot(T,X)
title('1m处水温延半径的分布'); xlabel('温度'); ylabel('半径'); %平均水温
x=linspace(0,10);
Tb=30+60*exp(-Lambda.*x/3.7037e+5); figure(2); plot(x,Tb)
title('平均水温延x的分布'); xlabel('圆管长度'); ylabel('平均水温'); 输入节点个数:100 可得:h =29.8341
x=1m截面上水温沿半径方向的分布图如下:
平均水温沿x方向的分布
因篇幅问题不能全部显示,请点此查看更多更全内容