一维导热计算

非稳态

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
open "a1.dat" for output as #1
ni=22:nim1=ni-1
dim x(ni),sew(ni),dxep(ni),dxpw(ni)
dim ae(ni),aw(ni),ap0(ni),ap1(ni),ap(ni)
dim t0(ni),t1(ni),su(ni),sp(ni),den0(ni),cp0(ni)
dim lam(ni),den1(ni),cp1(ni),p(ni),q(ni)
dx=0.5
x(1)=-dx/2
for i=2 to ni
x(i)=x(i-1)+dx
next
print #1,ni
for i=1 to ni
print #1,using"##.####";x(i)
next
dxpw(1)=0.0
dxep(ni)=0.0
for i=1 to nim1
dxep(i)=x(i+1)-x(i)
dxpw(i+1)=dxep(i)
next 
sew(1)=0.0
sew(ni)=0.0
for i=2 to nim1
sew(i)=0.5*(dxep(i)+dxpw(i))
next
tf=1000:t0=20:lamda0=20:alafa=2320:b=0.0:eps=0.01
densit=7000:cp0=15:dt=1:sel=1 //sel=1时为隐式，sel=2时为显式
for i=2 to nim1
t0(i)=t0
next i
time=0
10 time=time+dt
for i=2 to nim1
lam(i)=lamda0*(1+b*t0(i))
den0(i)=densit
cp0(i)=cp0
den1(i)=densit
cp1(i)=cp0
next i
for i=2 to nim1
ae(i)=2*lam(i)*lam(i+1)/(lam(i)+lam(i+1))/dxep(i)
aw(i)=2*lam(i)*lam(i-1)/(lam(i)+lam(i-1))/dxpw(i)
ap0(i)=den0(i)*cp0(i)/dt*sew(i)
ap1(i)=den1(i)*cp1(i)/dt*sew(i)
su(i)=0
sp(i)=0
next i
i=2
aw(i)=0
t1(1)=t1(2)
t0(i)=t0(2)
i=nim1
ae(i)=alfa
t1(ni)=tf
t0(ni)=tf
if sel=1 then
for i=2 to nim1
ap(i)=ap1(i)+ae(i)+aw(i)-sp(i)*sew(i)
next i
p(2)=ae(2)/ap(2)
q(2)=(aw(2)*tf+ap0(2)*t0(2)+su(2)*sew(2))/ap(2)
for i=3 to nim1
p(i)=ae(i)/(ap(i)-aw(i)*p(i-1))
q(i)=(aw(i)*q(i-1)+ap0(i)*t0(i)+su(i)*sew(i))/(ap(i)-aw(i)*p(i-1))
next i
t1(nim1)=q(nim1)
for i=ni-2 to 2 step -1
t1(i)=p(i)*t1(i+1)+q(i)
next i
end if
if sel=2 then
for i=2 to nim1
ap(i)=ap0(i)-ae(i)-aw(i)+sp(i)*sew(i)
if ap(i)<0 then stop
t1(i)=(ae(i)*t0(i+1)+aw(i)*t0(i-1)+su(i)*sew(i)+ap(i)*t0(i))/ap1(i)
next i
end if
print using"####.##";time,t1(2),t1(nim1)
print #1,time
for i=2 to nim1
print #1,using"####.##";t1(i)
next i
print #1,
for i=2 to nim1
t0(i)=t1(i)
next i
if time<3600 then 10
print time
for i=2 to nim1
print #1,using"####.##";t1(i)
next i
end

稳态

1
参考二维稳态导热

二维导热计算

稳态

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
open "a1.dat" for output as #1
open "a2.dat" for output as #2
ni=42:nj=42:nim1=ni-1:njm1=nj-1
dim x(ni),y(nj),sew(ni),sns(nj),dxep(ni),dxpw(ni),dynp(nj),dyps(nj)
dim ae(ni,nj),aw(ni,nj),as1(ni,nj),an(ni,nj),ap(ni,nj)
dim t0(ni,nj),t1(ni,nj),su(ni,nj),sp(ni,nj),lam(ni,nj)
dx1=0.2:dx2=0.3:dx3=0.01:dy1=0.05:dy2=0.05
x(1)=-dx1/2
for i=2 to 10
x(i)=x(i-1)+dx1
next i
for i=11 to 18
x(i)=x(i-1)+dx2
next i
for i=19 to ni
x(i)=x(i-1)+dx3
next i
y(1)=-dy1/2
for j=2 to 5
y(j)=y(j-1)+dy1
next j
for j=6 to nj
y(j)=y(j-1)+dy2
next j
print #1,ni,nj
for i=1 to ni
print using"##.####";x(i)
next i
for j=1 to nj
print #1,using"##.####";y(j)
next j
for i=2 to nim1
dxep(i)=x(i+1)-x(i)
dxpw(i)=x(i)-x(i-1)
next i
for j=2 to njm1
dynp(j)=y(j+1)-y(j)
dyps(j)=y(j)-y(j-1)
next j
for i=2 to nim1
sew(i)=0.5*(dxep(i)+dxpw(i))
next i
for j=2 to njm1
sns(j)=0.5*(dynp(j)+dyps(j))
next j
tf=1000:t0=100
lamda0=20:alfa=200:b=0.1:eps=0.01
niter=0
for i=2 to nim1
for j=2 to njm1
t0(i,j)=t0
next j,i
10 niter=niter+1
for i=2 to nim1
for j=2 to njm1
lam(i,j)=lamda0*(1+b*t0(i,j))
mext j,i
ae(i,j)=2*lam(i,j)*lam(i+1,j)/(lam(i,j)+lam(i+1,j))*sns(j)/dxep(i)
aw(i,j)=2*lam(i,j)*lam(i-1,j)/(lam(i,j)+lam(i-1,j))*sns(j)/dxpw(i)
as1(i,j)=2*lam(i,j)*lam(i,j-1)/(lam(i,j)+lam(i,j-1))*sew(i)/dyps(j)
an(i,j)=2*lam(i,j)*lam(i,j+1)/(lam(i,j)+lam(i,j+1))*sew(i)/dynp(j)
su(i,j)=0
sp(i,j)=0
next j,i
i=2
for j=2 to njm1
su(i,j)=1.0E30*100
sp(i,j)=-1.0E30
next j
i=nim1
for j=2 to njm1
su(i,j)=1.0E30*300
sp(i,j)=-1.0E30
next j
j=2
for i=2 to nim1
as1(i,j)=0
t0(i,1)=t0(i,j)
next i
j=njm1
for i=2 to nim1
an(i,j)=alfa*sew(i)
t0(1,nj)=tf
next i
for i=2 to nim1
for j=2 to njm1
ap(i,j)=ae(i,j)+aw(i,j)+as1(i,j)+an(i,j)-sp(i,j)*sew(i)*sns(j)
t1=ae(i,j)*t0(i+1,j)+aw(i,j)*t0(i-1,j)+as1(i,j)*t0(i,j-1)
t1=t1+an(i,j)*t0(i,j+1)+su(i,j)*sew(i)*sns(j)
t1(i,j)=t1/ap(i,j)
next j,i
num=0
for i=2 to nim1
for j=2 to njm1
if abs(t1(i,j)-t0(i,j))>eps then num=num+1
t0(i,j)=t1(i,j)
next j,i
print suing"###################";niter;num
if num>10 then goto 10
for j=njm1 to 2 step -1
for i=2 to nim1
print using"####.#";t1(i,j)
print #2,using"####.####";y(j),x(i),t1(i,j)
next 
print
next
end

非稳态-显式

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
open "a1.dat" for output as #1
open "a2.dat" for output as #2
ni=22:nj=12:nim1=ni-1:njm1=nj-1
dim x(ni),y(nj),sew(ni),sns(nj),dxep(ni),dxpw(ni),dynp(nj),dyps(nj)
dim ae(ni,nj),aw(ni,nj),as1(ni,nj),an(ni,nj),ap(ni,nj),ap0(ni,nj)
dim t0(ni,nj),t1(ni,nj),su(ni,nj),sp(ni,nj),ap1(ni,nj)
dim lam(ni,nj),den0(ni,nj),cp0(ni,nj),den1(ni,nj),cp1(ni,nj)
dx1=0.5:dx2=0.3:dx3=0.01:dy1=0.1:dy2=0.05
x(1)=-dx1/2
for i=2 to 10
x(i)=x(i-1)+dx1
next 
for i=11 to 18
x(i)=x(i-1)+dx2
next 
for i=19 to ni
x(i)=x(i-1)+dx3
next 
y(1)=-dy1/2
for j=2 to 5
y(j)=y(j-1)+dy1
next 
for j=6 to nj
y(j)=y(j-1)+dy2
next 
print #2 ni,nj
for i=1 to ni
print #2,using"##.####";x(i)
next 
for j=1 to nj
print #2,using"##.####";y(j)
next 
dxpw(1)=0.0
dxep(ni)=0.0
for i=1 to nim1
dxep(i)=x(i+1)-x(i)
dxpw(i+1)=dxep(i)
next 
dyps(1)=0.0
dynp(nj)=0.0
for j=1 to njm1
dynp(j)=y(j+1)-y(j)
dyps(j+1)=dynp(j)
next 
sew(1)=0.0
sew(ni)=0.0
for i=2 to nim1
sew(i)=0.5*(dxep(i)+dxpw(i))
next 
sns(1)=0.0
sns(nj)=0.0
for j=2 to njm1
sns(j)=0.5*(dynp(j)+dyps(j))
next j
tf=1000:t0=100
lamda0=20:alfa=200:b=0.0:eps=0.01
densit=7000:cp0=15:dt=0.01
niter=0
time=0
for i=2 to nim1
for j=2 to njm1
t0(i,j)=t0
next j,i
10 time=time+dt
for i=2 to nim1
for j=2 to njm1
lam(i,j)=lamda0*(1+b*t0(i,j))
den0(i,j)=densit
cp0(i,j)=cp0
den1(i,j)=densit
cp1(i,j)=cp0
next j,i
for i=2 to nim1
for j=2 to njm1
ae(i,j)=2*lam(i,j)*lam(i+1,j)/(lam(i,j)+lam(i+1,j))*sns(j)/dxep(i)
aw(i,j)=2*lam(i,j)*lam(i-1,j)/(lam(i,j)+lam(i-1,j))*sns(j)/dxpw(i)
as1(i,j)=2*lam(i,j)*lam(i,j-1)/(lam(i,j)+lam(i,j-1))*sew(i)/dyps(j)
an(i,j)=2*lam(i,j)*lam(i,j+1)/(lam(i,j)+lam(i,j+1))*sew(i)/dynp(j)
ap0(i,j)=den0(i,j)*cp0(i,j)/dt*sew(i)*sns(j)
ap1(i,j)=den1(i,j)*cp1(i,j)/dt*sew(i)*sns(j)
su(i,j)=0
sp(i,j)=0
next j,i
i=2
for j=2 to njm1
aw(i,j)=0
t0(1,j)=t0(2,j)
next j
i=nim1
for j=2 to njm1
ae(i,j)=alfa*sns(j)
t0(ni,j)=tf
next j
j=2
for i=2 to nim1
as1(i,j)=0
t0(i,1)=t0(i,j)
next i
j=njm1
for i=2 to nim1
an(i,j)=alfa*sew(i)
t0(i,nj)=tf
next i
for i=2 to nim1
for j=2 to njm1
ap(i,j)=ap0(i,j)-(ae(i,j)+aw(i,j)+as1(i,j)+an(i,j)-sp(i,j)*sew(i)*sns(j))
if ap(i,j)<0 then stop
t1=ae(i,j)*t0(i+1,j)+aw(i,j)*t0(i-1,j)+as1(i,j)*t0(i,j-1)+an(i,j)*t0(i,j+1)
t1(i,j)=(t1+su(i,j)*sew(i)*sns(j)+ap(i,j)*t0(i,j))/ap1(i,j)
next j,i
print using"#####.##";time,t1(2,2),t1(2,njm1),t1(nim1,2),t1(nim1,njm1)
print #1,using"#####.##";time,t1(2,2),t1(2,njm1),t1(nim1,2),t1(nim1,njm1)
for i=2 to nim1
for j=2 to njm1
t0(i,j)=t1(i,j)
next j,i
if time<3600 then goto 10
for i=2 to nim1
for j=2 to njm1
print #2,using"####.#";t1(i,j)
next j
print #2,
next i
end

非稳态-隐式

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
open "a1.dat" for output as #1
open "a2.dat" for output as #2
ni=22:nj=12:nim1=ni-1:njm1=nj-1
dim x(ni),y(nj),sew(ni),sns(nj),dxep(ni),dxpw(ni),dynp(nj),dyps(nj)
dim ae(ni,nj),aw(ni,nj),as1(ni,nj),an(ni,nj),ap(ni,nj),ap0(ni,nj)
dim t0(ni,nj),t1(ni,nj),su(ni,nj),sp(ni,nj),ap1(ni,nj),t10(ni,nj)
dim lam(ni,nj),den0(ni,nj),cp0(ni,nj),den1(ni,nj),cp1(ni,nj)
dx1=0.5:dx2=0.3:dx3=0.01:dy1=0.1:dy2=0.05
x(1)=-dx1/2
for i=2 to 10
x(i)=x(i-1)+dx1
next 
for i=11 to 18
x(i)=x(i-1)+dx2
next 
for i=19 to ni
x(i)=x(i-1)+dx3
next 
y(1)=-dy1/2
for j=2 to 5
y(j)=y(j-1)+dy1
next 
for j=6 to nj
y(j)=y(j-1)+dy2
next 
print #2 ni,nj
for i=1 to ni
print #2,using"##.####";x(i)
next 
for j=1 to nj
print #2,using"##.####";y(j)
next 
dxpw(1)=0.0
dxep(ni)=0.0
for i=1 to nim1
dxep(i)=x(i+1)-x(i)
dxpw(i+1)=dxep(i)
next 
dyps(1)=0.0
dynp(nj)=0.0
for j=1 to njm1
dynp(j)=y(j+1)-y(j)
dyps(j+1)=dynp(j)
next 
sew(1)=0.0
sew(ni)=0.0
for i=2 to nim1
sew(i)=0.5*(dxep(i)+dxpw(i))
next 
sns(1)=0.0
sns(nj)=0.0
for j=2 to njm1
sns(j)=0.5*(dynp(j)+dyps(j))
next j
tf=1000:t0=100
lamda0=20:alfa=200:b=0.0:eps=0.01
densit=7000:cp0=15:dt=1
niter=0
time=0
for i=2 to nim1
for j=2 to njm1
t0(i,j)=t0
t10(i,j)=t0
next j,i
10 time=time+dt
for i=2 to nim1
for j=2 to njm1
lam(i,j)=lamda0*(1+b*t0(i,j))
den0(i,j)=densit
cp0(i,j)=cp0
den1(i,j)=densit
cp1(i,j)=cp0
next j,i
for i=2 to nim1
for j=2 to njm1
ae(i,j)=2*lam(i,j)*lam(i+1,j)/(lam(i,j)+lam(i+1,j))*sns(j)/dxep(i)
aw(i,j)=2*lam(i,j)*lam(i-1,j)/(lam(i,j)+lam(i-1,j))*sns(j)/dxpw(i)
as1(i,j)=2*lam(i,j)*lam(i,j-1)/(lam(i,j)+lam(i,j-1))*sew(i)/dyps(j)
an(i,j)=2*lam(i,j)*lam(i,j+1)/(lam(i,j)+lam(i,j+1))*sew(i)/dynp(j)
ap0(i,j)=den0(i,j)*cp0(i,j)/dt*sew(i)*sns(j)
ap1(i,j)=den1(i,j)*cp1(i,j)/dt*sew(i)*sns(j)
su(i,j)=0
sp(i,j)=0
next j,i
i=2
for j=2 to njm1
aw(i,j)=0
t10(1,j)=t10(2,j)
next j
i=nim1
for j=2 to njm1
ae(i,j)=alfa*sns(j)
t10(ni,j)=tf
next j
j=2
for i=2 to nim1
as1(i,j)=0
t10(i,1)=t10(i,j)
next i
j=njm1
for i=2 to nim1
an(i,j)=alfa*sew(i)
t10(i,nj)=tf
next i
niter=0
20 niter=niter+1
num=0
for i=2 to nim1
for j=2 to njm1
ap(i,j)=ap1(i,j)+ae(i,j)+aw(i,j)+as1(i,j)+an(i,j)-sp(i,j)*sew(i)*sns(j)
t1=ae(i,j)*t10(i+1,j)+aw(i,j)*t10(i-1,j)+as1(i,j)*t10(i,j-1)+an(i,j)*t10(i,j+1)
t1(i,j)=(t1+su(i,j)*sew(i)*sns(j)+ap0(i,j)*t0(i,j)/ap(i,j)
if abs(t1(i,j)-t10(i,j))>eps then num=num+1
t10(i,j)=t1(i,j)
next j,i
print niter,num
print using"#####.##";time,t1(2,2),t1(2,njm1),t1(nim1,2),t1(nim1,njm1)
print #1,using"#####.##";time,t1(2,2),t1(2,njm1),t1(nim1,2),t1(nim1,njm1)
for i=2 to nim1
for j=2 to njm1
t0(i,j)=t1(i,j)
next j,i
if time<3600 then goto 10
for i=2 to nim1
for j=2 to njm1
print #2,using"####.#";t1(i,j)
next j
print #2,
next i
end

2个重要的理论基础

  • 相似转换-方程分析法

在解决工程问题的过程当中，常常拿实际问题去套用现有模型，在误差降低可接受范围时，这种方法是很实用的。而相似原理是这类方法的核心，推导相似准数是通过数学方法解决这类问题的关键，一般有2种方法：方程分析法和量纲分析法，考虑到最近传输原理课程当中的实际问题，就拿方程分析法来做一个算法总结。（所以相似原理的说明我就直接跳过了，和几何学中的原理类似）

给定一个方程，如v=l/t

在第一个系统里，为：v’=l’/t’

在第二个系统里，为：v”=l”/t”

这时，对同名物理量提出他们的相似常数（常数取值由自己给定），针对上式，分别为Cv、Cl、Ct

则用相似常数、各系统中方程，表示第二个系统中的方程。这里示例为：

v”=Cvv’,l”=Cll’,t”=Ct*t’

Cvv’=Cll’/Ct*t’

上式对比v’=l’/t’，并使之相等，可以得出相似条件：Cv=Cl/Ct，即Cv*Ct/Cl=1

因此，由于相似常数的特定关系为定值（常数），故在相似时，所对应的同名物理量也满足这个关系式。在这个例子当中为：

v’t’/l’=v”t”/l”

在某一类模型当中总是存在这样的相似关系，从而这个关系式（数，即相似准数）可视为这一模型中的一个影响因子（自变量），故可以给予一个物理学意义。例如相似时的雷诺数、努塞尔数都可以通过这种方法确定。

方法总结：在给定关系式的前提下，分别以2个相似系统，引入对应的相似常数并利用消元法去确定相似关系式（相似准数）。其推导过程单纯遵循数学推导

  • 相似转换-量纲分析法

量纲分析法也是相似变换当中常用的方法，其最重要的支持理论是相似第三定律（π定理）：现象的各量（影响因素）之间关系可表示为相似准数之间的函数关系

F（π1，π2，π3…，πn）=0

该方程称为准数方程。

由此，从现象的n个影响因素（物理量）中选取m个相互独立的物理量，则其余n-m个物理量均可由这m个物理量导出，其具体方法则为量纲分析法，是基于对各量（物理量）消除量纲（即分式当中，同名物理量的幂比值为1）而计算的。举这样一个例子：

一个现象的影响因素有：速度、长度、密度、压力、粘性系数、重力、时间，则他们的一般关系式为

f(v,l,p,P,u,g,t)=0

这里取出其中的v，l，p，去导出其他物理量之间的关系系数（即π，相似准数），此时准数方程应为

F（π1，π2，π3，π4）=0

根据剩余的物理量可以分别带入推出：

π[1]=P/v^a1l^b1p^c1

π[2]=u/v^a2l^b2p^c2

π[3]=g/v^a3l^b3p^c3

π[4]=t/v^a4l^b4p^c4

因为相似准数是无量纲的数（常数），从而以π[1]为例，有

[π1]=[ML^-1T^2]/([LT^-1]^a1)( [L]^b1)( [ML^-3]^c1)

对[M]:1=c1

对[L]:-1=a1+b1-3c1

对[T]:-2=-a1

解得a1=2,b1=0,c1=1，于是有π[1]=(P/pu^2 )=Eu

同理可得其余物理量分别为1/Re,Fr,Ho

则准数方程为F（Ho，Re，Fr，Eu）=0

（在结果上，量纲分析法和方程分析法的结果是相同的，原因是其计算过程中的主要目的相同——消去单位；同时，由π定理，准数方程中独立的相似准数个数=n-m。但是实际计算需要各物理量的量纲表达式！！）

常见物理量的量纲表达式：https://wenku.baidu.com/view/c43ca0f6f705cc17552709f3.html