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

Question

Andrea Somma 2022-11-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1840928-1d-simulation-diffusion-term-polar-coordinates-moving-boundary-condition

评论： Andrea Somma 2022-12-19

Paper A.pdf

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Sudarshan 2022-11-4

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?

请先登录，再进行评论。

请先登录，再回答此问题。