5.一厚度为250mm无限大平壁,其导热系数λ=43+0.08tw/(m·k),平壁一侧温度为250℃,另一侧温度为46℃,试用数值方法确定平壁内的温度分布,并确定通过该平壁的热流密度。解:(1)建立模型本题属于非常物性,无热源的一维稳态导热过程将平壁沿厚度方向(x方向)划分为N个均匀相等的间距。节点布置如图所示本题给出N=10。(2)通过热流法建立离散方程a)内节点离散方程对节点P(i)所代表的微元体,在x方向上与节点P相邻的节点分别为L(i-1)和R(i+1)。由于节点之间的间距很小,可以认为相邻节点间的温度分布是线性的。节点P所代表的网格单元与它周围各网格单元之间的导热量可根据傅里叶定律直接写为:11iiLPLPttxl--Φ=×Δ11iiRPRPttxl+-Φ=×Δ其中1()430.082iiRPttl++=+×1()430.082iiLPttl-+=+×对节点P(i)所代表的微元写热平衡式,即可得节点P(i)温度的离散方程0LPRPΦ+Φ=11110iiiiLPRPttttxxll-+--×+×=ΔΔ11LPiLPiiLPRPtttllll-++=+b)边界节点离散方程由于本题的壁面温度属于第一类边界条件,因此1250t=146nt+=c)热流密度的计算11()nRPiiittqtlldd+=-=Δ=∑N+1N4321RPL(3)编写流程图计算机程序中使用的变量标识符:i节点的坐标变量t(i)节点i的温度tt前一次算出的节点温度k迭代次数计算机程序中输入数据:n沿x方向的网格划分数e控制迭代过程终止的误差m允许的最大迭代次数开始N=10,t1=250,tn+1=46,ti=1,e=0.01,k=10000输入n,e,k,t1,tn+1,ti迭代次数k=0tti=tik=k+1km打印m,ti,q打印“不收敛”停机NoNoYesYes(4)编写程序本题使用matlab软件,所编写的程序如下:clear;clc;t=ones(1,11);%设定各项初始值q=ones(1,11);t(1)=250;t(11)=46;e=0.01;k=1;while1%迭代程序tt=t;fori=2:10a=43+0.08*(t(i-1)+t(i))/2;b=43+0.08*(t(i)+t(i+1))/2;t(i)=(a*t(i-1)+b*t(i+1))/(a+b);endk=k+1;ttt=abs(t(5)-tt(5));if(ttte)break;endendfori=2:11%计算热流密度a=43+0.08*(t(i-1)+t(i))/2;q(i)=q(i-1)+a*(t(i)-t(i-1));endktq=q(i)/0.25w=0:25:250;v=w/25+1;y=t(v);%绘制温度分布图plot(w,y),xlabel('位置(mm)'),ylabel('温度(℃)')legend('平壁内的温度分布',0),grid(5)计算结果迭代次数k=75平壁温度分布见下表节点i1234567891011温度t(℃)250.00232.01213.59194.71175.35155.46134.98113.8792.0669.4746.00热流密度q=-4.47E+04w/m2(6)温度分布图10.一砖墙厚240mm,内、外表面的表面传热系数分别为6.0w/(m2·k)和15w/(m2·k),墙体材料的导热系数λ=0.43w/(m·k),密度ρ=1668kg/m3,比热容c=0.75kJ/(kg·K),室内空气温度保持不变为20℃,室外空气温度周期性变化,中午12点温度最高为3℃,晚上12点温度最低为-9℃,试用数值计算方法确定内、外墙壁面温度在一天内的变化。解:(1)建立模型本题属于常物性,无热源的一维非稳态导热过程将平壁沿厚度方向(x方向)划分为N个均匀相等的间距,时间从τ=0开始,按Δτ分割为k段,令i表示内节点位置,k表示kΔτ时刻,节点布置如图所示本题给出N=10。(2)通过热流法建立离散方程(采用显式差分格式)a)内节点离散方程导热微分方程式为22ttaxt∂∂=∂∂(1)温度对时间的一阶导数,采用向前差分,则1,kkiiikttttt+-∂⎛⎞=⎜⎟∂Δ⎝⎠(2)温度对x的二阶导数,采用中心差分表达式应为21122,2kkkiiiikttttxx-+⎛⎞-+∂=⎜⎟∂Δ⎝⎠(3)将式(3)和式(2)代入式(1),就得到内节点(i,k)的节点离散方程11122kkkkkiiiiitttttaxt+-+--+=ΔΔ将上式移项整理得()1112212kkkkiiiiaattttxxtt+-+ΔΔ⎛⎞=++-⎜⎟ΔΔ⎝⎠上式可写作()()11112kkkkiiiitFottFot+-+=++-kii=1i=n+1xτ012Fo≤b)边界节点离散方程由于本题的壁面温度属于第三类边界条件,因此()1121112kkkkkkafttttxhttcxlrt+--Δ--=ΔΔ整理上式,可得()()212111112kkkkkkafhxcxttttttrllt+ΔΔ-+-=-Δ上式中aahxBilΔ=,21cxForltΔ=Δ,于是()()12111112kkkkkkafttBittttFo+-+-=-移项整理得()()11212122kkkkafatFotBitBiFoFot+=++--,122aFoBi≤+同理,在室外的壁面上()()1112122kkkknnbfbntFotBitBiFoFot+++=++--,122bFoBi≤+(3)编写流程图本题中,室外温度的周期变化函数和墙体的初始需要假设设室外温度按余弦变化,墙体初始温度均匀为1℃。则()tf36cos2bktp=--Δ×℃1211nttt+====L℃计算机程序中使用的变量标识符:i节点的坐标变量t(i)节点i的温度k时间间隔变量计算机程序中输入数据:n沿x方向的网格划分数np控制打印间隔tm时间间隔总数tfa室内温度tfb室外温度No开始n=10,tm=480,tfa=20,Bia=6*0.24/(n-1)/0.43,Bib=15*0.24/(n-1)/0.43;Fo=0.43/1668/750*86400/tm/(0.24/(n-1))^2输入n,np,tm,Bia,Bib,Fok=1tfb=-3-6*cos((s-1)/tm*2*pi)k=k+1控制打印NP打印Bia,Bib,Fo打印“不稳定”停机NoYesYesktm打印tn,twNoYes(4)编写程序本题使用matlab软件,所编写的程序如下:clear;clc;n=10%x节点数tm=480%时间节点tfa=20;%室内温度np=10;Bia=6*0.24/(n-1)/0.43;%计算Bi及Fo准则数初始值Bib=15*0.24/(n-1)/0.43;Fo=0.43/1668/750*86400/tm/(0.24/(n-1))^2;p=(2*Bib+2)*Fo%显式差分格式的稳定性条件P1t=ones(1,n);%输入墙体内部温度的初始值,各节点均设为1tw=ones(1,tm+1);tn=ones(1,tm+1);if(p1)%迭代程序开始disp('不稳定')elsefors=1:tm+1tfb=-3-6*cos((s-1)/tm*2*pi);%室外温度fork=2:n-1t(1)=2*Fo*(t(2)+tfa*Bia)+(1-2*Fo-2*Bia*Fo)*t(1);t(n)=2*Fo*(t(n-1)+tfb*Bib)+(1-2*Fo-2*Bib*Fo)*t(n-1);t(k)=Fo*(t(k-1)+t(k+1))+(1-2*Fo)*t(k);endtn(s)=t(1);tw(s)=t(n);end%迭代程序结束Bia%输出计算结果BibFofork=1:tm/nptnn(1,k)=tn(1,np*k);tww(1,k)=tw(1,np*k);endtnnt=0:tm;%绘制温度变化图plot(24*w/tm,tn(w+1),'-r',24*w/tm,tw(w+1)),gridaxis([0,24,-10,22])set(gca,'xtick',[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24])set(gca,'ytick',[-10,-8,-6,-4,-2,0,2,4,6,8,10,12,14,16,18,20,22,24]);xlabel('时间'),ylabel('温度(℃)')legend('内壁面温度','外壁面温度',0)end(5)计算结果准则数初始值:Bia=0.3721,Bib=0.9302,Fo=0.0870时间内壁面温度(℃)外壁面温度(℃)时间内壁面温度(℃)外壁面温度(℃)0:308.22-1.3912:3014.553.891:009.45-1.8613:0014.654.061:3010.28-2.1413:3014.744.192:0010.89-2.3014:0014.834.272:3011.37-2.3614:3014.934.313:0011.76-2.3215:0015.014.303:3012.09-2.2015:3015.104.264:0012.36-2.0116:0015.184.184:3012.59-1.7616:3015.264.065:0012.79-1.4717:0015.343.915:3012.97-1.1317:3015.413.736:0013.13-0.7618:0015.483.536:3013.27-0.3718:3015.543.317:0013.400.0419:0015.603.077:3013.520.4619:3015.652.828:0013.640.8920:0015.702.578:3013.751.3120:3015.742.339:0013.861.7221:0015.772.089:3013.962.1121:3015.801.8510:0014.062.4922:0015.831.6410:3014.162.8322:3015.851.4511:0014.263.1523:0015.861.2811:3014.363.4423:3015.871.1412:0014.463.680:0015.881.03(6)温度变化图clear;clc;t=ones(1,11);q=ones(1,11);t(1)=250;t(11)=46;e=0.05;k=1;while1tt=t;fori=2:10a=43+0.08*(t(i-1)+t(i))/2;b=43+0.08*(t(i)+t(i+1))/2;t(i)=(a*t(i-1)+b*t(i+1))/(a+b);endk=k+1;ttt=abs(t(5)-tt(5));if(ttte)break;endendfori=2:11a=43+0.08*(t(i-1)+t(i))/2;q(i)=q(i-1)+a*(t(i)-t(i-1));endktq=q(i)/0.25终止误差0.05迭代次数59热流密度q=4.4745e+004温度分布节点i温度t(℃)12502231.88193213.35724194.40465174.99496155.08987134.64058113.5867991.85581069.36121146clear;clc;m=10tm=480tfa=20;Bia=6*0.24/(m-1)/0.43;Bib=15*0.24/(m-1)/0.43;Fo=0.43/1668/750*86400/tm/(0.24/(m-1))^2;p=(2*Bib+2)*FoBia=6*0.24/(m-1)/0.43;Bib=15*0.24/(m-1)/0.43;Fo=0.43/1668/750*86400/tm/(0.24/(m-1))^2;p=(2*Bib+