%   Fluid flow through a heated pipe
%	
%   The PDE that models the heated pipe is
%
%      T = -v*T  + (4*hc)/(D*rho*Cp)*(Tw - T)               (1)
%       t      z
%
%   with the initial and boundary conditions
%
%      T(z,0) = T0, T(0,t) = Te                          (2)(3)
%
%   where
%
%      T     fluid temperature
%
%      t     time
%
%      z     position along the pipe
%
%      v     fluid velocity
%
%      hc    fluid to pipe wall heat transfer coefficient
%
%      D     pipe diameter
%
%      rho   fluid density
%
%      Cp    fluid specific heat
%
%      T0    initial fluid temperature
%
%      Te    entering fluid temperature
%
%      Tw    pipe wall temperature
%
%      L     pipe length
%
%   Of particular interest is the fluid exiting temperature, T(L,t).
%
%   The method of lines (MOL) solution of eq. (1) is coded below.
%   Specifically, the spatial derivative in eq. (1) is replaced
%   by one of four approximations. The resulting system of ODEs
%   in t, defined on a 21-point grid in z, is then integrated by
%   the Euler method.
%
%   The analytical solution to eqs. (1) to (3) is also programmed
%   for comparison with the numerical solution.
%
%   The following code is run for the four special cases:
%
%      Case 1: hc = 2, T0 = 25, Te = 25, Tw = 100
%
%      Case 2: hc = 0.2, T0 = 25, Te = 25, Tw = 100
%
%      Case 3: hc = 0, T0 = 25, Te = 25, Tw (NA)
%
%      Case 4: hc = 0, T0 = 25, Te = 100, TW (NA)
%
% Open output files
  fid1=fopen('htex1.out','w');
  fid2=fopen('htex1.plt','w');
%
% Model parameters (defined in the comments above)
  cp=1.0;
  rho=1.0;
  v=50.0;
  d=1.0;
  zl=100.0;
  n=21;
  dz=zl/(n-1);
%
% Select finite difference (FD) approximation of the spatial
% derivative
%
%   ifd = 1: Centered approximations
%
%   ifd = 2: Two point upwind approximation
%
%   ifd = 3: Five point, biased upwind approximation
%
%   ifd = 4: van Leer flux limiter
%
  ifd=1;
%
% Cases listed above
  for ncase=1:4
%
% Initial, final times, integration interval, number of Euler
% steps for each output
  t=0.0; tf=5.0; h=0.001; nout=250;
  if ncase==1
     hc=2.0;
     T0=25.0;
     Te=25.0;
     Tw=100.0;
  end
  if ncase==2
     hc=0.2;
     T0=25.0;
     Te=25.0;
     Tw=100.0;
  end
  if ncase==3
     hc=0.0;
     T0=25.0;
     Te=25.0;
     Tw=100.0;
  end
  if ncase==4
     hc=0.0;
     T0=25.0;
     Te=100.0;
     Tw=100.0;
  end
%
% Initial conditions
  for i=1:n
    T(i)=T0;
  end
%
% Write h
  fprintf(fid1,'\n\n ncase = %5d    h = %10.3e\n\n',ncase,h);
%
% Write heading
  fprintf(fid1,'   t     t-zl/v   T(zl,t)  Ta(zl,t)      diff\n');
%
% Integrate until t = tf
  while t < tf*1.0001
%
% Monitor solution by displaying t
  t
%
% Write selected output, including analytical solution and difference
% between numerical and analytical solutions
  alpha=4.0*hc/(d*rho*cp);
  Ta=Tw+(T0-Tw)*exp(-alpha*t);
  if(t > zl/v)Ta=Ta+exp(-alpha/v*zl)*((Te-Tw)-(T0-Tw)*...
                      exp(-alpha*(t-zl/v))); end
  diff=T(n)-Ta;
  fprintf(fid1,'%5.2f%10.2f%10.2f%10.2f%10.3f\n',...
          t,t-zl/v,T(n),Ta,diff);
  fprintf(fid2,'%5.2f%10.2f%10.2f%10.2f%10.3f\n',...
          t,t-zl/v,T(n),Ta,diff);
%
% Take nout Euler steps
  for iout=1:nout
%
%   Temporal derivatives
%
%     Centered approximations
      if(ifd==1)
%
%       Boundary condition at z = 0
        T(1)=Te;
%
%       Spatial derivative
        [Tz]=dss002(0.0,zl,n,T);
%
%       Temporal derivative
        Tt(1)=0.0;
        for i=2:n
          Tt(i)=-v*Tz(i)+alpha*(Tw-T(i));
        end
%
%     End of three point centered approximation
      end
%
%     Two point upwind approximation
      if(ifd==2)
%
%       Boundary condition at z = 0
        T(1)=Te;
%
%       Spatial derivative
        [Tz]=dss012(0.0,zl,n,T,v);
%
%       Temporal derivative
        Tt(1)=0.0;
        for i=2:n
          Tt(i)=-v*Tz(i)+alpha*(Tw-T(i));
        end
%
%     End of two point upwind approximation
      end
%
%     Five point, biased upwind approximation
      if(ifd==3)
%
%       Boundary condition at z = 0
        T(1)=Te;
%
%       Spatial derivative
        [Tz]=dss020(0.0,zl,n,T,v);
%
%       Temporal derivative
        Tt(1)=0.0;
        for i=2:n
          Tt(i)=-v*Tz(i)+alpha*(Tw-T(i));
        end
%
%     End of five point, biased upwind approximation
      end
%
%     van Leer flux limiTer
      if(ifd==4)
%
%       Boundary condition at z = 0
        T(1)=Te;
%
%       Spatial derivative
        [Tz]=vanl2(0.0,zl,n,T,v);
%
%       Temporal derivative
        Tt(1)=0.0;
        for i=2:n
          Tt(i)=-v*Tz(i)+alpha*(Tw-T(i));
        end
%
%     End of van Leer flux limiter
      end
%
%   Take Euler step
    for i=1:n
      T(i)=T(i)+Tt(i)*h;
    end
    t=t+h;
%
% Next Euler step
  end
%
% Next output
  end
%
% Next case
  end