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)=tridiagonal(b,a,c,d); where there error lies, Can someone help me with this.
3 次查看(过去 30 天)
显示 更早的评论
%% 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;
2 个评论
Alan Stevens
2020-12-14
You haven't supplied your tridiagonal function, so it's not possible to check.
回答(1 个)
Alan Stevens
2020-12-14
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!).
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Operating on Diagonal Matrices 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!