For loop for iteration.
6 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Rik
am 2 Sep. 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.
Antworten (0)
Siehe auch
Kategorien
Mehr zu Computational Fluid Dynamics (CFD) finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!