传热学编程实习何鹏举.docx
《传热学编程实习何鹏举.docx》由会员分享,可在线阅读,更多相关《传热学编程实习何鹏举.docx(20页珍藏版)》请在冰豆网上搜索。
![传热学编程实习何鹏举.docx](https://file1.bdocx.com/fileroot1/2023-2/4/beb56eec-3f95-46ff-8d02-18be13ca435a/beb56eec-3f95-46ff-8d02-18be13ca435a1.gif)
传热学编程实习何鹏举
传热学数值解法
学号:
**********
姓名:
***
1练习题一:
一维稳态导热的数值计算
1.1物理问题
图1示出了一个等截面直肋,处于温度
=80
的流体中。
肋表面与流体之间的对流换热系数为
肋基处温度
,肋端绝热。
肋片由铝合金制成,其导热系数为
,肋片厚度为
,高度为H=0.1m。
试计算肋内的温度分布及肋的总换热量。
1.2数学描述及其解析解
引入无量纲过余温度
,则无量纲温度
描述的肋片导热微分方程及其边界条件如下:
其中(其中符号含义与教科书一致)
上述数学模型的解析解为:
按(1-4)计算得到的在肋内的各点温度有表1给出。
表一等截面直肋内各点温度
坐标x(m)
温度t(
)
1.3数值离散
1.3.1区域离散
在对方程(1-1)~(1-3)进行数值离散之前,应首先进行计算区域的离散。
计算区域的离散如图所示,总节点数取N。
1.3.2微分方程的离散
由于方程(1-1)在计算区域内部处处成立,因而对图1所示的各离散点亦成立。
对任一借点i有:
用θ在节点i的二阶差分代替θ在节点i的二阶导数,得:
整理成迭代形式:
1.3.3边界条件离散
上面得到的离散方程式(1-5),对所有内部节点都成立,因此每个内部节点都可得出一个类似的方程。
事实上,式(1-5)表达的是一个代数方程组。
但这个方程组的个数少于未知数θi(i=1,2,……,N)的个数。
因此,还需要根据边界条件补充两个方程后代数方程组才封闭。
左边界(x=0)为第一类边界条件,温度为已知,因此可以根据式(1-2)直接补充一个方程为:
右边界为第二类边界条件,有图1中边界节点N的向后差分来代替式(1-3)中的导数,得:
将此式整理为迭代形式,得:
1.3.4最终离散格式
1.3.5代数方程组的求解及其程序
方法:
高斯-赛德尔迭代方法
实施步骤:
首先假定一个温度场的初始分布,即给出各点的温度初值,将这些初值带入方程组(1-6)中进行迭代计算,直至收敛。
C程序内容如下:
#include
#include
#defineN11
main()
{
inti;
floatcha;/*cha含义下面用到时会提到*/
floatt[N],a[N],b[N];
floath,t1,t0,r,D,H,x,m,A,p;/*r代表λ,x代表Δx,D代表δ*/
printf("\t\t\t一维稳态导热问题\t\t");
printf("\n\t\t\t\t\t\t----何鹏举\n");
printf("\n题目:
补充材料练习题一\n");
printf("已知:
h=45,t1=80,t0=200,r=110,D=0.01,H=0.1(ISO)\n");
/*下面根据题目赋值*/
h=45.0;t1=80.0;t0=300.0;r=110.0;D=0.01;H=0.1;
x=H/(N-1);A=3.1415926*D*D/4;p=3.1415926*D;m=sqrt((h*p)/(r*A));
/*x代表步长,p代表周长,A代表面积*/
printf("\n请首先假定一个温度场的初始分布,即给出各节点的温度初值:
\n");
for(i=0;i{
scanf("%f",&t[i]);
a[i]=(t[i]-t1)/(t0-t1);
b[i]=a[i];/*这里b[i]用记录一下a[i],后面迭代条件及二阶采用温度初场要用到*/
}
/*采用一阶精度的向后差分法数值离散*/
cha=1;
while(cha>0.0001)
{
a[0]=1;
for(i=1;ia[i]=(a[i+1]+a[i-1])/(2+m*m*x*x);
a[N-1]=a[N-2];
cha=0;
for(i=0;icha=cha+abs(a[i]-b[i]);
cha=cha/N;/*cha代表每次迭代后与上次迭代各点温度差值的平均值*/
}
for(i=0;it[i]=a[i]*(t0-t1)+t1;
printf("\n\n经数值离散(一阶精度的向后差分法)计算得肋片的温度分布为:
\n");
for(i=0;iprintf("%4.2f\t",t[i]);
printf("\n\n");
getchar();
/*采用二阶精度的元体平衡法数值离散(温度初值还用设定的初场,便于比较)*/
for(i=0;ia[i]=b[i];
cha=1;
while(cha>0.0001)
{
a[0]=1;
for(i=1;ia[i]=(a[i+1]+a[i-1])/(2+m*m*x*x);
a[N-1]=a[N-2]/(1+0.5*m*m*x*x);
cha=0;
for(i=0;icha=cha+a[i]-b[i];
cha=cha/N;
}
for(i=0;it[i]=a[i]*(t0-t1)+t1;
printf("\n\n经数值离散(二阶精度的元体平衡法)计算得肋片的温度分布为:
\n");
for(i=0;iprintf("%4.2f\t",t[i]);
printf("\n\n");
getchar();
}
2练习题二:
二维稳态导热的数值计算
2.1物理问题
图2示出了一矩形区域,其边长L=W=1,假设区域内无内热源,导热系数为常数,三个边温度为T1=0,一个边温度为T2=1,求该矩形区域内的温度分布。
2.2数学描述
对上述问题的微分方程及其边界条件为:
作为参考,一下给出该问题的解析解:
表2列出了由式(2-3)计算得到的在平面区域内各不同位置的温度指。
表2平面区域内温度分布(用式2-3计算得到)
坐标
温度
T
坐标
温度
T
x
y
x
y
2.3数值离散
2.3.1区域离散
区域离散如图2所示,x方向总节点数为N,y方向总节点数为M,区域内任一节点用I,j表示。
2.3.2方程的离散
对于图2中所有的内部节点方程(2-1)都适用,因此可写为:
用I,j节点的二阶中心差分代替上式中的二阶导数,得:
上式整理成迭代形式:
补充四个边界上的第一类边界条件得:
2.4计算程序
#include
#include
#defineN10
#defineM10
main()
{
chars;
inti,j,l;
floatcha,x,y;
floatt[N][M],a[N][M];
/*打印出题目*/
printf("\t\t\t二维稳态导热问题\t\t");
printf("\n\t\t\t\t\t\t----何鹏举\n");
printf("\n题目:
补充材料练习题二\n");
printf("\n矩形区域,边长L=W=1,假设区域内无内热源,导热系数为常熟,三个边温度为T1=0,一个边温度为T2=1,求该矩形区域内的温度分布。
\n");
printf("\n是否要手动对温度场赋予初值?
(Y/N):
");
scanf("%c",&s);
if(s=='y'||s=='Y')
/*手动赋予温度初场*/
{
printf("\n请首先假定一个温度场的初始分布,即给出各节点的温度初值(一行一行进行):
\n");
for(i=0;ifor(j=0;jscanf("%f",&t[i][j]);
}
else
/*自动赋予温度初场*/
{
for(i=0;ifor(j=0;jt[i][j]=0.5;
}
/*四个边界上的第一类边界条件*/
for(j=0;j{
t[0][j]=0;
t[M-1][j]=0;
}
for(i=0;i{
t[i][0]=0;
t[i][N-1]=1;
}
/*步长计算*/
x=1.0/(N-1);
y=1.0/(M-1);
/*迭代循环*/
cha=1;
while(cha>0.0001)
{
for(i=0;ifor(j=0;ja[i][j]=t[i][j];
for(i=1;ifor(j=1;jt[i][j]=0.5*y*y*(t[i+1][j]+t[i-1][j])/(x*x+y*y)+0.5*x*x*(t[i][j+1]+t[i][j-1])/(x*x+y*y);
cha=0;
for(i=0;ifor(j=0;jcha=cha+abs(a[i][j]-t[i][j]);
cha=cha/(N*M);
}
/*输出温度分布,其中l控制输出值的排列;这个结果是按照笛卡尔坐标系下平面从左上角开始依次的*/
printf("\n经数值离散计算的该矩形区域内温度分布为:
\n");
l=0;
for(j=M-1;j>=0;j--)
for(i=0;i{
printf("%4.3f",t[i][j]);
l=l+1;
if(l==N)
{
printf("\n");
l=0;
}
}
/*为了是生成的exe文件结果算的后不会立即退出,方便观看*/
getchar();getchar();
/*其中第一个getchar读取了回车键,第二个getchar读取任意键*/
}
3练习题三:
一维非稳态导热的数值计算
非稳态导热问题由于有时间变量,其数值计算出现了一些新的特点。
在非稳态导热微分方程中,与时间因素相关的非稳态项是温度对时间的一阶导数,这给差分离散带来了新的特点。
由于这个特点,可以采用不同的方法构造差分方程,从而得到几种不同的差分格式,即所谓的显式、隐式和半显式。
我们仍从一个具体问题出发来研究非稳态导热问题的数值计算。
3.1问题
一块无限大平板(如图3所示),其一半厚度为L=0.1m,初始温度T0=1000
,突然将其插入温度T
=20
的流体介质中。
平板的导热系数
=34.89W/(m
),密度
=7800kg/m3,比热c=712J/(kg
),平板与介质的对流换热系数为h=233W/(m2
),求平板内各点的温度分布。
3.2数学描述
由于平板换热关于中心线是对称的,仅对平板一半区域进行计算即可。
坐标x的原点选在平板中心线上,因而一半区域的非稳态导热的数学描述为:
该数学模型的解析解为:
其中,为方程的根,。
表3给出了在平板表面(x=L)处由式(3-5)计算得到的不同时刻的温度值。
表3平板表面各不同时刻温度值
时间
(S)
1
2
3
4
5
6
7
8
9
10
温度
(
)
3.3数值离散
3.3.1计算区域的离散
一维非稳态导热指的是空间坐标是一维的。
若考虑时间坐标,则所谓的一维非稳态导热实际上是二维问题(见图4),即:
有时间坐标T和空间坐标X两个变量。
但要注意,时间坐标是单向的,就是说,前一时刻的状态会对后一时刻的状态有影响,但后一时刻的状态却不能影响到前一时刻,图4示出了以X和T为坐标的计算区域的离散,时间从T=0开始,经过一个个时层增加到K时层和K+1时层。
3.3.2微分方程的离散
对于I节点,在K和K+1时刻可将微分方程(3-1)写成下面式子:
将式(3-6)、(3-7)的左端温度对时间的偏导数进行差分离散为:
观察是(3-8)和(3-9),这两个式子的左端差分式完全相同,但在两个式子中却有不同含义。
对式(3-8),右端项相对i点在K时刻的导数是向前差分。
而在式(3-9)中,右端项是i点在K+1时刻的导数的向后差分。
将式(3-8)和(3-9)分别代入(3-6)和(3-7),并将式(3-6)和(3-7)右端关于x的二阶导数用相应的差分代替,则可得到下列显式和隐式两种不同的差分格式:
显式:
全隐式:
从式(3-10)可见,其右端只涉及K时刻的温度,当从K=0(即
=0时刻)开始计算时,在K=0时等号右端都是已知值,因而直接可计算出K=1时刻各点的温度。
由K=1时刻的各点的温度值,又可以直接利用是(3-10)计算K=2时刻的各点的温度值,这样一个时层一个时层的往下推,各时层的温度都能利用事(3-10)直接计算出来,不要求解代数方程组。
而对于式(3-11)等号右端包含了与等号左端同一时刻但不同节点的温度,因而必须通过求解代数方程组才能求得这些节点的温度值。
3.3.3边界条件的离散
对于式(3-3)和(3-4)所给出的边界条件,可以直接用差分代替微分,也可以用元体平衡法给出相应的边界条件,亦有显式和隐式之分。
通常,当内部节点采用显式时,边界节点也用显式离散;内部节点用隐式时,边界节点亦用隐式。
边界节点的差分格式是显式还是隐式,取决于如何与内部节点的差分方程组合。
用K+1时刻相应节点的差分,代替式(3-3)和(3-4)中的微分,可得到边界节点的差分方程:
3.3.4最终离散格式
显式:
其中K=0,1,2……。
当采用二阶精度的元体平衡法离散时,与式(3-15)和(3-16)对应的离散格式为:
其中。
隐式:
其中K=0,1,2……
在用隐式差分计算时,每个时层都需要迭代求解代数方程组(3-17)~(3-20)。
在每个时层计算时,都要首先假定一个温度场(一般取上一时层的温度场为本时层的初始场),然后迭代计算直至收敛。
3.4程序
#include
#include
#defineN10
#defineK11
main()
{
inti,j,l;
floatcha;
floata,x,y,Fo,Bi;
floatt[N][K],b[N][K];
/*打印出题目*/
printf("\t\t\t一维非稳态导热问题\t\t");
printf("\n\t\t\t\t\t\t----何鹏举\n");
printf("\n题目:
补充材料练习题三\n");
y=1;/*y代表Δτ*/
x=0.05/(N-1);
a=34.89/(7800*712);
Fo=(a*y)/(x*x);
Bi=233*x/34.89;
printf("\n显示格式条件:
");
printf("\n1、Fo=%3.1f<0.5\t",Fo);
printf("\t2、1-2Fo*Bi-2Fo=%4.2f>0\n\n",1-2*Fo*Bi-2*Fo);
/*时刻为零时,赋予初场温度*/
for(i=0;it[i][0]=1000;
/*循环开始,每次计算一个时刻*/
for(j=0;j{
for(i=0;ib[i][j]=t[i][j];
/*下面对每一个时刻进行迭代求解对应的温度分布,公式按传热学课本P178页公式*/
cha=1;
while(cha>0.001)
{
for(i=0;i{
if(i==0)
t[i][j+1]=Fo*(t[i+1][j]+t[i+1][j])+(1-2*Fo)*t[i][j];/*当计算t[0]时,要用到t[-1],其中t[-1]=t[2]的(对称分布)*/
else
t[i][j+1]=Fo*(t[i+1][j]+t[i-1][j])+(1-2*Fo)*t[i][j];
t[N-1][j+1]=t[N-2][j]*(1-2*Fo*Bi-2*Fo)+2*Fo*t[N-1][j]+2*Fo*Bi*20;/*边界点温度用热平衡法推导出公式*/
}
cha=0;
for(i=0;icha=cha+abs(t[i][j]-b[i][j]);
cha=cha/N;
}
}
/*输出温度分布,其中l控制输出值的排列;这个结果是横轴为x,纵轴为τ的直角坐标下从左上角开始依次的*/
printf("\n经数值离散计算的温度分布为:
\n");
l=0;
for(j=K-1;j>=0;j--)
for(i=0;i{
if(t[i][j]>999.99)
printf("%6.1f",t[i][j]);
else
printf("%6.2f",t[i][j]);
l=l+1;
if(l==N)
{
printf("\n");
l=0;
}
}
getchar();/*为了是生成的exe文件结果算的后不会立即退出,方便观看*/
}
4总结
热传导问题的数值解法的基本流程如下:
1、建立物理模型
2、建立控制方程及定解条件
3、区域离散化
4、建立其节点的物理量的代数方程
5、设立迭代初场
6、求解代数方程组
7、解的分析
在这个题目的解决分析中,由于书中已给出了大部分步骤,我们只需要看懂并编辑好C程序即可,因此得到了极大的方便。
至此,传热学数值解法全部结束。
何鹏举2010年11月19日