1D simulation diffusion term, polar coordinates, moving boundary condition

3 visualizaciones (últimos 30 días)
I am trying to integrate the system I am showing in the page 3 of the attached pdf but I am having some problems with Cw and the radius because they should reach a minimum and a maximum respectively but in my simulation they go towards +inf and -inf, and Cw is decreasing too fast. Cw0, K, Cm are parameters, I also attach the whole paper if something is not clear enough. (I am using matlab). I am using forward euler and finite difference method. It is probably better to integrate with a variable space step but I cant figure out how to do it, so I just made a big enough domain to allow the radius to increase. I have started recently cfd so be kind ahahah.
clc; clear all; close all
%% DATA
K = 1;
Diff = 4e-9;
Cw0 = 55e-3; %mol/m3
r0 = 1e-3;
L = 12e-3;
N = 200;
tf = 100;
Cm = 30e-3; %mol/m3
Cinf = 0;
%% preprocessing
h = L/N;
grid1D = linspace(0,L,N+1);
%% SOLUTION
index = find(grid1D >= r0);
dt = 0.01;
tsteps = tf/dt;
C = [Cw0.*ones(1,index(1)-1)./K,Cm*ones(1,index(end)-index(1))]';
Cplot = zeros(size(grid1D));
%loop
for t = 1:tsteps
C0 = C;
% bc 6: eq 6 integration
if t ~= 1
indexrd = find(grid1D >= rd);
for counter_bc6 = index(1):indexrd(1) - 1
int6(counter_bc6) = (-4*pi*C(counter_bc6)*grid1D(counter_bc6)^2 -4*pi*C(counter_bc6 + 1)*grid1D(counter_bc6+1)^2)*h/2;
end
Cw = sum(int6)/(4/3*pi*r0^3) + Cw0;
else
Cw = Cw0;
end
% eq 3
C(index(1)) = Cw/K;
% eq 5
if t == 1
rd = r0 + h;
else
d = -Cm + 8*Cm -8*C(indexrd(1)- 2) + C(indexrd(1)- 3);
rd = -Diff*(d)/12/h/Cm*dt + rd;
end
if (mod(rd,h)~=0)
rd = rd + (h - mod(rd,h));
end
%eq 4
indexrd_new = find(grid1D>=rd);
%expl
for i = index(1):N-1
C(i) = C0(i) + dt*(Diff/h^2*(C0(i+1) + C0(i-1) - 2*C0(i)) + 2*Diff/i/h*(C0(i+1)-C0(i-1))/h );
end
for i = indexrd_new(1) :length(C)
C(i) = Cm;
end
for i = 1:index
C(i) = Cw;
end
if (mod(t,100)==0) % => Every 50 time steps
for i=1:N-1
Cplot(i) = 0.5*(C(i+1) + C(i));
end
plot(grid1D,Cplot);
hold on
plot(grid1D(indexrd_new(1))*ones(2,1),[0.028 0.052])
%scatter(t*dt,Cw)
hold off
ylim([2.8e-2 5.2e-2])
xlim([0 L/2])
xlabel("length [m]")
ylabel("Concentration [mol/m3]")
message = sprintf('t=%d', ceil(t*dt));
time = annotation('textbox',[0.15 0.8 0.1 0.1],'String',message,'EdgeColor','k','BackgroundColor','w');
drawnow;
end
end
  1 comentario
Sudarshan
Sudarshan el 4 de Nov. de 2022
Hi Andrea,
I tried going through the paper and the code. I needed some more information on what is the result that you want to achieve as I couldnt fully understand the paper. Also, could you point out if there is any specific place in the code where I can look into?

Iniciar sesión para comentar.

Respuesta aceptada

Andrea Somma
Andrea Somma el 4 de Nov. de 2022
I managed to solve it with a complete different approach and non-dimensional coordinates
  3 comentarios
Andrea Somma
Andrea Somma el 4 de Nov. de 2022
I will include the whole script after I will publish it for my project, I promise.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Ordinary Differential Equations en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by