Hi! I’m trying to simulate Cahn-Hilliard equation using explicit Euler scheme. However, no matter what changes I make and what parameters I choose I can’t get it to work. Can I use this scheme or does the non-linearity of the equation strictly require Newton iterations or some similar method? I also tried Euler semi-implicit method and implicit but nothing works. Any help is greatly appreciated.
real radius = 0.05;
real M=1;
real centerX = 0.35;
real centerX2 = 0.5;
real centerY = 0.5;
real Inside = 1;
real Outside = -1;
real Pe= 62;
real Ch=0.01;
func real initialCondition(real x, real y)
{
if (sqrt(pow(x - centerX, 2) + pow(y - centerY, 2)) <= radius ||
sqrt(pow(x - centerX2, 2) + pow(y - centerY, 2)) <= radius ||
sqrt(pow(x - centerX, 2) + pow(y - 0.8, 2)) <= radius ||
sqrt(pow(x - centerX2, 2) + pow(y - 0.8, 2)) <= radius)
return Inside;
else
return Outside;
}
mesh Th=square(72,72);
plot(Th,wait=1);
fespace Vh(Th,P1);
fespace Vh2(Th,P1);
real dt=0.01;
int i,k;
Vh2 u,w,phi,oldv,oldz,psi;
u=initialCondition(x,y);
plot (u,wait=1,cmm="Cahn-Hilliard",value=true);
for (int i=0;i<30;i++)
{oldv[]=u[];
solve CahnHill1(w,phi) =
int2d(Th)(Ch*Ch*(dx(oldv)*dx(phi)+dy(oldv)*dy(phi)) )+ int2d(Th)((oldv^3-oldv)*phi) -int2d(Th)(w*phi);
solve CahnHill2(u,phi)= int2d(Th)(u*phi) + int2d(Th)(dt*1/Pe*(dx(w)*dx(phi)+dy(w)*dy(phi)) ) -int2d(Th)(oldv*phi);
plot (u,wait=1,fill=true, cmm="Cahn-Hilliard_"+i,value=true);
}
Probably the biggest problem is the time step is too big. There may be others.
The attached seems to run just slightly out of bounds as min and max exceed
}!}. You can do w and u si,multaneously unless your analysis suggested they
may be better in sequence. Of course generally you want to do more things
at once. Uinsg “old” is a concession to the linearlity constrtain. You may be able
to replace some of the “oldv” things with “u” but again you should look at the
discretized form details.
See if this is any more of what you expect I changes the initial conditions to a
random thing and you can watch it segregate somewhat. I thought there
may be an issue with length scale and see also the slightly different
formulation which produce similar “u” as your’s except I probably dropped
some scale factors.
Hi! Thank you for taking a look at my code and making the changes you did, it’s really working now! I tried putting everything into one solve command instead of two and fixing the min and max bounds but I still got a solution that diverged so I though the main problem was somewhere in the formulation. Now I see from your suggestions that it was all of those issues + time step that was too large+ some other issues you fixed. Thank you once again for your help!
If you are doing any serious work you should probably actually
try to find some analytical cases to test. Should the solution stay bounded?
When do the time derivatives go to zero? etc.
Once I add Stokes equation to this one I can compare the results to an article I found online. So I’ll see how that goes. But maybe I can try to find one with just the C-H.
Just skimming it one comment was about the interface thickness
and this comes up in a lot of places like with electrochemistry
and etching a solid in a liquid as a lot goes on in a few nm
of the interface.