传热学编程实习何鹏举.docx

1、传热学编程实习何鹏举传热学数值解法学号：*姓名：*1练习题一：一维稳态导热的数值计算1.1物理问题图1示出了一个等截面直肋，处于温度=80的流体中。肋表面与流体之间的对流换热系数为,肋基处温度，肋端绝热。肋片由铝合金制成，其导热系数为，肋片厚度为，高度为H=0.1m。试计算肋内的温度分布及肋的总换热量。1.2数学描述及其解析解引入无量纲过余温度，则无量纲温度描述的肋片导热微分方程及其边界条件如下：其中 （其中符号含义与教科书一致）上述数学模型的解析解为：按（1-4）计算得到的在肋内的各点温度有表1给出。表一 等截面直肋内各点温度坐标x（m）温度t（）1.3数值离散1.3.1区域离散在对方程（1

2、-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）为第一类边界条件，温度为已知，因此

3、可以根据式（1-2）直接补充一个方程为：右边界为第二类边界条件，有图1中边界节点N的向后差分来代替式（1-3）中的导数，得：将此式整理为迭代形式，得：1.3.4最终离散格式1.3.5代数方程组的求解及其程序方法：高斯-赛德尔迭代方法实施步骤：首先假定一个温度场的初始分布，即给出各点的温度初值，将这些初值带入方程组（1-6）中进行迭代计算，直至收敛。C程序内容如下：#include#include#define N 11main() int i; float cha;/*cha含义下面用到时会提到*/ float tN,aN,bN; float h,t1,t0,r,D,H,x,m,A,p; /*

4、r代表,x代表x，D代表*/ printf(ttt一维稳态导热问题tt); printf(ntttttt-何鹏举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(

5、n请首先假定一个温度场的初始分布，即给出各节点的温度初值：n); for(i=0;i0.0001) a0=1; for(i=1;iN-1;i+) ai=(ai+1+ai-1)/(2+m*m*x*x); aN-1=aN-2; cha=0; for(i=0;iN;i+) cha=cha+abs(ai-bi); cha=cha/N;/*cha代表每次迭代后与上次迭代各点温度差值的平均值*/ for(i=0;iN;i+) ti=ai*(t0-t1)+t1; printf(nn经数值离散(一阶精度的向后差分法)计算得肋片的温度分布为:n); for(i=0;iN;i+) printf(%4.2ft,ti

6、); printf(nn); getchar(); /*采用二阶精度的元体平衡法数值离散(温度初值还用设定的初场，便于比较)*/ for(i=0;i0.0001) a0=1; for(i=1;iN-1;i+) ai=(ai+1+ai-1)/(2+m*m*x*x); aN-1=aN-2/(1+0.5*m*m*x*x); cha=0; for(i=0;iN;i+) cha=cha+ai-bi; cha=cha/N; for(i=0;iN;i+) ti=ai*(t0-t1)+t1; printf(nn经数值离散(二阶精度的元体平衡法)计算得肋片的温度分布为:n); for(i=0;iN;i+) pr

7、intf(%4.2ft,ti); printf(nn); getchar();2练习题二：二维稳态导热的数值计算2.1物理问题图2示出了一矩形区域，其边长L=W=1，假设区域内无内热源，导热系数为常数，三个边温度为T1=0，一个边温度为T2=1，求该矩形区域内的温度分布。2.2 数学描述对上述问题的微分方程及其边界条件为：作为参考，一下给出该问题的解析解：表2列出了由式（2-3）计算得到的在平面区域内各不同位置的温度指。表2 平面区域内温度分布（用式2-3计算得到）坐标温度T坐标温度Txyxy2.3数值离散2.3.1区域离散 区域离散如图2所示，x方向总节点数为N，y方向总节点数为M，区域内任

8、一节点用I,j表示。2.3.2方程的离散对于图2中所有的内部节点方程（2-1）都适用，因此可写为：用I,j节点的二阶中心差分代替上式中的二阶导数，得：上式整理成迭代形式：补充四个边界上的第一类边界条件得：2.4计算程序#include#include#define N 10#define M 10main() char s; int i,j,l; float cha,x,y; float tNM,aNM; /*打印出题目*/ printf(ttt二维稳态导热问题tt); printf(ntttttt-何鹏举n); printf(n题目：补充材料练习题二n); printf(n矩形区域，边长L=

9、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;iN;i+) for(j=0;jM;j+) scanf(%f,&tij); else /*自动赋予温度初场*/ for(i=0;iN;i+) for(j=0;jM;j+) tij=0.5; /*四个边界上的第一类边界条件*/

10、for(j=0;jM;j+) t0j=0; tM-1j=0; for(i=0;i0.0001) for(i=0;iN;i+) for(j=0;jM;j+) aij=tij; for(i=1;iN-1;i+) for(j=1;jM-1;j+) tij=0.5*y*y*(ti+1j+ti-1j)/(x*x+y*y)+0.5*x*x*(tij+1+tij-1)/(x*x+y*y); cha=0; for(i=0;iN;i+) for(j=0;j=0;j-) for(i=0;iN;i+) printf(%4.3f ,tij); l=l+1; if(l=N) printf(n); l=0; /*为了是生

11、成的exe文件结果算的后不会立即退出，方便观看*/ getchar();getchar(); /*其中第一个getchar读取了回车键，第二个getchar读取任意键*/3练习题三：一维非稳态导热的数值计算 非稳态导热问题由于有时间变量，其数值计算出现了一些新的特点。在非稳态导热微分方程中，与时间因素相关的非稳态项是温度对时间的一阶导数，这给差分离散带来了新的特点。由于这个特点，可以采用不同的方法构造差分方程，从而得到几种不同的差分格式，即所谓的显式、隐式和半显式。我们仍从一个具体问题出发来研究非稳态导热问题的数值计算。3.1问题一块无限大平板（如图3所示），其一半厚度为L=0.1m，初始温度

12、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）12345678910温度（）3.3数值离散3.3.1计算区域的离散一

13、维非稳态导热指的是空间坐标是一维的。若考虑时间坐标，则所谓的一维非稳态导热实际上是二维问题（见图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），这两个式子的左端差分式完全相同，但在两个式子中却有不同含义。对式（

14、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）直

15、接计算出来，不要求解代数方程组。而对于式（3-11）等号右端包含了与等号左端同一时刻但不同节点的温度，因而必须通过求解代数方程组才能求得这些节点的温度值。3.3.3边界条件的离散对于式（3-3）和（3-4）所给出的边界条件，可以直接用差分代替微分，也可以用元体平衡法给出相应的边界条件，亦有显式和隐式之分。通常，当内部节点采用显式时，边界节点也用显式离散；内部节点用隐式时，边界节点亦用隐式。边界节点的差分格式是显式还是隐式，取决于如何与内部节点的差分方程组合。用K+1时刻相应节点的差分，代替式（3-3）和（3-4）中的微分，可得到边界节点的差分方程：3.3.4最终离散格式显式：其中K=0，1，2

16、。当采用二阶精度的元体平衡法离散时，与式（3-15）和（3-16）对应的离散格式为：其中 。隐式：其中K=0，1，2在用隐式差分计算时，每个时层都需要迭代求解代数方程组（3-17）（3-20）。在每个时层计算时，都要首先假定一个温度场（一般取上一时层的温度场为本时层的初始场），然后迭代计算直至收敛。3.4程序#include#include#define N 10#define K 11main() int i,j,l; float cha; float a,x,y,Fo,Bi; float tNK,bNK; /*打印出题目*/ printf(ttt一维非稳态导热问题tt); printf(n

17、tttttt-何鹏举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.1f0nn,1-2*Fo*Bi-2*Fo); /*时刻为零时，赋予初场温度*/ for(i=0;iN;i+) ti0=1000; /*循环开始，每次计算一个时刻*/ for(j=0;jK-1;j+) for(i=0;i0.001) for(i=0;iN-1;i+) if(i=0) tij+1=Fo*(

18、ti+1j+ti+1j)+(1-2*Fo)*tij;/*当计算t0时，要用到t-1,其中t-1=t2的(对称分布)*/ else tij+1=Fo*(ti+1j+ti-1j)+(1-2*Fo)*tij; tN-1j+1=tN-2j*(1-2*Fo*Bi-2*Fo)+2*Fo*tN-1j+2*Fo*Bi*20;/*边界点温度用热平衡法推导出公式*/ cha=0; for(i=0;i=0;j-) for(i=0;i999.99) printf(%6.1f ,tij); else printf(%6.2f ,tij); l=l+1; if(l=N) printf(n); l=0; getchar();/*为了是生成的exe文件结果算的后不会立即退出，方便观看*/4总结热传导问题的数值解法的基本流程如下：1、建立物理模型2、建立控制方程及定解条件3、区域离散化4、建立其节点的物理量的代数方程5、设立迭代初场6、求解代数方程组7、解的分析在这个题目的解决分析中，由于书中已给出了大部分步骤，我们只需要看懂并编辑好C程序即可，因此得到了极大的方便。至此，传热学数值解法全部结束。何鹏举 2010年11月19日