Matlab编写的流体计算和传热程序.zip

% Twodimensional heat conduction % Finite Volume Method % Line by Line TDMA with SOR clear all; x=[];y=[];T=[];Told=[];Su=[];Sp=[];ap=[];ae=[];aw=[];as=[];an=[]; Q=[];P=[]; great = 1e20; lambda = 10; % thermal conductivity alfa = 10; % heat transfer coefficient dt = great; % Time step. If great steady state density = 6000;% density cp = 500;% heat capacity Lx = 0.12; % length x-direction Ly = 0.12; % length y -direction Tfluid = 20; % Fluid temperature Tinit = 50; % Initial guess and top- and bottom tempearature %cv_x = input('Number of CVs in x-direction = ') %cv_y = input('Number of CVs in y-direction = ') cv_x=10;cv_y=10; ni = cv_x+2; % grid points x-direction nj = cv_y+2; % grid points y-direction dx = Lx/cv_x; dy = Ly/cv_y; x(1) = 0; x(2)=dx/2; for i = 3:ni-1 x(i)=x(i-1)+dx; end; x(ni)=Lx; y(1) = 0; y(2)=dy/2; for j = 3:nj-1 y(j)=y(j-1)+dy; end y(nj)=Ly; % Initial values and coefficients for i = 1:ni for j = 1:nj T(i,j) = Tinit; %Initial temperature Told(i,j) = Tinit; T(i,1) = 50; T(i,nj) = 50 Su(i,j)=0; %Initial indendendent source term Sp(i,j)=0; %Initial dependent source term ae(i,j) = lambda*dy/dx; aw(i,j) = lambda*dy/dx; an(i,j) = lambda*dx/dy; as(i,j) = lambda*dx/dy; dV = dx*dy; ap0 = density*cp*dV/dt; if i==2 % convective heat transfer boundary Su(i,j) = Tfluid/(1/alfa+dx/(2*lambda))*dy/dV; Sp(i,j) = -1/(1/alfa+dx/(2*lambda))*dy/dV; aw(i,j) = 0; end; if i==ni-1 % insulated boundary ae(i,j) = 0; end if j==2 % bottom boundary, given temperature as(i,j)=2*lambda*dx/dy; end if j==nj-1 % top boundary, given temperature an(i,j)=2*lambda*dx/dy; end ap(i,j) = ae(i,j)+aw(i,j)+an(i,j)+as(i,j)-Sp(i,j)*dV+ap0; end; end; %%%%%%%%%%% maxres = 1.0e-6; maxit = 100; time=0; maxtime = 100; omega=1; while (time < (maxtime+dt/2)) Told=T; sumres = 1; counter = 0; while (sumres>maxres&counter<maxit) %Sweep in j-direction for i = 2:ni-1 P(1) = 0; Q(1) = T(1); for j = 2:nj-1 a=ap(i,j)/omega;b=an(i,j);c=as(i,j); d=ae(i,j)*T(i+1,j)+aw(i,j)*T(i-1,j)+Su(i,j)*dV +ap0*Told(i,j)... +(1-omega)*a*T(i,j); P(j) = b/(a-c*P(j-1)); Q(j) = (c*Q(j-1)+d)/(a-c*P(j-1)); end; for j = nj-1:-1:2 T(i,j) = P(j)*T(i,j+1)+Q(j); end; end; %Sweep in i-direction for j = 2:nj-1 P(1) = 0; Q(1) = T(1); for i = 2:ni-1 a=ap(i,j)/omega;b=ae(i,j);c=aw(i,j); d=an(i,j)*T(i,j+1)+as(i,j)*T(i,j-1)+Su(i,j)*dV+ap0*Told(i,j)... +(1-omega)*a*T(i,j); P(i) = b/(a-c*P(i-1)); Q(i) = (c*Q(i-1)+d)/(a-c*P(i-1)); end; for i = ni-1:-1:2 T(i,j) = P(i)*T(i+1,j)+Q(i); end; end; sumres = 0; % Calculate residual for i = 2:ni-1 for j = 2:nj-1 res = abs(ap(i,j)*T(i,j)-(ae(i,j)*T(i+1,j)+aw(i,j)*T(i-1,j)+... an(i,j)*T(i,j+1)+as(i,j)*T(i,j-1)+Su(i,j)*dV+ap0*Told(i,j))); sumres=sumres+res; end; end sumerr=sumres counter = counter + 1 end;%end while time = time+dt; end; %end time loop % Calculate boundary values for j = 2:nj-1 T(1,j)=(alfa*Tfluid+lambda/(dx/2)*T(2,j))/(alfa+lambda/(dx/2)); T(ni,j) = T(ni-1,j); end; % pcolor(x,y,T');shading interp;xlabel('x');ylabel('y');title('Temperature distribution');colorbar; %