一维导热计算

非稳态

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

稳态

参考二维稳态导热

二维导热计算

稳态

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

非稳态-显式

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

非稳态-隐式

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