0
$\begingroup$

I'm having a hard time refiguring this code we used in class to do a 1D finite-difference solution. If anyone could help it would be greatly appreciated. The problem set is listed below as well as the code.

enter image description here

clear all
close all
clc
tic

m=161; % number of y nodes 
n=81; % number of x nodes

Lx=0.5;     % length of plate in x-direction [m]
Ly=1;       % length of plate in y-direction [m]
kappa=0.6;   % thermal conductivity [W/m/K]

dx=Lx/(n-1);
dy=Ly/(m-1);

if dx~=dy
error("dx and dy must be equal")
end

A=zeros((m-2)*(n-2),(m-2)*(n-2));
B=zeros((m-2)*(n-2),1);
%row=1;

TL=zeros(n,1);
TR=zeros(n,1);
TB=zeros(1,m);
TT=zeros(1,m);

H=zeros(m,n);
H((m-1)/2+1,(n-1)/2+1)=1e4;
H(5,(n-1)/2+1)=1e4;

rowMap=[];
for i=1:m
    for j=1:n
        if i>=2 && i<=m-1
            if j>=2 && j<=n-1

                row=(n-2)*(i-2)+(j-2)+1;
                rowMap=[rowMap;row i j];
                B(row)=-dx^2/kappa*H(i,j);

                A(row,row)=-4;

                if(i-1>=2 && i-1<=m-1)
                    col=(n-2)*(i-1-2)+(j-2)+1;
                    A(row,col)=1;
                elseif i-1==1
                    B(row)=B(row)-TL(j);
                elseif i-1==m
                    B(row)=B(row)-TR(j);
                end

                if(i+1>=2 && i+1<=m-1)
                    col=(n-2)*(i+1-2)+(j-2)+1;
                    A(row,col)=1;
                elseif i+1==1
                    B(row)=B(row)-TL(j);
                elseif i+1==m
                    B(row)=B(row)-TR(j);

                end

                if(j-1>=2 && j-1<=n-1)
                    col=(n-2)*(i-2)+(j-1-2)+1;
                    A(row,col)=1;
                elseif j-1==1
                    B(row)=B(row)-TB(i);
                elseif j-1==n
                    B(row)=B(row)-TT(i);
                end


                if(j+1>=2 && j+1<=n-1)
                    col=(n-2)*(i-2)+(j+1-2)+1;
                    A(row,col)=1;
                elseif j+1==1
                    B(row)=B(row)-TB(i);
                elseif j+1==n
                    B(row)=B(row)-TT(i);
                end

%                row=row+1;

            end

        end
    end
end
disp("Assembly time")
toc

%A;
%eig(A)
S=sparse(A);
fillRatio=nnz(S)/prod(size(A))
Tv=S\B;
disp("Sparse Solve time")
toc

% Solve 
Tv=A\B;
disp("Dense Solve time")
toc

% Construct T matrix for plotting
T=zeros(m,n);
T(1,:)=TL;
T(end,:)=TR;
T(:,1)=TB;
T(:,end)=TT;

for k=1:size(rowMap,1)
    i=rowMap(k,2);
        j=rowMap(k,3);
    T(i,j)=Tv(k);
end

%% Plot
surf([1:n]/n*Lx,[1:m]/m*Ly,T)
xlabel('X [m]')
ylabel('Y [m]')
zlabel('\Delta T [K]')
axis equal
$\endgroup$
  • $\begingroup$ What have you tried? what are you stuck on? Do you have a specific question? $\endgroup$ – Daniel K Dec 5 '19 at 14:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.