传热学数值计算.docx

传热学数值计算.docx

《传热学数值计算.docx》由会员分享，可在线阅读，更多相关《传热学数值计算.docx（47页珍藏版）》请在冰豆网上搜索。

传热学数值计算

传热学数值计算作业

数值解程序：

tw1=40%三边温度

tw2=100%一边温度正弦变化幅度

l1=40%板长L1:

40厘米

l2=20%板宽L2:

20厘米

m=41%分划成40*20的网格

n=21

k=2

dx=l1/（m-1）

c=ones（n,m）

fori=1:

m

a2（i）=tw1+tw2*sin（pi*dx*（i-1）/l1）

c（1,i）=tw1,c（n,i）=a2（i）

end

forj=1:

n

c（j,1）=tw1

c（j,m）=tw1

end

while（abs（c（j,i）-k）>0.0001）

k=c（j,i）

fori=2:

m-1

forj=2:

n-1

c（j,i）=0.25*（c（j,i-1）+c（j,i+1）+c（j-1,i）+c（j+1,i））

end

end

end

数值解中各网格点的温度值：

数值二维温度分布图像：

解析解程序：

tw1=40

tw2=100

l1=40

l2=20

p=40

q=20

x

（1）=0

fori=1:

p

x（i+1）=x（i）+1

end

y

（1）=0

forj=1:

q

y（j+1）=y（j）+1

end

fori=1:

p+1

forj=1:

q+1

n（j,i）=tw1+tw2*sinh（pi*y（j）/l1）*sin（pi*x（i）/l1）/sinh（pi*l2/l1）

end

end

各网格点用解析式得到的温度值：

解析二维温度分布图像：

误差分析：

取x=21，即位于板长一半处，温度随y（宽度）的变化曲线。

c1（:

1）取自于数值解，c1（:

2）取自于解析解

c1（:

1）c1（:

2）

40.000040.0000

43.310643.4164

46.646546.8538

50.031350.3335

53.488953.8771

57.043057.5062

60.717861.2434

64.537665.1117

68.527369.1350

72.712273.3381

77.118777.7470

81.773682.3888

86.705087.2922

91.942392.4875

97.516298.0068

103.4592103.8840

109.8058110.1555

116.5925116.8600

123.8586124.0388

131.6461131.7363

140.0000140.0000

误差曲线：

由相对误差公式：

d1=（c1（:

2）-c1（:

1））./c1（:

2）可得：

d1=

0

0.0024

0.0044

0.0060

0.0072

0.0081

0.0086

0.0088

0.0088

0.0085

0.0081

0.0075

0.0067

0.0059

0.0050

0.0041

0.0032

0.0023

0.0015

0.0007

0

结论：

数值解与解析解吻合很好。

就x=21这一列，相对误差较小且在1%以内，但数值解较解析解偏小，且在平板中心附近的网格点的数值解较平板边缘数值解的相对误差大。

数值解程序：

tw1=50%三边温度

tw2=100%一边温度

l1=20%板长20厘米

l2=10%板宽10厘米

m=41%划分成40*20的网格

n=21

k=1

c=zeros（n,m）

c（n,1）=（tw1+tw2）/2

c（n,m）=（tw1+tw2）/2

c（1,1）=tw1

c（1,m）=tw1

fori=2:

（m-1）

c（1,i）=tw1

c（n,i）=tw2

end

forj=2:

（n-1）

c（j,1）=tw1

c（j,m）=tw1

end

while（abs（k-c（j,i））>0.0001）

k=c（j,i）

fori=2:

m-1

forj=2:

n-1

c（j,i）=0.25*（c（j,i-1）+c（j,i+1）+c（j-1,i）+c（j+1,i））

end

end

end

数值解中各网格点的温度值：

数值二维温度分布图像：

解析解程序：

tw1=50

tw2=100

m=40

n=20

l1=20

l2=10

tx=ones（20,40）

fori=1:

m+1

forj=1:

n+1

x=（i-1）*0.5

y=（j-1）*0.5

k=sym（'k'）

d=（（（-1）^（k+1）+1）/k）*sin（k*pi*x/l1）*sinh（k*pi*y/l1）/sinh（k*pi*l2/l1）

h=symsum（d,k,1,200）

tx（j,i）=2*h*（tw2-tw1）/pi+tw1

end

end

各网格点用解析式得到的温度值：

解析二维温度分布图像：

误差分析：

取x=21，即位于板长一半处，温度随y（宽度）的变化曲线。

c1（:

1）取自于数值解c1（:

2）取自于解析解

c1（:

1）=c1（:

2）=

50.000050.0000

51.979952.0881

53.973154.1851

55.990656.2998

58.043558.4409

60.141860.6165

62.295162.8344

64.511665.1016

66.798767.4243

69.162269.8077

71.606472.2558

74.133774.7710

76.744777.3546

79.437780.0056

82.209082.7214

85.052685.4977

87.960288.3278

90.921491.2035

93.924494.1151

96.955497.0513

100.000099.8408

误差分析曲线：

由相对误差公式：

d2=abs（c1（:

2）-c1（:

1））./c1（:

2）可得：

d2=

0

0.0021

0.0039

0.0055

0.0068

0.0078

0.0086

0.0091

0.0093

0.0092

0.0090

0.0085

0.0079

0.0071

0.0062

0.0052

0.0042

0.0031

0.0020

0.0010

0.0016

误差结论：

数值解与解析解吻合很好。

就x=21这一列，相对误差较小且在1%以内，但数值解较解析解普遍偏小（最后一个数除外），且在平板中心附近的网格点的数值解较平板边缘数值解的相对误差大。

数值解程序：

t0=200;%肋基温度为200度

tf=0;%环境温度为0度

L1=20;%肋片长20厘米

L2=2;%肋片高2厘米

h=0.01;%对流换热系数单位：

w/（cm2*K）

a=11;

b=101;%11*101的网格

dx=L1/（b-1）;

k=2;%导热系数单位：

w/（cm*K）

Bi=h*dx/k;

ti=ones（a,b）*10;

m1=ones（a,b）*3;

m1（2:

a-1,1）=zeros（a-2,1）;

m1（a,2:

b-1）=ones（1,b-2）;

m1（1,2:

b-1）=ones（1,b-2）*6;

m1（2:

a-1,b）=ones（a-2,1）*2;

m1（1,b）=ones（1,1）*4;

m1（a,b）=ones（1,1）*5;

m1（1,1）=7;

m1（a,1）=8;

tn=ti;

max1=1.0;

w=0;

while（max1>1e-6）

w=w+1;

max1=0;

fori=1:

a

forj=1:

b

m=m1（i,j）;

n=tn（i,j）;

switchm

case0

tn（i,j）=t0;

case1

tn（i,j）=（2*tn（i-1,j）+tn（i,j-1）+tn（i,j+1）-4*tf）/（4+2*Bi）+tf;

case2

tn（i,j）=（2*tn（i,j-1）+tn（i-1,j）+tn（i+1,j）-4*tf）/（4+2*Bi）+tf;

case3

tn（i,j）=0.25*（tn（i,j-1）+tn（i,j+1）+tn（i-1,j）+tn（i+1,j））;

case4

tn（i,j）=（tn（i,j-1）+tn（i+1,j）-2*tf）/（2*Bi+2）+tf;

case5

tn（i,j）=（tn（i,j-1）+tn（i-1,j）-2*tf）/（2*Bi+2）+tf;

case6

tn（i,j）=（2*tn（i+1,j）+tn（i,j-1）+tn（i,j+1）-4*tf）/（4+2*Bi）+tf;

case7

tn（i,j）=t0;

case8

tn（i,j）=t0;

end

er=abs（tn（i,j）-n）;

ifer>max1

max1=er;

end

end

end

ti=tn;

end

ti

%假设垂直于纸面方向肋宽为1cm

q1=0;

fori=1:

a

qi=（tn（i,b）-tf）*h*0.2;

q1=q1+qi;

end

q1%x=20cm处肋片的散热量

q2=0;

forj=1:

b-1

q2j=（tn（a,j）-tf）*h*0.2;

q2=q2+q2j;

end

q2=q2*2;

q2%y=1cm与y=-1cm处肋片的散热量

q=q1+q2;

q

数值解中各网格点的温度值：

数值二维温度分布图像：

数值解肋片换热量：

q1=1.9022w/cm肋端沿y方向

q2=49.4938w/cm肋端沿x方向

q=51.3960w/cmq=q1+q2

解析解程序：

t0=200;%肋基温度为200度

tf=0;%环境温度为0度

L1=20;%肋片长20厘米

L2=2;%肋片高2厘米

h=0.01;%对流换热系数单位：

w/（cm2*K）

a=11;

b=101;%11*101的网格

dx=L1/（b-1）;

k=2;%导热系数单位：

w/（cm*K）

Bi=h*1/k;

ta=ones（a,b）;

fori=1:

1:

a

forj=1:

1:

b

ifi>（a+1）/2

y=-（i-（a+1）/2）*dx;

elsey=（（a+1）/2-i）*dx;

end

x=dx*（j-1）;

ta（i,j）=（cosh（Bi^0.5*（L1-x）/1）+Bi^0.5*sinh（Bi^0.5*（L1-x）/1））*（t0-tf）/（cosh（Bi^0.5*L1/1）+Bi^0.5*sinh（Bi^0.5*L1/1））+tf;

end

end

ta

%假设垂直于纸面方向肋宽为1cm

q1=0;

fori=1:

1:

a

qi=（ta（i,b）-tf）*h*0.2;

q1=q1+qi;

end

q1%x=20cm处肋片的散热量

q2=0;

forj=1:

1:

b-1

q2j=（ta（a,j）-tf）*h*0.2;

q2=q2+q2j;

end

q2=q2*2;

q2%y=1cm与y=-1cm处肋片的散热量

q=q1+q2;

q

各网格点用解析式得到的温度值：

解析解二维温度图像：

解析解肋片换热量：

q1=1.9006w/cm肋端沿y方向

q2=49.5481w/cm肋端沿x方向

q=51.4488w/cmq=q1+q2

误差分析：

取y=0，即位于板宽一半处，温度随x（长度）的变化曲线。

b1（:

1）取自于数值解b1（:

2）取自于解析解

b1（:

1）=b1（:

2）=

200.0000200.0000

197.5383197.4675

195.0883194.9745

192.6596192.5205

190.2587190.1051

187.8895187.7276

185.5543185.3877

183.2542183.0848

180.9896180.8186

178.7605178.5885

176.5667176.3942

174.4081174.2351

172.2842172.1109

170.1946170.0211

168.1389167.9654

166.1169165.9432

164.1280163.9542

162.1718161.9979

160.2480160.0741

158.3563158.1823

156.4961156.3222

154.6672154.4933

152.8692152.6953

151.1017150.9278

149.3644149.1906

147.6569147.4831

145.9788145.8052

144.3300144.1564

142.7099142.5365

141.1184140.9451

139.5550139.3818

138.0195137.8464

136.5115136.3386

135.0308134.8581

133.5771133.4046

132.1500131.9777

130.7494130.5772

129.3748129.2028

128.0261127.8543

126.7029126.5314

125.4051125.2337

124.1322123.9611

122.8842122.7133

121.6607121.4901

120.4615120.2911

119.2863119.1162

118.1350117.9651

117.0073116.8376

115.9029115.7335

114.8216114.6525

113.7633113.5945

112.7277112.5592

111.7147111.5464

110.7239110.5559

109.7552109.5875

108.8085108.6410

107.8835107.7163

106.9800106.8131

106.0979105.9312

105.2369105.0706

104.3970104.2309

103.5780103.4122

102.7796102.6140

102.0017101.8365

101.2442101.0793

100.5069100.3423

99.789799.6253

99.092598.9283

98.415098.2511

97.757197.5935

97.118896.9555

96.499996.3368

95.900295.7374

95.319795.1572

94.758394.5960

94.215794.0537

93.692093.5302

93.187093.0254

92.700692.5393

92.232792.0716

91.783291.6224

91.352091.1915

90.939190.7788

90.544390.3842

90.167790.0078

89.809089.6494

89.468389.3088

89.145488.9862

88.840388.6813

88.553088.3942

88.283488.1248

88.031487.8729

87.796987.6387

87.580087.4220

87.380687.2228

87.198787.0410

87.034186.8766

86.886986.7296

86.757186.6000

86.644686.4876

86.549486.3926

误差分析曲线：

曲线适当放大：

由相对误差公式：

d2=abs（b1（:

2）-b1（:

1））./b1（:

2）可得：

d2=

0

0.0004

0.0006

0.0007

0.0008

0.0009

0.0009

0.0009

0.0009

0.0010

0.0010

0.0010

0.0010

0.0010

0.0010

0.0010

0.0011

0.0011

0.0011

0.0011

0.0011

0.0011

0.0011

0.0012

0.0012

0.0012

0.0012

0.0012

0.0012

0.0012

0.0012

0.0013

0.0013

0.0013

0.0013

0.0013

0.0013

0.0013

0.0013

0.0014

0.0014

0.0014

0.0014

0.0014

0.0014

0.0014

0.0014

0.0015

0.0015

0.0015

0.0015

0.0015

0.0015

0.0015

0.0015

0.0015

0.0016

0.0016

0.0016

0.0016

0.0016

0.0016

0.0016

0.0016

0.0016

0.0016

0.0017

0.0017

0.0017

0.0017

0.0017

0.0017

0.0017

0.0017

0.0017

0.0017

0.0017

0.0017

0.0017

0.0017

0.0018

0.0018

0.0018

0.0018

0.0018

0.0018

0.0018

0.0018

0.0018

0.0018

0.0018

0.0018

0.0018

0.0018

0.0018

0.0018

0.0018

0.0018

0.0018

0.0018

0.0018

肋片总散热量数值解与解析解的相对误差为

q（解析解）=51.4488w/cm

q（数值解）=51.3960w/cm

t=abs（q（解析解）-q（数值解））/q（解析解）=0.001026

误差结论：

数值解与解析解吻合很好。

就y=0这行，相对误差较小且在0.2%以内，但数值解较解析解普遍偏大，且离肋基较近的网格点的数值解比离肋基较远的数值解相对误差大。

解析解的肋片总散热量与数值解吻合较好，但总散热量比数值解大。

数值解程序：

tw=0;%环境温度0度

L2=10;%板宽的一半10厘米

c=0.5;

L1=L2/c;%板长的一半20厘米

p2=40;%划分成81*41的网格

p1=p2/c;

dx=2*L2/p2;

a=p2+1;

b=p1+1;

k=2;%导热系数单位w/（m2*K）

q=3;%单位体积的产热率w/m3

ti=ones（a,b）*5;

m1=ones（a,b）;

m1（1,:

）=zeros（1,b）;

m1（2:

a,b）=zeros（a-1,1）;

m1（2:

a,1）=zeros（a-1,1）;

m1（a,2:

b-1）=zeros（1,b-2）;

tn=ti;

max1=1.0;

y=0;

while（max1>1e-6）

max1=0;

y=y+1;

fori=1:

1:

a

forj=1:

1:

b

m=m1（i,j）;

n=tn（i,j）;

switchm

case0

tn（i,j）=tw;

case1

tn（i,j）=0.25*（tn（i-1,j）+tn（i+1,j）+tn（i,j-1）+tn（i,j+1）+q*（dx^2）/k）;

end

er=abs（tn（i,j）-n）;

ifer>max1

max1=er;

end

end

end

ti=tn;

end

ti

a1=ti（21:

41,41:

81）

数值解中各网格点的温度值：

数值二维温度分布图像：

解析解程序：

tw1=0;

L1=20;

L2=10;

m=41;

n=21;

dx=L1/（m-1）;

k=2;

q=3;

tx=ones（n,m）;

fori=1:

m

forj=1:

n

x=（i-1）*0.5;

y=（j-1）*0.5;

p=sym（'p'）;

xi=（2*p-1）*pi/2;

d=（（-1）^p）/（xi^3）*（cosh（xi*y/L1）/cosh（xi*L2/L1））*cos（xi*x/L1）;

h=symsum（d,p,1,100）;

tx（j,i）=2*q*L1^2/k*h+q*（L1^2-x^2）/（2*k）+tw1;

end

end

tx

各网格点用解析式得到的温度值：

解析二维温度分布图像：

误差分析：

取x=1，即位于板长一半处，温度随y（宽度）的变化曲线。

c1（:

1）取自于数值解c1（:

2）取自于解析解

c1（:

1）=c1（:

2）=

68.312168.3231

68.145268.1562

67.644467.6553

66.809366.8200

65.6393