fin02

 Program 8.5 Finite difference solver for 2D Poisson equation 
   with Dirichlet boundary conditions on a rectangle 
 Input: rectangle domain [xl,xr]x[yb,yt] with MxN space steps 
 Output: matrix w holding solution values 
 Example usage: w=poisson(0,1,1,2,4,4)

0001 % Program 8.5 Finite difference solver for 2D Poisson equation
0002 %   with Dirichlet boundary conditions on a rectangle
0003 % Input: rectangle domain [xl,xr]x[yb,yt] with MxN space steps
0004 % Output: matrix w holding solution values
0005 % Example usage: w=poisson(0,1,1,2,4,4)
0006 function w=fin02(xl,xr,yb,yt,z,M,N) 
0007 m=M+1; n=N+1; mn=m*n; 
0008 h=(xr-xl)/M; h2=h^2;
0009 k=(yt-yb)/N; k2=k^2; 
0010 x=xl+(0:M)*h; % set mesh values
0011 y=yb+(0:N)*k; 
0012 A=zeros(mn,mn);
0013 b=zeros(mn,1);
0014 L = 2; % length of input [ cm ]
0015 P = 5; % total power [ W ]
0016 H = 0.005; % convective heat transfer coefficient [ W / (cm^2 * degC) ]
0017 K = 1.68; % thermal conductivity [ W / (cm^2 * degC) ]
0018 delta = z; % depth in z-direction [ cm ]
0019 for i=2:m-1, % interior points
0020   for j=2:n-1, 
0021     A(i+(j-1)*m,i-1+(j-1)*m)=1/h2;      % u_i-1,j
0022     A(i+(j-1)*m,i+1+(j-1)*m)=1/h2;      % u_i+1,j
0023     A(i+(j-1)*m,i+(j-1)*m)=-2/h2-2/k2-((2*H)/(K*delta));  % u_i,j
0024     A(i+(j-1)*m,i+(j-2)*m)=1/k2;        % u_i,j-1
0025     A(i+(j-1)*m,i+j*m)=1/k2;            % u_i,j+1
0026     b(i+(j-1)*m)= 0;
0027   end 
0028 end 
0029 for i=1:m % bottom and top boundary points
0030   j=1;
0031   A(i+(j-1)*m,i+(j-1)*m) = -3/(2*k)+(H/K);  % u_i,1
0032   A(i+(j-1)*m,i+(j+1-1)*m) = 2/k;           % u_i,2
0033   A(i+(j-1)*m,i+(j+2-1)*m) = -1/(2*k);      % u_i,3
0034   b(i+(j-1)*m)=0;
0035   j=n;
0036   A(i+(j-1)*m,i+(j-1)*m) = -3/(2*k)+(H/K);  % u_i,n
0037   A(i+(j-1)*m,i+(j-1-1)*m) = 2/k;           % u_i,n-1
0038   A(i+(j-1)*m,i+(j-2-1)*m) = -1/(2*k);      % u_i,n-2
0039   b(i+(j-1)*m)=0;
0040 end 
0041 for j=2:n-1 % left and right boundary points
0042   i=1;
0043   A(i+(j-1)*m,i+(j-1)*m) = -3/(2*h)+(H/K);  % u_1,j
0044   A(i+(j-1)*m,i+1+(j-1)*m) = 2/h;           % u_2,j
0045   A(i+(j-1)*m,i+2+(j-1)*m) = -1/(2*h);        % u_3,j
0046   if j < n/2,
0047       b(i+(j-1)*m)= -P/(L*K*delta); % power input on left side
0048   else
0049       b(i+(j-1)*m)=0;
0050   end
0051   i=m;
0052   A(i+(j-1)*m,i+(j-1)*m) = -3/(2*h)+(H/K);  % u_m,j
0053   A(i+(j-1)*m,i-1+(j-1)*m) = 2/h;           % u_m-1,j
0054   A(i+(j-1)*m,i-2+(j-1)*m) = -1/(2*h);        % u_m-2,j
0055   b(i+(j-1)*m)=0;
0056 end 
0057 v=A\b; % solve for solution in v labeling
0058 w=reshape(v(1:mn),m,n); %translate from v to w
0059 w=w+20;
0060 T=max(max(w));
0061 figure(2);
0062 mesh(x,y,w');
0063 title(['Cooling Fin - Max Temp=' num2str(T,4) ' Power (W)=' num2str(P)]);
0064 zlabel('Temperature (C)');
0065 ylabel('Fin Length (cm)');
0066 xlabel('Fin Length (cm)');