tem.zip_matlab 温度_matlab 温度场_变化边界_温度变化_热TEM

function [temx]=tem(lmt,ro,cp,d,pretem,q,time,x) %一维导热问题温度场的显式有限差分求解 %考虑材料的变热物性 %输入参数： % lmt,ro,cp为热导率、密度、比热 % pretem为初始温度 % q热流密度 % time作用时间 % x测点位置 %返回参数可选： % temx测点处的温度历程 % temt在t=time时刻的温度场 % temxt测点位置在time时刻的温度值 timegrid=0.1; %sym timesteps; timesteps=min(max(501,fix(time/timegrid)+1),10001);%时间节点数 timegrid=time/(timesteps-1);%时间步长度 dgrid=0.01; dsteps=min(max(10,fix(d/dgrid)),50);%dsteps取值不能过大，否则计算不收敛 dgrid=d/dsteps; T=zeros(timesteps,dsteps+4); T(1,:)=pretem; for i=1:timesteps+1 T(i,1)=i-1; end for i=1:timesteps qq=get(q,(i-1)*timegrid,1); lmti=get(lmt,T(i,3),1); %T(i,2)=2*dgrid*qq/lmt+T(i,4); %不考虑表面辐射 e=0.8; zerotem=-273.15; xgm=5.67e-8; qe=e*xgm*((T(i,3)-zerotem)^4-(pretem-zerotem)^4); T(i,2)=2*dgrid*(qq-qe)/lmti+T(i,4); %考虑表面辐射 for j=3:dsteps+3 lmti=get(lmt,T(i,j),1); roi=get(ro,T(i,j),1); cpi=get(cp,T(i,j),1); alf=lmti/roi/cp; r=alf*timegrid/dgrid^2; T(i+1,j)=(1-2*r)*T(i,j)+r*(T(i,j-1)+T(i,j+1)); end %qq=get(q,i*timegrid,1); T(i+1,dsteps+4)=T(i+1,dsteps+2);%背壁绝热条件 lmti=get(lmt,T(i,3),1); %T(i+1,dsteps+4)=T(i+1,dsteps+2)-2*dgrid*qh/lmti;%背壁换热条件 end t=time; %%%温度场和温度历程输出%%% dgrids=x/dgrid+3; intd=fix(dgrids); timesteps=round(t/timegrid+1); Temxt=T(timesteps,intd)*(intd+1-dgrids)+T(timesteps,intd+1)*(dgrids-intd); %%%Txt2=T(intt+1,intd)*(intd+1-dgrids)+T(intt+1,intd+1)*(dgrids-intd); %%%Temxt=Txt1*(intt+1-tgrids)+Txt2*(tgrids-intt);%返回x位置在t时刻的温度 %tem=T(timesteps+1,intd)*(intd-dgrids+1)+T(timesteps+1,intd+1)*(dgrids-intd); time=timegrid*T(:,1); %TT=T(1:timesteps,3:dsteps+3);%各节点温度历程 temx=zeros(timesteps,2); temx(:,1)=time(1:timesteps); if x~=d%返回x位置的温度历程 temx(:,2)=T(1:timesteps,intd)*(intd-dgrids+1)+T(1:timesteps,intd+1)*(dgrids-intd); else temx(:,2)=T(1:timesteps,dgrids); end plot(temx(:,1),temx(:,2)); depth=zeros(dsteps+1); for i=1:dsteps+1 depth(i)=dgrid*(i-1); end Temt=T(timesteps,3:dsteps+3);%t时刻温度分布 plot(depth,Temt); end function [y]=get(data,x,f) [m,n]=size(data); if m==1 y=data; elseif f==3 y=interp1(data(:,1),data(:,2),x,'cublic'); elseif f==4 y=interp1(data(:,1),data(:,2),x,'spline'); else y=interp1(data(:,1),data(:,2),x); end end