%                                                        
%   BLOCK2D_SOV1     2-D Heat Conduction in a Rectangular Block
%                                                        
%   Analytical Soln using Separation of Variables (SOV) to the following heat 
%   conduction problem:
%     uxx(x,y) + uyy(x,y)  = 0    with ux(0,y) = 0   u(Lx/2,y) = 0            
%                                       u(x,0) = 0  -kuy(x,Ly) = h(u - Te)
%     where  u(x,y) = T(x,y) - Tw   and   Te = Tinf - Tw  
%   
%   The goal here is to plot the analytical solution a variety of ways to
%   help visualize the 2-D temperature distribution.  We also want to compute
%   the heat loss or gain across each face of the rectangular block.  Only the 
%   right half of the block is modeled  --  a symmetry condition at x = 0 
%   handles this (note that the heat flows given are only for the portion of 
%   the block modeled).
%
%   File prepared by J. R. White, UMass-Lowell (Aug. 2003)
%

%
%   getting started
      clear all,   close all,   nfig = 0;
%                                                                     
%   problem data                                                     
      Tw = 50;         % wall temperature on left, bottom, and right sides (C)
      Tinf = 100;      % air temperature on top of block (C)
      Te = Tinf-Tw;    % difference between air and wall temps (C)
      h = 100;         % heat transfer coeff (W/m^2-C)
      k = 1;           % thermal conductivity of block (W/m-C)
      Lx = 0.06;       % width of block in x direction (m)
      Ly = 0.09;       % height of block in y direction (m)
      maxt = 50;       % max number of nonzero terms
      tol = 0.01;      % convergence limit for series expansion (absolute temp diff)
%
%   points to evaluate function (domain of interest is a rectangle)
      Nx = 31;    x = linspace(0,Lx/2,Nx);  
      Ny = 37;    y = linspace(0,Ly,Ny);
      [xx,yy] = meshgrid(x,y);                                       
%
%   calc eigenvalues and expansion coeffs
%   Note: Since the bn values become so small so quickly, we can expect that 
%         numerical errors may become a problem (small # * large # in series expansion).  
%         Because of this I really should implement the series expansion differently
%         than given here --- use recursive form.  However, our focus here is
%         simply to get a reasonable visualization of the temperature and heat flux
%         distributions, so we will keep the implementation on the simple side (and
%         just be careful to make sure that the numerical errors do not affect the 
%         qualitative observations)
      dd = pi/Lx;     cc = 4*h*Te/Lx;                     
      for n = 1:maxt                            
        ln(n) = (2*n-1)*dd;
        bn(n) = (-1)^(n+1)*cc/ln(n)/(h*sinh(ln(n)*Ly) + k*ln(n)*cosh(ln(n)*Ly));
      end                                        
%
%   determine spatial temperature distribution              
      n = 0;   merr = 1.0;   u = zeros(Ny,Nx); 
      while merr > tol   &   n < maxt
        n = n+1;   un = bn(n)*cos(ln(n)*xx).*sinh(ln(n)*yy);
        u = u + un;    merr = max(max(abs(un)));  
      end
      fprintf(1,'\n  Needed %3i terms for convergence of u(x,y) \n\n',n)
%
%   now form actual temperature distribution of interest  --  T(x,y)
      T = u + Tw;
%
%   plot T(x,y) for a few y values              
      nfig = nfig+1;  figure(nfig)
      plot(x*100,T(Ny:-8:13,:),'LineWidth',2),grid
      title('Block2D\_SOV1:  T(x) profiles for various y values')
      xlabel('x location from block centerline (cm)');    
      ylabel('Temperature (C)')                                 
      for j = Ny:-8:13,  gtext(['y = ',num2str(y(j)*100),' cm']);   end   
%
%   various 3-D plots                           
      nfig = nfig+1;  figure(nfig)
      surf(xx*100,yy*100,T), shading interp, colorbar
      axis([0 Lx*100/2 0 Ly*100 50 90])
      title('Block2D\_SOV1:  T(x,y) in 2-D block')
      xlabel('x location (cm)');   ylabel('y location (cm)')
      zlabel('Temperature (C)')                                 
%
      nfig = nfig+1;  figure(nfig)
      surf(xx*100,yy*100,T), shading interp, colorbar, hold on
      [cc,hh] = contour3(x*100,y*100,T, [80 70 60 55 52.5]);
      clabel(cc), set(hh,'EdgeColor','k')
      axis([0 Lx*100/2 0 Ly*100 50 90])
      title('Block2D\_SOV1:  T(x,y) in 2-D block')
      xlabel('x location (cm)');   ylabel('y location (cm)')
      zlabel('Temperature (C)'),   hold off                                 
%
      nfig = nfig+1;  figure(nfig)
      surf(xx*100,yy*100,T), shading interp, view(2), colorbar, hold on
      [cc,hh] = contour3(x*100,y*100,T, [80 70 60 55 52.5]);
      clabel(cc), set(hh,'EdgeColor','k')
      axis([0 Lx*100/2 0 Ly*100 50 90])
      title('Block2D\_SOV1:  T(x,y) in 2-D block')
      xlabel('x location (cm)');   ylabel('y location (cm)')
      zlabel('Temperature (C)'),   hold off                                 
%
      nfig = nfig+1;  figure(nfig)
      surf(xx*100,yy*100,T), shading interp, view(2), colorbar, hold on
      [cc,hh] = contour3(x*100,y*100,T, [80 70 60 55 52.5]);
      clabel(cc), set(hh,'EdgeColor','k')
      axis equal
      axis([0 Lx*100/2 0 Ly*100 50 90])
      title('Block2D\_SOV1:  T(x,y) in 2-D block')
      xlabel('x location (cm)');   ylabel('y location (cm)')
      zlabel('Temperature (C)'),   hold off                                 
%
%   determine spatial heat flux distributions 
%   x-directed heat flux             
      n = 0;   merr = 1.0;   qx = zeros(Ny,Nx); 
      while merr > tol   &   n < maxt
        n = n+1;   qxn = k*bn(n)*ln(n)*sin(ln(n)*xx).*sinh(ln(n)*yy);
        qx = qx + qxn;   merr = max(max(abs(qxn)));
      end
      fprintf(1,'\n  Needed %3i terms for convergence of qx(x,y) \n\n',n)
%   y-directed heat flux             
      n = 0;   merr = 1.0;   qy = zeros(Ny,Nx); 
      while merr > tol   &   n < maxt
        n = n+1;   qyn = -k*bn(n)*ln(n)*cos(ln(n)*xx).*cosh(ln(n)*yy);
        qy = qy + qyn;   merr = max(max(abs(qyn)));
      end
      fprintf(1,'\n  Needed %3i terms for convergence of qy(x,y) \n\n',n)
%
%   plot the heat flux as a vector plot
      nfig = nfig+1;  figure(nfig)
      quiver(xx*100,yy*100,qx,qy),grid
      axis([0 Lx*100/2 0 Ly*100])
      title('Block2D\_SOV1:  Heat flux vector plot for 2-D block')
      xlabel('x location (cm)');    
      ylabel('y location (cm)')                                 
%
%   plot magnitude of the heat flux
      qmag = sqrt(qx.*qx + qy.*qy);
      nfig = nfig+1;  figure(nfig)
      surf(xx*100,yy*100,qmag/1000), shading interp, colorbar, hold on
      [cc,hh] = contour3(x*100,y*100,qmag/1000);
      clabel(cc), set(hh,'EdgeColor','k')
      axis([0 Lx*100/2 0 Ly*100 0 14])
      title('Block2D\_SOV1:  Distribution of heat flux magnitude (kW/m^2)')
      xlabel('x location (cm)');   ylabel('y location (cm)')
      zlabel('Temperature (C)'),   hold off                                 
%
%   now compute the heat flows across each surface (use all terms)
      qtop1 = -k*sum(bn.*sin(ln*Lx/2).*cosh(ln*Ly));
      qtop2 = h*sum(bn.*sin(ln*Lx/2).*sinh(ln*Ly)./ln) - h*Te*Lx/2;
      qrt = k*sum(bn.*sin(ln*Lx/2).*(cosh(ln*Ly)-1));
      qbot = -k*sum(bn.*sin(ln*Lx/2));
      qlt = 0;
%
      fprintf(1,'\n  Block2D_SOV1 -- Heat transfer (W) across boundaries for the right half of a 2-D block \n')
      fprintf(1,'    in top boundary (conduction):   %8.2f  \n',-qtop1)
      fprintf(1,'    in top boundary (convection):   %8.2f  \n',-qtop2)
      fprintf(1,'    out right boundary:             %8.2f  \n',qrt)
      fprintf(1,'    out bottom boundary:            %8.2f  \n',-qbot)
      fprintf(1,'    out left boundary:              %8.2f  \n',-qlt)
%
%   end of problem
