fin04

 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 z thickness, MxN space steps
 and P power input [Watts]
 Output: matrix w holding solution values 
 Example usage: w=fin04(0,4,0,4,.1,10,10,5)

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 z thickness, MxN space steps
0004 % and P power input [Watts]
0005 % Output: matrix w holding solution values
0006 % Example usage: w=fin04(0,4,0,4,.1,10,10,5)
0007 function [w,T]=fin04(xl,xr,yb,yt,z,M,N,P) 
0008 m=M+1; n=N+1; mn=m*n; 
0009 h=(xr-xl)/M; h2=h^2;
0010 k=(yt-yb)/N; k2=k^2; 
0011 x=xl+(0:M)*h; % set mesh values
0012 y=yb+(0:N)*k; 
0013 A=zeros(mn,mn);
0014 b=zeros(mn,1);
0015 L = 2; % length of input [ cm ]
0016 %P = 5; % total power [ W ]
0017 H = 0.005; % convective heat transfer coefficient [ W / (cm^2 * degC) ]
0018 K = 3.85; % thermal conductivity of COPPER[ W / (cm^2 * degC) ]
0019 delta = z; % depth in z-direction [ cm ]
0020 for i=2:m-1, % interior points
0021   for j=2:n-1, 
0022     A(i+(j-1)*m,i-1+(j-1)*m)=1/h2;      % u_i-1,j
0023     A(i+(j-1)*m,i+1+(j-1)*m)=1/h2;      % u_i+1,j
0024     A(i+(j-1)*m,i+(j-1)*m)=-2/h2-2/k2-((2*H)/(K*delta));  % u_i,j
0025     A(i+(j-1)*m,i+(j-2)*m)=1/k2;        % u_i,j-1
0026     A(i+(j-1)*m,i+j*m)=1/k2;            % u_i,j+1
0027     b(i+(j-1)*m)=0;
0028   end 
0029 end 
0030 for i=1:m % bottom and top boundary points
0031   j=1;
0032   A(i+(j-1)*m,i) = -3/(2*k)+(H/K);  % u_i,1
0033   A(i+(j-1)*m,i+j*m) = 2/k;       % u_i,2
0034   A(i+(j-1)*m,i+(j+1)*m) = -1/(2*k);    % u_i,3
0035   b(i+(j-1)*m)=0;
0036   j=n;
0037   A(i+(j-1)*m,i+(j-1)*m) = -3/(2*k)+(H/K);  % u_i,n
0038   A(i+(j-1)*m,i+(j-2)*m) = 2/k;       % u_i,n-1
0039   A(i+(j-1)*m,i+(j-3)*m) = -1/(2*k);    % u_i,n-2
0040   b(i+(j-1)*m)=0;
0041 end 
0042 for j=2:n-1 % left and right boundary points
0043   i=1;
0044   A(i+(j-1)*m,i+(j-1)*m) = -3/(2*h)+(H/K);  % u_1,j
0045   A(i+(j-1)*m,i+1+(j-1)*m) = 2/h;       % u_2,j
0046   A(i+(j-1)*m,i+2+(j-1)*m) = -1/(2*h);    % u_3,j
0047   
0048   if j < n/2, 
0049       b(i+(j-1)*m)= -P/(L*K*delta); % power input on left side
0050   else
0051       b(i+(j-1)*m)=0;
0052   end
0053   
0054   i=m;
0055   A(i+(j-1)*m,i+(j-1)*m) = -3/(2*h)+(H/K);  % u_m,j
0056   A(i+(j-1)*m,i-1+(j-1)*m) = 2/h;       % u_m-1,j
0057   A(i+(j-1)*m,i-2+(j-1)*m) = -1/(2*h);    % u_m-2,j
0058   b(i+(j-1)*m)=0;
0059 end 
0060 v=A\b; % solve for solution in v labeling
0061 w=reshape(v(1:mn),m,n); %translate from v to w
0062 w=w+20;
0063 T = max(max(w));
0064 %figure(1);
0065 %mesh(x,y,w');
0066 %ylabel('Heat Input (W)');