2
$\begingroup$

I'm using NDSolve to simulate the following model of phase separation: enter image description here

where rho_T=rho_1+rho_2. The term with chi causes rho_1 and rho_2 to repel each other. However, there is a numerical instability at the interface which causes the solution to blow up: enter image description here

Refining the mesh narrows the cusps, but they still diverge and cause the integration to halt. How can I use NDSolve to solve this system?

Here is the code in Mathematica 13.1. Note that decreasing chi allows NDSolve to converge, but you end up with a uniform solution instead of two domains.

rhoT = rho1[t, x] + rho2[t, x];

chi = 10.5;

vars = {{rho1[t, x], rho2[t, x]}, t, {x}};
op = D[vars[[1]], t] + 
   DiffusionPDETerm[
    vars, {{alpha1/rhoT, 
      alpha1*(rho1[t, x]/rhoT)*chi}, {alpha2*(rho2[t, x]/rhoT)*chi, 
      alpha2/rhoT}}];

(*no-flux boundary conditions*)

bcConserved = NeumannValue[0, x == 0] + NeumannValue[0, x == l];

alpha1 = 1;
alpha2 = 1;

tF = 1;
e0 = 0.1;
l = 1;

(*initial conditions*)

f1[x_] := (e0/
     l)*(20/(Log[20*(1 + E^(10*l))] - Log[20*(1 + E^(-10*l))]))/(1 + 
      E^(20*(x - l/2)));

f2[x_] := (e0/
     l)*(-(20/(Log[20*E^(10*l) (1 + E^(10*l))] - 
           Log[20*(1 + E^(10*l))]))/(1 + E^(20*(x - l/2))) + 2);

ic1 = rho1[0, x] == f1[x];

ic2 = rho2[0, x] == f2[x];


sol = NDSolveValue[{op == {bcConserved, bcConserved}, ic1, 
    ic2}, {rho1, rho2}, {t, 0, tF}, {x, 0, l}];

{r1, r2} = sol;

Plot[{r1[tF, x], r2[tF, x]}, {x, 0, l}, PlotLegends -> {"1", "2"}, 
 PlotRange -> {Automatic, Full}, LabelStyle -> Directive[Black, 10]]
Improve this question
$\endgroup$
8
  • $\begingroup$ This is a typical problem with pure Neumann boundary conditions discussed on mathematica.stackexchange.com/questions/191476/… $\endgroup$
    – Alex Trounev
    Jan 13 at 4:33
  • $\begingroup$ @AlexTrounev, how do get to this conclusion? This is a time dependent problem, while the one you link to is a stationary problem. $\endgroup$
    – user21
    Jan 13 at 5:52
  • $\begingroup$ The second equation is independent from the first, is that right? Where does this equation come from? $\endgroup$
    – user21
    Jan 13 at 5:57
  • $\begingroup$ @user21 There is a typo in the Latex form equation for $\rho_2$. Equations should be symmetrical. But in the code there is the symmetry. $\endgroup$
    – Alex Trounev
    Jan 13 at 6:48
  • 1
    $\begingroup$ @user21 It's a mass-transfer equation where rho is the mass density and the flux J=-alpha * concentration * gradient of chemical potential. The concentration is rho/rhoT. The chemical potential is mu_i =ln(rho_i)+chi * rho_j. (I'm playing fast and loose with units here, but that's not the source of the problem.) It's very similar to Eq. 10 in this paper, although I'm working on a slightly different problem. pubs.rsc.org/en/content/articlelanding/2019/sm/c8sm02045k $\endgroup$
    – ben.w
    Jan 13 at 16:37

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy