フィルターのクリア

Iteration coded weird behaviour

2 ビュー (過去 30 日間)
Dejan Cvijanovic
Dejan Cvijanovic 2012 年 5 月 8 日
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 件)

カテゴリ

Help Center および File ExchangeLogical についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by