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 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;
    return Outside;

mesh Th=square(72,72);
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)

varf vhJ([v,z],[phi,psi]) = int2d(Th)


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 );
          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);
          v[] -= vv[];

      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 Th=square(100,100);
//border C(t=0,2pi){x=cos(t);y=sin(t);};
//mesh Th=buildmesh (C(50));
macro Grad(u)[dx(u),dy(u)]//
macro uex(t)(exp(-2
t)(sin(pix)^2)(sin(piy)^2)) //

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

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)

varf vhJ([v,z],[phi,psi])=int2d(Th)



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

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?