fin07

 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,T]=fin07(xl,xr,yb,yt,z,M,N,P,K,H,offl,offr)
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) = -3/(2*k)+(H/K);  % u_i,1
0032   A(i+(j-1)*m,i+1*m) = 2/k;       % u_i,2
0033   A(i+(j-1)*m,i+2*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-2)*m) = 2/k;       % u_i,n-1
0038   A(i+(j-1)*m,i+(j-3)*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 < ceil((n+1)/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 % create a notch in the cooling fin
0058 % by setting a grid of zeros in A
0059 % offl = 0.75;
0060 % offr = 1.0;
0061 for i = floor(offl*m):floor(offr*m),
0062     for j = floor(0.75*n):n,
0063         A(i+(j-1)*m,i-1+(j-1)*m)=0;      % u_i-1,j
0064         A(i+(j-1)*m,i+1+(j-1)*m)=0;      % u_i+1,j
0065         A(i+(j-1)*m,i+(j-1)*m)=0;  % u_i,j
0066         A(i+(j-1)*m,i+(j-2)*m)=0;        % u_i,j-1
0067         A(i+(j-1)*m,i+j*m)=0;            % u_i,j+1
0068         b(i+(j-1)*m)=0;
0069     end
0070 end
0071 % now fix the boundaries conditions for our
0072 % modified fin
0073 % RIGHT SIDE
0074 for j = floor(0.75*n):n,
0075     i = ceil(offl*m);
0076     A(i+(j-1)*m,i+(j-1)*m) = -3/(2*h)+(H/K);  % u_m,j
0077     A(i+(j-1)*m,i-1+(j-1)*m) = 2/h;       % u_m-1,j
0078     A(i+(j-1)*m,i-2+(j-1)*m) = -1/(2*h);    % u_m-2,j
0079 end
0080 % TOP SIDE
0081 for i = floor(offl*m):floor(offr*m),
0082 %     j = floor(0.75*n);
0083     j = n;
0084     A(i+(j-1)*m,i+(j-1)*m) = -3/(2*k)+(H/K);  % u_i,n
0085     A(i+(j-1)*m,i+(j-2)*m) = 2/k;       % u_i,n-1
0086     A(i+(j-1)*m,i+(j-3)*m) = -1/(2*k);    % u_i,n-2
0087 end
0088 v=A\b; % solve for solution in v labeling
0089 w=reshape(v(1:mn),m,n); %translate from v to w
0090 w=w+20;
0091 T=max(max(w));
0092 figure(2);
0093 mesh(x,y,w');
0094 view(180,0);
0095 title(['Cooling Fin - ' num2str(P,4) ' Watts (RIGHT)']);
0096 zlabel('Temperature (C)');
0097 xlabel('Fin Length (cm)');
0098 ylabel('Fin Length (cm)');
0099 figure(4);
0100 mesh(x,y,w');
0101 view(0,90);
0102 title(['Cooling Fin - ' num2str(P,4) ' Watts (TOP)']);
0103 xlabel('Fin Length (cm)');
0104 ylabel('Fin Length (cm)');