Iteration coded weird behaviour
1 visualizzazione (ultimi 30 giorni)
Mostra commenti meno recenti
Sorry to bring up my old code again. The problem I'm having is the radius [a(i)] decreases to 0 (as it should), but the concentration [c(i,j)] decreases, in general, but unfortunately increase on some individual time steps (e.g; 0.5 0.4 0.45 0.3 0.32 0.12 0. 18 ...), when it should be strictly decreasing (e.g; 0.5 0.4 0.3 ....)
...I'm thinking the problem is either in c_ah or the definition for c(i,j+1). Any thoughts?
alpha = 0.5; beta = 1; a_0 = 1; dr = 0.1; dt = input('Value of dt:'); h = 0.0001; c_inf = 0; c_a = 1; c_s = 2; mesh_points = 3; rho = dt/(dr)^2; R= 10; N_i = R/dr; t_diss = 0.50; N_t = t_diss/dt;
if rho > 0.5; disp('Invalid dt value'); return; end
r = (0:N_i) * dr; c = zeros(N_i, N_t); c(:, 1) = c_inf; c(r < a_0, 1) = c_s;
a = zeros(1,N_t); a(1) = a_0;
dr3 = dr ^ 3;
for j = 1:N_t -1
if a(j) <= 0
a(j) = 0;
ajp1 = 0;
a(j+1) = 0;
else
b_1 = floor (a(j)/dr) + 3;
b_2 = b_1 + 1;
b_1s = (b_1 - 1)*dr;
b_2s = (b_2 - 1)*dr;
% c_ah = (a(j)+h-b_1s)*(a(j)+h-b_2s)*(c_a)/((a(j)-b_1s)*(a(j)-b_2s)) + ...
% h*(a(j)+h-b_2s)*c(b_1,j)/((b_1s-a(j))*(-1)) + ...
% h*(a(j)+h-b_1s)*c(b_2,j)/((b_2s-a(j))*(1));
c_ah = (a(j)+h-b_1s)*(a(j)+h-b_2s)*(c_a)/((a(j)-b_1s)*(a(j)-b_2s)) + ...
h*(a(j)+h-b_2s)*c(b_1,j)/((b_1s-a(j))*(b_1s-b_2s)) + ...
h*(a(j)+h-b_1s)*c(b_2,j)/((b_2s-a(j))*(b_2s-b_1s));
dc_dr = ((c_ah)-c_a)/h;
ajp1 = a(j) +dt*beta*dc_dr;
a(j+1) = ajp1;
end
tmp1 = (a(j)^2) * (alpha-1) * (ajp1-a(j)) / (2 * dr3);
for i = 2:N_i - 1
if r(i) > ajp1
tmp2 = rho / (i-1) + tmp1 / ((i-1)^2);
d = rho - tmp2;
e = 1 - 2 * rho;
f = rho + tmp2;
c(i, j+1) = d*c(i-1,j) + e*c(i,j) + f*c(i+1, j);
elseif r(i) == ajp1
c(i, j+1) = c_a;
else
c(i, j+1) = c_s;
end
disp(['c(' num2str(i) ',' num2str(j) ') = ' num2str(c(i,j))]);
disp(['a(j) = ' num2str(a(j))]);
end
c(1, j+1) = c(2, j+1);
c(N_i, j+1) = c_inf;
end
0 Commenti
Risposte (0)
Vedere anche
Categorie
Scopri di più su Logical in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!