传热学编程实习何鹏举Word文档格式.docx
- 文档编号:22138936
- 上传时间:2023-02-02
- 格式:DOCX
- 页数:20
- 大小:103.26KB
传热学编程实习何鹏举Word文档格式.docx
《传热学编程实习何鹏举Word文档格式.docx》由会员分享,可在线阅读,更多相关《传热学编程实习何鹏举Word文档格式.docx(20页珍藏版)》请在冰豆网上搜索。
计算区域的离散如图所示,总节点数取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<
stdio.h>
math.h>
#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"
);
\n\t\t\t\t\t\t----何鹏举\n"
\n题目:
补充材料练习题一\n"
已知:
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代表面积*/
\n请首先假定一个温度场的初始分布,即给出各节点的温度初值:
\n"
for(i=0;
i<
N;
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;
N-1;
a[i]=(a[i+1]+a[i-1])/(2+m*m*x*x);
a[N-1]=a[N-2];
cha=0;
for(i=0;
cha=cha+abs(a[i]-b[i]);
cha=cha/N;
/*cha代表每次迭代后与上次迭代各点温度差值的平均值*/
t[i]=a[i]*(t0-t1)+t1;
\n\n经数值离散(一阶精度的向后差分法)计算得肋片的温度分布为:
printf("
%4.2f\t"
t[i]);
\n\n"
getchar();
/*采用二阶精度的元体平衡法数值离散(温度初值还用设定的初场,便于比较)*/
a[i]=b[i];
a[N-1]=a[N-2]/(1+0.5*m*m*x*x);
cha=cha+a[i]-b[i];
\n\n经数值离散(二阶精度的元体平衡法)计算得肋片的温度分布为:
}
2练习题二:
二维稳态导热的数值计算
2.1物理问题
图2示出了一矩形区域,其边长L=W=1,假设区域内无内热源,导热系数为常数,三个边温度为T1=0,一个边温度为T2=1,求该矩形区域内的温度分布。
2.2数学描述
对上述问题的微分方程及其边界条件为:
作为参考,一下给出该问题的解析解:
表2列出了由式(2-3)计算得到的在平面区域内各不同位置的温度指。
表2平面区域内温度分布(用式2-3计算得到)
坐标
温度
T
x
y
2.3数值离散
2.3.1区域离散
区域离散如图2所示,x方向总节点数为N,y方向总节点数为M,区域内任一节点用I,j表示。
2.3.2方程的离散
对于图2中所有的内部节点方程(2-1)都适用,因此可写为:
用I,j节点的二阶中心差分代替上式中的二阶导数,得:
上式整理成迭代形式:
补充四个边界上的第一类边界条件得:
2.4计算程序
#defineN10
#defineM10
chars;
inti,j,l;
floatcha,x,y;
floatt[N][M],a[N][M];
/*打印出题目*/
\t\t\t二维稳态导热问题\t\t"
补充材料练习题二\n"
\n矩形区域,边长L=W=1,假设区域内无内热源,导热系数为常熟,三个边温度为T1=0,一个边温度为T2=1,求该矩形区域内的温度分布。
\n是否要手动对温度场赋予初值?
(Y/N):
"
scanf("
%c"
s);
if(s=='
y'
||s=='
Y'
/*手动赋予温度初场*/
\n请首先假定一个温度场的初始分布,即给出各节点的温度初值(一行一行进行):
for(j=0;
j<
M;
j++)
scanf("
t[i][j]);
else
/*自动赋予温度初场*/
for(j=0;
t[i][j]=0.5;
/*四个边界上的第一类边界条件*/
for(j=0;
t[0][j]=0;
t[M-1][j]=0;
t[i][0]=0;
t[i][N-1]=1;
/*步长计算*/
x=1.0/(N-1);
y=1.0/(M-1);
/*迭代循环*/
a[i][j]=t[i][j];
for(j=1;
M-1;
t[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=cha+abs(a[i][j]-t[i][j]);
cha=cha/(N*M);
/*输出温度分布,其中l控制输出值的排列;
这个结果是按照笛卡尔坐标系下平面从左上角开始依次的*/
\n经数值离散计算的该矩形区域内温度分布为:
l=0;
for(j=M-1;
j>
=0;
j--)
{
printf("
%4.3f"
t[i][j]);
l=l+1;
if(l==N)
{
printf("
l=0;
}
}
/*为了是生成的exe文件结果算的后不会立即退出,方便观看*/
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程序
#defineK11
floata,x,y,Fo,Bi;
floatt[N][K],b[N][K];
\t\t\t一维非稳态导热问题\t\t"
补充材料练习题三\n"
y=1;
/*y代表Δτ*/
x=0.05/(N-1);
a=34.89/(7800*712);
Fo=(a*y)/(x*x);
Bi=233*x/34.89;
\n显示格式条件:
\n1、Fo=%3.1f<
0.5\t"
Fo);
\t2、1-2Fo*Bi-2Fo=%4.2f>
0\n\n"
1-2*Fo*Bi-2*Fo);
/*时刻为零时,赋予初场温度*/
t[i][0]=1000;
/*循环开始,每次计算一个时刻*/
K-1;
b[i][j]=t[i][j];
/*下面对每一个时刻进行迭代求解对应的温度分布,公式按传热学课本P178页公式*/
0.001)
{
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;
cha=cha+abs(t[i][j]-b[i][j]);
这个结果是横轴为x,纵轴为τ的直角坐标下从左上角开始依次的*/
\n经数值离散计算的温度分布为:
for(j=K-1;
if(t[i][j]>
999.99)
%6.1f"
%6.2f"
}
/*为了是生成的exe文件结果算的后不会立即退出,方便观看*/
4总结
热传导问题的数值解法的基本流程如下:
1、建立物理模型
2、建立控制方程及定解条件
3、区域离散化
4、建立其节点的物理量的代数方程
5、设立迭代初场
6、求解代数方程组
7、解的分析
在这个题目的解决分析中,由于书中已给出了大部分步骤,我们只需要看懂并编辑好C程序即可,因此得到了极大的方便。
至此,传热学数值解法全部结束。
何鹏举2010年11月19日
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 传热学 编程 实习