相变传热与流体流动数值分析作业3.docx
《相变传热与流体流动数值分析作业3.docx》由会员分享，可在线阅读，更多相关《相变传热与流体流动数值分析作业3.docx（22页珍藏版）》请在冰豆网上搜索。
相变传热与流体流动数值分析作业3
相变传热与流体流动数值分析
作业3
学院（系）：
能源与动力学院
专业：
能源与动力工程
********
学号：
********
指导教师：
宋永臣教授
完成日期：
2012.3.6
大连理工大学
DalianUniversityofTechnology
TheFiniteVolumeMethodforConvection-DiffusionProblems
Subject:
ApropertyΦistransportedbymeansofconvectionanddiffusionthroughtheone-dimensionaldomainsketchedinFigure1.Thegoverningequationis
;boundaryconditionsareΦ0=1atx=0andΦL=0atx=L.UsingfiveequallyspacedcellsforconvectionanddiffusioncalculatethedistributionofΦafunctionofxforcase:
（i）Case1:
u=0.1m/s;
（ii）Case2:
u=2.5m/s;
（iii）Case3:
u=2.5m/swith20gridnodes;
Thefollowingdataapply:
LengthL=1.0m,ρ=1.0kg/m3,Γ=0.1kg/m/s.
Solution:
（I）Thecentraldifferencingscheme:
//王佳琪-作业-中心差分.cpp:
定义控制台应用程序的入口点。
//
#include
#include
#include
#include
#include
#include
#include
#include
#include
#defineN5
usingnamespacestd;
inti;
doubleaw[N],b[N],ae[N],f[N],x[N];
/*-------追赶法求解数组-------*/
voidtdma（）
{
doublel[N],u[N],y[N];
for（i=1;i{
u[0]=b[0];l[0]=0;
l[i]=aw[i]/u[i-1];
u[i]=b[i]-l[i]*ae[i-1];
}
y[0]=f[0];
for（i=1;iy[i]=f[i]-l[i]*y[i-1];
x[N-1]=y[N-1]/u[N-1];
for（i=N-2;i>=0;i--）
x[i]=（y[i]-ae[i]*x[i+1]）/u[i];
}
voidmain（）
{
voidOutput（）;
/*---------定义变量及边界条件---------*/
doubleF,D,u,ρ,Γ,x,L,φA,φB,Sp[N];
u=0.1;ρ=1;Γ=0.1;L=1;φA=1;φB=0;
x=L/N;F=ρ*u;D=Γ/x;
/*---------网格离散---------*/
for（i=0;i{
if（i==0）
{
aw[i]=0;ae[i]=-（D-F/2）;Sp[i]=-（2*D+F）;f[i]=（2*D+F）*φA;
}
elseif（i==N-1）
{
ae[i]=0;aw[i]=-（D+F/2）;Sp[i]=-（2*D-F）;f[i]=（2*D-F）*φB;
}
else
{
aw[i]=-（D+F/2）;ae[i]=-（D-F/2）;Sp[i]=0;f[i]=0;
}
b[i]=-aw[i]-ae[i]-Sp[i];
}
tdma（）;
Output（）;
}
voidOutput（）
{
/*---------后处理文件生成---------*/
ostringstreamname;
name.str（""）;
name<<"central_.plt";
ofstreamout（name.str（）.c_str（））;
out<<"Title=central"<for（i=0;i{
out<}
}
Results:
（i）case1:
u=0.1m/s;N=5;
（ii）Case2:
u=2.5m/s;N=5;
（iii）case3:
u=2.5m/s;N=20;
Discussion:
Fromthefigure,wecanseetheagreementswiththeanalyticalsolutionareverygoodincase1andcase3.Whileincase2,becauseofthefinitegridwithahighspeed,thenumericalsolutionappearstooscillate,called‘wiggles’.Inoneword,thecentraldifferencingschemedoesn’tadapttothehighpecletnumber（F/D>2）.
（II）Theupwinddifferencingscheme:
//王佳琪-作业3-迎风.cpp:
定义控制台应用程序的入口点。
//
#include"stdafx.h"
#include
#include
#include
#include
#include
#include
#include
#include
#defineN5
usingnamespacestd;
inti;
doubleaw[N],b[N],ae[N],f[N],x[N];
/*-------追赶法求解数组-------*/
voidtdma（）
{
doublel[N],u[N],y[N];
for（i=1;i{
u[0]=b[0];l[0]=0;
l[i]=aw[i]/u[i-1];
u[i]=b[i]-l[i]*ae[i-1];}
y[0]=f[0];
for（i=1;iy[i]=f[i]-l[i]*y[i-1];
x[N-1]=y[N-1]/u[N-1];
for（i=N-2;i>=0;i--）x[i]=（y[i]-ae[i]*x[i+1]）/u[i];
}
voidmain（）
{
voidOutput（）;
doubleFw1,Fw,Fe1,Fe,Dw,De,u,ρ,Γ,x,L,φA,φB,Sp[N];
u=0.1;ρ=1;Γ=0.1;L=1;φA=1;φB=0;
x=L/N;Fw=Fe=ρ*u;Dw=De=Γ/x;
Fw1=0>Fw?
0:
Fw;
Fe1=0>（-Fe）?
0:
（-Fe）;
/*---------网格离散---------*/
for（i=0;i{
if（i==0）
{aw[i]=0;ae[i]=-（De-Fe1）;Sp[i]=-（2*Dw+Fw1）;f[i]=（2*Dw+Fw1）*φA;}
elseif（i==N-1）
{aw[i]=-（Dw+Fw1）;ae[i]=0;Sp[i]=-（2*De+Fe1）;f[i]=（2*De-Fe1）*φB;}
else
{aw[i]=-（Dw+Fw1）;ae[i]=-（De+Fe1）;Sp[i]=0;f[i]=0;}
b[i]=-aw[i]-ae[i]-Sp[i];
}
tdma（）;
Output（）;
}
voidOutput（）
{
ostringstreamname;
name.str（""）;
name<<"upwind.plt";
ofstreamout（name.str（）.c_str（））;
out<<"Title=upwind"<for（i=0;iout<}
Results:
（i）case1:
u=0.1m/s;N=5;
（ii）Case2:
u=2.5m/s;N=5;
Discussion:
Theupwinddifferencingschemeproducesamuchmorerealisticsolutionthanthecentraldifferencingscheme,however,notveryclosetotheexactsolutionnearboundaryB.Thatisbecausethatitsaccuracyisonlyfirstorder,whilethecentraldifferencingissecondorder.
（III）Thehybriddifferencingscheme:
//王佳琪-作业-混合.cpp:
定义控制台应用程序的入口点。
//
#include"stdafx.h"
#include
#include
#include
#include
#include
#include
#include
#include
#defineN5
usingnamespacestd;
inti;
doubleaw[N],b[N],ae[N],f[N],x[N];
/*-------追赶法求解数组-------*/
voidtdma（）
{
doublel[N],u[N],y[N];
for（i=1;i{
u[0]=b[0];l[0]=0;
l[i]=aw[i]/u[i-1];
u[i]=b[i]-l[i]*ae[i-1];}
y[0]=f[0];
for（i=1;iy[i]=f[i]-l[i]*y[i-1];
x[N-1]=y[N-1]/u[N-1];
for（i=N-2;i>=0;i--）x[i]=（y[i]-ae[i]*x[i+1]）/u[i];
}
voidmain（）
{
voidOutput（）;
doubleFw1,Fw,Fe1,Fe,Dw,De,u,ρ,Γ,x,L,φA,φB,Sp[N];
u=2.5;ρ=1;Γ=0.1;L=1;φA=1;φB=0;
x=L/N;Fw=Fe=ρ*u;Dw=De=Γ/x;
if（Fw<0&&（Dw+Fw/2）<0）Fw1=0;
elseif（Fw<0&&（Dw+Fw/2）>0）Fw1=（Dw+Fw/2）;
elseif（Fw>0&&（Dw+Fw/2）<0）Fw1=Fw;
elseif（0elseif（0<（Dw+Fw/2）&&（Dw+Fw/2）if（（-Fe）<0&&（De-Fe/2）<0）Fe1=0;
elseif（（-Fe）<0&&（De-Fe/2）>0）Fe1=（De-Fe/2）;
elseif（（-Fe）>0&&（De-Fe/2）<0）Fe1=（-Fe）;
elseif（0<（-Fe）&&（-Fe）<（De-Fe/2））Fe1=（De-Fe/2）;
elseif（0<（De-Fe/2）&&（De-Fe/2）<（-Fe））Fe1=（-Fe）;
/*---------网格离散---------*/
for（i=0;i{
if（i==0）
{aw[i]=0;ae[i]=-Fe1;Sp[i]=-（2*Dw+Fw1）;f[i]=（2*Dw+Fw1）*φA;}
elseif（i==N-1）
{aw[i]=-Fw1;ae[i]=0;Sp[i]=-（2*De+Fe1）;f[i]=（2*De-Fe1）*φB;}
else
{aw[i]=-Fw1;ae[i]=-Fe1;Sp[i]=0;f[i]=0;}
b[i]=-aw[i]-ae[i]-Sp[i];
}
tdma（）;
Output（）;
}
voidOutput（）
{doublem=25;
ostringstreamname;
name.str（""）;
name<<"hybrid.plt";
ofstreamout（name.str（）.c_str（））;
out<<"Title=hybrid"<for（i=0;iout<}
Results:
（i）Case1:
u=2.5m/s;N=5;
（ii）Case2:
u=2.5m/s;N=25;
Discussion:
Onacombinationofcentralandupwinddifferencingschemes,thehybridschemeprovidesamuchmorerealisticsolutionandagreestheexactsolutionverywell.Themorenodes,themoreaccuratetheresultis.
（IV）Thepower-lawscheme:
//王佳琪-作业-指数.cpp:
定义控制台应用程序的入口点。
//
#include"stdafx.h"
#include
#include
#include
#include
#include
#include
#include
#include
#defineN5
usingnamespacestd;
inti;
doubleaw[N],b[N],ae[N],f[N],x[N];
/*-------追赶法求解数组-------*/
voidtdma（）
{
doublel[N],u[N],y[N];
for（i=1;i{
u[0]=b[0];l[0]=0;
l[i]=aw[i]/u[i-1];
u[i]=b[i]-l[i]*ae[i-1];}
y[0]=f[0];
for（i=1;iy[i]=f[i]-l[i]*y[i-1];
x[N-1]=y[N-1]/u[N-1];
for（i=N-2;i>=0;i--）x[i]=（y[i]-ae[i]*x[i+1]）/u[i];
}
voidmain（）
{
voidOutput（）;
ofstreamnote;
doubleFw,Fw1,Fe,Fe1,Dw,De,u,ρ,Γ,x,L,φA,φB,Sp[N],p,p1,Pew,Pee,q,q1;
u=2.5;ρ=1;Γ=0.1;L=1;φA=1;φB=0;
x=L/N;
Fw=Fe=ρ*u;
Dw=De=Γ/x;
Pew=Fw/Dw;
Pee=Fe/De;
Fw1=0>Fw?
0:
Fw;
Fe1=0>（-Fe）?
0:
（-Fe）;
p=pow（（1-0.1*fabs（Pew））,5）;
p1=0>p?
0:
p;
q=pow（（1-0.1*fabs（Pee））,5）;
q1=0>q?
0:
q;
/*---------网格离散---------*/
for（i=0;i{
if（i==0）
{
aw[i]=0;ae[i]=-（Fe1+De*q1）;Sp[i]=-（2*Dw+Fw1）;f[i]=（2*Dw+Fw1）*φA;
}
elseif（i==N-1）
{
aw[i]=-（Fw1+Dw*p1）;ae[i]=0;Sp[i]=-（2*De+Fe1）;f[i]=（2*De-Fe1）*φB;
}
else
{
aw[i]=-（Fw1+Dw*p1）;ae[i]=-（Fe1+De*q1）;Sp[i]=0;f[i]=0;
}
b[i]=-aw[i]-ae[i]-Sp[i];
}
tdma（）;
Output（）;
}
voidOutput（）
{doublem=25;
ostringstreamname;
name.str（""）;
name<<"power-law.plt";
ofstreamout（name.str（）.c_str（））;
out<<"power-low"<for（i=0;i{
out<}
}
Result:
（i）Case1:
u=2.5m/s;N=5;
（ii）Case2:
u=2.5m/s;N=25;
Discussion:
Thisresultismoreaccuratethanthethreeformerschemes.
升级会员 