传热学编程实习何鹏举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＋１时层。
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日
升级会员 