For loop for iteration.
6 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Hi all.
I try to solve 1D diffusion problem for diffusion through a rod. ∇⋅(-k∇T) = 0
k is dependent of the temperature
k=k0+alpha*T
Since k (k1) isn't constant and thereby the equation isn't linear I will have to make iterations.
My issue is though that my outer iteration forloop doesn't seem to work. It's like its not overwriting the initial temperature.
Can someone please help me out on this?
clear all
close all
% Simple 1D diffusion problem for Dirichlet BC's
% Steps of solution:
% 1-Write individual equations for control volumes
% 2-Correct for boundary conditions if necessary
% 3-Assemble Matrix
% 4-Solve Matrix
% 5-Plot and analyze results
p=4;
% Inputs
Ta=100; %left end boundary
Tb=500; %right end boundary
Area=10*10^(-3);
k=1000;
L=0.5;
N=10; %amount of elements
alpha=10;
%amount of iterations
x=linspace(0,L,N);
T1=(Ta+(Tb-Ta)/L*x)';
% Processing
dx=L/N;
a=zeros(N,N);
b=zeros(N,1);
for i=1:p
k1(:,i)=k+T1(:,i)*alpha;
for j=2:1:N-1
ae=k1(j,i)*Area/dx;
aw=k1(j,i)*Area/dx;
ap=ae+aw;
a(j,j-1)=-aw;
a(j,j)=ap;
a(j,j+1)=-ae;
b(j,1)=0;
end
% In this code the boundaries are hardwired.
j=1;
aw=2*k1(j,i)*Area/dx;
ae=k1(j,i)*Area/dx;
ap=ae+aw;
a(j,j)=ap;
a(j,j+1)=-ae;
b(j,1)=aw*Ta;
% In this code the boundaries are hardwired.
j=N;
aw=k1(j,i)*Area/dx;
ae=2*k1(j,i)*Area/dx;
ap=ae+aw;
a(j,j-1)=-aw;
a(j,j)=ap;
b(j,1)=ae*Tb;
% Assembly complete
% spy(a)
T2=(a^-1*b)';
T1(:,i+1)=T2;
end
xsol=dx/2:dx:L-dx/2;
xexact=0:L/100:L;
Texact=100+800*xexact;
plot(xsol,T1(:,i),'bs',xexact,Texact,'r')
xlabel('x-axis')
ylabel('Temperature')
legend('Numerical','Exact')
2 comentarios
Rik
el 2 de Sept. de 2019
Start by removing clear all from your code. If you want to avoid variables persisting across runs you should use clearvars during debugging and functions for normal use.
Then you can use breakpoints to step through your code and see when something unexpected happens.
Respuestas (0)
Ver también
Categorías
Más información sobre Computational Fluid Dynamics (CFD) en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!