计算传热学课程设计报告-对空气在有泡沫金属介质管内流动与传热的研究.doc

计算传热学课程设计计算传热学基础对空气在有泡沫金属介质管内流动与传热的研究热能与动力工程系10-1班张皓威 10123113雒飞 10123112陈诚 10123115白代立 10123122指导教师： 黄善波 巩亮 徐会金2013年7月目录题目1一、问题分析3二、解题过程4（一）对各个模型的流动和换热进行无量纲化41、对各个模型的换热进行无量纲化42、对各个模型的流动进行无量纲化5（二）Darcy模型的温度分布7（三）Brinkman模型的速度分布和温度分布9（四）Forchheimer模型的速度和温度分布13（五）Brinkman模型和Forchheimer模型的速度分布和温度分布进行对比。17（六）Brinkman模型和Forchheimer模型的f, fRe, Nu值18 （七）总结20附录21附录1 计算Darcy模型的温度程序21附录2 计算Brinkman模型的速度和温度及Nu程序23附录3 计算Forchheimer模型的速度和温度及Nu程序26附录4 计算f和fRe的程序30参考文献31问题三十三（难度：5.0） 一根完全填充多孔介质管外表面为恒热流边界条件（），管内径为，平均流速为的空气在管内流动，其内部层流充分发展流动模型通常有Darcy模型、Brinkman模型和forchheimer模型，管内填充孔隙率为的多孔介质，渗透率表示为：惯性系数表示为：有效导热系数表示为：充分发展的Darcy流动模型： （1）充分发展的Brinkman模型： （2）充分发展的forchheimer模型： （3）质量守恒方程：动量方程： （4）式中，J=u/u是沿坐标轴方向的单位速度矢量；和f分别为流体的密度和动力粘度；V为速度矢量；K，CF分别为渗透率、孔隙率和惯性系数。动量方程右边的4项分别为：压力梯度、Brinkman项、Darcy项和forchheimer项。能量守恒方程： （5）参数：试通过数值迭代计算：1. 将该问题的流动和换热进行无量纲化；2. Darcy模型的温度分布；3. Brinkman模型的速度分布和温度分布，计算该模型下的f、fRe、Nu数；4. forchheimer模型的速度分布和温度分布，计算该模型下的f、fRe、Nu数，并与Brinkman模型的速度分布和温度分布进行对比。组员分工：1、资料搜集：白代立2、能量方程无量纲化离散：张皓威3、Brinkman模型无量纲化离散：雒飞4、Forchheimer模型无量纲化离散：陈诚、白代立5、编程：雒飞、陈诚、张皓威6、word录入：白代立 陈诚7、ppt总结：张皓威一 问题分析（1）分析动量方程 题目的研究对象是充满泡沫金属管内充分发展段的空气，边界条件为恒热流密度，求在不同模型下的速度和温度分布。由于在垂直于轴向的圆形截面内，整个速度和温度分布呈轴对称，所以问题就变成求在圆形横截面内径向方向上的速度和温度分布。 由于流动过程为充分发展段，所以动量方程（4）等号左侧项等于0，可简化为： 根据上式可以得到不同模型的动量方程，即上面给到的式（1）、（2）、（3），所以在求每个模型的速度分布时可以直接对（1）、（2）、（3）式进行处理即可。（2） 分析温度的边界条件0 1 2 N-2 N-1 R N 用方法A在径向从截面圆心到外边界取点如上图所示。由元体能量法处理第二类边界条件(无内热源)二 解题过程（一） 对各个模型的流动和换热进行无量纲化 1对各个模型的换热进行无量纲化： 由于三种模型的能量方程均一样，温度解均用同一能量方程进行编程计算，唯一不同的是：在darcy模型中U=um=1,在Brinkman和Forchheimer模型中计算各点温度需要将已计算出的各点速度带入U即可，均为无量纲量。对能量守恒方程（5）处理： 其中温度的偏导为轴向的，可有以下得出： 在管子轴向取微元段： TT+dT2Hdx r X能量守恒得： 所以: 所以（5）式可化为： (6)定义无量纲： 带入（6）式得：化简得无量纲式如下： 2对各个模型的流动进行无量纲化（1）Darcy模型由于该模型下的速度为平均速度，所以流动模型无需进行无量纲化。（2）Brinkman模型 动量方程的处理： 该模型中的压降部分和Darcy模型中的压降相等，由式(1)可得:将上式带入到（2）式中去得： Brinkman模型无量纲化定义无量纲： u=Uum r=RH 代入得整理得该模型下的无量纲方程（3）Forchheimer模型 动量方程的处理： 该模型中的压降与前两种的有所不同，根据下式计算：将上式带入到（3）式得：整理得：定义无量纲： R=r/H 代入得 整理得无量纲式子：（二） Darcy模型的温度分布 无量纲式如下：离散后得到：由于该模型中的速度为已知的平均速度，取U=Um=1 整理得： （7） 边界条件：，（由第二类边界条件可得）由（7）式和边界条件可得计算方程组： 。 。 。 。 计算程序在附录1给出。取距离壁面9H/10处的点进行独立性检验： 所以取60段，61个节点。 由程序得到的结果如下图所示: 在图1中，R表示径向，由图可知，流体和固体温度均随离开壁面距离的增大而迅速降低,且该数值和文献中的的相近，证明本程序是正确的。（三） Brinkman模型的速度分布和温度分布 整个过程的流程图如下： 在计算温度的时候需要用到每个点的速度值，所以两个量的计算在一个程序中完成。1、 求速度分布该模型下的无量纲方程为：离散整理得 （8）速度边界条件 由式（8）和边界条件可得速度方程组为： （2） 求温度分布： 无量纲式如下：以下为对其离散过程： 整理得: （9）边界条件： （由第二类边界条件可得）由式（9）和边界条件可得温度方程组为： 计算程序在附录2给出。 取距离壁面9H/10处的点进行独立性检验 所以所以取60段，61个节点。 由程序得到的速度分布如下图所示: 图4为金属泡沫管截面的归一化速度分布曲线, 其中U为截面无量纲速度.，R为无量纲半径。由上图可知金属泡沫管内的速度分布很均匀, 这与光管内的抛物线分布不同. 金属泡沫的存在使边界层大大减薄, 边界层的厚度随着泡沫孔径与泡沫管半径之比d/R 的减小而减薄. 对于d/R 相同的泡沫管, 孔隙率越小,边界层就越薄. 该数值和文献中的的相近，证明本程序是正确的。 由程序得到的温度分布如下图所示: 由上图5可知，随着距离壁面的距离增加，流体和泡沫的无量纲温度和温度都不断减小，而且变化率不断减小。这是因为在离散方程中，温度分布受速度的影响很大，所以温度和速度的变化率相近。另由图可知，孔隙率越小，流体温度越高，而孔径比d/H则基本不影响温度的分布。该数值和文献中的的相近，证明本程序是正确的。（四） Forchheimer模型的速度和温度分布 对于速度分布求解，由于系数中存在，计算时需用到上次的值，故采用雅克比迭代法。对于温度分布求解，与Brinkman模型一样，需要用到每个点的速度值，直接带入最终求得的速度以tdma求解,所以两个量的计算在一个程序中完成。整个计算过程流程图如下：（1）、求速度分布该模型下的无量纲方程为：离散：整理得 源项的处理： 令 用最佳处理方式： 得最佳离散方程为：（10）速度边界条件 由边界条件和式（10）可得到速度方程组为： .（2）求温度分布：无量纲式如下： 以下为对其离散过程： 整理得: （7）边界条件： （由第二类边界条件可得）由式（9）和边界条件可得温度方程组为： .程序在附录3给出A 由程序得到的速度分布如下图所示: 图6为金属泡沫管截面的无量纲速度分布曲线, 其中U为截面无量纲速度.，R为无量纲半径。由上图可知金属泡沫管内的速度分布很均匀, 这与光管内的抛物线分布不同. 金属泡沫的存在使边界层大大减薄, 边界层的厚度随着泡沫孔径与泡沫管半径之比d/R 的减小而减薄. 对于d/R 相同的泡沫管, 孔隙率越小,边界层就越薄. 该数值和文献中的的相近，证明本程序是正确的。 B 由程序得到的温度分布如下图所示: 图7为金属泡沫管截面的无量纲温度分布曲线, 其中为截面无量纲速度，R为无量纲半径。随着距离壁面的距离增加，流体和泡沫的无量纲温度和温度都不断减小，而且变化率不断减小。这是因为在离散方程中，温度分布受速度的影响很大，所以温度和速度的变化率相近。另由图可知，孔隙率越大，流体温度越高，而孔径比d/H则基本不影响温度的分布。该数值和文献中的的相近，证明本程序是正确的。（五） Brinkman模型和Forchheimer模型的速度分布和温度分布进行对比。 A 速度分布的比较 图八为相同参数下两种模型的速度分布比较，由图可知，两种模型下速度分布的规律完全相同，在边界层处几乎重合，由计算结果可知雷诺数相同，在管内的相同位置处，Brinkman模型的边界层相对较厚，因为它没有考虑到流体的惯性力。两者雷诺数均小于临界雷诺数，处于层流区，速度边界层较薄，边界层内速度变化非常剧烈。 B 温度分布的比较 图9 对于两种模型的温度分布，由温度离散方程可知温度分布受速度分布的影响很大，故两者的温度分布也相似（六） Brinkman模型和Forchheimer模型的f, fRe, Nu值f 摩擦因子Re 雷诺数fRe 摩擦因子与雷诺数的乘积Nu 努塞尔数它们之间的关系及计算公式如下：对于Brinkman模型对于forchheimer模型 有以下方式求kse: e=0.339 当时， 当时， 为了算法的通用性，f, fRe, Nu 的值均由程序计算，其中f和fRe的值的计算程序在附录4中给出，而Nu得值需要截面平均温度值，所以将它的计算程序直接写入Brinkman模型和Forchheimer模型的程序中，与速度、温度的值一起输出。由程序得：Brinkman模型： f=66.759,fRe=88888.889Forchheimer模型： f=127.005,fRe=169106.48 孔隙率对努塞尔数有很大的影响，由图可知，随着孔隙率的减小，努塞尔数急剧增大。这是因为孔隙率减小意味着金属泡沫密度越大，固体骨架增多，导热大大增加，故努塞尔数增加。对于brinkman和forchheimer两种模型，由于平均无量纲温度接近，故努塞尔数也接近。（七）总结；本文为恒热流边界条件下，空气在泡沫金属管内的强制对流换热，分别在Darcy模型、Brinkman模型和Forchheimer模型下的数值计算，得出了三种模型下，管子横截面上的速度和温度分布情况。另外还把Brinkman模型和Forchheimer模型的速度和温度进行了分别比较，最后还算出了后两种模型的f、fRe、Nu值，还进一步探究了孔隙率对Nu数的影响。并把该结果和相关文献的结果相比较，证明本文的结果还是相对正确的。通过小组成员间的合作和几位老师的耐心帮助，我们在完成任务的同时，在已学知识的基础上对新知识进行了探索和领悟。在分析问题，解方程，编程，得出结果，分析结果，得出结论的过程中，我们收获很多。在此，还要感谢黄善波老师，巩亮老师和徐会金老师的指导和帮助。附录1：计算Darcy模型的温度程序#include #include #define N 60void tdma(double a,double b,double c,double f,double m);void main( )int i,qw;double mN,TN+1,x,s,Ke,Kf,Ks,d,H;double aN,bN,cN,fN;double K,H;double ke,kf,ks,k,e,RA,RB,RC,RD,tm;double kfe;/空气有效导热系数double kse;/铜的有效导热系数tm=0.36;/无量纲平均温度 e=0.339;x=1/(double)N; /步长qw=500; /热流密度Kf=0.0259; /空气导热Ks=398; /泡沫导热系数s=0.6; /孔隙率Ke=s*Kf+(1-s)*Ks;/有效导热系数d=0.001; /泡沫孔的直径H=0.02; /管直径 H=H/2; /定义管半径Hkf=Kf;ks=Ks;/令kf=0,得到kse kf=0;k=sqrt(2*sqrt(2)-1.25*e*e*e-2*s*sqrt(2)/sqrt(3.14*(3-4*e*sqrt(2)-e);RA=4*k/(2*e*e+3.14*k-3.14*k*e)*ks+(4-2*e*e-3.14*k+3.14*k*e)*kf);RB=(e-2*k)*(e-2*k)/(e-2*k)*e*e*ks+(2*e-4*k-e*e*e+2*k*e*e)*kf); RC=(sqrt(2)-2*e)*(sqrt(2)-2*e)/(2*3.14*k*k*(1-2*e*sqrt(2)*ks+2*kf*(sqrt(2)-2*e-3.14*k*k*(1-2*e*sqrt(2); RD=2*e/(e*e*ks+kf*(4-e*e); ke=sqrt(2)/(2*(RA+RB+RC+RD);kse=ke;/系数矩阵赋值for(i=1;i=N-2;i+) ai=1.0-0.5/i; ci=1.0+0.5/i; bi=2.0; fi=-2*x*x*(kse/Ke); b0=1.0; bN-1=1.0; a0=0.0; aN-1=0.0; c0=1.0; cN-1=0.0; f0=0.0; fN-1=-x*(kse/Ke);tdma(a,b,c,f,m); /调用追赶法for(i=0;iN;i+) Ti=mi;TN=0.0;for(i=0;iN+1;i+) /输出结果printf(T%d=%fn,i,Ti);void tdma(double a,double b,double c,double f,double m)int i;double PN,QN;/追的过程P0=c0/b0;Q0=f0/b0;for(i=1;i=0;i-)mi=Pi*mi+1+Qi;return;附录2 计算Brinkman模型的速度和温度及Nu程序#include #include #define N 60void tdma(double a,double b,double c,double f,double m);void main( )int i,qw;double UN+1,mN,TN+1,x,s,Ke,Kf,Ks,d,H;double aN,bN,cN,fN;double K,H,L;double ke,kf,ks,k,e,RA,RB,RC,RD,tm,Nu;double kse; /铜的有效导热系数tm=0.0; /无量纲平均温度 e=0.339;x=1/(double)N; /步长qw=500; /热流密度Kf=0.0259; /空气导热Ks=398; /泡沫导热系数s=0.6; /孔隙率Ke=s*Kf+(1-s)*Ks;/有效导热系数d=0.001; /泡沫孔的直径H=0.02; /管直径 K=d*d*s*s*s/(150*(1-s)*(1-s); /定义渗透率KH=H/2; /定义管半径HL=x*x*s*H*H/K; /源项整理后的表达式kf=Kf;ks=Ks;/令kf=0,得到kse kf=0;k=sqrt(2*sqrt(2)-1.25*e*e*e-2*s*sqrt(2)/sqrt(3.14*(3-4*e*sqrt(2)-e);RA=4*k/(2*e*e+3.14*k-3.14*k*e)*ks+(4-2*e*e-3.14*k+3.14*k*e)*kf);RB=(e-2*k)*(e-2*k)/(e-2*k)*e*e*ks+(2*e-4*k-e*e*e+2*k*e*e)*kf); RC=(sqrt(2)-2*e)*(sqrt(2)-2*e)/(2*3.14*k*k*(1-2*e*sqrt(2)*ks+2*kf*(sqrt(2)-2*e-3.14*k*k*(1-2*e*sqrt(2); RD=2*e/(e*e*ks+kf*(4-e*e); ke=sqrt(2)/(2*(RA+RB+RC+RD);kse=ke;/系数矩阵赋值for(i=1;i=N-2;i+) ai=1.0-0.5/i; ci=1.0+0.5/i; bi=2+L; fi=L; b0=1.0; bN-1=2+L; a0=0.0; aN-1=1.0-0.5/(N-1); c0=1.0; cN-1=0.0; f0=0.0; fN-1=L; tdma(a,b,c,f,m); /调用追赶法for(i=0;iN;i+)Ui=mi;UN=0.0;for(i=0;iN+1;i+)printf(U%d=%fn,i,Ui);printf(n);for(i=1;i=N-2;i+) ai=1.0-0.5/i; ci=1.0+0.5/i; bi=2.0; fi=-2*Ui*x*x*(kse/Ke); b0=1.0; bN-1=1.0; a0=0.0; aN-1=0.0; c0=1.0; cN-1=0.0; f0=0.0; fN-1=-(kse/Ke)*x;tdma(a,b,c,f,m); /调用追赶法for(i=0;iN;i+) Ti=mi;TN=0.0;for(i=0;iN+1;i+) /输出结果printf(T%d=%fn,i,Ti);/显示 kseprintf(kse=%fn,kse);/计算无量纲平均温度for(i=0;iN+1;i+)tm=tm+Ui*Ti*x;printf(tm=%fn,tm); /计算努塞尔数Nu=-2/(Kf/kse)*tm);printf(Nu=%fn,Nu);void tdma(double a,double b,double c,double f,double m)int i;double PN,QN;/追的过程P0=c0/b0;Q0=f0/b0;for(i=1;i=0;i-)mi=Pi*mi+1+Qi;return;附录3 计算Forchheimer模型的速度和温度及Nu程序#include #include #define N 60void tdma(double a,double b,double c,double f,double m);void main()int i,j,k0;double mN,tN,Um,Ke,Kf,Ks;double aN,bN,cN,fN;double Eps=0.000001,ERR,UN,U0N,sum;double yNN,zN;double kf,k,e,K,RA,RB,RC,RD,tm,Nu;double kse;/铜的有效导热系数/定义参数double H,d,s,p,v,x,cF,a1,b1,c1; tm=0.0; /无量纲平均温度 Um=1; /平均速度 H=0.01; /半径 d=0.001; /空半径 s=0.6; /孔隙率 p=1.205; /空气密度 v=0.0000181; /空气动力粘度 Kf=0.0259;/空气导热系数 Ks=398;/泡沫介质导热系数 Ke=s*Kf+(1-s)*Ks;/有效导热系数 e=0.339;/固定参量 K=d*d*s*s*s/(150*(1-s)*(1-s);/渗透率 kf=Kf;/空气导热系数 k0=0;/迭代次数 /计算参数 x=1/(double)(N-1); /步长 cF=1.75/sqrt(150)/pow(s,2/3); /惯性系数 /令kf=0,得到kse kf=0; k=sqrt(2*sqrt(2)-1.25*e*e*e-2*s*sqrt(2)/sqrt(3.14*(3-4*e*sqrt(2)-e); RA=4*k/(2*e*e+3.14*k-3.14*k*e)*Ks+(4-2*e*e-3.14*k+3.14*k*e)*kf); RB=(e-2*k)*(e-2*k)/(e-2*k)*e*e*Ks+(2*e-4*k-e*e*e+2*k*e*e)*kf); RC=(sqrt(2)-2*e)*(sqrt(2)-2*e)/(2*3.14*k*k*(1-2*e*sqrt(2)*Ks+2*kf*(sqrt(2)-2*e-3.14*k*k*(1-2*e*sqrt(2); RD=2*e/(e*e*Ks+kf*(4-e*e); kse=sqrt(2)/(2*(RA+RB+RC+RD);/雅可比迭代计算速度分布/定义中间变量 a1=K/(s*H*H*x*x); b1=p*cF*sqrt(K)/v; c1=K/(2*s*H*H*x);/定义初始值 for(i=0;iN;i+)U0i=0.85;/定义系数矩阵for(i=0;iN;i+)for(j=0;jN;j+)yij=0.0;for(i=0;iN;i+)yii=1; zi=(1+b1*Um+b1*Um*U0i*U0i)/(1+2*a1+2*b1*Um*U0i);for(i=1;iN-1;i+)yii+1=-(a1+c1/(i*x)/(1+2*a1+2*b1*Um*U0i);for(i=1;iN-1;i+)yii-1=-(a1-c1/(i*x)/(1+2*a1+2*b1*Um*U0i);z0=0;y01=-1;/速度迭代计算loop:for(i=0;iN;i+)sum=0.0;for(j=0;ji;j+)sum=sum+yij*U0j;for(j=i+1;jN;j+)sum=sum+yij*U0j;Ui=zi-sum;/速度误差分析ERR=0.0;for(i=0;iN;i+)ERR=ERR+(Ui-U0i);k0=k0+1;printf(迭代次数k0=%d,Er=%fn,k0,ERR);for(i=0;iEps)for(i=0;iN;i+) U0i=Ui;goto loop;/计算温度/定义系数矩阵for(i=1;i=N-2;i+) ai=1.0-0.5/i; ci=1.0+0.5/i; bi=2.0; fi=-2*x*x*Ui*(kse/Ke); b0=1.0; bN-1=1.0; a0=0.0; aN-1=0.0; c0=1.0; cN-1=0.0; f0=0.0; fN-1=-x*(kse/Ke);/追赶法计算tdma(a,b,c,f,m);/输出/最后一次的速度输出printf(迭代次数k0=%d,Er=%fn,k0,ERR);for(i=0;iN;i+)printf(U%d=%fn,i,Ui); printf(U%d=%fn,N,0);printf(n);/输出温度for(i=0;iN;i+)ti=mi;for(i=0;iN;i+)printf(t%d=%fn,i,ti);printf(t%d=%fn,N,0); /输出Kseprintf(kse=%fn,kse); /计算并输出平均温度 for(i=0;i=N;i+)tm=tm+Ui*ti*x; printf(tm=%fn,tm); /计算并输出努塞尔数Nu=-2/(Kf/kse)*tm);printf(Nu=%fn,Nu);/子程序，追赶法求温度分布void tdma(double a,double b,double c,double f,double m)int i;double PN,QN;P0=c0/b0;Q0=f0/b0;for(i=1;i=0;i-)mi=Pi*mi+1+Qi;return;附录4 计算f和fRe的程序#include#includevoid main()double f,s,d,Re,v,p,K,H,H,fRe,Da,cF,Um;double p1; /无量纲压强 double p2,p3; /轴向压降 p=1.205;/空气密度 v=0.0000181;/ Um=1;d=0.001; /泡沫孔的直径s=0.6; /孔隙率 K=d*d*s*s*s/(150*(1-s)*(1-s); /定义渗透率KcF=1.75/(sqrt(150)*pow(s,2/3);/H=0.02; /管直径 H=H/2; /定义管半径H Re=Um*H*p/v;/Da=K/H/H; p2=v*Um/K; /p3=p2+p*cF*Um*Um/sqrt(K);/ p1=K*p2/Um/v;f=8*p1/Re/Da;fRe=f*Re; printf(输出Brinkman的f,fRe:n);printf(f=%fn,f);print