Hi! I’m trying to simulate the mixing of two liquids of different viscosities using Cahn-Hilliard and Stokes equations. Due to the non-linear nature of C-H equation I’m using Newton iterations. If the viscosity ratio (lambda) is 1 the code works without any problems but for any other ratio after around 100 iterations I start getting these minimums and maximums that spread and take over the entire domain as shown in the picture. Any ideas what’s wrong with my code? Thanks!
real radius = 0.05;
real M=1;
real centerX = 0.35;
real centerX2 = 0.5;
real centerY = 0.5;
real Inside = 1.0;
real Outside = -1.0;
real Pe= 10000;
real Ch=0.01;
func real initialCondition(real x, real y)
{
if (y<0.5)
return Inside;
else
return Outside;
}
mesh Th=square(72,72);
plot(Th,wait=1);
fespace Vh(Th,P1);
fespace Vh2(Th,[P1,P1]);
real lambda=10;
real a=0.01;
real dt=0.01;
int i,k;
func real df(real u) { return u^3-u; }
func real ddf(real u) { return 3*u^2-1; }
Vh2 [u,w],[v,z],[oldu,oldz],[auxv,auxz],[vv,zz];
Vh dfalpha; // f'(u)
Vh ddfalpha; // f"(u)
fespace Uh(Th, P1);
Uh U, V;
Uh UU, VV;
fespace Ph(Th, P1);
Ph p, pp;
varf vdJ([v,z],[phi,psi]) = int2d(Th)
((u-oldu)*phi+dt*(dx(u)*U+dy(u)*V)*phi-dt*1/Pe*(dx(oldz)*dx(phi)+dy(oldz)*dy(phi))+oldz*psi-
Ch*Ch*(dx(oldu)*dx(psi)+dy(oldu)*dy(psi))-dfalpha*psi);
varf vhJ([v,z],[phi,psi]) = int2d(Th)
(-v*phi-dt*1/Pe*(dx(z)*dx(phi)+dy(z)*dy(phi))+z*psi-Ch*Ch*(dx(v)*dx(psi)
+dy(v)*dy(psi))-ddfalpha*v*psi);
[u,w]=[initialCondition(x,y),0];
plot (u,wait=1,cmm="Cahn-Hilliard",value=true);
macro SGrad(u,v) [[dx(u),0.5*(dx(v)+dy(u))],[0.5*(dx(v)+dy(u)),dy(v)]] //
macro div(u1,u2) (dx(u1)+dy(u2)) //
macro mu(a) (0.5*(a+1)-1/lambda*0.5*(a-1)) //
for (int i=0;i<150;i++)
{
oldu[ ]=u[ ];
solve Stokes( [U,V,p],[UU,VV,pp]) =
int2d(Th)( 2* mu(oldu)*(SGrad(U,V):SGrad(UU,VV)) - div(U,V)*pp)
- int2d(Th)( div(UU,VV)*p)
+ on(3, 2, 4, U=0, V=0)
+ on(1, U=16*(x)^2*(1-x)^2, V=0)
;
for (int k=0;k<10;k++) // boucle de Newton
{
dfalpha = df( oldu ) ;
ddfalpha = ddf( oldu );
auxv[]=vdJ(0,Vh2);
real res= auxv[]'*auxv[];
cout << i << " residu^2 = " << res << endl;
matrix H;
if( res< 1e-12) break;
H= vhJ(Vh2,Vh2,factorize=1,solver=LU);
vv[]=H^-1*auxv[];
v[] -= vv[];
oldu[]=v[];
}
u[]=v[];
plot (u,wait=1,fill=true, cmm="Cahn-Hilliard"+i,value=true);
}