%                                                        
%   CHIP2D_SOV1     2-D Heat Conduction in a Rectangular Silicon Chip
%                                                        
%   Analytical Soln using Separation of Variables to the following heat 
%   conduction problem:
%         Txx(x,y) + Tyy(x,y)  = -Q/k  (Q = internal volumetric heat generation rate)
%   with   left  BC:    Tx(0,y) = 0   (symmetry)
%         bottom BC:    Ty(x,0) = 0   (insulated)
%          right BC: -k*Tx(Lc/2,y) = hs*(T(Lc/2,y) - Tinf)  (convection)           
%           top  BC: -k*Ty(x,Hc) = ht*(T(x,Hc) - Tinf)  (convection) 
%   where the domain of interest is defined as 0 < x < Lc/2  and  0 < y < Hc and 
%   hs and ht are the side and top heat transfer coefficients, respectively.          
%   
%   The goal here is to compute and plot the analytical solution and to determine
%   the peak chip temperature.  
%
%   File prepared by J. R. White, UMass-Lowell (Nov. 2003)
%

%
%   getting started
      clear all,   close all,   nfig = 0;
      global hs kc Lc Q EIGN EIGM
%                                                                     
%   problem data                                                     
      Tinf = 20;         % air temperature on top and sides of chip (C)
      hs = 50;           % heat transfer coeff on chip sides (W/m^2-C)
      ht = 450;          % heat transfer coeff on top of chip (W/m^2-C)
      kc = 25;           % thermal conductivity of chip (W/m-C)
      Q = 1e+7;          % volumetric heat generation rate (W/m^3)
      Lc = 0.009;        % width of chip in x direction (m)
      Hc = 0.003;        % height of chip in y direction (m)
      maxt = 15;         % max number of nonzero terms (converges quickly)
      tol = 0.001;       % convergence limit for series expansion
%
%   points to evaluate function (domain of interest is a rectangle)
      Nx = 37;    x = linspace(0,Lc/2,Nx);  
      Ny = 31;    y = linspace(0,Hc,Ny);
      [xx,yy] = meshgrid(x,y);   
%
%   evaluate W(x)  (part of total solution)
      W = Q*Lc/(2*hs) + (Q/2/kc)*((Lc/2)^2 - x.^2);                                    
%   
%   calc eigenvalues (where guesses come from knowledge of the tan(p) function)
      ln = zeros(maxt,1);
      for n = 1:maxt
        g1 = (2*n-2)*pi/2;  g2 = (2*n-1)*pi/2;
        pn = fzero('chip2d_sov1a',[g1 g2]);
        ln(n) = 2*pn/Lc;
      end
%
%   evaluate integral In
      In = zeros(maxt,1);
      for n = 1:maxt
        EIGN = ln(n);
        In(n) = quadl('chip2d_sov1b',0,Lc/2);
      end
%
%   evaluate the norm of the eigenfunctions (also check orthogonality)
      norme = zeros(maxt,maxt);
      for m = 1:maxt
        EIGM = ln(m);
        for n = 1:maxt
          EIGN = ln(n);
          norme(m,n) = quadl('chip2d_sov1c',0,Lc/2);
        end
      end
%
%   let's do a quick SPY plot to check orthogonality
      nfig = nfig+1;  figure(nfig)
      maxnorm = max(max(norme));
      spy(abs(norme/maxnorm) > 0.001), axis tight
      title('CHIP2D\_SOV1:  SPY Plot of NormE(m,n) as an Orthogonality Check')
      xlabel('n value');   ylabel('m value')
%
%   calc expansion coeffs (use 'dot' arithmetic)
      bn = -(ht./diag(norme)).*In./(ht*cosh(ln*Hc) + kc*ln.*sinh(ln*Hc));
%
%   determine v(x,y) spatial temperature distribution              
      n = 0;   mrerr = 1.0;   v = zeros(Ny,Nx); 
      while mrerr > tol   &   n < maxt
        n = n+1;   vn = bn(n)*cos(ln(n)*xx).*cosh(ln(n)*yy);
        v = v + vn;  mrerr = max(max(abs(vn)));
      end
      fprintf(1,'\n  Needed %3i terms for convergence of v(x,y) \n\n',n)
%
%   now form desired temperature distribution  -->  T(x,y) = v(x,y) + W(x) + Tinf
      T = zeros(Ny,Nx);
      for j = 1:Ny
        T(j,:) = v(j,:) + W + Tinf;
      end
%
%   plot T(x,y) for a few y values              
      nfig = nfig+1;  figure(nfig)
      plot(x*100,T(Ny:-10:1,:),'LineWidth',2),grid
      title('CHIP2D\_SOV1:  T(x) profiles for various y values')
      xlabel('x location from chip centerline (cm)');    
      ylabel('Temperature (C)')                                 
      for j = Ny:-10:1,  gtext(['y = ',num2str(y(j)*100),' cm']);   end   
%
%  use a surface plot (top view) with contours
      Tmax = max(max(T));
      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, [82 82.5 83 83.5 84]);
      clabel(cc), set(hh,'EdgeColor','k')
      axis equal
      rr = axis; rr(1:4) = [0 Lc*100/2 0 Hc*100];  axis(rr)
      title(['CHIP2D\_SOV1:  T(x,y) in 2-D electronic chip (T_{max} = ', ...
            num2str(Tmax,'%5.2f'),' ^oC)'])
      xlabel('x location from chip centerline (cm)');   ylabel('y location (cm)')
      zlabel('Temperature (C)'),   hold off                                 
%
%   end of problem
