# Problem with C-H +Stokes + viscosity coupling

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 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 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)( 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);

}``````

Your Peclet number is huge. Try with a smaller Pe value or dramatically increase the resolution of your mesh. Also, you should use an inf-sup stable finite element discretization like P2-P1.

Thank you so much for the reply! I’ll try out those suggestions. If I can ask one more thing, did I add the advection correctly considering I’m using Newton iterations? I don’t have much experience with it so I’m not sure if I should add anything to vhJ?

EDIT: P2-P1 discretization and bigger mesh resolution fixed the problem. Thank you!

1 Like

//Mesh
mesh Th=square(100,100);
//border C(t=0,2pi){x=cos(t);y=sin(t);};
//mesh Th=buildmesh (C(50));
plot(Th,wait=1);
macro uex(t)(exp(-2
t)(sin(pix)^2)(sin(piy)^2)) //

//Fespace
fespace Vh(Th,P1);
fespace Vh2(Th,[P1,P1]);
//Vh uex;

//Parameter
real a=0.01;
real dt=0.10;
int i,k;
//real t;

// Function
func real df(real u){return u^3-u; }
func real ddf(real u){return 3*u^2-1; }

Vh2 [u,w],[v,z],[oldv,oldz],[auxv,auxz],[vv,zz];
Vh dfalpha;
Vh ddfalpha;
//Vh uex;

//macro uex(t) (exp(-2t)(sin(pix)^2)(sin(pi*y)^2)//
// Problem
varf vdJ([v,z],[phi,psi])=int2d(Th)
((oldv-u)phi+dt(dx(oldz)dx(phi)+dy(oldz)dy(phi))+oldzpsi-
a
(dx(oldv)*dx(psi)+dy(oldv)dy(psi))-dfalphapsi);

varf vhJ([v,z],[phi,psi])=int2d(Th)
(vphi+dt(dx(z)dx(phi)+dy(z)dy(phi))+zpsi-a(dx(v)dx(psi)
+dy(v)dy(psi))-ddfalphav
psi);

//u=uex(t);

[u,w]=[(sin(pix)^2)(sin(piy)^2),-2a*(pi^2)(sin(pi^2(x^2+y^2))
-2*(pi^2)(sin(pix)^2)(sin(piy)^2)(x^2+y^2))+((sin(pix)^2)^3)((sin(piy)^2)^3)
-(sin(pix)^2)(sin(pi*y)^2)];
//u=uex;

plot(u,wait=1,cmm=“Cahn-Hilliard” ,value=true);
//u=uex(t);
//macro uex(t) (exp(-2t)(sin(pix)^2)(sin(piy)^2))//
// u=uex(t);
//uex(t)=exp(-2
t)(sin(pix)^2)(sin(piy)^2);
//real t=0;// start from t=0
real t;
//t=0.1;
for (int i=0;i<1./dt; i++) // Time loop
{
t=t+dt;
//u=uex(t);
//=exp(-2t)((sin(pix)^2)(sin(piy)^2));
{
oldv[]=u[];
for (int k=0;k<10;k++)// Newton loop
{
dfalpha=df( oldv );
ddfalpha=ddf( oldv );
auxv[]=vdJ(0,Vh2);
real res= auxv[]'auxv[]; // Residual
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[]; // Newton
v[] -=vv[];
oldv[]=v[];
}
u[]=v[];
//u=uex(t);
//uex=exp
sin(pix)^2sin(piy)^2;
plot(u,wait=1,cmm=“Cahn-Hilliard”,value=true);
//uex(t)=exp(-2
t)sin(pix)^2sin(piy)^2;
//plot(uex(t),wait=1,cmm=“Cahn-Hilliard exact solution”,value=true);

``````}
cout << "t" << t << "L^2-error=" << sqrt(int2d(Th)((u-uex(t))^2))<<endl;
//cout << "t"<< t << "L^2-error=" << sqrt(int2d(Th)
//(abs(u-(exp(-2*t)*(sin(pi*x)^2)*(sin(pi*y)^2)))^2))<<endl;
}
``````

Can you tell me Please what is wrong in this code. Boundary of the approximate solution not matching with exact solution. Please help. It is very urgent.

Can you tell me how i post a queries in this platform?