I'm getting the error "Unable to perform assignment because the size of the left side is 1-by-1 and the size of the right side is 1-by-31." The function Tn1(:,i)=t​ridiagonal​(b,a,c,d); where there error lies, Can someone help me with this.

%% Geometry information

L= 0.10; %Length of Fin, meters

w= 0.20; %Width of Fin, meters

tt=0.01; %Full thickness of the Fin, meters

t= @(y)(tt/L)*y; %Thickness of Triangular Fin is a function of the y direction, meters

%% Property information

k= 200; %thermal conductivity, W/m-K

%% Boundary condition information

BCx0=0;%BC Flag at x=0,1 = constant temperature, 0= Specified Heat Transfer Coefficient(Convection), or Heat Flux

BCy0=0;%BC Flag at y=0,1 = Constant temperature, 0= Specified Heat Transfer Coeffieient(Convection), or Heat Flux

Tb= 120; %base temperature at y=0, deg C

%% Source term information

h=@(x)7*x^(-1/2); %convection coefficient equation

Tinf= 30; %freestream temperature, deg C

%% Discretization information

M= 41; %number of nodes in the x direction

N= 31; %number of nodes in the y direction

dx=w/(M-1); %spacing between nodes, m

x=[0:dx:w];

dy=L/(N-1); %spacing between nodes, m

y=[0:dy:L];

%% Determine cross section and surface area for each CV

ds=zeros(length(y),1);

dAs=zeros(length(y),1);

Ac=zeros(length(y),1);

for j=2:N-1

t1=(t(y(j))+t(y(j+1)))/2; %average thickness of CV on side closer to y=0

t2=(t(y(j-1))+t(y(j)))/2; %avgerage thickness of CV on side closer to y=L

ds(j)=sqrt((t(y(j+1))-t(t(j-1))).^2+(2*dy).^2)/2; %using central difference

dAs(j)= (w+((t1+t2)/L))*ds(j)*dx; %surface area

Ac(j)=(t(y(j)))*dx; %Cross sectional area

end

%Surface Area and Area for the j = 1 location of CV (Different because it's area at the domains edge)

t1=t(y(1));

t2=(t(y(1))+t(y(2)))/2;

ds(1)=sqrt((t(y(2))-t(y(1))).^2+(dy/2).^2);

dAs(1)=(w+((t1+t2)/L))*ds(1);

Ac(1)=(t(y(1)))*dx;

%Surface Area and Area for the j = N location of CV (Different because they area at the domains Tip)

t1=(t(y(end-1))+t(y(end)))/2;

t2=t(y(end));

ds(end)=sqrt((t(y(end))-t(y(end-1))).^2+(dx/2).^2);

dAs(end)=(w+((t1+t2)/L))*ds(end);

Ac(end)=(t(y(end)))*dx;

%% Start solution

Tn=Tinf; %initial guess (n) for domain temperature

Tn1=Tn; %initialize new iteration ("n+1") temperature

% Initialize tridiagonal vectors

a=zeros(N-1,1);

b=zeros(N,1);

c=zeros(N-1,1);

d=zeros(N,1);

% Loop through iterations, indexed by n

for n=1:100 %upper limit on iterations

% Sweep in positive x direction first, solving all j values

% simultaneously using tridiagonal solver

for i=1:M

%fill out first element of tridiagonal vectors

if BCy0==1

b(1)=1;

c(1)=0;

d(1)=Tb;

else

% a(1) = is absent because there is no j-1 term in the tridiagonal vector

b(1)= -((1/dy)*2*k*Ac(1+1)*k*Ac(1)/(k*Ac(1+1)+k*Ac(1))+(1/dx)*(k*(1/2*(t(1-1/2)+t(1+1/2))*dy))+(1/dx)*(k*(1/2*(t(1-1/2)+t(1+1/2))*dy))+h(x(i))*dAs(1));%Anorth+Awest+Aeast + Convection wrt x @base

c(1)= (1/dx)*(k*(1/2*(t(1-1/2)+t(1+1/2))*dy));%Aeast @base

d(1)=-((1/dx)*(k*(1/2*(t(1-1/2)+t(1+1/2))*dy))*(i-1)+(1/dx)*(k*(1/2*(t(1-1/2)+t(1+1/2))*dy))*(i+1)+(h(x(i))*dAs(1)*Tinf));%Awest+Aeast + Convection wrt x @base

end

%fill out middle of tridiagonal vectors( **Boundary Nodes are already found from the boundary conditions, and the

%for loop stops at j = N-1 because there is a boundary at N, and since the area of the fin tip is basically nothing

%I'm going to assume that the temperature is equal to the free stream temperature directly at the fin tip**)

for j=2:N-1

a(j-1)= (1/dy)*2*k*Ac(j)*k*Ac(j-1)/(k*Ac(j)+k*Ac(j-1)); %Asouth, Harmonic mean of Areas

c(j)= (1/dy)*2*k*Ac(j+1)*k*Ac(j)/(k*Ac(j+1)+k*Ac(j)); %Anorth, harmonic mean of Areas

b(j)= -(a(j-1)+c(j)+(1/dx)*(k*(1/2*(t(j-1/2)+t(j+1/2))*dy))+(1/dx)*(k*(1/2*(t(j-1/2)+t(j+1/2))*dy))+h(x(i))*dAs(j)); %Asouth+Anorth+Awest+Aeast + Convection wrt x

d(j)= -((1/dx)*(k*(1/2*(t(j-1/2)+t(j+1/2))*dy))*(i-1)+(1/dx)*(k*(1/2*(t(j-1/2)+t(j+1/2))*dy))*(i+1)+(h(x(i))*dAs(j)*Tinf)); %Aeast+Awest + Convection wrt x

end

%fill out last element of tridiagonal vectors (**Pretty Sure this is the

%Fin Tip Node FDE**)

%c(end) = is absent because there is no j+1 term in the tridiagonal vector

a(end)=(1/dx)*(k*(1/2*(t(end-1/2)+t(end))*dy));%Awest @Fin tip

b(end)= -((1/dy)*2*k*Ac(end)*k*Ac(end-1)/(k*Ac(end)+k*Ac(end-1))+(1/dx)*(k*(1/2*(t(end-1/2)+t(end))*dy))+(1/dx)*(k*(1/2*(t(end-1/2)+t(end))*dy))+(h(x(i))*dAs(end)*Tinf)); % Asouth+Awest+Aeast + Convection wrt x @Fin tip

d(end)= -((1/dx)*(k*(1/2*(t(end-1/2)+t(end))*dy))*(i-1)+(1/dx)*(k*(1/2*(t(end-1/2)+t(end))*dy))*(i+1)+(h(x(i))*dAs(end)*Tinf));% Awest+Aeast + Convection wrt x @fin tip

%solve tridiagonal matrix to update Tn+1 at current i location

Tn1(:,i)=tridiagonal(b,a,c,d);

end

%Initialize past solution at n with recently solved values at n+1 (** N+1 step is done

% after every i,j node is visited, it becomes the new Tn**)

Tn=Tn1;

[Alan Stevens]

You haven't supplied your tridiagonal function, so it's not possible to check.

[Matthew Williams]

Here is the tridiagonal function that I used in my code formulation

function x = tridiagonal(diag,lower,upper,rhs)

% tridiagonal: solves a tridiagonal matrix.

%

% Synopsis: x = tridiagonal(diag,lower,upper,rhs)

%

% Input: diag = diagonal vector of length N

% lower = lower diagonal vector of length N-1

% upper = upper diagonal vector of length N-1

% rhs = right hand side vector of length N

%

% Output: x = solution vector of length N

%

% --- determine the number of equations

N = length(diag);

% --- eliminate the lower diagonal

for i = 2:N

factor = lower(i-1)/diag(i-1);

diag(i) = diag(i)-upper(i-1)*factor;

rhs(i) = rhs(i)-rhs(i-1)*factor;

end

% --- back substitute

x(N) = rhs(N)/diag(N);

for i = N-1:-1:1

x(i) = ( rhs(i)-upper(i)*x(i+1) ) / diag(i);

end

end

Answer 1

[Alan Stevens]

You should probably define

Tn=Tinf*ones(N,M);

and then set

Tn1(:,i)=tridiagonal(b,a,c,d)';

However, there are other problems after that (the code works, but the results seem to be nonsense!).