Not sure why my loop continues? Solving a diffusion equation for context. Continues to print out time but should stop when conditions are met.
% Define constraints
alpha = 0.001; % cm^2/s
height = 2; % cm
dx = 0.1; % cm
dt = 0.1; % s
% Boundary conditions
Tbottom = 130; % C degrees
Tcondition=60; %degrees C
% Grid spacing
npoints = int32(height/dx)+1;
% Mesh
T=ones(1,npoints)*4.0; %initial temp is 4
disp(T)
time=0;
condition=true;
while condition
time=time+dt;
fprintf('Current Time is: %f\n', time);
Tnew=ones(1,npoints);
T(1)=Tbottom;
for i=2:npoints-1
Tnew(i)=(T(i))+alpha*dt/dx/dx * (T(i+1)+T(i-1)-2*T(i));
end
i=npoints-1;
Tnew(i)=T(i)-alpha*dt/dx/dx * (T(i)-T(i-1));
Tmin=inf;
for i=1:npoints
T(i)=Tnew(i);
if T(i)<Tmin
Tmin=T(i);
end
end
if Tmin>Tcondition
condition=false;
end
end
disp(T);

3 commentaires

Torsten
Torsten le 8 Mar 2024
Modifié(e) : Torsten le 8 Mar 2024
Which boundary condition did you try to set at x = height ?
In your code, your work with T = 4, and you set a strange update at x = height - dx:
i=npoints-1;
Tnew(i)=T(i)-alpha*dt/dx/dx * (T(i)-T(i-1));
And note that you never set Tnew(1) to a reasonable value in the while-loop.
Jalyn-Rose
Jalyn-Rose le 8 Mar 2024
I used the x=height for my grid spacing not for the boundary. My only boundaries are the T=130 and 60. T=4 is the initial temp across right side of my diffusion.
Torsten
Torsten le 8 Mar 2024
To answer your question:
Your while loop never stops because Tnew(1) is always 1 from the setting Tnew=ones(1,npoints);, thus < Tcondition.

Connectez-vous pour commenter.

 Réponse acceptée

Torsten
Torsten le 8 Mar 2024
Modifié(e) : Torsten le 8 Mar 2024
The initial condition at t = 0 seems to be T = 4 for all x in your code.
But you must set boundary conditions for T at x = 0 and x = 2.
At x = 0, this seems to be T = 130 for all t. But what about x = 2 ?
% Define constraints
alpha = 0.001; % cm^2/s
height = 2; % cm
dx = 0.1; % cm
nx = height/dx+1;
% Boundary conditions
Tbottom = 130; % C degrees
Tcondition=60; %degrees C
dt = 0.1; % s
nt = 5000;
T = 4*ones(nx,1);
T(1) = Tbottom;
Tnew = zeros(size(T));
condition = true;
counter = 0;
hold on
while condition
counter = counter + 1;
Tnew(1) = T(1); % T = 130 at x = 0
Tnew(2:nx-1) = T(2:nx-1)+dt*alpha/dx^2*(T(1:nx-2)-2*T(2:nx-1)+T(3:nx));
Tnew(nx) = Tnew(nx-1); %dT/dx = 0 at x = 2
T = Tnew;
condition = min(T) < Tcondition;
if counter == 500
counter = 0;
plot((0:nx-1)*dx,T)
end
end
hold off
grid on

7 commentaires

Jalyn-Rose
Jalyn-Rose le 8 Mar 2024
ok i meant to have At x=2 I want the T=4
Torsten
Torsten le 8 Mar 2024
But then 4=Tmin < Tcondition for all t, and you will never quit the while loop...
Jalyn-Rose
Jalyn-Rose le 8 Mar 2024
Hmmm I want the while loop to continue until the Tcondition is met so I think i have my sign switched. Starting at x=0,T=130 until x=2, T=60 and and display the time it takes to get there from x=2, T=4.
Torsten
Torsten le 8 Mar 2024
Modifié(e) : Torsten le 8 Mar 2024
Starting at x=0,T=130 until x=2, T=60 and and display the time it takes to get there from x=2, T=4.
That's what my code does. The only thing you have to do is add the dt's in the while loop.
Jalyn-Rose
Jalyn-Rose le 8 Mar 2024
I see, thank you for the help!
Jalyn-Rose
Jalyn-Rose le 8 Mar 2024
Im uncertain of how I would see the time?
time = 0;
while ...
time = time + dt;
...
end

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Loops and Conditional Statements dans Centre d'aide et File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by