I am trying to simulate the phase separation of a binary mixture. If the free energy F is known as a function of the concentration $c$, the dynamical equation is:
$ \frac{\partial c(x,t)}{\partial t}=\frac{d^2}{dx^2} \frac{\delta F[c]}{\delta c} $
For the Flory-Huggins free energy we have:
$ \frac{\delta F[c]}{\delta c}=\log\left(\frac{c}{1-c}\right) + \chi c^2 + \gamma (c')^2 $
First term is entropy, second term is attraction between particles ($\chi<0$) and third term is similar as a surface tension.
I use time forward and space centered differentiation. Even for very small $dt$ I get numerical instabilities. I first thought the term $\gamma(c'^2)$ was responsible, but instabilities remain even without.
Here is my Matlab code for no-flux boundary conditions, clear and simplified as much as I can. I am aware the boundary condition implementation may not be correct but I don't think the problem comes from this.
What should I do?
function phaseSep ()
clear all;
clf;
N = 101; %number of grid points in x
dx = 1/(N - 1);
x = 0:dx:1; %vector of x values
T = 1e3; %number of time steps
dt = 1e-8;
% Second derivative
function y=mylaplace(f,i)
y = f(i+1)-2*f(i)+f(i-1);
end
function y=dfdc(C)
p=-0.01;
g=0.01;
for i=2:N-1
y(i) = log(C(i)/(1-C(i)))+p*C(i)+g*mylaplace(C,i)/dx^2;
end
% No-flux Boundary conditions
y(1)=y(2);
y(N)=y(N-1);
end
%Initial concentration
for i=1:N
C(i)=0.2+0.1*tanh(10*(x(i)-0.5));
end
% Plot initial concentration
hold;
plot(x,C);
% iterate
for t=1:T
Z=dfdc(C);
for i=2:N-1
C(i) = C(i) + (mylaplace(Z,i)/dx^2)*dt;
end
%No flux boundary conditions
C(1) = C(2);
C(N) = C(N-1);
end
% Plot final concentration
plot(x,C);
end