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

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

[Sudarshan]

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?

Answer 1

[Andrea Somma]

I managed to solve it with a complete different approach and non-dimensional coordinates

[Torsten]

If you give this as an answer, you also have to include the code ... :-)

[Andrea Somma]

I will include the whole script after I will publish it for my project, I promise.

[Andrea Somma]

https://github.com/sommaa/poly_droplet